Bug Summary

File:libraries/ANALYSIS/DParticleComboBlueprint_factory.cc
Location:line 202, column 5
Description:Value stored to 'locSkimMissingFlag' is never read

Annotated Source Code

1#ifdef VTRACE
2#include "vt_user.h"
3#endif
4
5#include "DParticleComboBlueprint_factory.h"
6
7//------------------
8// init
9//------------------
10jerror_t DParticleComboBlueprint_factory::init(void)
11{
12 MAX_DParticleComboBlueprintStepPoolSize = 100;
13
14 dDebugLevel = 0;
15
16 //BEWARE: IF THIS IS CHANGED, CHANGE IN THE ANALYSIS UTILITIES AND THE EVENT WRITER ALSO!!
17 dShowerSelectionTag = "PreSelect";
18 dTrackSelectionTag = "PreSelect";
19
20 dMinProtonMomentum = pair<bool, double>(false, -1.0);
21
22 dMaxExtraGoodTracks = pair<bool, size_t>(false, 4);
23
24 vector<int> hypotheses;
25 hypotheses.push_back(PiPlus);
26 hypotheses.push_back(KPlus);
27 hypotheses.push_back(Proton);
28 hypotheses.push_back(PiMinus);
29 hypotheses.push_back(KMinus);
30
31 ostringstream locMassStream;
32 for(size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i)
33 {
34 locMassStream << hypotheses[loc_i];
35 if(loc_i != (hypotheses.size() - 1))
36 locMassStream << ",";
37 }
38
39 string HYPOTHESES = locMassStream.str();
40 gPARMS->SetDefaultParameter("TRKFIT:HYPOTHESES", HYPOTHESES);
41
42 dMaxNumNeutralShowers = 20;
43 gPARMS->SetDefaultParameter("COMBO:MAX_NEUTRALS", dMaxNumNeutralShowers);
44
45 // Parse MASS_HYPOTHESES strings to make list of masses to try
46 hypotheses.clear();
47 SplitString(HYPOTHESES, hypotheses, ",");
48 for(size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i)
49 dAvailablePIDs.insert(Particle_t(hypotheses[loc_i]));
50
51 return NOERROR;
52}
53
54//------------------
55// brun
56//------------------
57jerror_t DParticleComboBlueprint_factory::brun(jana::JEventLoop* locEventLoop, int32_t runnumber)
58{
59 //BE CAREFUL: DON'T DO ANYTHING THAT REQUIRES THE brun() METHOD OF THIS FACTORY TO BE CALLED!!!!
60 dTrackTimeBasedFactory_Combo = dynamic_cast<DTrackTimeBased_factory_Combo*>(locEventLoop->GetFactory("DTrackTimeBased", "Combo"));
61
62 gPARMS->SetDefaultParameter("COMBOBLUEPRINTS:DEBUG_LEVEL", dDebugLevel);
63
64 gPARMS->SetDefaultParameter("COMBO:TRACK_SELECT_TAG", dTrackSelectionTag);
65 gPARMS->SetDefaultParameter("COMBO:SHOWER_SELECT_TAG", dShowerSelectionTag);
66
67 if(gPARMS->Exists("COMBO:MIN_PROTON_MOMENTUM"))
68 {
69 dMinProtonMomentum.first = true;
70 gPARMS->GetParameter("COMBO:MIN_PROTON_MOMENTUM", dMinProtonMomentum.second);
71 }
72
73 if(gPARMS->Exists("COMBO:MAX_EXTRA_GOOD_TRACKS"))
74 {
75 dMaxExtraGoodTracks.first = true;
76 gPARMS->GetParameter("COMBO:MAX_EXTRA_GOOD_TRACKS", dMaxExtraGoodTracks.second);
77 }
78
79 vector<const DReaction*> locReactions;
80 Get_Reactions(locEventLoop, locReactions);
81 MAX_DParticleComboBlueprintStepPoolSize = 100*locReactions.size();
82
83 return NOERROR;
84}
85
86void DParticleComboBlueprint_factory::Get_Reactions(JEventLoop *locEventLoop, vector<const DReaction*>& locReactions) const
87{
88 // Get list of factories and find all the ones producing
89 // DReaction objects. (A simpler way to do this would be to
90 // just use locEventLoop->Get(...), but then only one plugin could
91 // be used at a time.)
92 locReactions.clear();
93 vector<JFactory_base*> locFactories = locEventLoop->GetFactories();
94 for(size_t loc_i = 0; loc_i < locFactories.size(); ++loc_i)
95 {
96 JFactory<DReaction>* locFactory = dynamic_cast<JFactory<DReaction>* >(locFactories[loc_i]);
97 if(locFactory == NULL__null)
98 continue;
99 if(string(locFactory->Tag()) == "Thrown")
100 continue;
101 // Found a factory producing DReactions. The reaction objects are
102 // produced at the init stage and are persistent through all event
103 // processing so we can grab the list here and append it to our
104 // overall list.
105 vector<const DReaction*> locReactionsSubset;
106 locFactory->Get(locReactionsSubset);
107 locReactions.insert(locReactions.end(), locReactionsSubset.begin(), locReactionsSubset.end());
108 }
109}
110
111void DParticleComboBlueprint_factory::Check_ReactionNames(vector<const DReaction*>& locReactions) const
112{
113 set<string> locReactionNames;
114 set<string> locDuplicateReactionNames;
115 for(auto locReaction : locReactions)
116 {
117 string locReactionName = locReaction->Get_ReactionName();
118 if(locReactionNames.find(locReactionName) == locReactionNames.end())
119 locReactionNames.insert(locReactionName);
120 else
121 locDuplicateReactionNames.insert(locReactionName);
122 }
123
124 if(locDuplicateReactionNames.empty())
125 return;
126
127 cout << "ERROR: MULTIPLE DREACTIONS WITH THE SAME NAME(S): " << endl;
128 for(auto locReactionName : locDuplicateReactionNames)
129 cout << locReactionName << ", ";
130 cout << endl;
131 cout << "ABORTING" << endl;
132 abort();
133}
134
135//------------------
136// evnt
137//------------------
138jerror_t DParticleComboBlueprint_factory::evnt(JEventLoop *locEventLoop, uint64_t eventnumber)
139{
140#ifdef VTRACE
141 VT_TRACER("DParticleComboBlueprint_factory::evnt()");
142#endif
143
144 //CHECK TRIGGER TYPE
145 const DTrigger* locTrigger = NULL__null;
146 locEventLoop->GetSingle(locTrigger);
147 if(!locTrigger->Get_IsPhysicsEvent())
148 return NOERROR;
149
150 Reset_Pools();
151
152 dBlueprintStepMap.clear();
153 dSavedBlueprintSteps.clear();
154 dChargedTracks.clear();
155 dNeutralShowers.clear();
156
157 vector<const DReaction*> locReactions;
158 Get_Reactions(locEventLoop, locReactions);
159 Check_ReactionNames(locReactions);
160
161 locEventLoop->GetSingle(dVertex);
162 locEventLoop->GetSingle(dDetectorMatches);
163
164 vector<const DESSkimData*> locESSkimDataVector;
165 locEventLoop->Get(locESSkimDataVector);
166 const DESSkimData* locESSkimData = locESSkimDataVector.empty() ? NULL__null : locESSkimDataVector[0];
167
168 locEventLoop->Get(dChargedTracks, dTrackSelectionTag.c_str());
169 locEventLoop->Get(dNeutralShowers, dShowerSelectionTag.c_str());
170
171 if(dNeutralShowers.size() > dMaxNumNeutralShowers)
172 return NOERROR; //don't even try
173
174 //sort charged particles into +/-
175 //Note that a DChargedTrack object can sometimes contain both positively and negatively charged hypotheses simultaneously: sometimes the tracking flips the sign of the track
176 dPositiveChargedTracks.clear();
177 dNegativeChargedTracks.clear();
178 for(size_t loc_i = 0; loc_i < dChargedTracks.size(); ++loc_i)
179 {
180 const DChargedTrack* locChargedTrack = dChargedTracks[loc_i];
181 if(locChargedTrack->Contains_Charge(1))
182 dPositiveChargedTracks.push_back(locChargedTrack);
183 if(locChargedTrack->Contains_Charge(-1))
184 dNegativeChargedTracks.push_back(locChargedTrack);
185 }
186 if(dDebugLevel > 0)
187 cout << "#+, #-, #0 particles = " << dPositiveChargedTracks.size() << ", " << dNegativeChargedTracks.size() << ", " << dNeutralShowers.size() << endl;
188
189 //build the combos for each DReaction, IF the event satisfies the DReaction skim requirements
190 for(size_t loc_i = 0; loc_i < locReactions.size(); ++loc_i)
191 {
192 if(locESSkimData != NULL__null)
193 {
194 string locReactionSkimString = locReactions[loc_i]->Get_EventStoreSkims();
195 vector<string> locReactionSkimVector;
196 SplitString(locReactionSkimString, locReactionSkimVector, ",");
197 bool locSkimMissingFlag = false;
198 for(size_t loc_j = 0; loc_j < locReactionSkimVector.size(); ++loc_j)
199 {
200 if(locESSkimData && locESSkimData->Get_IsEventSkim(locReactionSkimVector[loc_j]))
201 continue; //ok so far
202 locSkimMissingFlag = true;
Value stored to 'locSkimMissingFlag' is never read
203 break;
204 }
205 // if(locSkimMissingFlag)
206 // continue; //no blueprints for this reaction!
207 }
208
209 Build_ParticleComboBlueprints(locReactions[loc_i]);
210 }
211
212 return NOERROR;
213}
214
215void DParticleComboBlueprint_factory::Build_ParticleComboBlueprints(const DReaction* locReaction)
216{
217 if(dDebugLevel > 0)
218 cout << "Reaction Name, # Reaction steps = " << locReaction->Get_ReactionName() << ", " << locReaction->Get_NumReactionSteps() << endl;
219 if(locReaction->Get_NumReactionSteps() == 0)
220 return;
221
222 //make sure not more than one missing particle
223 size_t locNumMissingParticles = 0;
224 for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
225 {
226 if(locReaction->Get_ReactionStep(loc_i)->Get_MissingParticleIndex() != -1)
227 ++locNumMissingParticles;
228 }
229 if(locNumMissingParticles > 1)
230 {
231 cout << "ERROR: Too many missing particles in DReaction. No DParticleComboBlueprints generated." << endl;
232 return;
233 }
234
235 //set up combo loop
236 deque<deque<int> > locResumeAtIndexDeque; //1st index is step, 2nd is particle (the initial particle, then the final particles)
237 deque<deque<int> > locNumPossibilitiesDeque; //1st index is step, 2nd is particle (the initial particle, then the final particles)
238 map<int, int> locInitialParticleStepFromIndexMap; //ints are: step index, production-step index
239 map<pair<int, int>, int> locFinalStateDecayStepIndexMap; //ints are: step index, particle index, decay step index
240 if(!Setup_ComboLoop(locReaction, locResumeAtIndexDeque, locNumPossibilitiesDeque, locInitialParticleStepFromIndexMap, locFinalStateDecayStepIndexMap))
241 {
242 if(dDebugLevel > 0)
243 cout << "not enough detected particles with the correct charges for the event: no combos found." << endl;
244 return;
245 }
246
247 if(dDebugLevel > 10)
248 {
249 cout << "locResumeAtIndexDeque: ";
250 for(size_t loc_i = 0; loc_i < locResumeAtIndexDeque.size(); ++loc_i)
251 {
252 for(size_t loc_j = 0; loc_j < locResumeAtIndexDeque[loc_i].size(); ++loc_j)
253 cout << locResumeAtIndexDeque[loc_i][loc_j] << ", ";
254 cout << ";";
255 }
256 cout << endl;
257
258 cout << "locNumPossibilitiesDeque: ";
259 for(size_t loc_i = 0; loc_i < locNumPossibilitiesDeque.size(); ++loc_i)
260 {
261 for(size_t loc_j = 0; loc_j < locNumPossibilitiesDeque[loc_i].size(); ++loc_j)
262 cout << locNumPossibilitiesDeque[loc_i][loc_j] << ", ";
263 cout << ";";
264 }
265 cout << endl;
266
267 cout << "locInitialParticleStepFromIndexMap: ";
268 map<int, int>::iterator locIterator = locInitialParticleStepFromIndexMap.begin();
269 for(; locIterator != locInitialParticleStepFromIndexMap.end(); ++locIterator)
270 cout << locIterator->first << ", " << locIterator->second << endl;
271
272 cout << "locFinalStateDecayStepIndexMap: ";
273 map<pair<int, int>, int>::iterator locPairIterator = locFinalStateDecayStepIndexMap.begin();
274 for(; locPairIterator != locFinalStateDecayStepIndexMap.end(); ++locPairIterator)
275 cout << locPairIterator->first.first << ", " << locPairIterator->first.second << ": " << locPairIterator->second << endl;
276 }
277
278 //find the combos!!
279 dCurrentComboSourceObjects.clear();
280 Find_Combos(locReaction, locResumeAtIndexDeque, locNumPossibilitiesDeque, locInitialParticleStepFromIndexMap, locFinalStateDecayStepIndexMap);
281
282 if(dDebugLevel > 10)
283 {
284 cout << "print pointers: " << endl;
285 for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
286 {
287 cout << "COMBO " << loc_i << endl;
288 for(size_t loc_j = 0; loc_j < _data[loc_i]->Get_NumParticleComboBlueprintSteps(); ++loc_j)
289 {
290 cout << "Step " << loc_j << " pointers: ";
291 for(size_t loc_k = 0; loc_k < _data[loc_i]->Get_ParticleComboBlueprintStep(loc_j)->Get_NumFinalParticleSourceObjects(); ++loc_k)
292 cout << _data[loc_i]->Get_ParticleComboBlueprintStep(loc_j)->Get_FinalParticle_SourceObject(loc_k) << ", ";
293 cout << endl;
294 }
295 }
296 }
297
298 return;
299}
300
301bool DParticleComboBlueprint_factory::Setup_ComboLoop(const DReaction* locReaction, deque<deque<int> >& locResumeAtIndexDeque, deque<deque<int> >& locNumPossibilitiesDeque, map<int, int>& locInitialParticleStepFromIndexMap, map<pair<int, int>, int>& locFinalStateDecayStepIndexMap)
302{
303 //setup locResumeAtIndexDeque, & locNumPossibilitiesDeque
304 Particle_t locAnalysisPID;
305 int locMissingParticleIndex;
306 unsigned int locNumSteps = locReaction->Get_NumReactionSteps();
307 locResumeAtIndexDeque.clear();
308 locResumeAtIndexDeque.clear();
309 int locCharge;
310 int locNumNeededChargedParticles = 0, locNumNeededPositiveParticles = 0, locNumNeededNegativeParticles = 0, locNumNeededNeutralParticles = 0;
311
312 locInitialParticleStepFromIndexMap[0] = -1;
313 for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)
314 {
315 const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_i);
316 size_t locNumFinalParticles = locReactionStep->Get_NumFinalParticleIDs();
317
318 //setup final state particles num possibilities & resume-at index deque
319 deque<int> locTempDeque(locNumFinalParticles, 0);
320 locResumeAtIndexDeque.push_back(locTempDeque);
321 locMissingParticleIndex = locReactionStep->Get_MissingParticleIndex();
322 for(size_t loc_j = 0; loc_j < locNumFinalParticles; ++loc_j)
323 {
324 locAnalysisPID = locReactionStep->Get_FinalParticleID(loc_j);
325
326 if(locMissingParticleIndex == int(loc_j))
327 {
328 locTempDeque[loc_j] = 1;
329 continue;
330 }
331
332 // check to see if this particle has a decay that is represented in a future step
333 // e.g. on Lambda in g, p -> K+, Lambda; where a later step is Lambda -> p, pi-
334 int locDecayStepIndex = locReaction->Get_DecayStepIndex(loc_i, loc_j);
335 if(locDecayStepIndex >= 0)
336 {
337 if(dDebugLevel > 10)
338 cout << ParticleType(locAnalysisPID) << " decays later in the reaction, at step index " << locDecayStepIndex << endl;
339 locTempDeque[loc_j] = 1;
340 locFinalStateDecayStepIndexMap[pair<int, int>(loc_i, loc_j)] = locDecayStepIndex; //store step where this particle decays
341 locInitialParticleStepFromIndexMap[locDecayStepIndex] = loc_i; //store step where this particle is produced
342 continue;
343 }
344
345 //else use detected particles
346 locCharge = ParticleCharge(locAnalysisPID);
347 if(locCharge > 0)
348 {
349 ++locNumNeededPositiveParticles;
350 ++locNumNeededChargedParticles;
351 locTempDeque[loc_j] = dPositiveChargedTracks.size();
352 }
353 else if(locCharge < 0)
354 {
355 ++locNumNeededNegativeParticles;
356 ++locNumNeededChargedParticles;
357 locTempDeque[loc_j] = dNegativeChargedTracks.size();
358 }
359 else
360 {
361 ++locNumNeededNeutralParticles;
362 locTempDeque[loc_j] = dNeutralShowers.size();
363 }
364 }
365 locNumPossibilitiesDeque.push_back(locTempDeque);
366 }
367
368 //cut if too-few particles
369 if((locNumNeededPositiveParticles > int(dPositiveChargedTracks.size())) || (locNumNeededNegativeParticles > int(dNegativeChargedTracks.size())))
370 return false; //not enough particles of a given charge for the event
371 if((locNumNeededChargedParticles > int(dChargedTracks.size())) || (locNumNeededNeutralParticles > int(dNeutralShowers.size())))
372 return false; //not enough particles of a given charge for the event //#charged can fail here if a charged track has hypotheses with different charges
373
374 //cut if too-many particles
375 pair<bool, double> locMaxExtraGoodTracks = dMaxExtraGoodTracks.first ? dMaxExtraGoodTracks : locReaction->Get_MaxExtraGoodTracks();
376 if(locMaxExtraGoodTracks.first)
377 {
378 size_t locMaxNumTracks = locNumNeededPositiveParticles + locNumNeededNegativeParticles + locMaxExtraGoodTracks.second;
379 if(dChargedTracks.size() > locMaxNumTracks)
380 return false;
381 }
382
383 //make sure decaying particles are valid: one entry per step except the first one
384 for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)
385 {
386 if(locInitialParticleStepFromIndexMap.find(loc_i) == locInitialParticleStepFromIndexMap.end())
387 return false;
388 }
389
390 return true;
391}
392
393void DParticleComboBlueprint_factory::Find_Combos(const DReaction* locReaction, deque<deque<int> >& locResumeAtIndexDeque, const deque<deque<int> >& locNumPossibilitiesDeque, map<int, int>& locInitialParticleStepFromIndexMap, map<pair<int, int>, int>& locFinalStateDecayStepIndexMap)
394{
395#ifdef VTRACE
396 VT_TRACER("DParticleComboBlueprint_factory::Find_Combos()");
397#endif
398 DParticleComboBlueprint* locParticleComboBlueprint = new DParticleComboBlueprint();
399 locParticleComboBlueprint->Set_Reaction(locReaction);
400
401 int locStepIndex = locReaction->Get_NumReactionSteps() - 1;
402 int locParticleIndex = 0; //final = 0 -> (#final - 1)
403 DParticleComboBlueprintStep* locParticleComboBlueprintStep = Get_ParticleComboBlueprintStepResource();
404 locParticleComboBlueprintStep->Set_ReactionStep(locReaction->Get_ReactionStep(locStepIndex));
405 locParticleComboBlueprintStep->Set_InitialParticleDecayFromStepIndex(locInitialParticleStepFromIndexMap[locStepIndex]);
406
407 do
408 {
409 if(dDebugLevel > 10)
410 cout << "do loop: step & particle indices = " << locStepIndex << ", " << locParticleIndex << endl;
411 if(locParticleIndex == int(locNumPossibilitiesDeque[locStepIndex].size()))
412 {
413 if(dDebugLevel > 10)
414 cout << "end of reaction step" << endl;
415 if(!Handle_EndOfReactionStep(locReaction, locParticleComboBlueprint, locParticleComboBlueprintStep, locStepIndex, locParticleIndex, locResumeAtIndexDeque, locNumPossibilitiesDeque, locInitialParticleStepFromIndexMap))
416 break;
417 continue;
418 }
419
420 //get Analysis PID & the resume-at index
421 int locMissingParticleIndex = locReaction->Get_ReactionStep(locStepIndex)->Get_MissingParticleIndex();
422 Particle_t locAnalysisPID = locReaction->Get_ReactionStep(locStepIndex)->Get_FinalParticleID(locParticleIndex);
423 int& locResumeAtIndex = locResumeAtIndexDeque[locStepIndex][locParticleIndex];
424
425 if(dDebugLevel > 10)
426 cout << "do loop: locAnalysisPID, locMissingParticleIndex, locResumeAtIndex = " << ParticleType(locAnalysisPID) << ", " << locMissingParticleIndex << ", " << locResumeAtIndex << endl;
427
428 //handle if this is a missing particle
429 if(locMissingParticleIndex == locParticleIndex)
430 {
431 if(dDebugLevel > 10)
432 cout << "missing particle" << endl;
433 // e.g. on neutron in g, p -> pi+, (n)
434 //only one possibility ("missing"), so just set NULL and advance
435 locParticleComboBlueprintStep->Add_FinalParticle_SourceObject(NULL__null, -1); //missing
436 locResumeAtIndex = 1;
437 ++locParticleIndex;
438 continue;
439 }
440
441 // check to see if this particle has a decay that is represented in a future step
442 // e.g. on Lambda in g, p -> K+, Lambda; where a later step is Lambda -> p, pi-
443 pair<int, int> locParticlePair(locStepIndex, locParticleIndex);
444 if(locFinalStateDecayStepIndexMap.find(locParticlePair) != locFinalStateDecayStepIndexMap.end())
445 {
446 int locDecayStepIndex = locFinalStateDecayStepIndexMap[locParticlePair];
447 if(dDebugLevel > 10)
448 cout << "decaying particle" << endl;
449 locParticleComboBlueprintStep->Add_FinalParticle_SourceObject(NULL__null, locDecayStepIndex); //decaying
450 locResumeAtIndex = 1;
451 ++locParticleIndex;
452 continue;
453 }
454
455 //if two detected particles of same type in a step: locResumeAtIndex must always be >= the previous locResumeAtIndex (prevents duplicates) e.g. g, d -> p, p, pi-
456 //search for same pid previously in this step (and non-missing)
457 for(int loc_i = locParticleIndex - 1; loc_i >= 0; --loc_i)
458 {
459 if(loc_i == locMissingParticleIndex)
460 continue;
461 if(locReaction->Get_ReactionStep(locStepIndex)->Get_FinalParticleID(loc_i) == locAnalysisPID)
462 {
463 if(locResumeAtIndex < locResumeAtIndexDeque[locStepIndex][loc_i])
464 locResumeAtIndex = locResumeAtIndexDeque[locStepIndex][loc_i];
465 if(dDebugLevel > 10)
466 cout << "dupe type in step; locResumeAtIndex = " << locResumeAtIndex << endl;
467 break;
468 }
469 }
470
471 // else grab a detected track
472 const JObject* locSourceObject = Grab_DetectedParticle(locReaction, locAnalysisPID, locResumeAtIndex);
473 if(locSourceObject == NULL__null)
474 {
475 if(dDebugLevel > 10)
476 cout << "can't find detected particle" << endl;
477 if(!Handle_Decursion(locParticleComboBlueprint, locResumeAtIndexDeque, locNumPossibilitiesDeque, locParticleIndex, locStepIndex, locParticleComboBlueprintStep))
478 break;
479 continue;
480 }
481
482 if(dDebugLevel > 10)
483 cout << "detected track found, locResumeAtIndex now = " << locResumeAtIndex << endl;
484 locParticleComboBlueprintStep->Add_FinalParticle_SourceObject(locSourceObject, -2); //detected
485 dCurrentComboSourceObjects.insert(locSourceObject);
486 ++locParticleIndex;
487 }
488 while(true);
489 delete locParticleComboBlueprint; //delete the last, extra one
490}
491
492bool DParticleComboBlueprint_factory::Handle_EndOfReactionStep(const DReaction* locReaction, DParticleComboBlueprint*& locParticleComboBlueprint, DParticleComboBlueprintStep*& locParticleComboBlueprintStep, int& locStepIndex, int& locParticleIndex, deque<deque<int> >& locResumeAtIndexDeque, const deque<deque<int> >& locNumPossibilitiesDeque, map<int, int>& locInitialParticleStepFromIndexMap)
493{
494 //end of step
495
496 //if all of the particle types in this step are identical to all of the particle types in a previously-done step (and none of either are missing), regardless of the order they are listed:
497 //the locResumeAtIndex system of the new step's particles must be >= the locResumeAtIndex system of the particles in the LAST OF the matching previously-done steps (e.g. could be 3x pi0s)
498 //in other words: if 2x pi0 -> g, g; CD, AB ok and AB, CD bad; but also BC, AD ok (BC system occurs after AD system)
499 if(Check_IfDuplicateStepCombo(locParticleComboBlueprint, locParticleComboBlueprintStep, locStepIndex, locResumeAtIndexDeque, locNumPossibilitiesDeque)) //make sure none of the dupe particles are missing
500 {
501 if(dDebugLevel > 10)
502 cout << "duplicate step combo" << endl;
503 if(!Handle_Decursion(locParticleComboBlueprint, locResumeAtIndexDeque, locNumPossibilitiesDeque, locParticleIndex, locStepIndex, locParticleComboBlueprintStep))
504 return false;
505 return true;
506 }
507
508 // cut on invariant mass if desired
509 Particle_t locStepInitialPID = locParticleComboBlueprintStep->Get_InitialParticleID();
510 double locMinInvariantMass = 1.0, locMaxInvariantMass = -1.0;
511 if(locReaction->Get_InvariantMassCut(locStepInitialPID, locMinInvariantMass, locMaxInvariantMass))
512 {
513 DLorentzVector locP4;
514 if(Calc_FinalStateP4(locReaction->Get_NumReactionSteps(), locParticleComboBlueprint, locParticleComboBlueprintStep, -1, locP4))
515 {
516 double locInvariantMass = locP4.M();
517 if((locInvariantMass < locMinInvariantMass) || (locInvariantMass > locMaxInvariantMass))
518 {
519 if(dDebugLevel > 10)
520 cout << "bad invariant mass" << endl;
521 if(!Handle_Decursion(locParticleComboBlueprint, locResumeAtIndexDeque, locNumPossibilitiesDeque, locParticleIndex, locStepIndex, locParticleComboBlueprintStep))
522 return false;
523 return true;
524 }
525 }
526 }
527
528 //step is good: advance to next step
529
530 //first check to see if identical to a previous saved step; if so, just save the old step and recycle the current one
531 map<DParticleComboBlueprintStep, DParticleComboBlueprintStep*>::iterator locStepIterator = dBlueprintStepMap.find(*locParticleComboBlueprintStep);
532 if(locStepIterator != dBlueprintStepMap.end())
533 {
534 //identical step found, recycle current one
535 Recycle_ParticleComboBlueprintStep(locParticleComboBlueprintStep);
536 locParticleComboBlueprintStep = locStepIterator->second;
537 }
538 locParticleComboBlueprint->Prepend_ParticleComboBlueprintStep(locParticleComboBlueprintStep);
539
540 --locStepIndex;
541 if(dDebugLevel > 10)
542 cout << "handle end: new step index, #steps = " << locStepIndex << ", " << locReaction->Get_NumReactionSteps() << endl;
543 if(locStepIndex != -1)
544 {
545 // did not complete the chain yet
546 locParticleIndex = 0;
547
548 locParticleComboBlueprintStep = Get_ParticleComboBlueprintStepResource();
549 locParticleComboBlueprintStep->Set_ReactionStep(locReaction->Get_ReactionStep(locStepIndex));
550 locParticleComboBlueprintStep->Set_InitialParticleDecayFromStepIndex(locInitialParticleStepFromIndexMap[locStepIndex]);
551 return true;
552 }
553
554 if(dDebugLevel > 10)
555 cout << "save combo" << endl;
556 _data.push_back(locParticleComboBlueprint);
557
558 //register steps so they won't accidentally be recycled later, and so that they can be
559 for(size_t loc_i = 0; loc_i < locParticleComboBlueprint->Get_NumParticleComboBlueprintSteps(); ++loc_i)
560 {
561 const DParticleComboBlueprintStep* locParticleComboBlueprintStep = locParticleComboBlueprint->Get_ParticleComboBlueprintStep(loc_i);
562 if(dSavedBlueprintSteps.find(locParticleComboBlueprintStep) != dSavedBlueprintSteps.end())
563 continue;
564 dSavedBlueprintSteps.insert(locParticleComboBlueprintStep);
565 dBlueprintStepMap[*locParticleComboBlueprintStep] = const_cast<DParticleComboBlueprintStep*>(locParticleComboBlueprintStep);
566 }
567
568 locParticleComboBlueprint = new DParticleComboBlueprint(*locParticleComboBlueprint); //clone so don't alter saved object
569 locParticleComboBlueprintStep = NULL__null;
570
571 //if true, once one is found: bail on search
572 if(locReaction->Get_AnyBlueprintFlag())
573 return false; //skip to the end
574
575 if(!Handle_Decursion(locParticleComboBlueprint, locResumeAtIndexDeque, locNumPossibilitiesDeque, locParticleIndex, locStepIndex, locParticleComboBlueprintStep))
576 return false;
577 return true;
578}
579
580bool DParticleComboBlueprint_factory::Handle_Decursion(DParticleComboBlueprint* locParticleComboBlueprint, deque<deque<int> >& locResumeAtIndexDeque, const deque<deque<int> >& locNumPossibilitiesDeque, int& locParticleIndex, int& locStepIndex, DParticleComboBlueprintStep*& locParticleComboBlueprintStep)
581{
582 do
583 {
584 if(dDebugLevel > 50)
585 cout << "decursion: step, particle indices = " << locStepIndex << ", " << locParticleIndex << endl;
586 if(locStepIndex == -1) //just saved a blueprint
587 {
588 if(dDebugLevel > 50)
589 cout << "just saved blueprint" << endl;
590 ++locStepIndex;
591
592 //this looks like it loses the pointer during the pop step, but it doesn't: it was just saved in dSavedBlueprintSteps
593 locParticleComboBlueprintStep = Get_ParticleComboBlueprintStepResource();
594 *locParticleComboBlueprintStep = *(locParticleComboBlueprint->Pop_ParticleComboBlueprintStep());
595 locParticleIndex = locResumeAtIndexDeque[locStepIndex].size() - 1;
596 dCurrentComboSourceObjects.erase(locParticleComboBlueprintStep->Pop_FinalParticle_SourceObject());
597 if(dDebugLevel > 50)
598 cout << "step index, particle index, resume at, #possible = " << locStepIndex << ", " << locParticleIndex << ", " << locResumeAtIndexDeque[locStepIndex][locParticleIndex] << ", " << locNumPossibilitiesDeque[locStepIndex][locParticleIndex] << endl;
599 continue;
600 }
601 else if(locParticleIndex == int(locNumPossibilitiesDeque[locStepIndex].size())) //end of a step: step was found to be a duplicate OR failed a mass cut
602 {
603 if(dDebugLevel > 50)
604 cout << "failed at end of a step" << endl;
605 --locParticleIndex;
606 dCurrentComboSourceObjects.erase(locParticleComboBlueprintStep->Pop_FinalParticle_SourceObject());
607 if(dDebugLevel > 50)
608 cout << "step index, particle index, resume at, #possible = " << locStepIndex << ", " << locParticleIndex << ", " << locResumeAtIndexDeque[locStepIndex][locParticleIndex] << ", " << locNumPossibilitiesDeque[locStepIndex][locParticleIndex] << endl;
609 continue;
610 }
611
612 //locParticleIndex will represent the particle it failed to find a combo for
613 locResumeAtIndexDeque[locStepIndex][locParticleIndex] = 0;
614 --locParticleIndex;
615 if(locParticleIndex >= 0) //else is initial particle (will pop the entire step)
616 dCurrentComboSourceObjects.erase(locParticleComboBlueprintStep->Pop_FinalParticle_SourceObject());
617 else
618 {
619 if(dDebugLevel > 50)
620 cout << "pop the step" << endl;
621 ++locStepIndex;
622 if(locStepIndex >= int(locResumeAtIndexDeque.size()))
623 return false; //end of finding all blueprints
624
625 const DParticleComboBlueprintStep* locPoppedStep = locParticleComboBlueprint->Pop_ParticleComboBlueprintStep();
626 if(dSavedBlueprintSteps.find(locPoppedStep) == dSavedBlueprintSteps.end())
627 {
628 //popped step is not in a saved object: just use it, and recycle the step object we were just editing
629 dParticleComboBlueprintStepPool_Available.push_back(locParticleComboBlueprintStep);
630 locParticleComboBlueprintStep = const_cast<DParticleComboBlueprintStep*>(locPoppedStep);
631 }
632 else //popped step is in a saved object: cannot edit it: use the current object
633 *locParticleComboBlueprintStep = *locPoppedStep;
634
635 locParticleIndex = locResumeAtIndexDeque[locStepIndex].size() - 1;
636 dCurrentComboSourceObjects.erase(locParticleComboBlueprintStep->Pop_FinalParticle_SourceObject());
637 }
638 if(dDebugLevel > 50)
639 cout << "resume at, #possible = " << locResumeAtIndexDeque[locStepIndex][locParticleIndex] << ", " << locNumPossibilitiesDeque[locStepIndex][locParticleIndex] << endl;
640 }
641 while(locResumeAtIndexDeque[locStepIndex][locParticleIndex] == locNumPossibilitiesDeque[locStepIndex][locParticleIndex]);
642 return true;
643}
644
645bool DParticleComboBlueprint_factory::Check_IfDuplicateStepCombo(const DParticleComboBlueprint* locParticleComboBlueprint, const DParticleComboBlueprintStep* locCurrentStep, int locStepIndex, deque<deque<int> >& locResumeAtIndexDeque, const deque<deque<int> >& locNumPossibilitiesDeque) const
646{
647#ifdef VTRACE
648 VT_TRACER("DParticleComboBlueprint_factory::Check_IfDuplicateStepCombo()");
649#endif
650 //note that final particle ids could be rearranged in a different order
651 map<Particle_t, unsigned int> locParticleTypeCount_CurrentStep;
652 bool locIsAParticleDetected = false;
653 for(size_t loc_i = 0; loc_i < locCurrentStep->Get_NumFinalParticleSourceObjects(); ++loc_i)
654 {
655 Particle_t locPID = locCurrentStep->Get_FinalParticleID(loc_i);
656 if(locParticleTypeCount_CurrentStep.find(locPID) == locParticleTypeCount_CurrentStep.end())
657 locParticleTypeCount_CurrentStep[locPID] = 1;
658 else
659 ++(locParticleTypeCount_CurrentStep[locPID]);
660 if(locCurrentStep->Is_FinalParticleDetected(loc_i))
661 locIsAParticleDetected = true;
662 }
663 if(!locIsAParticleDetected)
664 return false; //dupes of this sort only occur when dealing with at least some detected particles
665
666
667 //if all of the particle types in this step are identical to all of the particle types in a previously-done step (and none of either are missing), regardless of the order they are listed:
668 //the locResumeAtIndex system of the new step's particles must be >= the locResumeAtIndex system of the particles in the LAST OF the matching previously-done steps (e.g. could be 3x pi0s)
669 //in other words: if 2x pi0 -> g, g; CD, AB ok and AB, CD bad; but also BC, AD ok (BC system occurs after AD system)
670 //this works regardless of how many particles or what their PIDs are, or if it's a long decay chain or not
671
672 //search future steps for a match (identical particles, none missing, may be in a different order)
673 int locNumFutureSteps = locParticleComboBlueprint->Get_NumParticleComboBlueprintSteps();
674 for(int loc_i = 0; loc_i < locNumFutureSteps; ++loc_i)
675 {
676 int locFutureStepIndex = locResumeAtIndexDeque.size() - locNumFutureSteps + loc_i;
677 const DParticleComboBlueprintStep* locFutureStep = locParticleComboBlueprint->Get_ParticleComboBlueprintStep(loc_i);
678 if(!Check_IfStepsAreIdentical(locParticleComboBlueprint, locCurrentStep, locFutureStep))
679 continue;
680
681 //step is a match
682 if(dDebugLevel > 10)
683 cout << "step at index = " << locStepIndex << " matches a future step at index " << locFutureStepIndex << endl;
684
685 //the resume-at index of the first detected particle listed in the future step MUST be less than that of the first particle listed in the current step that has the same pid
686 //else the current particle combo was already previously used for the future step in a previous combo
687 //there must be at least one detected particle because this was required earlier
688 size_t locFirstDetectedParticleIndex = 0;
689 for(size_t loc_j = 0; loc_j < locFutureStep->Get_NumFinalParticleSourceObjects(); ++loc_j)
690 {
691 if(!locFutureStep->Is_FinalParticleDetected(loc_j))
692 continue;
693 locFirstDetectedParticleIndex = loc_j;
694 break;
695 }
696
697 Particle_t locFutureFirstPID = locFutureStep->Get_FinalParticleID(locFirstDetectedParticleIndex);
698 int locFutureResumeAtIndex = locResumeAtIndexDeque[locFutureStepIndex][locFirstDetectedParticleIndex];
699 for(size_t loc_j = 0; loc_j < locCurrentStep->Get_NumFinalParticleSourceObjects(); ++loc_j)
700 {
701 if(locCurrentStep->Get_FinalParticleID(loc_j) != locFutureFirstPID)
702 continue;
703 int locCurrentResumeAtIndex = locResumeAtIndexDeque[locStepIndex][loc_j];
704 if(dDebugLevel > 10)
705 cout << "future first pid, future resume index, current resume index = " << locFutureFirstPID << ", " << locFutureResumeAtIndex << ", " << locCurrentResumeAtIndex << endl;
706 if(locCurrentResumeAtIndex < locFutureResumeAtIndex)
707 {
708 if(loc_j == (locCurrentStep->Get_NumFinalParticleSourceObjects() - 1))
709 {
710 //on the last particle index: advance resume-at index to the smallest-possible, non-duplicate value to speed up the search
711 locResumeAtIndexDeque[locStepIndex][loc_j] = locFutureResumeAtIndex;
712 if(dDebugLevel > 10)
713 cout << "resume-at index updated to " << locFutureResumeAtIndex << endl;
714 }
715 if(dDebugLevel > 10)
716 cout << "duplicate step, force abort" << endl;
717 return true;
718 }
719
720 if(dDebugLevel > 10)
721 cout << "combo is ok (step is duplicate, but combo is not)" << endl;
722 return false; //else combo is ok (step is duplicate, but combo is not)
723 }
724 }
725
726 if(dDebugLevel > 10)
727 cout << "step at index = " << locStepIndex << " has no similar steps." << endl;
728
729 return false;
730}
731
732bool DParticleComboBlueprint_factory::Check_IfStepsAreIdentical(const DParticleComboBlueprint* locParticleComboBlueprint, const DParticleComboBlueprintStep* locCurrentStep, const DParticleComboBlueprintStep* locPreviousStep) const
733{
734 if(locCurrentStep->Get_MissingParticleIndex() != -1)
735 return false; //missing particle somewhere in the step: cannot be a duplicate (already have checked not more than one missing particle)
736 if(locPreviousStep->Get_MissingParticleIndex() != -1)
737 return false; //missing particle somewhere in the step: cannot be a duplicate (already have checked not more than one missing particle)
738
739 if(locPreviousStep->Get_InitialParticleID() != locCurrentStep->Get_InitialParticleID())
740 return false;
741 if(locPreviousStep->Get_TargetParticleID() != locCurrentStep->Get_TargetParticleID())
742 return false;
743
744 if(locPreviousStep->Get_NumFinalParticleSourceObjects() != locCurrentStep->Get_NumFinalParticleSourceObjects())
745 return false;
746
747 //a step is identical only if the decay chains of it's final state particles are also identical
748 //note that final particle ids could be rearranged in a different order
749 map<Particle_t, unsigned int> locParticleTypeCount_CurrentStep;
750 map<Particle_t, deque<int> > locDecayingParticleStepIndices_CurrentStep;
751 for(size_t loc_i = 0; loc_i < locCurrentStep->Get_NumFinalParticleSourceObjects(); ++loc_i)
752 {
753 Particle_t locPID = locCurrentStep->Get_FinalParticleID(loc_i);
754 if(locParticleTypeCount_CurrentStep.find(locPID) == locParticleTypeCount_CurrentStep.end())
755 locParticleTypeCount_CurrentStep[locPID] = 1;
756 else
757 ++(locParticleTypeCount_CurrentStep[locPID]);
758 if(locCurrentStep->Is_FinalParticleDecaying(loc_i))
759 locDecayingParticleStepIndices_CurrentStep[locPID].push_back(locCurrentStep->Get_DecayStepIndex(loc_i));
760 }
761
762 //compare against the previous step
763 map<Particle_t, deque<int> > locDecayingParticleStepIndices_PreviousStep;
764 for(size_t loc_j = 0; loc_j < locPreviousStep->Get_NumFinalParticleSourceObjects(); ++loc_j)
765 {
766 Particle_t locPID = locPreviousStep->Get_FinalParticleID(loc_j);
767 if(locParticleTypeCount_CurrentStep.find(locPID) == locParticleTypeCount_CurrentStep.end())
768 return false;
769 if(locParticleTypeCount_CurrentStep[locPID] == 1)
770 locParticleTypeCount_CurrentStep.erase(locPID); //if another one, will fail .find() check
771 else
772 --(locParticleTypeCount_CurrentStep[locPID]);
773 if(locPreviousStep->Is_FinalParticleDecaying(loc_j))
774 locDecayingParticleStepIndices_PreviousStep[locPID].push_back(locPreviousStep->Get_DecayStepIndex(loc_j));
775 }
776 if(!locParticleTypeCount_CurrentStep.empty())
777 return false;
778
779 //all of the particle types are the same, now check the decays
780 map<Particle_t, deque<int> >::iterator locIterator;
781 for(locIterator = locDecayingParticleStepIndices_CurrentStep.begin(); locIterator != locDecayingParticleStepIndices_CurrentStep.end(); ++locIterator)
782 {
783 Particle_t locPID = locIterator->first;
784 deque<int>& locDecayStepIndices_Current = locIterator->second;
785 deque<int> locDecayStepIndices_Previous = locDecayingParticleStepIndices_PreviousStep[locPID];
786 deque<int>::iterator locDequeIterator;
787 for(size_t loc_i = 0; loc_i < locDecayStepIndices_Current.size(); ++loc_i)
788 {
789 bool locMatchFoundFlag = false;
790 for(locDequeIterator = locDecayStepIndices_Previous.begin(); locDequeIterator != locDecayStepIndices_Previous.end(); ++locDequeIterator)
791 {
792 const DParticleComboBlueprintStep* locNewCurrentStep = locParticleComboBlueprint->Get_ParticleComboBlueprintStep(locDecayStepIndices_Current[loc_i]);
793 const DParticleComboBlueprintStep* locNewPreviousStep = locParticleComboBlueprint->Get_ParticleComboBlueprintStep(*locDequeIterator);
794 if(!Check_IfStepsAreIdentical(locParticleComboBlueprint, locNewCurrentStep, locNewPreviousStep))
795 continue;
796 locMatchFoundFlag = true;
797 locDecayStepIndices_Previous.erase(locDequeIterator);
798 break;
799 }
800 if(!locMatchFoundFlag)
801 return false; //no matches found
802 }
803 }
804
805 //step is a match
806 return true;
807}
808
809const JObject* DParticleComboBlueprint_factory::Grab_DetectedParticle(const DReaction* locReaction, Particle_t locAnalysisPID, int& locResumeAtIndex)
810{
811 int locAnalysisCharge = ParticleCharge(locAnalysisPID);
812 if(dDebugLevel > 10)
813 cout << "Grab_DetectedTrack: PID, Charge = " << ParticleType(locAnalysisPID) << ", " << locAnalysisCharge << endl;
814 if(locAnalysisCharge == 0)
815 return Grab_NeutralShower(dNeutralShowers, locResumeAtIndex);
816 else if(locAnalysisCharge > 0)
817 return Grab_DetectedTrack(locReaction, locAnalysisPID, dPositiveChargedTracks, locResumeAtIndex);
818 else
819 return Grab_DetectedTrack(locReaction, locAnalysisPID, dNegativeChargedTracks, locResumeAtIndex);
820}
821
822const JObject* DParticleComboBlueprint_factory::Grab_NeutralShower(vector<const DNeutralShower*>& locNeutralShowers, int& locResumeAtIndex)
823{
824 if(dDebugLevel > 10)
825 cout << "Grab_NeutralShower: resume at, #possible = " << locResumeAtIndex << ", " << locNeutralShowers.size() << endl;
826 if(locResumeAtIndex >= int(locNeutralShowers.size()))
827 return NULL__null;
828 do
829 {
830 const JObject* locObject = static_cast<const JObject*>(locNeutralShowers[locResumeAtIndex]);
831 ++locResumeAtIndex;
832
833 //make sure not used currently
834 if(dCurrentComboSourceObjects.find(locObject) != dCurrentComboSourceObjects.end())
835 {
836 if(dDebugLevel > 20)
837 cout << "Source object already in use for locResumeAtIndex = " << locResumeAtIndex - 1 << endl;
838 continue;
839 }
840 return locObject; //not charged
841 }
842 while(locResumeAtIndex < int(locNeutralShowers.size()));
843 return NULL__null;
844}
845
846const JObject* DParticleComboBlueprint_factory::Grab_DetectedTrack(const DReaction* locReaction, Particle_t locAnalysisPID, vector<const DChargedTrack*>& locTracks, int& locResumeAtIndex) const
847{
848 if(dDebugLevel > 10)
849 cout << "Grab_DetectedTrack: resume at, #possible = " << locResumeAtIndex << ", " << locTracks.size() << endl;
850 if(locResumeAtIndex >= int(locTracks.size()))
851 return NULL__null;
852
853 do
854 {
855 const DChargedTrack* locChargedTrack = locTracks[locResumeAtIndex];
856 const JObject* locObject = static_cast<const JObject*>(locChargedTrack);
857 ++locResumeAtIndex;
858
859 //make sure not used currently
860 if(dCurrentComboSourceObjects.find(locObject) != dCurrentComboSourceObjects.end())
861 {
862 if(dDebugLevel > 20)
863 cout << "Source object already in use for locResumeAtIndex = " << locResumeAtIndex - 1 << endl;
864 continue;
865 }
866
867 //if charged, will make further cuts
868
869 const DChargedTrackHypothesis* locChargedTrackHypothesis = Get_ChargedHypothesisToUse(locChargedTrack, locAnalysisPID);
870 if(locChargedTrackHypothesis == NULL__null)
871 continue; //PID hypothesis failed a cut
872
873 //check to make sure the track momentum isn't too low (e.g. testing a 100 MeV pion to be a proton)
874 pair<bool, double> locMinProtonMomentum = dMinProtonMomentum.first ? dMinProtonMomentum : locReaction->Get_MinProtonMomentum();
875 if(locMinProtonMomentum.first && (ParticleMass(locChargedTrackHypothesis->PID()) < ParticleMass(Proton)) && (ParticleMass(locAnalysisPID) >= (ParticleMass(Proton) - 0.001)))
876 {
877 if(dDebugLevel > 20)
878 cout << "Proton candidate, momentum = " << locChargedTrackHypothesis->momentum().Mag() << endl;
879 if(locChargedTrackHypothesis->momentum().Mag() < locMinProtonMomentum.second)
880 {
881 if(dDebugLevel > 20)
882 cout << "Track momentum too low to be " << ParticleType(locAnalysisPID) << endl;
883 continue;
884 }
885 }
886
887 return locObject;
888 }
889 while(locResumeAtIndex < int(locTracks.size()));
890 return NULL__null;
891}
892
893const DChargedTrackHypothesis* DParticleComboBlueprint_factory::Get_ChargedHypothesisToUse(const DChargedTrack* locChargedTrack, Particle_t locAnalysisPID) const
894{
895 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(locAnalysisPID);
896 if(locChargedTrackHypothesis != NULL__null)
897 return locChargedTrackHypothesis;
898
899 if(dAvailablePIDs.find(locAnalysisPID) != dAvailablePIDs.end())
900 return NULL__null; //This PID was available: it must have been cut
901
902 //pid not found for this track: loop over other possible pids
903 deque<Particle_t> locPIDsToTry = dTrackTimeBasedFactory_Combo->Get_ParticleIDsToTry(locAnalysisPID);
904 for(size_t loc_i = 0; loc_i < locPIDsToTry.size(); ++loc_i)
905 {
906 locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(locPIDsToTry[loc_i]);
907 if(locChargedTrackHypothesis != NULL__null)
908 return locChargedTrackHypothesis;
909 }
910
911 //still none found, will have to take a new one: for now, take the one with the best FOM
912 return locChargedTrack->Get_BestFOM();
913}
914
915DParticleComboBlueprintStep* DParticleComboBlueprint_factory::Get_ParticleComboBlueprintStepResource(void)
916{
917 DParticleComboBlueprintStep* locParticleComboBlueprintStep;
918 if(dParticleComboBlueprintStepPool_Available.empty())
919 {
920 locParticleComboBlueprintStep = new DParticleComboBlueprintStep;
921 dParticleComboBlueprintStepPool_All.push_back(locParticleComboBlueprintStep);
922 }
923 else
924 {
925 locParticleComboBlueprintStep = dParticleComboBlueprintStepPool_Available.back();
926 locParticleComboBlueprintStep->Reset();
927 dParticleComboBlueprintStepPool_Available.pop_back();
928 }
929 return locParticleComboBlueprintStep;
930}
931
932void DParticleComboBlueprint_factory::Reset_Pools(void)
933{
934 // delete pool sizes if too large, preventing memory-leakage-like behavor.
935 if(dParticleComboBlueprintStepPool_All.size() > MAX_DParticleComboBlueprintStepPoolSize){
936 for(size_t loc_i = MAX_DParticleComboBlueprintStepPoolSize; loc_i < dParticleComboBlueprintStepPool_All.size(); ++loc_i)
937 delete dParticleComboBlueprintStepPool_All[loc_i];
938 dParticleComboBlueprintStepPool_All.resize(MAX_DParticleComboBlueprintStepPoolSize);
939 }
940 dParticleComboBlueprintStepPool_Available = dParticleComboBlueprintStepPool_All;
941}
942
943bool DParticleComboBlueprint_factory::Calc_FinalStateP4(size_t locTotalNumSteps, const DParticleComboBlueprint* locParticleComboBlueprint, const DParticleComboBlueprintStep* locNewParticleComboBlueprintStep, int locStepIndexToGrab, DLorentzVector& locFinalStateP4) const
944{
945 //The input locParticleComboBlueprint is under construction: it does not have all of the steps yet
946 //locNewParticleComboBlueprintStep is the step that will be added next, IF it passes this invariant mass cut
947 //This is a recursive function, so locStepIndexToGrab is the index for the step that should be used for the p4 calculation at this stage
948 //However, because the locParticleComboBlueprint is under construction, locStepIndexToGrab must be converted into the index that the step is actually located at
949 const DParticleComboBlueprintStep* locCurrentStep = locNewParticleComboBlueprintStep;
950 if(locStepIndexToGrab != -1)
951 {
952 int locStepFromBackIndex = locTotalNumSteps - locStepIndexToGrab;
953 int locActualStepIndex = locParticleComboBlueprint->Get_NumParticleComboBlueprintSteps() - locStepFromBackIndex;
954 locCurrentStep = locParticleComboBlueprint->Get_ParticleComboBlueprintStep(locActualStepIndex);
955 if(locCurrentStep == NULL__null)
956 return false;
957 }
958
959 locFinalStateP4.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
960 for(size_t loc_i = 0; loc_i < locCurrentStep->Get_NumFinalParticleSourceObjects(); ++loc_i)
961 {
962 if(locCurrentStep->Is_FinalParticleMissing(loc_i))
963 return false;
964 if(locCurrentStep->Is_FinalParticleDecaying(loc_i))
965 {
966 DLorentzVector locDecayP4;
967 if(!Calc_FinalStateP4(locTotalNumSteps, locParticleComboBlueprint, locNewParticleComboBlueprintStep, locCurrentStep->Get_DecayStepIndex(loc_i), locDecayP4))
968 return false;
969 locFinalStateP4 += locDecayP4;
970 }
971 else if(locCurrentStep->Is_FinalParticleCharged(loc_i))
972 {
973 const DChargedTrack* locChargedTrack = dynamic_cast<const DChargedTrack*>(locCurrentStep->Get_FinalParticle_SourceObject(loc_i));
974 Particle_t locPID = locCurrentStep->Get_FinalParticleID(loc_i);
975 const DChargedTrackHypothesis* locChargedTrackHypothesis = Get_ChargedHypothesisToUse(locChargedTrack, locPID);
976 DVector3 locMomentum(locChargedTrackHypothesis->momentum());
977 locFinalStateP4 += DLorentzVector(locMomentum, sqrt(locMomentum.Mag2() + ParticleMass(locPID)*ParticleMass(locPID)));
978 }
979 else //neutral
980 {
981 const DNeutralShower* locNeutralShower = dynamic_cast<const DNeutralShower*>(locCurrentStep->Get_FinalParticle_SourceObject(loc_i));
982 DVector3 locHitPoint = locNeutralShower->dSpacetimeVertex.Vect();
983 DVector3 locMomentum(locHitPoint - dVertex->dSpacetimeVertex.Vect());
984 Particle_t locPID = locCurrentStep->Get_FinalParticleID(loc_i);
985 if(locPID != Gamma)
986 {
987 double locDeltaT = locNeutralShower->dSpacetimeVertex.T() - dVertex->dSpacetimeVertex.T();
988 double locBeta = locMomentum.Mag()/(locDeltaT*29.9792458); //path length is locMomentum.Mag() (for now)
989 if(locBeta >= 1.0)
990 locBeta = 0.9999;
991 if(locBeta < 0.0)
992 locBeta = 0.0;
993 double locGamma = 1.0/sqrt(1.0 - locBeta*locBeta);
994 double locMass = ParticleMass(locPID);
995 double locPMag = locGamma*locBeta*locMass;
996 locMomentum.SetMag(locPMag);
997 locFinalStateP4 += DLorentzVector(locMomentum, sqrt(locMass*locMass + locPMag*locPMag));
998 }
999 else
1000 {
1001 locMomentum.SetMag(locNeutralShower->dEnergy);
1002 locFinalStateP4 += DLorentzVector(locMomentum, locNeutralShower->dEnergy);
1003 }
1004 }
1005 }
1006 return true;
1007}
1008
1009//------------------
1010// erun
1011//------------------
1012jerror_t DParticleComboBlueprint_factory::erun(void)
1013{
1014 return NOERROR;
1015}
1016
1017//------------------
1018// fini
1019//------------------
1020jerror_t DParticleComboBlueprint_factory::fini(void)
1021{
1022 for(size_t loc_i = 0; loc_i < dParticleComboBlueprintStepPool_All.size(); ++loc_i)
1023 delete dParticleComboBlueprintStepPool_All[loc_i];
1024
1025 return NOERROR;
1026}
1027