Bug Summary

File:/home/sdobbs/work/clang/halld_recon/src/libraries/DIRC/DDIRCLut.cc
Warning:line 297, column 16
Value stored to 'lenz' is never read

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -main-file-name DDIRCLut.cc -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0 -D HAVE_CCDB -D HAVE_RCDB -D HAVE_EVIO -D HAVE_TMVA=1 -D RCDB_MYSQL=1 -D RCDB_SQLITE=1 -D SQLITE_USE_LEGACY_STRUCT=ON -I .Linux_CentOS7.7-x86_64-gcc4.8.5/libraries/DIRC -I libraries/DIRC -I . -I libraries -I libraries/include -I /w/halld-scifs17exp/home/sdobbs/clang/halld_recon/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I external/xstream/include -I /usr/include/tirpc -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/root/root-6.08.06/include -I /w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/ccdb/ccdb_1.06.06/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/rcdb/rcdb_0.06.00/cpp/include -I /usr/include/mysql -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlitecpp/SQLiteCpp-2.2.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlite/sqlite-3.13.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/hdds/hdds-4.9.0/Linux_CentOS7.7-x86_64-gcc4.8.5/src -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/xerces-c/xerces-c-3.1.4/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/evio/evio-4.4.6/Linux-x86_64/include -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5 -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/x86_64-redhat-linux -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/backward -internal-isystem /usr/local/include -internal-isystem /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -std=c++11 -fdeprecated-macro -fdebug-compilation-dir /home/sdobbs/work/clang/halld_recon/src -ferror-limit 19 -fgnuc-version=4.2.1 -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -o /tmp/scan-build-2021-01-21-110224-160369-1 -x c++ libraries/DIRC/DDIRCLut.cc
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.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
65bool 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
101bool 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
153bool 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
241vector<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
415vector<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
423double 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
430double DDIRCLut::CalcAngle(double locP, double locMass) const {
431 return acos(sqrt(locP*locP + locMass*locMass)/locP/dIndex);
432}
433
434map<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}