19 dMinTrackingFOM = 0.0;
22 gPARMS->SetDefaultParameter(
"EVENTRFBUNCH:USE_TAG", OVERRIDE_TAG,
"Use a particular tag for the general RF bunch calculation instead of the default calculation.");
34 vector<double> locBeamPeriodVector;
35 locEventLoop->GetCalib(
"PHOTON_BEAM/RF/beam_period", locBeamPeriodVector);
36 dBeamBunchPeriod = locBeamPeriodVector[0];
38 double locTargetCenterZ;
40 dTargetCenter.SetXYZ(0.0, 0.0, locTargetCenterZ);
42 locEventLoop->GetSingle(dParticleID);
55 if(OVERRIDE_TAG !=
"") {
56 vector<const DEventRFBunch *> rfBunchVec;
57 locEventLoop->Get(rfBunchVec, OVERRIDE_TAG.c_str());
59 if(rfBunchVec.size() > 0) {
61 _data.push_back(rfBunchCopy);
64 jerr <<
"Could not find DEventRFBunch objects with tag " << OVERRIDE_TAG
65 <<
", falling back to default calculation ..." << endl;
71 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
72 Select_GoodTracks(locEventLoop, locTrackTimeBasedVector);
75 vector<const DRFTime*> locRFTimes;
76 locEventLoop->Get(locRFTimes);
77 if(!locRFTimes.empty())
78 return Select_RFBunch(locEventLoop, locTrackTimeBasedVector, locRFTimes[0]);
80 return Select_RFBunch_NoRFTime(locEventLoop, locTrackTimeBasedVector);
89 locSelectedTimeBasedTracks.clear();
91 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
92 locEventLoop->Get(locTrackTimeBasedVector);
95 map<JObject::oid_t, const DTrackTimeBased*> locBestTrackTimeBasedMap;
96 for(
size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
98 JObject::oid_t locCandidateID = locTrackTimeBasedVector[loc_i]->candidateid;
99 if(locBestTrackTimeBasedMap.find(locCandidateID) == locBestTrackTimeBasedMap.end())
100 locBestTrackTimeBasedMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
101 else if(locTrackTimeBasedVector[loc_i]->FOM > locBestTrackTimeBasedMap[locCandidateID]->FOM)
102 locBestTrackTimeBasedMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
106 map<JObject::oid_t, const DTrackTimeBased*>::iterator locIterator;
107 for(locIterator = locBestTrackTimeBasedMap.begin(); locIterator != locBestTrackTimeBasedMap.end(); ++locIterator)
110 if(locTrackTimeBased->
FOM >= dMinTrackingFOM)
111 locSelectedTimeBasedTracks.push_back(locTrackTimeBased);
134 locEventLoop->GetSingle(locDetectorMatches);
136 vector<pair<double, const JObject*> > locTimes;
137 int locBestRFBunchShift = 0, locHighestNumVotes = 0;
139 if(Find_TrackTimes_All(locDetectorMatches, locTrackTimeBasedVector, locTimes))
140 locBestRFBunchShift = Conduct_Vote(locEventLoop, locRFTime->
dTime, locTimes,
true, locHighestNumVotes);
141 else if(Find_NeutralTimes(locEventLoop, locTimes))
142 locBestRFBunchShift = Conduct_Vote(locEventLoop, locRFTime->
dTime, locTimes,
false, locHighestNumVotes);
144 locBestRFBunchShift = 0;
147 locEventRFBunch->
dTime = locRFTime->
dTime + dBeamBunchPeriod*double(locBestRFBunchShift);
151 _data.push_back(locEventRFBunch);
158 map<int, vector<const JObject*> > locNumBeamBucketsShiftedMap;
159 set<int> locBestRFBunchShifts;
161 locHighestNumVotes = Find_BestRFBunchShifts(locRFTime, locTimes, locNumBeamBucketsShiftedMap, locBestRFBunchShifts);
162 if(locBestRFBunchShifts.size() == 1)
163 return *locBestRFBunchShifts.begin();
166 vector<const DBeamPhoton*> locBeamPhotons;
167 locEventLoop->Get(locBeamPhotons);
168 if(Break_TieVote_BeamPhotons(locBeamPhotons, locRFTime, locNumBeamBucketsShiftedMap, locBestRFBunchShifts, locHighestNumVotes))
169 return *locBestRFBunchShifts.begin();
172 if(locUsedTracksFlag)
175 return Break_TieVote_Tracks(locNumBeamBucketsShiftedMap, locBestRFBunchShifts);
180 return Break_TieVote_Neutrals(locNumBeamBucketsShiftedMap, locBestRFBunchShifts);
187 for(
size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
189 const DTrackTimeBased* locTrackTimeBased = locTrackTimeBasedVector[loc_i];
190 double locP = locTrackTimeBased->
momentum().Mag();
196 shared_ptr<const DTOFHitMatchParams> locTOFHitMatchParams;
197 if((locP > 0.5) && dParticleID->Get_BestTOFMatchParams(locTrackTimeBased, locDetectorMatches, locTOFHitMatchParams))
199 double locPropagatedTime = locTOFHitMatchParams->dHitTime - locTOFHitMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->
z())/29.9792458;
200 locTimes.push_back(pair<double, const JObject*>(locPropagatedTime, locTrackTimeBased));
205 shared_ptr<const DSCHitMatchParams> locSCHitMatchParams;
206 if(!dParticleID->Get_BestSCMatchParams(locTrackTimeBased, locDetectorMatches, locSCHitMatchParams))
209 double locPropagatedTime = locSCHitMatchParams->dHitTime - locSCHitMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->
z())/29.9792458;
210 locTimes.push_back(pair<double, const JObject*>(locPropagatedTime, locTrackTimeBased));
213 return (!locTimes.empty());
223 for(
size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
225 const DTrackTimeBased* locTrackTimeBased = locTrackTimeBasedVector[loc_i];
228 shared_ptr<const DSCHitMatchParams> locSCHitMatchParams;
229 if(dParticleID->Get_BestSCMatchParams(locTrackTimeBased, locDetectorMatches, locSCHitMatchParams))
231 double locPropagatedTime = locSCHitMatchParams->dHitTime - locSCHitMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->
z())/29.9792458;
232 locTimes.push_back(pair<double, const JObject*>(locPropagatedTime, locTrackTimeBased));
242 shared_ptr<const DTOFHitMatchParams> locTOFHitMatchParams;
243 if(dParticleID->Get_BestTOFMatchParams(locTrackTimeBased, locDetectorMatches, locTOFHitMatchParams))
245 double locPropagatedTime = locTOFHitMatchParams->dHitTime - locTOFHitMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->
z())/29.9792458;
246 locTimes.push_back(pair<double, const JObject*>(locPropagatedTime, locTrackTimeBased));
251 double locP = locTrackTimeBased->
momentum().Mag();
254 if((locP < 0.25) && (dBeamBunchPeriod < 2.005))
257 shared_ptr<const DBCALShowerMatchParams> locBCALShowerMatchParams;
258 if(dParticleID->Get_BestBCALMatchParams(locTrackTimeBased, locDetectorMatches, locBCALShowerMatchParams))
260 const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
261 double locPropagatedTime = locBCALShower->
t - locBCALShowerMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->
z())/29.9792458;
262 locTimes.push_back(pair<double, const JObject*>(locPropagatedTime, locTrackTimeBased));
267 shared_ptr<const DFCALShowerMatchParams> locFCALShowerMatchParams;
268 if(dParticleID->Get_BestFCALMatchParams(locTrackTimeBased, locDetectorMatches, locFCALShowerMatchParams))
270 const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
271 double locPropagatedTime = locFCALShower->
getTime() - locFCALShowerMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->
z())/29.9792458;
272 locTimes.push_back(pair<double, const JObject*>(locPropagatedTime, locTrackTimeBased));
279 return (!locTimes.empty());
286 vector< const DFCALShower* > fcalShowers;
287 locEventLoop->Get( fcalShowers );
289 for(
size_t i = 0; i < fcalShowers.size(); ++i ){
291 DVector3 locHitPoint = fcalShowers[i]->getPosition();
292 DVector3 locPath = locHitPoint - dTargetCenter;
293 double locPathLength = locPath.Mag();
294 if(!(locPathLength > 0.0))
297 double locFlightTime = locPathLength/29.9792458;
298 double locHitTime = fcalShowers[i]->getTime();
299 locTimes.push_back( pair< double, const JObject*>( locHitTime - locFlightTime, fcalShowers[i] ) );
302 vector< const DBCALShower* > bcalShowers;
303 locEventLoop->Get( bcalShowers );
305 for(
size_t i = 0; i < bcalShowers.size(); ++i ){
307 DVector3 locHitPoint( bcalShowers[i]->
x, bcalShowers[i]->
y, bcalShowers[i]->z );
308 DVector3 locPath = locHitPoint - dTargetCenter;
309 double locPathLength = locPath.Mag();
310 if(!(locPathLength > 0.0))
313 double locFlightTime = locPathLength/29.9792458;
314 double locHitTime = bcalShowers[i]->t;
315 locTimes.push_back( pair< double, const JObject*>( locHitTime - locFlightTime, bcalShowers[i] ) );
318 return (locTimes.size() > 1);
324 int locHighestNumVotes = 0;
325 locNumBeamBucketsShiftedMap.clear();
326 locBestRFBunchShifts.clear();
328 for(
unsigned int loc_i = 0; loc_i < locTimes.size(); ++loc_i)
330 double locDeltaT = locTimes[loc_i].first - locRFHitTime;
331 int locNumBeamBucketsShifted = (locDeltaT > 0.0) ? int(locDeltaT/dBeamBunchPeriod + 0.5) : int(locDeltaT/dBeamBunchPeriod - 0.5);
332 locNumBeamBucketsShiftedMap[locNumBeamBucketsShifted].push_back(locTimes[loc_i].second);
334 int locNumVotesThisShift = locNumBeamBucketsShiftedMap[locNumBeamBucketsShifted].size();
335 if(locNumVotesThisShift > locHighestNumVotes)
337 locHighestNumVotes = locNumVotesThisShift;
338 locBestRFBunchShifts.clear();
339 locBestRFBunchShifts.insert(locNumBeamBucketsShifted);
341 else if(locNumVotesThisShift == locHighestNumVotes)
342 locBestRFBunchShifts.insert(locNumBeamBucketsShifted);
345 return locHighestNumVotes;
353 set<int> locInputRFBunchShifts = locBestRFBunchShifts;
354 for(
size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
356 double locDeltaT = locBeamPhotons[loc_i]->time() - locRFTime;
357 int locNumBeamBucketsShifted = (locDeltaT > 0.0) ? int(locDeltaT/dBeamBunchPeriod + 0.5) : int(locDeltaT/dBeamBunchPeriod - 0.5);
358 if(locInputRFBunchShifts.find(locNumBeamBucketsShifted) == locInputRFBunchShifts.end())
361 locNumBeamBucketsShiftedMap[locNumBeamBucketsShifted].push_back(locBeamPhotons[loc_i]);
363 int locNumVotesThisShift = locNumBeamBucketsShiftedMap[locNumBeamBucketsShifted].size();
364 if(locNumVotesThisShift > locHighestNumVotes)
366 locHighestNumVotes = locNumVotesThisShift;
367 locBestRFBunchShifts.clear();
368 locBestRFBunchShifts.insert(locNumBeamBucketsShifted);
370 else if(locNumVotesThisShift == locHighestNumVotes)
371 locBestRFBunchShifts.insert(locNumBeamBucketsShifted);
374 return (locBestRFBunchShifts.size() == 1);
383 int locBestRFBunchShift = 0;
384 pair<int, double> locBestTrackingTotals(-1, 9.9E99);
386 set<int>::const_iterator locSetIterator = locBestRFBunchShifts.begin();
387 for(; locSetIterator != locBestRFBunchShifts.end(); ++locSetIterator)
389 int locRFBunchShift = *locSetIterator;
390 int locTotalTrackingNDF = 0;
391 double locTotalTrackingChiSq = 0.0;
393 const vector<const JObject*>& locVoters = locNumBeamBucketsShiftedMap[locRFBunchShift];
394 for(
size_t loc_i = 0; loc_i < locVoters.size(); ++loc_i)
397 if(locTrackTimeBased == NULL)
399 locTotalTrackingNDF += locTrackTimeBased->
Ndof;
400 locTotalTrackingChiSq += locTrackTimeBased->
chisq;
403 if(locTotalTrackingNDF > locBestTrackingTotals.first)
405 locBestTrackingTotals = pair<int, double>(locTotalTrackingNDF, locTotalTrackingChiSq);
406 locBestRFBunchShift = locRFBunchShift;
408 else if((locTotalTrackingNDF == locBestTrackingTotals.first) && (locTotalTrackingChiSq < locBestTrackingTotals.second))
410 locBestTrackingTotals = pair<int, double>(locTotalTrackingNDF, locTotalTrackingChiSq);
411 locBestRFBunchShift = locRFBunchShift;
415 return locBestRFBunchShift;
422 int locBestRFBunchShift = 0;
423 double locHighestTotalEnergy = 0.0;
425 set<int>::const_iterator locSetIterator = locBestRFBunchShifts.begin();
426 for(; locSetIterator != locBestRFBunchShifts.end(); ++locSetIterator)
428 int locRFBunchShift = *locSetIterator;
429 double locTotalEnergy = 0.0;
431 const vector<const JObject*>& locVoters = locNumBeamBucketsShiftedMap[locRFBunchShift];
432 for(
size_t loc_i = 0; loc_i < locVoters.size(); ++loc_i)
438 if( fcalShower != NULL ){
440 locTotalEnergy += fcalShower->
getEnergy();
445 if( bcalShower != NULL ){
447 locTotalEnergy += bcalShower->
E;
452 if(locTotalEnergy > locHighestTotalEnergy)
454 locHighestTotalEnergy = locTotalEnergy;
455 locBestRFBunchShift = locRFBunchShift;
459 return locBestRFBunchShift;
472 locEventLoop->GetSingle(locDetectorMatches);
474 vector<pair<double, const JObject*> > locTimes;
475 if(!Find_TrackTimes_SCTOF(locDetectorMatches, locTrackTimeBasedVector, locTimes))
476 return Create_NaNRFBunch();
479 double locRFTimeGuess, locTimeVariance;
480 Get_RFTimeGuess(locTimes, locRFTimeGuess, locTimeVariance);
483 int locHighestNumVotes = 0;
484 int locBestRFBunchShift = Conduct_Vote(locEventLoop, locRFTimeGuess, locTimes,
true, locHighestNumVotes);
487 locEventRFBunch->
dTime = locRFTimeGuess + dBeamBunchPeriod*double(locBestRFBunchShift);
491 _data.push_back(locEventRFBunch);
504 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
505 Select_GoodTracks(locEventLoop, locTrackTimeBasedVector);
508 locEventLoop->GetSingle(locDetectorMatches);
511 vector<pair<double, const JObject*> > locTimes;
512 if(!Find_TrackTimes_SCTOF(locDetectorMatches, locTrackTimeBasedVector, locTimes))
516 Get_RFTimeGuess(locTimes, locRFTimeGuess, locRFVariance);
526 locRFTimeGuess = locTimes[0].first;
527 for(
size_t loc_i = 1; loc_i < locTimes.size(); ++loc_i)
529 double locDeltaT = locTimes[loc_i].first - locTimes[0].first;
530 double locNumBeamBucketsShifted = (locDeltaT > 0.0) ? floor(locDeltaT/dBeamBunchPeriod + 0.5) : floor(locDeltaT/dBeamBunchPeriod - 0.5);
531 locRFTimeGuess += locTimes[loc_i].first - locNumBeamBucketsShifted*dBeamBunchPeriod;
533 locRFTimeGuess /= double(locTimes.size());
535 locRFVariance = 0.3*0.3/double(locTimes.size());
541 locEventRFBunch->
dTime = numeric_limits<double>::quiet_NaN();
542 locEventRFBunch->
dTimeVariance = numeric_limits<double>::quiet_NaN();
545 _data.push_back(locEventRFBunch);
jerror_t Select_RFBunch_NoRFTime(JEventLoop *locEventLoop, vector< const DTrackTimeBased * > &locTrackTimeBasedVector)
jerror_t brun(jana::JEventLoop *locEventLoop, int32_t runnumber)
Called everytime a new run number is detected.
jerror_t init(void)
Called once at program start.
float chisq
Chi-squared for the track (not chisq/dof!)
bool Get_RFTimeGuess(JEventLoop *locEventLoop, double &locRFTimeGuess, double &locRFVariance, DetectorSystem_t &locTimeSource) const
bool Break_TieVote_BeamPhotons(vector< const DBeamPhoton * > &locBeamPhotons, double locRFTime, map< int, vector< const JObject * > > &locNumBeamBucketsShiftedMap, set< int > &locBestRFBunchShifts, int locHighestNumVotes)
int Find_BestRFBunchShifts(double locRFHitTime, const vector< pair< double, const JObject * > > &locTimes, map< int, vector< const JObject * > > &locNumBeamBucketsShiftedMap, set< int > &locBestRFBunchShifts)
jerror_t Create_NaNRFBunch(void)
bool Find_NeutralTimes(JEventLoop *locEventLoop, vector< pair< double, const JObject * > > &locTimes)
jerror_t evnt(jana::JEventLoop *locEventLoop, uint64_t eventnumber)
Called every event.
jerror_t fini(void)
Called after last event of last event source has been processed.
DGeometry * GetDGeometry(unsigned int run_number)
int Conduct_Vote(JEventLoop *locEventLoop, double locRFTime, vector< pair< double, const JObject * > > &locTimes, bool locUsedTracksFlag, int &locHighestNumVotes)
int Ndof
Number of degrees of freedom in the fit.
int Break_TieVote_Neutrals(map< int, vector< const JObject * > > &locNumBeamBucketsShiftedMap, set< int > &locBestRFBunchShifts)
bool Find_TrackTimes_SCTOF(const DDetectorMatches *locDetectorMatches, const vector< const DTrackTimeBased * > &locTrackTimeBasedVector, vector< pair< double, const JObject * > > &locTimes) const
int Break_TieVote_Tracks(map< int, vector< const JObject * > > &locNumBeamBucketsShiftedMap, set< int > &locBestRFBunchShifts)
const DVector3 & momentum(void) const
void Select_GoodTracks(JEventLoop *locEventLoop, vector< const DTrackTimeBased * > &locSelectedTimeBasedTracks) const
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
DetectorSystem_t dTimeSource
unsigned int dNumParticleVotes
bool GetTargetZ(double &z_target) const
z-location of center of target
jerror_t Select_RFBunch(JEventLoop *locEventLoop, vector< const DTrackTimeBased * > &locTrackTimeBasedVector, const DRFTime *locRFTime)
bool Find_TrackTimes_All(const DDetectorMatches *locDetectorMatches, const vector< const DTrackTimeBased * > &locTrackTimeBasedVector, vector< pair< double, const JObject * > > &locTimes)