30 pthread_mutex_init(&mutex, NULL);
48 TDirectory *
dir = (TDirectory*)gROOT->FindObject(
"TRACKING");
49 if(!dir)dir =
new TDirectoryFile(
"TRACKING",
"TRACKING");
53 trkeff =
new TTree(
"trkeff2",
"Tracking Efficiency");
54 trkeff->Branch(
"E",
"track2",&trk_ptr);
58 JParameterManager *parms = app->GetJParameterManager();
62 parms->SetDefaultParameter(
"TRKEFF:DEBUG",
DEBUG);
80 double cdc_axial_length;
82 CDCZmin = CDCZmax-cdc_axial_length;
112 vector<const DParticleID *> locPIDAlgorithms;
113 loop->Get(locPIDAlgorithms);
114 if(locPIDAlgorithms.size() < 1){
115 _DBG_<<
"Unable to get a DParticleID object! NO PID will be done!"<<endl;
116 return RESOURCE_UNAVAILABLE;
119 locPIDAlgorithm =
const_cast<DParticleID*
>(locPIDAlgorithms[0]);
122 if(!locPIDAlgorithm){
123 _DBG_<<
"Unable to get a DParticleID object! NO PID will be done!"<<endl;
124 return RESOURCE_UNAVAILABLE;
128 vector<const DCDCTrackHit*> cdctrackhits;
129 loop->Get(cdctrackhits);
133 vector<const DFDCHit*> fdchits;
137 vector<const DMCThrown*> mcthrowns;
138 vector<const DMCTrajectoryPoint*> mctrajpoints;
139 loop->Get(mctrajpoints);
140 loop->Get(mcthrowns);
143 japp->RootWriteLock();
146 for(
unsigned int i=0; i<mcthrowns.size(); i++){
147 const DMCThrown *mcthrown = mcthrowns[i];
154 if(fabs(mcthrowns[i]->charge())==0.0)
continue;
158 trk.pthrown = pthrown;
159 trk.q_thrown = mcthrown->
charge();
160 trk.PID_thrown = mcthrown->
type;
163 trk.pfit.SetXYZ(0,0,0);
164 trk.pfit_wire.SetXYZ(0,0,0);
165 trk.pcan.SetXYZ(0,0,0);
166 trk.trk_chisq=1.0E20;
168 trk.trk_chisq_wb=1.0E20;
170 trk.delta_pt_over_pt=1.0E20;
171 trk.delta_theta=1.0E20;
172 trk.delta_phi=1.0E20;
173 trk.delta_pt_over_pt_wire=1.0E20;
174 trk.delta_theta_wire=1.0E20;
175 trk.delta_phi_wire=1.0E20;
176 trk.delta_pt_over_pt_can=1.0E20;
177 trk.delta_theta_can=1.0E20;
178 trk.delta_phi_can=1.0E20;
179 trk.isreconstructable = isReconstructable(mcthrown, mctrajpoints);
183 trk.NLR_bad_stereo = 0;
185 trk.event = eventnumber;
186 trk.dTrackReconstructedFlag_Candidate =
false;
187 trk.dTrackReconstructedFlag_WireBased =
false;
188 trk.dTrackReconstructedFlag_TimeBased =
false;
189 trk.q_candidate = 0.0;
190 trk.q_wirebased = 0.0;
191 trk.q_timebased = 0.0;
192 trk.PID_candidate = 0;
193 trk.PID_wirebased = 0;
194 trk.PID_timebased = 0;
195 trk.PID_hypothesized = 0;
196 trk.q_hypothesized = 0.0;
197 trk.FOM_hypothesized = 0.0;
214 locFoundFlag = Search_ChargedTrackHypotheses(loop, eventnumber, mcthrown);
215 if(locFoundFlag ==
false)
216 trk.dTrackReconstructedFlag_Candidate =
false;
218 locFoundFlag = Search_WireBasedTracks(loop, eventnumber, mcthrown);
219 if(locFoundFlag ==
false)
220 trk.dTrackReconstructedFlag_WireBased =
false;
222 locFoundFlag = Search_TrackCandidates(loop, eventnumber, mcthrown);
223 if(locFoundFlag ==
false)
224 trk.dTrackReconstructedFlag_Candidate =
false;
235 vector<const DChargedTrackHypothesis*> locChargedTrackHypotheses;
236 vector<const DMCTrajectoryPoint*> mctrajpoints;
237 vector<const DCDCTrackHit*> cdctrackhits;
238 vector<const DFDCHit*> fdchits;
240 loop->Get(cdctrackhits);
242 loop->Get(locChargedTrackHypotheses);
243 loop->Get(mctrajpoints);
247 bool locFoundFlag =
false;
248 double fom_best = 1.0E8;
250 trk.num_timebased = locChargedTrackHypotheses.size();
253 for(
unsigned int j=0; j<locChargedTrackHypotheses.size(); j++){
255 const DTrackTimeBased *locTimeBasedTrack = locChargedTrackHypothesis->dTrackTimeBased;
258 vector<const DTrackWireBased*> tracks;
259 locTimeBasedTrack->Get(tracks);
261 vector<const DTrackCandidate*> trackcandidates;
262 if(track)track->Get(trackcandidates);
268 double delta_pt_over_pt = (pfit.Perp() - pthrown.Perp())/pthrown.Perp();
269 double delta_theta = (pfit.Theta() - pthrown.Theta())*1000.0;
270 double delta_phi = (pfit.Phi() - pthrown.Phi())*1000.0;
275 double fom = pow(delta_pt_over_pt/0.02, 2.0) + pow(delta_theta/20.0, 2.0) + pow(delta_phi/20.0, 2.0);
279 trk.PID_hypothesized = int(locChargedTrackHypothesis->dPID);
280 trk.q_hypothesized = locChargedTrackHypothesis->
charge();
281 trk.FOM_hypothesized = locChargedTrackHypothesis->dFOM;
283 trk.dTrackReconstructedFlag_TimeBased =
true;
285 trk.q_timebased = locTimeBasedTrack->
charge();
287 trk.PID_timebased = int(locPIDAlgorithm->IDTrack(locTimeBasedTrack->
charge(), locTimeBasedTrack->
mass()));
290 trk.trk_chisq = locTimeBasedTrack->
chisq;
291 trk.trk_Ndof = locTimeBasedTrack->
Ndof;
292 trk.delta_pt_over_pt = delta_pt_over_pt;
293 trk.delta_theta = delta_theta;
294 trk.delta_phi = delta_phi;
297 vector<const DCDCTrackHit*> cdchits;
298 locTimeBasedTrack->Get(cdchits);
300 for(
unsigned int k=0; k<cdchits.size(); k++)
if(cdchits[k]->wire->stereo!=0.0)trk.Nstereo++;
301 trk.Ncdc = cdchits.size();
302 vector<const DFDCPseudo*> fdchits;
303 locTimeBasedTrack->Get(fdchits);
304 trk.Nfdc = fdchits.size();
309 vector<int> LRthrown;
311 vector<const DCoordinateSystem*> wires;
312 for(
unsigned int k=0; k<cdchits.size(); k++)wires.push_back(cdchits[k]->wire);
313 for(
unsigned int k=0; k<fdchits.size(); k++)wires.push_back(fdchits[k]->wire);
315 FindLR(wires, rt_thrown, LRthrown);
317 FindLR(wires, mctrajpoints, LRthrown);
319 FindLR(wires, locTimeBasedTrack->rt, LRfit);
322 if(LRfit.size()!=LRthrown.size() || LRfit.size()!=wires.size()){
323 _DBG_<<
"LR vector sizes do not match! LRfit.size()="<<LRfit.size()<<
" LRthrown.size()="<<LRthrown.size()<<
" wires.size()="<<wires.size()<<endl;
328 trk.NLR_bad_stereo = trk.NLR_bad = 0;
329 for(
unsigned int k=0; k<wires.size(); k++){
330 if(LRfit[k] == LRthrown[k])
continue;
332 bool is_stereo = (wires[k]->udir.Theta()*57.3)>2.0;
333 if(is_stereo)trk.NLR_bad_stereo++;
341 vector<const DMCTrajectoryPoint*> mctrajpoints;
342 vector<const DCDCTrackHit*> cdctrackhits;
343 vector<const DFDCHit*> fdchits;
344 vector<const DTrackWireBased*> tracks;
346 loop->Get(cdctrackhits);
349 loop->Get(mctrajpoints);
352 double fom_best = 1.0E8;
354 trk.num_wirebased = tracks.size();
356 bool locFoundFlag =
false;
358 for(
unsigned int j=0; j<tracks.size(); j++){
365 double delta_pt_over_pt_wire = (pfit_wire.Perp() - pthrown.Perp())/pthrown.Perp();
366 double delta_theta_wire = (pfit_wire.Theta() - pthrown.Theta())*1000.0;
367 double delta_phi_wire = (pfit_wire.Phi() - pthrown.Phi())*1000.0;
372 double fom = pow(delta_pt_over_pt_wire/0.02, 2.0) + pow(delta_theta_wire/20.0, 2.0) + pow(delta_phi_wire/20.0, 2.0);
376 trk.dTrackReconstructedFlag_WireBased =
true;
378 trk.q_wirebased = track->
charge();
380 trk.PID_wirebased = int(locPIDAlgorithm->IDTrack(track->
charge(), track->
mass()));
382 trk.pfit_wire = pfit_wire;
383 trk.trk_chisq_wb = track->
chisq;
384 trk.trk_Ndof_wb = track->
Ndof;
385 trk.delta_pt_over_pt_wire = delta_pt_over_pt_wire;
386 trk.delta_theta_wire = delta_theta_wire;
387 trk.delta_phi_wire = delta_phi_wire;
390 vector<const DCDCTrackHit*> cdchits;
393 for(
unsigned int k=0; k<cdchits.size(); k++)
if(cdchits[k]->wire->stereo!=0.0)trk.Nstereo++;
394 trk.Ncdc = cdchits.size();
395 vector<const DFDCPseudo*> fdchits;
397 trk.Nfdc = fdchits.size();
402 vector<int> LRthrown;
404 vector<const DCoordinateSystem*> wires;
405 for(
unsigned int k=0; k<cdchits.size(); k++)wires.push_back(cdchits[k]->wire);
406 for(
unsigned int k=0; k<fdchits.size(); k++)wires.push_back(fdchits[k]->wire);
408 FindLR(wires, rt_thrown, LRthrown);
410 FindLR(wires, mctrajpoints, LRthrown);
412 FindLR(wires, track->rt, LRfit);
415 if(LRfit.size()!=LRthrown.size() || LRfit.size()!=wires.size()){
416 _DBG_<<
"LR vector sizes do not match! LRfit.size()="<<LRfit.size()<<
" LRthrown.size()="<<LRthrown.size()<<
" wires.size()="<<wires.size()<<endl;
421 trk.NLR_bad_stereo = trk.NLR_bad = 0;
422 for(
unsigned int k=0; k<wires.size(); k++){
423 if(LRfit[k] == LRthrown[k])
continue;
425 bool is_stereo = (wires[k]->udir.Theta()*57.3)>2.0;
426 if(is_stereo)trk.NLR_bad_stereo++;
434 vector<const DMCTrajectoryPoint*> mctrajpoints;
435 vector<const DCDCTrackHit*> cdctrackhits;
436 vector<const DFDCHit*> fdchits;
437 vector<const DTrackCandidate*> trackcandidates;
439 loop->Get(cdctrackhits);
441 loop->Get(mctrajpoints);
442 loop->Get(trackcandidates);
445 double fom_best = 1.0E8;
447 trk.num_candidates = trackcandidates.size();
449 bool locFoundFlag =
false;
451 for(
unsigned int j=0; j<trackcandidates.size(); j++){
458 double delta_pt_over_pt_can = (pcan.Perp() - pthrown.Perp())/pthrown.Perp();
459 double delta_theta_can = (pcan.Theta() - pthrown.Theta())*1000.0;
460 double delta_phi_can = (pcan.Phi() - pthrown.Phi())*1000.0;
465 double fom = pow(delta_pt_over_pt_can/0.02, 2.0) + pow(delta_theta_can/20.0, 2.0) + pow(delta_phi_can/20.0, 2.0);
469 trk.dTrackReconstructedFlag_Candidate =
true;
470 trk.q_candidate = trackcandidate->
charge();
471 trk.PID_candidate = int(locPIDAlgorithm->IDTrack(trackcandidate->
charge(), trackcandidate->
mass()));
473 trk.delta_pt_over_pt_can = delta_pt_over_pt_can;
474 trk.delta_theta_can = delta_theta_can;
475 trk.delta_phi_can = delta_phi_can;
478 vector<const DCDCTrackHit*> cdchits;
479 trackcandidate->Get(cdchits);
481 for(
unsigned int k=0; k<cdchits.size(); k++)
if(cdchits[k]->wire->stereo!=0.0)trk.Nstereo++;
482 trk.Ncdc = cdchits.size();
483 vector<const DFDCPseudo*> fdchits;
484 trackcandidate->Get(fdchits);
485 trk.Nfdc = fdchits.size();
516 for(
unsigned int i=0; i<mctrajpoints.size(); i++){
518 if(mctraj->
track != mcthrown->
myid)
continue;
520 if(mctraj_first==NULL)mctraj_first = mctraj;
521 mctraj_last = mctraj;
524 if(mctraj_last!=NULL){
525 double r =
sqrt(pow((
double)mctraj_last->
x,2.0) + pow((
double)mctraj_last->
y,2.0));
526 if(r>45.0)
return true;
527 if(mctraj_last->
z<17.0)
return true;
528 if(mctraj_last->
z<17.0)
return true;
547 for(
unsigned int i=0; i<wires.size(); i++){
556 DVector3 pos_diff = pos_doca-pos_wire;
558 LRhits.push_back(shift.Dot(pos_diff)<0.0 ? -1:1);
bool isReconstructable(const DMCThrown *mcthrown, vector< const DMCTrajectoryPoint * > &mctrajpoints)
float chisq
Chi-squared for the track (not chisq/dof!)
float chisq
Chi-squared for the track (not chisq/dof!)
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
DVector3 GetLastDOCAPoint(void) const
const DVector3 & position(void) const
double GetLastDistAlongWire(void) const
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
bool GetCDCAxialLength(double &cdc_axial_length) const
length of CDC axial wires in cm
DMagneticFieldMap * GetBfield(void) const
double DistToRT(double x, double y, double z) const
bool Search_TrackCandidates(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown)
~DEventProcessor_trackeff_hists2()
DGeometry * GetDGeometry(unsigned int run_number)
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
<A href="index.html#legend"> <IMG src="CORE.png" width="100"> </A>
int Ndof
Number of degrees of freedom in the fit.
jerror_t init(void)
Invoked via DEventProcessor virtual method.
const DVector3 & momentum(void) const
void FindLR(vector< const DCoordinateSystem * > &wires, const DReferenceTrajectory *crt, vector< int > &LRhits)
DEventProcessor_trackeff_hists2()
bool Search_WireBasedTracks(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown)
jerror_t brun(JEventLoop *loop, int32_t runnumber)
bool GetCDCEndplate(double &z, double &dz, double &rmin, double &rmax) const
int type
GEANT particle ID.
bool Search_ChargedTrackHypotheses(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown)
int myid
id of this particle from original generator