1 | |
2 | |
3 | |
4 | |
5 | |
6 | #include <cassert> |
7 | #include <math.h> |
8 | using namespace std; |
9 | |
10 | #include "DDIRCLut.h" |
11 | #include "DANA/DApplication.h" |
12 | #include <JANA/JCalibration.h> |
13 | |
14 | |
15 | |
16 | |
17 | DDIRCLut::DDIRCLut() |
18 | { |
19 | DIRC_DEBUG_HISTS = false; |
20 | gPARMS->SetDefaultParameter("DIRC:DEBUG_HISTS",DIRC_DEBUG_HISTS); |
21 | |
22 | DIRC_TRUTH_BARHIT = false; |
23 | gPARMS->SetDefaultParameter("DIRC:TRUTH_BARHIT",DIRC_TRUTH_BARHIT); |
24 | |
25 | DIRC_TRUTH_PIXELTIME = false; |
26 | gPARMS->SetDefaultParameter("DIRC:TRUTH_PIXELTIME",DIRC_TRUTH_PIXELTIME); |
27 | |
28 | |
29 | DIRC_CUT_TDIFFD = 2; |
30 | gPARMS->SetDefaultParameter("DIRC:CUT_TDIFFD",DIRC_CUT_TDIFFD); |
31 | DIRC_CUT_TDIFFR = 3; |
32 | gPARMS->SetDefaultParameter("DIRC:CUT_TDIFFR",DIRC_CUT_TDIFFR); |
33 | |
34 | |
35 | DIRC_LIGHT_V = 19.65; |
36 | gPARMS->SetDefaultParameter("DIRC:LIGHT_V",DIRC_LIGHT_V); |
37 | |
38 | |
39 | DIRC_SIGMA_THETAC = 0.008; |
40 | gPARMS->SetDefaultParameter("DIRC:SIGMA_THETAC",DIRC_SIGMA_THETAC); |
41 | |
42 | |
43 | DIRC_ROTATE_TRACK = true; |
44 | gPARMS->SetDefaultParameter("DIRC:ROTATE_TRACK",DIRC_ROTATE_TRACK); |
45 | |
46 | |
47 | DIRC_THETAC_OFFSET = false; |
48 | gPARMS->SetDefaultParameter("DIRC:THETAC_OFFSET",DIRC_THETAC_OFFSET); |
49 | |
50 | |
51 | DIRC_LUT_CORR = true; |
52 | gPARMS->SetDefaultParameter("DIRC:LUT_CORR",DIRC_LUT_CORR); |
53 | |
54 | |
55 | DIRC_CHROMATIC_CORR = true; |
56 | gPARMS->SetDefaultParameter("DIRC:CHROMATIC_CORR",DIRC_CHROMATIC_CORR); |
57 | DIRC_CHROMATIC_CONST = 0.0025; |
58 | gPARMS->SetDefaultParameter("DIRC:CHROMATIC_CONST",DIRC_CHROMATIC_CONST); |
59 | |
60 | |
61 | DIRC_CUT_ANGLE = 0.03; |
62 | gPARMS->SetDefaultParameter("DIRC:CUT_ANGLE",DIRC_CUT_ANGLE); |
63 | |
64 | |
65 | dFinalStatePIDs.push_back(Positron); |
66 | dFinalStatePIDs.push_back(PiPlus); |
67 | dFinalStatePIDs.push_back(KPlus); |
68 | dFinalStatePIDs.push_back(Proton); |
69 | |
70 | dMaxChannels = DDIRCGeometry::kPMTs*DDIRCGeometry::kPixels; |
71 | |
72 | dCriticalAngle = asin(1.00028/1.47125); |
73 | dIndex = 1.473; |
74 | |
75 | if(DIRC_DEBUG_HISTS) |
76 | CreateDebugHistograms(); |
77 | } |
78 | |
79 | bool DDIRCLut::brun(JEventLoop *loop) { |
80 | |
81 | dapp = dynamic_cast<DApplication*>(loop->GetJApplication()); |
82 | dDIRCLutReader = dapp->GetDIRCLut(loop->GetJEvent().GetRunNumber()); |
83 | |
84 | |
85 | vector<const DDIRCGeometry*> locDIRCGeometry; |
86 | loop->Get(locDIRCGeometry); |
87 | dDIRCGeometry = locDIRCGeometry[0]; |
88 | |
89 | |
90 | vector <double> thetac_offsets_bar(5184); |
91 | if(loop->GetCalib("/DIRC/thetac_offsets_bar", thetac_offsets_bar)) |
92 | jout << "Can't find requested /DIRC/thetac_offsets_bar in CCDB for this run!" << endl; |
93 | |
94 | for(int ibar=0; ibar<DDIRCGeometry::kBars; ibar++) { |
95 | for(int ipmt=0; ipmt<DDIRCGeometry::kPMTs; ipmt++) { |
96 | int index = ipmt + ibar*DDIRCGeometry::kPMTs; |
97 | dThetaCOffset[ibar][ipmt] = thetac_offsets_bar.at(index); |
98 | } |
99 | } |
100 | |
101 | |
102 | vector< map<string, double> > bar_rotation(DDIRCGeometry::kBars); |
103 | if(loop->GetCalib("/DIRC/bar_rotation", bar_rotation)) |
104 | jout << "Can't find requested /DIRC/bar_rotation in CCDB for this run!" << endl; |
105 | |
106 | for(int ibar=0; ibar<DDIRCGeometry::kBars; ibar++) { |
107 | dRotationX[ibar] = bar_rotation[ibar].at("rotationX")*TMath::DegToRad(); |
108 | dRotationY[ibar] = bar_rotation[ibar].at("rotationY")*TMath::DegToRad(); |
109 | dRotationZ[ibar] = bar_rotation[ibar].at("rotationZ")*TMath::DegToRad(); |
110 | } |
111 | |
112 | return true; |
113 | } |
114 | |
115 | bool DDIRCLut::CreateDebugHistograms() { |
116 | |
117 | |
118 | |
119 | |
120 | |
121 | |
122 | { |
123 | |
124 | |
125 | |
126 | |
127 | hDiff = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiff"); |
128 | hDiff_Pixel[0] = (TH2I*)gROOT(ROOT::GetROOT())->FindObject("hDiff_Pixel_N"); |
129 | hDiff_Pixel[1] = (TH2I*)gROOT(ROOT::GetROOT())->FindObject("hDiff_Pixel_S"); |
130 | hDiffT = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiffT"); |
131 | hDiffD = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiffD"); |
132 | hDiffR = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiffR"); |
133 | hTime = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hTime"); |
134 | hCalc = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hCalc"); |
135 | hNph = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hNph"); |
136 | if(!hDiff) hDiff = new TH1I("hDiff",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
137 | if(!hDiff_Pixel[0]) hDiff_Pixel[0] = new TH2I("hDiff_Pixel_N","; Channel ID; t_{calc}-t_{measured} [ns];entries [#]", dMaxChannels, 0, dMaxChannels, 400,-20,20); |
138 | if(!hDiff_Pixel[1]) hDiff_Pixel[1] = new TH2I("hDiff_Pixel_S","; Channel ID; t_{calc}-t_{measured} [ns];entries [#]", dMaxChannels, 0, dMaxChannels, 400,-20,20); |
139 | if(!hDiffT) hDiffT = new TH1I("hDiffT",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
140 | if(!hDiffD) hDiffD = new TH1I("hDiffD",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
141 | if(!hDiffR) hDiffR = new TH1I("hDiffR",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
142 | if(!hTime) hTime = new TH1I("hTime",";propagation time [ns];entries [#]", 1000,0,200); |
143 | if(!hCalc) hCalc = new TH1I("hCalc",";calculated time [ns];entries [#]", 1000,0,200); |
144 | if(!hNph) hNph = new TH1I("hNph",";detected photons [#];entries [#]", 150,0,150); |
145 | |
146 | |
147 | for(uint loc_i=0; loc_i<dFinalStatePIDs.size(); loc_i++) { |
148 | Particle_t locPID = dFinalStatePIDs[loc_i]; |
149 | string locParticleName = ParticleType(locPID); |
150 | string locParticleROOTName = ParticleName_ROOT(locPID); |
151 | |
152 | hDeltaThetaC[locPID] = (TH1I*)gROOT(ROOT::GetROOT())->FindObject(Form("hDeltaThetaC_%s",locParticleName.data())); |
153 | if(!hDeltaThetaC[locPID]) |
154 | hDeltaThetaC[locPID] = new TH1I(Form("hDeltaThetaC_%s",locParticleName.data()), "cherenkov angle; #Delta#theta_{C} [rad]", 200,-0.5,0.5); |
155 | hDeltaThetaC_Pixel[locPID] = (TH2I*)gROOT(ROOT::GetROOT())->FindObject(Form("hDeltaThetaC_Pixel_%s",locParticleName.data())); |
156 | if(!hDeltaThetaC_Pixel[locPID]) |
157 | hDeltaThetaC_Pixel[locPID] = new TH2I(Form("hDeltaThetaC_Pixel_%s",locParticleName.data()), "cherenkov angle; Pixel ID; #Delta#theta_{C} [rad]", 10000, 0, 10000, 200,-0.5,0.5); |
158 | } |
159 | |
160 | |
161 | } |
162 | |
163 | |
164 | return true; |
165 | } |
166 | |
167 | bool DDIRCLut::CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector<const DDIRCPmtHit*> locDIRCHits, double locFlightTime, double locMass, shared_ptr<DDIRCMatchParams>& locDIRCMatchParams, const vector<const DDIRCTruthBarHit*> locDIRCBarHits, map<shared_ptr<const DDIRCMatchParams>, vector<const DDIRCPmtHit*> >& locDIRCTrackMatchParams) const |
168 | { |
169 | |
170 | TVector3 momInBar = locProjMom; |
171 | TVector3 posInBar = locProjPos; |
172 | |
173 | |
174 | |
175 | |
176 | if(DIRC_TRUTH_BARHIT && locDIRCBarHits.size() > 0) { |
177 | |
178 | TVector3 bestMatchPos, bestMatchMom; |
179 | double bestFlightTime = 999.; |
180 | double bestMatchDist = 999.; |
181 | for(int i=0; i<(int)locDIRCBarHits.size(); i++) { |
182 | TVector3 locDIRCBarHitPos(locDIRCBarHits[i]->x, locDIRCBarHits[i]->y, locDIRCBarHits[i]->z); |
183 | TVector3 locDIRCBarHitMom(locDIRCBarHits[i]->px, locDIRCBarHits[i]->py, locDIRCBarHits[i]->pz); |
184 | if((posInBar - locDIRCBarHitPos).Mag() < bestMatchDist) { |
185 | bestMatchDist = (posInBar - locDIRCBarHitPos).Mag(); |
186 | bestMatchPos = locDIRCBarHitPos; |
187 | bestMatchMom = locDIRCBarHitMom; |
188 | bestFlightTime = locDIRCBarHits[i]->t; |
189 | } |
190 | } |
191 | |
192 | momInBar = bestMatchMom; |
193 | posInBar = bestMatchPos; |
194 | locFlightTime = bestFlightTime; |
195 | } |
196 | |
197 | |
198 | if(locDIRCMatchParams == nullptr) |
199 | locDIRCMatchParams = std::make_shared<DDIRCMatchParams>(); |
200 | |
201 | |
202 | double locAngle = CalcAngle(momInBar.Mag(), locMass); |
203 | map<Particle_t, double> locExpectedAngle = CalcExpectedAngles(momInBar.Mag()); |
204 | map<Particle_t, double> logLikelihoodSum; |
205 | Particle_t locHypothesisPID = PiPlus; |
206 | for(uint loc_i = 0; loc_i<dFinalStatePIDs.size(); loc_i++) { |
207 | logLikelihoodSum[dFinalStatePIDs[loc_i]]=0; |
208 | if(fabs(ParticleMass(dFinalStatePIDs[loc_i])-locMass) < 0.01) locHypothesisPID = dFinalStatePIDs[loc_i]; |
209 | } |
210 | |
211 | |
212 | int nPhotons = 0; |
213 | int nPhotonsThetaC = 0; |
214 | double meanThetaC = 0.; |
215 | double meanDeltaT = 0.; |
216 | for (unsigned int loc_i = 0; loc_i < locDIRCHits.size(); loc_i++){ |
217 | |
218 | const DDIRCPmtHit* locDIRCHit = locDIRCHits[loc_i]; |
219 | bool locIsGood = false; |
220 | bool locIsReflected = false; |
221 | CalcPhoton(locDIRCHit, locFlightTime, posInBar, momInBar, locExpectedAngle, locAngle, locHypothesisPID, locIsReflected, logLikelihoodSum, nPhotonsThetaC, meanThetaC, meanDeltaT, locIsGood); |
222 | if(locIsGood) { |
223 | |
224 | nPhotons++; |
225 | locDIRCTrackMatchParams[locDIRCMatchParams].push_back(locDIRCHit); |
226 | } |
227 | } |
228 | |
229 | if(DIRC_DEBUG_HISTS) { |
230 | dapp->RootWriteLock(); |
231 | hNph->Fill(nPhotons); |
232 | dapp->RootUnLock(); |
233 | } |
234 | |
235 | |
236 | if(nPhotons<5) |
237 | return false; |
238 | |
239 | |
240 | locDIRCMatchParams->dThetaC = meanThetaC/(double)nPhotonsThetaC; |
241 | locDIRCMatchParams->dDeltaT = meanDeltaT/(double)nPhotonsThetaC; |
242 | locDIRCMatchParams->dExpectedThetaC = locAngle; |
243 | locDIRCMatchParams->dLikelihoodElectron = logLikelihoodSum[Positron]; |
244 | locDIRCMatchParams->dLikelihoodPion = logLikelihoodSum[PiPlus]; |
245 | locDIRCMatchParams->dLikelihoodKaon = logLikelihoodSum[KPlus]; |
246 | locDIRCMatchParams->dLikelihoodProton = logLikelihoodSum[Proton]; |
247 | locDIRCMatchParams->dNPhotons = nPhotons; |
248 | locDIRCMatchParams->dExtrapolatedPos = posInBar; |
249 | locDIRCMatchParams->dExtrapolatedMom = momInBar; |
250 | locDIRCMatchParams->dExtrapolatedTime = locFlightTime; |
251 | |
252 | return true; |
253 | } |
254 | |
255 | vector<pair<double,double>> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, double locFlightTime, TVector3 posInBar, TVector3 momInBar, map<Particle_t, double> locExpectedAngle, double locAngle, Particle_t locPID, bool &isReflected, map<Particle_t, double> &logLikelihoodSum, int &nPhotonsThetaC, double &meanThetaC, double &meanDeltaT, bool &isGood) const |
256 | { |
257 | |
258 | pair<double, double> locDIRCPhoton(-999., -999.); |
259 | vector<pair<double, double>> locDIRCPhotons; |
260 | |
261 | TVector3 fnX1 = TVector3 (1,0,0); |
262 | TVector3 fnY1 = TVector3 (0,1,0); |
263 | TVector3 fnZ1 = TVector3 (0,0,1); |
264 | |
265 | double tangle,luttheta,evtime; |
266 | int64_t pathid; |
267 | TVector3 dir,dird; |
268 | |
269 | |
270 | int bar = dDIRCGeometry->GetBar(posInBar.Y()); |
271 | if(bar < 0 || bar > 47) return locDIRCPhotons; |
272 | |
273 | |
274 | int channel = locDIRCHit->ch; |
275 | |
276 | |
277 | int box = (bar < 24) ? 1 : 0; |
278 | if((box == 0 && channel < dMaxChannels) || (box == 1 && channel >= dMaxChannels)) |
279 | return locDIRCPhotons; |
280 | |
281 | int box_channel = channel%dMaxChannels; |
282 | int box_pmt = box_channel/DDIRCGeometry::kPixels; |
283 | |
284 | double rotX = 0, rotY = 0, rotZ=0; |
285 | uint xbinSize = 200.0 / dDIRCLutReader->GetLutCorrNbins(); |
286 | int bin = (int)( (posInBar.X()+100.) / xbinSize ); |
287 | if(bar>=0 && bar<=47 && bin>=0 && bin<=39) { |
288 | rotX = dRotationX[bar]; |
289 | rotY = dRotationY[bar]; |
290 | rotZ = dRotationZ[bar]; |
291 | } |
292 | |
293 | |
294 | double hitTime = locDIRCHit->t - locFlightTime; |
295 | |
296 | |
297 | vector<const DDIRCTruthPmtHit*> locTruthDIRCHits; |
298 | locDIRCHit->Get(locTruthDIRCHits); |
299 | if(DIRC_TRUTH_PIXELTIME && locTruthDIRCHits.size() > 0) { |
300 | double locTruthTime = locTruthDIRCHits[0]->t - locFlightTime; |
301 | hitTime = locTruthTime; |
302 | } |
303 | |
304 | |
305 | bool reflected = hitTime>35; |
306 | |
307 | |
308 | double radiatorL = dDIRCGeometry->GetBarLength(bar); |
309 | double barend = dDIRCGeometry->GetBarEnd(bar); |
310 | double lenz = fabs(barend - posInBar.X()); |
311 | |
312 | |
313 | double rlenz = 2*radiatorL - lenz; |
314 | double dlenz = lenz; |
315 | if(reflected) lenz = 2*radiatorL - lenz; |
| Value stored to 'lenz' is never read |
316 | |
317 | |
318 | if(dDIRCLutReader->GetLutPixelAngleSize(bar, box_channel) == 0) |
319 | return locDIRCPhotons; |
320 | |
321 | |
322 | for(uint i = 0; i < dDIRCLutReader->GetLutPixelAngleSize(bar, box_channel); i++){ |
323 | |
324 | dird = dDIRCLutReader->GetLutPixelAngle(bar, box_channel, i); |
325 | evtime = dDIRCLutReader->GetLutPixelTime(bar, box_channel, i); |
326 | pathid = dDIRCLutReader->GetLutPixelPath(bar, box_channel, i); |
327 | |
328 | |
329 | bool samepath(false); |
330 | if(!locTruthDIRCHits.empty() && fabs(pathid - locTruthDIRCHits[0]->path)<0.0001) |
331 | samepath=true; |
332 | |
333 | for(int r=0; r<2; r++){ |
334 | if(!reflected && r==1) continue; |
335 | |
336 | if(r) lenz = rlenz; |
337 | else lenz = dlenz; |
338 | |
339 | for(int u = 0; u < 4; u++){ |
340 | if(u == 0) dir = dird; |
341 | if(u == 1) dir.SetXYZ( dird.X(),-dird.Y(), dird.Z()); |
342 | if(u == 2) dir.SetXYZ( dird.X(), dird.Y(), -dird.Z()); |
343 | if(u == 3) dir.SetXYZ( dird.X(),-dird.Y(), -dird.Z()); |
344 | if(r) dir.SetXYZ( -dir.X(), dir.Y(), dir.Z()); |
345 | if(dir.Angle(fnY1) < dCriticalAngle || dir.Angle(fnZ1) < dCriticalAngle) continue; |
346 | |
347 | dir = dir.Unit(); |
348 | |
349 | luttheta = dir.Angle(TVector3(-1,0,0)); |
350 | if(luttheta > TMath::PiOver2()) luttheta = TMath::Pi()-luttheta; |
351 | |
352 | double bartime = lenz/cos(luttheta)/DIRC_LIGHT_V; |
353 | double totalTime = bartime+evtime; |
354 | |
355 | |
356 | if(DIRC_LUT_CORR){ |
357 | if (reflected) totalTime -= dDIRCLutReader->GetLutCorrTimeReflected(bar,box_pmt,bin); |
358 | else totalTime -= dDIRCLutReader->GetLutCorrTimeDirect(bar,box_pmt,bin); |
359 | } |
360 | |
361 | |
362 | double locDeltaT = totalTime-hitTime; |
363 | |
364 | TVector3 trackMom = momInBar; |
365 | if(DIRC_ROTATE_TRACK) { |
366 | trackMom.RotateX(rotX); |
367 | trackMom.RotateY(rotY); |
368 | trackMom.RotateZ(rotZ); |
369 | } |
370 | tangle = trackMom.Angle(dir); |
371 | |
372 | if(DIRC_THETAC_OFFSET) { |
373 | tangle -= dThetaCOffset[bar][box_pmt]; |
374 | } |
375 | |
376 | |
377 | if(DIRC_LUT_CORR){ |
378 | if (r) tangle += dDIRCLutReader->GetLutCorrAngleReflected(bar,box_pmt,bin); |
379 | else tangle += dDIRCLutReader->GetLutCorrAngleDirect(bar,box_pmt,bin); |
380 | } |
381 | |
382 | |
383 | if (fabs(locDeltaT) < 2.0 && DIRC_CHROMATIC_CORR) |
384 | tangle += DIRC_CHROMATIC_CONST * locDeltaT; |
385 | |
386 | if(DIRC_DEBUG_HISTS) { |
387 | dapp->RootWriteLock(); |
388 | hTime->Fill(hitTime); |
389 | hCalc->Fill(totalTime); |
390 | |
391 | if(fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))<0.2){ |
392 | hDiff->Fill(locDeltaT); |
393 | hDiff_Pixel[box]->Fill(channel%dMaxChannels, locDeltaT); |
394 | if(samepath){ |
395 | hDiffT->Fill(locDeltaT); |
396 | if(r) hDiffR->Fill(locDeltaT); |
397 | else hDiffD->Fill(locDeltaT); |
398 | } |
399 | } |
400 | dapp->RootUnLock(); |
401 | } |
402 | |
403 | |
404 | if(fabs(locDeltaT) < 100.0 && fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))<0.2) { |
405 | locDIRCPhoton.first = totalTime; |
406 | locDIRCPhoton.second = tangle; |
407 | locDIRCPhotons.push_back(locDIRCPhoton); |
408 | } |
409 | |
410 | |
411 | if(!r && fabs(locDeltaT)>DIRC_CUT_TDIFFD) continue; |
412 | if( r && fabs(locDeltaT)>DIRC_CUT_TDIFFR) continue; |
413 | |
414 | if(DIRC_DEBUG_HISTS) { |
415 | dapp->RootWriteLock(); |
416 | |
417 | |
418 | dapp->RootUnLock(); |
419 | } |
420 | |
421 | |
422 | if(locPID==Proton || locPID==AntiProton) { |
423 | if(fabs(tangle-locExpectedAngle[PiPlus])>DIRC_CUT_ANGLE && fabs(tangle-locExpectedAngle[KPlus])>DIRC_CUT_ANGLE && fabs(tangle-locExpectedAngle[Proton])>DIRC_CUT_ANGLE) continue; |
424 | } |
425 | else if(fabs(tangle-locExpectedAngle[PiPlus])>DIRC_CUT_ANGLE && fabs(tangle-locExpectedAngle[KPlus])>DIRC_CUT_ANGLE) continue; |
426 | |
427 | isReflected = r; |
428 | |
429 | |
430 | isGood = true; |
431 | |
432 | |
433 | nPhotonsThetaC++; |
434 | |
435 | |
436 | meanThetaC += tangle; |
437 | meanDeltaT += locDeltaT; |
438 | |
439 | |
440 | |
441 | for(uint loc_j = 0; loc_j<dFinalStatePIDs.size(); loc_j++) { |
442 | logLikelihoodSum[dFinalStatePIDs[loc_j]] += TMath::Log( CalcLikelihood(locExpectedAngle[dFinalStatePIDs[loc_j]], tangle)); |
443 | } |
444 | |
445 | } |
446 | } |
447 | } |
448 | |
449 | return locDIRCPhotons; |
450 | } |
451 | |
452 | |
453 | vector<pair<double,double>> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, double locFlightTime, TVector3 posInBar, TVector3 momInBar, map<Particle_t, double> locExpectedAngle, double locAngle, Particle_t locPID, bool &isReflected, map<Particle_t, double> &logLikelihoodSum) const |
454 | { |
455 | int nPhotonsThetaC=0; |
456 | double meanThetaC=0.0, meanDeltaT=0.0; |
457 | bool isGood=false; |
458 | return CalcPhoton(locDIRCHit, locFlightTime, posInBar, momInBar, locExpectedAngle, locAngle, locPID, isReflected, logLikelihoodSum, nPhotonsThetaC, meanThetaC, meanDeltaT, isGood); |
459 | } |
460 | |
461 | double DDIRCLut::CalcLikelihood(double locExpectedThetaC, double locThetaC) const { |
462 | |
463 | double locLikelihood = TMath::Exp(-0.5*( (locExpectedThetaC-locThetaC)/DIRC_SIGMA_THETAC * (locExpectedThetaC-locThetaC)/DIRC_SIGMA_THETAC ) ) + 0.00001; |
464 | |
465 | return locLikelihood; |
466 | } |
467 | |
468 | double DDIRCLut::CalcAngle(double locP, double locMass) const { |
469 | return acos(sqrt(locP*locP + locMass*locMass)/locP/dIndex); |
470 | } |
471 | |
472 | map<Particle_t, double> DDIRCLut::CalcExpectedAngles(double locP) const { |
473 | |
474 | map<Particle_t, double> locExpectedAngles; |
475 | for(uint loc_i = 0; loc_i<dFinalStatePIDs.size(); loc_i++) { |
476 | locExpectedAngles[dFinalStatePIDs[loc_i]] = acos(sqrt(locP*locP + ParticleMass(dFinalStatePIDs[loc_i])*ParticleMass(dFinalStatePIDs[loc_i]))/locP/dIndex); |
477 | } |
478 | return locExpectedAngles; |
479 | } |