1 | #ifdef VTRACE |
2 | #include "vt_user.h" |
3 | #endif |
4 | |
5 | #include "DParticleComboBlueprint_factory.h" |
6 | |
7 | |
8 | |
9 | |
10 | jerror_t DParticleComboBlueprint_factory::init(void) |
11 | { |
12 | MAX_DParticleComboBlueprintStepPoolSize = 100; |
13 | |
14 | dDebugLevel = 0; |
15 | |
16 | |
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 | |
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 | |
56 | |
57 | jerror_t DParticleComboBlueprint_factory::brun(jana::JEventLoop* locEventLoop, int32_t runnumber) |
58 | { |
59 | |
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 | |
86 | void DParticleComboBlueprint_factory::Get_Reactions(JEventLoop *locEventLoop, vector<const DReaction*>& locReactions) const |
87 | { |
88 | |
89 | |
90 | |
91 | |
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 | |
102 | |
103 | |
104 | |
105 | vector<const DReaction*> locReactionsSubset; |
106 | locFactory->Get(locReactionsSubset); |
107 | locReactions.insert(locReactions.end(), locReactionsSubset.begin(), locReactionsSubset.end()); |
108 | } |
109 | } |
110 | |
111 | void 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 | |
137 | |
138 | jerror_t DParticleComboBlueprint_factory::evnt(JEventLoop *locEventLoop, uint64_t eventnumber) |
139 | { |
140 | #ifdef VTRACE |
141 | VT_TRACER("DParticleComboBlueprint_factory::evnt()"); |
142 | #endif |
143 | |
144 | |
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; |
173 | |
174 | |
175 | |
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 | |
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; |
202 | locSkimMissingFlag = true; |
| Value stored to 'locSkimMissingFlag' is never read |
203 | break; |
204 | } |
205 | |
206 | |
207 | } |
208 | |
209 | Build_ParticleComboBlueprints(locReactions[loc_i]); |
210 | } |
211 | |
212 | return NOERROR; |
213 | } |
214 | |
215 | void 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 | |
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 | |
236 | deque<deque<int> > locResumeAtIndexDeque; |
237 | deque<deque<int> > locNumPossibilitiesDeque; |
238 | map<int, int> locInitialParticleStepFromIndexMap; |
239 | map<pair<int, int>, int> locFinalStateDecayStepIndexMap; |
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 | |
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 | |
301 | bool 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 | |
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 | |
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 | |
333 | |
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; |
341 | locInitialParticleStepFromIndexMap[locDecayStepIndex] = loc_i; |
342 | continue; |
343 | } |
344 | |
345 | |
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 | |
369 | if((locNumNeededPositiveParticles > int(dPositiveChargedTracks.size())) || (locNumNeededNegativeParticles > int(dNegativeChargedTracks.size()))) |
370 | return false; |
371 | if((locNumNeededChargedParticles > int(dChargedTracks.size())) || (locNumNeededNeutralParticles > int(dNeutralShowers.size()))) |
372 | return false; |
373 | |
374 | |
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 | |
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 | |
393 | void 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; |
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 | |
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 | |
429 | if(locMissingParticleIndex == locParticleIndex) |
430 | { |
431 | if(dDebugLevel > 10) |
432 | cout << "missing particle" << endl; |
433 | |
434 | |
435 | locParticleComboBlueprintStep->Add_FinalParticle_SourceObject(NULL__null, -1); |
436 | locResumeAtIndex = 1; |
437 | ++locParticleIndex; |
438 | continue; |
439 | } |
440 | |
441 | |
442 | |
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); |
450 | locResumeAtIndex = 1; |
451 | ++locParticleIndex; |
452 | continue; |
453 | } |
454 | |
455 | |
456 | |
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 | |
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); |
485 | dCurrentComboSourceObjects.insert(locSourceObject); |
486 | ++locParticleIndex; |
487 | } |
488 | while(true); |
489 | delete locParticleComboBlueprint; |
490 | } |
491 | |
492 | bool 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 | |
495 | |
496 | |
497 | |
498 | |
499 | if(Check_IfDuplicateStepCombo(locParticleComboBlueprint, locParticleComboBlueprintStep, locStepIndex, locResumeAtIndexDeque, locNumPossibilitiesDeque)) |
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 | |
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 | |
529 | |
530 | |
531 | map<DParticleComboBlueprintStep, DParticleComboBlueprintStep*>::iterator locStepIterator = dBlueprintStepMap.find(*locParticleComboBlueprintStep); |
532 | if(locStepIterator != dBlueprintStepMap.end()) |
533 | { |
534 | |
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 | |
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 | |
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); |
569 | locParticleComboBlueprintStep = NULL__null; |
570 | |
571 | |
572 | if(locReaction->Get_AnyBlueprintFlag()) |
573 | return false; |
574 | |
575 | if(!Handle_Decursion(locParticleComboBlueprint, locResumeAtIndexDeque, locNumPossibilitiesDeque, locParticleIndex, locStepIndex, locParticleComboBlueprintStep)) |
576 | return false; |
577 | return true; |
578 | } |
579 | |
580 | bool 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) |
587 | { |
588 | if(dDebugLevel > 50) |
589 | cout << "just saved blueprint" << endl; |
590 | ++locStepIndex; |
591 | |
592 | |
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())) |
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 | |
613 | locResumeAtIndexDeque[locStepIndex][locParticleIndex] = 0; |
614 | --locParticleIndex; |
615 | if(locParticleIndex >= 0) |
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; |
624 | |
625 | const DParticleComboBlueprintStep* locPoppedStep = locParticleComboBlueprint->Pop_ParticleComboBlueprintStep(); |
626 | if(dSavedBlueprintSteps.find(locPoppedStep) == dSavedBlueprintSteps.end()) |
627 | { |
628 | |
629 | dParticleComboBlueprintStepPool_Available.push_back(locParticleComboBlueprintStep); |
630 | locParticleComboBlueprintStep = const_cast<DParticleComboBlueprintStep*>(locPoppedStep); |
631 | } |
632 | else |
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 | |
645 | bool 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 | |
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; |
665 | |
666 | |
667 | |
668 | |
669 | |
670 | |
671 | |
672 | |
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 | |
682 | if(dDebugLevel > 10) |
683 | cout << "step at index = " << locStepIndex << " matches a future step at index " << locFutureStepIndex << endl; |
684 | |
685 | |
686 | |
687 | |
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 | |
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; |
723 | } |
724 | } |
725 | |
726 | if(dDebugLevel > 10) |
727 | cout << "step at index = " << locStepIndex << " has no similar steps." << endl; |
728 | |
729 | return false; |
730 | } |
731 | |
732 | bool DParticleComboBlueprint_factory::Check_IfStepsAreIdentical(const DParticleComboBlueprint* locParticleComboBlueprint, const DParticleComboBlueprintStep* locCurrentStep, const DParticleComboBlueprintStep* locPreviousStep) const |
733 | { |
734 | if(locCurrentStep->Get_MissingParticleIndex() != -1) |
735 | return false; |
736 | if(locPreviousStep->Get_MissingParticleIndex() != -1) |
737 | return false; |
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 | |
748 | |
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 | |
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); |
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 | |
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; |
802 | } |
803 | } |
804 | |
805 | |
806 | return true; |
807 | } |
808 | |
809 | const 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 | |
822 | const 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 | |
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; |
841 | } |
842 | while(locResumeAtIndex < int(locNeutralShowers.size())); |
843 | return NULL__null; |
844 | } |
845 | |
846 | const 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 | |
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 | |
868 | |
869 | const DChargedTrackHypothesis* locChargedTrackHypothesis = Get_ChargedHypothesisToUse(locChargedTrack, locAnalysisPID); |
870 | if(locChargedTrackHypothesis == NULL__null) |
871 | continue; |
872 | |
873 | |
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 | |
893 | const 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; |
901 | |
902 | |
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 | |
912 | return locChargedTrack->Get_BestFOM(); |
913 | } |
914 | |
915 | DParticleComboBlueprintStep* 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 | |
932 | void DParticleComboBlueprint_factory::Reset_Pools(void) |
933 | { |
934 | |
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 | |
943 | bool DParticleComboBlueprint_factory::Calc_FinalStateP4(size_t locTotalNumSteps, const DParticleComboBlueprint* locParticleComboBlueprint, const DParticleComboBlueprintStep* locNewParticleComboBlueprintStep, int locStepIndexToGrab, DLorentzVector& locFinalStateP4) const |
944 | { |
945 | |
946 | |
947 | |
948 | |
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 |
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); |
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 | |
1011 | |
1012 | jerror_t DParticleComboBlueprint_factory::erun(void) |
1013 | { |
1014 | return NOERROR; |
1015 | } |
1016 | |
1017 | |
1018 | |
1019 | |
1020 | jerror_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 | |