16 #define TOF_SIGMA 0.080 // TOF resolution in ns
50 for(
unsigned int i=0; i<a.size(); i++){
51 for(
unsigned int j=0; j<b.size(); j++){
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);
79 gPARMS->SetDefaultParameter(
"TRKFIT:DEBUG_HISTS",DEBUG_HISTS);
80 gPARMS->SetDefaultParameter(
"TRKFIT:DEBUG_LEVEL",DEBUG_LEVEL);
82 vector<int> hypotheses;
84 hypotheses.push_back(
PiPlus);
85 hypotheses.push_back(
KPlus);
86 hypotheses.push_back(
Proton);
89 hypotheses.push_back(
KMinus);
92 ostringstream locMassStream;
93 for(
size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i)
95 locMassStream << hypotheses[loc_i];
96 if(loc_i != (hypotheses.size() - 1))
100 string HYPOTHESES = locMassStream.str();
101 gPARMS->SetDefaultParameter(
"TRKFIT:HYPOTHESES", HYPOTHESES);
106 for(
size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i)
109 mass_hypotheses_positive.push_back(hypotheses[loc_i]);
111 mass_hypotheses_negative.push_back(hypotheses[loc_i]);
114 if(mass_hypotheses_positive.empty()){
115 static once_flag pwarn_flag;
116 call_once(pwarn_flag, [](){
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;
125 if(mass_hypotheses_negative.empty()){
126 static once_flag nwarn_flag;
127 call_once(nwarn_flag, [](){
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;
137 mNumHypPlus=mass_hypotheses_positive.size();
138 mNumHypMinus=mass_hypotheses_negative.size();
141 PID_FORCE_TRUTH =
false;
142 gPARMS->SetDefaultParameter(
"TRKFIT:PID_FORCE_TRUTH", PID_FORCE_TRUTH);
145 gPARMS->SetDefaultParameter(
"TRKFIT:USE_SC_TIME",USE_SC_TIME);
148 gPARMS->SetDefaultParameter(
"TRKFIT:USE_TOF_TIME",USE_TOF_TIME);
151 gPARMS->SetDefaultParameter(
"TRKFIT:USE_FCAL_TIME",USE_FCAL_TIME);
154 gPARMS->SetDefaultParameter(
"TRKFIT:USE_BCAL_TIME",USE_BCAL_TIME);
175 SetFactoryFlag(NOT_OBJECT_OWNER);
178 ClearFactoryFlag(NOT_OBJECT_OWNER);
182 vector<const DTrackFitter *> 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;
194 _DBG_<<
"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
195 return RESOURCE_UNAVAILABLE;
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;
206 pid_algorithm = pid_algorithms[0];
210 _DBG_<<
"Unable to get a DParticleID object! NO PID will be done!"<<endl;
211 return RESOURCE_UNAVAILABLE;
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");
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.);
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);
250 if(!
fitter)
return NOERROR;
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]));
266 vector<const DTrackWireBased*> tracks;
268 if (tracks.size()==0)
return NOERROR;
271 vector<const DSCHit*>sc_hits;
277 vector<const DTOFPoint*> tof_points;
279 loop->Get(tof_points);
283 vector<const DBCALShower*>bcal_showers;
285 loop->Get(bcal_showers);
287 vector<const DFCALShower*>fcal_showers;
289 loop->Get(fcal_showers);
292 vector<const DMCThrown*> mcthrowns;
293 loop->Get(mcthrowns,
"FinalState");
296 for(
unsigned int i=0; i<tracks.size(); i++){
299 unsigned int num=_data.size();
302 vector<DTrackTimeBased::DStartTime_t>start_times;
303 CreateStartTimeList(track,sc_hits,tof_points,bcal_showers,fcal_showers,start_times);
306 DoFit(track,start_times,loop,track->
mass());
309 if (PID_FORCE_TRUTH && _data.size()>num) {
313 _data[_data.size()-1]->FOM=GetTruthMatchingFOM(i,_data[_data.size()-1],
322 if (INSERT_MISSING_HYPOTHESES){
323 InsertMissingHypotheses(loop);
327 for(
size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
329 if(!mcthrowns.empty())
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);
338 _data[loc_i]->dMCThrownMatchMyID = -1;
339 _data[loc_i]->dNumHitsMatchedToThrown = 0;
370 if(_data.size()==0)
return;
372 if(DEBUG_LEVEL>2)
_DBG_<<
"Looking for clones of time-based tracks ..."<<endl;
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++){
390 vector<const DCDCTrackHit*> cdchits1;
391 vector<const DFDCPseudo*> fdchits1;
392 dtrack1->Get(cdchits1);
393 dtrack1->Get(fdchits1);
395 unsigned int num_cdc1=cdchits1.size();
396 unsigned int num_fdc1=fdchits1.size();
397 unsigned int total1 = num_cdc1+num_fdc1;
400 for(
unsigned int j=i+1; j<_data.size(); j++){
404 vector<const DCDCTrackHit*> cdchits2;
405 vector<const DFDCPseudo*> fdchits2;
406 dtrack2->Get(cdchits2);
407 dtrack2->Get(fdchits2);
410 unsigned int num_cdc2=cdchits2.size();
411 unsigned int num_fdc2=fdchits2.size();
412 unsigned int total2 = num_cdc2+num_fdc2;
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;
423 unsigned int total = Ncdc + Nfdc;
426 if (total<=1)
continue;
431 if (Ncdc>0 && (num_fdc1>0 || num_fdc2>0)
432 && (num_fdc1*num_fdc2)==0)
continue;
437 if (Nfdc>0 && (num_cdc1>0 || num_cdc2>0)
438 && (num_cdc1*num_cdc2)==0)
continue;
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;
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;
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);
461 candidates_to_delete.push_back(cand1);
462 candidates_to_keep.push_back(dtrack2->
candidateid);
467 if(DEBUG_LEVEL>2)
_DBG_<<
"Found "<<candidates_to_delete.size()<<
" time-based clones"<<endl;
471 if(candidates_to_delete.size()==0)
return;
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);
486 vector<DTrackTimeBased*> new_data;
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]){
494 _DBG_<<
"Deleting clone time-based fitted result "<<i
495 <<
" in event " << myevt << endl;
501 new_data.push_back(_data[i]);
503 else delete _data[i];
516 vector<DLorentzVector> gen_fourMom(mcthrowns.size());
517 for(
unsigned int i=0; i<mcthrowns.size(); i++){
518 gen_fourMom[i] = mcthrowns[i]->lorentzMomentum();
523 int thrownIndex = GetThrownIndex(mcthrowns, track, f);
524 if(thrownIndex<=0 || f<=0.5)
return 0.;
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;
528 double delta_phi = (fourMom.Phi()-gen_fourMom[thrownIndex-1].Phi())*1000.0;
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);
531 if (fabs(track->
mass()-mcthrowns[thrownIndex-1]->mass())<0.01 && track->
charge()==mcthrowns[thrownIndex-1]->charge())
534 double trk_chi2=track->
chisq;
535 unsigned int ndof=track->
Ndof;
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));
558 return (match) ? TMath::Prob(chisq,3) : 0.0;
568 vector<const DCDCTrackHit*> cdctrackhits;
569 kd->Get(cdctrackhits);
570 vector<const DFDCPseudo*> fdcpseudos;
573 int locTotalNumHits = cdctrackhits.size() + fdcpseudos.size();
574 if(locTotalNumHits == 0)
584 map<int, int> locHitMatches;
585 for(
size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
587 if(fabs(locMCThrowns[loc_i]->charge()) > 0.9)
588 locHitMatches[locMCThrowns[loc_i]->myid] = 0;
592 for(
size_t loc_i = 0; loc_i < cdctrackhits.size(); ++loc_i)
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;
600 int itrack = locTruthCDCHits[0]->itrack;
601 if(locHitMatches.find(itrack) == locHitMatches.end())
603 ++locHitMatches[itrack];
607 for(
size_t loc_i = 0; loc_i < fdcpseudos.size(); ++loc_i)
609 vector<const DMCTrackHit*> mctrackhits;
610 fdcpseudos[loc_i]->Get(mctrackhits);
611 if(mctrackhits.empty())
613 if(!mctrackhits[0]->primary)
616 int itrack = mctrackhits[0]->itrack;
617 if(locHitMatches.find(itrack) == locHitMatches.end())
619 ++locHitMatches[itrack];
623 map<int, int>::iterator locIterator = locHitMatches.begin();
624 int locBestMyID = -1;
625 int locBestNumHits = 0;
626 for(; locIterator != locHitMatches.end(); ++locIterator)
628 if(locIterator->second <= locBestNumHits)
630 locBestNumHits = locIterator->second;
631 locBestMyID = locIterator->first;
635 f = 1.0*locBestNumHits/locTotalNumHits;
636 if(DEBUG_HISTS)hitMatchFOM->Fill(f);
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){
654 double locStartTimeVariance = 0.0;
655 double track_t0=track->
t0();
656 double locStartTime = track_t0;
660 start_time.
t0=locStartTime;
664 start_times.push_back(start_time);
667 locStartTime = track_t0;
670 start_time.
t0=locStartTime;
674 start_times.push_back(start_time);
677 locStartTime = track_t0;
680 start_time.
t0=locStartTime;
684 start_times.push_back(start_time);
687 locStartTime=track_t0;
690 start_time.
t0=locStartTime;
694 start_times.push_back(start_time);
697 start_time.
t0=track_t0;
700 start_times.push_back(start_time);
704 mStartTime=start_times[0].t0;
705 mStartDetector=start_times[0].system;
713 vector<DTrackTimeBased::DStartTime_t>&start_times,
716 if(DEBUG_LEVEL>1){
_DBG__;
_DBG_<<
"---- Starting time based fit with mass: "<<mass<<endl;}
718 vector<const DFDCPseudo*>myfdchits;
719 track->GetT(myfdchits);
720 vector<const DCDCTrackHit *>mycdchits;
721 track->GetT(mycdchits);
725 if (USE_HITS_FROM_WIREBASED_FIT) {
733 track->
charge(),mass,mStartTime,mStartDetector);
740 mycdchits.size()+2*myfdchits.size(),
741 mStartTime,mStartDetector);
754 track->
charge(),mass,mStartTime,mStartDetector);
772 if (myfdchits.size()>3 && mycdchits.size()>3){
774 if (TMath::Prob(track->
chisq,track->
Ndof)>
789 *
static_cast<DTrackingData*
>(timebased_track) = *static_cast<const DTrackingData*>(track);
792 timebased_track->
Ndof = track->
Ndof;
795 timebased_track->
trackid = track->id;
797 timebased_track->
FOM=track->
FOM;
801 timebased_track->
start_times.assign(start_times.begin(),
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]);
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);
818 timebased_track->
ddEdx_FDC = locdEdx_FDC;
819 timebased_track->
ddx_FDC = locdx_FDC;
821 timebased_track->
ddEdx_CDC = locdEdx_CDC;
823 timebased_track->
ddx_CDC = locdx_CDC;
830 timebased_track->AddAssociatedObject(track);
831 _data.push_back(timebased_track);
842 timebased_track->
setTime(mStartTime);
848 timebased_track->
trackid = track->id;
853 timebased_track->
setT0(mStartTime,start_times[0].t0_sigma, mStartDetector);
854 timebased_track->
start_times.assign(start_times.begin(), start_times.end());
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;
864 Hstart_time->Fill(start_times[0].t0,
id);
879 if (cdchits.size()>0){
881 temp.
outer_layer=cdchits[cdchits.size()-1]->wire->ring;
889 if (fdchits.size()>0){
891 temp.
outer_layer=fdchits[fdchits.size()-1]->wire->layer;
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]);
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);
909 timebased_track->
ddEdx_FDC = locdEdx_FDC;
910 timebased_track->
ddx_FDC = locdx_FDC;
912 timebased_track->
ddEdx_CDC = locdEdx_CDC;
914 timebased_track->
ddx_CDC = locdx_CDC;
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);
930 timebased_track->AddAssociatedObject(track);
933 timebased_track->
FOM = TMath::Prob(timebased_track->
chisq, timebased_track->
Ndof);
936 _data.push_back(timebased_track);
960 *
static_cast<DTrackingData*
>(timebased_track) = *static_cast<const DTrackingData*>(src_track);
963 vector<const DCDCTrackHit *>src_cdchits;
964 src_track->GetT(src_cdchits);
965 vector<const DFDCPseudo *>src_fdchits;
966 src_track->GetT(src_fdchits);
971 timebased_track->
Ndof = src_track->
Ndof;
974 timebased_track->
trackid = src_track->id;
976 timebased_track->
FOM=src_track->
FOM;
990 vector<const DTrackWireBased*>wire_based_track;
991 src_track->GetT(wire_based_track);
992 timebased_track->AddAssociatedObject(wire_based_track[0]);
1009 if (momentum.Mag()>0.){
1010 CorrectForELoss(position,momentum,q,my_mass);
1020 src_cdchits.size()+2*src_fdchits.size(),
1021 timebased_track->
t0(),
1033 delete timebased_track;
1057 if (cdchits.size()>0){
1059 temp.
outer_layer=cdchits[cdchits.size()-1]->wire->ring;
1067 if (fdchits.size()>0){
1069 temp.
outer_layer=fdchits[fdchits.size()-1]->wire->layer;
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]);
1082 timebased_track->
FOM = TMath::Prob(timebased_track->
chisq, timebased_track->
Ndof);
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]);
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);
1103 timebased_track->
ddEdx_FDC = locdEdx_FDC;
1104 timebased_track->
ddx_FDC = locdx_FDC;
1106 timebased_track->
ddEdx_CDC = locdEdx_CDC;
1108 timebased_track->
ddx_CDC = locdx_CDC;
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);
1123 tracks_to_add.push_back(timebased_track);
1130 if (_data.size()==0)
return false;
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)
1146 AddMissingTrackHypotheses(mass_bits,tracks_to_add,myhypotheses,q,
1147 flipped_charge,loop);
1151 myhypotheses.clear();
1153 q=_data[i]->charge();
1154 flipped_charge=
false;
1156 mass_bits = 1<<_data[i]->PID();
1159 myhypotheses.push_back(_data[i]);
1162 myhypotheses.push_back(_data[i]);
1165 mass_bits |= 1<< _data[i]->PID();
1168 if (_data[i]->charge()!=q) flipped_charge=
true;
1171 old_id=_data[i]->candidateid;
1174 int num_hyp=myhypotheses.size();
1175 if ((q<0 && num_hyp!=mNumHypMinus)||(q>0 && num_hyp!=mNumHypPlus)
1177 AddMissingTrackHypotheses(mass_bits,tracks_to_add,myhypotheses,q,
1178 flipped_charge,loop);
1182 if (tracks_to_add.size()>0){
1184 for (
size_t i=0;i<tracks_to_add.size();i++){
1185 _data.push_back(tracks_to_add[i]);
1201 rt.SetMass(my_mass);
1206 rt.FastSwim(position,momentum,last_pos,last_mom,q,origin,dir,300.);
1213 vector<DTrackTimeBased*>&tracks_to_add,
1214 vector<DTrackTimeBased *>&myhypotheses,
1216 bool flipped_charge,
1221 unsigned int last_index=myhypotheses.size()-1;
1223 if (flipped_charge){
1225 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1228 for (
int i=0;i<3;i++){
1229 if ((mass_bits & (1<<negative_particles[i]))==0){
1230 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1236 if ((mass_bits & (1<<
Proton))==0){
1237 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1240 for (
int i=0;i<3;i++){
1241 if ((mass_bits & (1<<positive_particles[i]))==0){
1242 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1250 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1253 for (
int i=0;i<3;i++){
1254 if ((mass_bits & (1<<negative_particles[i]))==0){
1255 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
1260 if (flipped_charge){
1261 if ((mass_bits & (1<<
Proton))==0){
1262 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index],
1265 for (
int i=0;i<3;i++){
1266 if ((mass_bits & (1<<positive_particles[i]))==0){
1267 AddMissingTrackHypothesis(tracks_to_add,myhypotheses[0],
map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > extrapolations
void setMomentum(const DVector3 &aMomentum)
void setTime(double locTime)
vector< pull_t > & GetPulls(void)
float chisq
Chi-squared for the track (not chisq/dof!)
float chisq
Chi-squared for the track (not chisq/dof!)
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.
static unsigned int count_common_members(vector< T > &a, vector< T > &b)
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
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)
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)
oid_t candidateid
id of DTrackCandidate corresponding to this track
void SetDGeometry(const DGeometry *geom)
bool DTrackTimeBased_T0_cmp(DTrackTimeBased::DStartTime_t a, DTrackTimeBased::DStartTime_t b)
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<0 means get it from starting_params
hit_usage_t cdc_hit_usage
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
double GetChisq(void) const
void SetFitType(fit_type_t type)
void SplitString(string str, vector< T > &vec, const string &delim=" ")
DGeometry * GetDGeometry(unsigned int run_number)
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.
const vector< const DCDCTrackHit * > & GetCDCFitHits(void) const
unsigned int GetNumPotentialCDCHits(void) const
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
static double ParticleMass(Particle_t p)
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
const DTrackFitter * fitter
oid_t candidateid
which DTrackCandidate this came from
double GetTruthMatchingFOM(int trackIndex, DTrackTimeBased *dtrack, vector< const DMCThrown * >mcthrowns)
void FilterDuplicates(void)
map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > extrapolations
const DMagneticFieldMap * GetDMagneticFieldMap(void) const
void setPosition(const DVector3 &aPosition)
int GetThrownIndex(vector< const DMCThrown * > &locMCThrowns, const DKinematicData *kd, double &f)
static Particle_t IDTrack(float locCharge, float locMass)
unsigned int dNumHitsUsedFordEdx_FDC
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
vector< DStartTime_t > start_times
void CorrectForELoss(DVector3 &position, DVector3 &momentum, double q, double my_mass)
const DTrackingData & GetFitParameters(void) const