Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DTrackTimeBased_factory.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DTrackTimeBased_factory.cc
4 // Created: Thu Sep 4 14:02:44 EDT 2008
5 // Creator: davidl (on Darwin harriet.jlab.org 8.11.1 i386)
6 //
7 
8 
9 #include <iostream>
10 #include <iomanip>
11 #include <set>
12 #include <mutex>
13 #include <TMath.h>
14 using namespace std;
15 
16 #define TOF_SIGMA 0.080 // TOF resolution in ns
17 
18 #include <TROOT.h>
19 
23 #include <TRACKING/DTrackFitter.h>
25 #include <TRACKING/DMCTrackHit.h>
26 #include <SplitString.h>
28 #include <deque>
29 
30 using namespace jana;
31 
32 // Routine for sorting start times
35  return (a.system>b.system);
36 }
37 
39  if (a->candidateid==b->candidateid) return a->mass()<b->mass();
40  return a->candidateid<b->candidateid;
41 }
42 
43 
44 // count_common_members
45 //------------------
46 template<typename T>
47 static unsigned int count_common_members(vector<T> &a, vector<T> &b)
48 {
49  unsigned int n=0;
50  for(unsigned int i=0; i<a.size(); i++){
51  for(unsigned int j=0; j<b.size(); j++){
52  if(a[i]==b[j])n++;
53  }
54  }
55 
56  return n;
57 }
58 
59 
60 
61 //------------------
62 // init
63 //------------------
65 {
66  fitter = NULL;
67 
68  DEBUG_HISTS = false;
69  //DEBUG_HISTS = true;
70  DEBUG_LEVEL = 0;
71 
72  USE_HITS_FROM_WIREBASED_FIT=false;
73  gPARMS->SetDefaultParameter("TRKFIT:USE_HITS_FROM_WIREBASED_FIT",
74  USE_HITS_FROM_WIREBASED_FIT);
75  INSERT_MISSING_HYPOTHESES=true;
76  gPARMS->SetDefaultParameter("TRKFIT:INSERT_MISSING_HYPOTHESES",
77  INSERT_MISSING_HYPOTHESES);
78 
79  gPARMS->SetDefaultParameter("TRKFIT:DEBUG_HISTS",DEBUG_HISTS);
80  gPARMS->SetDefaultParameter("TRKFIT:DEBUG_LEVEL",DEBUG_LEVEL);
81 
82  vector<int> hypotheses;
83  hypotheses.push_back(Positron);
84  hypotheses.push_back(PiPlus);
85  hypotheses.push_back(KPlus);
86  hypotheses.push_back(Proton);
87  hypotheses.push_back(Electron);
88  hypotheses.push_back(PiMinus);
89  hypotheses.push_back(KMinus);
90  hypotheses.push_back(AntiProton);
91 
92  ostringstream locMassStream;
93  for(size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i)
94  {
95  locMassStream << hypotheses[loc_i];
96  if(loc_i != (hypotheses.size() - 1))
97  locMassStream << ",";
98  }
99 
100  string HYPOTHESES = locMassStream.str();
101  gPARMS->SetDefaultParameter("TRKFIT:HYPOTHESES", HYPOTHESES);
102 
103  // Parse MASS_HYPOTHESES strings to make list of masses to try
104  hypotheses.clear();
105  SplitString(HYPOTHESES, hypotheses, ",");
106  for(size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i)
107  {
108  if(ParticleCharge(Particle_t(hypotheses[loc_i])) > 0)
109  mass_hypotheses_positive.push_back(hypotheses[loc_i]);
110  else if(ParticleCharge(Particle_t(hypotheses[loc_i])) < 0)
111  mass_hypotheses_negative.push_back(hypotheses[loc_i]);
112  }
113 
114  if(mass_hypotheses_positive.empty()){
115  static once_flag pwarn_flag;
116  call_once(pwarn_flag, [](){
117  jout << endl;
118  jout << "############# WARNING !! ################ " <<endl;
119  jout << "There are no mass hypotheses for positive tracks!" << endl;
120  jout << "Be SURE this is what you really want!" << endl;
121  jout << "######################################### " <<endl;
122  jout << endl;
123  });
124  }
125  if(mass_hypotheses_negative.empty()){
126  static once_flag nwarn_flag;
127  call_once(nwarn_flag, [](){
128  jout << endl;
129  jout << "############# WARNING !! ################ " <<endl;
130  jout << "There are no mass hypotheses for negative tracks!" << endl;
131  jout << "Be SURE this is what you really want!" << endl;
132  jout << "######################################### " <<endl;
133  jout << endl;
134  });
135  }
136 
137  mNumHypPlus=mass_hypotheses_positive.size();
138  mNumHypMinus=mass_hypotheses_negative.size();
139 
140  // Forces correct particle id (when available)
141  PID_FORCE_TRUTH = false;
142  gPARMS->SetDefaultParameter("TRKFIT:PID_FORCE_TRUTH", PID_FORCE_TRUTH);
143 
144  USE_SC_TIME=true;
145  gPARMS->SetDefaultParameter("TRKFIT:USE_SC_TIME",USE_SC_TIME);
146 
147  USE_TOF_TIME=true;
148  gPARMS->SetDefaultParameter("TRKFIT:USE_TOF_TIME",USE_TOF_TIME);
149 
150  USE_FCAL_TIME=true;
151  gPARMS->SetDefaultParameter("TRKFIT:USE_FCAL_TIME",USE_FCAL_TIME);
152 
153  USE_BCAL_TIME=true;
154  gPARMS->SetDefaultParameter("TRKFIT:USE_BCAL_TIME",USE_BCAL_TIME);
155 
156  return NOERROR;
157 }
158 
159 //------------------
160 // brun
161 //------------------
162 jerror_t DTrackTimeBased_factory::brun(jana::JEventLoop *loop, int32_t runnumber)
163 {
164  // Get the geometry
165  DApplication* dapp=dynamic_cast<DApplication*>(loop->GetJApplication());
166  geom = dapp->GetDGeometry(runnumber);
167  // Check for magnetic field
168  dIsNoFieldFlag = (dynamic_cast<const DMagneticFieldMapNoField*>(dapp->GetBfield(runnumber)) != NULL);
169 
170  if(dIsNoFieldFlag){
171  //Setting this flag makes it so that JANA does not delete the objects in
172  //_data. This factory will manage this memory.
173  //This is all of these pointers are just copied from the "StraightLine"
174  //factory, and should not be re-deleted.
175  SetFactoryFlag(NOT_OBJECT_OWNER);
176  }
177  else{
178  ClearFactoryFlag(NOT_OBJECT_OWNER); //This factory will create it's own obje
179  }
180 
181  // Get pointer to TrackFitter object that actually fits a track
182  vector<const DTrackFitter *> fitters;
183  loop->Get(fitters);
184  if(fitters.size()<1){
185  _DBG_<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
186  return RESOURCE_UNAVAILABLE;
187  }
188 
189  // Drop the const qualifier from the DTrackFitter pointer (I'm surely going to hell for this!)
190  fitter = const_cast<DTrackFitter*>(fitters[0]);
191 
192  // Warn user if something happened that caused us NOT to get a fitter object pointer
193  if(!fitter){
194  _DBG_<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
195  return RESOURCE_UNAVAILABLE;
196  }
197 
198  // Get the particle ID algorithms
199  vector<const DParticleID *> pid_algorithms;
200  loop->Get(pid_algorithms);
201  if(pid_algorithms.size()<1){
202  _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
203  return RESOURCE_UNAVAILABLE;
204  }
205 
206  pid_algorithm = pid_algorithms[0];
207 
208  // Warn user if something happened that caused us NOT to get a pid_algorithm object pointer
209  if(!pid_algorithm){
210  _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
211  return RESOURCE_UNAVAILABLE;
212  }
213 
214 
215  if(DEBUG_HISTS){
216  dapp->Lock();
217 
218  // Histograms may already exist. (Another thread may have created them)
219  // Try and get pointers to the existing ones.
220  fom_chi2_trk = (TH1F*)gROOT->FindObject("fom_chi2_trk");
221  fom = (TH1F*)gROOT->FindObject("fom");
222  hitMatchFOM = (TH1F*)gROOT->FindObject("hitMatchFOM");
223  chi2_trk_mom = (TH2F*)gROOT->FindObject("chi2_trk_mom");
224 
225  if(!fom_chi2_trk)fom_chi2_trk = new TH1F("fom_chi2_trk","PID FOM: #chi^{2}/Ndf from tracking", 1000, 0.0, 100.0);
226  if(!fom)fom = new TH1F("fom","Combined PID FOM", 1000, 0.0, 1.01);
227  if(!hitMatchFOM)hitMatchFOM = new TH1F("hitMatchFOM","Total Fraction of Hit Matches", 101, 0.0, 1.01);
228  if(!chi2_trk_mom)chi2_trk_mom = new TH2F("chi2_trk_mom","Track #chi^{2}/Ndf versus Kinematic #chi^{2}/Ndf", 1000, 0.0, 100.0, 1000, 0.,100.);
229 
230 
231  Hstart_time=(TH2F*)gROOT->FindObject("Hstart_time");
232  if (!Hstart_time) Hstart_time=new TH2F("Hstart_time",
233  "vertex time source vs. time",
234  300,-50,50,9,-0.5,8.5);
235 
236  dapp->Unlock();
237 
238  }
239 
240  return NOERROR;
241 }
242 
243 //------------------
244 // evnt
245 //------------------
246 jerror_t DTrackTimeBased_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
247 {
248  // Save event number to help with debugging
249  myevt=eventnumber;
250  if(!fitter)return NOERROR;
251 
252  if(dIsNoFieldFlag){
253  //Clear previous objects:
254  //JANA doesn't do it because NOT_OBJECT_OWNER was set
255  //It DID delete them though, in the "StraightLine" factory
256  _data.clear();
257 
258  vector<const DTrackTimeBased*> locTimeBasedTracks;
259  loop->Get(locTimeBasedTracks, "StraightLine");
260  for(size_t loc_i = 0; loc_i < locTimeBasedTracks.size(); ++loc_i)
261  _data.push_back(const_cast<DTrackTimeBased*>(locTimeBasedTracks[loc_i]));
262  return NOERROR;
263  }
264 
265  // Get candidates and hits
266  vector<const DTrackWireBased*> tracks;
267  loop->Get(tracks);
268  if (tracks.size()==0) return NOERROR;
269 
270  // get start counter hits
271  vector<const DSCHit*>sc_hits;
272  if (USE_SC_TIME){
273  loop->Get(sc_hits);
274  }
275 
276  // Get TOF points
277  vector<const DTOFPoint*> tof_points;
278  if (USE_TOF_TIME){
279  loop->Get(tof_points);
280  }
281 
282  // Get BCAL and FCAL showers
283  vector<const DBCALShower*>bcal_showers;
284  if (USE_BCAL_TIME){
285  loop->Get(bcal_showers);
286  }
287  vector<const DFCALShower*>fcal_showers;
288  if (USE_FCAL_TIME){
289  loop->Get(fcal_showers);
290  }
291 
292  vector<const DMCThrown*> mcthrowns;
293  loop->Get(mcthrowns, "FinalState");
294 
295  // Loop over candidates
296  for(unsigned int i=0; i<tracks.size(); i++){
297  const DTrackWireBased *track = tracks[i];
298 
299  unsigned int num=_data.size();
300 
301  // Create vector of start times from various sources
302  vector<DTrackTimeBased::DStartTime_t>start_times;
303  CreateStartTimeList(track,sc_hits,tof_points,bcal_showers,fcal_showers,start_times);
304 
305  // Fit the track
306  DoFit(track,start_times,loop,track->mass());
307 
308  //_DBG_<< "eventnumber: " << eventnumber << endl;
309  if (PID_FORCE_TRUTH && _data.size()>num) {
310  // Add figure-of-merit based on difference between thrown and reconstructed momentum
311  // if more than half of the track's hits match MC truth hits and also (charge,mass)
312  // match; add FOM=0 otherwise
313  _data[_data.size()-1]->FOM=GetTruthMatchingFOM(i,_data[_data.size()-1],
314  mcthrowns);
315  }
316  } // loop over track candidates
317 
318  // Filter out duplicate tracks
319  FilterDuplicates();
320 
321  // Fill in track data for missing hypotheses
322  if (INSERT_MISSING_HYPOTHESES){
323  InsertMissingHypotheses(loop);
324  }
325 
326  // Set MC Hit-matching information
327  for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
328  {
329  if(!mcthrowns.empty())
330  {
331  double hitFraction;
332  int thrownIndex = GetThrownIndex(mcthrowns, (DKinematicData*)_data[loc_i], hitFraction);
333  _data[loc_i]->dMCThrownMatchMyID = thrownIndex;
334  _data[loc_i]->dNumHitsMatchedToThrown = int(hitFraction * float(Get_NumTrackHits(_data[loc_i])) + 0.01); // + 0.01 so that it rounds down properly
335  }
336  else
337  {
338  _data[loc_i]->dMCThrownMatchMyID = -1;
339  _data[loc_i]->dNumHitsMatchedToThrown = 0;
340  }
341  }
342 
343  return NOERROR;
344 }
345 
346 //------------------
347 // erun
348 //------------------
350 {
351  return NOERROR;
352 }
353 
354 //------------------
355 // fini
356 //------------------
358 {
359  return NOERROR;
360 }
361 
362 //------------------
363 // FilterDuplicates
364 //------------------
366 {
367  /// Look through all current DTrackTimeBased objects and remove any
368  /// that have most of their hits in common with another track
369 
370  if(_data.size()==0)return;
371 
372  if(DEBUG_LEVEL>2)_DBG_<<"Looking for clones of time-based tracks ..."<<endl;
373  // We want to remove duplicate tracks corresponding to actual particles,
374  // not just duplicate fitted tracks for certain mass hypotheses -- this
375  // is partly because at a later stage the holes in the list of mass
376  // hypotheses are filled in, thereby spoiling the whole point of this
377  // part of the code!
378  // Keep track of pairs of candididate id's, one for which we want to
379  // keep all the results of fitting with different mass hypotheses,
380  // the other for which we want to delete all the results of fitting.
381  // We need both vectors to take into account potential ambiguities:
382  // for one mass hypothesis starting with one candidate may be "better"
383  // than starting with a second clone candidate, whereas for a second
384  // mass hypothesis, the opposite may be true.
385  vector<unsigned int> candidates_to_keep;
386  vector<unsigned int> candidates_to_delete;
387  for(unsigned int i=0; i<_data.size()-1; i++){
388  DTrackTimeBased *dtrack1 = _data[i];
389 
390  vector<const DCDCTrackHit*> cdchits1;
391  vector<const DFDCPseudo*> fdchits1;
392  dtrack1->Get(cdchits1);
393  dtrack1->Get(fdchits1);
394  // Total number of hits in this candidate
395  unsigned int num_cdc1=cdchits1.size();
396  unsigned int num_fdc1=fdchits1.size();
397  unsigned int total1 = num_cdc1+num_fdc1;
398 
399  JObject::oid_t cand1=dtrack1->candidateid;
400  for(unsigned int j=i+1; j<_data.size(); j++){
401  DTrackTimeBased *dtrack2 = _data[j];
402  if (dtrack2->candidateid==cand1) continue;
403 
404  vector<const DCDCTrackHit*> cdchits2;
405  vector<const DFDCPseudo*> fdchits2;
406  dtrack2->Get(cdchits2);
407  dtrack2->Get(fdchits2);
408 
409  // Total number of hits in this candidate
410  unsigned int num_cdc2=cdchits2.size();
411  unsigned int num_fdc2=fdchits2.size();
412  unsigned int total2 = num_cdc2+num_fdc2;
413 
414  // Count number of cdc and fdc hits in common
415  unsigned int Ncdc = count_common_members(cdchits1, cdchits2);
416  unsigned int Nfdc = count_common_members(fdchits1, fdchits2);
417 
418  if(DEBUG_LEVEL>3){
419  _DBG_<<"cand1:"<<cand1<<" cand2:"<<dtrack2->candidateid<<endl;
420  _DBG_<<" Ncdc="<<Ncdc<<" num_cdc1="<<num_cdc1<<" num_cdc2="<<num_cdc2<<endl;
421  _DBG_<<" Nfdc="<<Nfdc<<" num_fdc1="<<num_fdc1<<" num_fdc2="<<num_fdc2<<endl;
422  }
423  unsigned int total = Ncdc + Nfdc;
424  // If the tracks share at most one hit, consider them
425  // to be separate tracks
426  if (total<=1) continue;
427 
428  // Deal with the case where there are cdc hits in
429  // common between the tracks but there were no fdc
430  // hits used in one of the tracks.
431  if (Ncdc>0 && (num_fdc1>0 || num_fdc2>0)
432  && (num_fdc1*num_fdc2)==0) continue;
433 
434  // Deal with the case where there are fdc hits in
435  // common between the tracks but no cdc hits used in
436  // one of the tracks.
437  if (Nfdc>0 && (num_cdc1>0 || num_cdc2>0)
438  && (num_cdc1*num_cdc2)==0) continue;
439 
440  // Look for tracks with many common hits in the CDC
441  if (num_cdc1>0 && num_cdc2>0){
442  if (double(Ncdc)/double(num_cdc1)<0.9) continue;
443  if (double(Ncdc)/double(num_cdc2)<0.9) continue;
444  }
445  // Look for tracks with many common hits in the FDC
446  if (num_fdc1>0 && num_fdc2>0){
447  if (double(Nfdc)/double(num_fdc1)<0.9) continue;
448  if (double(Nfdc)/double(num_fdc2)<0.9) continue;
449  }
450 
451  if(total1<total2){
452  candidates_to_delete.push_back(cand1);
453  candidates_to_keep.push_back(dtrack2->candidateid);
454  }else if(total2<total1){
455  candidates_to_delete.push_back(dtrack2->candidateid);
456  candidates_to_keep.push_back(cand1);
457  }else if(dtrack1->FOM > dtrack2->FOM){
458  candidates_to_delete.push_back(dtrack2->candidateid);
459  candidates_to_keep.push_back(cand1);
460  }else{
461  candidates_to_delete.push_back(cand1);
462  candidates_to_keep.push_back(dtrack2->candidateid);
463  }
464  }
465  }
466 
467  if(DEBUG_LEVEL>2)_DBG_<<"Found "<<candidates_to_delete.size()<<" time-based clones"<<endl;
468 
469 
470  // Return now if we're keeping everyone
471  if(candidates_to_delete.size()==0)return;
472 
473  // Deal with the ambiguity problem mentioned above
474  for (unsigned int i=0;i<candidates_to_keep.size();i++){
475  for (unsigned int j=0;j<candidates_to_delete.size();j++){
476  if (candidates_to_keep[i]==candidates_to_delete[j]){
477  candidates_to_delete.erase(candidates_to_delete.begin()+j);
478  break;
479  }
480  }
481 
482  }
483 
484  // Copy pointers that we want to keep to a new container and delete
485  // the clone objects
486  vector<DTrackTimeBased*> new_data;
487  sort(_data.begin(),_data.end(),DTrackTimeBased_cmp);
488  for (unsigned int i=0;i<_data.size();i++){
489  bool keep_track=true;
490  for (unsigned int j=0;j<candidates_to_delete.size();j++){
491  if (_data[i]->candidateid==candidates_to_delete[j]){
492  keep_track=false;
493  if(DEBUG_LEVEL>1){
494  _DBG_<<"Deleting clone time-based fitted result "<<i
495  << " in event " << myevt << endl;
496  }
497  break;
498  }
499  }
500  if (keep_track){
501  new_data.push_back(_data[i]);
502  }
503  else delete _data[i];
504  }
505  _data = new_data;
506 }
507 
508 // Returns a FOM based on difference between thrown and reconstructed momentum if track matches MC truth information,
509 // returns a FOM=0 otherwise;
510 // a match requires identical masses and charges, and that more than half of the track's hits match the truth hits
511 double DTrackTimeBased_factory::GetTruthMatchingFOM(int trackIndex,DTrackTimeBased *track,vector<const DMCThrown*>mcthrowns) {
512  bool match=false;
513 
514  DLorentzVector fourMom = track->lorentzMomentum();
515  //DLorentzVector gen_fourMom[mcthrowns.size()];
516  vector<DLorentzVector> gen_fourMom(mcthrowns.size());
517  for(unsigned int i=0; i<mcthrowns.size(); i++){
518  gen_fourMom[i] = mcthrowns[i]->lorentzMomentum();
519  }
520 
521  // Get info for thrown track
522  double f = 0.;
523  int thrownIndex = GetThrownIndex(mcthrowns, track, f);
524  if(thrownIndex<=0 || f<=0.5) return 0.;
525 
526  double delta_pt_over_pt = (fourMom.Pt()-gen_fourMom[thrownIndex-1].Pt())/gen_fourMom[thrownIndex-1].Pt();
527  double delta_theta = (fourMom.Theta()-gen_fourMom[thrownIndex-1].Theta())*1000.0; // in milliradians
528  double delta_phi = (fourMom.Phi()-gen_fourMom[thrownIndex-1].Phi())*1000.0; // in milliradians
529  double chisq = pow(delta_pt_over_pt/0.04, 2.0) + pow(delta_theta/20.0, 2.0) + pow(delta_phi/20.0, 2.0);
530 
531  if (fabs(track->mass()-mcthrowns[thrownIndex-1]->mass())<0.01 && track->charge()==mcthrowns[thrownIndex-1]->charge())
532  match = true;
533 
534  double trk_chi2=track->chisq;
535  unsigned int ndof=track->Ndof;
536 
537  if(DEBUG_HISTS&&match){
538  fom_chi2_trk->Fill(track->chisq);
539  chi2_trk_mom->Fill(chisq/3.,trk_chi2/ndof);
540  fom->Fill(TMath::Prob(chisq,3));
541  }
542 
543  /*_DBG_ << "f: " << f << endl;
544  _DBG_ << "trk_chi2: " << trk_chi2 << endl;
545  _DBG_ << "ndof: " << ndof << endl;
546  _DBG_ << "throwncharge: " << mcthrowns[thrownIndex-1]->charge() << endl;
547  _DBG_ << "trackcharge: " << track->charge() << endl;
548  _DBG_ << "chargediff: " << fabs(track->charge()-mcthrowns[thrownIndex-1]->charge()) << endl;
549  _DBG_ << "thrownmass: " << mcthrowns[thrownIndex-1]->mass() << endl;
550  _DBG_ << "trackmass: " << track->mass() << endl;
551  _DBG_ << "massdiff: " << fabs(track->mass()-mcthrowns[thrownIndex-1]->mass()) << endl;
552  _DBG_ << "chisq: " << chisq << endl;
553  _DBG_ << "match?: " << match << endl;
554  _DBG_ << "thrownIndex: " << thrownIndex << " trackIndex: " << trackIndex << endl;
555  _DBG_<< "track " << setprecision(4) << "Px: " << fourMom.Px() << " Py: " << fourMom.Py() << " Pz: " << fourMom.Pz() << " E: " << fourMom.E() << " M: " << fourMom.M() << " pt: " << fourMom.Pt() << " theta: " << fourMom.Theta() << " phi: " << fourMom.Phi() << endl;
556  _DBG_<< "thrown " << setprecision(4) << "Px: " << gen_fourMom[thrownIndex-1].Px() << " Py: " << gen_fourMom[thrownIndex-1].Py() << " Pz: " << gen_fourMom[thrownIndex-1].Pz() << " E: " << gen_fourMom[thrownIndex-1].E() << " M: " << gen_fourMom[thrownIndex-1].M() << " pt: " << gen_fourMom[thrownIndex-1].Pt() << " theta: " << gen_fourMom[thrownIndex-1].Theta() << " phi: " << gen_fourMom[thrownIndex-1].Phi() << endl;*/
557 
558  return (match) ? TMath::Prob(chisq,3) : 0.0;
559 }
560 
561 //------------------
562 // GetThrownIndex
563 //------------------
564 int DTrackTimeBased_factory::GetThrownIndex(vector<const DMCThrown*>& locMCThrowns, const DKinematicData *kd, double &f)
565 {
566  // The DKinematicData object should be a DTrackCandidate, DTrackWireBased, or DParticle which
567  // has associated objects for the hits
568  vector<const DCDCTrackHit*> cdctrackhits;
569  kd->Get(cdctrackhits);
570  vector<const DFDCPseudo*> fdcpseudos;
571  kd->Get(fdcpseudos);
572 
573  int locTotalNumHits = cdctrackhits.size() + fdcpseudos.size();
574  if(locTotalNumHits == 0)
575  {
576  f = 0;
577  return -1;
578  }
579 
580  // The track number is buried in the truth hit objects of type DMCTrackHit. These should be
581  // associated objects for the individual hit objects. We need to loop through them and
582  // keep track of how many hits for each track number we find
583 
584  map<int, int> locHitMatches; //first int is MC my id, second is num hits
585  for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
586  {
587  if(fabs(locMCThrowns[loc_i]->charge()) > 0.9)
588  locHitMatches[locMCThrowns[loc_i]->myid] = 0;
589  }
590 
591  // CDC hits
592  for(size_t loc_i = 0; loc_i < cdctrackhits.size(); ++loc_i)
593  {
594  const DCDCHit* locCDCHit = NULL;
595  cdctrackhits[loc_i]->GetSingle(locCDCHit);
596  vector<const DCDCHit*> locTruthCDCHits;
597  locCDCHit->Get(locTruthCDCHits);
598  if(locTruthCDCHits.empty()) continue;
599 
600  int itrack = locTruthCDCHits[0]->itrack;
601  if(locHitMatches.find(itrack) == locHitMatches.end())
602  continue;
603  ++locHitMatches[itrack];
604  }
605 
606  // FDC hits
607  for(size_t loc_i = 0; loc_i < fdcpseudos.size(); ++loc_i)
608  {
609  vector<const DMCTrackHit*> mctrackhits;
610  fdcpseudos[loc_i]->Get(mctrackhits);
611  if(mctrackhits.empty())
612  continue;
613  if(!mctrackhits[0]->primary)
614  continue;
615 
616  int itrack = mctrackhits[0]->itrack;
617  if(locHitMatches.find(itrack) == locHitMatches.end())
618  continue;
619  ++locHitMatches[itrack];
620  }
621 
622  // Find DMCThrown::myid with most wires hit
623  map<int, int>::iterator locIterator = locHitMatches.begin();
624  int locBestMyID = -1;
625  int locBestNumHits = 0;
626  for(; locIterator != locHitMatches.end(); ++locIterator)
627  {
628  if(locIterator->second <= locBestNumHits)
629  continue;
630  locBestNumHits = locIterator->second;
631  locBestMyID = locIterator->first;
632  }
633 
634  // total fraction of reconstructed hits that match truth hits
635  f = 1.0*locBestNumHits/locTotalNumHits;
636  if(DEBUG_HISTS)hitMatchFOM->Fill(f);
637 
638  return locBestMyID;
639 }
640 
641 
642 // Create a list of start (vertex) times from various sources, ordered by
643 // uncertainty.
646  vector<const DSCHit*>&sc_hits,
647  vector<const DTOFPoint*>&tof_points,
648  vector<const DBCALShower*>&bcal_showers,
649  vector<const DFCALShower*>&fcal_showers,
650  vector<DTrackTimeBased::DStartTime_t>&start_times){
652 
653  // Match to the start counter and the outer detectors
654  double locStartTimeVariance = 0.0;
655  double track_t0=track->t0();
656  double locStartTime = track_t0; // initial guess from tracking
657 
658  // Get start time estimate from Start Counter
659  if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_START),sc_hits,locStartTime)){
660  start_time.t0=locStartTime;
661  // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready
662  start_time.t0_sigma=sqrt(locStartTimeVariance);
663  start_time.system=SYS_START;
664  start_times.push_back(start_time);
665  }
666  // Get start time estimate from TOF
667  locStartTime = track_t0; // initial guess from tracking
668  if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_TOF),tof_points,locStartTime)){
669  // Fill in the start time vector
670  start_time.t0=locStartTime;
671  start_time.t0_sigma=sqrt(locStartTimeVariance);
672  // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready
673  start_time.system=SYS_TOF;
674  start_times.push_back(start_time);
675  }
676  // Get start time estimate from FCAL
677  locStartTime = track_t0; // initial guess from tracking
678  if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_FCAL),fcal_showers,locStartTime)){
679  // Fill in the start time vector
680  start_time.t0=locStartTime;
681  start_time.t0_sigma=sqrt(locStartTimeVariance);
682  // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready
683  start_time.system=SYS_FCAL;
684  start_times.push_back(start_time);
685  }
686  // Get start time estimate from BCAL
687  locStartTime=track_t0;
688  if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_BCAL),bcal_showers,locStartTime)){
689  // Fill in the start time vector
690  start_time.t0=locStartTime;
691  start_time.t0_sigma=0.5;
692  // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready
693  start_time.system=SYS_BCAL;
694  start_times.push_back(start_time);
695  }
696  // Add the t0 estimate from the tracking
697  start_time.t0=track_t0;
698  start_time.t0_sigma=5.;
699  start_time.system=track->t0_detector();
700  start_times.push_back(start_time);
701 
702  // Set t0 for the fit to the first entry in the list. Usually this will be
703  // from the start counter.
704  mStartTime=start_times[0].t0;
705  mStartDetector=start_times[0].system;
706 
707  // _DBG_ << mStartDetector << " " << mStartTime << endl;
708 
709 }
710 
711 // Create a list of start times and do the fit for a particular mass hypothesis
713  vector<DTrackTimeBased::DStartTime_t>&start_times,
714  JEventLoop *loop,
715  double mass){
716  if(DEBUG_LEVEL>1){_DBG__;_DBG_<<"---- Starting time based fit with mass: "<<mass<<endl;}
717  // Get the hits from the wire-based track
718  vector<const DFDCPseudo*>myfdchits;
719  track->GetT(myfdchits);
720  vector<const DCDCTrackHit *>mycdchits;
721  track->GetT(mycdchits);
722 
723  // Do the fit
725  if (USE_HITS_FROM_WIREBASED_FIT) {
726  fitter->Reset();
728 
729  fitter->AddHits(myfdchits);
730  fitter->AddHits(mycdchits);
731 
732  status=fitter->FitTrack(track->position(),track->momentum(),
733  track->charge(),mass,mStartTime,mStartDetector);
734  }
735  else{
736  fitter->Reset();
738  status = fitter->FindHitsAndFitTrack(*track, track->extrapolations,loop,
739  mass,
740  mycdchits.size()+2*myfdchits.size(),
741  mStartTime,mStartDetector);
742 
743  // If the status is kFitNotDone, then not enough hits were attached to this
744  // track using the hit-gathering algorithm. In this case get the hits
745  // from the wire-based track
746  if (status==DTrackFitter::kFitNotDone){
747  //_DBG_ << " Using wire-based hits " << endl;
748  fitter->Reset();
750  fitter->AddHits(myfdchits);
751  fitter->AddHits(mycdchits);
752 
753  status=fitter->FitTrack(track->position(),track->momentum(),
754  track->charge(),mass,mStartTime,mStartDetector);
755  }
756 
757  }
758 
759  // if the fit returns chisq=-1, something went terribly wrong. We may still
760  // have a usable track from the wire-based pass. In this case set
761  // kFitNoImprovement so we can save the wire-based results.
762  if (fitter->GetChisq()<0){
764  }
765 
766  // In the transition region between the CDC and the FDC where the track
767  // contains both CDC and FDC hits, sometimes too many hits are discarded in
768  // the time-based phase and the time-based fit result does not improve on the
769  // wire-based fit result. In this case set the status word to
770  // kFitNoImprovement and copy the wire-based parameters into the time-based
771  // class.
772  if (myfdchits.size()>3 && mycdchits.size()>3){
773  unsigned int ndof=fitter->GetNdof();
774  if (TMath::Prob(track->chisq,track->Ndof)>
775  TMath::Prob(fitter->GetChisq(),ndof)&&ndof<5)
777  }
778 
779  // Check the status value from the fit
780  switch(status){
782  //_DBG_<<"Fitter returned kFitNotDone. This should never happen!!"<<endl;
784  break;
786  {
787  // Create a new time-based track object
788  DTrackTimeBased *timebased_track = new DTrackTimeBased();
789  *static_cast<DTrackingData*>(timebased_track) = *static_cast<const DTrackingData*>(track);
790 
791  timebased_track->chisq = track->chisq;
792  timebased_track->Ndof = track->Ndof;
793  timebased_track->pulls = track->pulls;
794  timebased_track->extrapolations = track->extrapolations;
795  timebased_track->trackid = track->id;
796  timebased_track->candidateid=track->candidateid;
797  timebased_track->FOM=track->FOM;
799 
800  // add the list of start times
801  timebased_track->start_times.assign(start_times.begin(),
802  start_times.end());
803 
804  for(unsigned int m=0; m<myfdchits.size(); m++)
805  timebased_track->AddAssociatedObject(myfdchits[m]);
806  for(unsigned int m=0; m<mycdchits.size(); m++)
807  timebased_track->AddAssociatedObject(mycdchits[m]);
808 
809  timebased_track->measured_cdc_hits_on_track = mycdchits.size();
810  timebased_track->measured_fdc_hits_on_track = myfdchits.size();
811 
812  // dEdx
813  double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp;
814  double locdx_CDC_amp,locdx_CDC;
815  unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC;
816  pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp, locdx_CDC, locdx_CDC_amp,locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC);
817 
818  timebased_track->ddEdx_FDC = locdEdx_FDC;
819  timebased_track->ddx_FDC = locdx_FDC;
820  timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC;
821  timebased_track->ddEdx_CDC = locdEdx_CDC;
822  timebased_track->ddEdx_CDC_amp = locdEdx_CDC_amp;
823  timebased_track->ddx_CDC = locdx_CDC;
824  timebased_track->ddx_CDC_amp = locdx_CDC_amp;
825  timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC;
826 
829 
830  timebased_track->AddAssociatedObject(track);
831  _data.push_back(timebased_track);
832 
833  return true;
834  break;
835  }
837  {
838  // Create a new time-based track object
839  DTrackTimeBased *timebased_track = new DTrackTimeBased();
840  *static_cast<DTrackingData*>(timebased_track) = fitter->GetFitParameters();
841 
842  timebased_track->setTime(mStartTime);
843  timebased_track->chisq = fitter->GetChisq();
844  timebased_track->Ndof = fitter->GetNdof();
845  timebased_track->pulls = std::move(fitter->GetPulls());
846  timebased_track->extrapolations=std::move(fitter->GetExtrapolations());
847  timebased_track->IsSmoothed = fitter->GetIsSmoothed();
848  timebased_track->trackid = track->id;
849  timebased_track->candidateid=track->candidateid;
850  timebased_track->flags=DTrackTimeBased::FLAG__GOODFIT;
851 
852  // Set the start time and add the list of start times
853  timebased_track->setT0(mStartTime,start_times[0].t0_sigma, mStartDetector);
854  timebased_track->start_times.assign(start_times.begin(), start_times.end());
855 
856  if (DEBUG_HISTS){
857  int id=0;
858  if (mStartDetector==SYS_CDC) id=1;
859  else if (mStartDetector==SYS_FDC) id=2;
860  else if (mStartDetector==SYS_BCAL) id=3;
861  else if (mStartDetector==SYS_FCAL) id=4;
862  else if (mStartDetector==SYS_TOF) id=5;
863 
864  Hstart_time->Fill(start_times[0].t0,id);
865  }
866 
867 
868  // Add hits used as associated objects
869  const vector<const DCDCTrackHit*> &cdchits = fitter->GetCDCFitHits();
870  const vector<const DFDCPseudo*> &fdchits = fitter->GetFDCFitHits();
871 
872  unsigned int num_fdc_potential=fitter->GetNumPotentialFDCHits();
873  unsigned int num_cdc_potential=fitter->GetNumPotentialCDCHits();
874 
876  temp.inner_layer=0;
877  temp.outer_layer=0;
878  temp.total_hits=num_cdc_potential;
879  if (cdchits.size()>0){
880  temp.inner_layer=cdchits[0]->wire->ring;
881  temp.outer_layer=cdchits[cdchits.size()-1]->wire->ring;
882  }
883  timebased_track->cdc_hit_usage=temp;
884 
885  // Reset the structure
886  temp.inner_layer=0;
887  temp.outer_layer=0;
888  temp.total_hits=num_fdc_potential;
889  if (fdchits.size()>0){
890  temp.inner_layer=fdchits[0]->wire->layer;
891  temp.outer_layer=fdchits[fdchits.size()-1]->wire->layer;
892  }
893  timebased_track->fdc_hit_usage=temp;
894 
895  for(unsigned int m=0; m<cdchits.size(); m++)
896  timebased_track->AddAssociatedObject(cdchits[m]);
897  for(unsigned int m=0; m<fdchits.size(); m++)
898  timebased_track->AddAssociatedObject(fdchits[m]);
899 
900  timebased_track->measured_cdc_hits_on_track = cdchits.size();
901  timebased_track->measured_fdc_hits_on_track = fdchits.size();
902 
903  // dEdx
904  double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp;
905  double locdx_CDC,locdx_CDC_amp;
906  unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC;
907  pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp,locdx_CDC,locdx_CDC_amp, locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC);
908 
909  timebased_track->ddEdx_FDC = locdEdx_FDC;
910  timebased_track->ddx_FDC = locdx_FDC;
911  timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC;
912  timebased_track->ddEdx_CDC = locdEdx_CDC;
913  timebased_track->ddEdx_CDC_amp= locdEdx_CDC_amp;
914  timebased_track->ddx_CDC = locdx_CDC;
915  timebased_track->ddx_CDC_amp= locdx_CDC_amp;
916  timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC;
917 
918  // Set CDC ring & FDC plane hit patterns before candidate and wirebased tracks are associated
919  vector<const DCDCTrackHit*> tempCDCTrackHits;
920  vector<const DFDCPseudo*> tempFDCPseudos;
921  timebased_track->Get(tempCDCTrackHits);
922  timebased_track->Get(tempFDCPseudos);
923  timebased_track->dCDCRings = pid_algorithm->Get_CDCRingBitPattern(tempCDCTrackHits);
924  timebased_track->dFDCPlanes = pid_algorithm->Get_FDCPlaneBitPattern(tempFDCPseudos);
925 
928 
929  // Add DTrack object as associate object
930  timebased_track->AddAssociatedObject(track);
931 
932  // Compute the figure-of-merit based on tracking
933  timebased_track->FOM = TMath::Prob(timebased_track->chisq, timebased_track->Ndof);
934  //_DBG_<< "FOM: " << timebased_track->FOM << endl;
935 
936  _data.push_back(timebased_track);
937 
938  return true;
939  break;
940 
941  }
942  default:
943  break;
944  }
945  return false;
946 }
947 
948 
949 // Create a track with a mass hypothesis that was not present in the list of
950 // fitted tracks from an existing fitted track.
951 void DTrackTimeBased_factory::AddMissingTrackHypothesis(vector<DTrackTimeBased*>&tracks_to_add,
952  const DTrackTimeBased *src_track,
953  double my_mass,
954  double q,
955  JEventLoop *loop){
956 
957 
958  // Create a new time-based track object
959  DTrackTimeBased *timebased_track = new DTrackTimeBased();
960  *static_cast<DTrackingData*>(timebased_track) = *static_cast<const DTrackingData*>(src_track);
961 
962  // Get the hits used in the fit
963  vector<const DCDCTrackHit *>src_cdchits;
964  src_track->GetT(src_cdchits);
965  vector<const DFDCPseudo *>src_fdchits;
966  src_track->GetT(src_fdchits);
967 
968  // Copy over DKinematicData part from the result of a successful fit
969  timebased_track->setPID(IDTrack(q, my_mass));
970  timebased_track->chisq = src_track->chisq;
971  timebased_track->Ndof = src_track->Ndof;
972  timebased_track->pulls = src_track->pulls;
973  timebased_track->extrapolations = src_track->extrapolations;
974  timebased_track->trackid = src_track->id;
975  timebased_track->candidateid=src_track->candidateid;
976  timebased_track->FOM=src_track->FOM;
977  timebased_track->cdc_hit_usage=src_track->cdc_hit_usage;
978  timebased_track->fdc_hit_usage=src_track->fdc_hit_usage;
980 
981  // Add list of start times
982  timebased_track->start_times.assign(src_track->start_times.begin(),
983  src_track->start_times.end());
984  // Set the start time we used
985  timebased_track->setT0(timebased_track->start_times[0].t0,
986  timebased_track->start_times[0].t0_sigma,
987  timebased_track->start_times[0].system);
988 
989  // Add DTrack object as associate object
990  vector<const DTrackWireBased*>wire_based_track;
991  src_track->GetT(wire_based_track);
992  timebased_track->AddAssociatedObject(wire_based_track[0]);
993 
994  // (Partially) compensate for the difference in energy loss between the
995  // source track and a particle of mass my_mass
996  DVector3 position,momentum;
997  if (timebased_track->extrapolations.at(SYS_CDC).size()>0){
998  unsigned int index=timebased_track->extrapolations.at(SYS_CDC).size()-1;
999  position=timebased_track->extrapolations[SYS_CDC][index].position;
1000  momentum=timebased_track->extrapolations[SYS_CDC][index].momentum;
1001  }
1002  else if (timebased_track->extrapolations.at(SYS_FDC).size()>0){
1003  unsigned int index=timebased_track->extrapolations.at(SYS_FDC).size()-1;
1004  position=timebased_track->extrapolations[SYS_FDC][index].position;
1005  momentum=timebased_track->extrapolations[SYS_FDC][index].momentum;
1006  }
1007 
1009  if (momentum.Mag()>0.){
1010  CorrectForELoss(position,momentum,q,my_mass);
1011  timebased_track->setMomentum(momentum);
1012  timebased_track->setPosition(position);
1013 
1014  // Redo the fit with the new position and momentum as initial guesses
1015  fitter->Reset();
1017  status = fitter->FindHitsAndFitTrack(*timebased_track,
1018  timebased_track->extrapolations,loop,
1019  my_mass,
1020  src_cdchits.size()+2*src_fdchits.size(),
1021  timebased_track->t0(),
1022  timebased_track->t0_detector());
1023  // if the fit returns chisq=-1, something went terribly wrong. Do not
1024  // update the parameters for the track...
1025  if (fitter->GetChisq()<0) status=DTrackFitter::kFitFailed;
1026 
1027  // if the fit flips the charge of the track, then this is bad as well
1028  if(q != fitter->GetFitParameters().charge())
1029  status=DTrackFitter::kFitFailed;
1030 
1031  // if we can't refit the track, it is likely of poor quality, so stop here and do not add the hypothesis
1032  if(status == DTrackFitter::kFitFailed) {
1033  delete timebased_track;
1034  return;
1035  }
1036 
1037  if (status==DTrackFitter::kFitSuccess){
1038  timebased_track->chisq = fitter->GetChisq();
1039  timebased_track->Ndof = fitter->GetNdof();
1040  timebased_track->pulls = std::move(fitter->GetPulls());
1041  timebased_track->extrapolations=std::move(fitter->GetExtrapolations());
1042  timebased_track->IsSmoothed = fitter->GetIsSmoothed();
1043  *static_cast<DTrackingData*>(timebased_track) = fitter->GetFitParameters();
1044  timebased_track->flags=DTrackTimeBased::FLAG__GOODFIT;
1045 
1046  // Add hits used as associated objects
1047  const vector<const DCDCTrackHit*> &cdchits = fitter->GetCDCFitHits();
1048  const vector<const DFDCPseudo*> &fdchits = fitter->GetFDCFitHits();
1049 
1050  unsigned int num_fdc_potential=fitter->GetNumPotentialFDCHits();
1051  unsigned int num_cdc_potential=fitter->GetNumPotentialCDCHits();
1052 
1054  temp.inner_layer=0;
1055  temp.outer_layer=0;
1056  temp.total_hits=num_cdc_potential;
1057  if (cdchits.size()>0){
1058  temp.inner_layer=cdchits[0]->wire->ring;
1059  temp.outer_layer=cdchits[cdchits.size()-1]->wire->ring;
1060  }
1061  timebased_track->cdc_hit_usage=temp;
1062 
1063  // Reset the structure
1064  temp.inner_layer=0;
1065  temp.outer_layer=0;
1066  temp.total_hits=num_fdc_potential;
1067  if (fdchits.size()>0){
1068  temp.inner_layer=fdchits[0]->wire->layer;
1069  temp.outer_layer=fdchits[fdchits.size()-1]->wire->layer;
1070  }
1071  timebased_track->fdc_hit_usage=temp;
1072 
1073  for(unsigned int m=0; m<cdchits.size(); m++)
1074  timebased_track->AddAssociatedObject(cdchits[m]);
1075  for(unsigned int m=0; m<fdchits.size(); m++)
1076  timebased_track->AddAssociatedObject(fdchits[m]);
1077 
1078  timebased_track->measured_cdc_hits_on_track = cdchits.size();
1079  timebased_track->measured_fdc_hits_on_track = fdchits.size();
1080 
1081  // Compute the figure-of-merit based on tracking
1082  timebased_track->FOM = TMath::Prob(timebased_track->chisq, timebased_track->Ndof);
1083 
1084  }
1085  }
1086 
1087  if (status!=DTrackFitter::kFitSuccess){
1088  for(unsigned int m=0; m<src_fdchits.size(); m++)
1089  timebased_track->AddAssociatedObject(src_fdchits[m]);
1090  for(unsigned int m=0; m<src_cdchits.size(); m++)
1091  timebased_track->AddAssociatedObject(src_cdchits[m]);
1092 
1093  timebased_track->measured_cdc_hits_on_track = src_cdchits.size();
1094  timebased_track->measured_fdc_hits_on_track = src_fdchits.size();
1095  }
1096 
1097  // dEdx
1098  double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp;
1099  double locdx_CDC,locdx_CDC_amp;
1100  unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC;
1101  pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp,locdx_CDC,locdx_CDC_amp, locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC);
1102 
1103  timebased_track->ddEdx_FDC = locdEdx_FDC;
1104  timebased_track->ddx_FDC = locdx_FDC;
1105  timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC;
1106  timebased_track->ddEdx_CDC = locdEdx_CDC;
1107  timebased_track->ddEdx_CDC_amp = locdEdx_CDC_amp;
1108  timebased_track->ddx_CDC = locdx_CDC;
1109  timebased_track->ddx_CDC_amp = locdx_CDC_amp;
1110  timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC;
1111 
1112  // Set CDC ring & FDC plane hit patterns before candidate and wirebased tracks are associated
1113  vector<const DCDCTrackHit*> tempCDCTrackHits;
1114  vector<const DFDCPseudo*> tempFDCPseudos;
1115  timebased_track->Get(tempCDCTrackHits);
1116  timebased_track->Get(tempFDCPseudos);
1117  timebased_track->dCDCRings = pid_algorithm->Get_CDCRingBitPattern(tempCDCTrackHits);
1118  timebased_track->dFDCPlanes = pid_algorithm->Get_FDCPlaneBitPattern(tempFDCPseudos);
1119 
1122 
1123  tracks_to_add.push_back(timebased_track);
1124 }
1125 
1126 
1127 // If the fit failed for certain hypotheses, fill in the gaps using data from
1128 // successful fits for each candidate.
1130  if (_data.size()==0) return false;
1131 
1132  // Make sure the tracks are ordered by candidate id
1133  sort(_data.begin(),_data.end(),DTrackTimeBased_cmp);
1134 
1135  JObject::oid_t old_id=_data[0]->candidateid;
1136  unsigned int mass_bits=0;
1137  double q=_data[0]->charge();
1138  bool flipped_charge=false;
1139  vector<DTrackTimeBased*>myhypotheses;
1140  vector<DTrackTimeBased*>tracks_to_add;
1141  for (size_t i=0;i<_data.size();i++){
1142  if (_data[i]->candidateid!=old_id){
1143  int num_hyp=myhypotheses.size();
1144  if ((q<0 && num_hyp!=mNumHypMinus)||(q>0 && num_hyp!=mNumHypPlus)
1145  || flipped_charge){
1146  AddMissingTrackHypotheses(mass_bits,tracks_to_add,myhypotheses,q,
1147  flipped_charge,loop);
1148  }
1149 
1150  // Clear the myhypotheses vector for the next track
1151  myhypotheses.clear();
1152  // Reset flags and charge
1153  q=_data[i]->charge();
1154  flipped_charge=false;
1155  // Set the bit for this mass hypothesis
1156  mass_bits = 1<<_data[i]->PID();
1157 
1158  // Add the data to the myhypotheses vector
1159  myhypotheses.push_back(_data[i]);
1160  }
1161  else{
1162  myhypotheses.push_back(_data[i]);
1163 
1164  // Set the bit for this mass hypothesis
1165  mass_bits |= 1<< _data[i]->PID();
1166 
1167  // Check if the sign of the charge has flipped
1168  if (_data[i]->charge()!=q) flipped_charge=true;
1169  }
1170 
1171  old_id=_data[i]->candidateid;
1172  }
1173  // Deal with last track candidate
1174  int num_hyp=myhypotheses.size();
1175  if ((q<0 && num_hyp!=mNumHypMinus)||(q>0 && num_hyp!=mNumHypPlus)
1176  || flipped_charge){
1177  AddMissingTrackHypotheses(mass_bits,tracks_to_add,myhypotheses,q,
1178  flipped_charge,loop);
1179  }
1180 
1181  // Add the new list of tracks to the output list
1182  if (tracks_to_add.size()>0){
1183  //_DBG_ << "Adding tracks " <<endl;
1184  for (size_t i=0;i<tracks_to_add.size();i++){
1185  _data.push_back(tracks_to_add[i]);
1186  }
1187  // Make sure the tracks are ordered by candidate id
1188  sort(_data.begin(),_data.end(),DTrackTimeBased_cmp);
1189  }
1190 
1191  return true;
1192 }
1193 
1194 // Use the FastSwim method in DReferenceTrajectory to propagate back to the
1195 // POCA to the beam line, adding a bit of energy at each step that would have
1196 // been lost had the particle emerged from the target.
1197 void DTrackTimeBased_factory::CorrectForELoss(DVector3 &position,DVector3 &momentum,double q,double my_mass){
1199  DReferenceTrajectory rt(fitter->GetDMagneticFieldMap(),q,&dummy_step);
1200  rt.SetDGeometry(geom);
1201  rt.SetMass(my_mass);
1202  rt.SetPLossDirection(DReferenceTrajectory::kBackward);
1203  DVector3 last_pos,last_mom;
1204  DVector3 origin(0.,0.,65.);
1205  DVector3 dir(0.,0.,1.);
1206  rt.FastSwim(position,momentum,last_pos,last_mom,q,origin,dir,300.);
1207  position=last_pos;
1208  momentum=last_mom;
1209 }
1210 
1211 // Fill in all missing hypotheses for a given track candidate
1213  vector<DTrackTimeBased*>&tracks_to_add,
1214  vector<DTrackTimeBased *>&myhypotheses,
1215  double q,
1216  bool flipped_charge,
1217  JEventLoop *loop){
1218  Particle_t negative_particles[3]={KMinus,PiMinus,Electron};
1219  Particle_t positive_particles[3]={KPlus,PiPlus,Positron};
1220 
1221  unsigned int last_index=myhypotheses.size()-1;
1222  if (q>0){
1223  if (flipped_charge){
1224  if ((mass_bits & (1<<AntiProton))==0){
1225  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1226  ParticleMass(Proton),-1.,loop);
1227  }
1228  for (int i=0;i<3;i++){
1229  if ((mass_bits & (1<<negative_particles[i]))==0){
1230  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1231  ParticleMass(negative_particles[i]),
1232  -1.,loop);
1233  }
1234  }
1235  }
1236  if ((mass_bits & (1<<Proton))==0){
1237  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1238  ParticleMass(Proton),+1.,loop);
1239  }
1240  for (int i=0;i<3;i++){
1241  if ((mass_bits & (1<<positive_particles[i]))==0){
1242  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1243  ParticleMass(positive_particles[i]),
1244  +1.,loop);
1245  }
1246  }
1247  }
1248  else{
1249  if ((mass_bits & (1<<AntiProton))==0){
1250  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1251  ParticleMass(Proton),-1.,loop);
1252  }
1253  for (int i=0;i<3;i++){
1254  if ((mass_bits & (1<<negative_particles[i]))==0){
1255  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1256  ParticleMass(negative_particles[i]),
1257  -1.,loop);
1258  }
1259  }
1260  if (flipped_charge){
1261  if ((mass_bits & (1<<Proton))==0){
1262  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1263  ParticleMass(Proton),+1.,loop);
1264  }
1265  for (int i=0;i<3;i++){
1266  if ((mass_bits & (1<<positive_particles[i]))==0){
1267  AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1268  ParticleMass(positive_particles[i]),
1269  +1.,loop);
1270  }
1271  }
1272  }
1273  }
1274 }
map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > extrapolations
int GetNdof(void) const
Definition: DTrackFitter.h:155
Definition: track.h:16
void setMomentum(const DVector3 &aMomentum)
void setTime(double locTime)
DApplication * dapp
vector< pull_t > & GetPulls(void)
Definition: DTrackFitter.h:162
float chisq
Chi-squared for the track (not chisq/dof!)
float chisq
Chi-squared for the track (not chisq/dof!)
double t0(void) const
Definition: DTrackingData.h:23
The DTrackFitter class is a base class for different charged track fitting algorithms. It does not actually fit the track itself, but provides the interface and some common support features most algorthims will need to implement.
Definition: DTrackFitter.h:61
static unsigned int count_common_members(vector< T > &a, vector< T > &b)
TVector3 DVector3
Definition: DVector3.h:14
double total
Definition: FitGains.C:151
bool DoFit(const DTrackWireBased *track, vector< DTrackTimeBased::DStartTime_t > &start_times, JEventLoop *loop, double mass)
oid_t trackid
id of DTrackWireBased corresponding to this track
void AddHits(vector< const DCDCTrackHit * > cdchits)
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
jerror_t init(void)
Called once at program start.
bool DTrackTimeBased_cmp(DTrackTimeBased *a, DTrackTimeBased *b)
const DVector3 & position(void) const
bool GetIsSmoothed(void) const
Definition: DTrackFitter.h:160
unsigned int potential_cdc_hits_on_track
jerror_t brun(jana::JEventLoop *loop, int32_t runnumber)
Called everytime a new run number is detected.
void setT0(double at0, double at0_err, DetectorSystem_t at0_detector)
Definition: DTrackingData.h:64
jerror_t fini(void)
Called after last event of last event source has been processed.
TLorentzVector DLorentzVector
static int ParticleCharge(Particle_t p)
static char index(char c)
Definition: base64.cpp:115
oid_t candidateid
id of DTrackCandidate corresponding to this track
void SetDGeometry(const DGeometry *geom)
Definition: GlueX.h:17
Definition: GlueX.h:19
TF1 * f
Definition: FitGains.C:21
bool DTrackTimeBased_T0_cmp(DTrackTimeBased::DStartTime_t a, DTrackTimeBased::DStartTime_t b)
unsigned int dFDCPlanes
fit_status_t FindHitsAndFitTrack(const DKinematicData &starting_params, const DReferenceTrajectory *rt, JEventLoop *loop, double mass=-1.0, int N=0, double t0=QuietNaN, DetectorSystem_t t0_det=SYS_NULL)
mass&lt;0 means get it from starting_params
hit_usage_t cdc_hit_usage
void Reset(void)
Definition: DTrackFitter.cc:94
void AddMissingTrackHypotheses(unsigned int mass_bits, vector< DTrackTimeBased * > &tracks_to_add, vector< DTrackTimeBased * > &hypotheses, double q, bool flipped_charge, JEventLoop *loop)
size_t Get_NumTrackHits(const DTrackTimeBased *locTrackTimeBased)
vector< DTrackFitter::pull_t > pulls
Holds pulls used in chisq calc. (not including off-diagonals)
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
const map< DetectorSystem_t, vector< Extrapolation_t > > & GetExtrapolations(void) const
Definition: DTrackFitter.h:163
double GetChisq(void) const
Definition: DTrackFitter.h:154
Definition: GlueX.h:20
Definition: GlueX.h:18
void SetFitType(fit_type_t type)
Definition: DTrackFitter.h:170
void SplitString(string str, vector< T > &vec, const string &delim=" ")
Definition: SplitString.h:7
DGeometry * GetDGeometry(unsigned int run_number)
Definition: GlueX.h:22
#define _DBG_
Definition: HDEVIO.h:12
unsigned int potential_fdc_hits_on_track
void CreateStartTimeList(const DTrackWireBased *track, vector< const DSCHit * > &sc_hits, vector< const DTOFPoint * > &tof_points, vector< const DBCALShower * > &bcal_showers, vector< const DFCALShower * > &fcal_showers, vector< DTrackTimeBased::DStartTime_t > &start_times)
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
jerror_t evnt(jana::JEventLoop *loop, uint64_t eventnumber)
Called every event.
static double start_time
Definition: mc2coda.c:26
#define _DBG__
Definition: HDEVIO.h:13
const vector< const DCDCTrackHit * > & GetCDCFitHits(void) const
Definition: DTrackFitter.h:139
unsigned int GetNumPotentialCDCHits(void) const
Definition: DTrackFitter.h:157
int Ndof
Number of degrees of freedom in the fit.
void setPID(Particle_t locPID)
DLorentzVector lorentzMomentum(void) const
hit_usage_t fdc_hit_usage
const vector< const DFDCPseudo * > & GetFDCFitHits(void) const
Definition: DTrackFitter.h:140
static double ParticleMass(Particle_t p)
double sqrt(double)
unsigned int dCDCRings
vector< DTrackFitter::pull_t > pulls
Holds pulls used in chisq calc. (not including off-diagonals)
bool InsertMissingHypotheses(JEventLoop *loop)
unsigned int measured_cdc_hits_on_track
const DVector3 & momentum(void) const
fit_status_t FitTrack(const DVector3 &pos, const DVector3 &mom, double q, double mass, double t0=QuietNaN, DetectorSystem_t t0_det=SYS_NULL)
DetectorSystem_t t0_detector(void) const
Definition: DTrackingData.h:25
const DTrackFitter * fitter
oid_t candidateid
which DTrackCandidate this came from
double GetTruthMatchingFOM(int trackIndex, DTrackTimeBased *dtrack, vector< const DMCThrown * >mcthrowns)
map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > extrapolations
const DMagneticFieldMap * GetDMagneticFieldMap(void) const
Definition: DTrackFitter.h:168
void setPosition(const DVector3 &aPosition)
TDirectory * dir
Definition: bcal_hist_eff.C:31
int GetThrownIndex(vector< const DMCThrown * > &locMCThrowns, const DKinematicData *kd, double &f)
static Particle_t IDTrack(float locCharge, float locMass)
unsigned int dNumHitsUsedFordEdx_FDC
TH2F * temp
unsigned int measured_fdc_hits_on_track
unsigned int dNumHitsUsedFordEdx_CDC
void AddMissingTrackHypothesis(vector< DTrackTimeBased * > &tracks_to_add, const DTrackTimeBased *src_track, double my_mass, double q, JEventLoop *loop)
unsigned int GetNumPotentialFDCHits(void) const
Definition: DTrackFitter.h:156
vector< DStartTime_t > start_times
double mass(void) const
void CorrectForELoss(DVector3 &position, DVector3 &momentum, double q, double my_mass)
Particle_t
Definition: particleType.h:12
const DTrackingData & GetFitParameters(void) const
Definition: DTrackFitter.h:153