Bug Summary

File:libraries/PID/DParticleID.cc
Location:line 571, column 7
Description:Value stored to 'myphi' is never read

Annotated Source Code

1// $Id$
2//
3// File: DParticleID.cc
4// Created: Mon Feb 28 14:48:56 EST 2011
5// Creator: staylor (on Linux ifarml1 2.6.18-128.el5 x86_64)
6//
7
8#include "DParticleID.h"
9#include <TRACKING/DTrackFitter.h>
10#include <TMath.h>
11#include "FCAL/DFCALGeometry.h"
12
13#define C_EFFECTIVE15. 15. // start counter light propagation speed
14#define OUT_OF_TIME_CUT200. 200.
15
16// Routine for sorting dEdx data
17bool static DParticleID_dedx_cmp(DParticleID::dedx_t a,DParticleID::dedx_t b){
18 return a.dEdx < b.dEdx;
19}
20
21// Routine for sorting hypotheses according to FOM
22bool static DParticleID_hypothesis_cmp(const DTrackTimeBased *a,
23 const DTrackTimeBased *b){
24 return (a->FOM>b->FOM);
25}
26
27
28//---------------------------------
29// DParticleID (Constructor)
30//---------------------------------
31DParticleID::DParticleID(JEventLoop *loop)
32{
33
34 DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
35 if(!dapp){
36 _DBG_std::cerr<<"DParticleID.cc"<<":"<<36<<
" "<<"Cannot get DApplication from JEventLoop! (are you using a JApplication based program?)"<<endl;
37 return;
38 }
39 const DRootGeom *RootGeom = dapp->GetRootGeom();
40 bfield = dapp->GetBfield();
41 stepper= new DMagneticFieldStepper(bfield);
42
43 // Get material properties for chamber gas
44 double rho_Z_over_A_LnI=0,radlen=0;
45
46 RootGeom->FindMat("CDchamberGas", dRhoZoverA_CDC, rho_Z_over_A_LnI, radlen);
47 dLnI_CDC = rho_Z_over_A_LnI/dRhoZoverA_CDC;
48 dKRhoZoverA_CDC = 0.1535E-3*dRhoZoverA_CDC;
49
50 RootGeom->FindMat("FDchamberGas", dRhoZoverA_FDC, rho_Z_over_A_LnI, radlen);
51 dLnI_FDC = rho_Z_over_A_LnI/dRhoZoverA_FDC;
52 dKRhoZoverA_FDC = 0.1535E-3*dRhoZoverA_FDC;
53
54 // Get the geometry
55 geom = dapp->GetDGeometry(loop->GetJEvent().GetRunNumber());
56
57 vector<double>sc_origin;
58 geom->Get("//posXYZ[@volume='StartCntr']/@X_Y_Z",sc_origin);
59
60 //vector<double>sc_light_guide(3);
61 //geom->Get("//tubs[@name='STLG']/@Rio_Z",sc_light_guide);
62 //sc_light_guide_length=sc_light_guide[2];
63
64 vector<vector<double> > sc_rioz;
65 geom->GetMultiple("//pgon[@name='STRC']/polyplane/@Rio_Z", sc_rioz);
66
67 for (unsigned int k=0;k<sc_rioz.size()-1;k++){
68 DVector3 pos((sc_rioz[k][0]+sc_rioz[k][1])/2.,0.,sc_rioz[k][2]+sc_origin[2]);
69 DVector3 dir(sc_rioz[k+1][2]-sc_rioz[k][2],0,
70 -sc_rioz[k+1][0]+sc_rioz[k][0]);
71 dir.SetMag(1.);
72
73 sc_pos.push_back(pos);
74 sc_norm.push_back(dir);
75 }
76 // sc_leg_tcor=(sc_light_guide[2]-sc_pos[0].z())/C_EFFECTIVE;
77 sc_leg_tcor=-sc_pos[0].z()/C_EFFECTIVE15.;
78 double theta=sc_norm[sc_norm.size()-1].Theta();
79 sc_angle_cor=1./cos(M_PI3.14159265358979323846-theta);
80
81 //Get calibration constants
82 map<string, double> locPIDParams;
83 if ( !loop->GetCalib("PID/photon_track_matching", locPIDParams)){
84 cout<<"DParticleID: loading values from PID data base"<<endl;
85 DELTA_R_BCAL = locPIDParams["DELTA_R_BCAL"];
86 DELTA_R_FCAL = locPIDParams["DELTA_R_FCAL"];
87 } else {
88 cout << "DParticleID: Error loading values from PID data base" <<endl;
89 DELTA_R_BCAL = 15.0;
90 DELTA_R_FCAL = 15.0;
91 }
92
93 dTargetZCenter = 0.0;
94 geom->GetTargetZ(dTargetZCenter);
95}
96
97//---------------------------------
98// ~DParticleID (Destructor)
99//---------------------------------
100DParticleID::~DParticleID()
101{
102 if (stepper) delete stepper;
103}
104
105// Group fitted tracks according to candidate id
106jerror_t DParticleID::GroupTracks(vector<const DTrackTimeBased *> &tracks,
107 vector<vector<const DTrackTimeBased*> >&grouped_tracks){
108 if (tracks.size()==0) return RESOURCE_UNAVAILABLE;
109
110 JObject::oid_t old_id=tracks[0]->candidateid;
111 vector<const DTrackTimeBased *>hypotheses;
112 for (unsigned int i=0;i<tracks.size();i++){
113 const DTrackTimeBased *track=tracks[i];
114
115 if (old_id != track->candidateid){
116 // Sort according to FOM
117 sort(hypotheses.begin(),hypotheses.end(),DParticleID_hypothesis_cmp);
118
119 // Add to the grouped_tracks vector
120 grouped_tracks.push_back(hypotheses);
121
122 // Clear the hypothesis list for the next track
123 hypotheses.clear();
124 }
125 hypotheses.push_back(track);
126 old_id=track->candidateid;
127 }
128 // Final set
129 sort(hypotheses.begin(),hypotheses.end(),DParticleID_hypothesis_cmp);
130 grouped_tracks.push_back(hypotheses);
131
132 return NOERROR;
133}
134
135// Compute the energy losses and the path lengths in the chambers for each hit
136// on the track. Returns a list of dE and dx pairs with the momentum at the
137// hit.
138jerror_t DParticleID::GetDCdEdxHits(const DTrackTimeBased *track, vector<dedx_t>& dEdxHits_CDC, vector<dedx_t>& dEdxHits_FDC){
139 // Position and momentum
140 DVector3 pos,mom;
141
142 //dE and dx pairs
143 pair<double,double>de_and_dx;
144
145 // We cast away the const-ness of the reference trajectory so that we can use the DisToRT method
146 DReferenceTrajectory *my_rt=const_cast<DReferenceTrajectory*>(track->rt);
147
148 //Get the list of cdc hits used in the fit
149 vector<const DCDCTrackHit*>cdchits;
150 track->GetT(cdchits);
151
152 // Loop over cdc hits
153 for (unsigned int i=0;i<cdchits.size();i++){
154 double locReturnValue = my_rt->DistToRT(cdchits[i]->wire);
155 if(!((locReturnValue >= 0.0) || (locReturnValue <= 0.0)))
156 continue; //NaN
157
158 my_rt->GetLastDOCAPoint(pos, mom);
159
160 // Create the dE,dx pair from the position and momentum using a helical approximation for the path
161 // in the straw and keep track of the momentum in the active region of the detector
162 if (CalcdEdxHit(mom,pos,cdchits[i],de_and_dx)==NOERROR)
163 dEdxHits_CDC.push_back(dedx_t(de_and_dx.first, de_and_dx.second, mom.Mag()));
164 }
165
166 //Get the list of fdc hits used in the fit
167 vector<const DFDCPseudo*>fdchits;
168 track->GetT(fdchits);
169
170 // loop over fdc hits
171 for (unsigned int i=0;i<fdchits.size();i++){
172 double locReturnValue = my_rt->DistToRT(fdchits[i]->wire);
173 if(!((locReturnValue >= 0.0) || (locReturnValue <= 0.0)))
174 continue; //NaN
175 my_rt->GetLastDOCAPoint(pos, mom);
176
177 double gas_thickness = 1.0; // cm
178 dEdxHits_FDC.push_back(dedx_t(fdchits[i]->dE, gas_thickness/cos(mom.Theta()), mom.Mag()));
179 }
180
181 // Sort the dEdx entries from smallest to largest
182 sort(dEdxHits_FDC.begin(),dEdxHits_FDC.end(),DParticleID_dedx_cmp);
183 sort(dEdxHits_CDC.begin(),dEdxHits_CDC.end(),DParticleID_dedx_cmp);
184
185 return NOERROR;
186}
187
188jerror_t DParticleID::CalcDCdEdx(const DTrackTimeBased *locTrackTimeBased, double& locdEdx_FDC, double& locdx_FDC, double& locdEdx_CDC, double& locdx_CDC, unsigned int& locNumHitsUsedFordEdx_FDC, unsigned int& locNumHitsUsedFordEdx_CDC){
189 vector<dedx_t> locdEdxHits_CDC, locdEdxHits_FDC;
190 jerror_t locReturnStatus = GetDCdEdxHits(locTrackTimeBased, locdEdxHits_CDC, locdEdxHits_FDC);
191 if(locReturnStatus != NOERROR){
192 locdEdx_FDC = NaNstd::numeric_limits<double>::quiet_NaN();
193 locdx_FDC = NaNstd::numeric_limits<double>::quiet_NaN();
194 locNumHitsUsedFordEdx_FDC = 0;
195 locdEdx_CDC = NaNstd::numeric_limits<double>::quiet_NaN();
196 locdx_CDC = NaNstd::numeric_limits<double>::quiet_NaN();
197 locNumHitsUsedFordEdx_CDC = 0;
198 return locReturnStatus;
199 }
200 return CalcDCdEdx(locTrackTimeBased, locdEdxHits_CDC, locdEdxHits_FDC, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdx_CDC, locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC);
201}
202
203jerror_t DParticleID::CalcDCdEdx(const DTrackTimeBased *locTrackTimeBased, const vector<dedx_t>& locdEdxHits_CDC, const vector<dedx_t>& locdEdxHits_FDC, double& locdEdx_FDC, double& locdx_FDC, double& locdEdx_CDC, double& locdx_CDC, unsigned int& locNumHitsUsedFordEdx_FDC, unsigned int& locNumHitsUsedFordEdx_CDC){
204 locdx_CDC = 0.0;
205 locdEdx_CDC = 0.0;
206 locNumHitsUsedFordEdx_CDC = locdEdxHits_CDC.size()/2;
207 if(locNumHitsUsedFordEdx_CDC > 0){
208 for(unsigned int loc_i = 0; loc_i < locNumHitsUsedFordEdx_CDC; ++loc_i){
209 locdEdx_CDC += locdEdxHits_CDC[loc_i].dE; //weight is ~ #e- (scattering sites): dx!
210 locdx_CDC += locdEdxHits_CDC[loc_i].dx;
211 }
212 locdEdx_CDC /= locdx_CDC;
213 }
214
215 locdx_FDC = 0.0;
216 locdEdx_FDC = 0.0;
217 locNumHitsUsedFordEdx_FDC = locdEdxHits_FDC.size()/2;
218 if(locNumHitsUsedFordEdx_FDC > 0){
219 for(unsigned int loc_i = 0; loc_i < locNumHitsUsedFordEdx_FDC; ++loc_i){
220 locdEdx_FDC += locdEdxHits_FDC[loc_i].dE; //weight is ~ #e- (scattering sites): dx!
221 locdx_FDC += locdEdxHits_FDC[loc_i].dx;
222 }
223 locdEdx_FDC /= locdx_FDC; //weight is dx/dx_total
224 }
225 return NOERROR;
226}
227
228// Calculate the path length for a single hit in a straw and return ds and the
229// energy deposition in the straw. It returns dE as the first element in the
230// dedx pair and ds as the second element in the dedx pair.
231jerror_t DParticleID::CalcdEdxHit(const DVector3 &mom,
232 const DVector3 &pos,
233 const DCDCTrackHit *hit,
234 pair <double,double> &dedx){
235 if (hit==NULL__null || hit->wire==NULL__null) return RESOURCE_UNAVAILABLE;
236
237 // Track direction parameters
238 double phi=mom.Phi();
239 double lambda=M_PI_21.57079632679489661923-mom.Theta();
240 double cosphi=cos(phi);
241 double sinphi=sin(phi);
242 double tanl=tan(lambda);
243
244 //Position relative to wire origin
245 double dz=pos.z()-hit->wire->origin.z();
246 double dx=pos.x()-hit->wire->origin.x();
247 double dy=pos.y()-hit->wire->origin.y();
248
249 // square of straw radius
250 double rs2=0.776*0.776;
251
252 // Useful temporary variables related to the direction of the wire
253 double ux=hit->wire->udir.x();
254 double uy=hit->wire->udir.y();
255 double uz=hit->wire->udir.z();
256 double A=1.-ux*ux;
257 double B=-2.*ux*uy;
258 double C=-2.*ux*uz;
259 double D=-2.*uy*uz;
260 double E=1.-uy*uy;
261 double F=1.-uz*uz;
262
263 // The path length in the straw is given by s=sqrt(b*b-4*a*c)/a/cosl.
264 // a, b, and c follow.
265 double a=A*cosphi*cosphi+B*cosphi*sinphi+C*cosphi*tanl+D*sinphi*tanl
266 +E*sinphi*sinphi+F*tanl*tanl;
267 double b=2.*A*dx*cosphi+B*dx*sinphi+B*dy*cosphi+C*dx*tanl+C*cosphi*dz
268 +D*dy*tanl+D*sinphi*dz+2.*E*dy*sinphi+2.*F*dz*tanl;
269 double c=A*dx*dx+B*dx*dy+C*dx*dz+D*dy*dz+E*dy*dy+F*dz*dz-rs2;
270
271 // Check for valid arc length and compute dEdx
272 double temp=b*b-4.*a*c;
273 if (temp>0){
274 double cosl=fabs(cos(lambda));
275 // double gas_density=0.0018;
276
277 // arc length and energy deposition
278 //dedx.second=gas_density*sqrt(temp)/a/cosl; // g/cm^2
279 dedx.second=sqrt(temp)/a/cosl;
280 dedx.first=hit->dE; //GeV
281
282 return NOERROR;
283 }
284
285 return VALUE_OUT_OF_RANGE;
286}
287
288//------------------
289// MatchToTOF
290//------------------
291// Loop over TOF points, looking for minimum distance of closest approach
292// of track to a point in the TOF and using this to check for a match.
293//
294// NOTE: an initial guess for tproj is expected as input so that out-of-time
295// hits can be skipped
296jerror_t DParticleID::MatchToTOF(const DReferenceTrajectory *rt, DTrackFitter::fit_type_t fit_type, vector<const DTOFPoint*>&tof_points, double &tproj, unsigned int &tof_match_id, double &locPathLength, double &locFlightTime){
297 //tproj=NaN;
298 tof_match_id=0;
299 if (tof_points.size()==0){
300 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
301 return RESOURCE_UNAVAILABLE;
302 }
303
304 double dmin=10000.;
305 // loop over tof points
306 double locTempPathLength;
307 for (unsigned int k=0;k<tof_points.size();k++){
308
309 // Get the TOF cluster position and normal vector for the TOF plane
310 DVector3 tof_pos=tof_points[k]->pos;
311 DVector3 norm(0,0,1);
312 DVector3 proj_pos,dir;
313
314 // Find the distance of closest approach between the track trajectory
315 // and the tof cluster position, looking for the minimum
316 double locTempFlightTime=0.;
317 rt->GetIntersectionWithPlane(tof_pos,norm,proj_pos,dir,&locTempPathLength,&locTempFlightTime);
318 double d=(tof_pos-proj_pos).Mag();
319
320 // Check that the hit is not out of time with respect to the track
321 if (fabs(tof_points[k]->t - locTempFlightTime - tproj) > OUT_OF_TIME_CUT200.) continue;
322
323 if (d<dmin){
324 dmin=d;
325 locPathLength = locTempPathLength;
326 locFlightTime = locTempFlightTime;
327 tof_match_id=k;
328 }
329 }
330
331 // Check for a match
332 // double p=rt->swim_steps[0].mom.Mag();
333 double match_cut=0.;
334 match_cut = 6.15; //current dPositionMatchCut_DoubleEnded variable in DTOFPoint_factory.cc
335// if (fit_type==DTrackFitter::kTimeBased) match_cut=3.624+0.488/p;
336// else match_cut=4.0+0.488/p;
337 if (dmin<match_cut){
338 // Projected time at the target
339 tproj=tof_points[tof_match_id]->t - locFlightTime;
340
341 return NOERROR;
342 }
343
344 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
345 return VALUE_OUT_OF_RANGE;
346}
347
348//------------------
349// MatchToBCAL
350//------------------
351// Loop over bcal clusters, looking for minimum distance of closest approach
352// of track to a cluster and using this to check for a match.
353//
354// NOTE: an initial guess for tproj is expected as input so that out-of-time
355// hits can be skipped
356jerror_t DParticleID::MatchToBCAL(const DReferenceTrajectory *rt, const vector<const DBCALShower*>& locInputBCALShowers, vector<const DBCALShower*>& locMatchedBCALShowers, double& locProjectedTime, double& locPathLength, double& locFlightTime){
357 if (locInputBCALShowers.size() == 0){
358 locProjectedTime = NaNstd::numeric_limits<double>::quiet_NaN();
359 return RESOURCE_UNAVAILABLE;
360 }
361
362 //Loop over bcal showers
363 double locTempPathLength, dphi, dz;
364 double locLargestShowerEnergy = -1.0;
365 int locBestShowerMatchIndex = -1, locBestShowerMatchIndexInMatchVector = -1;
366 for (unsigned int k = 0; k < locInputBCALShowers.size(); ++k){
367 // Get the BCAL cluster position and normal
368 const DBCALShower *shower = locInputBCALShowers[k];
369 DVector3 bcal_pos(shower->x, shower->y, shower->z);
370 DVector3 proj_pos;
371 // and the bcal cluster position, looking for the minimum
372 double locTempFlightTime = 0.0;
373 double d = rt->DistToRTwithTime(bcal_pos, &locTempPathLength, &locTempFlightTime);
374 if(!finite(d))
375 continue;
376 // Check that the hit is not out of time with respect to the track
377 if (fabs(locInputBCALShowers[k]->t - locTempFlightTime - locProjectedTime) > OUT_OF_TIME_CUT200.) continue;
378
379 proj_pos = rt->GetLastDOCAPoint();
380
381 dz = proj_pos.z() - bcal_pos.z();
382 dphi = proj_pos.Phi() - bcal_pos.Phi();
383 double p = rt->swim_steps[0].mom.Mag();
384 dphi += 0.002 + 8.314e-3/(p + 0.3788)/(p + 0.3788);
385 while(dphi > TMath::Pi())
386 dphi -= 2.0*TMath::Pi();
387 while(dphi < -1.0*TMath::Pi())
388 dphi += 2.0*TMath::Pi();
389 double phi_sigma = 0.025 + 5.8e-4/(p*p*p);
390 if (fabs(dz) < 10. && fabs(dphi) < 3.*phi_sigma){
391 locMatchedBCALShowers.push_back(shower);
392 if(shower->E > locLargestShowerEnergy){
393 locLargestShowerEnergy = shower->E;
394 locPathLength = locTempPathLength;
395 locFlightTime = locTempFlightTime;
396 locBestShowerMatchIndex = k;
397 locBestShowerMatchIndexInMatchVector = locMatchedBCALShowers.size() - 1;
398 }
399 }
400 }
401
402 if(locBestShowerMatchIndexInMatchVector > 0){ //move highest energy shower to the front of the list
403 locMatchedBCALShowers.erase(locMatchedBCALShowers.begin() + locBestShowerMatchIndexInMatchVector);
404 locMatchedBCALShowers.insert(locMatchedBCALShowers.begin(), locInputBCALShowers[locBestShowerMatchIndex]);
405 }
406
407 if(locMatchedBCALShowers.size() > 0){
408 // Projected time at the target
409 locProjectedTime = locMatchedBCALShowers[0]->t - locFlightTime;
410 // locProjectedTime += 0.03; // correction determined from observation of simulated data //Paul Mattione
411 return NOERROR;
412 }
413
414 locProjectedTime = NaNstd::numeric_limits<double>::quiet_NaN();
415 return VALUE_OUT_OF_RANGE;
416}
417
418//------------------
419// MatchToFCAL
420//------------------
421// Loop over fcal clusters, looking for minimum distance of closest approach
422// of track to a cluster and using this to check for a match.
423//
424// NOTE: an initial guess for locProjectedTime is expected as input so that out-of-time
425// hits can be skipped
426jerror_t DParticleID::MatchToFCAL(const DReferenceTrajectory *rt, const vector<const DFCALShower*>& locInputFCALShowers, vector<const DFCALShower*>& locMatchedFCALShowers, double& locProjectedTime, double& locPathLength, double& locFlightTime){
427 if (locInputFCALShowers.size() == 0){
428 locProjectedTime = NaNstd::numeric_limits<double>::quiet_NaN();
429 return RESOURCE_UNAVAILABLE;
430 }
431
432 // Set minimum matching distance to a large default value
433 double locLargestShowerEnergy = -1.0;
434 // loop over clusters
435 double locTempPathLength;
436 int locBestShowerMatchIndex = -1, locBestShowerMatchIndexInMatchVector = -1;
437 for (unsigned int k = 0;k < locInputFCALShowers.size(); k++){
438
439 // Get the FCAL cluster position and normal vector for the FCAL plane
440 DVector3 fcal_pos = locInputFCALShowers[k]->getPosition();
441 // This is a bit of a kludge...
442 //fcal_pos.SetZ(DFCALGeometry::fcalFaceZ());
443 DVector3 norm(0,0,1);
444 DVector3 proj_pos, dir;
445
446 // Find the distance of closest approach between the track trajectory
447 // and the tof cluster position, looking for the minimum
448 double locTempFlightTime = 0.;
449 rt->GetIntersectionWithPlane(fcal_pos, norm, proj_pos, dir, &locTempPathLength, &locTempFlightTime);
450 double d = (fcal_pos - proj_pos).Mag();
451 // Check that the hit is not out of time with respect to the track
452 if (fabs(locInputFCALShowers[k]->getTime() - locTempFlightTime - locProjectedTime) > OUT_OF_TIME_CUT200.) continue;
453 if(d < DELTA_R_FCAL){
454 locMatchedFCALShowers.push_back(locInputFCALShowers[k]);
455 if(locInputFCALShowers[k]->getEnergy() > locLargestShowerEnergy){
456 locLargestShowerEnergy = locInputFCALShowers[k]->getEnergy();
457 locPathLength = locTempPathLength;
458 locFlightTime = locTempFlightTime;
459 locBestShowerMatchIndex = k;
460 locBestShowerMatchIndexInMatchVector = locMatchedFCALShowers.size() - 1;
461 }
462 }
463 }
464
465 if(locBestShowerMatchIndexInMatchVector > 0){ //move highest energy shower to the front of the list
466 locMatchedFCALShowers.erase(locMatchedFCALShowers.begin() + locBestShowerMatchIndexInMatchVector);
467 locMatchedFCALShowers.insert(locMatchedFCALShowers.begin(), locInputFCALShowers[locBestShowerMatchIndex]);
468 }
469
470 // Check for a match
471 if(locMatchedFCALShowers.size() > 0){
472 locProjectedTime = locMatchedFCALShowers[0]->getTime() - locFlightTime;
473 locProjectedTime -= 2.218; // correction determined from fit to simulated data
474 return NOERROR;
475 }
476
477 locProjectedTime = NaNstd::numeric_limits<double>::quiet_NaN();
478 return VALUE_OUT_OF_RANGE;
479}
480
481//------------------
482// MatchToSC
483//------------------
484// Match track to the start counter paddles with hits. If a match
485// is found, use the z-position of the track projection to the start counter
486// planes to correct for the light propagation within the scintillator and
487// estimate the "vertex" time.
488//
489// Unlike the next method, this method does not use the reference trajectory.
490//
491// NOTE: an initial guess for tproj is expected as input so that out-of-time
492// hits can be skipped
493jerror_t DParticleID::MatchToSC(const DKinematicData &parms,
494 vector<const DSCHit*>&sc_hits,
495 double &tproj,unsigned int &sc_match_id){
496 sc_match_id=0;
497 if (sc_hits.size()==0){
498 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
499 return RESOURCE_UNAVAILABLE;
500 }
501 double myz=0.;
502 double dphi_min=10000.,myphi=0.;
503 unsigned int num=sc_norm.size()-1;
504
505 DVector3 pos(parms.position());
506 DVector3 mom(parms.momentum());
507 stepper->SetCharge(parms.charge());
508
509 // Swim to barrel representing the start counter straight portion
510 double ds=0.,myds=0;
511 if (stepper->SwimToRadius(pos,mom,sc_pos[1].x(),&ds)){
512 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
513 return VALUE_OUT_OF_RANGE;
514 }
515 // Change the sign of ds -- most of the time we will need to add a small
516 // amount of time to the time calculated at the start counter
517 ds*=-1.;
518
519 // Position along z and phi at intersection
520 double proj_z=pos.z();
521 double proj_phi=pos.Phi();
522 if (proj_phi<0) proj_phi+=2.*M_PI3.14159265358979323846;
523
524 // Position of transition between start counter nose and leg
525 double sc_pos1=sc_pos[1].z();
526 // position in z at the end of the nose
527 double sc_posn=sc_pos[num].z();
528
529 // loop over sc hits, looking for the one with closest phi value to the
530 // projected phi value
531 for (unsigned int i=0;i<sc_hits.size();i++){
532 // Check that the hit is not out of time with respect to the track
533 if (fabs(tproj-sc_hits[i]->t)>OUT_OF_TIME_CUT200.) continue;
534
535 double phi=(6.0+12.*(sc_hits[i]->sector-1))*M_PI3.14159265358979323846/180.;
536 double dphi=phi-proj_phi;
537
538 // If the z position is in the nose region, match to the appropriate start
539 // counter plane
540 if (proj_z>sc_pos1){
541 double cosphi=cos(phi);
542 double sinphi=sin(phi);
543 for (unsigned int i=1;i<num;i++){
544 DVector3 mymom=(-1.)*mom;
545 DVector3 mypos=pos;
546 double xhat=sc_norm[i].x();
547 double r=sc_pos[i].X();
548 double x=r*cosphi;
549 double y=r*sinphi;
550 double z=sc_pos[i].z();
551 DVector3 norm(x*xhat,y*xhat,z);
552 DVector3 plane(x,y,z);
553 double ds2=0.;
554 if (stepper->SwimToPlane(mypos,mymom,plane,norm,&ds2)) continue;
555
556 proj_z=mypos.z();
557 if (proj_z<sc_pos[i+1].z()){
558 proj_phi=mypos.Phi();
559 if (proj_phi<0) proj_phi+=2.*M_PI3.14159265358979323846;
560 dphi=proj_phi-phi;
561 ds+=ds2;
562
563 break;
564 }
565 }
566 }
567
568 // Look for smallest difference in phi
569 if (fabs(dphi)<dphi_min){
570 dphi_min=fabs(dphi);
571 myphi=phi;
Value stored to 'myphi' is never read
572 myz=proj_z;
573 myds=ds;
574 sc_match_id=i;
575 }
576 }
577
578 // Look for a match in phi
579 if (dphi_min<0.16){
580 // Find the time at the start counter, before correcting for position
581 // along the start counter
582 tproj=sc_hits[sc_match_id]->t-sc_leg_tcor;
583
584 // Check position along z
585 if (myz<sc_pos[0].z()) myz=sc_pos[0].z();
586 if (myz<sc_pos1){
587 // Leg region
588 tproj-=myz/C_EFFECTIVE15.;
589 }
590 else if (myz<sc_posn){
591 // Nose region
592 tproj-=((myz-sc_pos1)*sc_angle_cor+sc_pos1)/C_EFFECTIVE15.;
593 }
594 else{
595 // Assume that the particle hit the most downstream z position of the
596 // start counter
597 tproj-=((sc_posn-sc_pos1)*sc_angle_cor+sc_pos1)/C_EFFECTIVE15.;
598 }
599
600 // Adjust for flight time to the position pos
601 double mass=parms.mass();
602 double p2=mom.Mag2();
603 double one_over_beta=sqrt(1.+mass*mass/p2);
604 tproj += myds*one_over_beta/SPEED_OF_LIGHT29.9792;
605
606 return NOERROR;
607 }
608
609 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
610 return VALUE_OUT_OF_RANGE;
611
612}
613
614
615//------------------
616// MatchToSC
617//------------------
618// Match track to the start counter paddles with hits. If a match
619// is found, use the z-position of the track projection to the start counter
620// planes to correct for the light propagation within the scintillator and
621// estimate the "vertex" time.
622//
623// NOTE: an initial guess for tproj is expected as input so that out-of-time
624// hits can be skipped
625jerror_t DParticleID::MatchToSC(const DReferenceTrajectory *rt, DTrackFitter::fit_type_t fit_type, vector<const DSCHit*>&sc_hits, double &tproj,unsigned int &sc_match_id, double &locPathLength, double &locFlightTime){
626
627 //tproj=NaN;
628 sc_match_id=0;
629 if (sc_hits.size()==0){
630 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
631 return RESOURCE_UNAVAILABLE;
632 }
633 double myz=0.;
634 double dphi_min=10000.,myphi=0.;
635 double locTempPathLength, locTempFlightTime = 0.0;
636 DVector3 proj_pos;
637
638 // Find intersection with a "barrel" approximation for the start counter
639 rt->GetIntersectionWithRadius(sc_pos[1].x(),proj_pos,&locTempPathLength,&locTempFlightTime);
640 double proj_phi=proj_pos.Phi();
641 if (proj_phi<0) proj_phi+=2.*M_PI3.14159265358979323846;
642
643 // loop over sc hits, looking for the one with closest phi value to the
644 // projected phi value
645 for (unsigned int i=0;i<sc_hits.size();i++){
646 // Check that the hit is not out of time with respect to the track
647 if (fabs(tproj-sc_hits[i]->t)>OUT_OF_TIME_CUT200.) continue;
648
649 double phi=(6.0+12.*(sc_hits[i]->sector-1))*M_PI3.14159265358979323846/180.;
650 double dphi=phi-proj_phi;
651
652 if (fabs(dphi)<dphi_min){
653 dphi_min=fabs(dphi);
654 myphi=phi;
655 myz=proj_pos.z();
656 sc_match_id=i;
657 }
658 }
659 // Look for a match in phi
660 if (dphi_min<0.21){
661 // Now check to see if the intersection is in the nose region and find the
662 // start time
663 tproj=sc_hits[sc_match_id]->t-sc_leg_tcor;
664 if (myz<sc_pos[0].z()) myz=sc_pos[0].z();
665 if (myz>sc_pos[1].z()){
666 unsigned int num=sc_norm.size()-1;
667 for (unsigned int i=1;i<num;i++){
668 double xhat=sc_norm[i].x();
669 DVector3 norm(cos(myphi)*xhat,sin(myphi)*xhat,sc_norm[i].z());
670 double r=sc_pos[i].X();
671 DVector3 pos(r*cos(myphi),r*sin(myphi),sc_pos[i].z());
672 rt->GetIntersectionWithPlane(pos,norm,proj_pos,&locTempPathLength,&locTempFlightTime);
673 myz=proj_pos.z();
674 if (myz<sc_pos[i+1].z())
675 break;
676 }
677 double sc_pos1=sc_pos[1].z();
678 if (myz<sc_pos1){
679 tproj-=locTempFlightTime+sc_pos1/C_EFFECTIVE15.;
680 locFlightTime = locTempFlightTime;
681 locPathLength = locTempPathLength;
682 }else if (myz>sc_pos[num].z()){
683 // Assume that the particle hit the most downstream z position of the
684 // start counter
685 double costheta=rt->swim_steps[0].mom.CosTheta();
686 double s=(sc_pos[num].z()-rt->swim_steps[0].origin.z())/costheta;
687 double mass=rt->GetMass();
688 double p2=rt->swim_steps[0].mom.Mag2();
689 double one_over_beta=sqrt(1.+mass*mass/p2);
690
691 tproj -= s*one_over_beta/SPEED_OF_LIGHT29.9792 + ((sc_pos[num].z()-sc_pos1)*sc_angle_cor+sc_pos1)/C_EFFECTIVE15.;
692 locFlightTime = s*one_over_beta/SPEED_OF_LIGHT29.9792;
693 locPathLength = s;
694 }else{
695 tproj-=locTempFlightTime+((myz-sc_pos1)*sc_angle_cor+sc_pos1)/C_EFFECTIVE15.;
696 locFlightTime = locTempFlightTime;
697 locPathLength = locTempPathLength;
698 }
699 }else{
700 tproj-=locTempFlightTime+myz/C_EFFECTIVE15.;
701 locFlightTime = locTempFlightTime;
702 locPathLength = locTempPathLength;
703 }
704 return NOERROR;
705 }
706
707 tproj=NaNstd::numeric_limits<double>::quiet_NaN();
708 return VALUE_OUT_OF_RANGE;
709}
710
711void DParticleID::Calc_TimingChiSq(DChargedTrackHypothesis* locChargedTrackHypothesis, double locRFTime, double locRFBunchFrequency)
712{
713 if((locChargedTrackHypothesis->t0_detector() == SYS_NULL) || (locChargedTrackHypothesis->t1_detector() == SYS_NULL) || (locChargedTrackHypothesis->t1_detector() == SYS_START))
714 {
715 //uncertainty so huge on SYS_START that for t1() it won't help distinguish PID anyway
716 locChargedTrackHypothesis->dChiSq_Timing = 0.0;
717 locChargedTrackHypothesis->dNDF_Timing = 0;
718 return;
719 }
720
721 // Use ST hit to select RF beam bucket
722 double locPropagatedRFTime = locRFTime + (locChargedTrackHypothesis->z() - dTargetZCenter)/SPEED_OF_LIGHT29.9792;
723 double locSTRFTimeDifference = locChargedTrackHypothesis->t0() - locPropagatedRFTime;
724 while(fabs(locSTRFTimeDifference) > locRFBunchFrequency/2.0){
725 locPropagatedRFTime += (locSTRFTimeDifference > 0.0) ? locRFBunchFrequency : -1.0*locRFBunchFrequency;
726 locSTRFTimeDifference = locChargedTrackHypothesis->t0() - locPropagatedRFTime;
727 }
728
729 // Compare time difference between RF & TOF/BCAL/FCAL times at the vertex
730 double locTimeDifference = locPropagatedRFTime - locChargedTrackHypothesis->dProjectedStartTime;
731
732 // Calculate ChiSq, FOM
733 double locTUncertainty = locChargedTrackHypothesis->dProjectedStartTimeUncertainty;
734 double locTimingChiSq = locTimeDifference*locTimeDifference/(locTUncertainty*locTUncertainty);
735 locChargedTrackHypothesis->dChiSq_Timing = locTimingChiSq;
736 locChargedTrackHypothesis->dNDF_Timing = 1;
737}
738
739Particle_t DParticleID::IDTrack(float locCharge, float locMass){
740 float locMassTolerance = 0.010;
741 if (locCharge > 0.0){ // Positive particles
742 if (fabs(locMass - ParticleMass(Proton)) < locMassTolerance) return Proton;
743 if (fabs(locMass - ParticleMass(PiPlus)) < locMassTolerance) return PiPlus;
744 if (fabs(locMass - ParticleMass(KPlus)) < locMassTolerance) return KPlus;
745 if (fabs(locMass - ParticleMass(Positron)) < locMassTolerance) return Positron;
746 if (fabs(locMass - ParticleMass(MuonPlus)) < locMassTolerance) return MuonPlus;
747 } else if (locCharge < 0.0){ // Negative particles
748 if (fabs(locMass - ParticleMass(PiMinus)) < locMassTolerance) return PiMinus;
749 if (fabs(locMass - ParticleMass(KMinus)) < locMassTolerance) return KMinus;
750 if (fabs(locMass - ParticleMass(MuonMinus)) < locMassTolerance) return MuonMinus;
751 if (fabs(locMass - ParticleMass(Electron)) < locMassTolerance) return Electron;
752 if (fabs(locMass - ParticleMass(AntiProton)) < locMassTolerance) return AntiProton;
753 } else { //Neutral Track
754 if (fabs(locMass - ParticleMass(Gamma)) < locMassTolerance) return Gamma;
755 if (fabs(locMass - ParticleMass(Neutron)) < locMassTolerance) return Neutron;
756 }
757 return Unknown;
758}
759