Bug Summary

File:/volatile/halld/gluex/nightly/2024-04-23/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/DIRC/DDIRCLut.cc
Location:line 315, column 16
Description:Value stored to 'lenz' is never read

Annotated Source Code

1// $Id$
2//
3// File: DDIRCLut.cc
4//
5
6#include <cassert>
7#include <math.h>
8using namespace std;
9
10#include "DDIRCLut.h"
11#include "DANA/DApplication.h"
12#include <JANA/JCalibration.h>
13
14//---------------------------------
15// DDIRCLut (Constructor)
16//---------------------------------
17DDIRCLut::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.65; // 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 radians
39 DIRC_SIGMA_THETAC = 0.008;
40 gPARMS->SetDefaultParameter("DIRC:SIGMA_THETAC",DIRC_SIGMA_THETAC);
41
42 // Rotate tracks angle based on bar box survey data
43 DIRC_ROTATE_TRACK = true;
44 gPARMS->SetDefaultParameter("DIRC:ROTATE_TRACK",DIRC_ROTATE_TRACK);
45
46 // PMT angle offsets for each bar from CCDB table (deprecated)
47 DIRC_THETAC_OFFSET = false;
48 gPARMS->SetDefaultParameter("DIRC:THETAC_OFFSET",DIRC_THETAC_OFFSET);
49
50 // PMT angle and time offsets to correct LUT from CCDB resource
51 DIRC_LUT_CORR = true;
52 gPARMS->SetDefaultParameter("DIRC:LUT_CORR",DIRC_LUT_CORR);
53
54 // CHROMATIC CORRECTION
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 // CHERENKOV ANGLE CUT
61 DIRC_CUT_ANGLE = 0.03;
62 gPARMS->SetDefaultParameter("DIRC:CUT_ANGLE",DIRC_CUT_ANGLE);
63
64 // set PID for different passes in debuging histograms
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); // n_quarzt = 1.47125; //(1.47125 <==> 390nm)
73 dIndex = 1.473;
74
75 if(DIRC_DEBUG_HISTS)
76 CreateDebugHistograms();
77}
78
79bool DDIRCLut::brun(JEventLoop *loop) {
80
81 dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
82 dDIRCLutReader = dapp->GetDIRCLut(loop->GetJEvent().GetRunNumber());
83
84 // get DIRC geometry
85 vector<const DDIRCGeometry*> locDIRCGeometry;
86 loop->Get(locDIRCGeometry);
87 dDIRCGeometry = locDIRCGeometry[0];
88
89 // get cherenkov angle corrections from CCDB
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 // get track rotation corrections from CCDB
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
115bool DDIRCLut::CreateDebugHistograms() {
116
117 ////////////////////////////////////
118 // LUT histograms and diagnostics //
119 ////////////////////////////////////
120
121 //dapp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
122 {
123 //TDirectory *mainDir = gDirectory;
124 //TDirectory *dircDir = gDirectory->mkdir("DIRC_debug");
125 //dircDir->cd();
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 // DeltaThetaC for each particle type (e, pi, K, p) and per pixel
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 //mainDir->cd();
161 }
162 //dapp->RootUnLock(); //REMOVE ROOT LOCK!!
163
164 return true;
165}
166
167bool 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 // get bar and track position/momentum from extrapolation
170 TVector3 momInBar = locProjMom;
171 TVector3 posInBar = locProjPos;
172
173 ////////////////////////////////////////////
174 // option to cheat and use truth position //
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 // clear object to store good photons
198 if(locDIRCMatchParams == nullptr)
199 locDIRCMatchParams = std::make_shared<DDIRCMatchParams>();
200
201 // initialize variables for LUT summary
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 // loop over DIRC hits
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 // count good photons and add hits to associated objects
224 nPhotons++;
225 locDIRCTrackMatchParams[locDIRCMatchParams].push_back(locDIRCHit);
226 }
227 }// end loop over hits
228
229 if(DIRC_DEBUG_HISTS) {
230 dapp->RootWriteLock();
231 hNph->Fill(nPhotons);
232 dapp->RootUnLock();
233 }
234
235 // skip tracks without enough photons
236 if(nPhotons<5)
237 return false;
238
239 // set DIRCMatchParameters contents
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
255vector<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 // initialize photon pairs for time and thetaC
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 // get bar number from geometry
270 int bar = dDIRCGeometry->GetBar(posInBar.Y());
271 if(bar < 0 || bar > 47) return locDIRCPhotons;
272
273 // get channel information for LUT
274 int channel = locDIRCHit->ch;
275
276 // get box from bar number (North/Upper = 0 and South/Lower = 1)
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 // use hit time to determine if reflected or not
294 double hitTime = locDIRCHit->t - locFlightTime;
295
296 // option to use truth hit time rather than smeared hit time
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 // needs to be X dependent choice for reflection cut (from CCDB?)
305 bool reflected = hitTime>35; // try only some photons as reflected for now
306
307 // get position along bar for calculated time
308 double radiatorL = dDIRCGeometry->GetBarLength(bar);
309 double barend = dDIRCGeometry->GetBarEnd(bar);
310 double lenz = fabs(barend - posInBar.X());
311
312 // get length for reflected and direct photons
313 double rlenz = 2*radiatorL - lenz; // reflected
314 double dlenz = lenz; // direct
315 if(reflected) lenz = 2*radiatorL - lenz;
Value stored to 'lenz' is never read
316
317 // check for pixel before going through loop
318 if(dDIRCLutReader->GetLutPixelAngleSize(bar, box_channel) == 0)
319 return locDIRCPhotons;
320
321 // loop over LUT table for this bar/pixel to calculate thetaC
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 // in MC we can check if the path of the LUT and measured photon are the same
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 // LUT time corrections
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 // calculate time difference
362 double locDeltaT = totalTime-hitTime;
363
364 TVector3 trackMom = momInBar;
365 if(DIRC_ROTATE_TRACK) { // rotate tracks to bar plane from survey data
366 trackMom.RotateX(rotX);
367 trackMom.RotateY(rotY);
368 trackMom.RotateZ(rotZ);
369 }
370 tangle = trackMom.Angle(dir);
371
372 if(DIRC_THETAC_OFFSET) { // ad-hoc correction per-PMT
373 tangle -= dThetaCOffset[bar][box_pmt];
374 }
375
376 // LUT angle corrections
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 // chromatic correction
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 // save hits array which pass some lose time and angle criteria
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 // reject photons that are too far out of time
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 //hDeltaThetaC[locPID]->Fill(tangle-locAngle);
417 //hDeltaThetaC_Pixel[locPID]->Fill(channel, tangle-locAngle);
418 dapp->RootUnLock();
419 }
420
421 // remove photon candidates not used in likelihood
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 // save good photons to matched list
430 isGood = true;
431
432 // count good photons
433 nPhotonsThetaC++;
434
435 // calculate average ThetaC and DeltaT
436 meanThetaC += tangle;
437 meanDeltaT += locDeltaT;
438
439
440 // calculate likelihood for each mass hypothesis
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 } // end loop over reflections
447 } // end loop over nodes
448
449 return locDIRCPhotons;
450}
451
452// overloaded function when calculating outside LUT factory
453vector<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
461double 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
468double DDIRCLut::CalcAngle(double locP, double locMass) const {
469 return acos(sqrt(locP*locP + locMass*locMass)/locP/dIndex);
470}
471
472map<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}