Particle physics - exercise 1b solution

Repository

<p dir="auto"><span><a href="https://github.com/BFuks/mad5-utopian-exercises" target="_blank" rel="noreferrer noopener" title="This link will take you away from hive.blog" class="external_link">https://github.com/BFuks/mad5-utopian-exercises <h3>Introduction <p dir="auto">This is a description of my solution provided for <a href="https://steemit.com/utopian-io/@lemouth/particle-physics-utopian-implementing-an-lhc-analysis-on-a-computer-the-physics-objects" target="_blank" rel="noreferrer noopener" title="This link will take you away from hive.blog" class="external_link">particle physics exercise<span> by <a href="/@lemouth">@lemouth <p dir="auto">The exercise code is located in two files: <pre><code>ex1b-mactro.h ex1b-mactro.cpp <p dir="auto">The goal of the exercise was to implement a simplified mechanism of identifying electrons, muons and jets useful for further study. This process consists of three parts: <ol> <li>Identifying baseline objects <li>Removal of overlapping objects <li>Identifying signal objects <p dir="auto">It is described in length in <a href="https://arxiv.org/pdf/1711.11520.pdf" target="_blank" rel="noreferrer noopener" title="This link will take you away from hive.blog" class="external_link">ATLAS research paper<br /><br /> [ Image credits <a href="https://home.cern/topics/large-hadron-collider" target="_blank" rel="noreferrer noopener" title="This link will take you away from hive.blog" class="external_link">CERN ]<span> <img src="https://images.hive.blog/768x0/https://home.cern/sites/home.web.cern.ch/files/image/inline-images/old/lhc_long_1.jpg" srcset="https://images.hive.blog/768x0/https://home.cern/sites/home.web.cern.ch/files/image/inline-images/old/lhc_long_1.jpg 1x, https://images.hive.blog/1536x0/https://home.cern/sites/home.web.cern.ch/files/image/inline-images/old/lhc_long_1.jpg 2x" /> <h3>Implementation <p dir="auto">This exercise was definitely more fun then the previous one, and required more work. I implemented 7 helper functions to get it done: <pre><code> std::vector<RecLeptonFormat> GetBaselineElectrons(const std::vector<RecLeptonFormat>& candidate_electrons); std::vector<RecLeptonFormat> GetBaselineMuons(const std::vector<RecLeptonFormat>& candidate_muons); std::vector<RecJetFormat> GetBaselineJets(const std::vector<RecJetFormat>& candidate_jets); std::vector<RecLeptonFormat> GetSignalElectrons(const std::vector<RecLeptonFormat>& baseline_electrons); std::vector<RecLeptonFormat> GetSignalMuons(const std::vector<RecLeptonFormat>& baseline_muons); std::vector<RecJetFormat> GetSignalJets(const std::vector<RecJetFormat>& baseline_jets); void RemoveOverlap(std::vector<RecJetFormat>& jets, std::vector<RecLeptonFormat>& electrons, std::vector<RecLeptonFormat>& muons); <p dir="auto">All functions except <code>RemoveOverlap() create a new vectors of Lepton or Jet format. I saw no real value in keeping vectors containing overlapping objects, so <code>RemoveOverlap() works on the vectors provided as its arguments. <p dir="auto">Bodies of all <code>GetSignal... and <code>GetBaseline.. functions look pretty similar: <pre><code>std::vector<RecLeptonFormat> test_analysis::GetBaselineElectrons(const std::vector<RecLeptonFormat> &candidate_electrons) { std::vector<RecLeptonFormat> result; for (auto& candidate : candidate_electrons) { if(candidate.pt() <= 5.7) continue; if(candidate.abseta() >= 2.47) continue; result.push_back(candidate); } return result; } <p dir="auto">Inside the <code>for loop that iterates over all provided object, we check their properties, and reject those that doesn't match. If all checks are passed, the candidate object is added to result vector. That process could be rewritten as: <pre><code> if(candidate.pt() > 5.7 && candidate.abseta() < 2.47) result.push_back(candidate); <p dir="auto">Which is shorter, but in case there are many checks to be done (which may be the case if we need to implement full procedure in some of the next exercises), then IMO it is easier to read and debug the longer version. <p dir="auto">Removing overlap is also pretty straightforward: <pre><code> // Second step from Table 3. for(auto& electron: electrons) { for(auto jet_it = jets.begin(); jet_it != jets.end(); jet_it++) { if(jet_it->dr(electron) < 0.2) { jets.erase(jet_it); // we have to move iterator back after removing the element jet_it--; } } } <p dir="auto">Outer loop iterates over elements that have a precedence (need to stay in case of overlap), while inner loop checks the elements that may be potentially removed. <p dir="auto">The next overlap removal step is implemented in a very similar way, and the rest of the code is printing of the results. <h3>Results <pre><code>=> progress: [===> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 6 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 5 => progress: [======> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 6 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 6 => progress: [==========> ] Initial electrons count: 2, initial muons count: 1, initial jets count: 9 Baseline electrons count: 2, baseline muons count: 1, baseline jets count: 9 After overlap removal: Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 9 signal electrons count: 1, signal muons count: 0, signal jets count: 8 => progress: [=============> ] Initial electrons count: 0, initial muons count: 1, initial jets count: 10 Baseline electrons count: 0, baseline muons count: 1, baseline jets count: 10 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 10 signal electrons count: 0, signal muons count: 0, signal jets count: 8 => progress: [=================> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 8 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 8 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 8 signal electrons count: 0, signal muons count: 0, signal jets count: 7 => progress: [====================> ] Initial electrons count: 0, initial muons count: 1, initial jets count: 6 Baseline electrons count: 0, baseline muons count: 1, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 6 => progress: [========================> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 4 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 4 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 4 signal electrons count: 0, signal muons count: 0, signal jets count: 4 => progress: [===========================> ] Initial electrons count: 1, initial muons count: 0, initial jets count: 6 Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 6 After overlap removal: Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 5 signal electrons count: 0, signal muons count: 0, signal jets count: 5 => progress: [===============================> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 10 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 10 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 10 signal electrons count: 0, signal muons count: 0, signal jets count: 8 => progress: [==================================>] Initial electrons count: 1, initial muons count: 0, initial jets count: 7 Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 7 After overlap removal: Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 6 signal electrons count: 1, signal muons count: 0, signal jets count: 6 => progress: [===================================] <h3>Possible improvements <ol> <li>I don't know with what number of objects such analysis should deal. Is it tens? Thousands? Or maybe billions? If the number is big, more care should be taken to optimization of the code - the overlap removal procedure seems to be really costly right now. <li>"Magic numbers" in the code should be replaced by named constants - I actually noticed that writing this article. <li>Maybe, for better readability, results could be printed only once instead for each data set. <li><span>Not an improvement, but question about requirements: <a href="/@lemouth">@lemouth Are we limited to any specific C++ standard? I'm using C++11, but maybe should give it up? <h3>Previous exercises <p dir="auto"><span><a href="https://steemit.com/utopian-io/@mactro/particle-physics-utopian-project-first-exercise" target="_blank" rel="noreferrer noopener" title="This link will take you away from hive.blog" class="external_link">https://steemit.com/utopian-io/@mactro/particle-physics-utopian-project-first-exercise <h3>EDIT <p dir="auto"><span>I did some fixes according to the comments by <a href="/@irelandscape">@irelandscape and <a href="/@lemouth">@lemouth. Both fixes are related to overlap removal: <pre><code> for(auto jet_it = jets.begin(); jet_it != jets.end();) { if(jet_it->dr(electron) < 0.2 && jet_it->true_btag()) { jet_it = jets.erase(jet_it); // Iterator now points to element after the removed one, and it will be further // incremented by the for loop, so we need to move it back. jet_it--; } } <p dir="auto">I added the <code>jet_it->true_btag() check that was missing which caused wrong results. I also assign vector erase result to the iterator <code>jet_it = jets.erase(jet_it);<span>. As pointed out by <a href="/@irelandscape">@irelandscape, not doing so is undefined behavior. <p dir="auto">Results after the correction: <pre><code>=> progress: [===> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 6 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 5 => progress: [======> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 6 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 6 => progress: [==========> ] Initial electrons count: 2, initial muons count: 1, initial jets count: 9 Baseline electrons count: 2, baseline muons count: 1, baseline jets count: 9 After overlap removal: Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 9 signal electrons count: 1, signal muons count: 0, signal jets count: 8 => progress: [=============> ] Initial electrons count: 0, initial muons count: 1, initial jets count: 10 Baseline electrons count: 0, baseline muons count: 1, baseline jets count: 10 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 10 signal electrons count: 0, signal muons count: 0, signal jets count: 8 => progress: [=================> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 8 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 8 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 8 signal electrons count: 0, signal muons count: 0, signal jets count: 7 => progress: [====================> ] Initial electrons count: 0, initial muons count: 1, initial jets count: 6 Baseline electrons count: 0, baseline muons count: 1, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 6 => progress: [========================> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 4 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 4 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 4 signal electrons count: 0, signal muons count: 0, signal jets count: 4 => progress: [===========================> ] Initial electrons count: 1, initial muons count: 0, initial jets count: 6 Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 6 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 6 signal electrons count: 0, signal muons count: 0, signal jets count: 6 => progress: [===============================> ] Initial electrons count: 0, initial muons count: 0, initial jets count: 10 Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 10 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 10 signal electrons count: 0, signal muons count: 0, signal jets count: 8 => progress: [==================================>] Initial electrons count: 1, initial muons count: 0, initial jets count: 7 Baseline electrons count: 1, baseline muons count: 0, baseline jets count: 7 After overlap removal: Baseline electrons count: 0, baseline muons count: 0, baseline jets count: 7 signal electrons count: 0, signal muons count: 0, signal jets count: 7 => progress: [===================================]