File: | /home/sdobbs/work/clang/halld_recon/src/libraries/DIRC/DDIRCLut.cc |
Warning: | line 297, column 16 Value stored to 'lenz' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
1 | // $Id$ |
2 | // |
3 | // File: DDIRCLut.cc |
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 | // DDIRCLut (Constructor) |
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 | // timing cuts for photons |
29 | DIRC_CUT_TDIFFD = 2; // direct cut in ns |
30 | gPARMS->SetDefaultParameter("DIRC:CUT_TDIFFD",DIRC_CUT_TDIFFD); |
31 | DIRC_CUT_TDIFFR = 3; // reflected cut in ns |
32 | gPARMS->SetDefaultParameter("DIRC:CUT_TDIFFR",DIRC_CUT_TDIFFR); |
33 | |
34 | // Gives DeltaT = 0, but shouldn't it be v=20.3767 [cm/ns] for 1.47125 |
35 | DIRC_LIGHT_V = 19.80; // v=19.8 [cm/ns] for 1.5141 |
36 | gPARMS->SetDefaultParameter("DIRC:LIGHT_V",DIRC_LIGHT_V); |
37 | |
38 | // sigma (thetaC for single photon) in radiams |
39 | DIRC_SIGMA_THETAC = 0.01; |
40 | gPARMS->SetDefaultParameter("DIRC:SIGMA_THETAC",DIRC_SIGMA_THETAC); |
41 | |
42 | // Rotate tracks angle based on bar box survey data |
43 | DIRC_ROTATE_TRACK = false; |
44 | gPARMS->SetDefaultParameter("DIRC:ROTATE_TRACK",DIRC_ROTATE_TRACK); |
45 | |
46 | // Rotate tracks angle based on bar box survey data |
47 | DIRC_THETAC_OFFSET = true; |
48 | gPARMS->SetDefaultParameter("DIRC:THETAC_OFFSET",DIRC_THETAC_OFFSET); |
49 | |
50 | // set PID for different passes in debuging histograms |
51 | dFinalStatePIDs.push_back(Positron); |
52 | dFinalStatePIDs.push_back(PiPlus); |
53 | dFinalStatePIDs.push_back(KPlus); |
54 | dFinalStatePIDs.push_back(Proton); |
55 | |
56 | dMaxChannels = DDIRCGeometry::kPMTs*DDIRCGeometry::kPixels; |
57 | |
58 | dCriticalAngle = asin(1.00028/1.47125); // n_quarzt = 1.47125; //(1.47125 <==> 390nm) |
59 | dIndex = 1.473; |
60 | |
61 | if(DIRC_DEBUG_HISTS) |
62 | CreateDebugHistograms(); |
63 | } |
64 | |
65 | bool DDIRCLut::brun(JEventLoop *loop) { |
66 | |
67 | dapp = dynamic_cast<DApplication*>(loop->GetJApplication()); |
68 | dDIRCLutReader = dapp->GetDIRCLut(loop->GetJEvent().GetRunNumber()); |
69 | |
70 | // get DIRC geometry |
71 | vector<const DDIRCGeometry*> locDIRCGeometry; |
72 | loop->Get(locDIRCGeometry); |
73 | dDIRCGeometry = locDIRCGeometry[0]; |
74 | |
75 | // get cherenkov angle corrections from CCDB |
76 | vector <double> thetac_offsets_bar(5184); |
77 | if(loop->GetCalib("/DIRC/thetac_offsets_bar", thetac_offsets_bar)) |
78 | jout << "Can't find requested /DIRC/thetac_offsets_bar in CCDB for this run!" << endl; |
79 | |
80 | for(int ibar=0; ibar<DDIRCGeometry::kBars; ibar++) { |
81 | for(int ipmt=0; ipmt<DDIRCGeometry::kPMTs; ipmt++) { |
82 | int index = ipmt + ibar*DDIRCGeometry::kPMTs; |
83 | dThetaCOffset[ibar][ipmt] = thetac_offsets_bar.at(index); |
84 | } |
85 | } |
86 | |
87 | // get track rotation corrections from CCDB |
88 | vector< map<string, double> > bar_rotation(DDIRCGeometry::kBars); |
89 | if(loop->GetCalib("/DIRC/bar_rotation", bar_rotation)) |
90 | jout << "Can't find requested /DIRC/bar_rotation in CCDB for this run!" << endl; |
91 | |
92 | for(int ibar=0; ibar<DDIRCGeometry::kBars; ibar++) { |
93 | dRotationX[ibar] = bar_rotation[ibar].at("rotationX")*TMath::DegToRad(); |
94 | dRotationY[ibar] = bar_rotation[ibar].at("rotationY")*TMath::DegToRad(); |
95 | dRotationZ[ibar] = bar_rotation[ibar].at("rotationZ")*TMath::DegToRad(); |
96 | } |
97 | |
98 | return true; |
99 | } |
100 | |
101 | bool DDIRCLut::CreateDebugHistograms() { |
102 | |
103 | //////////////////////////////////// |
104 | // LUT histograms and diagnostics // |
105 | //////////////////////////////////// |
106 | |
107 | //dapp->RootWriteLock(); //ACQUIRE ROOT LOCK!! |
108 | { |
109 | //TDirectory *mainDir = gDirectory; |
110 | //TDirectory *dircDir = gDirectory->mkdir("DIRC_debug"); |
111 | //dircDir->cd(); |
112 | |
113 | hDiff = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiff"); |
114 | hDiff_Pixel[0] = (TH2I*)gROOT(ROOT::GetROOT())->FindObject("hDiff_Pixel_N"); |
115 | hDiff_Pixel[1] = (TH2I*)gROOT(ROOT::GetROOT())->FindObject("hDiff_Pixel_S"); |
116 | hDiffT = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiffT"); |
117 | hDiffD = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiffD"); |
118 | hDiffR = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hDiffR"); |
119 | hTime = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hTime"); |
120 | hCalc = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hCalc"); |
121 | hNph = (TH1I*)gROOT(ROOT::GetROOT())->FindObject("hNph"); |
122 | if(!hDiff) hDiff = new TH1I("hDiff",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
123 | 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); |
124 | 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); |
125 | if(!hDiffT) hDiffT = new TH1I("hDiffT",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
126 | if(!hDiffD) hDiffD = new TH1I("hDiffD",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
127 | if(!hDiffR) hDiffR = new TH1I("hDiffR",";t_{calc}-t_{measured} [ns];entries [#]", 400,-20,20); |
128 | if(!hTime) hTime = new TH1I("hTime",";propagation time [ns];entries [#]", 1000,0,200); |
129 | if(!hCalc) hCalc = new TH1I("hCalc",";calculated time [ns];entries [#]", 1000,0,200); |
130 | if(!hNph) hNph = new TH1I("hNph",";detected photons [#];entries [#]", 150,0,150); |
131 | |
132 | // DeltaThetaC for each particle type (e, pi, K, p) and per pixel |
133 | for(uint loc_i=0; loc_i<dFinalStatePIDs.size(); loc_i++) { |
134 | Particle_t locPID = dFinalStatePIDs[loc_i]; |
135 | string locParticleName = ParticleType(locPID); |
136 | string locParticleROOTName = ParticleName_ROOT(locPID); |
137 | |
138 | hDeltaThetaC[locPID] = (TH1I*)gROOT(ROOT::GetROOT())->FindObject(Form("hDeltaThetaC_%s",locParticleName.data())); |
139 | if(!hDeltaThetaC[locPID]) |
140 | hDeltaThetaC[locPID] = new TH1I(Form("hDeltaThetaC_%s",locParticleName.data()), "cherenkov angle; #Delta#theta_{C} [rad]", 200,-0.5,0.5); |
141 | hDeltaThetaC_Pixel[locPID] = (TH2I*)gROOT(ROOT::GetROOT())->FindObject(Form("hDeltaThetaC_Pixel_%s",locParticleName.data())); |
142 | if(!hDeltaThetaC_Pixel[locPID]) |
143 | 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); |
144 | } |
145 | |
146 | //mainDir->cd(); |
147 | } |
148 | //dapp->RootUnLock(); //REMOVE ROOT LOCK!! |
149 | |
150 | return true; |
151 | } |
152 | |
153 | 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 |
154 | { |
155 | // get bar and track position/momentum from extrapolation |
156 | TVector3 momInBar = locProjMom; |
157 | TVector3 posInBar = locProjPos; |
158 | |
159 | //////////////////////////////////////////// |
160 | // option to cheat and use truth position // |
161 | //////////////////////////////////////////// |
162 | if(DIRC_TRUTH_BARHIT && locDIRCBarHits.size() > 0) { |
163 | |
164 | TVector3 bestMatchPos, bestMatchMom; |
165 | double bestFlightTime = 999.; |
166 | double bestMatchDist = 999.; |
167 | for(int i=0; i<(int)locDIRCBarHits.size(); i++) { |
168 | TVector3 locDIRCBarHitPos(locDIRCBarHits[i]->x, locDIRCBarHits[i]->y, locDIRCBarHits[i]->z); |
169 | TVector3 locDIRCBarHitMom(locDIRCBarHits[i]->px, locDIRCBarHits[i]->py, locDIRCBarHits[i]->pz); |
170 | if((posInBar - locDIRCBarHitPos).Mag() < bestMatchDist) { |
171 | bestMatchDist = (posInBar - locDIRCBarHitPos).Mag(); |
172 | bestMatchPos = locDIRCBarHitPos; |
173 | bestMatchMom = locDIRCBarHitMom; |
174 | bestFlightTime = locDIRCBarHits[i]->t; |
175 | } |
176 | } |
177 | |
178 | momInBar = bestMatchMom; |
179 | posInBar = bestMatchPos; |
180 | locFlightTime = bestFlightTime; |
181 | } |
182 | |
183 | // clear object to store good photons |
184 | if(locDIRCMatchParams == nullptr) |
185 | locDIRCMatchParams = std::make_shared<DDIRCMatchParams>(); |
186 | |
187 | // initialize variables for LUT summary |
188 | double locAngle = CalcAngle(momInBar.Mag(), locMass); |
189 | map<Particle_t, double> locExpectedAngle = CalcExpectedAngles(momInBar.Mag()); |
190 | map<Particle_t, double> logLikelihoodSum; |
191 | Particle_t locHypothesisPID = PiPlus; |
192 | for(uint loc_i = 0; loc_i<dFinalStatePIDs.size(); loc_i++) { |
193 | logLikelihoodSum[dFinalStatePIDs[loc_i]]=0; |
194 | if(fabs(ParticleMass(dFinalStatePIDs[loc_i])-locMass) < 0.01) locHypothesisPID = dFinalStatePIDs[loc_i]; |
195 | } |
196 | |
197 | // loop over DIRC hits |
198 | int nPhotons = 0; |
199 | int nPhotonsThetaC = 0; |
200 | double meanThetaC = 0.; |
201 | double meanDeltaT = 0.; |
202 | for (unsigned int loc_i = 0; loc_i < locDIRCHits.size(); loc_i++){ |
203 | |
204 | const DDIRCPmtHit* locDIRCHit = locDIRCHits[loc_i]; |
205 | bool locIsGood = false; |
206 | bool locIsReflected = false; |
207 | CalcPhoton(locDIRCHit, locFlightTime, posInBar, momInBar, locExpectedAngle, locAngle, locHypothesisPID, locIsReflected, logLikelihoodSum, nPhotonsThetaC, meanThetaC, meanDeltaT, locIsGood); |
208 | if(locIsGood) { |
209 | // count good photons and add hits to associated objects |
210 | nPhotons++; |
211 | locDIRCTrackMatchParams[locDIRCMatchParams].push_back(locDIRCHit); |
212 | } |
213 | }// end loop over hits |
214 | |
215 | if(DIRC_DEBUG_HISTS) { |
216 | dapp->RootWriteLock(); |
217 | hNph->Fill(nPhotons); |
218 | dapp->RootUnLock(); |
219 | } |
220 | |
221 | // skip tracks without enough photons |
222 | if(nPhotons<5) |
223 | return false; |
224 | |
225 | // set DIRCMatchParameters contents |
226 | locDIRCMatchParams->dThetaC = meanThetaC/(double)nPhotonsThetaC; |
227 | locDIRCMatchParams->dDeltaT = meanDeltaT/(double)nPhotonsThetaC; |
228 | locDIRCMatchParams->dExpectedThetaC = locAngle; |
229 | locDIRCMatchParams->dLikelihoodElectron = logLikelihoodSum[Positron]; |
230 | locDIRCMatchParams->dLikelihoodPion = logLikelihoodSum[PiPlus]; |
231 | locDIRCMatchParams->dLikelihoodKaon = logLikelihoodSum[KPlus]; |
232 | locDIRCMatchParams->dLikelihoodProton = logLikelihoodSum[Proton]; |
233 | locDIRCMatchParams->dNPhotons = nPhotons; |
234 | locDIRCMatchParams->dExtrapolatedPos = posInBar; |
235 | locDIRCMatchParams->dExtrapolatedMom = momInBar; |
236 | locDIRCMatchParams->dExtrapolatedTime = locFlightTime; |
237 | |
238 | return true; |
239 | } |
240 | |
241 | 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 |
242 | { |
243 | // initialize photon pairs for time and thetaC |
244 | pair<double, double> locDIRCPhoton(-999., -999.); |
245 | vector<pair<double, double>> locDIRCPhotons; |
246 | |
247 | TVector3 fnX1 = TVector3 (1,0,0); |
248 | TVector3 fnY1 = TVector3 (0,1,0); |
249 | TVector3 fnZ1 = TVector3 (0,0,1); |
250 | |
251 | double tangle,luttheta,evtime; |
252 | int64_t pathid; |
253 | TVector3 dir,dird; |
254 | |
255 | // get bar number from geometry |
256 | int bar = dDIRCGeometry->GetBar(posInBar.Y()); |
257 | if(bar < 0 || bar > 47) return locDIRCPhotons; |
258 | |
259 | // get channel information for LUT |
260 | int channel = locDIRCHit->ch; |
261 | |
262 | // get box from bar number (North/Upper = 0 and South/Lower = 1) |
263 | int box = (bar < 24) ? 1 : 0; |
264 | if((box == 0 && channel < dMaxChannels) || (box == 1 && channel >= dMaxChannels)) |
265 | return locDIRCPhotons; |
266 | |
267 | double rotX = 0, rotY = 0, rotZ=0; |
268 | int xbin = (int)(posInBar.X()/5.0) + 19; |
269 | if(bar>=0 && bar<=47 && xbin>=0 && xbin<=39) { |
270 | rotX = dRotationX[bar]; |
271 | rotY = dRotationY[bar]; |
272 | rotZ = dRotationZ[bar]; |
273 | } |
274 | |
275 | // use hit time to determine if reflected or not |
276 | double hitTime = locDIRCHit->t - locFlightTime; |
277 | |
278 | // option to use truth hit time rather than smeared hit time |
279 | vector<const DDIRCTruthPmtHit*> locTruthDIRCHits; |
280 | locDIRCHit->Get(locTruthDIRCHits); |
281 | if(DIRC_TRUTH_PIXELTIME && locTruthDIRCHits.size() > 0) { |
282 | double locTruthTime = locTruthDIRCHits[0]->t - locFlightTime; |
283 | hitTime = locTruthTime; |
284 | } |
285 | |
286 | // needs to be X dependent choice for reflection cut (from CCDB?) |
287 | bool reflected = hitTime>35; // try only some photons as reflected for now |
288 | |
289 | // get position along bar for calculated time |
290 | double radiatorL = dDIRCGeometry->GetBarLength(bar); |
291 | double barend = dDIRCGeometry->GetBarEnd(bar); |
292 | double lenz = fabs(barend - posInBar.X()); |
293 | |
294 | // get length for reflected and direct photons |
295 | double rlenz = 2*radiatorL - lenz; // reflected |
296 | double dlenz = lenz; // direct |
297 | if(reflected) lenz = 2*radiatorL - lenz; |
Value stored to 'lenz' is never read | |
298 | |
299 | // check for pixel before going through loop |
300 | int box_channel = channel%dMaxChannels; |
301 | int box_pmt = box_channel/DDIRCGeometry::kPixels; |
302 | if(dDIRCLutReader->GetLutPixelAngleSize(bar, box_channel) == 0) |
303 | return locDIRCPhotons; |
304 | |
305 | // loop over LUT table for this bar/pixel to calculate thetaC |
306 | for(uint i = 0; i < dDIRCLutReader->GetLutPixelAngleSize(bar, box_channel); i++){ |
307 | |
308 | dird = dDIRCLutReader->GetLutPixelAngle(bar, box_channel, i); |
309 | evtime = dDIRCLutReader->GetLutPixelTime(bar, box_channel, i); |
310 | pathid = dDIRCLutReader->GetLutPixelPath(bar, box_channel, i); |
311 | |
312 | // in MC we can check if the path of the LUT and measured photon are the same |
313 | bool samepath(false); |
314 | if(!locTruthDIRCHits.empty() && fabs(pathid - locTruthDIRCHits[0]->path)<0.0001) |
315 | samepath=true; |
316 | |
317 | for(int r=0; r<2; r++){ |
318 | if(!reflected && r==1) continue; |
319 | |
320 | if(r) lenz = rlenz; |
321 | else lenz = dlenz; |
322 | |
323 | for(int u = 0; u < 4; u++){ |
324 | if(u == 0) dir = dird; |
325 | if(u == 1) dir.SetXYZ( dird.X(),-dird.Y(), dird.Z()); |
326 | if(u == 2) dir.SetXYZ( dird.X(), dird.Y(), -dird.Z()); |
327 | if(u == 3) dir.SetXYZ( dird.X(),-dird.Y(), -dird.Z()); |
328 | if(r) dir.SetXYZ( -dir.X(), dir.Y(), dir.Z()); |
329 | if(dir.Angle(fnY1) < dCriticalAngle || dir.Angle(fnZ1) < dCriticalAngle) continue; |
330 | |
331 | TVector3 trackMom = momInBar; |
332 | if(DIRC_ROTATE_TRACK) { // rotate tracks to bar plane from survey data |
333 | trackMom.RotateX(rotX); |
334 | trackMom.RotateY(rotY); |
335 | trackMom.RotateZ(rotZ); |
336 | } |
337 | tangle = trackMom.Angle(dir); |
338 | |
339 | if(DIRC_THETAC_OFFSET) { // ad-hoc correction per-PMT |
340 | tangle -= dThetaCOffset[bar][box_pmt]; |
341 | } |
342 | |
343 | luttheta = dir.Angle(TVector3(-1,0,0)); |
344 | if(luttheta > TMath::PiOver2()) luttheta = TMath::Pi()-luttheta; |
345 | double bartime = lenz/cos(luttheta)/DIRC_LIGHT_V; |
346 | double totalTime = bartime+evtime; |
347 | |
348 | // calculate time difference |
349 | double locDeltaT = totalTime-hitTime; |
350 | |
351 | if(DIRC_DEBUG_HISTS) { |
352 | dapp->RootWriteLock(); |
353 | hTime->Fill(hitTime); |
354 | hCalc->Fill(totalTime); |
355 | |
356 | if(fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))<0.2){ |
357 | hDiff->Fill(locDeltaT); |
358 | hDiff_Pixel[box]->Fill(channel%dMaxChannels, locDeltaT); |
359 | if(samepath){ |
360 | hDiffT->Fill(locDeltaT); |
361 | if(r) hDiffR->Fill(locDeltaT); |
362 | else hDiffD->Fill(locDeltaT); |
363 | } |
364 | } |
365 | dapp->RootUnLock(); |
366 | } |
367 | |
368 | // save hits array which pass some lose time and angle criteria |
369 | if(fabs(locDeltaT) < 100.0 && fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))<0.2) { |
370 | locDIRCPhoton.first = totalTime; |
371 | locDIRCPhoton.second = tangle; |
372 | locDIRCPhotons.push_back(locDIRCPhoton); |
373 | } |
374 | |
375 | // reject photons that are too far out of time |
376 | if(!r && fabs(locDeltaT)>DIRC_CUT_TDIFFD) continue; |
377 | if( r && fabs(locDeltaT)>DIRC_CUT_TDIFFR) continue; |
378 | |
379 | if(DIRC_DEBUG_HISTS) { |
380 | dapp->RootWriteLock(); |
381 | //hDeltaThetaC[locPID]->Fill(tangle-locAngle); |
382 | //hDeltaThetaC_Pixel[locPID]->Fill(channel, tangle-locAngle); |
383 | dapp->RootUnLock(); |
384 | } |
385 | |
386 | // remove photon candidates not used in likelihood |
387 | if(fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))>0.03) continue; |
388 | |
389 | isReflected = r; |
390 | |
391 | // save good photons to matched list |
392 | isGood = true; |
393 | |
394 | // count good photons |
395 | nPhotonsThetaC++; |
396 | |
397 | // calculate average ThetaC and DeltaT |
398 | meanThetaC += tangle; |
399 | meanDeltaT += locDeltaT; |
400 | |
401 | |
402 | // calculate likelihood for each mass hypothesis |
403 | for(uint loc_j = 0; loc_j<dFinalStatePIDs.size(); loc_j++) { |
404 | logLikelihoodSum[dFinalStatePIDs[loc_j]] += TMath::Log( CalcLikelihood(locExpectedAngle[dFinalStatePIDs[loc_j]], tangle)); |
405 | } |
406 | |
407 | } |
408 | } // end loop over reflections |
409 | } // end loop over nodes |
410 | |
411 | return locDIRCPhotons; |
412 | } |
413 | |
414 | // overloaded function when calculating outside LUT factory |
415 | 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 |
416 | { |
417 | int nPhotonsThetaC=0; |
418 | double meanThetaC=0.0, meanDeltaT=0.0; |
419 | bool isGood=false; |
420 | return CalcPhoton(locDIRCHit, locFlightTime, posInBar, momInBar, locExpectedAngle, locAngle, locPID, isReflected, logLikelihoodSum, nPhotonsThetaC, meanThetaC, meanDeltaT, isGood); |
421 | } |
422 | |
423 | double DDIRCLut::CalcLikelihood(double locExpectedThetaC, double locThetaC) const { |
424 | |
425 | double locLikelihood = TMath::Exp(-0.5*( (locExpectedThetaC-locThetaC)/DIRC_SIGMA_THETAC * (locExpectedThetaC-locThetaC)/DIRC_SIGMA_THETAC ) ) + 0.00001; |
426 | |
427 | return locLikelihood; |
428 | } |
429 | |
430 | double DDIRCLut::CalcAngle(double locP, double locMass) const { |
431 | return acos(sqrt(locP*locP + locMass*locMass)/locP/dIndex); |
432 | } |
433 | |
434 | map<Particle_t, double> DDIRCLut::CalcExpectedAngles(double locP) const { |
435 | |
436 | map<Particle_t, double> locExpectedAngles; |
437 | for(uint loc_i = 0; loc_i<dFinalStatePIDs.size(); loc_i++) { |
438 | locExpectedAngles[dFinalStatePIDs[loc_i]] = acos(sqrt(locP*locP + ParticleMass(dFinalStatePIDs[loc_i])*ParticleMass(dFinalStatePIDs[loc_i]))/locP/dIndex); |
439 | } |
440 | return locExpectedAngles; |
441 | } |