Marcel Hohmann / DualTrackTagger

Commit c5732764d0f543773a6f3e94d35d94fe46639a1d

Authored by Marcel Hohmann
1 parent 21d80fe9

small changes

Showing 2 changed files with 227 additions and 215 deletions   Show diff stats
Inline Side-by-side
DualTrackTagger/include/DualTrackTaggerModule.h
  Show/Hide comments   View file @c573276
... ... @@ -98,6 +98,8 @@ namespace Belle2 {
98 98  
99 99 //histograms
100 100 TH1F *h_genParticleIndexMulti;
  101 + TH1F *h_taggedDuplicate;
  102 + TH1F *h_taggedSingle;
101 103 };
102 104 }
103 105 #endif
... ...
DualTrackTagger/src/DualTrackTaggerModule.cc
  Show/Hide comments   View file @c573276
1   -// Own Include
2   -#include <analysis/modules/DualTrackTagger/DualTrackTaggerModule.h>
3   -
4   -
5   -// -- Info --
6   -// Belle Cuts taken from BN#1079
7   -
8   -// ---- TODO ----
9   -/*
10   - * Update Belle Cuts (maybe train a BDT?)
11   - * Add Belle 2 Cuts (retrain same BDT as for Belle?)
12   - * Remove found duplicate pairs from candidate lists to avoid double tagging and
13   - increase speed
14   - * Find new way to tag duplicates to avoid having to add extraInfo to all particles.
15   -*/
16   -
17   -using namespace Belle2;
18   -
19   -REG_MODULE(DualTrackTagger)
20   -
21   -DualTrackTaggerModule::DualTrackTaggerModule() : Module()
22   -{
23   - setDescription("Dual Track Tagger Module");
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>() );
27   -
28   -}
29   -
30   -
31   -DualTrackTaggerModule::~DualTrackTaggerModule()
32   -{
33   -}
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   -}
43   -
44   -void DualTrackTaggerModule::initialize()
45   -{
46   - m_particles.isRequired();
47   - total_particles = 0;
48   - lowPt_particles = 0;
49   - total_events = 0;
50   -
51   -}
52   -
53   -void DualTrackTaggerModule::beginRun()
54   -{
55   - defineHisto();
56   -}
57   -
58   -void DualTrackTaggerModule::event()
59   -{
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) {
  1 + // Own Include
  2 + #include <analysis/modules/DualTrackTagger/DualTrackTaggerModule.h>
  3 +
  4 +
  5 + // -- Info --
  6 + // Belle Cuts taken from BN#1079
  7 +
  8 + // ---- TODO ----
  9 + /*
  10 + * Update Belle Cuts (maybe train a BDT?)
  11 + * Add Belle 2 Cuts (retrain same BDT as for Belle?)
  12 + * Remove found duplicate pairs from candidate lists to avoid double tagging and
  13 + increase speed
  14 + * Find new way to tag duplicates to avoid having to add extraInfo to all particles.
  15 + */
  16 +
  17 + using namespace Belle2;
  18 +
  19 + REG_MODULE(DualTrackTagger)
  20 +
  21 + DualTrackTaggerModule::DualTrackTaggerModule() : Module()
  22 + {
  23 + setDescription("Dual Track Tagger Module");
  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>() );
  27 +
  28 + }
  29 +
  30 +
  31 + DualTrackTaggerModule::~DualTrackTaggerModule()
  32 + {
  33 + }
  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_taggedDuplicate = new TH1F("h_taggedDuplicate","h_taggedDuplicate; 0 = Singles tagged as Duplicate, 1 = Duplicates tagged correctly",2,0,2);
  41 + h_taggedSingle = new TH1F("h_taggedSingle","h_taggedSingle; 0 = Duplicates tagged as Single, 1 = Singles tagged correctly",2,0,2);
  42 + //h_pt_dual = new TH1F
  43 + //h_pt_sig = new TH1F
  44 + }
  45 +
  46 + void DualTrackTaggerModule::initialize()
  47 + {
  48 + m_particles.isRequired();
  49 + total_particles = 0;
  50 + lowPt_particles = 0;
  51 + total_events = 0;
  52 +
  53 + }
  54 +
  55 + void DualTrackTaggerModule::beginRun()
  56 + {
  57 + defineHisto();
  58 + }
  59 +
  60 + void DualTrackTaggerModule::event()
  61 + {
  62 + for (auto& iList : m_strParticleLists) {
  63 + // reset for new particle list
  64 + m_candidateParticles.clear();
  65 + m_allParticlesInList.clear();
  66 +
  67 + StoreObjPtr<ParticleList> particlelist(iList);
  68 +
  69 + if (!particlelist) {
  70 + B2ERROR("ParticleList " << iList << " not found");
72 71 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()
85   -
86   -void DualTrackTaggerModule::endRun()
87   -{
88   -}
89   -
90   -void DualTrackTaggerModule::terminate()
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;
95   -}
96   -
97   -void DualTrackTaggerModule::identifyPotentialDuplicateTracks(Particle* p)
98   -{
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
  72 + } else {
  73 + if (particlelist->getListSize() == 0) {
  74 + continue;
  75 + }
  76 + for (unsigned int i = 0; i < particlelist->getListSize(); ++i) {
  77 + Particle* iParticle = particlelist->getParticle(i);
  78 + m_allParticlesInList.push_back(iParticle);
  79 + identifyPotentialDuplicateTracks(iParticle);
  80 + } // loop through particles
  81 + loopTracks();
  82 + fill_h_genParticleIndexMulti();
  83 + } // particle list exists
  84 + } // loop list of particle lists
  85 +
  86 + } // DualTrackTaggerModule::event()
  87 +
  88 + void DualTrackTaggerModule::endRun()
  89 + {
104 90 }
105 91  
106   - // select cut for Belle or Belle 2
107   - float ptCut;
108   - if (m_BelleFlag) {
109   - ptCut = 0.275;
110   - } else {
111   - ptCut = 0.275;
  92 + void DualTrackTaggerModule::terminate()
  93 + {
  94 + //std::cout << "FINAL STATISTICS:" << total_particles << " | " << lowPt_particles << std::endl;
  95 + //std::cout << "FINAL STATISTICS:" << total_events << std::endl;
  96 + //std::cout << "FINAL STATISTICS:" << (float)total_particles / total_events << " | " << (float)lowPt_particles / total_events << std::endl;
112 97 }
113 98  
114   - // add particles that pass selection to vector.
115   -
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);
131   -}
132   -
133   -void DualTrackTaggerModule::loopTracks()
134   -{
135   - for (Particle* p1 : m_candidateParticles) {
136   - for (Particle* p2 : m_candidateParticles) {
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   - }
  99 + void DualTrackTaggerModule::identifyPotentialDuplicateTracks(Particle* p)
  100 + {
  101 + //Initiate all particles as not being dual tracks
  102 + if (!p->hasExtraInfo("DualTrack")) {
  103 + p -> addExtraInfo("DualTrack",0);
  104 + } else {
  105 + p -> setExtraInfo("DualTrack",0); // This should not be needed
  106 + }
142 107  
143   - bool duplicates = compareTracks(p1, p2);
144   - if (duplicates) {
145   - tagDuplicate(p1 , p2);
146   - } // duplicats
147   - } //p2
148   - } //p1
149   -}
150   -
151   -bool DualTrackTaggerModule::compareTracks(Particle* p1, Particle* p2)
152   -{
153   - if (m_BelleFlag) {
154   - // momentum cut
155   - TVector3 deltaP = p1->getMomentum() - p2->getMomentum();
156   - if (deltaP.Mag() >= 0.1) {
157   - return false;
  108 + // select cut for Belle or Belle 2
  109 + float ptCut;
  110 + if (m_BelleFlag) {
  111 + ptCut = 0.275;
  112 + } else {
  113 + ptCut = 0.275;
158 114 }
159 115  
160   - //angle cut
161   - float chargeMultiple = p1->getCharge()*p2->getCharge();
  116 + // add particles that pass selection to vector.
  117 +
  118 + //remove particles with no tracking hits
  119 + if (Variable::trackNCDCHits(p) == 0 && Variable::trackNVXDHits(p) == 0) {
  120 + return;
  121 + } // no tracking
  122 + //remove neutral particles - should already be done by the tracker cut
  123 + if (p->getCharge() == 0) {
  124 + return;
  125 + } // neutral
  126 + //remove particles that dont meet ptCut
  127 + if (Variable::particlePt(p) >= ptCut) {
  128 + return;
  129 + } // ptCut
  130 +
  131 + m_candidateParticles.push_back(p);
  132 + }
  133 +
  134 + void DualTrackTaggerModule::loopTracks()
  135 + {
  136 + for (Particle* p1 : m_candidateParticles) {
  137 + for (Particle* p2 : m_candidateParticles) {
  138 + if (p1 == p2) continue; //dont compare particle against itself.
  139 +
  140 +
  141 +
  142 + bool duplicates = compareTracks(p1, p2);
  143 + if (duplicates) {
  144 + tagDuplicate(p1 , p2);
  145 + } // duplicats
  146 +
  147 + // Fill histograms
  148 + if (Variable::genParticleIndex(p1) == Variable::genParticleIndex(p2)) {
  149 + if (duplicates) {
  150 + h_taggedDuplicate -> Fill(1);
  151 + } else {
  152 + h_taggedDuplicate -> Fill(0);
  153 + }
  154 + } else {
  155 + if (duplicates) {
  156 + h_taggedSingle -> Fill(0);
  157 + } else {
  158 + h_taggedSingle -> Fill(1);
  159 + }
  160 + } // different genParticleIndex
  161 +
  162 + } //p2
  163 + } //p1
  164 + }
162 165  
163   - // same charge (or neutral but neutrals should be removed already)
164   - if (chargeMultiple >= 0) {
165   - if (p1->getMomentum().Angle(p2->getMomentum()) < 15*TMath::Pi()/180) {
166   - return true;
  166 + bool DualTrackTaggerModule::compareTracks(Particle* p1, Particle* p2)
  167 + {
  168 + if (m_BelleFlag) {
  169 + // momentum cut
  170 + TVector3 deltaP = p1->getMomentum() - p2->getMomentum();
  171 + if (deltaP.Mag() >= 0.1) {
  172 + return false;
167 173 }
168   - return false;
169 174  
170   - //opposite charge
171   - } else if (chargeMultiple < 0) {
172   - if (p1->getMomentum().Angle(p2->getMomentum()) > 165*TMath::Pi()/180) {
173   - return true;
  175 + //angle cut
  176 + float chargeMultiple = p1->getCharge()*p2->getCharge();
  177 +
  178 + // same charge (or neutral but neutrals should be removed already)
  179 + if (chargeMultiple >= 0) {
  180 + if (p1->getMomentum().Angle(p2->getMomentum()) < 15*TMath::Pi()/180) {
  181 + return true;
  182 + }
  183 + return false;
  184 +
  185 + //opposite charge
  186 + } else if (chargeMultiple < 0) {
  187 + if (p1->getMomentum().Angle(p2->getMomentum()) > 165*TMath::Pi()/180) {
  188 + return true;
  189 + }
  190 + return false;
174 191 }
175 192 return false;
176   - }
177   - return false;
178 193  
179   - } else { // BELLE 2
  194 + } else { // BELLE 2
  195 + return false;
  196 + } // if else
180 197 return false;
181   - } // if else
182   - return false;
183   -}
184   -
185   -void DualTrackTaggerModule::tagDuplicate(Particle* p1, Particle* p2)
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 198 }
195   - */
196 199  
  200 + void DualTrackTaggerModule::tagDuplicate(Particle* p1, Particle* p2)
  201 + {
  202 + //prevent double tagging ---------- Allow double tagging for cases where there are multiple (>2) tracks. At least 1 will always be tagged as a normal track.
  203 + //if (p1 -> getExtraInfo("DualTrack") == 1 || p2 -> getExtraInfo("DualTrack") == 1) {
  204 + // return;
  205 + //}
197 206  
198   - double p1chi2 = chi2(p1);
199   - double p2chi2 = chi2(p2);
200 207  
201   - if (p1chi2 > p2chi2) {
202   - p1 -> setExtraInfo("DualTrack",1);
203   - } else {
204   - p2 -> setExtraInfo("DualTrack",1);
  208 + double p1chi2 = chi2(p1);
  209 + double p2chi2 = chi2(p2);
  210 +
  211 + if (p1chi2 > p2chi2) {
  212 + p1 -> setExtraInfo("DualTrack",1);
  213 + } else {
  214 + p2 -> setExtraInfo("DualTrack",1);
  215 + }
205 216 }
206   -}
207 217  
208   -double DualTrackTaggerModule::chi2(Particle* p)
209   -{
210   - float gamma = 5;
211   - return TMath::Power(gamma*Variable::particleDRho(p),2) + TMath::Power(Variable::particleDZ(p),2);
212   -}
  218 + double DualTrackTaggerModule::chi2(Particle* p)
  219 + {
  220 + float gamma = 5;
  221 + return TMath::Power(gamma*Variable::particleDRho(p),2) + TMath::Power(Variable::particleDZ(p),2);
  222 + }
213 223  
214   -void DualTrackTaggerModule::fill_h_genParticleIndexMulti()
215   -{
216   - std::vector<std::pair<int,int>> v_MCIndex;
  224 + void DualTrackTaggerModule::fill_h_genParticleIndexMulti()
  225 + {
  226 + std::vector<std::pair<int,int>> v_MCIndex;
217 227  
218   - const unsigned int nParticles = m_allParticlesInList.size();
  228 + const unsigned int nParticles = m_allParticlesInList.size();
219 229  
220   - for (unsigned int iParticle = 0; iParticle < nParticles; iParticle++) {
221   - Particle* p = m_allParticlesInList[iParticle];
  230 + for (unsigned int iParticle = 0; iParticle < nParticles; iParticle++) {
  231 + Particle* p = m_allParticlesInList[iParticle];
222 232  
223 233  
224   - bool MCIndexExists = false;
  234 + bool MCIndexExists = false;
225 235  
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
  236 + for (unsigned int iPair = 0; iPair < v_MCIndex.size(); iPair++) {
  237 + if (Variable::genParticleIndex(p) == v_MCIndex[iPair].first) {
  238 + MCIndexExists = true;
  239 + v_MCIndex[iPair].second++;
  240 + }
  241 + } // loop through vec
232 242  
233 243  
234 244  
235   - if (!MCIndexExists) {
236   - std::pair<int,int> temp_pair = {Variable::genParticleIndex(p), 1};
237   - v_MCIndex.push_back(temp_pair);
238   - }
  245 + if (!MCIndexExists) {
  246 + std::pair<int,int> temp_pair = {Variable::genParticleIndex(p), 1};
  247 + v_MCIndex.push_back(temp_pair);
  248 + }
239 249  
240   - } // loop through StoreArray Particles
  250 + } // loop through StoreArray Particles
241 251  
242   - for (auto pair : v_MCIndex) {
243   - h_genParticleIndexMulti -> Fill(pair.second);
  252 + for (auto pair : v_MCIndex) {
  253 + h_genParticleIndexMulti -> Fill(pair.second);
  254 + }
244 255 }
245   -}
246 256  
247   -//void DualTrackTaggerModule::GroupParticles()
248   -//{
249   -//}
  257 + //void DualTrackTaggerModule::GroupParticles()
  258 + //{
  259 + //}
... ...