Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DTrackWireBased_factory_THROWN.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DTrackWireBased_factory_THROWN.cc
4 // Created: Mon Sep 3 19:57:11 EDT 2007
5 // Creator: davidl (on Darwin Amelia.local 8.10.1 i386)
6 //
7 
8 #include <cmath>
9 #include <limits>
10 using namespace std;
11 
12 #include <DANA/DApplication.h>
13 
14 #include <CDC/DCDCTrackHit.h>
15 #include <FDC/DFDCPseudo.h>
16 
18 #include "DMCThrown.h"
19 #include "DReferenceTrajectory.h"
20 #include "DRandom.h"
21 #include "DMatrix.h"
22 #include "DTrackHitSelector.h"
23 
24 
25 //------------------
26 // DTrackWireBased_factory_THROWN
27 //------------------
29 {
30  fitter = NULL;
31  hitselector=NULL;
32 
33  RootGeom = NULL;
34  geom = NULL;
35 }
36 
37 //------------------
38 // brun
39 //------------------
40 jerror_t DTrackWireBased_factory_THROWN::brun(jana::JEventLoop *loop, int32_t runnumber)
41 {
42  // Get pointer to DTrackFitter object that actually fits a track
43  vector<const DTrackFitter *> fitters;
44  loop->Get(fitters);
45  if(fitters.size()<1){
46  _DBG_<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl;
47  return RESOURCE_UNAVAILABLE;
48  }
49 
50  // Drop the const qualifier from the DTrackFitter pointer (I'm surely going to hell for this!)
51  fitter = const_cast<DTrackFitter*>(fitters[0]);
52 
53  // Warn user if something happened that caused us NOT to get a fitter object pointer
54  if(!fitter){
55  _DBG_<<"ERROR: Unable to get a DTrackFitter object! Chisq for DTrackWireBased:THROWN will NOT be calculated!"<<endl;
56  return RESOURCE_UNAVAILABLE;
57  }
58 
59  // Get pointer to DTrackHitSelector object
60  vector<const DTrackHitSelector *> hitselectors;
61  loop->Get(hitselectors);
62  if(hitselectors.size()<1){
63  _DBG_<<"ERROR: Unable to get a DTrackHitSelector object! NO DTrackWireBased:THROWN objects will be created!"<<endl;
64  return RESOURCE_UNAVAILABLE;
65  }
66  hitselector = hitselectors[0];
67 
68  // Set DGeometry pointers so it can be used by the DReferenceTrajectory class
69  DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
70  geom = dapp->GetDGeometry(runnumber);
71 
72  // Get the particle ID algorithms
73  loop->GetSingle(dParticleID);
74 
75  // Set magnetic field pointer
76  bfield = dapp->GetBfield(runnumber);
77 
78  return NOERROR;
79 }
80 
81 //------------------
82 // evnt
83 //------------------
84 jerror_t DTrackWireBased_factory_THROWN::evnt(JEventLoop *loop, uint64_t eventnumber)
85 {
86  vector<const DMCThrown*> mcthrowns;
87  vector<const DCDCTrackHit*> cdctrackhits;
88  vector<const DFDCPseudo*> fdcpseudos;
89  loop->Get(mcthrowns);
90  loop->Get(cdctrackhits);
91  loop->Get(fdcpseudos);
92 
93  for(unsigned int i=0; i< mcthrowns.size(); i++){
94  const DMCThrown *thrown = mcthrowns[i];
95 
96  if(fabs(thrown->charge())<1)continue;
97 
98  // First, copy over the DKinematicData part
100  *static_cast<DKinematicData*>(track) = *static_cast<const DKinematicData*>(thrown);
101 
102  // Add DMCThrown as associated object
103  track->AddAssociatedObject(thrown);
104 
105  // We need to swim a reference trajectory here. To avoid the overhead
106  // of allocating/deallocating them every event, we keep a pool and
107  // re-use them. If the pool is not big enough, then add one to the
108  // pool.
109  unsigned int locNumInitialReferenceTrajectories = rt_pool.size();
110  if(rt_pool.size()<=_data.size()){
111  // This is a little ugly, but only gets called a few times throughout the life of the process
112  // Note: these never get deleted, even at the end of process.
113  rt_pool.push_back(new DReferenceTrajectory(bfield));
114  }
115  DReferenceTrajectory *rt = rt_pool[_data.size()];
116  if(locNumInitialReferenceTrajectories == rt_pool.size()) //didn't create a new one
117  rt->Reset();
118  rt->q = track->charge();
119 
120  DVector3 pos = track->position();
121  DVector3 mom = track->momentum();
122  rt->SetMass(thrown->mass());
123  rt->SetDGeometry(geom);
124  rt->SetDRootGeom(RootGeom);
125  rt->Swim(pos, mom, track->charge());
126 
127  // Find hits that should be on this track and add them as associated objects
128  vector<const DCDCTrackHit*> cdchits;
129  vector<const DFDCPseudo*> fdchits;
130  if(hitselector && cdctrackhits.size()>0)hitselector->GetCDCHits(DTrackHitSelector::kHelical, rt, cdctrackhits, cdchits);
131  if(hitselector && fdcpseudos.size()>0)hitselector->GetFDCHits(DTrackHitSelector::kHelical, rt, fdcpseudos, fdchits);
132  for(unsigned int i=0; i<cdchits.size(); i++)track->AddAssociatedObject(cdchits[i]);
133  for(unsigned int i=0; i<fdchits.size(); i++)track->AddAssociatedObject(fdchits[i]);
134 
135  // We want to get chisq and Ndof values for this track using the hits from above.
136  // We do this using the DTrackFitter object. This more or less guarantees that the
137  // chisq calculation is done in the same way as it is for track fitting. Note
138  // that no fitting is actually done here so this should be reasonably fast
139  if(fitter){
140  fitter->Reset();
141  fitter->AddHits(cdchits);
142  fitter->AddHits(fdchits);
143  double chisq;
144  int Ndof;
145  vector<DTrackFitter::pull_t> pulls;
146  fitter->ChiSq(DTrackFitter::kWireBased, rt, &chisq, &Ndof, &pulls);
147  track->chisq = chisq;
148  track->Ndof = Ndof;
149  track->FOM = TMath::Prob(track->chisq, track->Ndof);
150  track->pulls = pulls;
151  }else{
152  track->chisq = 0.0;
153  track->Ndof = 0;
154  track->FOM = numeric_limits<double>::quiet_NaN();
155  }
156 
157  //
158  //
159  // These are to provide symmetry with DTrackTimeBased_factory_THROWN.cc
160  //
161  //
162  track->candidateid = thrown->id;
163  //track->trackid = thrown->id;
164  //track->FOM = 1.0;
165 
166  _data.push_back(track);
167  }
168 
169  // Set CDC ring & FDC plane hit patterns
170  for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i)
171  {
172  vector<const DCDCTrackHit*> locCDCTrackHits;
173  _data[loc_i]->Get(locCDCTrackHits);
174 
175  vector<const DFDCPseudo*> locFDCPseudos;
176  _data[loc_i]->Get(locFDCPseudos);
177 
178  _data[loc_i]->dCDCRings = dParticleID->Get_CDCRingBitPattern(locCDCTrackHits);
179  _data[loc_i]->dFDCPlanes = dParticleID->Get_FDCPlaneBitPattern(locFDCPseudos);
180  }
181 
182  return NOERROR;
183 }
184 
Definition: track.h:16
DApplication * dapp
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
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.
Definition: DTrackFitter.h:61
TVector3 DVector3
Definition: DVector3.h:14
void AddHits(vector< const DCDCTrackHit * > cdchits)
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
void SetMass(double mass)
const DVector3 & position(void) const
void Swim(const DVector3 &pos, const DVector3 &mom, double q=-1000.0, const TMatrixFSym *cov=NULL, double smax=2000.0, const DCoordinateSystem *wire=NULL)
void SetDGeometry(const DGeometry *geom)
virtual double ChiSq(fit_type_t fit_type, DReferenceTrajectory *rt, double *chisq_ptr=NULL, int *dof_ptr=NULL, vector< pull_t > *pulls_ptr=NULL)=0
void Reset(void)
Definition: DTrackFitter.cc:94
vector< DTrackFitter::pull_t > pulls
Holds pulls used in chisq calc. (not including off-diagonals)
DGeometry * GetDGeometry(unsigned int run_number)
#define _DBG_
Definition: HDEVIO.h:12
void SetDRootGeom(const DRootGeom *RootGeom)
double charge(void) const
int Ndof
Number of degrees of freedom in the fit.
const DVector3 & momentum(void) const
const DTrackFitter * fitter
oid_t candidateid
which DTrackCandidate this came from
jerror_t brun(JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
double mass(void) const