Commit 21d80fe99d5bd4cec3d27d0a59a2b702c5484e2c

Authored by Marcel Hohmann
1 parent 30e8cb57

Update to allow for multiple ParticleLists

DualTrackTagger/SConscript
1 1 Import('env')
2 2  
3   -env['LIBS'] = ['framework', 'analysis_dataobjects', 'analysis']
  3 +env['LIBS'] = ['framework', 'analysis_dataobjects', 'analysis','$ROOT_LIBS']
4 4  
5 5 Return('env')
... ...
DualTrackTagger/include/DualTrackTaggerModule.h
... ... @@ -3,23 +3,31 @@
3 3  
4 4 // Framework
5 5 #include <framework/datastore/StoreArray.h>
6   -#include <framework/core/Module.h>
  6 +#include <framework/datastore/StoreObjPtr.h>
  7 +//#include <framework/core/Module.h>
  8 +#include <framework/core/HistoModule.h>
7 9  
8 10 // Dataobjects
9 11 #include <analysis/dataobjects/Particle.h>
  12 +#include <analysis/dataobjects/ParticleList.h>
10 13 #include <analysis/dataobjects/EventExtraInfo.h>
11 14  
12 15 // Variables
13 16 #include <analysis/VariableManager/Variables.h>
  17 +#include <analysis/VariableManager/TrackVariables.h>
14 18  
15 19 // Std C++ includes
16 20 #include <iostream>
17 21 #include <vector>
  22 +#include <string>
18 23  
19 24 // Root includes
20 25 #include "TLorentzVector.h"
21 26 #include "TVector3.h"
22 27 #include "TMath.h"
  28 +#include "TH1F.h"
  29 +#include "TString.h"
  30 +
23 31  
24 32  
25 33  
... ... @@ -37,6 +45,8 @@ namespace Belle2 {
37 45  
38 46 virtual ~DualTrackTaggerModule();
39 47  
  48 + virtual void defineHisto(); //Define the histograms
  49 +
40 50 virtual void initialize() override;
41 51  
42 52 virtual void beginRun() override;
... ... @@ -52,7 +62,7 @@ namespace Belle2 {
52 62 private:
53 63 /** Loops through all particles in event adding those below a Pt threshold
54 64 to \code m_candidateParticles \endcode */
55   - virtual void identifyPotentialDuplicateTracks();
  65 + virtual void identifyPotentialDuplicateTracks(Particle* p);
56 66  
57 67 /** Compares two particles as per the cuts outlined in BN#1079 sec 4.1.
58 68 Returns true if the particles meet the cuts. False if they do not. */
... ... @@ -68,10 +78,26 @@ namespace Belle2 {
68 78 /** Calculates BN#1079 sec 4.1 chi2 */
69 79 virtual double chi2(Particle* p);
70 80  
71   - std::vector<Particle*> m_candidateParticles;
  81 + /** Used to fill histogram showing # of tracks associated with one MC particle */
  82 + virtual void fill_h_genParticleIndexMulti();
  83 +
  84 + std::vector<Particle*> m_candidateParticles; // particles with low pt
  85 + std::vector<Particle*> m_allParticlesInList; // used for statistics
72 86 StoreArray<Particle> m_particles;
73 87  
  88 + // statistics - currently unused.
  89 + long int total_particles;
  90 + long int lowPt_particles;
  91 + long int nParticlesVar;
  92 + int total_events;
  93 +
  94 + //params
74 95 bool m_BelleFlag;
  96 + bool m_MCFlag;
  97 + std::vector<std::string> m_strParticleLists;
  98 +
  99 + //histograms
  100 + TH1F *h_genParticleIndexMulti;
75 101 };
76 102 }
77 103 #endif
... ...
DualTrackTagger/src/DualTrackTaggerModule.cc
... ... @@ -21,7 +21,9 @@ REG_MODULE(DualTrackTagger)
21 21 DualTrackTaggerModule::DualTrackTaggerModule() : Module()
22 22 {
23 23 setDescription("Dual Track Tagger Module");
24   - addParam("BelleFlag",m_BelleFlag,"Used to differentiate Belle (b2bii) data from Belle 2");
  24 + addParam("BelleFlag",m_BelleFlag,"Used to differentiate Belle (b2bii) data from Belle 2", false);
  25 + addParam("MCFlag",m_MCFlag,"True = Data is MC, False = Experimental Data. If true outputs additional statistics based on truth information.", false);
  26 + addParam("particleLists", m_strParticleLists, "List of ParticleLists", std::vector<std::string>() );
25 27  
26 28 }
27 29  
... ... @@ -30,21 +32,56 @@ DualTrackTaggerModule::~DualTrackTaggerModule()
30 32 {
31 33 }
32 34  
  35 +void DualTrackTaggerModule::defineHisto()
  36 +{
  37 + h_genParticleIndexMulti = new TH1F("h_genParticleIndexMulti","h_genParticleIndexMulti; # of particles with same genParticleIndex",20,0,20);
  38 + //h_SameChargeMultAngle = new TH1F("h_SameChargeMultAngle", "h_SameChargeMultAngle; #Angle between two particles - both pos/neg",100,0,TMath::Pi());
  39 + //h_OppoChargeMultAngle = new TH1F("h_OppoChargeMultAngle", "h_OppoChargeMultAngle; #Angle between two particles - one pos, one neg",100,0,TMath::Pi());
  40 + //h_pt_dual = new TH1F
  41 + //h_pt_sig = new TH1F
  42 +}
33 43  
34 44 void DualTrackTaggerModule::initialize()
35 45 {
36 46 m_particles.isRequired();
  47 + total_particles = 0;
  48 + lowPt_particles = 0;
  49 + total_events = 0;
  50 +
37 51 }
38 52  
39 53 void DualTrackTaggerModule::beginRun()
40 54 {
  55 + defineHisto();
41 56 }
42 57  
43 58 void DualTrackTaggerModule::event()
44 59 {
45   - identifyPotentialDuplicateTracks();
46   - loopTracks();
47   -}
  60 + for (auto& iList : m_strParticleLists) {
  61 + // reset for new particle list
  62 + m_candidateParticles.clear();
  63 + m_allParticlesInList.clear();
  64 +
  65 + StoreObjPtr<ParticleList> particlelist(iList);
  66 +
  67 + if (!particlelist) {
  68 + B2ERROR("ParticleList " << iList << " not found");
  69 + continue;
  70 + } else {
  71 + if (particlelist->getListSize() == 0) {
  72 + continue;
  73 + }
  74 + for (unsigned int i = 0; i < particlelist->getListSize(); ++i) {
  75 + Particle* iParticle = particlelist->getParticle(i);
  76 + m_allParticlesInList.push_back(iParticle);
  77 + identifyPotentialDuplicateTracks(iParticle);
  78 + } // loop through particles
  79 + loopTracks();
  80 + fill_h_genParticleIndexMulti();
  81 + } // particle list exists
  82 + } // loop list of particle lists
  83 +
  84 +} // DualTrackTaggerModule::event()
48 85  
49 86 void DualTrackTaggerModule::endRun()
50 87 {
... ... @@ -52,13 +89,19 @@ void DualTrackTaggerModule::endRun()
52 89  
53 90 void DualTrackTaggerModule::terminate()
54 91 {
  92 + //std::cout << "FINAL STATISTICS:" << total_particles << " | " << lowPt_particles << std::endl;
  93 + //std::cout << "FINAL STATISTICS:" << total_events << std::endl;
  94 + //std::cout << "FINAL STATISTICS:" << (float)total_particles / total_events << " | " << (float)lowPt_particles / total_events << std::endl;
55 95 }
56 96  
57   -void DualTrackTaggerModule::identifyPotentialDuplicateTracks()
  97 +void DualTrackTaggerModule::identifyPotentialDuplicateTracks(Particle* p)
58 98 {
59   - // clear candidate particle vector from last event
60   - m_candidateParticles.clear();
61   -
  99 + //Initiate all particles as not being dual tracks
  100 + if (!p->hasExtraInfo("DualTrack")) {
  101 + p -> addExtraInfo("DualTrack",0);
  102 + } else {
  103 + p -> setExtraInfo("DualTrack",0); // This should not be needed
  104 + }
62 105  
63 106 // select cut for Belle or Belle 2
64 107 float ptCut;
... ... @@ -69,28 +112,34 @@ void DualTrackTaggerModule::identifyPotentialDuplicateTracks()
69 112 }
70 113  
71 114 // add particles that pass selection to vector.
72   - const unsigned int nParticles = m_particles.getEntries();
73   -
74   - for (unsigned int iParticle = 0; iParticle < nParticles; iParticle++) {
75   - Particle* p = m_particles[iParticle];
76   -
77   - if (!p->hasExtraInfo("DualTrack")) {
78   - p -> addExtraInfo("DualTrack",0);
79   - } else {
80   - p -> setExtraInfo("DualTrack",0); // This should not be needed
81   - }
82 115  
83   - if (p->get4Vector().Perp() < ptCut) {
84   - m_candidateParticles.push_back(p);
85   - } // ptCut
86   - } //iParticles
  116 + //remove particles with no tracking hits
  117 + if (Variable::trackNCDCHits(p) == 0 && Variable::trackNVXDHits(p) == 0) {
  118 + std::cout << "Particle has no tracking hits!" << std::endl;
  119 + return;
  120 + } // no tracking
  121 + //remove neutral particles - should already be done by the tracker cut
  122 + if (p->getCharge() == 0) {
  123 + return;
  124 + } // neutral
  125 + //remove particles that dont meet ptCut
  126 + if (Variable::particlePt(p) >= ptCut) {
  127 + return;
  128 + } // ptCut
  129 +
  130 + m_candidateParticles.push_back(p);
87 131 }
88 132  
89 133 void DualTrackTaggerModule::loopTracks()
90 134 {
91 135 for (Particle* p1 : m_candidateParticles) {
92 136 for (Particle* p2 : m_candidateParticles) {
93   - if (p1 == p2) continue;
  137 + if (p1 == p2) continue; //dont compare particle against itself.
  138 +
  139 + if (Variable::genParticleIndex(p1) == Variable::genParticleIndex(p2)) {
  140 + //std::cout << "Both particles have same MC particle " << Variable::genParticleIndex(p1) << " " << Variable::genParticleIndex(p2) << std::endl;
  141 + }
  142 +
94 143 bool duplicates = compareTracks(p1, p2);
95 144 if (duplicates) {
96 145 tagDuplicate(p1 , p2);
... ... @@ -111,7 +160,7 @@ bool DualTrackTaggerModule::compareTracks(Particle* p1, Particle* p2)
111 160 //angle cut
112 161 float chargeMultiple = p1->getCharge()*p2->getCharge();
113 162  
114   - // same charge or neutral
  163 + // same charge (or neutral but neutrals should be removed already)
115 164 if (chargeMultiple >= 0) {
116 165 if (p1->getMomentum().Angle(p2->getMomentum()) < 15*TMath::Pi()/180) {
117 166 return true;
... ... @@ -135,6 +184,16 @@ bool DualTrackTaggerModule::compareTracks(Particle* p1, Particle* p2)
135 184  
136 185 void DualTrackTaggerModule::tagDuplicate(Particle* p1, Particle* p2)
137 186 {
  187 + //prevent double tagging
  188 + /*
  189 + std:: cout << p1 -> getPDGCode() << " " << p2 -> getPDGCode() << std::endl;
  190 +
  191 + if (p1 -> getExtraInfo("DualTrack") == 1 || p2 -> getExtraInfo("DualTrack") == 1) {
  192 + std:: cout << p1 -> getExtraInfo("DualTrack") << " " << p2 -> getExtraInfo("DualTrack") << std::endl;
  193 + return;
  194 + }
  195 + */
  196 +
138 197  
139 198 double p1chi2 = chi2(p1);
140 199 double p2chi2 = chi2(p2);
... ... @@ -151,3 +210,40 @@ double DualTrackTaggerModule::chi2(Particle* p)
151 210 float gamma = 5;
152 211 return TMath::Power(gamma*Variable::particleDRho(p),2) + TMath::Power(Variable::particleDZ(p),2);
153 212 }
  213 +
  214 +void DualTrackTaggerModule::fill_h_genParticleIndexMulti()
  215 +{
  216 + std::vector<std::pair<int,int>> v_MCIndex;
  217 +
  218 + const unsigned int nParticles = m_allParticlesInList.size();
  219 +
  220 + for (unsigned int iParticle = 0; iParticle < nParticles; iParticle++) {
  221 + Particle* p = m_allParticlesInList[iParticle];
  222 +
  223 +
  224 + bool MCIndexExists = false;
  225 +
  226 + for (unsigned int iPair = 0; iPair < v_MCIndex.size(); iPair++) {
  227 + if (Variable::genParticleIndex(p) == v_MCIndex[iPair].first) {
  228 + MCIndexExists = true;
  229 + v_MCIndex[iPair].second++;
  230 + }
  231 + } // loop through vec
  232 +
  233 +
  234 +
  235 + if (!MCIndexExists) {
  236 + std::pair<int,int> temp_pair = {Variable::genParticleIndex(p), 1};
  237 + v_MCIndex.push_back(temp_pair);
  238 + }
  239 +
  240 + } // loop through StoreArray Particles
  241 +
  242 + for (auto pair : v_MCIndex) {
  243 + h_genParticleIndexMulti -> Fill(pair.second);
  244 + }
  245 +}
  246 +
  247 +//void DualTrackTaggerModule::GroupParticles()
  248 +//{
  249 +//}
... ...