Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventProcessor_trackeff_hists2.cc
Go to the documentation of this file.
1 // $Id: DEventProcessor_trackeff_hists2.cc 2774 2007-07-19 15:59:02Z davidl $
2 //
3 // File: DEventProcessor_trackeff_hists2.cc
4 // Created: Sun Apr 24 06:45:21 EDT 2005
5 // Creator: davidl (on Darwin Harriet.local 7.8.0 powerpc)
6 //
7 
9 
10 using namespace std;
11 
12 // Routine used to create our DEventProcessor
13 extern "C"{
14 void InitPlugin(JApplication *app){
15  InitJANAPlugin(app);
16  app->AddProcessor(new DEventProcessor_trackeff_hists2());
17 }
18 } // "C"
19 
20 
21 //------------------
22 // DEventProcessor_trackeff_hists2
23 //------------------
25 {
26  trk_ptr = &trk;
27 
28  //trkres = new DTrackingResolutionGEANT();
29 
30  pthread_mutex_init(&mutex, NULL);
31 
32 }
33 
34 //------------------
35 // ~DEventProcessor_trackeff_hists2
36 //------------------
38 {
39 
40 }
41 
42 //------------------
43 // init
44 //------------------
46 {
47  // Create TRACKING directory
48  TDirectory *dir = (TDirectory*)gROOT->FindObject("TRACKING");
49  if(!dir)dir = new TDirectoryFile("TRACKING","TRACKING");
50  dir->cd();
51 
52  // Create Trees
53  trkeff = new TTree("trkeff2","Tracking Efficiency");
54  trkeff->Branch("E","track2",&trk_ptr);
55 
56  dir->cd("../");
57 
58  JParameterManager *parms = app->GetJParameterManager();
59 
60  DEBUG = 1;
61 
62  parms->SetDefaultParameter("TRKEFF:DEBUG", DEBUG);
63 
64  return NOERROR;
65 }
66 
67 //------------------
68 // brun
69 //------------------
70 jerror_t DEventProcessor_trackeff_hists2::brun(JEventLoop *loop, int32_t runnumber)
71 {
72  DApplication *dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
73  const DGeometry *dgeom = dapp->GetDGeometry(runnumber);
74 
75  rt_thrown = new DReferenceTrajectory(dgeom->GetBfield(runnumber));
76 
77  double dz, rmin, rmax;
78  dgeom->GetCDCEndplate(CDCZmax, dz, rmin, rmax);
79 
80  double cdc_axial_length;
81  dgeom->GetCDCAxialLength(cdc_axial_length);
82  CDCZmin = CDCZmax-cdc_axial_length;
83 
84  use_rt_thrown = true; //mctrajpoints.size()<20;
85 
86  return NOERROR;
87 }
88 
89 //------------------
90 // erun
91 //------------------
93 {
94  return NOERROR;
95 }
96 
97 //------------------
98 // fini
99 //------------------
101 {
102  return NOERROR;
103 }
104 
105 //------------------
106 // evnt
107 //------------------
108 jerror_t DEventProcessor_trackeff_hists2::evnt(JEventLoop *loop, uint64_t eventnumber)
109 {
110 
111  // Get the particle ID algorithms
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;
117  }
118  // Drop the const qualifier from the DParticleID pointer (I'm surely going to hell for this!)
119  locPIDAlgorithm = const_cast<DParticleID*>(locPIDAlgorithms[0]);
120 
121  // Warn user if something happened that caused us NOT to get a locPIDAlgorithm object pointer
122  if(!locPIDAlgorithm){
123  _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
124  return RESOURCE_UNAVAILABLE;
125  }
126 
127  // Bail quick on events with too many or too few CDC hits
128  vector<const DCDCTrackHit*> cdctrackhits;
129  loop->Get(cdctrackhits);
130  //if(cdctrackhits.size()>30 || cdctrackhits.size()<6)return NOERROR;
131 
132  // Bail quick on events with too many FDC hits
133  vector<const DFDCHit*> fdchits;
134  loop->Get(fdchits);
135  //if(fdchits.size()>30)return NOERROR;
136 
137  vector<const DMCThrown*> mcthrowns;
138  vector<const DMCTrajectoryPoint*> mctrajpoints;
139  loop->Get(mctrajpoints);
140  loop->Get(mcthrowns);
141 
142  // Although we are only filling objects local to this plugin, TTree::Fill() periodically writes to file: Global ROOT lock
143  japp->RootWriteLock(); //ACQUIRE ROOT LOCK
144 
145  // Get hit list for all throwns
146  for(unsigned int i=0; i<mcthrowns.size(); i++){
147  const DMCThrown *mcthrown = mcthrowns[i];
148 
149  // If there aren't enough DMCTrajectoryPoint objects then we will need to
150  // get the LR information by swimming the thrown value ourself.
151  if(use_rt_thrown)rt_thrown->Swim(mcthrown->position(), mcthrown->momentum(), mcthrown->charge());
152 
153  // if this isn't a charged track, then skip it
154  if(fabs(mcthrowns[i]->charge())==0.0)continue;
155 
156  // Momentum of thrown particle
157  DVector3 pthrown = mcthrown->momentum();
158  trk.pthrown = pthrown;
159  trk.q_thrown = mcthrown->charge(); //make sure this value isn't bogus!!
160  trk.PID_thrown = mcthrown->type;
161 
162  // Initialize with the "not found" values
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;
167  trk.trk_Ndof=1;
168  trk.trk_chisq_wb=1.0E20;
169  trk.trk_Ndof_wb=1;
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);
180  trk.Nstereo = 0;
181  trk.Ncdc = 0;
182  trk.Nfdc = 0;
183  trk.NLR_bad_stereo = 0;
184  trk.NLR_bad = 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;
198 
199  bool locFoundFlag;
200 /*
201  locFoundFlag = Search_ChargedTrackHypotheses(loop, eventnumber, mcthrown);
202  if(locFoundFlag == false){
203  trk.dTrackReconstructedFlag_TimeBased = false;
204  locFoundFlag = Search_WireBasedTracks(loop, eventnumber, mcthrown);
205  if(locFoundFlag == false){
206  trk.dTrackReconstructedFlag_WireBased = false;
207  locFoundFlag = Search_TrackCandidates(loop, eventnumber, mcthrown);
208  if(locFoundFlag == false)
209  trk.dTrackReconstructedFlag_Candidate = false;
210  }
211  }
212 */
213 
214  locFoundFlag = Search_ChargedTrackHypotheses(loop, eventnumber, mcthrown);
215  if(locFoundFlag == false)
216  trk.dTrackReconstructedFlag_Candidate = false;
217 
218  locFoundFlag = Search_WireBasedTracks(loop, eventnumber, mcthrown);
219  if(locFoundFlag == false)
220  trk.dTrackReconstructedFlag_WireBased = false;
221 
222  locFoundFlag = Search_TrackCandidates(loop, eventnumber, mcthrown);
223  if(locFoundFlag == false)
224  trk.dTrackReconstructedFlag_Candidate = false;
225 
226  trkeff->Fill();
227  }
228 
229  japp->RootUnLock(); //RELEASE ROOT LOCK
230 
231  return NOERROR;
232 }
233 
234 bool DEventProcessor_trackeff_hists2::Search_ChargedTrackHypotheses(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown){
235  vector<const DChargedTrackHypothesis*> locChargedTrackHypotheses;
236  vector<const DMCTrajectoryPoint*> mctrajpoints;
237  vector<const DCDCTrackHit*> cdctrackhits;
238  vector<const DFDCHit*> fdchits;
239 
240  loop->Get(cdctrackhits);
241  loop->Get(fdchits);
242  loop->Get(locChargedTrackHypotheses);
243  loop->Get(mctrajpoints);
244 
245  DVector3 pthrown = trk.pthrown;
246 
247  bool locFoundFlag = false;
248  double fom_best = 1.0E8;
249 
250  trk.num_timebased = locChargedTrackHypotheses.size();
251 
252  // Loop over found/fit tracks
253  for(unsigned int j=0; j<locChargedTrackHypotheses.size(); j++){
254  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrackHypotheses[j];
255  const DTrackTimeBased *locTimeBasedTrack = locChargedTrackHypothesis->dTrackTimeBased;
256 
257  // Get DTrackWireBased and DTrackCandidate objects for this DTrackTimeBased
258  vector<const DTrackWireBased*> tracks;
259  locTimeBasedTrack->Get(tracks);
260  const DTrackWireBased *track = tracks.size()==1 ? tracks[0]:NULL;
261  vector<const DTrackCandidate*> trackcandidates;
262  if(track)track->Get(trackcandidates);
263 
264  // Copy momentum vectors to convenient local variables
265  DVector3 pfit = locTimeBasedTrack->momentum();
266 
267  // Calculate residuals from momentum parameters from DTrackTimeBased
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;
271 
272  // Formulate a figure of merit to decide if this fit track is closer to
273  // the thrown track than the best one we found so far. We hardwire
274  // dpt/pt=2%, dtheta=20mrad and dphi=20mrad for now.
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);
276  if(fom<fom_best){
277  fom_best = fom;
278  locFoundFlag = true;
279  trk.PID_hypothesized = int(locChargedTrackHypothesis->dPID);
280  trk.q_hypothesized = locChargedTrackHypothesis->charge();
281  trk.FOM_hypothesized = locChargedTrackHypothesis->dFOM;
282 
283  trk.dTrackReconstructedFlag_TimeBased = true;
284 
285  trk.q_timebased = locTimeBasedTrack->charge();
286 
287  trk.PID_timebased = int(locPIDAlgorithm->IDTrack(locTimeBasedTrack->charge(), locTimeBasedTrack->mass()));
288 
289  trk.pfit = pfit;
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;
295 
296  // Get Nstereo, Ncdc, and Nfdc
297  vector<const DCDCTrackHit*> cdchits;
298  locTimeBasedTrack->Get(cdchits);
299  trk.Nstereo = 0;
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();
305 
306  // Get the number LR signs for all hits used on this track. We have to
307  // do this for the thrown track here so we get the list for the same hits
308  // used in fitting this track.
309  vector<int> LRthrown;
310  vector<int> LRfit;
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);
314  if(use_rt_thrown){
315  FindLR(wires, rt_thrown, LRthrown);
316  }else{
317  FindLR(wires, mctrajpoints, LRthrown);
318  }
319  FindLR(wires, locTimeBasedTrack->rt, LRfit);
320 
321  // Make sure the number of entries is the same for both LR vectors
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;
324  continue;
325  }
326 
327  // count total number of incorrect LR choices and how many are stereo
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;
331  trk.NLR_bad++;
332  bool is_stereo = (wires[k]->udir.Theta()*57.3)>2.0;
333  if(is_stereo)trk.NLR_bad_stereo++;
334  }
335  }
336  }
337  return locFoundFlag;
338 }
339 
340 bool DEventProcessor_trackeff_hists2::Search_WireBasedTracks(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown){
341  vector<const DMCTrajectoryPoint*> mctrajpoints;
342  vector<const DCDCTrackHit*> cdctrackhits;
343  vector<const DFDCHit*> fdchits;
344  vector<const DTrackWireBased*> tracks;
345 
346  loop->Get(cdctrackhits);
347  loop->Get(fdchits);
348  loop->Get(tracks);
349  loop->Get(mctrajpoints);
350 
351  DVector3 pthrown = trk.pthrown;
352  double fom_best = 1.0E8;
353 
354  trk.num_wirebased = tracks.size();
355 
356  bool locFoundFlag = false;
357  // Loop over found/fit tracks
358  for(unsigned int j=0; j<tracks.size(); j++){
359  const DTrackWireBased *track = tracks[j];
360 
361  // Copy momentum vectors to convenient local variables
362  DVector3 pfit_wire = track->momentum();
363 
364  // Calculate residuals from momentum parameters from DTrackTimeBased
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;
368 
369  // Formulate a figure of merit to decide if this fit track is closer to
370  // the thrown track than the best one we found so far. We hardwire
371  // dpt/pt=2%, dtheta=20mrad and dphi=20mrad for now.
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);
373  if(fom<fom_best){
374  fom_best = fom;
375  locFoundFlag = true;
376  trk.dTrackReconstructedFlag_WireBased = true;
377 
378  trk.q_wirebased = track->charge();
379 
380  trk.PID_wirebased = int(locPIDAlgorithm->IDTrack(track->charge(), track->mass()));
381 
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;
388 
389  // Get Nstereo, Ncdc, and Nfdc
390  vector<const DCDCTrackHit*> cdchits;
391  track->Get(cdchits);
392  trk.Nstereo = 0;
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;
396  track->Get(fdchits);
397  trk.Nfdc = fdchits.size();
398 
399  // Get the number LR signs for all hits used on this track. We have to
400  // do this for the thrown track here so we get the list for the same hits
401  // used in fitting this track.
402  vector<int> LRthrown;
403  vector<int> LRfit;
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);
407  if(use_rt_thrown){
408  FindLR(wires, rt_thrown, LRthrown);
409  }else{
410  FindLR(wires, mctrajpoints, LRthrown);
411  }
412  FindLR(wires, track->rt, LRfit);
413 
414  // Make sure the number of entries is the same for both LR vectors
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;
417  continue;
418  }
419 
420  // count total number of incorrect LR choices and how many are stereo
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;
424  trk.NLR_bad++;
425  bool is_stereo = (wires[k]->udir.Theta()*57.3)>2.0;
426  if(is_stereo)trk.NLR_bad_stereo++;
427  }
428  }
429  }
430  return locFoundFlag;
431 }
432 
433 bool DEventProcessor_trackeff_hists2::Search_TrackCandidates(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown){
434  vector<const DMCTrajectoryPoint*> mctrajpoints;
435  vector<const DCDCTrackHit*> cdctrackhits;
436  vector<const DFDCHit*> fdchits;
437  vector<const DTrackCandidate*> trackcandidates;
438 
439  loop->Get(cdctrackhits);
440  loop->Get(fdchits);
441  loop->Get(mctrajpoints);
442  loop->Get(trackcandidates);
443 
444  DVector3 pthrown = trk.pthrown;
445  double fom_best = 1.0E8;
446 
447  trk.num_candidates = trackcandidates.size();
448 
449  bool locFoundFlag = false;
450  // Loop over found/fit tracks
451  for(unsigned int j=0; j<trackcandidates.size(); j++){
452  const DTrackCandidate *trackcandidate = trackcandidates[j];
453 
454  // Copy momentum vectors to convenient local variables
455  DVector3 pcan = trackcandidate->momentum();
456 
457  // Calculate residuals from momentum parameters from DTrackTimeBased
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;
461 
462  // Formulate a figure of merit to decide if this fit track is closer to
463  // the thrown track than the best one we found so far. We hardwire
464  // dpt/pt=2%, dtheta=20mrad and dphi=20mrad for now.
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);
466  if(fom<fom_best){
467  fom_best = fom;
468  locFoundFlag = true;
469  trk.dTrackReconstructedFlag_Candidate = true;
470  trk.q_candidate = trackcandidate->charge();
471  trk.PID_candidate = int(locPIDAlgorithm->IDTrack(trackcandidate->charge(), trackcandidate->mass()));
472  trk.pcan = pcan;
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;
476 
477  // Get Nstereo, Ncdc, and Nfdc
478  vector<const DCDCTrackHit*> cdchits;
479  trackcandidate->Get(cdchits);
480  trk.Nstereo = 0;
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();
486  }
487  }
488  return locFoundFlag;
489 }
490 
491 //------------------
492 // isReconstructable
493 //------------------
494 bool DEventProcessor_trackeff_hists2::isReconstructable(const DMCThrown *mcthrown, vector<const DMCTrajectoryPoint*> &mctrajpoints)
495 {
496  /// In order to test the efficiency of the finder/fitter, we must first
497  /// determine whether a track is "reconstructible" or not. (See
498  /// COMPASS-Note 2004-1 section 5.1 and
499  /// Mankel Rep. Prog. Phys. 67 (2004) 553-622 section 2.5.1)
500  ///
501  /// We do this here by checking if the "death" point of the
502  /// track is inside or outside of some inner volume of the detector.
503  /// Specifically, if the death point is outside of a cylinder
504  /// of radius 45cm and z-extent covering from the upstream end
505  /// of the CDC to the downstream end of the 2nd FDC package.
506  /// This ensures that it either passes through the first 2 layers
507  /// of the outermost axial superlayer of the CDC or, the first
508  /// 2 packages of the FDC.
509  ///
510  /// Note that the beam hole is not considered here. It is assumed
511  /// that efficiency plots will be made either with a cut on theta
512  /// or plotted against theta in order to accomodate that area.
513 
514  const DMCTrajectoryPoint *mctraj_first = NULL;
515  const DMCTrajectoryPoint *mctraj_last = NULL;
516  for(unsigned int i=0; i<mctrajpoints.size(); i++){
517  const DMCTrajectoryPoint *mctraj = mctrajpoints[i];
518  if(mctraj->track != mcthrown->myid)continue;
519  if(mctraj->primary_track != mcthrown->myid)continue; // redundant?
520  if(mctraj_first==NULL)mctraj_first = mctraj;
521  mctraj_last = mctraj;
522  }
523 
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;
529  }
530 
531  return false;
532 }
533 
534 //------------------
535 // FindLR
536 //------------------
537 void DEventProcessor_trackeff_hists2::FindLR(vector<const DCoordinateSystem*> &wires, const DReferenceTrajectory *crt, vector<int> &LRhits)
538 {
539  /// Fill the vector LRhits with +1 or -1 values indicating the side of each wire in the "wires"
540  /// vector the given reference trajectory passed on.
541 
542  LRhits.clear();
543 
544  // This first bit is shameful. In order to use the DistToRT methods of the reference trajectory,
545  // we have to cast away the const qualifier since those methods must modify the object
546  DReferenceTrajectory *rt = const_cast<DReferenceTrajectory*>(crt);
547  for(unsigned int i=0; i<wires.size(); i++){
548  const DCoordinateSystem *wire = wires[i];
549 
550  DVector3 pos_doca, mom_doca;
551  rt->DistToRT(wire);
552  rt->GetLastDOCAPoint(pos_doca, mom_doca);
553  DVector3 shift = wire->udir.Cross(mom_doca);
554  double u = rt->GetLastDistAlongWire();
555  DVector3 pos_wire = wire->origin + u*wire->udir;
556  DVector3 pos_diff = pos_doca-pos_wire;
557 
558  LRhits.push_back(shift.Dot(pos_diff)<0.0 ? -1:1);
559  }
560 }
561 
562 //------------------
563 // FindLR
564 //------------------
565 void DEventProcessor_trackeff_hists2::FindLR(vector<const DCoordinateSystem*> &wires, vector<const DMCTrajectoryPoint*> &trajpoints, vector<int> &LRhits)
566 {
567  /// Fill the vector LRhits with +1 or -1 values indicating the side of each wire in the "wires"
568  /// vector the particle swum by GEANT passed on according to the closest DMCTrajectoryPoint.
569 
570 }
571 
Definition: track.h:16
void trkeff(int pid)
Definition: trkeff.C:1
DApplication * dapp
bool isReconstructable(const DMCThrown *mcthrown, vector< const DMCTrajectoryPoint * > &mctrajpoints)
float chisq
Chi-squared for the track (not chisq/dof!)
double rmax
float chisq
Chi-squared for the track (not chisq/dof!)
TVector3 DVector3
Definition: DVector3.h:14
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
DVector3 GetLastDOCAPoint(void) const
const DVector3 & position(void) const
double GetLastDistAlongWire(void) const
double rmin
JApplication * japp
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
Definition: DGeometry.cc:1425
DMagneticFieldMap * GetBfield(void) const
Definition: DGeometry.cc:61
double DistToRT(double x, double y, double z) const
InitPlugin_t InitPlugin
bool Search_TrackCandidates(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown)
DGeometry * GetDGeometry(unsigned int run_number)
#define _DBG_
Definition: HDEVIO.h:12
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
&lt;A href=&quot;index.html#legend&quot;&gt; &lt;IMG src=&quot;CORE.png&quot; width=&quot;100&quot;&gt; &lt;/A&gt;
int Ndof
Number of degrees of freedom in the fit.
double sqrt(double)
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)
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
Definition: DGeometry.cc:1497
TDirectory * dir
Definition: bcal_hist_eff.C:31
int type
GEANT particle ID.
Definition: DMCThrown.h:20
bool Search_ChargedTrackHypotheses(JEventLoop *loop, uint64_t eventnumber, const DMCThrown *mcthrown)
union @6::@8 u
#define DEBUG
Definition: tw_corr.C:41
int myid
id of this particle from original generator
Definition: DMCThrown.h:22
double mass(void) const