11 gPARMS->SetDefaultParameter(
"DAnalysisUtilities:DEBUG_LEVEL",
DEBUG_LEVEL);
18 DGeometry *locGeometry = locApplication ? locApplication->
GetDGeometry(locEventLoop->GetJEvent().GetRunNumber()) : NULL;
36 VT_TRACER(
"DAnalysisUtilities::Check_IsBDTSignalEvent()");
53 vector<const DReaction*> locThrownReactions;
54 locEventLoop->Get(locThrownReactions,
"Thrown");
55 if(locThrownReactions.empty())
57 auto locActualThrownReaction = locThrownReactions[0];
60 size_t locStepIndex = 0;
61 vector<DReactionStep> locNewReactionSteps;
62 const DReaction* locCurrentReaction = locReaction;
79 for(
size_t loc_i = 0; loc_i < locStepIndex; ++loc_i)
85 bool locFoundFlag =
false;
86 for(
size_t loc_j = 0; loc_j < locPIDs.size(); ++loc_j)
89 if(locFinalStatePID != locInitialPID)
91 locPIDs.erase(locPIDs.begin() + loc_j);
92 locPIDs.insert(locPIDs.begin(), locDecayProducts.begin(), locDecayProducts.end());
104 for(
int loc_j = 0; loc_j < int(locPIDs.size()); ++loc_j)
105 locNewReactionStep.
Add_FinalParticleID(locPIDs[loc_j], (loc_j == locMissingParticleIndex));
106 locNewReactionSteps.push_back(locNewReactionStep);
115 else if(loc_j == locStepIndex)
120 locNewReaction = locReactionBuild;
121 locCurrentReaction = &locNewReaction;
125 while(locStepIndex < locCurrentReaction->Get_NumReactionSteps());
128 deque<pair<const DMCThrown*, deque<const DMCThrown*> > > locThrownSteps;
131 if(locThrownSteps.size() == 1)
135 map<Particle_t, size_t> locNumDecayingParticles_Thrown;
136 for(
size_t loc_i = 0; loc_i < locThrownSteps.size(); ++loc_i)
138 if(locThrownSteps[loc_i].first == NULL)
140 Particle_t locPID = locThrownSteps[loc_i].first->PID();
141 if(locNumDecayingParticles_Thrown.find(locPID) == locNumDecayingParticles_Thrown.end())
142 locNumDecayingParticles_Thrown[locPID] = 1;
144 ++locNumDecayingParticles_Thrown[locPID];
147 map<Particle_t, size_t> locNumDecayingParticles_Reaction;
153 if(locNumDecayingParticles_Reaction.find(locPID) == locNumDecayingParticles_Reaction.end())
154 locNumDecayingParticles_Reaction[locPID] = 1;
156 ++locNumDecayingParticles_Reaction[locPID];
160 map<Particle_t, size_t>::iterator locIterator = locNumDecayingParticles_Reaction.begin();
161 for(; locIterator != locNumDecayingParticles_Reaction.end(); ++locIterator)
164 if(locNumDecayingParticles_Thrown.find(locPID) == locNumDecayingParticles_Thrown.end())
173 pair<const DMCThrown*, deque<const DMCThrown*> > locStepPair = locThrownSteps[locStepIndex];
176 const DMCThrown* locThrownParent = locStepPair.first;
178 bool locInitialPIDFoundFlag = (locNumDecayingParticles_Reaction.find(locInitialPID) != locNumDecayingParticles_Reaction.end());
182 if(((!locInitialPIDFoundFlag) && locIncludeDecayingToReactionFlag) || (locInitialPID ==
omega) || (locInitialPID ==
phiMeson))
185 if(locNumDecayingParticles_Thrown[locInitialPID] == 1)
186 locNumDecayingParticles_Thrown.erase(locInitialPID);
188 --locNumDecayingParticles_Thrown[locInitialPID];
193 while(locStepIndex < locThrownSteps.size());
195 if(!locIncludeDecayingToReactionFlag)
200 bool locCheckResult =
Check_ThrownsMatchReaction(locActualThrownReaction, locThrownCombo, locCurrentReaction, locExclusiveMatchFlag);
204 return locCheckResult;
208 if(locExclusiveMatchFlag)
219 vector<Particle_t> locPIDVector;
220 vector<int> locResumeAtIndex;
221 for(locIterator = locNumDecayingParticles_Reaction.begin(); locIterator != locNumDecayingParticles_Reaction.end(); ++locIterator)
224 size_t locNumThrown = locNumDecayingParticles_Thrown[locPID];
225 if(locNumThrown < locIterator->second)
227 if(locNumThrown == locIterator->second)
229 for(
size_t loc_i = 0; loc_i < locNumThrown - locIterator->second; ++loc_i)
231 locPIDVector.push_back(locPID);
232 locResumeAtIndex.push_back(0);
237 if(locPIDVector.empty())
241 bool locCheckResult =
Check_ThrownsMatchReaction(locActualThrownReaction, locThrownCombo, locCurrentReaction, locExclusiveMatchFlag);
245 return locCheckResult;
251 deque<int> locDecayReplacementIndices;
252 int locParticleIndex = 0;
255 deque<deque<pair<const DMCThrown*, deque<const DMCThrown*> > > > locReplacementThrownSteps(1, locThrownSteps);
258 if(locParticleIndex == -1)
261 deque<pair<const DMCThrown*, deque<const DMCThrown*> > > locCurrentThrownSteps = locReplacementThrownSteps.back();
262 if(locParticleIndex ==
int(locPIDVector.size()))
267 bool locCheckResult =
Check_ThrownsMatchReaction(locActualThrownReaction, locThrownCombo, locCurrentReaction, locExclusiveMatchFlag);
274 locReplacementThrownSteps.pop_back();
279 Particle_t locPID = locPIDVector[locParticleIndex];
281 bool locFoundFlag =
false;
282 for(
size_t loc_i = locResumeAtIndex[locParticleIndex]; loc_i < locCurrentThrownSteps.size(); ++loc_i)
284 if(locCurrentThrownSteps[loc_i].first == NULL)
286 Particle_t locInitialPID = locCurrentThrownSteps[loc_i].first->PID();
287 if(locInitialPID != locPID)
292 locReplacementThrownSteps.push_back(locCurrentThrownSteps);
293 locResumeAtIndex[locParticleIndex] = loc_i + 1;
302 locResumeAtIndex[locParticleIndex] = 0;
303 locReplacementThrownSteps.pop_back();
314 if(locStepIndex >= locThrownSteps.size())
318 int locInitialParticleDecayFromStepIndex = -1;
319 const DMCThrown* locThrownParent = locThrownSteps[locStepIndex].first;
320 if(locThrownParent == NULL)
323 for(
size_t loc_j = 0; loc_j < locStepIndex; ++loc_j)
325 for(
size_t loc_k = 0; loc_k < locThrownSteps[loc_j].second.size(); ++loc_k)
327 if(locThrownParent != locThrownSteps[loc_j].second[loc_k])
329 locInitialParticleDecayFromStepIndex = loc_j;
332 if(locInitialParticleDecayFromStepIndex != -1)
336 if(locInitialParticleDecayFromStepIndex == -1)
338 cout <<
"ERROR: SOMETHING IS WRONG WITH DAnalysisUtilities::Replace_DecayingParticleWithProducts(). ABORTING." << endl;
343 deque<const DMCThrown*>& locProductionStepFinalParticles = locThrownSteps[locInitialParticleDecayFromStepIndex].second;
344 deque<const DMCThrown*>& locDecayStepFinalParticles = locThrownSteps[locStepIndex].second;
345 for(
size_t loc_j = 0; loc_j < locProductionStepFinalParticles.size(); ++loc_j)
347 if(locProductionStepFinalParticles[loc_j] != locThrownParent)
349 locProductionStepFinalParticles.erase(locProductionStepFinalParticles.begin() + loc_j);
350 locProductionStepFinalParticles.insert(locProductionStepFinalParticles.end(), locDecayStepFinalParticles.begin(), locDecayStepFinalParticles.end());
354 locThrownSteps.erase(locThrownSteps.begin() + locStepIndex);
367 vector<const DReaction*> locReactions;
368 locEventLoop->Get(locReactions,
"Thrown");
369 if(locReactions.empty())
380 VT_TRACER(
"DAnalysisUtilities::Check_ThrownsMatchReaction()");
388 if(locThrownCombo == NULL)
390 if(locExclusiveMatchFlag)
397 map<size_t, int> locReactionInitialParticleDecayFromStepIndexMap;
404 locReactionInitialParticleDecayFromStepIndexMap[0] = -1;
415 if(locReactionInitialParticleDecayFromStepIndexMap.find(loc_k) != locReactionInitialParticleDecayFromStepIndexMap.end())
417 locReactionInitialParticleDecayFromStepIndexMap[loc_k] = loc_i;
424 set<size_t> locMatchedInputStepIndices;
425 map<int, int> locStepMatching;
426 map<int, int> locReverseStepMatching;
432 size_t locStartSearchIndex = 0;
433 if(locReverseStepMatching.find(locInitialParticleDecayFromStepIndex) != locReverseStepMatching.end())
434 locStartSearchIndex = locReverseStepMatching[locInitialParticleDecayFromStepIndex] + 1;
437 bool locMatchFoundFlag =
false;
438 int locPossibleMatchIndex = -1;
442 if(locMatchedInputStepIndices.find(loc_j) != locMatchedInputStepIndices.end())
450 if(locReactionInitialParticleDecayFromStepIndexMap.find(loc_j) == locReactionInitialParticleDecayFromStepIndexMap.end())
454 locPossibleMatchIndex = loc_j;
458 int locReactionInitialParticleDecayFromStepIndex = locReactionInitialParticleDecayFromStepIndexMap[loc_j];
459 if(locReactionInitialParticleDecayFromStepIndex != -1)
461 if(locStepMatching.find(locReactionInitialParticleDecayFromStepIndex) == locStepMatching.end())
463 int locReactionInitialParticleDecayFromStepIndexMappedBackToThrown = locStepMatching[locReactionInitialParticleDecayFromStepIndex];
464 if(locInitialParticleDecayFromStepIndex != locReactionInitialParticleDecayFromStepIndexMappedBackToThrown)
467 else if(locInitialParticleDecayFromStepIndex != -1)
471 locMatchedInputStepIndices.insert(loc_j);
472 locStepMatching[loc_j] = loc_i;
473 locReverseStepMatching[loc_i] = loc_j;
474 locMatchFoundFlag =
true;
478 if((!locMatchFoundFlag) && (locPossibleMatchIndex != -1))
481 locMatchedInputStepIndices.insert(locPossibleMatchIndex);
482 locStepMatching[locPossibleMatchIndex] = loc_i;
483 locReverseStepMatching[loc_i] = locPossibleMatchIndex;
484 locMatchFoundFlag =
true;
486 if(locExclusiveMatchFlag && (!locMatchFoundFlag))
489 if(locExclusiveMatchFlag)
495 if(locMatchedInputStepIndices.find(loc_i) == locMatchedInputStepIndices.end())
504 locUnusedChargedTracks.clear();
505 locEventLoop->Get(locUnusedChargedTracks,
"Combo");
506 std::sort(locUnusedChargedTracks.begin(), locUnusedChargedTracks.end());
509 for(
size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
511 for(
auto locIterator = locUnusedChargedTracks.begin(); locIterator != locUnusedChargedTracks.end();)
513 if(locChargedSourceObjects[loc_i] != *locIterator)
516 locIterator = locUnusedChargedTracks.erase(locIterator);
523 locUnusedTimeBasedTracks.clear();
524 locEventLoop->Get(locUnusedTimeBasedTracks);
526 vector<const DTrackTimeBased*> locComboTimeBasedTracks;
527 locEventLoop->Get(locComboTimeBasedTracks,
"Combo");
528 locUnusedTimeBasedTracks.insert(locUnusedTimeBasedTracks.end(), locComboTimeBasedTracks.begin(), locComboTimeBasedTracks.end());
531 for(
size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
534 auto locChargedTrack =
static_cast<const DChargedTrack*
>(locChargedSourceObjects[loc_i]);
535 for(
auto locIterator = locUnusedTimeBasedTracks.begin(); locIterator != locUnusedTimeBasedTracks.end();)
537 if(locChargedTrack->candidateid != (*locIterator)->candidateid)
540 locIterator = locUnusedTimeBasedTracks.erase(locIterator);
547 locUnusedWireBasedTracks.clear();
548 locEventLoop->Get(locUnusedWireBasedTracks);
551 for(
size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
554 auto locChargedTrack =
static_cast<const DChargedTrack*
>(locChargedSourceObjects[loc_i]);
555 for(
auto locIterator = locUnusedWireBasedTracks.begin(); locIterator != locUnusedWireBasedTracks.end();)
557 if(locChargedTrack->candidateid != (*locIterator)->candidateid)
560 locIterator = locUnusedWireBasedTracks.erase(locIterator);
567 locUnusedTrackCandidates.clear();
568 locEventLoop->Get(locUnusedTrackCandidates);
571 set<unsigned int> locUsedCandidateIndices;
572 for(
size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
575 auto locChargedTrack =
static_cast<const DChargedTrack*
>(locChargedSourceObjects[loc_i]);
576 locUsedCandidateIndices.insert(locChargedTrack->candidateid - 1);
579 for(
int loc_i = locUnusedTrackCandidates.size() - 1; loc_i >= 0; --loc_i)
581 if(locUsedCandidateIndices.find(loc_i) != locUsedCandidateIndices.end())
582 locUnusedTrackCandidates.erase(locUnusedTrackCandidates.begin() + loc_i);
588 locUnusedNeutralShowers.clear();
592 for(
size_t loc_i = 0; loc_i < locNeutralSourceObjects.size(); ++loc_i)
594 auto locNeutralShower =
static_cast<const DNeutralShower*
>(locNeutralSourceObjects[loc_i]);
595 for(
auto locIterator = locUnusedNeutralShowers.begin(); locIterator != locUnusedNeutralShowers.end();)
597 if(locNeutralShower != *locIterator)
600 locIterator = locUnusedNeutralShowers.erase(locIterator);
607 locUnusedNeutralParticles.clear();
611 for(
size_t loc_i = 0; loc_i < locNeutralSourceObjects.size(); ++loc_i)
613 auto locNeutralShower =
static_cast<const DNeutralShower*
>(locNeutralSourceObjects[loc_i]);
614 for(
auto locIterator = locUnusedNeutralParticles.begin(); locIterator != locUnusedNeutralParticles.end();)
616 if(locNeutralShower != (*locIterator)->dNeutralShower)
619 locIterator = locUnusedNeutralParticles.erase(locIterator);
626 vector<const DMCThrown*> locMCThrowns;
627 locEventLoop->Get(locMCThrowns);
629 map<size_t, const DMCThrown*> locIDMap;
630 for(
size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
631 locIDMap[locMCThrowns[loc_i]->myid] = locMCThrowns[loc_i];
633 locThrownSteps.clear();
634 if(locMCThrowns.empty())
636 locThrownSteps.push_back(pair<
const DMCThrown*, deque<const DMCThrown*> >(NULL, deque<const DMCThrown*>() ) );
638 for(
size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
643 if(locMCThrowns[loc_i]->PID() ==
Unknown)
647 int locParentID = locMCThrowns[loc_i]->parentid;
650 locThrownSteps[0].second.push_back(locMCThrowns[loc_i]);
653 if(locIDMap.find(locParentID) == locIDMap.end())
657 Particle_t locParentPID = locIDMap[locParentID]->PID();
658 bool locDoneFlag =
false;
662 locParentID = locIDMap[locParentID]->parentid;
665 locThrownSteps[0].second.push_back(locMCThrowns[loc_i]);
668 else if(locIDMap.find(locParentID) == locIDMap.end())
671 locParentPID = locIDMap[locParentID]->PID();
678 bool locListedAsDecayingFlag =
false;
679 for(
size_t loc_j = 1; loc_j < locThrownSteps.size(); ++loc_j)
681 if(locThrownSteps[loc_j].first->myid != locParentID)
683 locThrownSteps[loc_j].second.push_back(locMCThrowns[loc_i]);
684 locListedAsDecayingFlag =
true;
687 if(locListedAsDecayingFlag)
692 const DMCThrown* locThrownParent = locIDMap[locParentID];
693 bool locFoundFlag =
false;
694 for(
size_t loc_j = 0; loc_j < locThrownSteps.size(); ++loc_j)
696 for(
size_t loc_k = 0; loc_k < locThrownSteps[loc_j].second.size(); ++loc_k)
698 if(locThrownSteps[loc_j].second[loc_k] != locThrownParent)
710 locThrownSteps.push_back(pair<
const DMCThrown*, deque<const DMCThrown*> >(locThrownParent, deque<const DMCThrown*>(1, locMCThrowns[loc_i]) ));
727 vector<const DMCThrown*> locMCThrowns;
728 locEventLoop->Get(locMCThrowns,
"FinalState");
729 deque<Particle_t> locDesiredPIDs_Copy = locDesiredPIDs;
731 bool locMissingPIDMatchedFlag =
false;
732 for(
size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
736 if((!locMissingPIDMatchedFlag) && (locMissingPID == locPID))
739 locMissingPIDMatchedFlag =
true;
743 bool locPIDFoundFlag =
false;
744 for(deque<Particle_t>::iterator locIterator = locDesiredPIDs_Copy.begin(); locIterator != locDesiredPIDs_Copy.end(); ++locIterator)
746 if(*locIterator != locPID)
748 locDesiredPIDs_Copy.erase(locIterator);
749 locPIDFoundFlag =
true;
756 return (locDesiredPIDs_Copy.empty());
761 set<pair<const JObject*, unsigned int> > locSourceObjects;
762 return Calc_MissingP4(locReaction, locParticleCombo, 0, -1, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
767 return Calc_MissingP4(locReaction, locParticleCombo, 0, -1, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
772 set<pair<const JObject*, unsigned int> > locSourceObjects;
773 return Calc_MissingP4(locReaction, locParticleCombo, locStepIndex, locUpToStepIndex, locUpThroughIndices, locSourceObjects, locUseKinFitDataFlag);
781 return Calc_MissingP4(locReaction, locParticleCombo, locStepIndex, locUpToStepIndex, locUpThroughIndices, locSourceObjects,
false);
788 if(locStepIndex == 0)
792 locSourceObjects.insert(pair<const JObject*, unsigned int>(locKinematicData, abs(
PDGtype(locKinematicData->
PID()))));
793 if(locUseKinFitDataFlag)
799 Particle_t locPID = locReactionStep->Get_TargetPID();
807 for(
size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
809 if(
int(loc_j) == locReactionStep->Get_MissingParticleIndex())
811 if((
int(locStepIndex) == locUpToStepIndex) && (locUpThroughIndices.find(loc_j) == locUpThroughIndices.end()))
815 if(locDecayStepIndex > 0)
818 locMissingP4 +=
Calc_MissingP4(locReaction, locParticleCombo, locDecayStepIndex, locUpToStepIndex, locUpThroughIndices, locSourceObjects, locUseKinFitDataFlag);
822 Particle_t locPID = locReactionStep->Get_FinalPID(loc_j);
823 auto locDetectedP4 = locParticles[loc_j]->lorentzMomentum();
824 locMissingP4 -= locDetectedP4;
849 TMatrixFSym locMissingCovarianceMatrix(3);
850 locMissingCovarianceMatrix.Zero();
853 if(locStepIndex == 0)
857 TMatrixFSym locParticleCovarianceMatrix = *(locKinematicData->
errorMatrix().get());
858 locParticleCovarianceMatrix.ResizeTo(3, 3);
859 locMissingCovarianceMatrix += locParticleCovarianceMatrix;
863 for(
size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
865 if(
int(loc_j) == locReactionStep->Get_MissingParticleIndex())
867 if((
int(locStepIndex) == locUpToStepIndex) && (locUpThroughIndices.find(loc_j) == locUpThroughIndices.end()))
871 if(locDecayStepIndex > 0)
872 locMissingCovarianceMatrix +=
Calc_MissingP3Covariance(locReaction, locParticleCombo, locDecayStepIndex, locUpToStepIndex, locUpThroughIndices);
875 TMatrixFSym locParticleCovarianceMatrix = *(locParticles[loc_j]->errorMatrix().get());
876 locParticleCovarianceMatrix.ResizeTo(3, 3);
877 locMissingCovarianceMatrix += locParticleCovarianceMatrix;
881 return locMissingCovarianceMatrix;
886 set<pair<const JObject*, unsigned int> > locSourceObjects;
887 return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
892 return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
897 set<pair<const JObject*, unsigned int> > locSourceObjects;
898 return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, locToIncludeIndices, locSourceObjects, locUseKinFitDataFlag);
904 return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, locToIncludeIndices, locSourceObjects,
false);
908 if(locParticleComboStep == NULL)
915 if(locStepIndex != 0)
917 Particle_t locPID = locReactionStep->Get_TargetPID();
922 bool locDoSubsetFlag = !locToIncludeIndices.empty();
923 for(
size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
925 if(locDoSubsetFlag && (locToIncludeIndices.find(loc_i) == locToIncludeIndices.end()))
928 if(locReactionStep->Get_MissingParticleIndex() == int(loc_i))
932 if(locDecayStepIndex >= 0)
935 if((!locUseKinFitDataFlag) || (!
IsFixedMass(locReactionStep->Get_FinalPID(loc_i))))
936 locFinalStateP4 +=
Calc_FinalStateP4(locReaction, locParticleCombo, locDecayStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
939 locFinalStateP4 += locParticles[loc_i]->lorentzMomentum();
941 Calc_FinalStateP4(locReaction, locParticleCombo, locDecayStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
946 locFinalStateP4 += locParticles[loc_i]->lorentzMomentum();
950 return locFinalStateP4;
957 double locRFTime = (locEventRFBunch != NULL) ? locEventRFBunch->
dTime : numeric_limits<double>::quiet_NaN();
959 vector<const DNeutralShower*> locUnusedNeutralShowers;
962 double locEnergy_UnusedShowers = 0.;
963 for(
size_t loc_i = 0; loc_i < locUnusedNeutralShowers.size(); ++loc_i) {
964 const DNeutralShower* locUnusedNeutralShower = locUnusedNeutralShowers[loc_i];
968 double locDeltaT = locUnusedNeutralShower->
dSpacetimeVertex.T() - locFlightTime - locRFTime;
969 double locDetectorTheta = (locUnusedNeutralShower->
dSpacetimeVertex.Vect()-locVertex).Theta()*180./TMath::Pi();
970 if(locDetectorTheta < 2.0 || fabs(locDeltaT) > 4.)
973 locEnergy_UnusedShowers += locUnusedNeutralShower->
dEnergy;
976 return locEnergy_UnusedShowers;
981 vector<const DChargedTrack*> locUnusedChargedTracks;
984 for(
size_t loc_i = 0; loc_i < locUnusedChargedTracks.size(); ++loc_i) {
985 const DChargedTrack* locUnusedChargedTrack = locUnusedChargedTracks[loc_i];
988 locSumPMag_UnusedTracks += locUnusedChargedTrackHypothesis->
pmag();
989 locSumP3_UnusedTracks += locUnusedChargedTrackHypothesis->
momentum();
992 return (
int)locUnusedChargedTracks.size();
1024 return Calc_DOCA(locUnitDir, locUnitDir, locPosition, locVertex, locPOCA1, locPOCA2);
1030 return Calc_DOCA(locUnitDir, locUnitDir, locPosition, locVertex, locPOCA, locPOCA2);
1039 locPOCA1,locPOCA2,locDOCA)==NOERROR){
1040 locDOCAVertex=0.5*(locPOCA1+locPOCA2);
1049 return Calc_DOCAVertex(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locDOCAVertex);
1058 locPOCA1,locPOCA2,locDOCA)==NOERROR){
1059 locDOCAVertex=0.5*(locPOCA1+locPOCA2);
1068 return Calc_DOCAVertex(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locDOCAVertex);
1074 double locDOCA =
Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1075 locDOCAVertex = 0.5*(locPOCA1 + locPOCA2);
1082 return Calc_DOCA(locKinFitParticle1, locKinFitParticle2, locPOCA1, locPOCA2);
1088 return Calc_DOCA(locKinematicData1, locKinematicData2, locPOCA1, locPOCA2);
1097 return Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1106 return Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1112 return Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1119 double locUnitDot = locUnitDir1.Dot(locUnitDir2);
1120 double locDenominator = locUnitDot*locUnitDot - 1.0;
1121 double locDistVertToInterDOCA1 = 0.0, locDistVertToInterDOCA2 = 0.0;
1123 if(fabs(locDenominator) < 1.0
e-15)
1125 locDistVertToInterDOCA1 = (locVertex2 - locVertex1).Dot(locUnitDir2)/locUnitDot;
1126 locDistVertToInterDOCA2 = (locVertex1 - locVertex2).Dot(locUnitDir1)/locUnitDot;
1130 double locA = (locVertex1 - locVertex2).Dot(locUnitDir1);
1131 double locB = (locVertex1 - locVertex2).Dot(locUnitDir2);
1132 locDistVertToInterDOCA1 = (locA - locUnitDot*locB)/locDenominator;
1133 locDistVertToInterDOCA2 = (locUnitDot*locA - locB)/locDenominator;
1136 locPOCA1 = locVertex1 + locDistVertToInterDOCA1*locUnitDir1;
1137 locPOCA2 = locVertex2 + locDistVertToInterDOCA2*locUnitDir2;
1138 return (locPOCA1 - locPOCA2).Mag();
1148 double &doca)
const{
1161 return Calc_DOCA(q1,q2,pos1_in,pos2_in,mom1_in,mom2_in,pos1_out,pos2_out,
1174 double q1=kinematicData1->
charge();
1175 double q2=kinematicData2->
charge();
1185 return Calc_DOCA(q1,q2,pos1_in,pos2_in,mom1_in,mom2_in,pos1_out,pos2_out,
1199 DVector3 avg_pos=0.5*(pos1_in+pos2_in);
1202 double kap1=-0.003*B*q1/(2.*mom1_in.Perp());
1203 double kap2=-0.003*B*q2/(2.*mom2_in.Perp());
1204 double dx=diff_pos.x(),dy=diff_pos.y(),dz=diff_pos.z();
1205 double tan1=tan(M_PI_2-mom1_in.Theta());
1206 double tan1sq=tan1*tan1;
1207 double phi1=mom1_in.Phi();
1208 double cos1=cos(phi1),sin1=
sin(phi1);
1209 double cos1sq=cos1*cos1,sin1sq=sin1*sin1;
1210 double tan2=tan(M_PI_2-mom2_in.Theta());
1211 double tan2sq=tan2*tan2;
1212 double phi2=mom2_in.Phi();
1213 double cos2=cos(phi2),sin2=
sin(phi2);
1214 double cos2sq=cos2*cos2,sin2sq=sin2*sin2;
1215 double dx_sq=dx*dx,dy_sq=dy*dy;
1216 double cos2sin1_minus_cos1sin2=cos2*sin1 - cos1*sin2;
1218 double B1=pow(2*dx*kap2*sin1sq*sin2 - 2*dx*kap1*sin1*sin2sq - sin1sq*sin2sq
1219 + tan1sq + 2*dx*kap2*sin2*tan1sq - 2*dx*kap1*sin1*tan2sq
1220 - 2*sin1*sin2*(2*dx_sq*kap1*kap2 + tan1*tan2)
1221 + sin1sq*(1 + tan2sq) + cos1sq*(-2*cos2*dy*kap2
1222 + 4*dx*kap2*sin2 + sin2sq
1224 + 2*cos2*(-2*dy*kap2*sin1sq + dy*kap1*sin1*sin2
1225 - dy*kap2*tan1sq + sin1*(2*dx*dy*kap1*kap2
1228 + 2*cos1*(cos2sq*dy*kap1 + dy*kap2*sin1*sin2 + dy*kap1*tan2sq
1229 + sin2*(2*dx*dy*kap1*kap2 + dz*kap2*tan1
1231 - cos2*(2*dy_sq*kap1*kap2 + dx*kap2*sin1
1232 + dx*kap1*sin2 + sin1*sin2 + tan1*tan2)),2)
1233 + 8*(cos2sin1_minus_cos1sin2)*(cos1*kap1*(cos2 - 2*dy*kap2)
1234 + kap2*(-1 + 2*dx*kap1*sin1 - tan1sq)
1235 + kap1*(sin1*sin2 + tan1*tan2))*
1236 (2*cos2*dy_sq*kap2*sin1 + cos2*dx*sin1*sin2 - dz*tan1
1237 - 2*dx*dz*kap2*sin2*tan1 + dz*sin1*sin2*tan2 + cos2*dx*tan1*tan2
1238 + dy*(-2*dx*kap2*sin1*sin2 + sin1*sin2sq + 2*cos2*dz*kap2*tan1
1239 + sin2*tan1*tan2 - sin1*(1 + tan2sq))
1240 + cos1*(cos2*(2*dx*dy*kap2 + dy*sin2 + dz*tan2)
1241 - dx*(2*dx*kap2*sin2 + sin2sq + tan2sq)));
1244 return VALUE_OUT_OF_RANGE;
1247 double B2=pow(2*dx*kap2*sin1sq*sin2 - 2*dx*kap1*sin1*sin2sq - sin1sq*sin2sq
1248 + tan1sq + 2*dx*kap2*sin2*tan1sq - 2*dx*kap1*sin1*tan2sq
1249 - 2*sin1*sin2*(2*dx_sq*kap1*kap2 + tan1*tan2)
1250 + sin1sq*(1 + tan2sq) + cos1sq*(-2*cos2*dy*kap2
1251 + 4*dx*kap2*sin2 + sin2sq
1253 + 2*cos2*(-2*dy*kap2*sin1sq + dy*kap1*sin1*sin2 -dy*kap2*tan1sq
1254 + sin1*(2*dx*dy*kap1*kap2 - dz*kap2*tan1
1256 + 2*cos1*(cos2sq*dy*kap1 + dy*kap2*sin1*sin2 + dy*kap1*tan2sq
1257 + sin2*(2*dx*dy*kap1*kap2 + dz*kap2*tan1
1259 - cos2*(2*dy_sq*kap1*kap2 + dx*kap2*sin1
1260 + dx*kap1*sin2 + sin1*sin2 + tan1*tan2)),2)
1261 + 8*(cos2sin1_minus_cos1sin2)*(cos1*kap1*(cos2 - 2*dy*kap2)
1262 + kap2*(-1 + 2*dx*kap1*sin1 - tan1sq)
1263 + kap1*(sin1*sin2 + tan1*tan2))
1264 * (2*cos2*dy_sq*kap2*sin1 + cos2*dx*sin1*sin2 - dz*tan1
1265 - 2*dx*dz*kap2*sin2*tan1 + dz*sin1*sin2*tan2 + cos2*dx*tan1*tan2
1266 + dy*(-2*dx*kap2*sin1*sin2 + sin1*sin2sq + 2*cos2*dz*kap2*tan1
1267 + sin2*tan1*tan2 - sin1*(1 + tan2sq))
1268 + cos1*(cos2*(2*dx*dy*kap2 + dy*sin2 + dz*tan2)
1269 - dx*(2*dx*kap2*sin2 + sin2sq + tan2sq)));
1272 return VALUE_OUT_OF_RANGE;
1275 double A1=2*cos1*cos2sq*dy*kap1 - 2*cos1sq*cos2*dy*kap2
1276 - 4*cos1*cos2*dy_sq*kap1*kap2 - 2*cos1*cos2*dx*kap2*sin1
1277 + 4*cos2*dx*dy*kap1*kap2*sin1 + cos2sq*sin1sq - 4*cos2*dy*kap2*sin1sq
1278 - 2*cos1*cos2*dx*kap1*sin2 + 4*cos1sq*dx*kap2*sin2
1279 + 4*cos1*dx*dy*kap1*kap2*sin2 - 2*cos1*cos2*sin1*sin2
1280 + 2*cos2*dy*kap1*sin1*sin2 + 2*cos1*dy*kap2*sin1*sin2
1281 - 4*dx_sq*kap1*kap2*sin1*sin2 + 2*dx*kap2*sin1sq*sin2 + cos1sq*sin2sq
1282 - 2*dx*kap1*sin1*sin2sq - 2*cos2*dz*kap2*sin1*tan1
1283 + 2*cos1*dz*kap2*sin2*tan1 + cos2sq*tan1sq - 2*cos2*dy*kap2*tan1sq
1284 + 2*dx*kap2*sin2*tan1sq + sin2sq*tan1sq + 2*cos2*dz*kap1*sin1*tan2
1285 - 2*cos1*dz*kap1*sin2*tan2 - 2*cos1*cos2*tan1*tan2 - 2*sin1*sin2*tan1*tan2
1286 + cos1sq*tan2sq + 2*cos1*dy*kap1*tan2sq - 2*dx*kap1*sin1*tan2sq
1288 double C1=-4.*(cos2sin1_minus_cos1sin2)
1289 *(cos1*kap1*(cos2 - 2*dy*kap2) + kap2*(-1 + 2*dx*kap1*sin1 - tan1sq)
1290 + kap1*(sin1*sin2 + tan1*tan2));
1292 double A2=2*cos1*cos2sq*dy*kap1 - 2*cos1sq*cos2*dy*kap2
1293 - 4*cos1*cos2*dy_sq*kap1*kap2 - 4*cos2sq*dx*kap1*sin1
1294 + 2*cos1*cos2*dx*kap2*sin1 + 4*cos2*dx*dy*kap1*kap2*sin1 + cos2sq*sin1sq
1295 + 2*cos1*cos2*dx*kap1*sin2 + 4*cos1*dx*dy*kap1*kap2*sin2
1296 - 2*cos1*cos2*sin1*sin2 - 2*cos2*dy*kap1*sin1*sin2
1297 - 2*cos1*dy*kap2*sin1*sin2 - 4*dx_sq*kap1*kap2*sin1*sin2
1298 + 2*dx*kap2*sin1sq*sin2 + cos1sq*sin2sq + 4*cos1*dy*kap1*sin2sq
1299 - 2*dx*kap1*sin1*sin2sq + 2*cos2*dz*kap2*sin1*tan1
1300 - 2*cos1*dz*kap2*sin2*tan1 + cos2sq*tan1sq - 2*cos2*dy*kap2*tan1sq
1301 + 2*dx*kap2*sin2*tan1sq + sin2sq*tan1sq - 2*cos2*dz*kap1*sin1*tan2
1302 + 2*cos1*dz*kap1*sin2*tan2 - 2*cos1*cos2*tan1*tan2 - 2*sin1*sin2*tan1*tan2
1303 + cos1sq*tan2sq + 2*cos1*dy*kap1*tan2sq - 2*dx*kap1*sin1*tan2sq
1305 double C2=4.*(cos2sin1_minus_cos1sin2)
1306 *(kap2*(cos1*cos2 + sin1*sin2 + tan1*tan2)
1307 + kap1*(-1 + 2*cos2*dy*kap2 - 2*dx*kap2*sin2 - tan2sq));
1310 double s1_minus=(A1-
sqrt(B1))/C1;
1311 double s2_minus=(A2-
sqrt(B2))/C2;
1312 double s1_plus=(A1+
sqrt(B1))/C1;
1313 double s2_plus=(A2+
sqrt(B2))/C2;
1316 double x1=pos1_in.x();
1317 double y1=pos1_in.y();
1318 double z1=pos1_in.z();
1319 double x2=pos2_in.x();
1320 double y2=pos2_in.y();
1321 double z2=pos2_in.z();
1324 if (fabs(s1_minus+s2_minus)<fabs(s1_plus+s2_plus)){
1325 pos1_out.SetXYZ(x1+cos1*s1_minus,y1+sin1*s1_minus,z1+tan1*s1_minus);
1326 pos2_out.SetXYZ(x2+cos2*s2_minus,y2+sin2*s2_minus,z2+tan2*s2_minus);
1329 pos1_out.SetXYZ(x1+cos1*s1_plus,y1+sin1*s1_plus,z1+tan1*s1_plus);
1330 pos2_out.SetXYZ(x2+cos2*s2_plus,y2+sin2*s2_plus,z2+tan2*s2_plus);
1333 doca=(pos1_out-pos2_out).Mag();
1336 DVector3 myvertex=0.5*(pos1_out+pos2_out);
1337 printf(
"Vertex position (x,y,z)=(%3.2f, %3.2f, %3.2f)\n",myvertex.x(),
1338 myvertex.y(),myvertex.z());
1352 double locAverageTime = 0.0;
1353 for(
size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
1356 locDeltaVertex = locPOCA - locParticles[loc_i]->position();
1357 locMomentum = locParticles[loc_i]->momentum();
1358 double locTime = locParticles[loc_i]->time() + locDeltaVertex.Dot(locMomentum)*locParticles[loc_i]->energy()/(29.9792458*locMomentum.Mag2());
1359 locAverageTime += locTime;
1361 return locAverageTime/(double(locParticles.size()));
1369 double locAverageTime = 0.0;
1370 for(
size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
1372 double locTime = 0.0;
1373 double locE = locParticles[loc_i]->Get_ShowerEnergy();
1374 if((locParticles[loc_i]->Get_Charge() == 0) && (locE > 0.0))
1376 double locMass = locParticles[loc_i]->Get_Mass();
1377 double locPMag =
sqrt(locE*locE - locMass*locMass);
1378 DVector3 locPosition = locParticles[loc_i]->Get_Position();
1379 DVector3 locDPosition(locPosition.X(), locPosition.Y(), locPosition.Z());
1380 DVector3 locDeltaVertex = locDPosition - locCommonVertex;
1381 locTime = locParticles[loc_i]->Get_Time() + locDeltaVertex.Mag()*locE/(29.9792458*locPMag);
1386 locDeltaVertex = locPOCA -
DVector3(locParticles[loc_i]->Get_Position().
X(),locParticles[loc_i]->Get_Position().Y(),locParticles[loc_i]->Get_Position().Z());
1387 DVector3 locMomentum(locParticles[loc_i]->Get_Momentum().
X(),locParticles[loc_i]->Get_Momentum().Y(),locParticles[loc_i]->Get_Momentum().Z());
1388 locTime = locParticles[loc_i]->Get_Time() + locDeltaVertex.Dot(locMomentum)*locParticles[loc_i]->Get_Energy()/(29.9792458*locMomentum.Mag2());
1390 locAverageTime += locTime;
1392 return locAverageTime/(double(locParticles.size()));
1401 if(locParticles.size() == 0)
1403 if(locParticles.size() == 1)
1404 return locParticles[0]->position();
1406 double locDOCA, locSmallestDOCA;
1409 locSmallestDOCA = 9.9E9;
1410 for(
int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
1412 for(
size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
1414 locDOCA =
Calc_DOCAVertex(locParticles[loc_j], locParticles[loc_k], locTempVertex);
1416 if(locDOCA < locSmallestDOCA)
1418 locSmallestDOCA = locDOCA;
1419 locVertex = locTempVertex;
1432 if(locParticles.size() == 0)
1434 if(locParticles.size() == 1)
1435 return locParticles[0]->position();
1437 double locDOCA, locSmallestDOCA;
1440 locSmallestDOCA = 9.9E9;
1441 for(
int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
1443 for(
size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
1445 locDOCA =
Calc_DOCAVertex(locParticles[loc_j], locParticles[loc_k], locTempVertex);
1446 if(locDOCA < locSmallestDOCA)
1448 locSmallestDOCA = locDOCA;
1449 locVertex = locTempVertex;
1462 if(locParticles.size() == 0)
1464 if(locParticles.size() == 1)
1465 return locParticles[0]->position();
1467 double locDOCA, locSmallestDOCA;
1470 locSmallestDOCA = 9.9E9;
1471 for(
int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
1473 for(
size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
1475 locDOCA =
Calc_DOCAVertex(locParticles[loc_j], locParticles[loc_k], locTempVertex);
1476 if(locDOCA < locSmallestDOCA)
1478 locSmallestDOCA = locDOCA;
1479 locVertex = locTempVertex;
1492 if(locParticles.size() == 0)
1495 if(locParticles.size() == 1)
1496 return DVector3(locParticles[0]->Get_Position().X(), locParticles[0]->Get_Position().Y(), locParticles[0]->Get_Position().Z());
1498 double locDOCA, locSmallestDOCA;
1501 locSmallestDOCA = 9.9E9;
1502 for(
int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
1504 for(
size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
1506 locDOCA =
Calc_DOCAVertex(locParticles[loc_j].
get(), locParticles[loc_k].
get(), locTempVertex);
1507 if(locDOCA < locSmallestDOCA)
1509 locSmallestDOCA = locDOCA;
1510 locVertex = locTempVertex;
1520 set<set<size_t> > locCombos;
1522 auto locReactionStepPIDs = locReactionStep->
Get_FinalPIDs();
1525 if(locToIncludePIDs.empty())
1527 set<size_t> locCombo;
1528 for(
size_t loc_i = 0; loc_i < locReactionStepPIDs.size(); ++loc_i)
1530 if(
int(loc_i) == locMissingParticleIndex)
1532 locCombo.insert(loc_i);
1534 locCombos.insert(locCombo);
1539 deque<deque<size_t> > locPossibilities(locToIncludePIDs.size(), deque<size_t>());
1540 deque<int> locResumeAtIndices(locToIncludePIDs.size(), 0);
1543 for(
size_t loc_i = 0; loc_i < locReactionStepPIDs.size(); ++loc_i)
1545 if(
int(loc_i) == locMissingParticleIndex)
1547 Particle_t locPID = locReactionStepPIDs[loc_i];
1550 for(
size_t loc_j = 0; loc_j < locToIncludePIDs.size(); ++loc_j)
1552 if(locToIncludePIDs[loc_j] == locPID)
1553 locPossibilities[loc_j].push_back(loc_i);
1558 int locParticleIndex = 0;
1559 deque<size_t> locComboDeque;
1562 if(locParticleIndex ==
int(locPossibilities.size()))
1564 set<size_t> locComboSet;
1565 for(
size_t loc_i = 0; loc_i < locComboDeque.size(); ++loc_i)
1566 locComboSet.insert(locComboDeque[loc_i]);
1567 locCombos.insert(locComboSet);
1569 if(!
Handle_Decursion(locParticleIndex, locComboDeque, locResumeAtIndices, locPossibilities))
1574 int& locResumeAtIndex = locResumeAtIndices[locParticleIndex];
1578 Particle_t locToIncludePID = locToIncludePIDs[locParticleIndex];
1579 for(
int loc_i = locParticleIndex - 1; loc_i >= 0; --loc_i)
1581 if(locToIncludePIDs[loc_i] == locToIncludePID)
1583 if(locResumeAtIndex < locResumeAtIndices[loc_i])
1584 locResumeAtIndex = locResumeAtIndices[loc_i];
1589 if(locResumeAtIndex >=
int(locPossibilities[locParticleIndex].
size()))
1591 if(!
Handle_Decursion(locParticleIndex, locComboDeque, locResumeAtIndices, locPossibilities))
1597 locComboDeque.push_back(locPossibilities[locParticleIndex][locResumeAtIndex]);
1609 if(locParticleIndex <
int(locResumeAtIndices.size()))
1610 locResumeAtIndices[locParticleIndex] = 0;
1613 if(locParticleIndex < 0)
1616 locComboDeque.pop_back();
1618 while(locResumeAtIndices[locParticleIndex] ==
int(locPossibilities[locParticleIndex].
size()));
1629 double locDistance = (locMeasuredX4.Vect() - locPropagateToPoint).Mag();
1634 auto locP3 = locMeasuredP4.Vect();
1635 Calc_DOCAToVertex(locP3.Unit(), locMeasuredX4.Vect(), locPropagateToPoint, locPOCA);
1636 locMeasuredX4.SetVect(locPOCA);
1637 auto locDistanceVector = locPOCA - locMeasuredX4.Vect();
1639 auto locDeltaPathLength = (locP3.Dot(locDistanceVector) > 0.0) ? -1.0*locDistanceVector.Mag() : locDistanceVector.Mag();
1640 locMeasuredX4.SetT(locMeasuredX4.T() - locDeltaPathLength/(locMeasuredP4.Beta()*
SPEED_OF_LIGHT));
1641 return locDeltaPathLength;
1645 double locTotalDeltaPathLength = 0.0;
1648 if(fabs(locMeasuredP4.Pz()) > 0.0)
1650 locTotalDeltaPathLength += (locPropagateToPoint.Z() - locMeasuredX4.Z())*locMeasuredP4.P()/locMeasuredP4.Pz();
1651 Propagate_Track(locTotalDeltaPathLength, locCharge, locTempPropagatedPosition, locTempPropagatedMomentum, NULL);
1655 locTotalDeltaPathLength +=
Calc_PathLength_Step(locCharge, locPropagateToPoint, locTempPropagatedPosition, locTempPropagatedMomentum);
1658 locTotalDeltaPathLength +=
Calc_PathLength_FineGrained(locCharge, locPropagateToPoint, locTempPropagatedPosition.Vect(), locTempPropagatedMomentum.Vect());
1662 Propagate_Track(locTotalDeltaPathLength, locCharge, locMeasuredX4, locMeasuredP4, locCovarianceMatrix);
1665 return -1.0*locTotalDeltaPathLength;
1673 double locA = -0.00299792458*(double(locCharge))*locBField.Mag();
1675 double locPMag = locMeasuredP4.P();
1676 double locPathLengthOneRotation = 2.0*TMath::Pi()*locPMag/locA;
1677 double locStepDistance = locPathLengthOneRotation/12.0;
1680 double locStepDeltaZ = locMeasuredP4.Pz()*locStepDistance/locPMag;
1681 double locTargetDeltaZ = 1.0;
1682 if(fabs(locStepDeltaZ) > 5.0)
1683 locStepDistance = locTargetDeltaZ*locPMag/locMeasuredP4.Pz();
1686 double locStepDirection = 1.0;
1687 double locTotalDeltaPathLength = 0.0;
1688 double locDeltaX3 = (locMeasuredX4.Vect() - locPropagateToPoint).Mag();
1689 double locRhoS = locStepDistance*locA/locPMag;
1690 double locRhoSCutOff = 0.1;
1691 while(locRhoS > locRhoSCutOff)
1693 double locDeltaPathLength = locStepDirection*locStepDistance;
1694 locTotalDeltaPathLength += locDeltaPathLength;
1696 Propagate_Track(locDeltaPathLength, locCharge, locMeasuredX4, locMeasuredP4, NULL);
1698 double locNewDeltaX3 = (locMeasuredX4.Vect() - locPropagateToPoint).Mag();
1699 bool locGettingCloserFlag = (locNewDeltaX3 < locDeltaX3);
1700 locDeltaX3 = locNewDeltaX3;
1702 if(locGettingCloserFlag)
1706 locStepDirection *= -1.0;
1707 locStepDistance /= 2.0;
1708 locRhoS = locStepDistance*locA/locPMag;
1711 if(locRhoS < locRhoSCutOff)
1713 locTotalDeltaPathLength -= locDeltaPathLength;
1714 Propagate_Track(-1.0*locDeltaPathLength, locCharge, locMeasuredX4, locMeasuredP4, NULL);
1718 return locTotalDeltaPathLength;
1739 double locA = -0.00299792458*(double(locCharge))*locBField.Mag();
1740 double locPMag = locMeasuredMomentum.Mag();
1741 double locRho = locA/locPMag;
1743 DVector3 locDeltaX3 = locMeasuredPosition - locPropagateToPoint;
1744 double locF = locMeasuredMomentum.Pz()/locPMag;
1745 double locC = locDeltaX3.Z();
1746 double locD = locRho*(locMeasuredMomentum.Px()*locDeltaX3.X() + locMeasuredMomentum.Py()*locDeltaX3.Y())/locA;
1747 double locPerpPSq = locMeasuredMomentum.Px()*locMeasuredMomentum.Px() + locMeasuredMomentum.Py()*locMeasuredMomentum.Py();
1748 double locE = locRho*(locPerpPSq/locA - locMeasuredMomentum.Py()*locDeltaX3.X() + locMeasuredMomentum.Px()*locDeltaX3.Y())/locA;
1749 return -1.0*(locF*locC + locD)/(locF*locF + locE*locRho);
1757 if(!(locBField.Mag() > 0.0))
1760 double locA = -0.00299792458*(double(locCharge))*locBField.Mag();
1762 double locPMag = locP4.P();
1763 double locRhoS = locDeltaPathLength*locA/locPMag;
1766 locDeltaX4.SetX(locP4.Px()*
sin(locRhoS)/locA - locP4.Py()*(1.0 - cos(locRhoS))/locA);
1767 locDeltaX4.SetY(locP4.Py()*
sin(locRhoS)/locA + locP4.Px()*(1.0 - cos(locRhoS))/locA);
1768 locDeltaX4.SetZ(locP4.Pz()*locDeltaPathLength/locPMag);
1769 locDeltaX4.SetT(locDeltaPathLength/(locP4.Beta()*
SPEED_OF_LIGHT));
1770 locX4 += locDeltaX4;
1772 if(locCovarianceMatrix == NULL)
1774 locP4.SetVect(locP4.Vect() - locDeltaX4.Vect().Cross(locA*locH));
1779 TMatrixF locTransformMatrix(7, 7);
1780 locTransformMatrix.Zero();
1782 double locPCubed = locPMag*locPMag*locPMag;
1783 double locSOverP3 = locDeltaPathLength/locPCubed;
1785 double locSOverP3_PxPx = locSOverP3*locP4.Px()*locP4.Px();
1786 double locSOverP3_PxPy = locSOverP3*locP4.Px()*locP4.Py();
1787 double locSOverP3_PxPz = locSOverP3*locP4.Px()*locP4.Pz();
1789 double locSOverP3_PyPy = locSOverP3*locP4.Py()*locP4.Py();
1790 double locSOverP3_PyPz = locSOverP3*locP4.Py()*locP4.Pz();
1791 double locSOverP3_PzPz = locSOverP3*locP4.Pz()*locP4.Pz();
1794 locTransformMatrix(0, 0) = (1.0 + locA*locSOverP3_PxPy)*cos(locRhoS) + locA*locSOverP3_PxPx*
sin(locRhoS);
1795 locTransformMatrix(0, 1) = (-1.0 + locA*locSOverP3_PxPy)*
sin(locRhoS) + locA*locSOverP3_PyPy*cos(locRhoS);
1796 locTransformMatrix(0, 2) = locA*locSOverP3_PyPz*cos(locRhoS) + locA*locSOverP3_PxPz*
sin(locRhoS);
1799 locTransformMatrix(1, 0) = -1.0*locA*locSOverP3_PxPx*cos(locRhoS) + (1.0 + locA*locSOverP3_PxPy)*
sin(locRhoS);
1800 locTransformMatrix(1, 1) = (1.0 - locA*locSOverP3_PxPy)*cos(locRhoS) + locA*locSOverP3_PyPy*
sin(locRhoS);
1801 locTransformMatrix(1, 2) = locA*locSOverP3_PyPz*
sin(locRhoS) - locA*locSOverP3_PxPz*cos(locRhoS);
1804 locTransformMatrix(2, 2) = 1.0;
1807 locTransformMatrix(3, 0) = (1.0/locA + locSOverP3_PxPy)*
sin(locRhoS) - locSOverP3_PxPx*cos(locRhoS);
1808 locTransformMatrix(3, 1) = -1.0/locA + (1.0/locA - locSOverP3_PxPy)*cos(locRhoS) + locSOverP3_PyPy*
sin(locRhoS);
1809 locTransformMatrix(3, 2) = locSOverP3_PyPz*
sin(locRhoS) - locSOverP3_PxPz*cos(locRhoS);
1810 locTransformMatrix(3, 3) = 1.0;
1813 locTransformMatrix(4, 0) = 1.0/locA - (1.0/locA + locSOverP3_PxPy)*cos(locRhoS) + locSOverP3_PxPx*
sin(locRhoS);
1814 locTransformMatrix(4, 1) = (1.0/locA - locSOverP3_PxPy)*
sin(locRhoS) - locSOverP3_PyPy*cos(locRhoS);
1815 locTransformMatrix(4, 2) = -1.0*locSOverP3_PyPz*cos(locRhoS) - locSOverP3_PxPz*
sin(locRhoS);
1816 locTransformMatrix(4, 4) = 1.0;
1819 locTransformMatrix(5, 0) = -1.0*locSOverP3_PxPz;
1820 locTransformMatrix(5, 1) = -1.0*locSOverP3_PyPz;
1821 locTransformMatrix(5, 2) = locDeltaPathLength/locPMag - locSOverP3_PzPz;
1822 locTransformMatrix(5, 5) = 1.0;
1825 locTransformMatrix(6, 0) = -1.0*locP4.M2()*locSOverP3*locP4.Px()/(
SPEED_OF_LIGHT*locP4.E());
1826 locTransformMatrix(6, 1) = -1.0*locP4.M2()*locSOverP3*locP4.Py()/(
SPEED_OF_LIGHT*locP4.E());
1827 locTransformMatrix(6, 2) = -1.0*locP4.M2()*locSOverP3*locP4.Pz()/(
SPEED_OF_LIGHT*locP4.E());
1828 locTransformMatrix(6, 6) = 1.0;
1831 locCovarianceMatrix->Similarity(locTransformMatrix);
1834 locP4.SetVect(locP4.Vect() - locDeltaX4.Vect().Cross(locA*locH));
vector< const DKinematicData * > Get_FinalParticles_Measured(void) const
double dMinDistanceForStraightTrack
vector< const JObject * > Get_FinalParticle_SourceObjects(Charge_t locCharge=d_AllCharges) const
char Get_Charge(void) const
void Add_ReactionStep(const DReactionStep *locReactionStep)
int Calc_Momentum_UnusedTracks(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, double &locSumPMag_UnusedTracks, TVector3 &locSumP3_UnusedTracks) const
pair< int, int > Get_InitialParticleDecayFromIndices(const DReaction *locReaction, int locStepIndex)
size_t Get_NumParticleComboSteps(void) const
bool Check_ChannelEquality(const DReaction *lhs, const DReaction *rhs, bool locSameOrderFlag=true, bool locRightSubsetOfLeftFlag=false)
DAnalysisUtilities(JEventLoop *loop)
void Get_UnusedTimeBasedTracks(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, vector< const DTrackTimeBased * > &locUnusedTimeBasedTracks) const
static unsigned short int IsResonance(Particle_t p)
const DChargedTrackHypothesis * Get_BestTrackingFOM(void) const
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
bool Get_IsBFieldNearBeamline(void) const
const DVector3 & position(void) const
bool Check_IsBDTSignalEvent(JEventLoop *locEventLoop, const DReaction *locReaction, bool locExclusiveMatchFlag, bool locIncludeDecayingToReactionFlag) const
const DParticleCombo * Build_ThrownCombo(JEventLoop *locEventLoop)
void Clear_ReactionSteps(void)
Particle_t Get_TargetPID(void) const
double Calc_Energy_UnusedShowers(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo) const
string dShowerSelectionTag
set< set< size_t > > Build_IndexCombos(const DReactionStep *locReactionStep, deque< Particle_t > locToIncludePIDs) const
double Calc_DOCAToVertex(const DVector3 &locUnitDir, const DVector3 &locPosition, const DVector3 &locVertex) const
const DKinematicData * Get_InitialParticle_Measured(void) const
TLorentzVector DLorentzVector
const DEventRFBunch * Get_EventRFBunch(void) const
void Recycle_Reaction(DReaction *locReaction)
static int PDGtype(Particle_t p)
double Calc_DOCAVertex(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2, DVector3 &locDOCAPoint) const
DVector3 Get_BField(const DVector3 &locPosition) const
double Calc_DOCA(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2) const
bool Handle_Decursion(int &locParticleIndex, deque< size_t > &locComboDeque, deque< int > &locResumeAtIndices, deque< deque< size_t > > &locPossibilities) const
DLorentzVector dSpacetimeVertex
DLorentzVector Calc_FinalStateP4(const DReaction *locReaction, const DParticleCombo *locParticleCombo, size_t locStepIndex, bool locUseKinFitDataFlag) const
void Replace_DecayingParticleWithProducts(deque< pair< const DMCThrown *, deque< const DMCThrown * > > > &locThrownSteps, size_t locStepIndex) const
Particle_t Get_InitialPID(void) const
void Get_UnusedChargedTracks(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, vector< const DChargedTrack * > &locUnusedChargedTracks) const
void Get_UnusedTrackCandidates(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, vector< const DTrackCandidate * > &locUnusedTrackCandidates) const
void Add_FinalParticleID(Particle_t locPID, bool locIsMissingFlag=false)
DVector3 Calc_CrudeVertex(const vector< const DKinematicData * > &locParticles) const
DLorentzVector Calc_MissingP4(const DReaction *locReaction, const DParticleCombo *locParticleCombo, bool locUseKinFitDataFlag) const
DGeometry * GetDGeometry(unsigned int run_number)
double Calc_PathLength_FineGrained(int locCharge, const DVector3 &locPropagateToPoint, DVector3 locMeasuredPosition, DVector3 locMeasuredMomentum) const
double Calc_CrudeTime(const vector< const DKinematicData * > &locParticles, const DVector3 &locCommonVertex) const
static unsigned short int IsFixedMass(Particle_t p)
vector< const DKinematicData * > Get_FinalParticles(void) const
Particle_t Get_FinalPID(size_t locIndex) const
bool Are_ThrownPIDsSameAsDesired(JEventLoop *locEventLoop, const deque< Particle_t > &locDesiredPIDs, Particle_t locMissingPID=Unknown) const
void Get_UnusedWireBasedTracks(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, vector< const DTrackWireBased * > &locUnusedWireBasedTracks) const
static int Is_FinalStateParticle(Particle_t locPID)
double charge(void) const
int Get_MissingParticleIndex(void) const
void Get_ThrownParticleSteps(JEventLoop *locEventLoop, deque< pair< const DMCThrown *, deque< const DMCThrown * > > > &locThrownSteps) const
DLorentzVector lorentzMomentum(void) const
static double ParticleMass(Particle_t p)
const JObject * Get_FinalParticle_SourceObject(size_t locFinalParticleIndex) const
virtual double GetBz(double x, double y, double z) const =0
const DKinematicData * Get_InitialParticle(void) const
size_t Get_NumFinalPIDs(void) const
void Set_InitialParticleID(Particle_t locPID, bool locIsMissingFlag=false)
DParticleComboCreator * dParticleComboCreator
const DVector3 & momentum(void) const
vector< Particle_t > Get_FinalPIDs(bool locIncludeMissingFlag=true, Charge_t locCharge=d_AllCharges, bool locIncludeDuplicatesFlag=true) const
double Propagate_Track(int locCharge, const DVector3 &locPropagateToPoint, DLorentzVector &locMeasuredX4, DLorentzVector &locMeasuredP4, TMatrixFSym *locCovarianceMatrix) const
DReaction * Build_ThrownReaction(JEventLoop *locEventLoop, deque< pair< const DMCThrown *, deque< const DMCThrown * > > > &locThrownSteps)
const DReactionStep * Get_ReactionStep(size_t locStepIndex) const
void Get_UnusedNeutralParticles(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, vector< const DNeutralParticle * > &locUnusedNeutralParticles) const
int Get_DecayStepIndex(const DReaction *locReaction, size_t locStepIndex, size_t locParticleIndex)
const DKinFitResults * Get_KinFitResults(void) const
TVector3 Get_Position(void) const
void Get_UnusedNeutralShowers(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo, vector< const DNeutralShower * > &locUnusedNeutralShowers) const
shared_ptr< const TMatrixFSym > errorMatrix(void) const
bool Check_ThrownsMatchReaction(JEventLoop *locEventLoop, const DReaction *locReaction, bool locExclusiveMatchFlag) const
const DParticleComboStep * Get_ParticleComboStep(size_t locStepIndex) const
size_t Get_NumReactionSteps(void) const
const DMagneticFieldMap * dMagneticFieldMap
printf("string=%s", string)
string dTrackSelectionTag
const DParticleID * dPIDAlgorithm
void Set_TargetParticleID(Particle_t locPID)
bool GetTargetZ(double &z_target) const
z-location of center of target
Particle_t PID(void) const
TVector3 Get_Momentum(void) const
double Calc_PathLength_Step(int locCharge, const DVector3 &locPropagateToPoint, DLorentzVector &locMeasuredX4, DLorentzVector &locMeasuredP4) const
TMatrixFSym Calc_MissingP3Covariance(const DReaction *locReaction, const DParticleCombo *locParticleCombo) const