1 | |
2 | |
3 | |
4 | |
5 | |
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 | |
17 | bool static DParticleID_dedx_cmp(DParticleID::dedx_t a,DParticleID::dedx_t b){ |
18 | return a.dEdx < b.dEdx; |
19 | } |
20 | |
21 | |
22 | bool static DParticleID_hypothesis_cmp(const DTrackTimeBased *a, |
23 | const DTrackTimeBased *b){ |
24 | return (a->FOM>b->FOM); |
25 | } |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | DParticleID::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 | |
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 | |
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 | |
61 | |
62 | |
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 | |
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 | |
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 | |
99 | |
100 | DParticleID::~DParticleID() |
101 | { |
102 | if (stepper) delete stepper; |
103 | } |
104 | |
105 | |
106 | jerror_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 | |
117 | sort(hypotheses.begin(),hypotheses.end(),DParticleID_hypothesis_cmp); |
118 | |
119 | |
120 | grouped_tracks.push_back(hypotheses); |
121 | |
122 | |
123 | hypotheses.clear(); |
124 | } |
125 | hypotheses.push_back(track); |
126 | old_id=track->candidateid; |
127 | } |
128 | |
129 | sort(hypotheses.begin(),hypotheses.end(),DParticleID_hypothesis_cmp); |
130 | grouped_tracks.push_back(hypotheses); |
131 | |
132 | return NOERROR; |
133 | } |
134 | |
135 | |
136 | |
137 | |
138 | jerror_t DParticleID::GetDCdEdxHits(const DTrackTimeBased *track, vector<dedx_t>& dEdxHits_CDC, vector<dedx_t>& dEdxHits_FDC){ |
139 | |
140 | DVector3 pos,mom; |
141 | |
142 | |
143 | pair<double,double>de_and_dx; |
144 | |
145 | |
146 | DReferenceTrajectory *my_rt=const_cast<DReferenceTrajectory*>(track->rt); |
147 | |
148 | |
149 | vector<const DCDCTrackHit*>cdchits; |
150 | track->GetT(cdchits); |
151 | |
152 | |
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; |
157 | |
158 | my_rt->GetLastDOCAPoint(pos, mom); |
159 | |
160 | |
161 | |
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 | |
167 | vector<const DFDCPseudo*>fdchits; |
168 | track->GetT(fdchits); |
169 | |
170 | |
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; |
175 | my_rt->GetLastDOCAPoint(pos, mom); |
176 | |
177 | double gas_thickness = 1.0; |
178 | dEdxHits_FDC.push_back(dedx_t(fdchits[i]->dE, gas_thickness/cos(mom.Theta()), mom.Mag())); |
179 | } |
180 | |
181 | |
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 | |
188 | jerror_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 | |
203 | jerror_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; |
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; |
221 | locdx_FDC += locdEdxHits_FDC[loc_i].dx; |
222 | } |
223 | locdEdx_FDC /= locdx_FDC; |
224 | } |
225 | return NOERROR; |
226 | } |
227 | |
228 | |
229 | |
230 | |
231 | jerror_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 | |
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 | |
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 | |
250 | double rs2=0.776*0.776; |
251 | |
252 | |
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 | |
264 | |
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 | |
272 | double temp=b*b-4.*a*c; |
273 | if (temp>0){ |
274 | double cosl=fabs(cos(lambda)); |
275 | |
276 | |
277 | |
278 | |
279 | dedx.second=sqrt(temp)/a/cosl; |
280 | dedx.first=hit->dE; |
281 | |
282 | return NOERROR; |
283 | } |
284 | |
285 | return VALUE_OUT_OF_RANGE; |
286 | } |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | |
295 | |
296 | jerror_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 | |
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 | |
306 | double locTempPathLength; |
307 | for (unsigned int k=0;k<tof_points.size();k++){ |
308 | |
309 | |
310 | DVector3 tof_pos=tof_points[k]->pos; |
311 | DVector3 norm(0,0,1); |
312 | DVector3 proj_pos,dir; |
313 | |
314 | |
315 | |
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 | |
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 | |
332 | |
333 | double match_cut=0.; |
334 | match_cut = 6.15; |
335 | |
336 | |
337 | if (dmin<match_cut){ |
338 | |
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 | |
350 | |
351 | |
352 | |
353 | |
354 | |
355 | |
356 | jerror_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 | |
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){ |
| 2 | | Loop condition is true. Entering loop body | |
|
| 7 | | Loop condition is true. Entering loop body | |
|
| 14 | | Loop condition is false. Execution continues on line 402 | |
|
367 | |
368 | const DBCALShower *shower = locInputBCALShowers[k]; |
369 | DVector3 bcal_pos(shower->x, shower->y, shower->z); |
370 | DVector3 proj_pos; |
371 | |
372 | double locTempFlightTime = 0.0; |
373 | double d = rt->DistToRTwithTime(bcal_pos, &locTempPathLength, &locTempFlightTime); |
374 | if(!finite(d)) |
| |
| |
375 | continue; |
376 | |
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()) |
| 5 | | Loop condition is false. Execution continues on line 387 | |
|
| 10 | | Loop condition is false. Execution continues on line 387 | |
|
386 | dphi -= 2.0*TMath::Pi(); |
387 | while(dphi < -1.0*TMath::Pi()) |
| 6 | | Loop condition is false. Execution continues on line 389 | |
|
| 11 | | Loop condition is false. Execution continues on line 389 | |
|
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){ |
| |
403 | locMatchedBCALShowers.erase(locMatchedBCALShowers.begin() + locBestShowerMatchIndexInMatchVector); |
404 | locMatchedBCALShowers.insert(locMatchedBCALShowers.begin(), locInputBCALShowers[locBestShowerMatchIndex]); |
405 | } |
406 | |
407 | if(locMatchedBCALShowers.size() > 0){ |
| |
408 | |
409 | locProjectedTime = locMatchedBCALShowers[0]->t - locFlightTime; |
| 17 | | Dereference of null pointer |
|
410 | |
411 | return NOERROR; |
412 | } |
413 | |
414 | locProjectedTime = NaNstd::numeric_limits<double>::quiet_NaN(); |
415 | return VALUE_OUT_OF_RANGE; |
416 | } |
417 | |
418 | |
419 | |
420 | |
421 | |
422 | |
423 | |
424 | |
425 | |
426 | jerror_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 | |
433 | double locLargestShowerEnergy = -1.0; |
434 | |
435 | double locTempPathLength; |
436 | int locBestShowerMatchIndex = -1, locBestShowerMatchIndexInMatchVector = -1; |
437 | for (unsigned int k = 0;k < locInputFCALShowers.size(); k++){ |
438 | |
439 | |
440 | DVector3 fcal_pos = locInputFCALShowers[k]->getPosition(); |
441 | |
442 | |
443 | DVector3 norm(0,0,1); |
444 | DVector3 proj_pos, dir; |
445 | |
446 | |
447 | |
448 | double locTempFlightTime = 0.; |
449 | rt->GetIntersectionWithPlane(fcal_pos, norm, proj_pos, dir, &locTempPathLength, &locTempFlightTime); |
450 | double d = (fcal_pos - proj_pos).Mag(); |
451 | |
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){ |
466 | locMatchedFCALShowers.erase(locMatchedFCALShowers.begin() + locBestShowerMatchIndexInMatchVector); |
467 | locMatchedFCALShowers.insert(locMatchedFCALShowers.begin(), locInputFCALShowers[locBestShowerMatchIndex]); |
468 | } |
469 | |
470 | |
471 | if(locMatchedFCALShowers.size() > 0){ |
472 | locProjectedTime = locMatchedFCALShowers[0]->getTime() - locFlightTime; |
473 | locProjectedTime -= 2.218; |
474 | return NOERROR; |
475 | } |
476 | |
477 | locProjectedTime = NaNstd::numeric_limits<double>::quiet_NaN(); |
478 | return VALUE_OUT_OF_RANGE; |
479 | } |
480 | |
481 | |
482 | |
483 | |
484 | |
485 | |
486 | |
487 | |
488 | |
489 | |
490 | |
491 | |
492 | |
493 | jerror_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 | |
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 | |
516 | |
517 | ds*=-1.; |
518 | |
519 | |
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 | |
525 | double sc_pos1=sc_pos[1].z(); |
526 | |
527 | double sc_posn=sc_pos[num].z(); |
528 | |
529 | |
530 | |
531 | for (unsigned int i=0;i<sc_hits.size();i++){ |
532 | |
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 | |
539 | |
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 | |
569 | if (fabs(dphi)<dphi_min){ |
570 | dphi_min=fabs(dphi); |
571 | myphi=phi; |
572 | myz=proj_z; |
573 | myds=ds; |
574 | sc_match_id=i; |
575 | } |
576 | } |
577 | |
578 | |
579 | if (dphi_min<0.16){ |
580 | |
581 | |
582 | tproj=sc_hits[sc_match_id]->t-sc_leg_tcor; |
583 | |
584 | |
585 | if (myz<sc_pos[0].z()) myz=sc_pos[0].z(); |
586 | if (myz<sc_pos1){ |
587 | |
588 | tproj-=myz/C_EFFECTIVE15.; |
589 | } |
590 | else if (myz<sc_posn){ |
591 | |
592 | tproj-=((myz-sc_pos1)*sc_angle_cor+sc_pos1)/C_EFFECTIVE15.; |
593 | } |
594 | else{ |
595 | |
596 | |
597 | tproj-=((sc_posn-sc_pos1)*sc_angle_cor+sc_pos1)/C_EFFECTIVE15.; |
598 | } |
599 | |
600 | |
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 | |
617 | |
618 | |
619 | |
620 | |
621 | |
622 | |
623 | |
624 | |
625 | jerror_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 | |
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 | |
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 | |
644 | |
645 | for (unsigned int i=0;i<sc_hits.size();i++){ |
646 | |
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 | |
660 | if (dphi_min<0.21){ |
661 | |
662 | |
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 | |
684 | |
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 | |
711 | void 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 | |
716 | locChargedTrackHypothesis->dChiSq_Timing = 0.0; |
717 | locChargedTrackHypothesis->dNDF_Timing = 0; |
718 | return; |
719 | } |
720 | |
721 | |
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 | |
730 | double locTimeDifference = locPropagatedRFTime - locChargedTrackHypothesis->dProjectedStartTime; |
731 | |
732 | |
733 | double locTUncertainty = locChargedTrackHypothesis->dProjectedStartTimeUncertainty; |
734 | double locTimingChiSq = locTimeDifference*locTimeDifference/(locTUncertainty*locTUncertainty); |
735 | locChargedTrackHypothesis->dChiSq_Timing = locTimingChiSq; |
736 | locChargedTrackHypothesis->dNDF_Timing = 1; |
737 | } |
738 | |
739 | Particle_t DParticleID::IDTrack(float locCharge, float locMass){ |
740 | float locMassTolerance = 0.010; |
741 | if (locCharge > 0.0){ |
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){ |
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 { |
754 | if (fabs(locMass - ParticleMass(Gamma)) < locMassTolerance) return Gamma; |
755 | if (fabs(locMass - ParticleMass(Neutron)) < locMassTolerance) return Neutron; |
756 | } |
757 | return Unknown; |
758 | } |
759 | |