36 gPARMS->SetDefaultParameter(
"MCMATCH:DEBUG_LEVEL",
dDebugLevel);
56 VT_TRACER(
"DMCThrownMatching_factory::evnt()");
59 vector<const DMCThrown*> locMCThrowns;
60 locEventLoop->Get(locMCThrowns);
62 if(locMCThrowns.empty())
66 locEventLoop->Get(locMCThrowns,
"FinalState");
68 vector<const DMCThrown*> locMCThrowns_Charged;
69 vector<const DMCThrown*> locMCThrowns_Neutral;
71 for(
size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
74 locMCThrowns_Neutral.push_back(locMCThrowns[loc_i]);
76 locMCThrowns_Charged.push_back(locMCThrowns[loc_i]);
79 cout <<
"input #thrown, ok charged # thrown, ok neutral # thrown = " << locMCThrowns.size() <<
", " << locMCThrowns_Charged.size() <<
", " << locMCThrowns_Neutral.size() << endl;
81 vector<const DChargedTrack*> locChargedTracks;
82 locEventLoop->Get(locChargedTracks);
84 vector<const DNeutralParticle*> locNeutralParticles;
85 locEventLoop->Get(locNeutralParticles);
87 vector<const DChargedTrackHypothesis*> locChargedTrackHypotheses;
88 locEventLoop->Get(locChargedTrackHypotheses);
90 vector<const DNeutralParticleHypothesis*> locNeutralParticleHypotheses;
91 locEventLoop->Get(locNeutralParticleHypotheses);
109 double locMatchFOM = 0.0;
110 cout <<
"Charged Track Matching Summary:" << endl;
111 for(
size_t loc_i = 0; loc_i < locChargedTrackHypotheses.size(); ++loc_i)
113 double locP = locChargedTrackHypotheses[loc_i]->momentum().Mag();
114 double locTheta = locChargedTrackHypotheses[loc_i]->momentum().Theta()*180.0/TMath::Pi();
115 double locPhi = locChargedTrackHypotheses[loc_i]->momentum().Phi()*180.0/TMath::Pi();
116 Particle_t locPID = locChargedTrackHypotheses[loc_i]->PID();
117 cout <<
"charged info: " << locChargedTrackHypotheses[loc_i]->Get_TrackTimeBased()->candidateid <<
", " <<
ParticleType(locPID) <<
", " << locP <<
", " << locTheta <<
", " << locPhi << endl;
118 cout <<
"matched info: ";
120 if(locMCThrown == NULL)
122 cout <<
"NO MATCH." << endl;
125 locP = locMCThrown->
momentum().Mag();
126 locTheta = locMCThrown->
momentum().Theta()*180.0/TMath::Pi();
127 locPhi = locMCThrown->
momentum().Phi()*180.0/TMath::Pi();
128 locPID = locMCThrown->
PID();
129 cout <<
ParticleType(locPID) <<
", " << locP <<
", " << locTheta <<
", " << locPhi <<
", " << locMatchFOM << endl;
131 cout <<
"Neutral Particle Matching Summary:" << endl;
132 for(
size_t loc_i = 0; loc_i < locNeutralParticleHypotheses.size(); ++loc_i)
134 double locP = locNeutralParticleHypotheses[loc_i]->momentum().Mag();
135 double locTheta = locNeutralParticleHypotheses[loc_i]->momentum().Theta()*180.0/TMath::Pi();
136 double locPhi = locNeutralParticleHypotheses[loc_i]->momentum().Phi()*180.0/TMath::Pi();
137 Particle_t locPID = locNeutralParticleHypotheses[loc_i]->PID();
138 cout <<
"neutral info: " << locNeutralParticleHypotheses[loc_i] <<
", " <<
ParticleType(locPID) <<
", " << locP <<
", " << locTheta <<
", " << locPhi << endl;
139 cout <<
"matched info: ";
141 if(locMCThrown == NULL)
143 cout <<
"NO MATCH." << endl;
146 locP = locMCThrown->
momentum().Mag();
147 locTheta = locMCThrown->
momentum().Theta()*180.0/TMath::Pi();
148 locPhi = locMCThrown->
momentum().Phi()*180.0/TMath::Pi();
149 locPID = locMCThrown->
PID();
150 cout <<
ParticleType(locPID) <<
", " << locP <<
", " << locTheta <<
", " << locPhi <<
", " << locMatchFOM << endl;
152 cout <<
"Unmatched Charged DMCThrowns:" << endl;
153 for(
size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
159 double locP = locMCThrowns[loc_i]->momentum().Mag();
160 double locTheta = locMCThrowns[loc_i]->momentum().Theta()*180.0/TMath::Pi();
161 double locPhi = locMCThrowns[loc_i]->momentum().Phi()*180.0/TMath::Pi();
162 Particle_t locPID = locMCThrowns[loc_i]->PID();
163 cout <<
"thrown info: " <<
ParticleType(locPID) <<
", " << locP <<
", " << locTheta <<
", " << locPhi << endl;
167 _data.push_back(locMCThrownMatching);
174 vector<const DBeamPhoton*> locBeamPhotons;
175 locEventLoop->Get(locBeamPhotons);
177 vector<const DBeamPhoton*> locBeamPhotons_MCGEN;
178 locEventLoop->Get(locBeamPhotons_MCGEN,
"MCGEN");
180 vector<const DBeamPhoton*> locBeamPhotons_TAGGEDMCGEN;
181 locEventLoop->Get(locBeamPhotons_TAGGEDMCGEN,
"TAGGEDMCGEN");
183 vector<const DBeamPhoton*> locBeamPhotons_TRUTH;
184 locEventLoop->Get(locBeamPhotons_TRUTH,
"TRUTH");
187 locMCThrownMatching->
Set_MCGENBeamPhoton(locBeamPhotons_MCGEN.empty() ?
nullptr : locBeamPhotons_MCGEN[0]);
188 locMCThrownMatching->
Set_TaggedMCGENBeamPhoton(locBeamPhotons_TAGGEDMCGEN.empty() ?
nullptr : locBeamPhotons_TAGGEDMCGEN[0]);
191 map<const DBeamPhoton*, deque<pair<const DBeamPhoton*, double> > > locPossibleMatches;
192 for(
size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
194 for(
size_t loc_j = 0; loc_j < locBeamPhotons_TRUTH.size(); ++loc_j)
197 locBeamPhotons[loc_i]->GetSingle(locTAGMHit);
199 locBeamPhotons[loc_i]->GetSingle(locTAGHHit);
201 const DTAGMHit* locTAGMHit_TRUTH = NULL;
202 locBeamPhotons_TRUTH[loc_j]->GetSingle(locTAGMHit_TRUTH);
203 const DTAGHHit* locTAGHHit_TRUTH = NULL;
204 locBeamPhotons_TRUTH[loc_j]->GetSingle(locTAGHHit_TRUTH);
206 if(locTAGMHit != NULL)
208 if(locTAGMHit_TRUTH == NULL)
210 if((locTAGMHit->
row != locTAGMHit_TRUTH->
row) || (locTAGMHit->
column != locTAGMHit_TRUTH->
column))
215 if(locTAGHHit_TRUTH == NULL)
221 double locDeltaT = fabs(locBeamPhotons[loc_i]->time() - locBeamPhotons_TRUTH[loc_j]->time());
222 pair<const DBeamPhoton*, double> locMatchPair(locBeamPhotons_TRUTH[loc_j], locDeltaT);
223 locPossibleMatches[locBeamPhotons[loc_i]].push_back(locMatchPair);
227 map<const DBeamPhoton*, const DBeamPhoton*> locBeamPhotonToTruthMap;
228 map<const DBeamPhoton*, const DBeamPhoton*> locBeamTruthToPhotonMap;
230 map<const DBeamPhoton*, deque<pair<const DBeamPhoton*, double> > >::iterator locIterator = locPossibleMatches.begin();
231 for(; locIterator != locPossibleMatches.end(); ++locIterator)
233 deque<pair<const DBeamPhoton*, double> > locPairs = locIterator->second;
234 double locBestDeltaT = 9.9E9;
236 for(
size_t loc_i = 0; loc_i < locPairs.size(); ++loc_i)
238 double locDeltaT = locPairs[loc_i].second;
239 if(locDeltaT >= locBestDeltaT)
241 locBestDeltaT = locDeltaT;
242 locBestBeamPhoton = locPairs[loc_i].first;
245 if((
dDebugLevel > 20) && (locBestBeamPhoton != NULL))
246 cout <<
"Photon match: delta-E, delta-t = " << fabs(locIterator->first->energy() - locBestBeamPhoton->
energy()) <<
", " << locBestDeltaT << endl;
247 locBeamPhotonToTruthMap[locIterator->first] = locBestBeamPhoton;
248 locBeamTruthToPhotonMap[locBestBeamPhoton] = locIterator->first;
257 vector<const DFCALTruthShower*> locFCALTruthShowers;
259 locEventLoop->Get(locFCALTruthShowers);
261 vector<const DFCALShower*> locFCALShowers;
263 locEventLoop->Get(locFCALShowers);
265 map<const DFCALShower*, pair<const DFCALTruthShower*, double> > locFCALShowerToTruthMap;
266 map<const DFCALTruthShower*, pair<const DFCALShower*, double> > locFCALTruthToShowerMap;
268 size_t locBestFCALShowerIndex = 0, locBestFCALTruthShowerIndex = 0;
269 double locBestMatchDistance = 0.0, locMatchDistance;
270 while((!locFCALShowers.empty()) && (!locFCALTruthShowers.empty()))
272 bool locMatchFoundFlag =
false;
273 double locLargestEnergy = 0.0;
274 for(
size_t loc_i = 0; loc_i < locFCALTruthShowers.size(); ++loc_i)
276 locFCALTruthShower = locFCALTruthShowers[loc_i];
277 for(
size_t loc_j = 0; loc_j < locFCALShowers.size(); ++loc_j)
279 locFCALShower = locFCALShowers[loc_j];
284 double locDeltaZ = locFCALShowerPosition.Z() - locFCALTruthShower->
z();
285 DVector3 locTrueMomentum(locFCALTruthShower->
px(), locFCALTruthShower->
py(), locFCALTruthShower->
pz());
286 double locDeltaPathLength = locDeltaZ/cos(locTrueMomentum.Theta());
287 double locTrueX = locFCALTruthShower->
x() + locDeltaPathLength*
sin(locTrueMomentum.Theta())*cos(locTrueMomentum.Phi());
288 double locTrueY = locFCALTruthShower->
y() + locDeltaPathLength*
sin(locTrueMomentum.Theta())*
sin(locTrueMomentum.Phi());
289 DVector3 locFCALTruthShowerPosition(locTrueX, locTrueY, locFCALShowerPosition.Z());
291 locMatchDistance = (locFCALShowerPosition - locFCALTruthShowerPosition).Mag();
296 if(locFCALShower->
getEnergy() > locLargestEnergy)
298 locBestMatchDistance = locMatchDistance;
299 locMatchFoundFlag =
true;
300 locLargestEnergy = locFCALShower->
getEnergy();
301 locBestFCALTruthShowerIndex = loc_i;
302 locBestFCALShowerIndex = loc_j;
307 if(!locMatchFoundFlag)
310 locFCALTruthShower = locFCALTruthShowers[locBestFCALTruthShowerIndex];
311 locFCALShower = locFCALShowers[locBestFCALShowerIndex];
313 locFCALShowerToTruthMap[locFCALShower] = pair<const DFCALTruthShower*, double>(locFCALTruthShower, locBestMatchDistance);
314 locFCALTruthToShowerMap[locFCALTruthShower] = pair<const DFCALShower*, double>(locFCALShower, locBestMatchDistance);
316 locFCALShowers.erase(locFCALShowers.begin() + locBestFCALShowerIndex);
317 locFCALTruthShowers.erase(locFCALTruthShowers.begin() + locBestFCALTruthShowerIndex);
326 vector<const DBCALTruthShower*> locBCALTruthShowers;
328 locEventLoop->Get(locBCALTruthShowers);
330 vector<const DBCALShower*> locBCALShowers;
332 locEventLoop->Get(locBCALShowers);
334 map<const DBCALShower*, pair<const DBCALTruthShower*, double> > locBCALShowerToTruthMap;
335 map<const DBCALTruthShower*, pair<const DBCALShower*, double> > locBCALTruthToShowerMap;
337 vector<const DMCThrown*> locMCThrowns;
338 locEventLoop->Get(locMCThrowns);
341 map<const DBCALTruthShower*, const DMCThrown*> locBCALTruthToMCThrownMap;
342 for(
size_t loc_i = 0; loc_i < locBCALTruthShowers.size(); ++loc_i)
344 for(
size_t loc_j = 0; loc_j < locMCThrowns.size(); ++loc_j)
346 if(locBCALTruthShowers[loc_i]->
track != locMCThrowns[loc_j]->myid)
348 locBCALTruthToMCThrownMap[locBCALTruthShowers[loc_i]] = locMCThrowns[loc_j];
353 size_t locBestBCALShowerIndex = 0, locBestBCALTruthShowerIndex = 0;
354 while((!locBCALShowers.empty()) && (!locBCALTruthShowers.empty()))
356 bool locMatchFoundFlag =
false;
357 double locLargestEnergy = 0.0;
358 double locBestUnitCircleArcLength = 0.0;
359 for(
size_t loc_i = 0; loc_i < locBCALTruthShowers.size(); ++loc_i)
361 locBCALTruthShower = locBCALTruthShowers[loc_i];
362 for(
size_t loc_j = 0; loc_j < locBCALShowers.size(); ++loc_j)
364 locBCALShower = locBCALShowers[loc_j];
366 DVector3 locBCALShowerPosition(locBCALShower->
x, locBCALShower->
y, locBCALShower->
z);
367 DVector3 locBCALTruthShowerPosition(locBCALTruthShower->
r*cos(locBCALTruthShower->
phi), locBCALTruthShower->
r*
sin(locBCALTruthShower->
phi), locBCALTruthShower->
z);
370 map<const DBCALTruthShower*, const DMCThrown*>::const_iterator locIterator = locBCALTruthToMCThrownMap.find(locBCALTruthShower);
372 if(locIterator != locBCALTruthToMCThrownMap.end())
373 locProductionVertex = locIterator->second->position();
374 double locTrueTheta = (locBCALTruthShowerPosition - locProductionVertex).Theta();
375 double locReconstructedTheta = (locBCALShowerPosition - locProductionVertex).Theta();
377 double locDeltaPhi = locBCALShowerPosition.Phi() - locBCALTruthShower->
phi;
378 while(locDeltaPhi > TMath::Pi())
379 locDeltaPhi -= 2.0*TMath::Pi();
380 while(locDeltaPhi < -1.0*TMath::Pi())
381 locDeltaPhi += 2.0*TMath::Pi();
383 double locUnitCircleArcLength = acos(
sin(locTrueTheta)*
sin(locReconstructedTheta)*cos(locDeltaPhi) + cos(locTrueTheta)*cos(locReconstructedTheta));
398 if(locBCALShower->
E > locLargestEnergy)
400 locBestUnitCircleArcLength = locUnitCircleArcLength;
401 locMatchFoundFlag =
true;
402 locLargestEnergy = locBCALShower->
E;
403 locBestBCALTruthShowerIndex = loc_i;
404 locBestBCALShowerIndex = loc_j;
409 if(!locMatchFoundFlag)
412 locBCALTruthShower = locBCALTruthShowers[locBestBCALTruthShowerIndex];
413 locBCALShower = locBCALShowers[locBestBCALShowerIndex];
415 locBCALShowerToTruthMap[locBCALShower] = pair<const DBCALTruthShower*, double>(locBCALTruthShower, locBestUnitCircleArcLength);
416 locBCALTruthToShowerMap[locBCALTruthShower] = pair<const DBCALShower*, double>(locBCALShower, locBestUnitCircleArcLength);
418 locBCALShowers.erase(locBCALShowers.begin() + locBestBCALShowerIndex);
419 locBCALTruthShowers.erase(locBCALTruthShowers.begin() + locBestBCALTruthShowerIndex);
428 vector<const DTOFTruth*> locTOFTruths;
430 locEventLoop->Get(locTOFTruths);
432 vector<const DTOFPoint*> locTOFPoints;
434 locEventLoop->Get(locTOFPoints);
436 map<const DTOFPoint*, pair<const DTOFTruth*, double> > locTOFPointToTruthMap;
437 map<const DTOFTruth*, pair<const DTOFPoint*, double> > locTOFTruthToPointMap;
439 size_t locBestTOFPointIndex = 0, locBestTOFTruthIndex = 0;
440 double locBestMatchDistance = 0.0, locMatchDistance;
441 while((!locTOFPoints.empty()) && (!locTOFTruths.empty()))
443 bool locMatchFoundFlag =
false;
444 double locLargestEnergy = 0.0;
445 for(
size_t loc_i = 0; loc_i < locTOFTruths.size(); ++loc_i)
447 locTOFTruth = locTOFTruths[loc_i];
448 for(
size_t loc_j = 0; loc_j < locTOFPoints.size(); ++loc_j)
450 locTOFPoint = locTOFPoints[loc_j];
454 double locDeltaZ = locTOFPointPosition.Z() - locTOFTruth->
z;
455 DVector3 locTrueMomentum(locTOFTruth->
px, locTOFTruth->
py, locTOFTruth->
pz);
456 double locDeltaPathLength = locDeltaZ/cos(locTrueMomentum.Theta());
457 double locTrueX = locTOFTruth->
x + locDeltaPathLength*
sin(locTrueMomentum.Theta())*cos(locTrueMomentum.Phi());
458 double locTrueY = locTOFTruth->
y + locDeltaPathLength*
sin(locTrueMomentum.Theta())*
sin(locTrueMomentum.Phi());
459 DVector3 locTOFTruthPosition(locTrueX, locTrueY, locTOFPointPosition.Z());
461 locMatchDistance = (locTOFTruthPosition - locTOFPointPosition).Mag();
466 if(locTOFTruth->
E > locLargestEnergy)
468 locBestMatchDistance = locMatchDistance;
469 locMatchFoundFlag =
true;
470 locLargestEnergy = locTOFTruth->
E;
471 locBestTOFTruthIndex = loc_i;
472 locBestTOFPointIndex = loc_j;
477 if(!locMatchFoundFlag)
480 locTOFTruth = locTOFTruths[locBestTOFTruthIndex];
481 locTOFPoint = locTOFPoints[locBestTOFPointIndex];
483 locTOFPointToTruthMap[locTOFPoint] = pair<const DTOFTruth*, double>(locTOFTruth, locBestMatchDistance);
484 locTOFTruthToPointMap[locTOFTruth] = pair<const DTOFPoint*, double>(locTOFPoint, locBestMatchDistance);
486 locTOFPoints.erase(locTOFPoints.begin() + locBestTOFPointIndex);
487 locTOFTruths.erase(locTOFTruths.begin() + locBestTOFTruthIndex);
497 map<const DChargedTrackHypothesis*, pair<const DMCThrown*, double> > locChargedHypoToThrownMap;
498 map<const DMCThrown*, pair<deque<const DChargedTrackHypothesis*>,
double> > locThrownToChargedHypoMap;
502 map<const DChargedTrack*, pair<const DMCThrown*, double> > locChargedToThrownMap;
503 map<const DMCThrown*, pair<const DChargedTrack*, double> > locThrownToChargedMap;
505 map<const DMCThrown*, pair<deque<const DChargedTrackHypothesis*>,
double> >::iterator locIterator;
506 for(locIterator = locThrownToChargedHypoMap.begin(); locIterator != locThrownToChargedHypoMap.end(); ++locIterator)
508 for(
size_t loc_j = 0; loc_j < locChargedTracks.size(); ++loc_j)
510 if(locChargedTracks[loc_j]->candidateid == (locIterator->second.first)[0]->Get_TrackTimeBased()->candidateid)
512 locChargedToThrownMap[locChargedTracks[loc_j]] = pair<const DMCThrown*, double>(locIterator->first, locIterator->second.second);
513 locThrownToChargedMap[locIterator->first] = pair<const DChargedTrack*, double>(locChargedTracks[loc_j], locIterator->second.second);
525 map<const DChargedTrackHypothesis*, pair<const DMCThrown*, double> > locChargedToThrownMap;
526 map<const DMCThrown*, pair<deque<const DChargedTrackHypothesis*>,
double> > locThrownToChargedMap;
528 map<int, const DMCThrown*> locMyIDToThrownMap;
529 for(
size_t loc_i = 0; loc_i < locMCThrownVector.size(); ++loc_i)
530 locMyIDToThrownMap[locMCThrownVector[loc_i]->myid] = locMCThrownVector[loc_i];
533 cout <<
"START IT!" << endl;
536 set<pair<double, pair<const DMCThrown*, const DChargedTrackHypothesis*> > > locParticleMatches;
537 for(
size_t loc_j = 0; loc_j < locChargedTrackHypothesisVector.size(); ++loc_j)
542 if(locMyIDToThrownMap.find(locTrackTimeBased->dMCThrownMatchMyID) == locMyIDToThrownMap.end())
545 const DMCThrown* locMCThrown = locMyIDToThrownMap[locTrackTimeBased->dMCThrownMatchMyID];
547 double locHitFraction = locTrackTimeBased->dNumHitsMatchedToThrown/locNumTrackHits;
552 double locMatchFOM = locTrackTimeBased->dNumHitsMatchedToThrown * locHitFraction;
556 cout <<
"MATCHING: MCTHROWN: ";
558 cout <<
"MATCHING: CHARGEDHYPO: ";
559 cout <<
ParticleType(locChargedTrackHypothesis->
PID()) <<
", " << locChargedTrackHypothesis->
momentum().Mag() <<
", " << locChargedTrackHypothesis->
momentum().Theta()*180.0/TMath::Pi() <<
", " << locChargedTrackHypothesis->
momentum().Phi()*180.0/TMath::Pi() << endl;
563 pair<const DMCThrown*, const DChargedTrackHypothesis*> locTrackPair(locMCThrown, locChargedTrackHypothesis);
564 pair<double, pair<const DMCThrown*, const DChargedTrackHypothesis*> > locMatchPair(locMatchFOM, locTrackPair);
565 locParticleMatches.insert(locMatchPair);
567 if(locParticleMatches.empty())
571 set<pair<double, pair<const DMCThrown*, const DChargedTrackHypothesis*> > >::iterator locIterator = locParticleMatches.end();
572 set<const DMCThrown*> locMatchedThrowns;
573 set<const DChargedTrackHypothesis*> locMatchedHypotheses;
574 for(--locIterator; locIterator != locParticleMatches.begin(); --locIterator)
576 double locMatchFOM = locIterator->first;
577 const DMCThrown* locMCThrown = locIterator->second.first;
580 if(locMatchedThrowns.find(locMCThrown) != locMatchedThrowns.end())
582 if(locMatchedHypotheses.find(locChargedTrackHypothesis) != locMatchedHypotheses.end())
585 locMatchedThrowns.insert(locMCThrown);
588 deque<const DChargedTrackHypothesis*> locMatchedChargedHypos;
589 for(
size_t loc_i = 0; loc_i < locChargedTrackHypothesisVector.size(); ++loc_i)
593 cout <<
"CHECKING: CHARGEDHYPO: ";
594 cout <<
ParticleType(locChargedTrackHypothesisVector[loc_i]->PID()) <<
", " << locChargedTrackHypothesisVector[loc_i]->momentum().Mag() <<
", " << locChargedTrackHypothesisVector[loc_i]->momentum().Theta()*180.0/TMath::Pi() <<
", " << locChargedTrackHypothesisVector[loc_i]->momentum().Phi()*180.0/TMath::Pi() << endl;
595 cout <<
"best id, test id = " << locChargedTrackHypothesis->
Get_TrackTimeBased()->
candidateid <<
", " << locChargedTrackHypothesisVector[loc_i]->Get_TrackTimeBased()->candidateid << endl;
597 if(locChargedTrackHypothesisVector[loc_i]->Get_TrackTimeBased()->candidateid == locChargedTrackHypothesis->
Get_TrackTimeBased()->
candidateid)
600 cout <<
"save!" << endl;
601 locChargedToThrownMap[locChargedTrackHypothesisVector[loc_i]] = pair<const DMCThrown*, double>(locMCThrown, locMatchFOM);
602 locMatchedChargedHypos.push_back(locChargedTrackHypothesisVector[loc_i]);
603 locMatchedHypotheses.insert(locChargedTrackHypothesisVector[loc_i]);
606 locThrownToChargedMap[locMCThrown] = pair<deque<const DChargedTrackHypothesis*>,
double>(locMatchedChargedHypos, locMatchFOM);
616 map<const DNeutralParticleHypothesis*, pair<const DMCThrown*, double> > locNeutralHypoToThrownMap;
617 map<const DMCThrown*, pair<deque<const DNeutralParticleHypothesis*>,
double> > locThrownToNeutralHypoMap;
621 map<const DNeutralParticle*, pair<const DMCThrown*, double> > locNeutralToThrownMap;
622 map<const DMCThrown*, pair<const DNeutralParticle*, double> > locThrownToNeutralMap;
624 map<const DMCThrown*, pair<deque<const DNeutralParticleHypothesis*>,
double> >::iterator locIterator;
625 for(locIterator = locThrownToNeutralHypoMap.begin(); locIterator != locThrownToNeutralHypoMap.end(); ++locIterator)
627 for(
size_t loc_j = 0; loc_j < locNeutralParticles.size(); ++loc_j)
629 if((locIterator->second.first)[0]->Get_NeutralShower() == locNeutralParticles[loc_j]->dNeutralShower)
631 locNeutralToThrownMap[locNeutralParticles[loc_j]] = pair<const DMCThrown*, double>(locIterator->first, locIterator->second.second);
632 locThrownToNeutralMap[locIterator->first] = pair<const DNeutralParticle*, double>(locNeutralParticles[loc_j], locIterator->second.second);
643 map<const DNeutralParticleHypothesis*, pair<const DMCThrown*, double> > locNeutralToThrownMap;
644 map<const DMCThrown*, pair<deque<const DNeutralParticleHypothesis*>,
double> > locThrownToNeutralMap;
647 cout <<
"START IT!" << endl;
650 map<const DNeutralParticleHypothesis*, TMatrixDSym> locInverseCovMatrixMap;
651 for(
size_t loc_i = 0; loc_i < locNeutralParticleHypothesisVector.size(); ++loc_i)
654 TMatrixDSym locInverse3x3Matrix(3);
656 locInverseCovMatrixMap.insert(pair<const DNeutralParticleHypothesis*, TMatrixFSym>(locNeutralParticleHypothesis, locInverse3x3Matrix));
660 set<pair<double, pair<const DMCThrown*, const DNeutralParticleHypothesis*> > > locParticleMatches;
661 for(
size_t loc_i = 0; loc_i < locMCThrownVector.size(); ++loc_i)
663 const DMCThrown* locMCThrown = locMCThrownVector[loc_i];
664 for(
size_t loc_j = 0; loc_j < locNeutralParticleHypothesisVector.size(); ++loc_j)
667 if(locInverseCovMatrixMap.find(locNeutralParticleHypothesis) == locInverseCovMatrixMap.end())
669 TMatrixDSym& locInverse3x3Matrix = locInverseCovMatrixMap[locNeutralParticleHypothesis];
675 cout <<
"MATCHING: MCTHROWN: ";
677 cout <<
"MATCHING: NEUTRALHYPO: ";
678 cout <<
ParticleType(locNeutralParticleHypothesis->
PID()) <<
", " << locNeutralParticleHypothesis->
momentum().Mag() <<
", " << locNeutralParticleHypothesis->
momentum().Theta()*180.0/TMath::Pi() <<
", " << locNeutralParticleHypothesis->
momentum().Phi()*180.0/TMath::Pi() << endl;
679 cout <<
"MATCHING: FOM: " << locMatchFOM << endl;
682 pair<const DMCThrown*, const DNeutralParticleHypothesis*> locTrackPair(locMCThrown, locNeutralParticleHypothesis);
683 pair<double, pair<const DMCThrown*, const DNeutralParticleHypothesis*> > locMatchPair(locMatchFOM, locTrackPair);
684 locParticleMatches.insert(locMatchPair);
688 if(locParticleMatches.empty())
692 set<pair<double, pair<const DMCThrown*, const DNeutralParticleHypothesis*> > >::iterator locIterator = locParticleMatches.end();
693 set<const DMCThrown*> locMatchedThrowns;
694 set<const DNeutralParticleHypothesis*> locMatchedHypotheses;
695 for(--locIterator; locIterator != locParticleMatches.begin(); --locIterator)
697 double locMatchFOM = locIterator->first;
698 const DMCThrown* locMCThrown = locIterator->second.first;
701 if(locMatchedThrowns.find(locMCThrown) != locMatchedThrowns.end())
703 if(locMatchedHypotheses.find(locNeutralParticleHypothesis) != locMatchedHypotheses.end())
706 locMatchedThrowns.insert(locMCThrown);
709 deque<const DNeutralParticleHypothesis*> locMatchedNeutralHypos(1, locNeutralParticleHypothesis);
710 for(
size_t loc_i = 0; loc_i < locNeutralParticleHypothesisVector.size(); ++loc_i)
712 if(locNeutralParticleHypothesis->
Get_NeutralShower() == locNeutralParticleHypothesisVector[loc_i]->Get_NeutralShower())
714 locNeutralToThrownMap[locNeutralParticleHypothesisVector[loc_i]] = pair<const DMCThrown*, double>(locMCThrown, locMatchFOM);
715 locMatchedNeutralHypos.push_back(locNeutralParticleHypothesisVector[loc_i]);
716 locMatchedHypotheses.insert(locNeutralParticleHypothesisVector[loc_i]);
719 locThrownToNeutralMap[locMCThrown] = pair<deque<const DNeutralParticleHypothesis*>,
double>(locMatchedNeutralHypos, locMatchFOM);
728 double locTotalError =
sqrt(locInputCovarianceMatrix(0, 0) + locInputCovarianceMatrix(1, 1) + locInputCovarianceMatrix(2, 2));
732 locInverse3x3Matrix.ResizeTo(3, 3);
733 for(
unsigned int loc_i = 0; loc_i < 3; ++loc_i)
735 for(
unsigned int loc_j = 0; loc_j < 3; ++loc_j)
736 locInverse3x3Matrix(loc_i, loc_j) = locInputCovarianceMatrix(loc_i, loc_j);
740 TDecompLU locDecompLU(locInverse3x3Matrix);
742 if((!locDecompLU.Decompose()) || (fabs(locInverse3x3Matrix.Determinant()) < 1.0E-300))
744 locInverse3x3Matrix.Invert();
750 DVector3 locDeltaP3 = locMomentum_Detected - locMomentum_Thrown;
751 TMatrix locDeltas(3, 1);
752 locDeltas(0, 0) = locDeltaP3.Px();
753 locDeltas(1, 0) = locDeltaP3.Py();
754 locDeltas(2, 0) = locDeltaP3.Pz();
756 double locChiSq = (locInverse3x3Matrix.SimilarityT(locDeltas))(0, 0);
757 double locFOM = TMath::Prob(locChiSq, 3);
759 cout <<
"delta pxyz, chisq, FOM = " << locDeltaP3.Px() <<
", " << locDeltaP3.Py() <<
", " << locDeltaP3.Pz() <<
", " << locChiSq <<
", " << locFOM << endl;
void Set_ThrownToChargedMap(const map< const DMCThrown *, pair< const DChargedTrack *, double > > &locThrownToChargedMap)
void Set_ThrownToNeutralMap(const map< const DMCThrown *, pair< const DNeutralParticle *, double > > &locThrownToNeutralMap)
jerror_t brun(jana::JEventLoop *locEventLoop, int32_t runnumber)
Called everytime a new run number is detected.
double dMaximumBCALMatchAngleDegrees
jerror_t init(void)
Called once at program start.
void Find_GenReconMatches_NeutralHypo(const vector< const DMCThrown * > &locInputMCThrownVector, const vector< const DNeutralParticleHypothesis * > &locInputNeutralParticleHypothesisVector, DMCThrownMatching *locMCThrownMatching) const
bool Calc_InverseMatrix(const TMatrixFSym &locInputCovarianceMatrix, TMatrixDSym &locInverse3x3Matrix) const
double dMaximumFCALMatchDistance
void Set_FCALTruthToShowerMap(map< const DFCALTruthShower *, pair< const DFCALShower *, double > > &locFCALTruthToShowerMap)
void Set_FCALShowerToTruthMap(map< const DFCALShower *, pair< const DFCALTruthShower *, double > > &locFCALShowerToTruthMap)
double dMinTrackMatchHitFraction
void Set_NeutralHypoToThrownMap(const map< const DNeutralParticleHypothesis *, pair< const DMCThrown *, double > > &locNeutralHypoToThrownMap)
double energy(void) const
double dMaximumTOFMatchDistance
void Find_GenReconMatches_TOFPoints(JEventLoop *locEventLoop, DMCThrownMatching *locMCThrownMatching) const
void Set_BCALTruthToShowerMap(map< const DBCALTruthShower *, pair< const DBCALShower *, double > > &locBCALTruthToShowerMap)
const DTrackTimeBased * Get_TrackTimeBased(void) const
void Find_GenReconMatches_NeutralParticle(const vector< const DNeutralParticle * > &locNeutralParticles, DMCThrownMatching *locMCThrownMatching) const
static char * ParticleType(Particle_t p)
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
void Set_TOFTruthToPointMap(const map< const DTOFTruth *, pair< const DTOFPoint *, double > > &locTOFTruthToPointMap)
static int ParticleCharge(Particle_t p)
oid_t candidateid
id of DTrackCandidate corresponding to this track
const DNeutralShower * Get_NeutralShower(void) const
void Set_TaggedMCGENBeamPhoton(const DBeamPhoton *locBeamPhoton)
void Set_ThrownToChargedHypoMap(const map< const DMCThrown *, pair< deque< const DChargedTrackHypothesis * >, double > > &locThrownToChargedHypoMap)
void Set_ThrownToNeutralHypoMap(const map< const DMCThrown *, pair< deque< const DNeutralParticleHypothesis * >, double > > &locThrownToNeutralHypoMap)
size_t Get_NumTrackHits(const DTrackTimeBased *locTrackTimeBased)
jerror_t evnt(jana::JEventLoop *locEventLoop, uint64_t eventnumber)
Called every event.
const DChargedTrack * Get_MatchingChargedTrack(const DMCThrown *locInputMCThrown, double &locMatchFOM) const
void Find_GenReconMatches_BeamPhotons(JEventLoop *locEventLoop, DMCThrownMatching *locMCThrownMatching) const
void Get_ChargedHypoToThrownMap(map< const DChargedTrackHypothesis *, pair< const DMCThrown *, double > > &locChargedHypoToThrownMap) const
void Find_GenReconMatches_ChargedTrack(const vector< const DChargedTrack * > &locChargedTracks, DMCThrownMatching *locMCThrownMatching) const
void Set_BeamTruthToPhotonMap(map< const DBeamPhoton *, const DBeamPhoton * > &locBeamTruthToPhotonMap)
double Calc_MatchFOM(const DVector3 &locMomentum_Thrown, const DVector3 &locMomentum_Detected, TMatrixDSym locInverse3x3Matrix) const
void Set_TOFPointToTruthMap(const map< const DTOFPoint *, pair< const DTOFTruth *, double > > &locTOFPointToTruthMap)
DGeometry * GetDGeometry(unsigned int run_number)
jerror_t fini(void)
Called after last event of last event source has been processed.
void Set_ChargedToThrownMap(const map< const DChargedTrack *, pair< const DMCThrown *, double > > &locChargedToThrownMap)
void Find_GenReconMatches_ChargedHypo(const vector< const DMCThrown * > &locInputMCThrownVector, const vector< const DChargedTrackHypothesis * > &locInputChargedTrackHypothesisVector, DMCThrownMatching *locMCThrownMatching) const
void Set_NeutralToThrownMap(const map< const DNeutralParticle *, pair< const DMCThrown *, double > > &locNeutralToThrownMap)
void Find_GenReconMatches_FCALShowers(JEventLoop *locEventLoop, DMCThrownMatching *locMCThrownMatching) const
void Get_NeutralHypoToThrownMap(map< const DNeutralParticleHypothesis *, pair< const DMCThrown *, double > > &locNeutralHypoToThrownMap) const
const DVector3 & momentum(void) const
void Set_MCGENBeamPhoton(const DBeamPhoton *locBeamPhoton)
void Set_ChargedHypoToThrownMap(const map< const DChargedTrackHypothesis *, pair< const DMCThrown *, double > > &locChargedHypoToThrownMap)
void Set_BeamPhotonToTruthMap(map< const DBeamPhoton *, const DBeamPhoton * > &locBeamPhotonToTruthMap)
void Find_GenReconMatches_BCALShowers(JEventLoop *locEventLoop, DMCThrownMatching *locMCThrownMatching) const
void Get_ThrownToChargedHypoMap(map< const DMCThrown *, pair< deque< const DChargedTrackHypothesis * >, double > > &locThrownToChargedHypoMap) const
void Set_BCALShowerToTruthMap(map< const DBCALShower *, pair< const DBCALTruthShower *, double > > &locBCALShowerToTruthMap)
shared_ptr< const TMatrixFSym > errorMatrix(void) const
double dMaxTotalParticleErrorForMatch
int type
GEANT particle ID.
DVector3 getPosition() const
bool GetTargetZ(double &z_target) const
z-location of center of target
Particle_t PID(void) const
const DMCThrown * Get_MatchingMCThrown(const DChargedTrackHypothesis *locChargedTrackHypothesis, double &locMatchFOM) const
void Get_ThrownToNeutralHypoMap(map< const DMCThrown *, pair< deque< const DNeutralParticleHypothesis * >, double > > &locThrownToNeutralHypoMap) const