Bug Summary

File:libraries/TRACKING/DTrackWireBased_factory.cc
Location:line 251, column 5
Description:Value stored to 'Ntracks_to_fit' is never read

Annotated Source Code

1// $Id: DTrackWireBased_factory.cc 5612 2009-10-15 20:51:25Z staylor $
2//
3// File: DTrackWireBased_factory.cc
4// Created: Wed Sep 3 09:33:40 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 <cmath>
13using namespace std;
14
15#include "DTrackWireBased_factory.h"
16#include <TRACKING/DTrackCandidate.h>
17#include <TRACKING/DReferenceTrajectory.h>
18#include <CDC/DCDCTrackHit.h>
19#include <FDC/DFDCPseudo.h>
20#include <SplitString.h>
21
22#include <TROOT.h>
23
24#define C_EFFECTIVE15. 15.
25
26using namespace jana;
27
28//------------------
29// CDCSortByRincreasing
30//------------------
31bool CDCSortByRincreasing(const DCDCTrackHit* const &hit1, const DCDCTrackHit* const &hit2) {
32 // use the ring number to sort by R(decreasing) and then straw(increasing)
33 if(hit1->wire->ring == hit2->wire->ring){
34 return hit1->wire->straw < hit2->wire->straw;
35 }
36 return hit1->wire->ring < hit2->wire->ring;
37}
38
39//------------------
40// FDCSortByZincreasing
41//------------------
42bool FDCSortByZincreasing(const DFDCPseudo* const &hit1, const DFDCPseudo* const &hit2) {
43 // use the layer number to sort by Z(decreasing) and then wire(increasing)
44 if(hit1->wire->layer == hit2->wire->layer){
45 return hit1->wire->wire < hit2->wire->wire;
46 }
47 return hit1->wire->layer < hit2->wire->layer;
48}
49
50
51//------------------
52// count_common_members
53//------------------
54template<typename T>
55static unsigned int count_common_members(vector<T> &a, vector<T> &b)
56{
57 unsigned int n=0;
58 for(unsigned int i=0; i<a.size(); i++){
59 for(unsigned int j=0; j<b.size(); j++){
60 if(a[i]==b[j])n++;
61 }
62 }
63
64 return n;
65}
66
67//------------------
68// init
69//------------------
70jerror_t DTrackWireBased_factory::init(void)
71{
72 fitter = NULL__null;
73 MAX_DReferenceTrajectoryPoolSize = 50;
74
75 //DEBUG_HISTS = true;
76 DEBUG_HISTS = false;
77 DEBUG_LEVEL = 0;
78
79 gPARMS->SetDefaultParameter("TRKFIT:DEBUG_LEVEL",DEBUG_LEVEL);
80 COSMICS=false;
81 gPARMS->SetDefaultParameter("TRKFIND:COSMICS",COSMICS);
82
83 return NOERROR;
84}
85
86//------------------
87// brun
88//------------------
89jerror_t DTrackWireBased_factory::brun(jana::JEventLoop *loop, int32_t runnumber)
90{
91 // Get the geometry
92 DApplication* dapp=dynamic_cast<DApplication*>(loop->GetJApplication());
93 geom = dapp->GetDGeometry(runnumber);
94 // Check for magnetic field
95 dIsNoFieldFlag = (dynamic_cast<const DMagneticFieldMapNoField*>(dapp->GetBfield(runnumber)) != NULL__null);
96
97 // Get pointer to DTrackFitter object that actually fits a track
98 vector<const DTrackFitter *> fitters;
99 loop->Get(fitters);
100 if(fitters.size()<1){
101 _DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<101<<" "<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
102 return RESOURCE_UNAVAILABLE;
103 }
104
105 // Drop the const qualifier from the DTrackFitter pointer (I'm surely going to hell for this!)
106 fitter = const_cast<DTrackFitter*>(fitters[0]);
107
108 // Warn user if something happened that caused us NOT to get a fitter object pointer
109 if(!fitter){
110 _DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<110<<" "<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
111 return RESOURCE_UNAVAILABLE;
112 }
113 SKIP_MASS_HYPOTHESES_WIRE_BASED=false;
114 gPARMS->SetDefaultParameter("TRKFIT:SKIP_MASS_HYPOTHESES_WIRE_BASED",
115 SKIP_MASS_HYPOTHESES_WIRE_BASED);
116
117 USE_HITS_FROM_CANDIDATE=false;
118 gPARMS->SetDefaultParameter("TRKFIT:USE_HITS_FROM_CANDIDATE",
119 USE_HITS_FROM_CANDIDATE);
120
121 MIN_FIT_P = 0.050; // GeV
122 gPARMS->SetDefaultParameter("TRKFIT:MIN_FIT_P", MIN_FIT_P, "Minimum fit momentum in GeV/c for fit to be considered successful");
123
124 if (SKIP_MASS_HYPOTHESES_WIRE_BASED==false){
125
126 ostringstream locMassStream_Positive, locMassStream_Negative;
127 locMassStream_Positive << ParticleMass(PiPlus) << "," << ParticleMass(KPlus) << "," << ParticleMass(Proton);
128 locMassStream_Negative << ParticleMass(PiMinus) << "," << ParticleMass(KMinus);
129 string MASS_HYPOTHESES_POSITIVE = locMassStream_Positive.str();
130 string MASS_HYPOTHESES_NEGATIVE = locMassStream_Negative.str();
131 gPARMS->SetDefaultParameter("TRKFIT:MASS_HYPOTHESES_POSITIVE", MASS_HYPOTHESES_POSITIVE);
132 gPARMS->SetDefaultParameter("TRKFIT:MASS_HYPOTHESES_NEGATIVE", MASS_HYPOTHESES_NEGATIVE);
133
134 // Parse MASS_HYPOTHESES strings to make list of masses to try
135 SplitString(MASS_HYPOTHESES_POSITIVE, mass_hypotheses_positive, ",");
136 SplitString(MASS_HYPOTHESES_NEGATIVE, mass_hypotheses_negative, ",");
137 if(mass_hypotheses_positive.size()==0)mass_hypotheses_positive.push_back(0.0); // If empty string is specified, assume they want massless particle
138 if(mass_hypotheses_negative.size()==0)mass_hypotheses_negative.push_back(0.0); // If empty string is specified, assume they want massless particle
139
140
141 }
142 if(DEBUG_HISTS){
143 dapp->Lock();
144
145 // Histograms may already exist. (Another thread may have created them)
146 // Try and get pointers to the existing ones.
147
148
149 dapp->Unlock();
150 }
151
152 // Get the particle ID algorithms
153 loop->GetSingle(dPIDAlgorithm);
154
155 //Pre-allocate memory for DReferenceTrajectory objects early
156 //The swim-step objects of these arrays take up a significant amount of memory, and it can be difficult to find enough free contiguous space for them.
157 //Therefore, allocate them at the beginning before the available memory becomes randomly populated
158 while(rtv.size() < MAX_DReferenceTrajectoryPoolSize)
159 rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap()));
160
161 return NOERROR;
162}
163
164//------------------
165// evnt
166//------------------
167jerror_t DTrackWireBased_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
168{
169 if(!fitter)return NOERROR;
170
171 if(rtv.size() > MAX_DReferenceTrajectoryPoolSize){
172 for(size_t loc_i = MAX_DReferenceTrajectoryPoolSize; loc_i < rtv.size(); ++loc_i)
173 delete rtv[loc_i];
174 rtv.resize(MAX_DReferenceTrajectoryPoolSize);
175 }
176
177 // Get candidates and hits
178 vector<const DTrackCandidate*> candidates;
179 if(COSMICS){
180 loop->Get(candidates,"CDCCOSMIC");
181 }
182 else{
183 loop->Get(candidates);
184 }
185 if (candidates.size()==0) return NOERROR;
186
187 if (dIsNoFieldFlag){
188 // Copy results over from the StraightLine or CDCCOSMIC candidate and add reference
189 // trajectory
190 for (unsigned int i=0;i<candidates.size();i++){
191 const DTrackCandidate *cand=candidates[i];
192
193 // Make a new wire-based track
194 DTrackWireBased *track = new DTrackWireBased;
195
196 // Copy over DKinematicData part
197 DKinematicData *track_kd = track;
198 *track_kd=*cand;
199
200 // Attach a reference trajectory -- make sure there are enough DReferenceTrajectory objects
201 unsigned int locNumInitialReferenceTrajectories = rtv.size();
202 while(rtv.size()<=num_used_rts){
203 //printf("Adding %d %d\n",rtv.size(),_data.size());
204 rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap()));
205 }
206 DReferenceTrajectory *rt = rtv[num_used_rts];
207 if(locNumInitialReferenceTrajectories == rtv.size()) //didn't create a new one
208 rt->Reset();
209 rt->SetDGeometry(geom);
210 rt->FastSwim(track->position(),track->momentum(),track->charge());
211 track->rt=rt;
212
213 // candidate id
214 track->candidateid=i+1;
215
216 // Track quality properties
217 track->Ndof=cand->Ndof;
218 track->chisq=cand->chisq;
219 track->FOM = TMath::Prob(track->chisq, track->Ndof);
220
221 // pull vector
222 track->pulls = cand->pulls;
223
224 // Lists of hits used in the previous pass
225 vector<const DCDCTrackHit *>cdchits;
226 cand->GetT(cdchits);
227 vector<const DFDCPseudo *>fdchits;
228 cand->GetT(fdchits);
229
230 for (unsigned int k=0;k<cdchits.size();k++){
231 track->AddAssociatedObject(cdchits[k]);
232 }
233 for (unsigned int k=0;k<fdchits.size();k++){
234 track->AddAssociatedObject(fdchits[k]);
235 }
236
237 _data.push_back(track);
238
239 }
240 return NOERROR;
241 }
242
243
244
245 // Reset the number of used reference trajectories from the pool
246 num_used_rts=0;
247
248 // Count the number of tracks we'll be fitting
249 unsigned int Ntracks_to_fit = 0;
250 if (SKIP_MASS_HYPOTHESES_WIRE_BASED){
251 Ntracks_to_fit=candidates.size();
Value stored to 'Ntracks_to_fit' is never read
252 }
253 else{
254 for(unsigned int i=0; i<candidates.size(); i++){
255 Ntracks_to_fit += candidates[i]->charge()<0.0 ? mass_hypotheses_negative.size():mass_hypotheses_positive.size();
256 }
257 }
258
259 // Loop over candidates
260 for(unsigned int i=0; i<candidates.size(); i++){
261 const DTrackCandidate *candidate = candidates[i];
262
263 // Skip candidates with momentum below some cutoff
264 if (candidate->momentum().Mag()<MIN_FIT_P){
265 continue;
266 }
267
268 if (SKIP_MASS_HYPOTHESES_WIRE_BASED){
269 // Make sure there are enough DReferenceTrajectory objects
270 unsigned int locNumInitialReferenceTrajectories = rtv.size();
271 while(rtv.size()<=num_used_rts){
272 //printf("Adding %d %d\n",rtv.size(),_data.size());
273 rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap()));
274 }
275 DReferenceTrajectory *rt = rtv[num_used_rts];
276 if(locNumInitialReferenceTrajectories == rtv.size()) //didn't create a new one
277 rt->Reset();
278 rt->SetDGeometry(geom);
279 rt->q = candidate->charge();
280
281 // Increment the number of used reference trajectories
282 num_used_rts++;
283
284 DoFit(i,candidate,rt,loop,0.13957);
285 }
286 else{
287 // Choose list of mass hypotheses based on charge of candidate
288 vector<double> mass_hypotheses;
289 if(candidate->charge()<0.0){
290 mass_hypotheses = mass_hypotheses_negative;
291 }else{
292 mass_hypotheses = mass_hypotheses_positive;
293 }
294
295 if ((!isfinite(candidate->momentum().Mag())) || (!isfinite(candidate->position().Mag())))
296 _DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<296<<" " << "Invalid seed data for event "<< eventnumber <<"..."<<endl;
297
298 // Loop over potential particle masses
299 for(unsigned int j=0; j<mass_hypotheses.size(); j++){
300 if(DEBUG_LEVEL>1){_DBG__std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<300<<std::endl;_DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<300<<" "<<"---- Starting wire based fit with mass: "<<mass_hypotheses[j]<<endl;}
301 // Make sure there are enough DReferenceTrajectory objects
302 unsigned int locNumInitialReferenceTrajectories = rtv.size();
303 while(rtv.size()<=num_used_rts){
304 //printf("Adding %d\n",rtv.size());
305 rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap()));
306 }
307 DReferenceTrajectory *rt = rtv[num_used_rts];
308 if(locNumInitialReferenceTrajectories == rtv.size()) //didn't create a new one
309 rt->Reset();
310 rt->SetDGeometry(geom);
311 rt->q = candidate->charge();
312
313 // Increment the number of used reference trajectories
314 num_used_rts++;
315
316 DoFit(i,candidate,rt,loop,mass_hypotheses[j]);
317 }
318
319 }
320 }
321
322 // Filter out duplicate tracks
323 FilterDuplicates();
324
325 // Set CDC ring & FDC plane hit patterns
326 for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
327 {
328 vector<const DCDCTrackHit*> locCDCTrackHits;
329 _data[loc_i]->Get(locCDCTrackHits);
330
331 vector<const DFDCPseudo*> locFDCPseudos;
332 _data[loc_i]->Get(locFDCPseudos);
333
334 _data[loc_i]->dCDCRings = dPIDAlgorithm->Get_CDCRingBitPattern(locCDCTrackHits);
335 _data[loc_i]->dFDCPlanes = dPIDAlgorithm->Get_FDCPlaneBitPattern(locFDCPseudos);
336 }
337
338 return NOERROR;
339}
340
341
342//------------------
343// erun
344//------------------
345jerror_t DTrackWireBased_factory::erun(void)
346{
347 return NOERROR;
348}
349
350//------------------
351// fini
352//------------------
353jerror_t DTrackWireBased_factory::fini(void)
354{
355 for(unsigned int i=0; i<rtv.size(); i++)delete rtv[i];
356 rtv.clear();
357
358 return NOERROR;
359}
360
361//------------------
362// FilterDuplicates
363//------------------
364void DTrackWireBased_factory::FilterDuplicates(void)
365{
366 /// Look through all current DTrackWireBased objects and remove any
367 /// that have all of their hits in common with another track
368
369 if(_data.size()==0)return;
370
371 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<371<<" "<<"Looking for clones of wire-based tracks ..."<<endl;
372
373 set<unsigned int> indexes_to_delete;
374 for(unsigned int i=0; i<_data.size()-1; i++){
375 DTrackWireBased *dtrack1 = _data[i];
376
377 vector<const DCDCTrackHit*> cdchits1;
378 vector<const DFDCPseudo*> fdchits1;
379 dtrack1->Get(cdchits1);
380 dtrack1->Get(fdchits1);
381
382 JObject::oid_t cand1=dtrack1->candidateid;
383 for(unsigned int j=i+1; j<_data.size(); j++){
384 DTrackWireBased *dtrack2 = _data[j];
385 if (dtrack2->candidateid==cand1) continue;
386 if (dtrack2->mass() != dtrack1->mass())continue;
387
388 vector<const DCDCTrackHit*> cdchits2;
389 vector<const DFDCPseudo*> fdchits2;
390 dtrack2->Get(cdchits2);
391 dtrack2->Get(fdchits2);
392
393 // Count number of cdc and fdc hits in common
394 unsigned int Ncdc = count_common_members(cdchits1, cdchits2);
395 unsigned int Nfdc = count_common_members(fdchits1, fdchits2);
396 unsigned int total = Ncdc + Nfdc;
397
398 if (total==0) continue;
399 if(Ncdc!=cdchits1.size() && Ncdc!=cdchits2.size())continue;
400 if(Nfdc!=fdchits1.size() && Nfdc!=fdchits2.size())continue;
401
402 unsigned int total1 = cdchits1.size()+fdchits1.size();
403 unsigned int total2 = cdchits2.size()+fdchits2.size();
404 if(total!=total1 && total!=total2)continue;
405
406 if(total1<total2){
407 // The two track candidates actually correspond to
408 // a single track. Set the candidate id for this
409 // track to the candidate id from the clone match to
410 // prevent multiple clone tracks confusing matters
411 // at a later stage of the reconstruction...
412 _data[j]->candidateid=cand1;
413 indexes_to_delete.insert(i);
414 }else{
415 indexes_to_delete.insert(j);
416 }
417
418 }
419 }
420
421 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<421<<" "<<"Found "<<indexes_to_delete.size()<<" wire-based clones"<<endl;
422
423 // Return now if we're keeping everyone
424 if(indexes_to_delete.size()==0)return;
425
426 // Copy pointers that we want to keep to a new container and delete
427 // the clone objects
428 vector<DTrackWireBased*> new_data;
429 for(unsigned int i=0; i<_data.size(); i++){
430 if(indexes_to_delete.find(i)==indexes_to_delete.end()){
431 new_data.push_back(_data[i]);
432 }else{
433 delete _data[i];
434 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<434<<" "<<"Deleting clone wire-based track "<<i<<endl;
435 }
436 }
437 _data = new_data;
438}
439
440// Routine to find the hits, do the fit, and fill the list of wire-based tracks
441void DTrackWireBased_factory::DoFit(unsigned int c_id,
442 const DTrackCandidate *candidate,
443 DReferenceTrajectory *rt,
444 JEventLoop *loop, double mass){
445 // Do the fit
446 DTrackFitter::fit_status_t status = DTrackFitter::kFitNotDone;
447 if (USE_HITS_FROM_CANDIDATE) {
448 fitter->Reset();
449 fitter->SetFitType(DTrackFitter::kWireBased);
450
451 // Get the hits from the track candidate
452 vector<const DFDCPseudo*>myfdchits;
453 candidate->GetT(myfdchits);
454 fitter->AddHits(myfdchits);
455 vector<const DCDCTrackHit *>mycdchits;
456 candidate->GetT(mycdchits);
457 fitter->AddHits(mycdchits);
458
459 status=fitter->FitTrack(candidate->position(),candidate->momentum(),
460 candidate->charge(),mass,0.);
461 }
462 else{
463 fitter->SetFitType(DTrackFitter::kWireBased);
464 // Swim a reference trajectory using the candidate starting momentum
465 // and position
466 rt->SetMass(mass);
467 //rt->Swim(candidate->position(),candidate->momentum(),candidate->charge());
468 rt->FastSwim(candidate->position(),candidate->momentum(),candidate->charge(),2000.0,0.,370.);
469
470 status=fitter->FindHitsAndFitTrack(*candidate,rt,loop,mass,candidate->Ndof+3);
471 if (/*false && */status==DTrackFitter::kFitNotDone){
472 if (DEBUG_LEVEL>1)_DBG_std::cerr<<"libraries/TRACKING/DTrackWireBased_factory.cc"
<<":"<<472<<" " << "Using hits from candidate..." << endl;
473 fitter->Reset();
474
475 // Get the hits from the candidate
476 vector<const DFDCPseudo*>myfdchits;
477 candidate->GetT(myfdchits);
478 fitter->AddHits(myfdchits);
479 vector<const DCDCTrackHit *>mycdchits;
480 candidate->GetT(mycdchits);
481 fitter->AddHits(mycdchits);
482
483 status=fitter->FitTrack(candidate->position(),candidate->momentum(),
484 candidate->charge(),mass,0.);
485 }
486
487
488 }
489
490 // Check the status of the fit
491 switch(status){
492 case DTrackFitter::kFitNotDone:
493 //_DBG_<<"Fitter returned kFitNotDone. This should never happen!!"<<endl;
494 case DTrackFitter::kFitFailed:
495 break;
496 case DTrackFitter::kFitNoImprovement:
497 case DTrackFitter::kFitSuccess:
498 if(!isfinite(fitter->GetFitParameters().position().X())) break;
499 {
500 // Make a new wire-based track
501 DTrackWireBased *track = new DTrackWireBased;
502
503 // Copy over DKinematicData part
504 DKinematicData *track_kd = track;
505 *track_kd = fitter->GetFitParameters();
506 track_kd->setPID(dPIDAlgorithm->IDTrack(track_kd->charge(), track_kd->mass()));
507 track_kd->setTime(track_kd->t0());
508
509 // Fill reference trajectory
510 rt->q = candidate->charge();
511 rt->SetMass(track_kd->mass());
512 //rt->Swim(track->position(), track->momentum(), track->charge());
513 rt->FastSwim(track->position(), track->momentum(), track->charge());
514
515 if(rt->Nswim_steps <= 1)
516 {
517 //Track parameters are bogus (e.g. track position closer to infinity than the beamline)
518 delete track;
519 return;
520 }
521
522 track->rt = rt;
523 track->chisq = fitter->GetChisq();
524 track->Ndof = fitter->GetNdof();
525 track->FOM = TMath::Prob(track->chisq, track->Ndof);
526 track->pulls = fitter->GetPulls();
527 track->candidateid = c_id+1;
528
529 // Add hits used as associated objects
530 vector<const DCDCTrackHit*> cdchits = fitter->GetCDCFitHits();
531 vector<const DFDCPseudo*> fdchits = fitter->GetFDCFitHits();
532 sort(cdchits.begin(), cdchits.end(), CDCSortByRincreasing);
533 sort(fdchits.begin(), fdchits.end(), FDCSortByZincreasing);
534 for(unsigned int m=0; m<cdchits.size(); m++)track->AddAssociatedObject(cdchits[m]);
535 for(unsigned int m=0; m<fdchits.size(); m++)track->AddAssociatedObject(fdchits[m]);
536
537 // Add DTrackCandidate as associated object
538 track->AddAssociatedObject(candidate);
539
540 _data.push_back(track);
541 break;
542 }
543 default:
544 break;
545 }
546}