Bug Summary

File:libraries/TRACKING/DTrackTimeBased_factory.cc
Location:line 499, column 7
Description:Dereference of null pointer

Annotated Source Code

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 <TMath.h>
13using namespace std;
14
15#define TOF_SIGMA0.080 0.080 // TOF resolution in ns
16
17#include <TROOT.h>
18
19#include "DTrackTimeBased_factory.h"
20#include <TRACKING/DTrackWireBased.h>
21#include <TRACKING/DReferenceTrajectory.h>
22#include <TRACKING/DTrackFitter.h>
23#include <TRACKING/DTrackHitSelector.h>
24#include <TRACKING/DMCTrackHit.h>
25#include <SplitString.h>
26
27using namespace jana;
28
29// Routine for sorting start times
30bool DTrackTimeBased_T0_cmp(DTrackTimeBased::DStartTime_t a,
31 DTrackTimeBased::DStartTime_t b){
32 return (a.system>b.system);
33}
34
35
36// count_common_members
37//------------------
38template<typename T>
39static unsigned int count_common_members(vector<T> &a, vector<T> &b)
40{
41 unsigned int n=0;
42 for(unsigned int i=0; i<a.size(); i++){
43 for(unsigned int j=0; j<b.size(); j++){
44 if(a[i]==b[j])n++;
45 }
46 }
47
48 return n;
49}
50
51
52
53//------------------
54// init
55//------------------
56jerror_t DTrackTimeBased_factory::init(void)
57{
58 fitter = NULL__null;
59 MAX_DReferenceTrajectoryPoolSize = 20;
60
61 //DEBUG_HISTS = false;
62 DEBUG_HISTS = true;
63 DEBUG_LEVEL = 0;
64 MOMENTUM_CUT_FOR_DEDX=0.5;
65 MOMENTUM_CUT_FOR_PROTON_ID=2.0;
66
67 MIN_CDC_HITS_FOR_TB_FORWARD_TRACKING=3;
68 BYPASS_TB_FOR_FORWARD_TRACKS=false;
69
70 //which BCAL reconstruction algorithm to use
71 USE_KLOE=1;
72
73 USE_HITS_FROM_WIREBASED_FIT=false;
74 gPARMS->SetDefaultParameter("TRKFIT:USE_HITS_FROM_WIREBASED_FIT",
75 USE_HITS_FROM_WIREBASED_FIT);
76
77 gPARMS->SetDefaultParameter("TRKFIT:BYPASS_TB_FOR_FORWARD_TRACKS",
78 BYPASS_TB_FOR_FORWARD_TRACKS);
79 gPARMS->SetDefaultParameter("TRKFIT:MIN_CDC_HITS_FOR_TB_FORWARD_TRACKING",
80 MIN_CDC_HITS_FOR_TB_FORWARD_TRACKING);
81
82
83 gPARMS->SetDefaultParameter("TRKFIT:DEBUG_HISTS", DEBUG_HISTS);
84 gPARMS->SetDefaultParameter("TRKFIT:DEBUG_LEVEL", DEBUG_LEVEL);
85 gPARMS->SetDefaultParameter("TRKFIT:MOMENTUM_CUT_FOR_DEDX",MOMENTUM_CUT_FOR_DEDX);
86 gPARMS->SetDefaultParameter("TRKFIT:MOMENTUM_CUT_FOR_PROTON_ID",MOMENTUM_CUT_FOR_PROTON_ID);
87
88 SKIP_MASS_HYPOTHESES_WIRE_BASED=false;
89 gPARMS->SetDefaultParameter("TRKFIT:SKIP_MASS_HYPOTHESES_WIRE_BASED",
90 SKIP_MASS_HYPOTHESES_WIRE_BASED);
91
92 if (SKIP_MASS_HYPOTHESES_WIRE_BASED){
93 string MASS_HYPOTHESES_POSITIVE = "0.13957,0.93827";
94 string MASS_HYPOTHESES_NEGATIVE = "0.13957";
95 gPARMS->SetDefaultParameter("TRKFIT:MASS_HYPOTHESES_POSITIVE", MASS_HYPOTHESES_POSITIVE);
96 gPARMS->SetDefaultParameter("TRKFIT:MASS_HYPOTHESES_NEGATIVE", MASS_HYPOTHESES_NEGATIVE);
97
98 // Parse MASS_HYPOTHESES strings to make list of masses to try
99 SplitString(MASS_HYPOTHESES_POSITIVE, mass_hypotheses_positive, ",");
100 SplitString(MASS_HYPOTHESES_NEGATIVE, mass_hypotheses_negative, ",");
101 if(mass_hypotheses_positive.size()==0)mass_hypotheses_positive.push_back(0.0); // If empty string is specified, assume they want massless particle
102 if(mass_hypotheses_negative.size()==0)mass_hypotheses_negative.push_back(0.0); // If empty string is specified, assume they want massless particle
103
104
105
106 }
107
108
109
110 // Forces correct particle id (when available)
111 PID_FORCE_TRUTH = false;
112 gPARMS->SetDefaultParameter("TRKFIT:PID_FORCE_TRUTH", PID_FORCE_TRUTH);
113
114 gPARMS->SetDefaultParameter("BCALRECON:USE_KLOE", USE_KLOE);
115
116 return NOERROR;
117}
118
119//------------------
120// brun
121//------------------
122jerror_t DTrackTimeBased_factory::brun(jana::JEventLoop *loop, int runnumber)
123{
124 // Get the geometry
125 DApplication* dapp=dynamic_cast<DApplication*>(loop->GetJApplication());
126 geom = dapp->GetDGeometry(runnumber);
127
128 // Get pointer to TrackFitter object that actually fits a track
129 vector<const DTrackFitter *> fitters;
130 loop->Get(fitters);
131 if(fitters.size()<1){
132 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
132<<" "<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
133 return RESOURCE_UNAVAILABLE;
134 }
135
136 // Drop the const qualifier from the DTrackFitter pointer (I'm surely going to hell for this!)
137 fitter = const_cast<DTrackFitter*>(fitters[0]);
138
139 // Warn user if something happened that caused us NOT to get a fitter object pointer
140 if(!fitter){
141 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
141<<" "<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
142 return RESOURCE_UNAVAILABLE;
143 }
144
145 // Get the particle ID algorithms
146 vector<const DParticleID *> pid_algorithms;
147 loop->Get(pid_algorithms);
148 if(pid_algorithms.size()<1){
149 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
149<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
150 return RESOURCE_UNAVAILABLE;
151 }
152 // Drop the const qualifier from the DParticleID pointer (I'm surely going to hell for this!)
153 pid_algorithm = const_cast<DParticleID*>(pid_algorithms[0]);
154
155 // Warn user if something happened that caused us NOT to get a pid_algorithm object pointer
156 if(!pid_algorithm){
157 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
157<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
158 return RESOURCE_UNAVAILABLE;
159 }
160
161
162 if(DEBUG_HISTS){
163 dapp->Lock();
164
165 // Histograms may already exist. (Another thread may have created them)
166 // Try and get pointers to the existing ones.
167 fom_chi2_trk = (TH1F*)gROOT->FindObject("fom_chi2_trk");
168 fom = (TH1F*)gROOT->FindObject("fom");
169 hitMatchFOM = (TH1F*)gROOT->FindObject("hitMatchFOM");
170 chi2_trk_mom = (TH2F*)gROOT->FindObject("chi2_trk_mom");
171
172 if(!fom_chi2_trk)fom_chi2_trk = new TH1F("fom_chi2_trk","PID FOM: #chi^{2}/Ndf from tracking", 1000, 0.0, 100.0);
173 if(!fom)fom = new TH1F("fom","Combined PID FOM", 1000, 0.0, 1.01);
174 if(!hitMatchFOM)hitMatchFOM = new TH1F("hitMatchFOM","Total Fraction of Hit Matches", 101, 0.0, 1.01);
175 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.);
176
177
178 Hstart_time=(TH2F*)gROOT->FindObject("Hstart_time");
179 if (!Hstart_time) Hstart_time=new TH2F("Hstart_time",
180 "vertex time source vs. time",
181 300,-50,50,9,-0.5,8.5);
182
183 dapp->Unlock();
184
185 }
186
187
188 return NOERROR;
189}
190
191//------------------
192// evnt
193//------------------
194jerror_t DTrackTimeBased_factory::evnt(JEventLoop *loop, int eventnumber)
195{
196 if(!fitter)return NOERROR;
197
198 if(rtv.size() > MAX_DReferenceTrajectoryPoolSize){
199 //printf("rtv Deleting\n");
200 for(size_t loc_i = MAX_DReferenceTrajectoryPoolSize; loc_i < rtv.size(); ++loc_i)
201 delete rtv[loc_i];
202 rtv.resize(MAX_DReferenceTrajectoryPoolSize);
203 }
204
205 // Get candidates and hits
206 vector<const DTrackWireBased*> tracks;
207 loop->Get(tracks);
208 if (tracks.size()==0) return NOERROR;
209
210 // get start counter hits
211 vector<const DSCHit*>sc_hits;
212 eventLoop->Get(sc_hits);
213
214 // Get TOF points
215 vector<const DTOFPoint*> tof_points;
216 eventLoop->Get(tof_points);
217
218 // Get BCAL and FCAL showers
219 vector<const DBCALShower*>bcal_showers;
220 if (USE_KLOE) {
221 eventLoop->Get(bcal_showers, "KLOE");
222 } else {
223 eventLoop->Get(bcal_showers);
224 }
225
226 vector<const DMCThrown*> mcthrowns;
227 loop->Get(mcthrowns);
228
229 // Loop over candidates
230 for(unsigned int i=0; i<tracks.size(); i++){
231 const DTrackWireBased *track = tracks[i];
232
233 if (SKIP_MASS_HYPOTHESES_WIRE_BASED){
234 // Choose list of mass hypotheses based on charge of candidate
235 vector<double> mass_hypotheses;
236 if(track->charge()<0.0){
237 mass_hypotheses = mass_hypotheses_negative;
238 }else{
239 mass_hypotheses = mass_hypotheses_positive;
240 }
241
242 for (unsigned int j=0;j<mass_hypotheses.size();j++){
243 if (mass_hypotheses[j]>0.9
244 && track->momentum().Mag()>MOMENTUM_CUT_FOR_PROTON_ID) continue;
245 DoFit(track,sc_hits,tof_points,bcal_showers,loop,mass_hypotheses[j]);
246 }
247 }
248 else{
249 if (BYPASS_TB_FOR_FORWARD_TRACKS){
250 // Lists of hits used in the previous pass
251 vector<const DCDCTrackHit *>cdchits;
252 track->GetT(cdchits);
253 vector<const DFDCPseudo *>fdchits;
254 track->GetT(fdchits);
255
256 if (fdchits.size()>0
257 && cdchits.size()<MIN_CDC_HITS_FOR_TB_FORWARD_TRACKING){
258 // Copy over the results of the wire-based fit to DTrackTimeBased
259 DTrackTimeBased *timebased_track = new DTrackTimeBased;
260
261 // Copy over DKinematicData part
262 DKinematicData *track_kd = timebased_track;
263 *track_kd = *track;
264
265 timebased_track->rt = track->rt;
266 timebased_track->chisq = track->chisq;
267 timebased_track->Ndof = track->Ndof;
268 timebased_track->pulls = track->pulls;
269 timebased_track->trackid = track->id;
270 timebased_track->candidateid=track->candidateid;
271
272 for(unsigned int m=0; m<fdchits.size(); m++)
273 timebased_track->AddAssociatedObject(fdchits[m]);
274 for(unsigned int m=0; m<cdchits.size(); m++)
275 timebased_track->AddAssociatedObject(cdchits[m]);
276
277 // Compute the figure-of-merit based on tracking
278 timebased_track->FOM = TMath::Prob(timebased_track->chisq, timebased_track->Ndof);
279
280 _data.push_back(timebased_track);
281 }
282 }
283 else{
284 unsigned int num=_data.size();
285 DoFit(track,sc_hits,tof_points,bcal_showers,loop,track->mass());
286
287 //_DBG_<< "eventnumber: " << eventnumber << endl;
288 if (PID_FORCE_TRUTH && _data.size()>num) {
289 // Add figure-of-merit based on difference between thrown and reconstructed momentum
290 // if more than half of the track's hits match MC truth hits and also (charge,mass)
291 // match; add FOM=0 otherwise
292 _data[_data.size()-1]->FOM=GetTruthMatchingFOM(i,_data[_data.size()-1],
293 mcthrowns);
294 }
295 }
296 }
297 }
298
299 // Filter out duplicate tracks
300 FilterDuplicates();
301
302 return NOERROR;
303}
304
305//------------------
306// erun
307//------------------
308jerror_t DTrackTimeBased_factory::erun(void)
309{
310 return NOERROR;
311}
312
313//------------------
314// fini
315//------------------
316jerror_t DTrackTimeBased_factory::fini(void)
317{
318 for(unsigned int i=0; i<rtv.size(); i++)delete rtv[i];
319 rtv.clear();
320
321 return NOERROR;
322}
323
324//------------------
325// FilterDuplicates
326//------------------
327void DTrackTimeBased_factory::FilterDuplicates(void)
328{
329 /// Look through all current DTrackTimeBased objects and remove any
330 /// that have all of their hits in common with another track
331
332 if(_data.size()==0)return;
333
334 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
334<<" "<<"Looking for clones of time-based tracks ..."<<endl;
335
336 set<unsigned int> indexes_to_delete;
337 for(unsigned int i=0; i<_data.size()-1; i++){
338 DTrackTimeBased *dtrack1 = _data[i];
339
340 vector<const DCDCTrackHit*> cdchits1;
341 vector<const DFDCPseudo*> fdchits1;
342 dtrack1->Get(cdchits1);
343 dtrack1->Get(fdchits1);
344
345 JObject::oid_t cand1=dtrack1->candidateid;
346 for(unsigned int j=i+1; j<_data.size(); j++){
347 DTrackTimeBased *dtrack2 = _data[j];
348 if (dtrack2->candidateid==cand1) continue;
349
350 // Particles with the same mass but from different
351 // candidates are filtered at the Wire-based level.
352 // Here, it's possible to have multiple tracks with
353 // different masses that are clones due to that.
354 // Hence, we cut different mass clones is appropriate.
355 //if (dtrack2->mass() != dtrack1->mass())continue;
356
357 vector<const DCDCTrackHit*> cdchits2;
358 vector<const DFDCPseudo*> fdchits2;
359 dtrack2->Get(cdchits2);
360 dtrack2->Get(fdchits2);
361
362 // Count number of cdc and fdc hits in common
363 unsigned int Ncdc = count_common_members(cdchits1, cdchits2);
364 unsigned int Nfdc = count_common_members(fdchits1, fdchits2);
365
366 if(DEBUG_LEVEL>3){
367 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
367<<" "<<"cand1:"<<cand1<<" cand2:"<<dtrack2->candidateid<<endl;
368 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
368<<" "<<" Ncdc="<<Ncdc<<" cdchits1.size()="<<cdchits1.size()<<" cdchits2.size()="<<cdchits2.size()<<endl;
369 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
369<<" "<<" Nfdc="<<Nfdc<<" fdchits1.size()="<<fdchits1.size()<<" fdchits2.size()="<<fdchits2.size()<<endl;
370 }
371
372 if(Ncdc!=cdchits1.size() && Ncdc!=cdchits2.size())continue;
373 if(Nfdc!=fdchits1.size() && Nfdc!=fdchits2.size())continue;
374
375 unsigned int total = Ncdc + Nfdc;
376 unsigned int total1 = cdchits1.size()+fdchits1.size();
377 unsigned int total2 = cdchits2.size()+fdchits2.size();
378 if(total!=total1 && total!=total2)continue;
379
380 if(total1<total2){
381 indexes_to_delete.insert(i);
382 }else if(total2<total1){
383 indexes_to_delete.insert(j);
384 }else if(dtrack1->FOM > dtrack2->FOM){
385 indexes_to_delete.insert(j);
386 }else{
387 indexes_to_delete.insert(i);
388 }
389 }
390 }
391
392 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
392<<" "<<"Found "<<indexes_to_delete.size()<<" time-based clones"<<endl;
393
394 // Return now if we're keeping everyone
395 if(indexes_to_delete.size()==0)return;
396
397 // Copy pointers that we want to keep to a new container and delete
398 // the clone objects
399 vector<DTrackTimeBased*> new_data;
400 for(unsigned int i=0; i<_data.size(); i++){
401 if(indexes_to_delete.find(i)==indexes_to_delete.end()){
402 new_data.push_back(_data[i]);
403 }else{
404 delete _data[i];
405 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
405<<" "<<"Deleting clone time-based track "<<i<<endl;
406 }
407 }
408 _data = new_data;
409}
410
411// Returns a FOM based on difference between thrown and reconstructed momentum if track matches MC truth information,
412// returns a FOM=0 otherwise;
413// a match requires identical masses and charges, and that more than half of the track's hits match the truth hits
414double DTrackTimeBased_factory::GetTruthMatchingFOM(int trackIndex,DTrackTimeBased *track,vector<const DMCThrown*>mcthrowns) {
415 bool match=false;
416
417 DLorentzVector fourMom = track->lorentzMomentum();
418 //DLorentzVector gen_fourMom[mcthrowns.size()];
419 vector<DLorentzVector> gen_fourMom(mcthrowns.size());
420 for(unsigned int i=0; i<mcthrowns.size(); i++){
421 gen_fourMom[i] = mcthrowns[i]->lorentzMomentum();
422 }
423
424 // Get info for thrown track
425 int MAX_TRACKS = (int)mcthrowns.size()+1, thrownIndex=-1; double f = 0.;
426 GetThrownIndex(track,MAX_TRACKS,f,thrownIndex);
427 if(thrownIndex<=0 || thrownIndex>=MAX_TRACKS || f<=0.5) return 0.;
428
429 double delta_pt_over_pt = (fourMom.Pt()-gen_fourMom[thrownIndex-1].Pt())/gen_fourMom[thrownIndex-1].Pt();
430 double delta_theta = (fourMom.Theta()-gen_fourMom[thrownIndex-1].Theta())*1000.0; // in milliradians
431 double delta_phi = (fourMom.Phi()-gen_fourMom[thrownIndex-1].Phi())*1000.0; // in milliradians
432 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);
433
434 if (fabs(track->mass()-mcthrowns[thrownIndex-1]->mass())<0.01 && track->charge()==mcthrowns[thrownIndex-1]->charge())
435 match = true;
436
437 double trk_chi2=track->chisq;
438 unsigned int ndof=track->Ndof;
439
440 if(DEBUG_HISTS&&match){
441 fom_chi2_trk->Fill(track->chisq);
442 chi2_trk_mom->Fill(chisq/3.,trk_chi2/ndof);
443 fom->Fill(TMath::Prob(chisq,3));
444 }
445
446 /*_DBG_ << "f: " << f << endl;
447 _DBG_ << "trk_chi2: " << trk_chi2 << endl;
448 _DBG_ << "ndof: " << ndof << endl;
449 _DBG_ << "throwncharge: " << mcthrowns[thrownIndex-1]->charge() << endl;
450 _DBG_ << "trackcharge: " << track->charge() << endl;
451 _DBG_ << "chargediff: " << fabs(track->charge()-mcthrowns[thrownIndex-1]->charge()) << endl;
452 _DBG_ << "thrownmass: " << mcthrowns[thrownIndex-1]->mass() << endl;
453 _DBG_ << "trackmass: " << track->mass() << endl;
454 _DBG_ << "massdiff: " << fabs(track->mass()-mcthrowns[thrownIndex-1]->mass()) << endl;
455 _DBG_ << "chisq: " << chisq << endl;
456 _DBG_ << "match?: " << match << endl;
457 _DBG_ << "thrownIndex: " << thrownIndex << " trackIndex: " << trackIndex << endl;
458 _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;
459 _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;*/
460
461 return (match) ? TMath::Prob(chisq,3) : 0.0;
462}
463
464//------------------
465// GetThrownIndex
466//------------------
467void DTrackTimeBased_factory::GetThrownIndex(const DKinematicData *kd, int &MAX_TRACKS, double &f, int &track)
468{
469 vector<const DCDCTrackHit*> cdctrackhits;
470 vector<const DFDCPseudo*> fdcpseudos;
471
472 // The DKinematicData object should be a DTrackCandidate, DTrackWireBased, or DParticle which
473 // has associated objects for the hits
474 kd->Get(cdctrackhits);
475 kd->Get(fdcpseudos);
476
477 // The track number is buried in the truth hit objects of type DMCTrackHit. These should be
478 // associated objects for the individual hit objects. We need to loop through them and
479 // keep track of how many hits for each track number we find
480
481 // CDC hits
482 vector<int> cdc_track_no(MAX_TRACKS, 0);
483 for(unsigned int i=0; i<cdctrackhits.size(); i++){
1
Loop condition is false. Execution continues on line 494
→
484 vector<const DMCTrackHit*> mctrackhits;
485 cdctrackhits[i]->Get(mctrackhits);
486 if(mctrackhits.size()==0)continue;
487 if(!mctrackhits[0]->primary)continue;
488 int track = mctrackhits[0]->track;
489 if(track>=0 && track<MAX_TRACKS)cdc_track_no[track]++;
490 //_DBG_ << "cdc:(i,trackhitssize,mchitssize,TrackNo,NhitsforTrackNo): " << "(" << i << "," << cdctrackhits.size() << "," << mctrackhits.size() << "," << track << "," << cdc_track_no[track] << ")" << endl;
491 //_DBG_ << "cdc:(system,ptype,r,phi,z): " << "(" << mctrackhits[0]->system << "," << mctrackhits[0]->ptype << "," << mctrackhits[0]->r << "," << mctrackhits[0]->phi << "," << mctrackhits[0]->z << ")" << endl;
492 }
493 // FDC hits
494 vector<int> fdc_track_no(MAX_TRACKS, 0);
495 for(unsigned int i=0; i<fdcpseudos.size(); i++){
2
←
Loop condition is true. Entering loop body
→
5
←
Loop condition is true. Entering loop body
→
8
←
Loop condition is true. Entering loop body
→
496 vector<const DMCTrackHit*> mctrackhits;
497 fdcpseudos[i]->Get(mctrackhits);
498 if(mctrackhits.size()==0)continue;
3
←
Taking true branch
→
4
←
 Execution continues on line 495
→
6
←
Taking true branch
→
7
←
 Execution continues on line 495
→
9
←
Taking false branch
→
499 if(!mctrackhits[0]->primary)continue;
10
←
Dereference of null pointer
500 int track = mctrackhits[0]->track;
501 if(track>=0 && track<MAX_TRACKS)fdc_track_no[track]++;
502 //_DBG_ << "fdc:(i,trackhitssize,mchitssize,TrackNo,NhitsforTrackNo): " << "(" << i << "," << fdcpseudos.size() << "," << mctrackhits.size() << "," << track << "," << fdc_track_no[track] << ")" << endl;
503 //_DBG_ << "fdc:(system,ptype,r,phi,z): " << "(" << mctrackhits[0]->system << "," << mctrackhits[0]->ptype << "," << mctrackhits[0]->r << "," << mctrackhits[0]->phi << "," << mctrackhits[0]->z << ")" << endl;
504 }
505
506 // Find track number with most wires hit
507 int track_with_max_hits = 0;
508 int tot_hits_max = cdc_track_no[0] + fdc_track_no[0];
509 for(int i=1; i<MAX_TRACKS; i++){
510 int tot_hits = cdc_track_no[i] + fdc_track_no[i];
511 if(tot_hits > tot_hits_max){
512 track_with_max_hits=i;
513 tot_hits_max = tot_hits;
514 }
515 //_DBG_ << "tot_hits_max: " << tot_hits_max << endl;
516 //_DBG_ << "track_with_max_hits: " << track_with_max_hits << endl;
517 }
518
519 int Ncdc = cdc_track_no[track_with_max_hits];
520 int Nfdc = fdc_track_no[track_with_max_hits];
521
522 // total fraction of reconstructed hits that match truth hits
523 if (cdctrackhits.size()+fdcpseudos.size()) f = 1.*(Ncdc+Nfdc)/(cdctrackhits.size()+fdcpseudos.size());
524 //_DBG_ << "(Ncdc(match),Nfdc(match),Ncdc(recon),Nfdc(recon)): " << "(" << Ncdc << "," << Nfdc << "," << cdctrackhits.size() << "," << fdcpseudos.size() << ")" << endl;
525 if(DEBUG_HISTS)hitMatchFOM->Fill(f);
526
527 // If there are no hits on this track, then we really should report
528 // a "non-track" (i.e. track=-1)
529 track = tot_hits_max>0 ? track_with_max_hits:-1;
530}
531
532
533// Create a list of start (vertex) times from various sources, ordered by
534// uncertainty.
535void DTrackTimeBased_factory
536 ::CreateStartTimeList(const DTrackWireBased *track,
537 vector<const DSCHit*>&sc_hits,
538 vector<const DTOFPoint*>&tof_points,
539 vector<const DBCALShower*>&bcal_showers,
540 vector<DTrackTimeBased::DStartTime_t>&start_times){
541 // Add the t0 estimate from the tracking
542 DTrackTimeBased::DStartTime_t start_time;
543 start_time.t0=track->t0();
544 start_time.t0_sigma=5.;
545 start_time.system=track->t0_detector();
546 start_times.push_back(start_time);
547
548 // Match to the start counter and the outer detectors
549 double tproj=track->t0(); // initial guess from tracking
550 unsigned int tof_id=0,sc_id=0;
551 double locPathLength, locFlightTime;
552
553 if (pid_algorithm->MatchToSC(track->rt,DTrackFitter::kWireBased,sc_hits,
554 tproj,sc_id, locPathLength, locFlightTime)==NOERROR){
555 // Fill in the start time vector
556 start_time.t0=tproj;
557 start_time.t0_sigma=0.3;
558 start_time.system=SYS_START;
559 start_times.push_back(start_time);
560 }
561
562 if (pid_algorithm->MatchToTOF(track->rt,DTrackFitter::kWireBased,tof_points,
563 tproj,tof_id, locPathLength, locFlightTime)==NOERROR){
564 // Fill in the start time vector
565 start_time.t0=tproj;
566 start_time.t0_sigma=0.1;
567 start_time.system=SYS_TOF;
568 start_times.push_back(start_time);
569 }
570 vector<const DBCALShower*> locMatchedBCALShowers;
571 if (pid_algorithm->MatchToBCAL(track->rt, bcal_showers, locMatchedBCALShowers, tproj, locPathLength, locFlightTime) == NOERROR){
572 // Fill in the start time vector
573 start_time.t0=tproj;
574 start_time.t0_sigma=0.5;
575 start_time.system=SYS_BCAL;
576 start_times.push_back(start_time);
577 }
578 // Sort the list of start times according to uncertainty and set
579 // t0 for the fit to the first entry
580 sort(start_times.begin(),start_times.end(),DTrackTimeBased_T0_cmp);
581 mStartTime=start_times[0].t0;
582 mStartDetector=start_times[0].system;
583
584 // for (unsigned int i=0;i<start_times.size();i++){
585 // printf("%d t0 %f sys %d\n",i,start_times[i].t0,start_times[i].system);
586 // }
587
588}
589
590// Create a list of start times and do the fit for a particular mass hypothesis
591void DTrackTimeBased_factory::DoFit(const DTrackWireBased *track,
592 vector<const DSCHit*>&sc_hits,
593 vector<const DTOFPoint*>&tof_points,
594 vector<const DBCALShower*>&bcal_showers,
595 JEventLoop *loop,
596 double mass){
597 // Create vector of start times from various sources
598 vector<DTrackTimeBased::DStartTime_t>start_times;
599 CreateStartTimeList(track,sc_hits,tof_points,bcal_showers,start_times);
600
601 // Make sure there are enough DReferenceTrajectory objects
602 unsigned int locNumInitialReferenceTrajectories = rtv.size();
603 while(rtv.size()<=_data.size()){
604 //printf("TB adding\n");
605 rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap()));
606 }
607 DReferenceTrajectory *rt = rtv[_data.size()];
608 if(locNumInitialReferenceTrajectories == rtv.size()) //didn't create a new one
609 rt->Reset();
610 rt->SetMass(mass);
611 rt->q = track->charge();
612
613 rt->SetDGeometry(geom);
614
615 if(DEBUG_LEVEL>1){_DBG__std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
615<<std::endl;_DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
615<<" "<<"---- Starting time based fit with mass: "<<mass<<endl;}
616
617 // Do the fit
618 DTrackFitter::fit_status_t status = DTrackFitter::kFitNotDone;
619 if (USE_HITS_FROM_WIREBASED_FIT) {
620 fitter->Reset();
621 fitter->SetFitType(DTrackFitter::kTimeBased);
622
623 // Get the hits from the wire-based track
624 vector<const DFDCPseudo*>myfdchits;
625 track->GetT(myfdchits);
626 fitter->AddHits(myfdchits);
627 vector<const DCDCTrackHit *>mycdchits;
628 track->GetT(mycdchits);
629 fitter->AddHits(mycdchits);
630
631 status=fitter->FitTrack(track->position(),track->momentum(),
632 track->charge(),mass,mStartTime,mStartDetector);
633 }
634 else{
635 fitter->SetFitType(DTrackFitter::kTimeBased);
636 status = fitter->FindHitsAndFitTrack(*track, track->rt,loop, mass, mStartTime,
637 mStartDetector);
638 // If the status is kFitNotDone, then no hits were attached to this track
639 // using the hit-gathering algorithm. In this case get the hits from the
640 // wire-based track
641 if (status==DTrackFitter::kFitNotDone){
642 //_DBG_ << " Using wire-based hits " << endl;
643
644 vector<const DFDCPseudo*>myfdchits;
645 track->GetT(myfdchits);
646 fitter->AddHits(myfdchits);
647 vector<const DCDCTrackHit *>mycdchits;
648 track->GetT(mycdchits);
649 fitter->AddHits(mycdchits);
650
651 status=fitter->FitTrack(track->position(),track->momentum(),
652 track->charge(),mass,mStartTime,mStartDetector);
653 }
654
655 }
656
657 // Check the status value from the fit
658 switch(status){
659 case DTrackFitter::kFitNotDone:
660 _DBG_std::cerr<<"DTrackTimeBased_factory.cc"<<":"<<
660<<" "<<"Fitter returned kFitNotDone. This should never happen!!"<<endl;
661 case DTrackFitter::kFitFailed:
662 break;
663 case DTrackFitter::kFitSuccess:
664 case DTrackFitter::kFitNoImprovement:
665 {
666 // Allocate a DReferenceTrajectory object if needed.
667 // These each have a large enough memory footprint that
668 // it causes noticable performance problems if we allocated
669 // and deallocated them every event. Therefore, we allocate
670 // when needed, but recycle them on the next event.
671 // They are deleted in the fini method.
672 locNumInitialReferenceTrajectories = rtv.size();
673 while(rtv.size()<=_data.size())rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap()));
674 DReferenceTrajectory *rt = rtv[_data.size()];
675 if(locNumInitialReferenceTrajectories == rtv.size()) //didn't create a new one
676 rt->Reset();
677
678 // Create a new time-based track object
679 DTrackTimeBased *timebased_track = new DTrackTimeBased;
680
681 // Copy over DKinematicData part
682 DKinematicData *track_kd = timebased_track;
683 *track_kd = fitter->GetFitParameters();
684 rt->SetMass(mass);
685 rt->SetDGeometry(geom);
686 rt->q = timebased_track->charge();
687 rt->Swim(timebased_track->position(), timebased_track->momentum(), timebased_track->charge());
688
689 timebased_track->rt = rt;
690 timebased_track->chisq = fitter->GetChisq();
691 timebased_track->Ndof = fitter->GetNdof();
692 timebased_track->pulls = fitter->GetPulls();
693 timebased_track->trackid = track->id;
694 timebased_track->candidateid=track->candidateid;
695
696 // Set the start time and add the list of start times
697 timebased_track->setT0(mStartTime,start_times[0].t0_sigma,
698 mStartDetector);
699 timebased_track->start_times.assign(start_times.begin(),
700 start_times.end());
701
702 if (DEBUG_HISTS){
703 if (start_times[0].system==SYS_START){
704 Hstart_time->Fill(start_times[0].t0,0.);
705 }
706 else{
707 Hstart_time->Fill(start_times[0].t0,start_times[0].system);
708 }
709
710 }
711
712
713 // Add hits used as associated objects
714 const vector<const DCDCTrackHit*> &cdchits = fitter->GetCDCFitHits();
715 const vector<const DFDCPseudo*> &fdchits = fitter->GetFDCFitHits();
716
717 for(unsigned int m=0; m<cdchits.size(); m++)
718 timebased_track->AddAssociatedObject(cdchits[m]);
719 for(unsigned int m=0; m<fdchits.size(); m++)
720 timebased_track->AddAssociatedObject(fdchits[m]);
721
722 // dEdx
723 double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdx_CDC;
724 unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC;
725 pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdx_CDC, locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC);
726
727 timebased_track->ddEdx_FDC = locdEdx_FDC;
728 timebased_track->ddx_FDC = locdx_FDC;
729 timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC;
730 timebased_track->ddEdx_CDC = locdEdx_CDC;
731 timebased_track->ddx_CDC = locdx_CDC;
732 timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC;
733 timebased_track->setdEdx((locNumHitsUsedFordEdx_CDC >= locNumHitsUsedFordEdx_FDC) ? locdEdx_CDC : locdEdx_FDC);
734
735 // Add DTrack object as associate object
736 timebased_track->AddAssociatedObject(track);
737
738 // Compute the figure-of-merit based on tracking
739 timebased_track->FOM = TMath::Prob(timebased_track->chisq, timebased_track->Ndof);
740 //_DBG_<< "FOM: " << timebased_track->FOM << endl;
741
742 _data.push_back(timebased_track);
743 break;
744
745 }
746 default:
747 break;
748 }
749}
750