Bug Summary

File:libraries/TRACKING/DTrackFitterALT1.cc
Location:line 1303, column 14
Description:Dereference of null pointer

Annotated Source Code

1// $Id: DTrackFitterALT1.cc 9024 2012-04-18 18:26:47Z pmatt $
2//
3// File: DTrackFitterALT1.cc
4// Created: Tue Sep 2 11:18:22 EDT 2008
5// Creator: davidl
6//
7namespace{const char* GetMyID(){return "$Id: DTrackFitterALT1.cc 9024 2012-04-18 18:26:47Z pmatt $";}}
8
9#include <math.h>
10
11#include <TROOT.h>
12
13#include <DVector3.h>
14#include <DMatrix.h>
15
16#include <JANA/JEventLoop.h>
17using namespace jana;
18
19#include "GlueX.h"
20#include "DANA/DApplication.h"
21#include "DMagneticFieldStepper.h"
22#include "DTrackCandidate.h"
23#include "DTrackFitterALT1.h"
24#include "CDC/DCDCTrackHit.h"
25#include "FDC/DFDCPseudo.h"
26#include "DReferenceTrajectory.h"
27#include "DMCThrown.h"
28#include "DMCTrackHit.h"
29#include "DTrackHitSelectorTHROWN.h"
30
31extern double GetCDCCovariance(int layer1, int layer2);
32extern double GetFDCCovariance(int layer1, int layer2);
33extern double GetFDCCathodeCovariance(int layer1, int layer2);
34
35
36#define NaNstd::numeric_limits<double>::quiet_NaN() std::numeric_limits<double>::quiet_NaN()
37
38// The GNU implementation of STL includes definitions of "greater" and "less"
39// but the SunOS implementation does not. Since it is a bit of a pain to
40// define this only for SunOS, we just define "greaterthan" and use it for
41// all platforms/compilers. Note that this is essentially the same as the
42// GNU definition from stl_function.h, except it does not derive from the
43// templated "binary_function" class.
44template<typename T>
45class greaterthan{
46 public: bool operator()(const T &a, const T &b) const {return a>b;}
47};
48
49
50//------------------
51// DTrackFitterALT1 (Constructor)
52//------------------
53DTrackFitterALT1::DTrackFitterALT1(JEventLoop *loop):DTrackFitter(loop)
54{
55 this->loop = loop;
56
57 // Define target center
58 target = new DCoordinateSystem();
59 target->origin.SetXYZ(0.0, 0.0, 65.0);
60 target->sdir.SetXYZ(1.0, 0.0, 0.0);
61 target->tdir.SetXYZ(0.0, 1.0, 0.0);
62 target->udir.SetXYZ(0.0, 0.0, 1.0);
63 target->L = 30.0;
64
65 // DReferenceTrajectory objects
66 rt = new DReferenceTrajectory(bfield);
67 tmprt = new DReferenceTrajectory(bfield);
68
69 DEBUG_HISTS = false;
70 DEBUG_LEVEL = 0;
71 MAX_CHISQ_DIFF = 1.0E-2;
72 MAX_FIT_ITERATIONS = 50;
73 SIGMA_CDC = 0.0150;
74 CDC_USE_PARAMETERIZED_SIGMA = true;
75 SIGMA_FDC_ANODE = 0.0200;
76 SIGMA_FDC_CATHODE = 0.0200;
77 CHISQ_GOOD_LIMIT = 2.0;
78 LEAST_SQUARES_DP = 0.0001;
79 LEAST_SQUARES_DX = 0.010;
80 LEAST_SQUARES_MIN_HITS = 3;
81 LEAST_SQUARES_MAX_E2NORM = 1.0E20;
82 DEFAULT_MASS = rt->GetMass(); // Get default mass from DReferenceTrajectory class itself
83 TARGET_CONSTRAINT = false;
84 LR_FORCE_TRUTH = false;
85 USE_MULS_COVARIANCE = true;
86 USE_FDC = true;
87 USE_CDC = true;
88 USE_FDC_CATHODE = true;
89 MATERIAL_MAP_MODEL = "DGeometry";
90
91 gPARMS->SetDefaultParameter("TRKFIT:DEBUG_HISTS", DEBUG_HISTS);
92 gPARMS->SetDefaultParameter("TRKFIT:DEBUG_LEVEL", DEBUG_LEVEL);
93 gPARMS->SetDefaultParameter("TRKFIT:MAX_CHISQ_DIFF", MAX_CHISQ_DIFF);
94 gPARMS->SetDefaultParameter("TRKFIT:MAX_FIT_ITERATIONS", MAX_FIT_ITERATIONS);
95 gPARMS->SetDefaultParameter("TRKFIT:SIGMA_CDC", SIGMA_CDC);
96 gPARMS->SetDefaultParameter("TRKFIT:CDC_USE_PARAMETERIZED_SIGMA", CDC_USE_PARAMETERIZED_SIGMA);
97 gPARMS->SetDefaultParameter("TRKFIT:SIGMA_FDC_ANODE", SIGMA_FDC_ANODE);
98 gPARMS->SetDefaultParameter("TRKFIT:SIGMA_FDC_CATHODE", SIGMA_FDC_CATHODE);
99 gPARMS->SetDefaultParameter("TRKFIT:CHISQ_GOOD_LIMIT", CHISQ_GOOD_LIMIT);
100 gPARMS->SetDefaultParameter("TRKFIT:LEAST_SQUARES_DP", LEAST_SQUARES_DP);
101 gPARMS->SetDefaultParameter("TRKFIT:LEAST_SQUARES_DX", LEAST_SQUARES_DX);
102 gPARMS->SetDefaultParameter("TRKFIT:LEAST_SQUARES_MIN_HITS", LEAST_SQUARES_MIN_HITS);
103 gPARMS->SetDefaultParameter("TRKFIT:LEAST_SQUARES_MAX_E2NORM", LEAST_SQUARES_MAX_E2NORM);
104 gPARMS->SetDefaultParameter("TRKFIT:DEFAULT_MASS", DEFAULT_MASS);
105 gPARMS->SetDefaultParameter("TRKFIT:TARGET_CONSTRAINT", TARGET_CONSTRAINT);
106 gPARMS->SetDefaultParameter("TRKFIT:LR_FORCE_TRUTH", LR_FORCE_TRUTH);
107 gPARMS->SetDefaultParameter("TRKFIT:USE_MULS_COVARIANCE", USE_MULS_COVARIANCE);
108 gPARMS->SetDefaultParameter("TRKFIT:USE_FDC", USE_FDC);
109 gPARMS->SetDefaultParameter("TRKFIT:USE_CDC", USE_CDC);
110 gPARMS->SetDefaultParameter("TRKFIT:USE_FDC_CATHODE", USE_FDC_CATHODE);
111 gPARMS->SetDefaultParameter("TRKFIT:MATERIAL_MAP_MODEL", MATERIAL_MAP_MODEL);
112
113 // Set default mass based on configuration parameter
114 rt->SetMass(DEFAULT_MASS);
115 tmprt->SetMass(DEFAULT_MASS);
116
117 // Set the geometry object pointers based on the material map
118 // specified via configuration parameter.
119 if(MATERIAL_MAP_MODEL=="DRootGeom"){
120 rt->SetDRootGeom(RootGeom);
121 tmprt->SetDRootGeom(RootGeom);
122 }else if(MATERIAL_MAP_MODEL=="DGeometry"){
123 rt->SetDGeometry(geom);
124 tmprt->SetDGeometry(geom);
125 }else if(MATERIAL_MAP_MODEL=="NONE"){
126 rt->SetDRootGeom(NULL__null);
127 tmprt->SetDRootGeom(NULL__null);
128 rt->SetDGeometry(NULL__null);
129 tmprt->SetDGeometry(NULL__null);
130 }else{
131 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<131<<
" "<<"ERROR: Invalid value for TRKFIT:MATERIAL_MAP_MODEL (=\""<<MATERIAL_MAP_MODEL<<"\")"<<endl;
132 exit(-1);
133 }
134
135 cout<<__FILE__"DTrackFitterALT1.cc"<<":"<<__LINE__135<<"-------------- Least Squares TRACKING --------------"<<endl;
136
137 DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
138 if(!dapp){
139 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<139<<
" "<<"Cannot get DApplication from JEventLoop! (are you using a JApplication based program?)"<<endl;
140 return;
141 }
142
143 if(DEBUG_HISTS){
144 dapp->Lock();
145
146 // Histograms may already exist. (Another thread may have created them)
147 // Try and get pointers to the existing ones.
148 cdcdoca_vs_dist = (TH2F*)gROOT->FindObject("cdcdoca_vs_dist");
149 cdcdoca_vs_dist_vs_ring = (TH3F*)gROOT->FindObject("cdcdoca_vs_dist_vs_ring");
150 fdcdoca_vs_dist = (TH2F*)gROOT->FindObject("fdcdoca_vs_dist");
151 fdcu_vs_s = (TH2F*)gROOT->FindObject("fdcu_vs_s");
152 dist_axial = (TH1F*)gROOT->FindObject("dist_axial");
153 doca_axial = (TH1F*)gROOT->FindObject("doca_axial");
154 dist_stereo = (TH1F*)gROOT->FindObject("dist_stereo");
155 doca_stereo = (TH1F*)gROOT->FindObject("doca_stereo");
156 chisq_final_vs_initial = (TH2F*)gROOT->FindObject("chisq_final_vs_initial");
157 nhits_final_vs_initial = (TH2F*)gROOT->FindObject("nhits_final_vs_initial");
158 Npasses = (TH1F*)gROOT->FindObject("Npasses");
159 ptotal = (TH1F*)gROOT->FindObject("ptotal");
160 residuals_cdc = (TH2F*)gROOT->FindObject("residuals_cdc");
161 residuals_fdc_anode = (TH2F*)gROOT->FindObject("residuals_fdc_anode");
162 residuals_fdc_cathode = (TH2F*)gROOT->FindObject("residuals_fdc_cathode");
163 residuals_cdc_vs_s = (TH3F*)gROOT->FindObject("residuals_cdc_vs_s");
164 residuals_fdc_anode_vs_s = (TH3F*)gROOT->FindObject("residuals_fdc_anode_vs_s");
165 residuals_fdc_cathode_vs_s = (TH3F*)gROOT->FindObject("residuals_fdc_cathode_vs_s");
166 initial_chisq_vs_Npasses = (TH2F*)gROOT->FindObject("initial_chisq_vs_Npasses");
167 chisq_vs_pass = (TH2F*)gROOT->FindObject("chisq_vs_pass");
168 dchisq_vs_pass = (TH2F*)gROOT->FindObject("dchisq_vs_pass");
169 cdc_single_hit_prob = (TH1F*)gROOT->FindObject("cdc_single_hit_prob");
170 cdc_double_hit_prob = (TH1F*)gROOT->FindObject("cdc_double_hit_prob");
171 fdc_single_hit_prob = (TH1F*)gROOT->FindObject("fdc_single_hit_prob");
172 fdc_double_hit_prob = (TH1F*)gROOT->FindObject("fdc_double_hit_prob");
173 cdc_can_resi = (TH1F*)gROOT->FindObject("cdc_can_resi");
174 fdc_can_resi = (TH1F*)gROOT->FindObject("fdc_can_resi");
175 fdc_can_resi_cath = (TH1F*)gROOT->FindObject("fdc_can_resi_cath");
176 chisq_vs_p_vs_theta = (TH2F*)gROOT->FindObject("chisq_vs_p_vs_theta");
177 lambda = (TH1F*)gROOT->FindObject("lambda");
178
179 if(!cdcdoca_vs_dist)cdcdoca_vs_dist = new TH2F("cdcdoca_vs_dist","DOCA vs. DIST",300, 0.0, 1.2, 300, 0.0, 1.2);
180 if(!cdcdoca_vs_dist_vs_ring)cdcdoca_vs_dist_vs_ring = new TH3F("cdcdoca_vs_dist_vs_ring","DOCA vs. DIST vs. ring",300, 0.0, 1.2, 300, 0.0, 1.2,23,0.5,23.5);
181 if(!fdcdoca_vs_dist)fdcdoca_vs_dist = new TH2F("fdcdoca_vs_dist","DOCA vs. DIST",500, 0.0, 2.0, 500, 0.0, 2.0);
182 if(!fdcu_vs_s)fdcu_vs_s = new TH2F("fdcu_vs_s","DOCA vs. DIST along wire",500, -60.0, 60.0, 500, -60.0, 60.0);
183 if(!dist_axial)dist_axial = new TH1F("dist_axial","Distance from drift time for axial CDC wires",300,0.0,3.0);
184 if(!doca_axial)doca_axial = new TH1F("doca_axial","DOCA of track for axial CDC wires",300,0.0,3.0);
185 if(!dist_stereo)dist_stereo = new TH1F("dist_stereo","Distance from drift time for stereo CDC wires",300,0.0,3.0);
186 if(!doca_stereo)doca_stereo = new TH1F("doca_stereo","DOCA of track for axial CDC wires",300,0.0,3.0);
187 if(!chisq_final_vs_initial)chisq_final_vs_initial = new TH2F("chisq_final_vs_initial","Final vs. initial chi-sq.",200, 0.0, 10.0,50, 0.0, 2.5);
188 if(!nhits_final_vs_initial)nhits_final_vs_initial = new TH2F("nhits_final_vs_initial","Final vs. initial nhits used in chi-sq.",30, -0.5, 29.5, 30, -0.5, 29.5);
189 if(!Npasses)Npasses = new TH1F("Npasses","Npasses", 101, -0.5, 100.5);
190 if(!ptotal)ptotal = new TH1F("ptotal","ptotal",1000, 0.1, 8.0);
191 if(!residuals_cdc)residuals_cdc = new TH2F("residuals_cdc","Residuals in CDC",1000,-2.0,2.0,24,0.5,24.5);
192 if(!residuals_fdc_anode)residuals_fdc_anode = new TH2F("residuals_fdc_anode","Residuals in FDC anodes",1000,-2.0,2.0,24,0.5,24.5);
193 if(!residuals_fdc_cathode)residuals_fdc_cathode = new TH2F("residuals_fdc_cathode","Residuals in FDC cathode",1000,-2.0,2.0,24,0.5,24.5);
194 if(!residuals_cdc_vs_s)residuals_cdc_vs_s = new TH3F("residuals_cdc_vs_s","Residuals in CDC vs. pathlength",1000,-2.0,2.0,24,0.5,24.5,100, 0.0, 800);
195 if(!residuals_fdc_anode_vs_s)residuals_fdc_anode_vs_s = new TH3F("residuals_fdc_anode_vs_s","Residuals in FDC anode vs. pathlength",1000,-2.0,2.0,24,0.5,24.5,100, 0.0, 800);
196 if(!residuals_fdc_cathode_vs_s)residuals_fdc_cathode_vs_s = new TH3F("residuals_fdc_cathode_vs_s","Residuals in FDC cathode vs. pathlength",1000,-2.0,2.0,24,0.5,24.5,100, 0.0, 800);
197 if(!initial_chisq_vs_Npasses)initial_chisq_vs_Npasses = new TH2F("initial_chisq_vs_Npasses","Initial chi-sq vs. number of iterations", 25, -0.5, 24.5, 400, 0.0, 40.0);
198 if(!chisq_vs_pass)chisq_vs_pass = new TH2F("chisq_vs_pass","Chi-sq vs. fit iteration", 25, -0.5, 24.5, 400, 0.0, 40.0);
199 if(!dchisq_vs_pass)dchisq_vs_pass = new TH2F("dchisq_vs_pass","Change in Chi-sq vs. fit iteration", 25, -0.5, 24.5, 400, 0.0, 8.0);
200 if(!cdc_single_hit_prob)cdc_single_hit_prob = new TH1F("cdc_single_hit_prob","Probability a CDC hit belongs to a track for all tracks",200,0.0,1.0);
201 if(!cdc_double_hit_prob)cdc_double_hit_prob = new TH1F("cdc_double_hit_prob","Probability a CDC hit belongs to two tracks",200,0.0,1.0);
202 if(!fdc_single_hit_prob)fdc_single_hit_prob = new TH1F("fdc_single_hit_prob","Probability a FDC hit belongs to a track for all tracks",200,0.0,1.0);
203 if(!fdc_double_hit_prob)fdc_double_hit_prob = new TH1F("fdc_double_hit_prob","Probability a FDC hit belongs to two tracks",200,0.0,1.0);
204 if(!cdc_can_resi)cdc_can_resi = new TH1F("cdc_can_resi","Residual of CDC hits with candidate tracks", 200, -1.0, 1.0);
205 if(!fdc_can_resi)fdc_can_resi = new TH1F("fdc_can_resi","Residual of FDC hits with candidate tracks", 200, -1.0, 1.0);
206 if(!fdc_can_resi_cath)fdc_can_resi_cath = new TH1F("fdc_can_resi_cath","Residual of FDC cathode hits with candidate tracks", 200, -1.0, 1.0);
207 if(!lambda)lambda = new TH1F("lambda","Scaling factor #lambda for Newton-Raphson calculated step", 2048, -2.0, 2.0);
208
209 dapp->Unlock();
210 }
211}
212
213//------------------
214// DTrackFitterALT1 (Destructor)
215//------------------
216DTrackFitterALT1::~DTrackFitterALT1()
217{
218 if(rt)delete rt;
219 if(tmprt)delete tmprt;
220 if(target)delete target;
221}
222
223//------------------
224// FitTrack
225//------------------
226DTrackFitter::fit_status_t DTrackFitterALT1::FitTrack(void)
227{
228 /// Fit a track candidate
229
230 // Debug message
231 if(DEBUG_LEVEL>2){
232 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<232<<
" "<<"cdchits.size="<<cdchits.size()<<" fdchits.size="<<fdchits.size()<<endl;
233 if(fit_type==kTimeBased)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<233<<
" "<<"============= Time-based ==============="<<endl;
234 if(fit_type==kWireBased)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<234<<
" "<<"============= Wire-based ==============="<<endl;
235 }
236
237 // Set mass for our DReferenceTrajectory objects
238 rt->SetMass(input_params.mass());
239 tmprt->SetMass(input_params.mass());
240
241 // Do material boundary checking only for time-based tracking
242 // and only if beta*gamma is less than 1.0
243 double betagamma = input_params.momentum().Mag()/input_params.mass();
244 bool check_material_boundaries = fit_type==kTimeBased && betagamma<=1.0;
245 rt->SetCheckMaterialBoundaries(check_material_boundaries);
246 tmprt->SetCheckMaterialBoundaries(check_material_boundaries);
247
248 // Swim reference trajectory
249 fit_status = kFitNotDone; // initialize to a safe default
250 DVector3 mom = input_params.momentum();
251 if(mom.Mag()>12.0)mom.SetMag(8.0); // limit crazy big momenta from candidates
252 rt->Swim(input_params.position(), mom, input_params.charge());
253
254 if(DEBUG_LEVEL>1){
255 double chisq;
256 int Ndof;
257 ChiSq(fit_type, rt, &chisq, &Ndof);
258 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<258<<
" "<<"starting parameters for fit: chisq="<<chisq<<" Ndof="<<Ndof<<" (chisq/Ndof="<<chisq/(double)Ndof<<") p="<<rt->swim_steps[0].mom.Mag()<<endl;
259 }
260
261 // Iterate over fitter until it converges or we reach the
262 // maximum number of iterations
263 double last_chisq = 1.0E6;
264 int last_Ndof=1;
265 int Niterations;
266 for(Niterations=0; Niterations<MAX_FIT_ITERATIONS; Niterations++){
267
268 hitsInfo hinfo;
269 GetHits(fit_type, rt, hinfo);
270
271 // Optionally use the truth information from the Monte Carlo
272 // to force the correct LR choice for each hit. This is
273 // for debugging (obviously).
274 if(LR_FORCE_TRUTH && fit_type==kTimeBased)ForceLRTruth(loop, rt, hinfo);
275
276 switch(fit_status = LeastSquaresB(hinfo, rt)){
277 case kFitSuccess:
278 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<278<<
" "<<"Fit succeeded ----- "<<endl;
279 break;
280 case kFitNoImprovement:
281 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<281<<
" "<<"Unable to improve on input parameters (chisq="<<chisq<<")"<<endl;
282 break;
283 case kFitFailed:
284 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<284<<
" "<<"Fit failed on iteration "<<Niterations<<endl;
285 if(Niterations>0){
286 if(DEBUG_LEVEL>2)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<286<<
" "<<"Number of iterations >4. Trying to keep fit from last iteration... "<<endl;
287 chisq=last_chisq;
288 Ndof = last_Ndof;
289 fit_status = kFitNoImprovement;
290 break;
291 }
292 return fit_status;
293 break;
294 case kFitNotDone: // avoid compiler warnings
295 break;
296 }
297
298 // Fill debug histos
299 if(DEBUG_HISTS){
300 chisq_vs_pass->Fill(Niterations+1, chisq/Ndof);
301 dchisq_vs_pass->Fill(Niterations+1, last_chisq/last_Ndof - chisq/Ndof);
302 }
303
304 // If the chisq is too large, then consider it a hopeless cause
305 if(chisq/Ndof>1.0E4 || !finite(chisq/Ndof)){
306 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<306<<
" "<<"Fit chisq too large on iteration "<<Niterations<<endl;
307 return fit_status=kFitFailed;
308 }
309
310 // Check if we converged
311 if(fabs(last_chisq/last_Ndof-chisq/Ndof) < MAX_CHISQ_DIFF){
312 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<312<<
" "<<"Fit chisq change below threshold fabs("<<last_chisq/last_Ndof<<" - "<<chisq/Ndof<<")<"<<MAX_CHISQ_DIFF<<endl;
313 break;
314 }
315 last_chisq = chisq;
316 last_Ndof = Ndof;
317 }
318
319 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<319<<
" "<<" Niterations="<<Niterations<<" chisq="<<chisq<<" Ndof="<<Ndof<<" (chisq/Ndof="<<chisq/(double)Ndof<<") p="<<rt->swim_steps[0].mom.Mag()<<endl;
320
321 // At this point we must decide whether the fit succeeded or not.
322 // We'll consider the fit a success if any of the following is true:
323 //
324 // 1. We got through at least one iteration in the above loop
325 // 2. The chi-sq is less than CHISQ_GOOD_LIMIT
326 // 3. MAX_FIT_ITERATIONS is zero (for debugging)
327 bool fit_succeeded = false;
328 if(Niterations>0)fit_succeeded = true;
329 if(chisq<=CHISQ_GOOD_LIMIT)fit_succeeded = true;
330 if(MAX_FIT_ITERATIONS==0){fit_succeeded = true; fit_status=kFitSuccess;}
331 if(!fit_succeeded)return fit_status=kFitFailed;
332
333 if(DEBUG_LEVEL>9){
334 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<334<<
" "<<"-------- Final Chisq for track = "<<this->chisq<<" Ndof="<<this->Ndof<<endl;
335 double chisq;
336 int Ndof;
337 ChiSq(fit_type, rt, &chisq, &Ndof);
338 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<338<<
" "<<"-------- Check chisq = "<<chisq<<" Ndof="<<Ndof<<endl;
339 }
340
341 if(DEBUG_LEVEL>2){
342 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<342<<
" "<<"start POCA pos: "; rt->swim_steps->origin.Print();
343 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<343<<
" "<<"start POCA mom: "; rt->swim_steps->mom.Print();
344 ChiSq(fit_type, rt, &chisq, &Ndof);
345 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<345<<
" "<<"start POCA chisq="<<chisq<<" Ndof="<<Ndof<<" (chisq/Ndof="<<chisq/(double)Ndof<<") p="<<rt->swim_steps[0].mom.Mag()<<endl;
346 }
347
348 // At this point, the first point on the rt is likely not the POCA to the
349 // beamline. To find that, we swim the particle in, towards the target from
350 // the swim step closest to the first wire hit. We make sure to set the tmprt to
351 // energy *gain* for the backwards swim.
352 DVector3 vertex_pos = rt->swim_steps->origin;
353 DVector3 vertex_mom = rt->swim_steps->mom;
354 const DCoordinateSystem *first_wire = NULL__null;
355 if(cdchits.size()>0) first_wire=cdchits[0]->wire;
356 else if(fdchits.size()>0) first_wire=fdchits[0]->wire;
357 if(first_wire){
358 DReferenceTrajectory::swim_step_t *step = rt->FindClosestSwimStep(first_wire);
359 if(step){
360#if 0
361 tmprt->SetPLossDirection(DReferenceTrajectory::kBackward);
362 tmprt->Swim(step->origin, -step->mom, -rt->q, 100.0, target);
363 tmprt->DistToRT(target);
364 tmprt->GetLastDOCAPoint(vertex_pos, vertex_mom);
365 vertex_mom = -vertex_mom;
366 tmprt->SetPLossDirection(DReferenceTrajectory::kForward);
367#endif
368 // Swim back towards target (note: we may already be past it!)
369 tmprt->SetPLossDirection(DReferenceTrajectory::kBackward);
370 tmprt->Swim(rt->swim_steps->origin, -rt->swim_steps->mom, -rt->q, 100.0, target);
371
372 // Swim swim in farward direction to target (another short swim!)
373 tmprt->SetPLossDirection(DReferenceTrajectory::kForward);
374 vertex_pos = tmprt->swim_steps[tmprt->Nswim_steps-1].origin;
375 vertex_mom = tmprt->swim_steps[tmprt->Nswim_steps-1].mom;
376 tmprt->Swim(vertex_pos, -vertex_mom, rt->q, 100.0, target);
377
378 double locReturnValue = tmprt->DistToRT(target);
379 if((locReturnValue >= -1.0) || (locReturnValue <= 1.0)){ //else NaN
380 tmprt->GetLastDOCAPoint(vertex_pos, vertex_mom);
381
382 // Now, swim out the rt one last time such that it starts at the POCA to
383 // the beamline. We need to do this so that we can calculate the
384 // chisq/Ndof based on the vertex parameters
385 rt->Swim(vertex_pos, vertex_mom);
386 ChiSq(fit_type, rt, &this->chisq, &this->Ndof);
387 }
388 }
389 }else{
390 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<390<<
" "<<"NO WIRES IN CANDIDATE!! (event "<<eventnumber<<")"<<endl;
391 }
392
393 if(DEBUG_LEVEL>2){
394 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<394<<
" "<<"end POCA pos: "; rt->swim_steps->origin.Print();
395 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<395<<
" "<<"end POCA mom: "; rt->swim_steps->mom.Print();
396 ChiSq(fit_type, rt, &chisq, &Ndof);
397 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<397<<
" "<<"end POCA chisq="<<chisq<<" Ndof="<<Ndof<<" (chisq/Ndof="<<chisq/(double)Ndof<<") p="<<rt->swim_steps[0].mom.Mag()<<endl;
398 }
399
400
401#if 0
402 // At this point, the first point on the rt is likely not the POCA to the
403 // beamline. To find that, we swim the particle in, towards the target from
404 // about 10 cm out along the trajectory. We make sure to set the tmprt to
405 // energy gain for the backwards swim.
406 DReferenceTrajectory::swim_step_t *step = rt->swim_steps;
407 for(int i=0; i<rt->Nswim_steps-1; i++, step++){if(step->s>10.0)break;}
408 tmprt->SetPLossDirection(DReferenceTrajectory::kBackward);
409 tmprt->Swim(step->origin, -step->mom, -rt->q, 100.0, target);
410 tmprt->SetPLossDirection(DReferenceTrajectory::kForward);
411 double s;
412 tmprt->DistToRT(target, &s);
413 DVector3 vertex_mom, vertex_pos;
414 tmprt->GetLastDOCAPoint(vertex_pos, vertex_mom);
415 vertex_mom = -vertex_mom;
416#endif
417
418 // Copy final fit parameters into TrackFitter classes data members. Note that the chisq and Ndof
419 // members are copied in during the ChiSq() method call in LeastSquaresB().
420 fit_params.setPosition(vertex_pos);
421 fit_params.setMomentum(vertex_mom);
422 fit_params.setCharge(rt->q);
423 fit_params.setMass(input_params.mass());
424 cdchits_used_in_fit = cdchits;
425 fdchits_used_in_fit = fdchits;
426
427 // Fill debugging histos if requested
428 if(DEBUG_HISTS){
429 Npasses->Fill(Niterations);
430 initial_chisq_vs_Npasses->Fill(Niterations, chisq);
431 FillDebugHists(rt, vertex_pos, vertex_mom);
432 }
433
434 // Debugging messages
435 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<435<<
" "<<" -- Fit succeeded: Final params after vertex adjustment: q="<<rt->q<<" p="<<vertex_mom.Mag()<<" theta="<<vertex_mom.Theta()*57.3<<" phi="<<vertex_mom.Phi()*57.3<<" chisq="<<chisq<<" Ndof="<<Ndof<<" (chisq/Ndof="<<chisq/(double)Ndof<<")"<<endl;
436
437 return fit_status;
438}
439
440//------------------
441// ChiSq
442//------------------
443double DTrackFitterALT1::ChiSq(fit_type_t fit_type, DReferenceTrajectory *rt, double *chisq_ptr, int *dof_ptr, vector<pull_t> *pulls_ptr)
444{
445 /// Calculate the chisq for the given reference trajectory based on the hits
446 /// currently registered through the DTrackFitter base class into the cdchits
447 /// and fdchits vectors. This does not get called by the core part of the
448 /// fitter, but is used, rather, to give an "independent" value of the
449 /// chisq based only on a reference trajectory and the input hits. It is
450 /// called after the fit to calculate the final chisq.
451
452 hitsInfo hinfo;
453 GetHits(fit_type, rt, hinfo);
454 if(LR_FORCE_TRUTH && fit_type==kTimeBased)ForceLRTruth(loop, rt, hinfo);
455
456 vector<resiInfo> residuals;
457 GetResiInfo(rt, hinfo, residuals);
458
459 double my_chisq = ChiSq(residuals, chisq_ptr, dof_ptr);
460 if(pulls_ptr)*pulls_ptr = pulls;
461
462 return my_chisq;
463}
464
465//------------------
466// ChiSq
467//------------------
468double DTrackFitterALT1::ChiSq(vector<resiInfo> &residuals, double *chisq_ptr, int *dof_ptr)
469{
470 // The residuals vector contains a list of only good hits, both
471 // anode and cathode, along with their measurment errors. Use these to fill
472 // the resiv and cov_meas matrices now. Fill in the cov_muls below.
473 int Nmeasurements = (int)residuals.size();
474 resiv.ResizeTo(Nmeasurements, 1);
475 cov_meas.ResizeTo(Nmeasurements, Nmeasurements);
476 cov_muls.ResizeTo(Nmeasurements, Nmeasurements);
477 pulls.clear();
478
479 // Return "infinite" chisq if an empty residuals vector was passed.
480 if(Nmeasurements<1){
481 if(chisq_ptr)*chisq_ptr=1.0E6;
482 if(dof_ptr)*dof_ptr=1;
483 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<483<<
" "<<"Chisq called with empty residuals vector!"<<endl;
484 return 1.0E6;
485 }
486
487 for(unsigned int i=0; i<residuals.size(); i++){
488 resiInfo &ri = residuals[i];
489
490 resiv[i][0] = ri.resi;
491 cov_meas[i][i] = pow(ri.err, 2.0);
492 }
493
494 // Fill in the cov_muls matrix.
495 cov_muls.Zero();
496 if(USE_MULS_COVARIANCE){
497 // Find the first sensitive detector hit. Ignore
498 // the target point if present by checking layer!=0
499 const swim_step_t *step0 = NULL__null;
500 for(unsigned int i=0; i<residuals.size(); i++){
501 resiInfo &ri = residuals[i];
502 if(ri.layer<1)continue;
503 if(!ri.step)continue;
504 if( (step0==NULL__null) || (ri.step->s < step0->s) ){
505 step0 = ri.step;
506 }
507 }
508
509 if(step0){
510 // Loop over all residuals
511 for(unsigned int i=0; i<residuals.size(); i++){
512 const swim_step_t *stepA = residuals[i].step;
513 if(!stepA)continue;
514 // resiInfo &riA = residuals[i]; // see 2010/07/12 note below
515
516 // Loop over all residuals
517 for(unsigned int j=i; j<residuals.size(); j++){
518 const swim_step_t *stepB = residuals[j].step;
519 if(!stepB)continue;
520 // resiInfo &riB = residuals[j]; // see 2010/07/12 note below
521
522 // Correlations are very hard to get correct given the left-right
523 // ambiguity which determines the sign of the correlation.
524 // Because of this, only include the diagonal elements of the
525 // covariance matrix. Attempts were made to do otherwise, but they
526 // failed. The code is left here in case we need to try it again
527 // at a later time.
528 if(i!=j)continue;
529
530 // Correlations between A and B are due only to material between
531 // the first detector and the most upstream of A or B.
532 double sA = stepA->s;
533 double sB = stepB->s;
534 const swim_step_t *step_end = sA<sB ? stepA:stepB;
535
536 if(step0->s>step_end->s)continue; // Bullet proof
537
538 // Need a comment here to explain this ...
539 double itheta02 = step_end->itheta02 - step0->itheta02;
540 double itheta02s = step_end->itheta02s - step0->itheta02s;
541 double itheta02s2 = step_end->itheta02s2 - step0->itheta02s2;
542 double sigmaAB = sA*sB*itheta02 - (sA+sB)*itheta02s + itheta02s2;
543
544 // sigmaAB represents the magnitude of the covariance, but not the
545 // sign.
546 //
547 // For wires generally in the same direction, the sign
548 // should be positive if they are on the same side of the wire
549 // but negative if they are on opposite sides.
550 //
551 // For wires that are orhogonal (like a CDC is to an FDC wire), the
552 // the procedure is not so clear.
553
554 // Note: 2010/07/12 DL
555 // The following code is unneeded since we are not implementing off
556 // diagonal elements of the covariance matrix. It was noted by Simon,
557 // however, that the value of "sign" was sometimes being calculated
558 // as negative, even for on diagonal elements leading to overall
559 // worse resolutions. Upon re-reading the comment below, I'm so sure
560 // I believe the method to be correct anyway. Hence, I'm commenting
561 // it and the corresponding code out, but leaving it here in case this
562 // issue gets revisited in the future.
563 //
564 // // Find LR of this hit.
565 // // Vector pointing from origin of wire to step position crossed into
566 // // wire direction gives a vector that will either be pointing
567 // // generally in the direction of the momentum or opposite to it.
568 // // Take dot product of above described vectors for each hit
569 // // use it to determine L or R.
570 // DVector3 crossA = (stepA->origin - riA.hit->wire->origin).Cross(riA.hit->wire->udir);
571 // DVector3 crossB = (stepB->origin - riB.hit->wire->origin).Cross(riB.hit->wire->udir);
572 // double sides_angle = crossA.Angle(crossB)*57.3;
573 // double sign = fabs(sides_angle) < 90.0 ? +1.0:-1.0;
574 double sign = +1.0;
575
576 cov_muls[i][j] = cov_muls[j][i] = sign*sigmaAB;
577 }
578 }
579 }
580 }
581
582 // Multiply resiv with inverse of total error matrix to get chisq = r^T (cov_meas+cov_muls)^-1 r)
583 DMatrix resiv_t(DMatrix::kTransposed, resiv);
584 DMatrix W(DMatrix::kInverted, cov_meas+cov_muls);
585 DMatrix chisqM(resiv_t*W*resiv);
586 double chisq = chisqM[0][0];
587 int Ndof = (int)residuals.size() - 5; // assume 5 fit parameters
588 weights.ResizeTo(W);
589 weights = W;
590
591 // Copy pulls into pulls vector
592 for(unsigned int i=0; i<residuals.size(); i++){
593 double err = sqrt(cov_meas[i][i]+cov_muls[i][i]);
594 pulls.push_back(pull_t(resiv[i][0], err, residuals[i].step->s));
595 }
596
597 if(DEBUG_LEVEL>100 || (DEBUG_LEVEL>10 && !finite(chisq)))PrintChisqElements(resiv, cov_meas, cov_muls, weights);
598
599 // If the caller supplied pointers to chisq and dof variables, copy the values into them
600 if(chisq_ptr)*chisq_ptr = chisq;
601 if(dof_ptr)*dof_ptr = Ndof; // assume 5 fit parameters
602
603 // If it turns out the dof is zero, return 1.0E6. Otherwise, return
604 // the chisq/dof
605 return Ndof<2 ? 1.0E6:chisq/(double)Ndof;
606}
607
608//------------------
609// GetHits
610//------------------
611void DTrackFitterALT1::GetHits(fit_type_t fit_type, DReferenceTrajectory *rt, hitsInfo &hinfo)
612{
613
614 // -- Target --
615 if(TARGET_CONSTRAINT){
616 hitInfo hi;
617 hi.wire = target;
618 hi.dist = 0.0;
619 hi.err = 0.1; // 1mm beam width
620 hi.u_dist = 0.0;
621 hi.u_err = 0.0;
622 hi.good = hi.good_u = false;
623 hinfo.push_back(hi);
624 }
625
626 // --- CDC ---
627 if(USE_CDC){
628 for(unsigned int i=0; i<cdchits.size(); i++){
629 const DCDCTrackHit *hit = cdchits[i];
630 const DCoordinateSystem *wire = hit->wire;
631 hitInfo hi;
632
633 hi.wire = wire;
634 hi.u_dist = 0.0;
635 hi.u_err = 0.0;
636 hi.good = hi.good_u = false;
637
638 // Fill in shifts and errs vectors based on whether we're doing
639 // hit-based or time-based tracking
640 if(fit_type==kWireBased){
641 // If we're doing hit-based tracking then only the wire positions
642 // are used and the drift time info is ignored.
643 // NOTE: The track quality itself goes into the residual resoultion
644 // and so we use something larger than the variance over an evenly
645 // illuminated cell size (commented out). The value of 0.31 is
646 // empirical from forward (2-40 degree) pi+ tracks when MULS was
647 // off. This will likely need to be higher when MULS is on...
648 hi.dist = 0.0;
649 hi.err = 0.45; // empirical. (see note above)
650 //hi.err = 0.8/sqrt(12.0); // variance for evenly illuminated straw
651 }else{
652 // Find the DOCA point for the RT and use the momentum direction
653 // there and the wire direction to determine the "shift".
654 // Note that whether the shift is to the left or to the right
655 // is not decided here. That ambiguity is left to be resolved later
656 // by applying a minus sign (or not) to the shift.
657 DVector3 pos_doca, mom_doca;
658 double s;
659 rt->DistToRT(wire, &s);
660 rt->GetLastDOCAPoint(pos_doca, mom_doca);
661
662 // The magnitude of the shift is based on the drift time. The
663 // value of the dist member of the DCDCTrackHit object does not
664 // subtract out the TOF. This can add 50-100 microns to the
665 // resolution in the CDC. Here, we actually can calculate the TOF
666 // (for a given mass hypothesis).
667 double mass = rt->GetMass();
668 double beta = 1.0/sqrt(1.0 + pow(mass/mom_doca.Mag(), 2.0))*2.998E10;
669 double tof = s/beta/1.0E-9; // in ns
670 hi.dist = hit->dist*((hit->tdrift-tof)/hit->tdrift);
671 hi.err = SIGMA_CDC;
672
673 // Default is to use parameterized sigma, by allow user to specify
674 // using constant sigma.
675 if(CDC_USE_PARAMETERIZED_SIGMA){
676 // The following is from a fit to Yves' numbers circa Aug 2009. The values fit were
677 // resolution (microns) vs. drift distance (mm).
678 // par[8] = {699.875, -559.056, 149.391, 25.6929, -22.0238, 4.75091, -0.452373, 0.0163858};
679 double x = hi.dist;
680 if(x<0.0)x=0.0; // this can be negative due to tof subtraction above.
681 //double sigma_d = (699.875) + x*((-559.056) + x*((149.391) + x*((25.6929) + x*((-22.0238) + x*((4.75091) + x*((-0.452373) + x*((0.0163858))))))));
682 double sigma_d = 108.55 + 7.62391*x + 556.176*exp(-(1.12566)*pow(x,1.29645));
683 hi.err = sigma_d/10000.0;
684 }
685 }
686 hinfo.push_back(hi);
687 }
688 } // USE_CDC
689
690 // --- FDC ---
691 if(USE_FDC){
692 for(unsigned int i=0; i<fdchits.size(); i++){
693 const DFDCPseudo *hit = fdchits[i];
694 const DCoordinateSystem *wire = hit->wire;
695 hitInfo hi;
696
697 hi.wire = wire;
698 hi.good = hi.good_u = false;
699
700 // Fill in shifts and errs vectors based on whether we're doing
701 // hit-based or time-based tracking
702 if(fit_type==kWireBased){
703 // If we're doing hit-based tracking then only the wire positions
704 // are used and the drift time info is ignored.
705 // NOTE: The track quality itself goes into the residual resoultion
706 // and so we use something larger than the variance over an evenly
707 // illuminated cell size (commented out). The value of 0.30 is
708 // empirical from forward (2-40 degree) pi+ tracks when MULS was
709 hi.dist = 0.0;
710 hi.err = 0.30; // emprical. (see note above)
711 //hi.err = 0.5/sqrt(12.0); // variance for evenly illuminated cell - anode
712 hi.u_dist = 0.0;
713 hi.u_err = 0.0; // variance for evenly illuminated cell - cathode
714 }else{
715 // Find the DOCA point for the RT and use the momentum direction
716 // there and the wire direction to determine the "shift".
717 // Note that whether the shift is to the left or to the right
718 // is not decided here. That ambiguity is left to be resolved later
719 // by applying a minus sign (or not) to the shift.
720 DVector3 pos_doca, mom_doca;
721 double s;
722 rt->DistToRT(wire, &s);
723 rt->GetLastDOCAPoint(pos_doca, mom_doca);
724
725 // The magnitude of the shift is based on the drift time. The
726 // value of the dist member of the DCDCTrackHit object does not
727 // subtract out the TOF. This can add 50-100 microns to the
728 // resolution in the CDC. Here, we actually can calculate the TOF
729 // (for a given mass hypothesis).
730 double mass = rt->GetMass();
731 double beta = 1.0/sqrt(1.0 + pow(mass/mom_doca.Mag(), 2.0))*2.998E10;
732 double tof = s/beta/1.0E-9; // in ns
733 hi.dist = 55E-4*(hit->time-tof);
734 hi.err = SIGMA_FDC_ANODE;
735
736 if(USE_FDC_CATHODE){
737 // Find whether the track is on the "left" or "right" of the wire
738 DVector3 shift = wire->udir.Cross(mom_doca);
739 if(shift.Mag()!=0.0)shift.SetMag(1.0);
740 double u = rt->GetLastDistAlongWire();
741 DVector3 pos_wire = wire->origin + u*wire->udir;
742 double LRsign = shift.Dot(pos_doca-pos_wire)<0.0 ? +1.0:-1.0;
743
744 // Lorentz corrected poisition along the wire is contained in x,y
745 // values, BUT, it is based on a left-right decision of the track
746 // segment. This may or may not be the same as the global track.
747 // As such, we have to determine the correction for our track.
748 //DVector3 wpos(hit->x, hit->y, wire->origin.Z());
749 //DVector3 wdiff = wpos - wire->origin;
750 //double u_corr = wire->udir.Dot(wdiff);
751 double alpha = mom_doca.Angle(DVector3(0,0,1));
752 hi.u_lorentz = LRsign*lorentz_def->GetLorentzCorrection(pos_doca.X(), pos_doca.Y(), pos_doca.Z(), alpha, hi.dist);
753 hi.u_dist = hit->s;
754 hi.u_err = SIGMA_FDC_CATHODE;
755 }else{
756 // User specified not to use FDC cathode information in the fit.
757 hi.u_dist = 0.0;
758 hi.u_err = 0.0; // setting u_err to zero means it's excluded from chi-sq
759 }
760 }
761 hinfo.push_back(hi);
762 }
763 } // USE_FDC
764}
765
766
767//------------------
768// GetResiInfo
769//------------------
770vector<bool> DTrackFitterALT1::GetResiInfo(DMatrix &state, const swim_step_t *start_step, DReferenceTrajectory *rt, hitsInfo &hinfo, vector<resiInfo> &residuals)
771{
772 /// Calculate the chi-squared for a track specified by state relative to the
773 /// given reference trajectory. This is just a wrapper for
774 /// ChiSq(DReferenceTrajectory *rt, vector<const DCoordinateSystem*> &wires, vector<DVector3> &shifts, vector<double> &errs, vector<double> &chisqv)
775 /// that accepts the state vector and re-swims the trajectory.
776
777 // In case we need to return early
778 int Nhits=0;
779 for(unsigned int i=0; i<hinfo.size(); i++){Nhits++; if(hinfo[i].u_err!=0.0)Nhits++;}
780 vector<bool> good_none(Nhits, false);
781
782 // "v" direction is perpendicular to both the rt direction and the
783 // x-direction. See LeastSquares() for more.
784 DVector3 vdir = start_step->sdir.Cross(start_step->mom);
785 if(vdir.Mag()!=0.0)vdir.SetMag(1.0);
786
787 DVector3 pos = start_step->origin
788 + state[state_x ][0]*start_step->sdir
789 + state[state_v ][0]*vdir;
790 DVector3 mom = state[state_px][0]*start_step->sdir
791 + state[state_py][0]*start_step->tdir
792 + state[state_pz][0]*start_step->udir;
793
794 if(!rt){
795 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<795<<
" "<<"NULL pointer passed for DReferenceTrajectory object!"<<endl;
796 return good_none;
797 }
798
799 if(pos.Mag()>200.0 || fabs(state[state_x ][0])>100.0 || fabs(state[state_v ][0])>100.0){
800 if(DEBUG_LEVEL>3)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<800<<
" "<<"state values out of range"<<endl;
801 if(DEBUG_LEVEL>6){
802 pos.Print();
803 mom.Print();
804 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<804<<
" "<<"state[state_x ][0]="<<state[state_x ][0]<<" state[state_v ][0]="<<state[state_x ][0]<<endl;
805 }
806 return good_none;
807 }
808
809 // Swim the trajectory with the specified state
810 rt->Swim(pos,mom);
811
812 // Sometimes, a bad state vector is passed that leads to referance trajectory with no steps.
813 if(rt->Nswim_steps<1){
814 residuals.clear();
815 return good_none;
816 }
817
818 return GetResiInfo(rt, hinfo, residuals);
819}
820
821//------------------
822// GetResiInfo
823//------------------
824vector<bool> DTrackFitterALT1::GetResiInfo(DReferenceTrajectory *rt, hitsInfo &hinfo, vector<resiInfo> &residuals)
825{
826 // Make a serialized list of "good" hits as we sparsely fill the residuals
827 // container with only the good ones.
828 vector<bool> good;
829
830 // Loop over wires hit. Make lists of finite residuals with layer numbers
831 // from which to build the actual matrices used to calculate chisq below.
832 residuals.clear();
833 for(unsigned int i=0; i<hinfo.size(); i++){
834 hitInfo &hi = hinfo[i];
835 hi.good = hi.good_u = false;
836
837 const DCoordinateSystem *wire = hi.wire;
838
839 // Figure out whether this is a CDC or FDC wire. Note that
840 // it could be neither if the target constraint is used.
841 const DFDCWire *fdcwire = dynamic_cast<const DFDCWire*>(wire);
842 const DCDCWire *cdcwire = dynamic_cast<const DCDCWire*>(wire);
843
844 // Get distance of the wire from the reference trajectory and the
845 // distance s along the track to the point of closest approach.
846 double s;
847 double d = rt->DistToRT(wire, &s);
848
849 // Residual. If we're on the correct side of the wire, then this is
850 // dist-doca. If we're on the incorrect side of the wire, then this
851 // is -dist-doca. Prior to calling us, the value of hi.dist will have
852 // a sign that has already been assigned to indicate the side of the wire
853 // the track is believed to be on.
854 double resi = hi.dist - d;
855 if(finite(resi)){
856 resiInfo ri;
857 ri.hit = &hi;
858 ri.layer = cdcwire ? cdcwire->ring:(fdcwire ? fdcwire->layer:0);
859 ri.resi_type = cdcwire ? resi_type_cdc_anode:(fdcwire ? resi_type_fdc_anode:resi_type_other);
860 ri.resi = resi;
861 ri.err = hi.err;
862 ri.step = rt->GetLastSwimStep();
863 hi.good = true;
864 residuals.push_back(ri);
865 good.push_back(true);
866 }else{
867 good.push_back(false);
868 }
869
870 // Also add residual along the wire. If the value of u_err is zero
871 // that indicates no measurement was made along the wire for this hit.
872 if(hi.u_err!=0.0){
873 // The sign of hi.dist indicates whether we want to treat this hit as
874 // being on the same side of the wire as the track(+) or the opposite
875 // side (-). In the latter case, we need to apply the Lorentz correction
876 // to the position along the wire in the opposite direction than we
877 // would otherwise. Set the sign of the Lorentz deflection based on the
878 // sign of hi.dist.
879 double LRsign = hi.dist<0.0 ? -1.0:1.0;
880
881 double u = rt->GetLastDistAlongWire();
882 double u_corrected = hi.u_dist + LRsign*hi.u_lorentz;
883 double resic = u - u_corrected;
884 if(finite(resic)){
885 resiInfo ri;
886 ri.hit = &hi;
887 ri.layer = fdcwire ? fdcwire->layer:0;
888 ri.resi_type =fdcwire ? resi_type_fdc_cathode:resi_type_other;
889 ri.resi = resic;
890 ri.err = hi.u_err;
891 ri.step = rt->GetLastSwimStep();
892 hi.good_u = true;
893 residuals.push_back(ri);
894 good.push_back(true);
895 }else{
896 good.push_back(false);
897 }
898 }
899 }
900
901 return good;
902}
903
904//------------------
905// LeastSquaresB
906//------------------
907DTrackFitterALT1::fit_status_t DTrackFitterALT1::LeastSquaresB(hitsInfo &hinfo, DReferenceTrajectory *rt)
908{
909 /// Fit the track with starting parameters given in the first step
910 /// of the reference trajectory rt. On return, the reference
911 /// trajectory rt will represent the final fit parameters and
912 /// chisq, Ndof, resiv, cov_meas, cov_muls, and cov_parm will be
913 /// filled based on the fit results.
914 ///
915 /// This determines the best fit of the track using the least squares method
916 /// described by R. Mankel Rep. Prog. Phys. 67 (2004) 553-622 pg 565.
917 /// Since it uses a linear approximation for the chisq dependance on
918 /// the fit parameters, several calls may be required for convergence.
919
920 // Initialize the chisq and Ndof data members in case we need to bail early
921 this->chisq = 1.0E6;
922 this->Ndof = 0;
923
924 // Make sure RT is not empty
925 if(rt->Nswim_steps<1)return kFitFailed;
926
927 // For fitting, we want to define a coordinate system very similar to the
928 // Reference Trajectory coordinate system. The difference is that we want
929 // the position to be defined in a plane perpendicular to the momentum.
930 // The RT x-direction is in this plane, but the momentum itself lies
931 // somewhere in the Y/Z plane so that neither Y nor Z makes a good choice
932 // for the second postion dimension. We will call the second coordinate in
933 // the perpendicular plane "v" which is the coordinate along the axis that
934 // is perpendicular to both the x-direction and the momentum direction.
935 swim_step_t &start_step = rt->swim_steps[0];
936 DVector3 pos = start_step.origin;
937 DVector3 mom = start_step.mom;
938
939 // Get the chi-squared vector for the initial reference trajectory
940 double initial_chisq;
941 int initial_Ndof;
942 vector<resiInfo> tmpresiduals;
943 vector<bool> good_initial = GetResiInfo(rt, hinfo, tmpresiduals);
944 double initial_chisq_per_dof = ChiSq(tmpresiduals, &initial_chisq, &initial_Ndof);
945 DMatrix resiv_initial(resiv);
946 DMatrix cov_meas_initial(cov_meas);
947 DMatrix cov_muls_initial(cov_muls);
948 DMatrix weights_initial(weights);
949
950 // Check that the initial chisq is reasonable before continuing
951 if(initial_Ndof<1){
952 if(DEBUG_LEVEL>0)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<952<<
" "<<"Initial Ndof<1. Aborting fit."<<endl;
953 return kFitFailed;
954 }
955
956 // Because we have a complicated non-linear system, we take the derivatives
957 // numerically.
958 //
959 // Note that in the calculations of the deltas below,
960 // the change in state should be set first and the value
961 // of deltas[...] calculated from that. See Numerical
962 // Recipes in C 2nd ed. section 5.7 ppg. 186-189.
963 const int Nparameters = 5;
964 double deltas[Nparameters];
965 DMatrix state(5,1);
966 switch(Nparameters){
967 case 5: state[state_v ][0] = 0.0;
968 case 4: state[state_x ][0] = 0.0;
969 case 3: state[state_pz ][0] = mom.Dot(start_step.udir);
970 case 2: state[state_py ][0] = mom.Dot(start_step.tdir);
971 case 1: state[state_px ][0] = mom.Dot(start_step.sdir);
972 }
973
974 // For the swimming below, we use a second reference trajectory so as
975 // to preserve the original. Set the charge here. The rest of the
976 // parameters (starting position and momentum) will be set using
977 // values from the state vector.
978 tmprt->q = rt->q;
979
980 // dpx : tweak by 0.0001
981 DMatrix state_dpx = state;
982 state_dpx[state_px][0] += LEAST_SQUARES_DP;
983 deltas[state_px] = state_dpx[state_px][0] - state[state_px][0];
984 vector<bool> good_px = GetResiInfo(state_dpx, &start_step, tmprt, hinfo, tmpresiduals);
985 double chisq_dpx = ChiSq(tmpresiduals);
986 DMatrix resiv_dpx_hi(resiv);
987 DMatrix &resiv_dpx_lo = resiv_initial;
988
989 // dpy : tweak by 0.0001
990 DMatrix state_dpy = state;
991 state_dpy[state_py][0] += LEAST_SQUARES_DP;
992 deltas[state_py] = state_dpy[state_py][0] - state[state_py][0];
993 vector<bool> good_py = GetResiInfo(state_dpy, &start_step, tmprt, hinfo, tmpresiduals);
994 double chisq_dpy = ChiSq(tmpresiduals);
995 DMatrix resiv_dpy_hi(resiv);
996 DMatrix &resiv_dpy_lo = resiv_initial;
997
998 // dpz : tweak by 0.0001
999 DMatrix state_dpz = state;
1000 state_dpz[state_pz][0] += LEAST_SQUARES_DP;
1001 deltas[state_pz] = state_dpz[state_pz][0] - state[state_pz][0];
1002 vector<bool> good_pz = GetResiInfo(state_dpz, &start_step, tmprt, hinfo, tmpresiduals);
1003 double chisq_dpz = ChiSq(tmpresiduals);
1004 DMatrix resiv_dpz_hi(resiv);
1005 DMatrix &resiv_dpz_lo = resiv_initial;
1006
1007 // dv : tweak by 0.01
1008 DMatrix state_dv = state;
1009 state_dv[state_v][0] += LEAST_SQUARES_DX;
1010 deltas[state_v] = state_dv[state_v][0] - state[state_v][0];
1011 vector<bool> good_v = GetResiInfo(state_dv, &start_step, tmprt, hinfo, tmpresiduals);
1012 double chisq_dv = ChiSq(tmpresiduals);
1013 DMatrix resiv_dv_hi(resiv);
1014 DMatrix &resiv_dv_lo = resiv_initial;
1015
1016 // dx : tweak by 0.01
1017 DMatrix state_dx = state;
1018 state_dx[state_x][0] += LEAST_SQUARES_DX;
1019 deltas[state_x] = state_dx[state_x][0] - state[state_x][0];
1020 vector<bool> good_x = GetResiInfo(state_dx, &start_step, tmprt, hinfo, tmpresiduals);
1021 double chisq_dx = ChiSq(tmpresiduals);
1022 DMatrix resiv_dx_hi(resiv);
1023 DMatrix &resiv_dx_lo = resiv_initial;
1024
1025 // Verify all of the "good" vectors are of the same length
1026 unsigned int size_good = good_initial.size();
1027 if( (good_px.size() != size_good)
1028 || (good_py.size() != size_good)
1029 || (good_pz.size() != size_good)
1030 || (good_v.size() != size_good)
1031 || (good_x.size() != size_good)){
1032 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1032
<<" "<<"Size of \"good\" vectors don't match! size_good="<<size_good<<endl;
1033 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1033
<<" "<<"good_px.size()="<<good_px.size()<<endl;
1034 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1034
<<" "<<"good_py.size()="<<good_py.size()<<endl;
1035 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1035
<<" "<<"good_pz.size()="<<good_pz.size()<<endl;
1036 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1036
<<" "<<" good_v.size()="<<good_v.size()<<endl;
1037 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1037
<<" "<<" good_x.size()="<<good_x.size()<<endl;
1038 return kFitFailed;
1039 }
1040
1041 // We need to get a list of hits that are good for all of the tweaked
1042 // tracks as well as the initial track.
1043 unsigned int Ngood = 0;
1044 vector<bool> good_all;
1045 for(unsigned int i=0; i<size_good; i++){
1046 bool isgood = true;
1047 isgood &= good_initial[i];
1048 isgood &= good_px[i];
1049 isgood &= good_py[i];
1050 isgood &= good_pz[i];
1051 isgood &= good_v[i];
1052 isgood &= good_x[i];
1053 good_all.push_back(isgood);
1054 if(isgood)Ngood++;
1055 }
1056
1057 // Make sure there is a minimum number of hits before attempting a fit
1058 if(Ngood<LEAST_SQUARES_MIN_HITS){
1059 if(DEBUG_LEVEL>0)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1059
<<" "<<"Not enough good hits to do fit. Aborting..."<<endl;
1060 return kFitFailed;
1061 }
1062
1063 // Use the good_all map to filter out hits for each residual vector
1064 // for which somebody else did not have a good hit.
1065 FilterGood(resiv_initial, good_initial, good_all);
1066 FilterGood(resiv_dpx_hi , good_px, good_all);
1067 FilterGood(resiv_dpy_hi , good_py, good_all);
1068 FilterGood(resiv_dpz_hi , good_pz, good_all);
1069 FilterGood(resiv_dv_hi , good_v , good_all);
1070 FilterGood(resiv_dx_hi , good_x , good_all);
1071 FilterGood(weights_initial , good_x , good_all);
1072
1073 // Here, we check that the the residual vectors are of the same
1074 // dimension for all tweaked tracks and the initial track so that we
1075 // can proceed with building the fit matrices below.
1076 int Nhits = resiv_initial.GetNrows();
1077 if( (resiv_dpx_hi.GetNrows() != Nhits)
1078 || (resiv_dpy_hi.GetNrows() != Nhits)
1079 || (resiv_dpz_hi.GetNrows() != Nhits)
1080 || (resiv_dx_hi.GetNrows() != Nhits)
1081 || (resiv_dv_hi.GetNrows() != Nhits)){
1082 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1082
<<" "<<"Size of residual vectors don't match! Nhits="<<Nhits<<endl;
1083 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1083
<<" "<<"resiv_dpx_hi.GetNrows()="<<resiv_dpx_hi.GetNrows()<<endl;
1084 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1084
<<" "<<"resiv_dpy_hi.GetNrows()="<<resiv_dpy_hi.GetNrows()<<endl;
1085 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1085
<<" "<<"resiv_dpz_hi.GetNrows()="<<resiv_dpz_hi.GetNrows()<<endl;
1086 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1086
<<" "<<" resiv_dx_hi.GetNrows()="<<resiv_dx_hi.GetNrows()<<endl;
1087 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1087
<<" "<<" resiv_dv_hi.GetNrows()="<<resiv_dv_hi.GetNrows()<<endl;
1088 return kFitFailed;
1089 }
1090
1091 // Print some debug messages
1092 if(DEBUG_LEVEL>4){
1093 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1093
<<" "<<"initial_chisq_per_dof="<<initial_chisq_per_dof<<endl;
1094 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1094
<<" "<<"chisq_dpx="<<chisq_dpx<<endl;
1095 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1095
<<" "<<"chisq_dpy="<<chisq_dpy<<endl;
1096 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1096
<<" "<<"chisq_dpz="<<chisq_dpz<<endl;
1097 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1097
<<" "<<"chisq_dv="<<chisq_dv<<endl;
1098 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1098
<<" "<<"chisq_dx="<<chisq_dx<<endl;
1099 if(DEBUG_LEVEL>10){
1100 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1100
<<" "<<"hit\tinitial\tpx \tpy \tpz \tx \tv"<<endl;
1101 for(int j=0; j<resiv_initial.GetNrows(); j++){
1102 cout<<j<<"\t";
1103 cout<<resiv_initial[j][0]/sqrt(cov_meas[j][j])<<"\t";
1104 cout<<resiv_dpx_hi[j][0]/sqrt(cov_meas[j][j])<<"\t";
1105 cout<<resiv_dpy_hi[j][0]/sqrt(cov_meas[j][j])<<"\t";
1106 cout<<resiv_dpz_hi[j][0]/sqrt(cov_meas[j][j])<<"\t";
1107 cout<<resiv_dv_hi[j][0]/sqrt(cov_meas[j][j])<<"\t";
1108 cout<<resiv_dx_hi[j][0]/sqrt(cov_meas[j][j])<<endl;
1109 }
1110 }
1111 }
1112
1113 // Build "F" matrix of derivatives
1114 DMatrix F(Nhits,Nparameters);
1115 for(int i=0; i<Nhits; i++){
1116 switch(Nparameters){
1117 // Note: This is a funny way to use a switch!
1118 case 5: F[i][state_v ] = (resiv_dv_hi[i][0]-resiv_dv_lo[i][0])/deltas[state_v];
1119 case 4: F[i][state_x ] = (resiv_dx_hi[i][0]-resiv_dx_lo[i][0])/deltas[state_x];
1120 case 3: F[i][state_pz] = (resiv_dpz_hi[i][0]-resiv_dpz_lo[i][0])/deltas[state_pz];
1121 case 2: F[i][state_py] = (resiv_dpy_hi[i][0]-resiv_dpy_lo[i][0])/deltas[state_py];
1122 case 1: F[i][state_px] = (resiv_dpx_hi[i][0]-resiv_dpx_lo[i][0])/deltas[state_px];
1123 }
1124 }
1125
1126 // Sometimes, "F" has lots of values like 1.44E+09 indicating a problem (I think comming
1127 // from some nan values floating around.) Anyway, in these cases, the E2Norm value is
1128 // quite large (>1E+18) which we can use to punt now. In reality, we do this to avoid
1129 // ROOT error messages about a matrix being singular when Ft*Vinv*F is inverted below.
1130 if(F.E2Norm()>1.0E18){
1131 if(DEBUG_LEVEL>1){
1132 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1132
<<" "<<" -- F matrix E2Norm out of range(E2Norm="<<F.E2Norm()<<" max="<<1.0E18<<")"<<endl;
1133 }
1134 return kFitFailed;
1135 }
1136
1137 // Transpose of "F" matrix
1138 DMatrix Ft(DMatrix::kTransposed, F);
1139
1140 // Calculate the "B" matrix using the weights from the initial chisq
1141 DMatrix &Vinv = weights_initial; // Stick with Mankel naming convention
1142 DMatrix B(DMatrix::kInverted, Ft*Vinv*F);
1143
1144 // If the inversion failed altogether then the invalid flag
1145 // will be set on the matrix. In these cases, we're dead.
1146 if(!B.IsValid()){
1147 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1147
<<" "<<" -- B matrix invalid"<<endl;
1148 return kFitFailed;
1149 }
1150
1151 // The "B" matrix happens to be the covariance matrix of the
1152 // state parameters. A problem sometimes occurs where one or
1153 // more elements of B are very large. This tends to happen
1154 // when a column of F is essentially all zeros making the
1155 // matrix un-invertable. What we should really do in these
1156 // cases is check beforehand and drop the bad column(s)
1157 // before trying to invert. That will add complication that
1158 // I don't want to introduce quite yet. What we do now
1159 // is check for it and punt rather than return a nonsensical
1160 // value.
1161 if(B.E2Norm() > LEAST_SQUARES_MAX_E2NORM){
1162 if(DEBUG_LEVEL>1){
1163 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1163
<<" "<<" -- B matrix E2Norm out of range(E2Norm="<<B.E2Norm()<<" max="<<LEAST_SQUARES_MAX_E2NORM<<")"<<endl;
1164 cout<<"--- B ---"<<endl;
1165 B.Print();
1166 }
1167 return kFitFailed;
1168 }
1169
1170 // Copy the B matrix into cov_parm to later copy into DTrack
1171 cov_parm.ResizeTo(B);
1172 cov_parm = B;
1173
1174 // Calculate step direction and magnitude
1175 DMatrix delta_state = B*Ft*Vinv*resiv_initial;
1176
1177 // The following is based on Numerical Recipes in C 2nd Ed.
1178 // ppg. 384-385.
1179 //
1180 // We now have the "direction" by which to step in the-d
1181 // parameter space as well as an amplitude by which to
1182 // step in "delta_state". A potential problem is that
1183 // we can over-step such that we end up at a worse
1184 // chi-squared value. To address this, we try taking
1185 // the full step, but if we end up at a worse chi-sq
1186 // then we backtrack and take a smaller step until
1187 // we see the chi-sq improve.
1188 double min_chisq_per_dof = initial_chisq_per_dof;
1189 double min_lambda = 0.0;
1190 double lambda = -1.0;
1191 int Ntrys = 0;
1192 int max_trys = 6;
1193 DMatrix new_state(5,1);
1194 for(; Ntrys<max_trys; Ntrys++){
1195
1196 // Scale the delta by lambda to take a partial step (except the 1st iteration where lambda is 1)
1197 for(int i=0; i<Nparameters; i++)new_state[i][0] = state[i][0] + delta_state[i][0]*lambda;
1198
1199 GetResiInfo(new_state, &start_step, tmprt, hinfo, tmpresiduals);
1200 double chisq_per_dof = ChiSq(tmpresiduals);
1201
1202 if(chisq_per_dof<min_chisq_per_dof){
1203 min_chisq_per_dof = chisq_per_dof;
1204 min_lambda = lambda;
1205 }
1206
1207 // If we're at a lower chi-sq then we're done
1208 if(DEBUG_LEVEL>4)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1208
<<" "<<" -- initial_chisq_per_dof="<<initial_chisq_per_dof<<" new chisq_per_dof="<<chisq_per_dof<<" nhits="<<resiv.GetNrows()<<" p="<<tmprt->swim_steps[0].mom.Mag()<<" lambda="<<lambda<<endl;
1209 if(chisq_per_dof-initial_chisq_per_dof < 0.1 && chisq_per_dof<2.0)break;
1210
1211 if(chisq_per_dof<initial_chisq_per_dof)break;
1212
1213 // Chi-sq was increased, try a smaller step on the next iteration
1214 lambda/=2.0;
1215 }
1216
1217 // If we failed to find a better Chi-Sq above, maybe we were looking
1218 // in the wrong direction(??) Try looking in the opposite direction.
1219 //if(Ntrys>=max_trys && (min_chisq_per_dof>=initial_chisq_per_dof || min_chisq_per_dof>1.0)){
1220 if(Ntrys>=max_trys){
1221 lambda = 1.0/4.0;
1222 for(int j=0; j<3; j++, Ntrys++){
1223
1224 // Scale the delta by lambda to take a partial step (except the 1st iteration where lambda is 1)
1225 for(int i=0; i<Nparameters; i++)new_state[i][0] = state[i][0] + delta_state[i][0]*lambda;
1226
1227 GetResiInfo(new_state, &start_step, tmprt, hinfo, tmpresiduals);
1228 double chisq_per_dof = ChiSq(tmpresiduals);
1229
1230 if(chisq_per_dof<min_chisq_per_dof){
1231 min_chisq_per_dof = chisq_per_dof;
1232 min_lambda = lambda;
1233 }
1234
1235 // If we're at a lower chi-sq then we're done
1236 if(DEBUG_LEVEL>4)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1236
<<" "<<" -- initial_chisq_per_dof="<<initial_chisq_per_dof<<" new chisq_per_dof="<<chisq_per_dof<<" nhits="<<resiv.GetNrows()<<" p="<<tmprt->swim_steps[0].mom.Mag()<<" lambda="<<lambda<<endl;
1237 if(chisq_per_dof-initial_chisq_per_dof < 0.1 && chisq_per_dof<2.0)break;
1238
1239 if(chisq_per_dof<initial_chisq_per_dof)break;
1240
1241 // Chi-sq was increased, try a smaller step on the next iteration
1242 lambda/=2.0;
1243 }
1244 }
1245
1246 // If we failed to make a step to a smaller chi-sq then signal
1247 // that we were unable to make any improvement.
1248 if(min_lambda==0.0){
1249 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1249
<<" "<<"Chisq only increased (both directions searched!)"<<endl;
1250 GetResiInfo(rt, hinfo, tmpresiduals);
1251 ChiSq(tmpresiduals, &this->chisq, &this->Ndof); // refill resiv, cov_meas, ...
1252 return kFitNoImprovement;
1253 }
1254
1255 // Re-create new_state using min_lambda
1256 for(int i=0; i<Nparameters; i++)new_state[i][0] = state[i][0] + delta_state[i][0]*min_lambda;
1257 if(DEBUG_HISTS)this->lambda->Fill(min_lambda);
1258
1259 // Re-swim reference trajectory using these parameters and re-calc chisq.
1260 // Note that here we have the chisq and Ndof members set.
1261 GetResiInfo(new_state, &start_step, rt, hinfo, tmpresiduals);
1262 ChiSq(tmpresiduals, &this->chisq, &this->Ndof);
1263
1264 if(DEBUG_LEVEL>3){
1265 DVector3 pos = start_step.origin;
1266 DVector3 mom = start_step.mom;
1267 double phi = mom.Phi();
1268 if(phi<0.0)phi+=2.0*M_PI3.14159265358979323846;
1269 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1269
<<" "<<"LeastSquaresB succeeded: p="<<mom.Mag()<<" theta="<<mom.Theta()<<" phi="<<phi<<" z="<<pos.Z()<<endl;
1270 }
1271
1272 return kFitSuccess;
1273}
1274
1275//------------------
1276// FilterGood
1277//------------------
1278void DTrackFitterALT1::FilterGood(DMatrix &my_resiv, vector<bool> &my_good, vector<bool> &good_all)
1279{
1280 /// Remove elements from my_resiv that are not "good" according to the good_all list.
1281 /// The my_good and good_all vectors should be the same size. The number of "true"
1282 /// entries in my_good should be the size of the my_resiv vector. For entries in the
1283 /// my_good vector that are true, but have an entry in good_all that is false, the
1284 /// corresponding entry in my_resiv will be removed. The my_resiv matrix (which should be
1285 /// N x 1) will be resized upon exit. The my_good vector will be set equal to good_all
1286 /// upon exit also so that my_resiv and my_good stay in sync.
1287
1288 // Make list of rows (and columns) we should keep
1289 vector<int> rows_to_keep;
1290 for(unsigned int i=0, n=0; i<my_good.size(); i++){
1
Loop condition is true. Entering loop body
→
3
←
Loop condition is true. Entering loop body
→
5
←
Loop condition is true. Entering loop body
→
8
←
Loop condition is false. Execution continues on line 1296
→
1291 if(my_good[i] && good_all[i])rows_to_keep.push_back(n);
6
←
Taking true branch
→
1292 if(my_good[i])n++;
2
←
Taking false branch
→
4
←
Taking false branch
→
7
←
Taking false branch
→
1293 }
1294
1295 // Copy my_resiv to a temporary matrix and resize my_resiv to the new size
1296 int Nrows = (int)rows_to_keep.size();
1297 int Ncols = my_resiv.GetNcols()>1 ? Nrows:1;
9
←
'?' condition is false
→
1298 DMatrix tmp(my_resiv);
1299 my_resiv.ResizeTo(Nrows, Ncols);
1300
1301 // Loop over rows and columns copying in the elements we're keeping
1302 for(int i=0; i<Nrows; i++){
10
←
Loop condition is true. Entering loop body
→
1303 int irow = rows_to_keep[i];
11
←
Dereference of null pointer
1304 for(int j=0; j<Ncols; j++){
1305 int icol = Ncols>1 ? rows_to_keep[j]:0;
1306 my_resiv[i][j] = tmp[irow][icol];
1307 }
1308 }
1309
1310 my_good = good_all;
1311}
1312
1313//------------------
1314// PrintChisqElements
1315//------------------
1316void DTrackFitterALT1::PrintChisqElements(DMatrix &resiv, DMatrix &cov_meas, DMatrix &cov_muls, DMatrix &weights)
1317{
1318 /// This is for debugging only.
1319 int Nhits = resiv.GetNrows();
1320 double chisq_diagonal = 0.0;
1321 for(int i=0; i<Nhits; i++){
1322 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1322
<<" "<<" r/sigma "<<i<<": "<<resiv[i][0]*sqrt(weights[i][i])
1323 <<" resi="<<resiv[i][0]
1324 <<" sigma="<<1.0/sqrt(weights[i][i])
1325 <<" cov_meas="<<cov_meas[i][i]
1326 <<" cov_muls="<<cov_muls[i][i]
1327 <<endl;
1328 chisq_diagonal += pow(resiv[i][0], 2.0)*weights[i][i];
1329 }
1330 DMatrix resiv_t(DMatrix::kTransposed, resiv);
1331 DMatrix chisqM(resiv_t*weights*resiv);
1332 int Ndof = Nhits - 5; // assume 5 fit parameters
1333
1334 _DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1334
<<" "<<" chisq/Ndof: "<<chisqM[0][0]/(double)Ndof<<" chisq/Ndof diagonal elements only:"<<chisq_diagonal/(double)Ndof<<endl;
1335}
1336
1337//------------------
1338// ForceLRTruth
1339//------------------
1340void DTrackFitterALT1::ForceLRTruth(JEventLoop *loop, DReferenceTrajectory *rt, hitsInfo &hinfo)
1341{
1342 /// This routine is called when the TRKFIT:LR_FORCE_TRUTH parameters is
1343 /// set to a non-zero value (e.g. -PTRKFIT:LR_FORCE_TRUTH=1 is passed
1344 /// on the command line). This is used only for debugging and only with
1345 /// Monte Carlo data.
1346 ///
1347 /// This routine will adjust the left-right choice of each hit based
1348 /// on the current fit track (represented by rt), and the truth information
1349 /// contained in the the DMCTruthPoint objects. It assumes the hits in
1350 /// hinfo correspond to the track represented by rt.
1351
1352 // Get Truth hits
1353 vector<const DMCTrackHit*> mctrackhits;
1354 loop->Get(mctrackhits);
1355
1356 // Loop over hits
1357 for(unsigned int i=0; i<hinfo.size(); i++){
1358 hitInfo &hi = hinfo[i];
1359 const DCoordinateSystem *wire = hi.wire;
1360 if(wire==target)continue; // ignore target
1361
1362 // Sometimes dist is NaN
1363 if(!finite(hi.dist)){
1364 hi.err = 100.0;
1365 hi.u_err = 0.0;
1366 continue;
1367 }
1368
1369 // Find the truth hit corresponding to this real hit
1370 const DMCTrackHit *mctrackhit = DTrackHitSelectorTHROWN::GetMCTrackHit(hi.wire, fabs(hi.dist), mctrackhits, 0);
1371
1372 // If no truth hit was found, then this may be a noise hit. Set
1373 // the error to something large so it is ignored and warn the
1374 // user.
1375 if(!mctrackhit){
1376 if(DEBUG_LEVEL>1)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1376
<<" "<<"No DMCTrackHit found corresponding to hit "<<i<<" in hinfo! (noise hit?)"<<endl;
1377 hi.err = 100.0;
1378 hi.u_err = 0.0;
1379 continue;
1380 }
1381
1382 // We do this by looking at the direction of the vector pointing from the
1383 // DOCA point on the wire to that of the fit track and comparing it to
1384 // a similar vector pointing to the truth point. If they both are generally
1385 // in the same direction (small angle) then the fit track is considered
1386 // as being on the correct side of the wire. If they have an angle somewhere
1387 // in the 180 degree range, the fit is considered to be on the wrong side
1388 // of the wire and it is "flipped".
1389 //
1390 // It can also be that the truth vector and the fit vector are at nearly
1391 // a right angle with respect to one another. This really only happens in
1392 // the CDC for tracks going in roughly the same direction as the wire.
1393 // In these cases, the truth point can't help us solve the LR so we set
1394 // the drift distance to zero and give it a larger error so this hit will
1395 // still be included, but appropriately weighted.
1396
1397 // Vector pointing from wire to the fit track
1398 DVector3 pos_doca, mom_doca;
1399 rt->DistToRT(wire);
1400 rt->GetLastDOCAPoint(pos_doca, mom_doca);
1401 DVector3 shift = wire->udir.Cross(mom_doca);
1402 double u = rt->GetLastDistAlongWire();
1403 DVector3 pos_wire = wire->origin + u*wire->udir;
1404
1405 // Vector pointing from wire to the truth point
1406 DVector3 pos_truth(mctrackhit->r*cos(mctrackhit->phi), mctrackhit->r*sin(mctrackhit->phi), mctrackhit->z);
1407 DVector3 pos_wire_truth = wire->origin + (pos_truth - wire->origin).Dot(wire->udir)*wire->udir;
1408
1409 if(DEBUG_LEVEL>8)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1409
<<" "<<" "<<i+1<<": (pos_doca-pos_wire).Angle(pos_truth-pos_wire_truth)="<<(pos_doca-pos_wire).Angle(pos_truth-pos_wire_truth)*57.3<<endl;
1410
1411 // Decide what to do based on angle between fit and truth vectors
1412 double angle_fit_truth = (pos_doca-pos_wire).Angle(pos_truth-pos_wire_truth);
1413 if(fabs( fabs(angle_fit_truth)-M_PI3.14159265358979323846/2.0) < M_PI3.14159265358979323846/3.0){
1414 // Fit and truth vector are nearly at 90 degrees. Fit this hit
1415 // only to wire position
1416 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1416
<<" "<<"Downgrading "<<i+1<<"th hit to wire-based (hi.u_err was:"<<hi.u_err<<")"<<endl;
1417 hi.dist = 0.0;
1418 hi.err = 0.35;
1419 hi.u_err = 0.0;
1420 }else{
1421 // Fit and truth vector are nearly (anti)parallel. Decide whether to "flip" hit
1422 if(fabs(angle_fit_truth) > M_PI3.14159265358979323846/2.0){
1423 if(DEBUG_LEVEL>5)_DBG_std::cerr<<"DTrackFitterALT1.cc"<<":"<<1423
<<" "<<"Flipping side "<<i+1<<"th hit "<<endl;
1424 hi.dist = -hi.dist;
1425 }
1426 }
1427 }
1428}
1429
1430//------------------
1431// FillDebugHists
1432//------------------
1433void DTrackFitterALT1::FillDebugHists(DReferenceTrajectory *rt, DVector3 &vertex_pos, DVector3 &vertex_mom)
1434{
1435 //vertex_mom.SetMagThetaPhi(6.0, 17.2*M_PI/180.0, 90.0*M_PI/180.0);
1436 //vertex_pos.SetXYZ(0.0,0.0,65.0);
1437 //rt->Swim(vertex_pos, vertex_mom);
1438 ptotal->Fill(vertex_mom.Mag());
1439
1440 // Calculate particle beta
1441 double beta = 1.0/sqrt(1.0+pow(rt->GetMass(), 2.0)/vertex_mom.Mag2()); // assume this is a pion for now. This should eventually come from outer detectors
1442
1443 for(unsigned int j=0; j<cdchits.size(); j++){
1444 const DCDCTrackHit *hit = cdchits[j];
1445 const DCDCWire *wire = hit->wire;
1446
1447 // Distance of closest approach for track to wire
1448 double s;
1449 double doca = rt->DistToRT(wire, &s);
1450
1451 // Distance from drift time. Hardwired for simulation for now
1452 double tof = s/(beta*3E10*1E-9);
1453 double dist = (hit->tdrift-tof)*55E-4;
1454
1455 // NOTE: Sometimes this could be nan
1456 double resi = dist - doca;
1457 if(!finite(resi))continue;
1458
1459 // Fill histos
1460 residuals_cdc->Fill(resi, wire->ring);
1461 residuals_cdc_vs_s->Fill(resi, wire->ring, s);
1462
1463 cdcdoca_vs_dist->Fill(dist, doca);
1464 cdcdoca_vs_dist_vs_ring->Fill(dist, doca, wire->ring);
1465 if(wire->stereo==0.0){
1466 dist_axial->Fill(dist);
1467 doca_axial->Fill(doca);
1468 }else{
1469 dist_stereo->Fill(dist);
1470 doca_stereo->Fill(doca);
1471 }
1472
1473 }
1474
1475 for(unsigned int j=0; j<fdchits.size(); j++){
1476 const DFDCPseudo *hit = fdchits[j];
1477 const DFDCWire *wire = hit->wire;
1478
1479 // Distance of closest approach for track to wire
1480 double s;
1481 double doca = rt->DistToRT(wire,&s);
1482
1483 double tof = s/(beta*3E10*1E-9);
1484 double dist = (hit->time-tof)*55E-4;
1485
1486 // NOTE: Sometimes this is nan
1487 double resi = dist - doca;
1488 if(finite(resi)){
1489 fdcdoca_vs_dist->Fill(dist, doca);
1490 residuals_fdc_anode->Fill(resi, wire->layer);
1491 residuals_fdc_anode_vs_s->Fill(resi, wire->layer,s);
1492 }
1493
1494 double u = rt->GetLastDistAlongWire();
1495 resi = u - hit->s;
1496 if(finite(resi)){
1497 fdcu_vs_s->Fill(u, hit->s);
1498 residuals_fdc_cathode->Fill(resi, wire->layer);
1499 residuals_fdc_cathode_vs_s->Fill(resi, wire->layer,s);
1500 }
1501 }
1502}
1503
1504