Bug Summary

File:/volatile/halld/gluex/nightly/2024-03-23/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/TRACKING/DTrackHitSelectorALT2.cc
Location:line 819, column 5
Description:Value stored to 'var_phi' is never read

Annotated Source Code

1// $Id$
2//
3// File: DTrackHitSelectorALT2.cc
4// Created: Fri Feb 6 08:22:58 EST 2009
5// Creator: davidl (on Darwin harriet.jlab.org 9.6.0 i386)
6//
7
8#include <cmath>
9using namespace std;
10
11#include <TROOT.h>
12
13#include <TRACKING/DReferenceTrajectory.h>
14
15#include "DTrackHitSelectorALT2.h"
16
17#ifndef ansi_escape((char)0x1b)
18#define ansi_escape((char)0x1b) ((char)0x1b)
19#define ansi_bold((char)0x1b)<<"[1m" ansi_escape((char)0x1b)<<"[1m"
20#define ansi_normal((char)0x1b)<<"[0m" ansi_escape((char)0x1b)<<"[0m"
21#define ansi_red((char)0x1b)<<"[31m" ansi_escape((char)0x1b)<<"[31m"
22#define ansi_green((char)0x1b)<<"[32m" ansi_escape((char)0x1b)<<"[32m"
23#define ansi_blue((char)0x1b)<<"[34m" ansi_escape((char)0x1b)<<"[34m"
24#endif // ansi_escape
25
26#define ONE_OVER_SQRT120.288675 0.288675
27#define ONE_OVER_120.08333333333333 0.08333333333333
28#define EPS1e-6 1e-6
29
30bool static DTrackHitSelector_cdchit_cmp(pair<double,const DCDCTrackHit *>a,
31 pair<double,const DCDCTrackHit *>b){
32 if (a.second->wire->ring!=b.second->wire->ring)
33 return (a.second->wire->ring>b.second->wire->ring);
34 return (a.first>b.first);
35}
36bool static DTrackHitSelector_fdchit_cmp(pair<double,const DFDCPseudo *>a,
37 pair<double,const DFDCPseudo *>b){
38 if (a.second->wire->layer!=b.second->wire->layer)
39 return (a.second->wire->layer>b.second->wire->layer);
40 return (a.first>b.first);
41}
42bool static DTrackHitSelector_gemhit_cmp(pair<double,const DGEMPoint *>a,
43 pair<double,const DGEMPoint *>b){
44 if (a.second->detector!=b.second->detector)
45 return (a.second->detector>b.second->detector);
46 return (a.first>b.first);
47}
48bool static DTrackHitSelector_trdhit_cmp(pair<double,const DTRDPoint *>a,
49 pair<double,const DTRDPoint *>b){
50 if (a.second->detector!=b.second->detector)
51 return (a.second->detector>b.second->detector);
52 return (a.first>b.first);
53}
54
55bool static DTrackHitSelector_cdchit_in_cmp(const DCDCTrackHit *a, const DCDCTrackHit *b){
56 if (a->wire->ring != b->wire->ring) return (a->wire->ring < b->wire->ring);
57 return (a->wire->straw < b->wire->straw);
58}
59
60bool static DTrackHitSelector_fdchit_in_cmp(const DFDCPseudo *a, const DFDCPseudo *b){
61 if (a->wire->layer != b->wire->layer) return (a->wire->layer < b->wire->layer);
62 return (a->wire->wire < b->wire->wire);
63}
64
65
66//---------------------------------
67// DTrackHitSelectorALT2 (Constructor)
68//---------------------------------
69DTrackHitSelectorALT2::DTrackHitSelectorALT2(jana::JEventLoop *loop, int32_t runnumber):DTrackHitSelector(loop)
70{
71 HS_DEBUG_LEVEL = 0;
72 MAKE_DEBUG_TREES = false;
73 MIN_HIT_PROB_CDC = 0.01;
74 MIN_HIT_PROB_FDC = 0.01;
75 MIN_HIT_PROB_TRD = 0.01;
76 MIN_HIT_PROB_GEM = 0.01;
77 MIN_FDC_SIGMA_ANODE_CANDIDATE = 0.1000;
78 MIN_FDC_SIGMA_CATHODE_CANDIDATE = 0.1000;
79 MIN_FDC_SIGMA_ANODE_WIREBASED = 0.0100;
80 MIN_FDC_SIGMA_CATHODE_WIREBASED = 0.0100;
81 MAX_DOCA=2.5;
82
83 gPARMS->SetDefaultParameter("TRKFIT:MAX_DOCA",MAX_DOCA,"Maximum doca for associating hit with track");
84
85 gPARMS->SetDefaultParameter("TRKFIT:HS_DEBUG_LEVEL", HS_DEBUG_LEVEL, "Debug verbosity level for hit selector used in track fitting (0=no debug messages)");
86 gPARMS->SetDefaultParameter("TRKFIT:MAKE_DEBUG_TREES", MAKE_DEBUG_TREES, "Create a TTree with debugging info on hit selection for the FDC and CDC");
87 gPARMS->SetDefaultParameter("TRKFIT:MIN_HIT_PROB_CDC", MIN_HIT_PROB_CDC, "Minimum probability a CDC hit may have to be associated with a track to be included in list passed to fitter");
88 gPARMS->SetDefaultParameter("TRKFIT:MIN_HIT_PROB_FDC", MIN_HIT_PROB_FDC, "Minimum probability a FDC hit may have to be associated with a track to be included in list passed to fitter");
89 gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_ANODE_CANDIDATE", MIN_FDC_SIGMA_ANODE_CANDIDATE, "Minimum sigma used for FDC anode hits on track candidates");
90 gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_CATHODE_CANDIDATE", MIN_FDC_SIGMA_CATHODE_CANDIDATE, "Minimum sigma used for FDC cathode hits on track candidates");
91 gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_ANODE_WIREBASED", MIN_FDC_SIGMA_ANODE_WIREBASED, "Minimum sigma used for FDC anode hits on wire-based tracks");
92 gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_CATHODE_WIREBASED", MIN_FDC_SIGMA_CATHODE_WIREBASED, "Minimum sigma used for FDC cathode hits on wire-based tracks");
93
94 cdchitsel = NULL__null;
95 fdchitsel = NULL__null;
96 if(MAKE_DEBUG_TREES){
97 loop->GetJApplication()->Lock();
98
99 cdchitsel= (TTree*)gROOT(ROOT::GetROOT())->FindObject("cdchitsel");
100 if(!cdchitsel){
101 cdchitsel = new TTree("cdchitsel", "CDC Hit Selector");
102 cdchitsel->Branch("H", &cdchitdbg, "fit_type/I:p/F:theta:mass:sigma:x:y:z:s:itheta02:itheta02s:itheta02s2:dist:doca:resi:chisq:prob:sig_phi:sig_lambda:sig_pt");
103 }else{
104 _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<104<<std::endl;
105 jerr<<" !!! WARNING !!!"<<endl;
106 jerr<<"It appears that the cdchitsel TTree is already defined."<<endl;
107 jerr<<"This is probably means you are running with multiple threads."<<endl;
108 jerr<<"To avoid complication, filling of the hit selector debug"<<endl;
109 jerr<<"trees will be disabled for this thread."<<endl;
110 _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<110<<std::endl;
111 cdchitsel = NULL__null;
112 }
113
114 fdchitsel= (TTree*)gROOT(ROOT::GetROOT())->FindObject("fdchitsel");
115 if(!fdchitsel){
116 fdchitsel = new TTree("fdchitsel", "FDC Hit Selector");
117 fdchitsel->Branch("H", &fdchitdbg, "fit_type/I:p/F:theta:mass:sigma_anode:sigma_cathode:x:y:z:s:itheta02:itheta02s:itheta02s2:dist:doca:resi:u:u_cathodes:resic:chisq:prob:sig_phi:sig_lambda:sig_pt");
118 }else{
119 _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<119<<std::endl;
120 jerr<<" !!! WARNING !!!"<<endl;
121 jerr<<"It appears that the fdchitsel TTree is already defined."<<endl;
122 jerr<<"This is probably means you are running with multiple threads."<<endl;
123 jerr<<"To avoid complication, filling of the hit selector debug"<<endl;
124 jerr<<"trees will be disabled for this thread."<<endl;
125 _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<125<<std::endl;
126 fdchitsel = NULL__null;
127 }
128
129 loop->GetJApplication()->Unlock();
130 }
131
132 DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
133 bfield = dapp->GetBfield(runnumber); // this should be run number based!
134
135}
136
137//---------------------------------
138// ~DTrackHitSelectorALT2 (Destructor)
139//---------------------------------
140DTrackHitSelectorALT2::~DTrackHitSelectorALT2()
141{
142
143}
144
145
146//---------------------------------
147// GetTRDHits
148//---------------------------------
149void DTrackHitSelectorALT2::GetTRDHits(const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DTRDPoint*> &trdhits_in, vector<const DTRDPoint*> &trdhits_out) const {
150 // Vector of pairs storing the hit with the probability it is on the track
151 vector<pair<double,const DTRDPoint*> >trdhits_tmp;
152
153 for (unsigned int k=0;k<extrapolations.size();k++){
154 DVector3 pos=extrapolations[k].position;
155 for (unsigned int j=0;j<trdhits_in.size();j++){
156 const DTRDPoint *hit=trdhits_in[j];
157 if (int(k)==hit->detector){
158 double dx=hit->x-pos.X();
159 double dy=hit->y-pos.Y();
160 double varx=1.,vary=1.;
161 double chisq=dx*dx/varx+dy*dy/vary;
162 // Use chi-sq probability function with Ndof=2 to calculate probability
163 double probability = TMath::Prob(chisq, 2);
164 if(probability>=MIN_HIT_PROB_GEM){
165 pair<double,const DTRDPoint*>myhit;
166 myhit.first=probability;
167 myhit.second=hit;
168 trdhits_tmp.push_back(myhit);
169 }
170 }
171 }
172 }
173 // Order according to layer number and probability,then put the hits in the
174 // output list with the following algorithm: put hits with the highest
175 // probability in a given layer in the output list.
176 sort(trdhits_tmp.begin(),trdhits_tmp.end(),DTrackHitSelector_trdhit_cmp);
177 int old_layer=1000;
178 for (unsigned int i=0;i<trdhits_tmp.size();i++){
179 if (trdhits_tmp[i].second->detector!=old_layer){
180 trdhits_out.push_back(trdhits_tmp[i].second);
181 }
182 old_layer=trdhits_tmp[i].second->detector;
183 }
184}
185
186//---------------------------------
187// GetGEMHits
188//---------------------------------
189void DTrackHitSelectorALT2::GetGEMHits(const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DGEMPoint*> &gemhits_in, vector<const DGEMPoint*> &gemhits_out) const {
190 // Vector of pairs storing the hit with the probability it is on the track
191 vector<pair<double,const DGEMPoint*> >gemhits_tmp;
192
193 for (unsigned int k=0;k<extrapolations.size();k++){
194 DVector3 pos=extrapolations[k].position;
195 for (unsigned int j=0;j<gemhits_in.size();j++){
196 const DGEMPoint *hit=gemhits_in[j];
197 if (int(k)==hit->detector){
198 double dx=hit->x-pos.X();
199 double dy=hit->y-pos.Y();
200 double varx=1.,vary=1.;
201 double chisq=dx*dx/varx+dy*dy/vary;
202 // Use chi-sq probability function with Ndof=2 to calculate probability
203 double probability = TMath::Prob(chisq, 2);
204 if(probability>=MIN_HIT_PROB_GEM){
205 pair<double,const DGEMPoint*>myhit;
206 myhit.first=probability;
207 myhit.second=hit;
208 gemhits_tmp.push_back(myhit);
209 }
210 }
211 }
212
213 }
214 // Order according to layer number and probability,then put the hits in the
215 // output list with the following algorithm: put hits with the highest
216 // probability in a given layer in the output list.
217 sort(gemhits_tmp.begin(),gemhits_tmp.end(),DTrackHitSelector_gemhit_cmp);
218 int old_layer=1000;
219 for (unsigned int i=0;i<gemhits_tmp.size();i++){
220 if (gemhits_tmp[i].second->detector!=old_layer){
221 gemhits_out.push_back(gemhits_tmp[i].second);
222 }
223 old_layer=gemhits_tmp[i].second->detector;
224 }
225
226}
227
228//---------------------------------
229// GetCDCHits
230//---------------------------------
231void DTrackHitSelectorALT2::GetCDCHits(double Bz,double q,
232 const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DCDCTrackHit*> &cdchits_in, vector<const DCDCTrackHit*> &cdchits_out, int N) const
233{
234 // Vector of pairs storing the hit with the probability it is on the track
235 vector<pair<double,const DCDCTrackHit*> >cdchits_tmp;
236
237 // Sort so innermost ring is first and outermost is last
238 vector<const DCDCTrackHit*> cdchits_in_sorted = cdchits_in;
239 sort(cdchits_in_sorted.begin(),cdchits_in_sorted.end(),DTrackHitSelector_cdchit_in_cmp);
240
241 // The variance on the residual due to measurement error.
242 double var=0.25*2.4336*ONE_OVER_120.08333333333333;
243
244 // To estimate the impact of errors in the track momentum on the variance
245 // of the residual, use a helical approximation.
246 double a=-0.003*Bz*q;
247 DVector3 mom=extrapolations[extrapolations.size()/2].momentum;
248 double p=mom.Mag();
249 double p_over_a=p/a;
250 double a_over_p=1./p_over_a;
251 double lambda=M_PI_21.57079632679489661923-mom.Theta();
252 double tanl=tan(lambda);
253 double tanl2=tanl*tanl;
254 double sinl=sin(lambda);
255 double cosl=cos(lambda);
256 double cosl2=cosl*cosl;
257 double sinl2=sinl*sinl;
258 double pt_over_a=cosl*p_over_a;
259
260 // variances
261 double var_lambda_res=0.;
262 double var_phi_radial=0.;
263 double var_x0=0.0,var_y0=0.0;
264 double var_k=0.;
265
266 // Loop over all the CDC hits looking for matches with the track
267 bool outermost_hit=true;
268 vector<const DCDCTrackHit*>::const_reverse_iterator iter;
269 for(iter=cdchits_in_sorted.rbegin(); iter!=cdchits_in_sorted.rend(); iter++){
270 const DCDCTrackHit *hit = *iter;
271 DVector3 origin=hit->wire->origin;
272 DVector3 dir=hit->wire->udir;
273 double ux=dir.x();
274 double uy=dir.y();
275 double uz=dir.z();
276 double z0=origin.z();
277 double d2=0.,d2_old=1.e6;
278 double s=0.;
279 DVector3 old_trackpos;
280 for (unsigned int i=0;i<extrapolations.size();i++){
281 DVector3 trackpos=extrapolations[i].position;
282 double dz=trackpos.z()-z0;
283 DVector3 wirepos=origin+dz/uz*dir;
284 d2=(trackpos-wirepos).Perp2();
285 if (d2>d2_old){
286 if (d2<4){ // has to be reasonably close to the wire in question
287 double phi=extrapolations[i].momentum.Phi();
288 double sinphi=sin(phi);
289 double cosphi=cos(phi);
290 // Variables relating wire direction and track direction
291 double my_ux=ux*sinl-cosl*cosphi;
292 double my_uy=uy*sinl-cosl*sinphi;
293 double denom=my_ux*my_ux+my_uy*my_uy;
294 // For simplicity make a linear approximation for the path increment
295 double ds=((old_trackpos.x()-origin.x()-ux*dz)*my_ux
296 +(old_trackpos.y()-origin.y()-uy*dz)*my_uy)/denom;
297 wirepos+=ds*dir;
298 double dx=old_trackpos.x()+mom.X()/mom.Mag()*ds-wirepos.x();
299 double dy=old_trackpos.y()+mom.Y()/mom.Mag()*ds-wirepos.y();
300 d2_old=d2;
301 d2=dx*dx+dy*dy;
302 if (d2>d2_old){
303 // linear approximation did not work...
304 d2=d2_old;
305 }
306 // Variance in dip angle due to multiple scattering
307 double var_lambda = extrapolations[i].theta2ms_sum/3.;
308 // variance in phi due to multiple scattering
309 double var_phi=var_lambda*(1.+tanl2);
310
311 if (outermost_hit){
312 // variance in the curvature k due to resolution and multiple scattering
313 double s_sq=s*s;
314 double s_4=s_sq*s_sq;
315 double sum_s_theta_ms=extrapolations[i].s_theta_ms_sum;
316 double var_k_ms=(4./3.)*sum_s_theta_ms*sum_s_theta_ms/s_4;
317 double var_k_res=var*0.0720/double(N+4)/s_4/(cosl2*cosl2);
318 var_k=var_k_ms+var_k_res;
319
320 // Variance in dip angle due to measurement error
321 var_lambda_res=12.0*var*double(N-1)/double(N*(N+1))
322 *sinl2*sinl2/s_sq;
323
324 // Variance in phi, crude approximation
325 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
326 var_phi_radial=tan_s_2R*tan_s_2R*pt_over_a*pt_over_a*var_k;
327
328 outermost_hit=false;
329 }
330 // Include error in lambda due to measurements
331 var_lambda+=var_lambda_res;
332
333 // Include uncertainty in phi due to uncertainty in the center of
334 // the circle
335 var_phi+=var_phi_radial;
336
337 // Variance in position due to multiple scattering
338 double var_pos_ms=extrapolations[i].s_theta_ms_sum/3.;
339
340 // Hit doca and residual
341 double doca=sqrt(d2);
342 double dist=0.39;
343 double resi=dist-doca;
344
345 // Variances in x and y due to uncertainty in track parameters
346 double as_over_p=s*a_over_p;
347 double sin_as_over_p=sin(as_over_p);
348 double cos_as_over_p=cos(as_over_p);
349 double one_minus_cos_as_over_p=1-cos_as_over_p;
350 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
351 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
352 double dx_dk=2.*pt_over_a*pt_over_a*(sinphi*diff2-cosphi*diff1);
353 double dx_dcosl=p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
354 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
355 double var_x=var_x0+dx_dk*dx_dk*var_k+var_pos_ms
356 +dx_dcosl*dx_dcosl*sinl*sinl*var_lambda+dx_dphi*dx_dphi*var_phi;
357
358 double dy_dk=-2.*pt_over_a*pt_over_a*(sinphi*diff1+cosphi*diff2);
359 double dy_dcosl=p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
360 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
361 double var_y=var_y0+dy_dk*dy_dk*var_k+var_pos_ms
362 +dy_dcosl*dy_dcosl*sinl*sinl*var_lambda+dy_dphi*dy_dphi*var_phi;
363
364 double dd_dx=dx/doca;
365 double dd_dy=dy/doca;
366 double var_d=dd_dx*dd_dx*var_x+dd_dy*dd_dy*var_y;
367 double chisq=resi*resi/(var+var_d);
368
369 // Use chi-sq probability function with Ndof=1 to calculate probability
370 double probability = TMath::Prob(chisq, 1);
371 if(probability>=MIN_HIT_PROB_CDC){
372 pair<double,const DCDCTrackHit*>myhit;
373 myhit.first=probability;
374 myhit.second=hit;
375 cdchits_tmp.push_back(myhit);
376 }
377
378 // Optionally fill debug tree
379 if(cdchitsel){
380 cdchitdbg.fit_type = kWireBased;
381 cdchitdbg.p = p;
382 cdchitdbg.theta = extrapolations[i].momentum.Theta();
383 //cdchitdbg.mass = mass;
384 cdchitdbg.sigma = sqrt(var+var_d);
385 cdchitdbg.x = old_trackpos.X();
386 cdchitdbg.y = old_trackpos.Y();
387 cdchitdbg.z = old_trackpos.Z();
388 cdchitdbg.s = s;
389 cdchitdbg.itheta02 = extrapolations[i].theta2ms_sum;
390 //cdchitdbg.itheta02s = last_step->itheta02s;
391 //cdchitdbg.itheta02s2 = last_step->itheta02s2;
392 cdchitdbg.dist = dist;
393 cdchitdbg.doca = doca;
394 cdchitdbg.resi = resi;
395 cdchitdbg.chisq = chisq;
396 cdchitdbg.prob = probability;
397 cdchitdbg.sig_phi=sqrt(var_phi);
398 cdchitdbg.sig_lambda=sqrt(var_lambda);
399 cdchitdbg.sig_pt=fabs(2.*pt_over_a)*sqrt(var_k);
400
401 cdchitsel->Fill();
402 }
403
404 if(HS_DEBUG_LEVEL>10){
405 _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<405<<" ";
406 if (probability>=MIN_HIT_PROB_CDC) jerr<<"---> ";
407 jerr<<"s="<<s<<" t=" <<hit->tdrift << " doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" sigma="/*<<sigma_total*/<<" prob="<<probability<<endl;
408 }
409 }
410 break;
411 }
412 d2_old=d2;
413 old_trackpos=trackpos;
414 s=extrapolations[i].s;
415 }
416
417 }
418
419 // Order according to ring number and probability, then put the hits in the
420 // output list with the following algorithm: hits with the highest
421 // probability in a given ring are automatically put in the output list,
422 // but if there is more than one hit in a given ring, only those hits
423 // that are within +/-1 of the straw # of the most probable hit are added
424 // to the list.
425 sort(cdchits_tmp.begin(),cdchits_tmp.end(),DTrackHitSelector_cdchit_cmp);
426 int old_straw=1000,old_ring=1000;
427 for (unsigned int i=0;i<cdchits_tmp.size();i++){
428 if (cdchits_tmp[i].second->wire->ring!=old_ring ||
429 abs(cdchits_tmp[i].second->wire->straw-old_straw)==1){
430 cdchits_out.push_back(cdchits_tmp[i].second);
431 }
432 old_straw=cdchits_tmp[i].second->wire->straw;
433 old_ring=cdchits_tmp[i].second->wire->ring;
434 }
435}
436
437//---------------------------------
438// GetCDCHits
439//---------------------------------
440void DTrackHitSelectorALT2::GetCDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector<const DCDCTrackHit*> &cdchits_in, vector<const DCDCTrackHit*> &cdchits_out, int N) const
441{
442 // Vector of pairs storing the hit with the probability it is on the track
443 vector<pair<double,const DCDCTrackHit*> >cdchits_tmp;
444
445 /// Determine the probability that for each CDC hit that it came from the
446 /// track with the given trajectory.
447 ///
448 /// This will calculate a probability for each CDC hit that
449 /// it came from the track represented by the given
450 /// DReference trajectory. The probability is based on
451 /// the residual between the distance of closest approach
452 /// of the trajectory to the wire and the drift time for
453 /// time-based tracks and the distance to the wire for
454 /// wire-based tracks.
455
456 // Sort so innermost ring is first and outermost is last
457 vector<const DCDCTrackHit*> cdchits_in_sorted = cdchits_in;
458 sort(cdchits_in_sorted.begin(),cdchits_in_sorted.end(),DTrackHitSelector_cdchit_in_cmp);
459
460 // Calculate beta of particle.
461 //double my_mass=rt->GetMass();
462 // double one_over_beta =sqrt(1.0+my_mass*my_mass/rt->swim_steps[0].mom.Mag2());
463
464 // The variance on the residual due to measurement error.
465 double var=0.25*2.4336*ONE_OVER_120.08333333333333;
466
467 // To estimate the impact of errors in the track momentum on the variance of the residual,
468 // use a helical approximation.
469 const DReferenceTrajectory::swim_step_t *my_step=&rt->swim_steps[0];
470 double Bz0=my_step->B.z();
471 double a=-0.003*Bz0*rt->q;
472 double p=my_step->mom.Mag();
473 double p_over_a=p/a;
474 double a_over_p=1./p_over_a;
475 double lambda=M_PI_21.57079632679489661923-my_step->mom.Theta();
476 double cosl=cos(lambda);
477 double sinl=sin(lambda);
478 double sinl2=sinl*sinl;
479 double cosl2=cosl*cosl;
480 double tanl=tan(lambda);
481 double tanl2=tanl*tanl;
482 double pt_over_a=cosl*p_over_a;
483 double phi=my_step->mom.Phi();
484 double cosphi=cos(phi);
485 double sinphi=sin(phi);
486
487 double var_lambda_res=0.;
488 double var_phi_radial=0.;
489 double var_x0=0.0,var_y0=0.0;
490 double mass=rt->GetMass();
491 double one_over_beta=sqrt(1.+mass*mass/(p*p));
492 double var_pt_factor=0.016*one_over_beta/(cosl*a);
493 double var_pt_over_pt_sq=0.;
494
495 // Keep track of straws and rings
496 int old_straw=1000,old_ring=1000;
497
498 // Loop over hits
499 bool outermost_hit=true;
500 vector<const DCDCTrackHit*>::const_reverse_iterator iter;
501 for(iter=cdchits_in_sorted.rbegin(); iter!=cdchits_in_sorted.rend(); iter++){
502 const DCDCTrackHit *hit = *iter;
503
504 // Skip hit if it is on the same wire as the previous hit
505 if (hit->wire->ring == old_ring && hit->wire->straw==old_straw){
506 //_DBG_ << "ring " << hit->wire->ring << " straw " << hit->wire->straw << endl;
507 continue;
508 }
509 old_ring=hit->wire->ring;
510 old_straw=hit->wire->straw;
511
512 // Find the DOCA to this wire
513 double s=0.;
514 double doca = rt->DistToRT(hit->wire, &s);
515
516 if(!isfinite(doca)) continue;
517 if(!isfinite(s))continue;
518 if (s<1.) continue;
519 if (doca>MAX_DOCA) continue;
520
521 const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep();
522
523 // Position along trajectory
524 DVector3 pos=rt->GetLastDOCAPoint();
525
526 // Compensate for the fact that the field at the "origin" of the
527 // track does not correspond to the average Bz used to compute pt
528 double Bratio=last_step->B.z()/Bz0;
529 double invBratio=1./Bratio;
530 pt_over_a*=invBratio;
531 p_over_a*=invBratio;
532 a_over_p*=Bratio;
533
534 // Variances in angles due to multiple scattering
535 double var_lambda = last_step->itheta02/3.;
536 double var_phi=var_lambda*(1.+tanl2);
537
538 if (outermost_hit){
539 // Fractional variance in the curvature k due to resolution and multiple scattering
540 double s_sq=s*s;
541 double var_k_over_k_sq_res=var*p_over_a*p_over_a*0.0720/double(N+4)/(s_sq*s_sq*sinl2)/cosl2;
542 double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->invX0/s;
543 // Fractional variance in pt
544 var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res;
545
546 // Variance in dip angle due to measurement error
547 var_lambda_res=12.0*var*double(N-1)/double(N*(N+1))*cosl2*cosl2/s_sq;
548
549 // Variance in phi, crude approximation
550 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
551 var_phi_radial=tan_s_2R*tan_s_2R*var_pt_over_pt_sq;
552
553 outermost_hit=false;
554 }
555 // Include error in lambda due to measurements
556 var_lambda+=var_lambda_res;
557
558 // Include uncertainty in phi due to uncertainty in the radius of
559 // the circle
560 var_phi+=var_phi_radial;
561
562 // Get "measured" distance to wire.
563 double dist=0.39;
564
565 // Residual
566 double resi = dist - doca;
567
568 // Variance in position due to multiple scattering
569 double var_pos_ms=last_step->itheta02s2/3.;
570
571 // Variances in x and y due to uncertainty in track parameters
572 double as_over_p=s*a_over_p;
573 double sin_as_over_p=sin(as_over_p);
574 double cos_as_over_p=cos(as_over_p);
575 double one_minus_cos_as_over_p=1-cos_as_over_p;
576 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
577 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
578 double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2);
579 double dx_dcosl=p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
580 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
581 double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms
582 +dx_dcosl*dx_dcosl*sinl*sinl*var_lambda+dx_dphi*dx_dphi*var_phi;
583
584 double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2);
585 double dy_dcosl=p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
586 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
587 double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms
588 +dy_dcosl*dy_dcosl*sinl*sinl*var_lambda+dy_dphi*dy_dphi*var_phi;
589
590
591 DVector3 origin=hit->wire->origin;
592 DVector3 dir=hit->wire->udir;
593 double uz=dir.z();
594 double z0=origin.z();
595 DVector3 wirepos=origin+(pos.z()-z0)/uz*dir;
596 double dd_dx=(pos.x()-wirepos.x())/doca;
597 double dd_dy=(pos.y()-wirepos.y())/doca;
598 double var_d=dd_dx*dd_dx*var_x+dd_dy*dd_dy*var_y;
599
600 double chisq=resi*resi/(var+var_d);
601
602 // Use chi-sq probability function with Ndof=1 to calculate probability
603 double probability = TMath::Prob(chisq, 1);
604 if(probability>=MIN_HIT_PROB_CDC){
605 pair<double,const DCDCTrackHit*>myhit;
606 myhit.first=probability;
607 myhit.second=hit;
608 cdchits_tmp.push_back(myhit);
609 }
610
611 // Optionally fill debug tree
612 if(cdchitsel){
613 cdchitdbg.fit_type = fit_type;
614 cdchitdbg.p = p;
615 cdchitdbg.theta = rt->swim_steps[0].mom.Theta();
616 cdchitdbg.mass = mass;
617 cdchitdbg.sigma = sqrt(var+var_d);
618 cdchitdbg.x = pos.X();
619 cdchitdbg.y = pos.Y();
620 cdchitdbg.z = pos.Z();
621 cdchitdbg.s = s;
622 cdchitdbg.itheta02 = last_step->itheta02;
623 cdchitdbg.itheta02s = last_step->itheta02s;
624 cdchitdbg.itheta02s2 = last_step->itheta02s2;
625 cdchitdbg.dist = dist;
626 cdchitdbg.doca = doca;
627 cdchitdbg.resi = resi;
628 cdchitdbg.chisq = chisq;
629 cdchitdbg.prob = probability;
630 cdchitdbg.sig_phi=sqrt(var_phi);
631 cdchitdbg.sig_lambda=sqrt(var_lambda);
632 cdchitdbg.sig_pt=sqrt(var_pt_over_pt_sq);
633
634 cdchitsel->Fill();
635 }
636
637 if(HS_DEBUG_LEVEL>10){
638 _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<638<<" ";
639 if (probability>=MIN_HIT_PROB_CDC) jerr<<"---> ";
640 jerr<<"s="<<s<<" t=" <<hit->tdrift << " doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" sigma="/*<<sigma_total*/<<" prob="<<probability<<endl;
641 }
642 }
643
644 // Order according to ring number and probability, then put the hits in the
645 // output list with the following algorithm: hits with the highest
646 // probability in a given ring are automatically put in the output list,
647 // but if there is more than one hit in a given ring, only those hits
648 // that are within +/-1 of the straw # of the most probable hit are added
649 // to the list.
650 sort(cdchits_tmp.begin(),cdchits_tmp.end(),DTrackHitSelector_cdchit_cmp);
651 old_straw=1000,old_ring=1000;
652 for (unsigned int i=0;i<cdchits_tmp.size();i++){
653 if (cdchits_tmp[i].second->wire->ring!=old_ring ||
654 abs(cdchits_tmp[i].second->wire->straw-old_straw)==1){
655 cdchits_out.push_back(cdchits_tmp[i].second);
656 }
657 old_straw=cdchits_tmp[i].second->wire->straw;
658 old_ring=cdchits_tmp[i].second->wire->ring;
659 }
660}
661
662//---------------------------------
663// GetFDCHits
664//---------------------------------
665void DTrackHitSelectorALT2::GetFDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector<const DFDCPseudo*> &fdchits_in, vector<const DFDCPseudo*> &fdchits_out,int N) const
666{
667 // Vector of pairs storing the hit with the probability it is on the track
668 vector<pair<double,const DFDCPseudo*> >fdchits_tmp;
669
670 /// Determine the probability that for each FDC hit that it came from the
671 /// track with the given trajectory.
672 ///
673 /// This will calculate a probability for each FDC hit that
674 /// it came from the track represented by the given
675 /// DReference trajectory. The probability is based on
676 /// the residual between the distance of closest approach
677 /// of the trajectory to the wire and the drift distance
678 /// and the distance along the wire.
679
680 // Sort so innermost ring is first and outermost is last
681 vector<const DFDCPseudo*> fdchits_in_sorted = fdchits_in;
682 sort(fdchits_in_sorted.begin(),fdchits_in_sorted.end(),DTrackHitSelector_fdchit_in_cmp);
683
684 // The variance on the residual due to measurement error.
685 double var_anode = 0.25*ONE_OVER_120.08333333333333;
686 const double VAR_CATHODE_STRIPS=0.000225;
687 double var_cathode = VAR_CATHODE_STRIPS;
688
689 // To estimate the impact of errors in the track momentum on the variance of the residual,
690 // use a helical approximation.
691 const DReferenceTrajectory::swim_step_t *my_step=&rt->swim_steps[0];
692 double Bz0=my_step->B.z();
693 double z0=my_step->origin.z();
694 double a=-0.003*Bz0*rt->q;
695 double p=my_step->mom.Mag();
696 double p_sq=p*p;
697 double p_over_a=p/a;
698 double a_over_p=1./p_over_a;
699 double lambda=M_PI_21.57079632679489661923-my_step->mom.Theta();
700 double cosl=cos(lambda);
701 double cosl2=cosl*cosl;
702 double sinl=sin(lambda);
703 double sinl2=sinl*sinl;
704 double tanl=tan(lambda);
705 double tanl2=tanl*tanl;
706 double pt_over_a=cosl*p_over_a;
707 double phi=my_step->mom.Phi();
708 double cosphi=cos(phi);
709 double sinphi=sin(phi);
710 double var_lambda=0.,var_phi=0.,var_lambda_res=0.;
711 double var_phi_radial=0.;
712 double mass=rt->GetMass();
713
714 double var_x0=0.0,var_y0=0.0;
715 double var_pt_over_pt_sq=0.;
716
717 // Loop over hits
718 bool most_downstream_hit=true;
719 vector<const DFDCPseudo*>::const_reverse_iterator iter;
720 for(iter=fdchits_in_sorted.rbegin(); iter!=fdchits_in_sorted.rend(); iter++){
721 const DFDCPseudo *hit = *iter;
722
723 // Find the DOCA to this wire
724 double s=0.;
725 double doca = rt->DistToRT(hit->wire, &s);
726
727 if(!isfinite(doca)) continue;
728 if(!isfinite(s))continue;
729 if (s<1.) continue;
730 if (doca>MAX_DOCA)continue;
731
732 const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep();
733
734 // Position along trajectory
735 DVector3 pos=rt->GetLastDOCAPoint();
736 double dz=pos.z()-z0;
737
738 // Direction variables for wire
739 double cosa=hit->wire->udir.y();
740 double sina=hit->wire->udir.x();
741
742 // Cathode Residual
743 double u=rt->GetLastDistAlongWire();
744 double u_cathodes = hit->s;
745 double resic = u - u_cathodes;
746
747 // Get "measured" distance to wire.
748 double dist=0.25;
749
750 // Take into account non-normal incidence to FDC plane
751 double pz=last_step->mom.z();
752 double tx=last_step->mom.x()/pz;
753 double ty=last_step->mom.y()/pz;
754 double tu=tx*cosa-ty*sina;
755 double alpha=atan(tu);
756 double cosalpha=cos(alpha);
757
758 // Compensate for the fact that the field at the "origin" of the
759 // track does not correspond to the average Bz used to compute pt
760 double Bz=last_step->B.z();
761 double Bratio=Bz/Bz0;
762 double invBratio=1./Bratio;
763 pt_over_a*=invBratio;
764 p_over_a*=invBratio;
765 a_over_p*=Bratio;
766
767 // Anode Residual
768 double resi = dist - doca/cosalpha;
769
770 // Initialize some probability-related variables needed later
771 double probability=0.,chisq=0.;
772
773 if (fit_type!=kHelical){
774 // Correct for deflection of avalanche position due to Lorentz force
775 double sign=(pos.x()*cosa-pos.y()*sina-hit->w)<0?1:-1.;
776 double ucor=0.1458*Bz*(1.-0.048*last_step->B.Perp())*sign*doca;
777 resic-=ucor;
778 }
779 else{
780 // Cathode variance due to Lorentz deflection
781 double max_deflection=0.1458*Bz*(1.-0.048*last_step->B.Perp())*0.5;
782 var_cathode=VAR_CATHODE_STRIPS+max_deflection*max_deflection/12.;
783 }
784 double var_tot=var_anode+var_cathode;
785
786 // Variance in angles due to multiple scattering
787 var_lambda = last_step->itheta02/3.;
788 var_phi=var_lambda*(1.+tanl2);
789
790 if (most_downstream_hit){
791 // Fractional variance in the curvature k due to resolution and multiple scattering
792 double s_sq=s*s;
793 double var_k_over_k_sq_res=var_tot*p_over_a*p_over_a
794 *0.0720/double(N+4)/(s_sq*s_sq)/cosl2;
795
796
797 double one_over_beta=sqrt(1.+mass*mass/p_sq);
798 double var_pt_factor=0.016*one_over_beta/(cosl*0.003*last_step->B.z());
799 double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->invX0/s;
800 // Fractional variance in pt
801 var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res;
802
803 // Variance in dip angle due to measurement error
804 var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1))
805 *sinl2*sinl2/s_sq;
806
807 // Variance in phi, crude approximation
808 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
809 var_phi_radial=tan_s_2R*tan_s_2R*var_pt_over_pt_sq;
810
811 most_downstream_hit=false;
812 }
813
814 // Include error in lambda due to measurements
815 var_lambda+=var_lambda_res;
816
817 // Include uncertainty in phi due to uncertainty in the center of
818 // the circle
819 var_phi+=var_phi_radial;
Value stored to 'var_phi' is never read
820
821 var_phi=0.;
822 var_lambda=0.;
823
824 // Variance in position due to multiple scattering
825 double var_pos_ms=last_step->itheta02s2/3.;
826
827 // Variances in x and y due to uncertainty in track parameters
828 double as_over_p=s*a_over_p;
829 double sin_as_over_p=sin(as_over_p);
830 double cos_as_over_p=cos(as_over_p);
831 double one_minus_cos_as_over_p=1-cos_as_over_p;
832 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
833 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
834 double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2);
835 double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p);
836 double ds_dcosl=dz*cosl/(sinl*sinl2);
837 double dx_dcosl
838 =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p)
839 +dx_ds*ds_dcosl;
840 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
841 double var_z0=2.*tanl2*(var_tot)*double(2*N-1)/double(N*(N+1));
842
843 double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms
844 +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi
845 +dx_ds*dx_ds*var_z0/sinl2;
846
847 double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2);
848 double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p);
849 double dy_dcosl
850 =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p)
851 +dy_ds*ds_dcosl;
852 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
853 double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms
854 +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi
855 +dy_ds*dy_ds*var_z0/sinl2;
856
857 // Rotate from global coordinate system into FDC local system
858 double cos2a=cosa*cosa;
859 double sin2a=sina*sina;
860 double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha);
861 double var_u=cos2a*var_y+sin2a*var_x;
862
863 if (fit_type!=kHelical){
864 // Factors take into account improved resolution after wire-based fit
865 var_d*=0.1;
866 var_u*=0.1;
867 }
868
869 // Calculate chisq
870 chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode);
871
872 // Probability of this hit being on the track
873 probability = TMath::Prob(chisq,2);
874
875 if(probability>=MIN_HIT_PROB_FDC){
876 pair<double,const DFDCPseudo*>myhit;
877 myhit.first=probability;
878 myhit.second=hit;
879 fdchits_tmp.push_back(myhit);
880 }
881
882 // Optionally fill debug tree
883 if(fdchitsel){
884 fdchitdbg.fit_type = fit_type;
885 fdchitdbg.p = p;
886 fdchitdbg.theta = rt->swim_steps[0].mom.Theta();
887 fdchitdbg.mass = mass;
888 fdchitdbg.sigma_anode = sqrt(var_d+var_anode);
889 fdchitdbg.sigma_cathode = sqrt(var_u+var_cathode);
890 fdchitdbg.x = pos.X();
891 fdchitdbg.y = pos.Y();
892 fdchitdbg.z = pos.Z();
893 fdchitdbg.s = s;
894 fdchitdbg.itheta02 = last_step->itheta02;
895 fdchitdbg.itheta02s = last_step->itheta02s;
896 fdchitdbg.itheta02s2 = last_step->itheta02s2;
897 fdchitdbg.dist = dist;
898 fdchitdbg.doca = doca;
899 fdchitdbg.resi = resi;
900 fdchitdbg.u = u;
901 fdchitdbg.u_cathodes = u_cathodes;
902 fdchitdbg.resic = resic;
903 fdchitdbg.chisq = chisq;
904 fdchitdbg.prob = probability;
905 fdchitdbg.sig_phi=sqrt(var_phi);
906 fdchitdbg.sig_lambda=sqrt(var_lambda);
907 fdchitdbg.sig_pt=sqrt(var_pt_over_pt_sq);
908
909 fdchitsel->Fill();
910 }
911
912 if(HS_DEBUG_LEVEL>10){
913 _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<913<<" ";
914 if(probability>=MIN_HIT_PROB_FDC)jerr<<"----> ";
915 jerr<<"s="<<s<<" doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" resic="<<resic<<" chisq="<<chisq<<" prob="<<probability<<endl;
916 }
917 }
918 // Order according to layer number and probability,then put the hits in the
919 // output list with the following algorithm: hits with the highest
920 // probability in a given layer are automatically put in the output list,
921 // but if there is more than one hit in a given layer, only those hits
922 // that are within +/-1 of the wire # of the most probable hit are added
923 // to the list.
924 sort(fdchits_tmp.begin(),fdchits_tmp.end(),DTrackHitSelector_fdchit_cmp);
925 int old_layer=1000,old_wire=1000;
926 for (unsigned int i=0;i<fdchits_tmp.size();i++){
927 if (fdchits_tmp[i].second->wire->layer!=old_layer ||
928 abs(fdchits_tmp[i].second->wire->wire-old_wire)==1){
929 fdchits_out.push_back(fdchits_tmp[i].second);
930 }
931 old_wire=fdchits_tmp[i].second->wire->wire;
932 old_layer=fdchits_tmp[i].second->wire->layer;
933 }
934}
935void DTrackHitSelectorALT2::GetFDCHits(double Bz,double q,
936 const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DFDCPseudo*> &fdchits_in, vector<const DFDCPseudo*> &fdchits_out,int N) const
937{
938 // Vector of pairs storing the hit with the probability it is on the track
939 vector<pair<double,const DFDCPseudo*> >fdchits_tmp;
940
941 /// Determine the probability that for each FDC hit that it came from the
942 /// track with the given trajectory. The probability is based on
943 /// the residual between the distance of closest approach
944 /// of the trajectory to the wire and the drift distance
945 /// and the distance along the wire.
946
947 // Sort so innermost plane is first and outermost is last
948 vector<const DFDCPseudo*> fdchits_in_sorted = fdchits_in;
949 sort(fdchits_in_sorted.begin(),fdchits_in_sorted.end(),DTrackHitSelector_fdchit_in_cmp);
950
951 // The variance on the residual due to measurement error.
952 double var_anode = 0.25*ONE_OVER_120.08333333333333;
953 const double VAR_CATHODE_STRIPS=0.000225;
954 double var_cathode = VAR_CATHODE_STRIPS+0.15*0.15*ONE_OVER_120.08333333333333;
955 double var_tot=var_anode+var_cathode;
956
957 // To estimate the impact of errors in the track momentum on the variance of the residual,
958 // use a helical approximation.
959 double a=-0.003*Bz*q;
960 DVector3 mom=extrapolations[extrapolations.size()/2].momentum;
961 double p=mom.Mag();
962 double p_over_a=p/a;
963 double a_over_p=1./p_over_a;
964 double lambda=M_PI_21.57079632679489661923-mom.Theta();
965 double tanl=tan(lambda);
966 double tanl2=tanl*tanl;
967 double sinl=sin(lambda);
968 double cosl=cos(lambda);
969 double cosl2=cosl*cosl;
970 double sinl2=sinl*sinl;
971 double pt_over_a=cosl*p_over_a;
972 double z0=extrapolations[0].position.z();
973 // variances
974 double var_lambda=0.,var_phi=0.,var_lambda_res=0.,var_k=0.;
975 double var_phi_radial=0.;
976 double var_x0=0.0,var_y0=0.0;
977 double var_z0=2.*tanl2*(var_tot)*double(2*N-1)/double(N*(N+1));
978
979 // Loop over hits
980 bool most_downstream_hit=true;
981 int last_extrapolation_index=extrapolations.size()-1;
982 vector<const DFDCPseudo*>::const_reverse_iterator iter;
983 for(iter=fdchits_in_sorted.rbegin(); iter!=fdchits_in_sorted.rend(); iter++){
984 const DFDCPseudo *hit = *iter;
985 for (int k=last_extrapolation_index;k>=0;k--){
986 // Position along trajectory
987 DVector3 pos=extrapolations[k].position;
988 double dz=pos.z()-z0;
989 double s=extrapolations[k].s;
990 if (fabs(pos.z()-hit->wire->origin.z())<0.5){
991 // Variance in dip angle due to multiple scattering
992 var_lambda = extrapolations[k].theta2ms_sum/3.;
993 // the above expression seems to lead to overestimation of the uncertainty in the dip angle after the wire-based pass..
994 var_lambda*=0.01;
995 // Variance in phi due to multiple scattering
996 var_phi=var_lambda*(1.+tanl2);
997
998 // azimuthal angle
999 double phi=extrapolations[k].momentum.Phi();
1000 double sinphi=sin(phi);
1001 double cosphi=cos(phi);
1002
1003 if (most_downstream_hit){
1004 // variance in the curvature k due to resolution and multiple scattering
1005 double s_sq=s*s;
1006 double s_4=s_sq*s_sq;
1007 double sum_s_theta_ms=extrapolations[k].s_theta_ms_sum;
1008 double var_k_ms=(4./3.)*sum_s_theta_ms*sum_s_theta_ms/s_4;
1009 double var_k_res=var_tot*0.0720/double(N+4)/s_4/(cosl2*cosl2);
1010 var_k=var_k_ms+var_k_res;
1011
1012 // Variance in dip angle due to measurement error
1013 var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1))
1014 *sinl2*sinl2/s_sq;
1015
1016 // Variance in phi, crude approximation
1017 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
1018 var_phi_radial=tan_s_2R*tan_s_2R*pt_over_a*pt_over_a*var_k;
1019
1020 most_downstream_hit=false;
1021 }
1022 // Include error in lambda due to measurements
1023 var_lambda+=var_lambda_res;
1024
1025 // Include uncertainty in phi due to uncertainty in the center of
1026 // the circle
1027 var_phi+=var_phi_radial;
1028
1029 // Variance in position due to multiple scattering
1030 double var_pos_ms=extrapolations[k].s_theta_ms_sum/3.;
1031
1032 // doca to wire
1033 double x=pos.x();
1034 double y=pos.y();
1035 double dx=hit->xy.X()-x;
1036 double dy=hit->xy.Y()-y;
1037 double doca=sqrt(dx*dx+dy*dy);
1038
1039 // Direction variables for wire
1040 double cosa=hit->wire->udir.y();
1041 double sina=hit->wire->udir.x();
1042
1043 // Cathode Residual
1044 double u=x*sina+y*cosa;
1045 double u_cathodes = hit->s;
1046 double resic = u - u_cathodes;
1047
1048 // Get "measured" distance to wire.
1049 double dist=0.25;
1050
1051 // Take into account non-normal incidence to FDC plane
1052 double pz=extrapolations[k].momentum.z();
1053 double tx=extrapolations[k].momentum.x()/pz;
1054 double ty=extrapolations[k].momentum.y()/pz;
1055 double tu=tx*cosa-ty*sina;
1056 double alpha=atan(tu);
1057 double cosalpha=cos(alpha);
1058
1059 // Anode Residual
1060 double resi = dist - doca/cosalpha;
1061
1062 // Initialize some probability-related variables needed later
1063 double probability=0.,chisq=0.;
1064
1065 // Variances in x and y due to uncertainty in track parameters
1066 double as_over_p=s*a_over_p;
1067 double sin_as_over_p=sin(as_over_p);
1068 double cos_as_over_p=cos(as_over_p);
1069 double one_minus_cos_as_over_p=1-cos_as_over_p;
1070 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
1071 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
1072 double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p);
1073 double ds_dcosl=dz*cosl/(sinl*sinl2);
1074 double dx_dcosl
1075 =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p)
1076 +dx_ds*ds_dcosl;
1077 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
1078 double dx_dk=2.*pt_over_a*pt_over_a*(sinphi*diff2-cosphi*diff1);
1079
1080 double var_x=var_x0+dx_dk*dx_dk*var_k+var_pos_ms
1081 +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi
1082 +dx_ds*dx_ds*var_z0/sinl2;
1083
1084 double dy_dk=-2.*pt_over_a*pt_over_a*(sinphi*diff1+cosphi*diff2);
1085 double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p);
1086 double dy_dcosl
1087 =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p)
1088 +dy_ds*ds_dcosl;
1089 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
1090 double var_y=var_y0+dy_dk*dy_dk*var_k+var_pos_ms
1091 +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi
1092 +dy_ds*dy_ds*var_z0/sinl2;
1093
1094 // Rotate from global coordinate system into FDC local system
1095 double cos2a=cosa*cosa;
1096 double sin2a=sina*sina;
1097 double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha);
1098 double var_u=cos2a*var_y+sin2a*var_x;
1099
1100 // Factors take into account improved resolution after wire-based fit
1101 //var_d*=0.1;
1102 //var_u*=0.1;
1103
1104 // Calculate chisq
1105 chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode);
1106
1107 // Probability of this hit being on the track
1108 probability = TMath::Prob(chisq,2);
1109
1110 if(probability>=MIN_HIT_PROB_FDC){
1111 pair<double,const DFDCPseudo*>myhit;
1112 myhit.first=probability;
1113 myhit.second=hit;
1114 fdchits_tmp.push_back(myhit);
1115 }
1116 // Optionally fill debug tree
1117 if(fdchitsel){
1118 fdchitdbg.fit_type = kWireBased;
1119 fdchitdbg.p = p;
1120 fdchitdbg.theta = extrapolations[k].momentum.Theta();
1121 //fdchitdbg.mass = mass;
1122 fdchitdbg.sigma_anode = sqrt(var_d+var_anode);
1123 fdchitdbg.sigma_cathode = sqrt(var_u+var_cathode);
1124 fdchitdbg.x = pos.X();
1125 fdchitdbg.y = pos.Y();
1126 fdchitdbg.z = pos.Z();
1127 fdchitdbg.s = s;
1128 fdchitdbg.itheta02 = extrapolations[k].theta2ms_sum;
1129 //fdchitdbg.itheta02s = last_step->itheta02s;
1130 //fdchitdbg.itheta02s2 = last_step->itheta02s2;
1131 fdchitdbg.dist = dist;
1132 fdchitdbg.doca = doca;
1133 fdchitdbg.resi = resi;
1134 fdchitdbg.u = u;
1135 fdchitdbg.u_cathodes = u_cathodes;
1136 fdchitdbg.resic = resic;
1137 fdchitdbg.chisq = chisq;
1138 fdchitdbg.prob = probability;
1139 fdchitdbg.sig_phi=sqrt(var_phi);
1140 fdchitdbg.sig_lambda=sqrt(var_lambda);
1141 fdchitdbg.sig_pt=fabs(2.*pt_over_a)*sqrt(var_k);
1142
1143 fdchitsel->Fill();
1144 }
1145
1146 if(HS_DEBUG_LEVEL>10){
1147 _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc"
<<":"<<1147<<" ";
1148 if(probability>=MIN_HIT_PROB_FDC)jerr<<"----> ";
1149 jerr<<"s="<<s<<" doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" resic="<<resic<<" chisq="<<chisq<<" prob="<<probability<<endl;
1150 }
1151
1152 }
1153 }
1154
1155 }
1156 // Order according to layer number and probability,then put the hits in the
1157 // output list with the following algorithm: hits with the highest
1158 // probability in a given layer are automatically put in the output list,
1159 // but if there is more than one hit in a given layer, only those hits
1160 // that are within +/-1 of the wire # of the most probable hit are added
1161 // to the list.
1162 sort(fdchits_tmp.begin(),fdchits_tmp.end(),DTrackHitSelector_fdchit_cmp);
1163 int old_layer=1000,old_wire=1000;
1164 for (unsigned int i=0;i<fdchits_tmp.size();i++){
1165 if (fdchits_tmp[i].second->wire->layer!=old_layer ||
1166 abs(fdchits_tmp[i].second->wire->wire-old_wire)==1){
1167 fdchits_out.push_back(fdchits_tmp[i].second);
1168 }
1169 old_wire=fdchits_tmp[i].second->wire->wire;
1170 old_layer=fdchits_tmp[i].second->wire->layer;
1171 }
1172}