Bug Summary

File:/home/sdobbs/work/clang/halld_recon/src/libraries/FCAL/DFCALShower_factory.cc
Warning:line 500, column 2
Value stored to 'Egamma' 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 DFCALShower_factory.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/FCAL -I libraries/FCAL -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/FCAL/DFCALShower_factory.cc
1//
2// File: DFCALShower_factory.cc
3// Created: Tue May 17 11:57:50 EST 2005
4// Creator: remitche (on Linux mantrid00 2.4.20-18.8smp i686)
5
6#include <thread>
7#include <math.h>
8#include <DVector3.h>
9#include "TH2F.h"
10#include "TROOT.h"
11#include "TDirectory.h"
12using namespace std;
13
14#include "FCAL/DFCALShower_factory.h"
15#include "FCAL/DFCALGeometry.h"
16#include "FCAL/DFCALCluster.h"
17#include "FCAL/DFCALHit.h"
18#include "TRACKING/DTrackWireBased.h"
19#include <JANA/JEvent.h>
20#include <JANA/JApplication.h>
21using namespace jana;
22
23//----------------
24// Constructor
25//----------------
26DFCALShower_factory::DFCALShower_factory()
27{
28 // should we use CCDB constants?
29 LOAD_NONLIN_CCDB = 1.;
30 LOAD_TIMING_CCDB = 1.;
31 // 29/03/2020 ijaegle@jlab.org decouple non linear and timing correction
32 gPARMS->SetDefaultParameter("FCAL:LOAD_NONLIN_CCDB", LOAD_NONLIN_CCDB);
33 gPARMS->SetDefaultParameter("FCAL:LOAD_TIMING_CCDB", LOAD_TIMING_CCDB);
34
35 SHOWER_ENERGY_THRESHOLD = 50*k_MeV;
36 gPARMS->SetDefaultParameter("FCAL:SHOWER_ENERGY_THRESHOLD", SHOWER_ENERGY_THRESHOLD);
37
38 // these need to come from database to ensure accuracy
39 // remove default value which might be close to the right solution,
40 // but not quite correct -- allow command line tuning
41
42 cutoff_energy= 0;
43 linfit_slope = 0;
44 linfit_intercept = 0;
45 expfit_param1 = 0;
46 expfit_param2 = 0;
47 expfit_param3 = 0;
48
49 timeConst0 = 0;
50 timeConst1 = 0;
51 timeConst2 = 0;
52 timeConst3 = 0;
53 timeConst4 = 0;
54
55 gPARMS->SetDefaultParameter("FCAL:cutoff_energy", cutoff_energy);
56 gPARMS->SetDefaultParameter("FCAL:linfit_slope", linfit_slope);
57 gPARMS->SetDefaultParameter("FCAL:linfit_intercept", linfit_intercept);
58 gPARMS->SetDefaultParameter("FCAL:expfit_param1", expfit_param1);
59 gPARMS->SetDefaultParameter("FCAL:expfit_param2", expfit_param2);
60 gPARMS->SetDefaultParameter("FCAL:expfit_param3", expfit_param3);
61
62 gPARMS->SetDefaultParameter("FCAL:P0", timeConst0);
63 gPARMS->SetDefaultParameter("FCAL:P1", timeConst1);
64 gPARMS->SetDefaultParameter("FCAL:P2", timeConst2);
65 gPARMS->SetDefaultParameter("FCAL:P3", timeConst3);
66 gPARMS->SetDefaultParameter("FCAL:P4", timeConst4);
67
68 // Parameters to make shower-depth correction taken from Radphi,
69 // slightly modifed to match photon-polar angle
70 FCAL_RADIATION_LENGTH = 0;
71 FCAL_CRITICAL_ENERGY = 0;
72 FCAL_SHOWER_OFFSET = 0;
73
74 gPARMS->SetDefaultParameter("FCAL:FCAL_RADIATION_LENGTH", FCAL_RADIATION_LENGTH);
75 gPARMS->SetDefaultParameter("FCAL:FCAL_CRITICAL_ENERGY", FCAL_CRITICAL_ENERGY);
76 gPARMS->SetDefaultParameter("FCAL:FCAL_SHOWER_OFFSET", FCAL_SHOWER_OFFSET);
77
78 VERBOSE = 0; ///< >0 once off info ; >2 event by event ; >3 everything
79 COVARIANCEFILENAME = ""; ///< Setting the filename will take precidence over the CCDB. Files must end in ij.txt, where i and j are integers corresponding to the element of the matrix
80 gPARMS->SetDefaultParameter("DFCALShower:VERBOSE", VERBOSE, "Verbosity level for DFCALShower objects and factories");
81 gPARMS->SetDefaultParameter("DFCALShower:COVARIANCEFILENAME", COVARIANCEFILENAME, "File name for covariance files");
82
83
84 log_position_const = 4.2;
85 gPARMS->SetDefaultParameter("FCAL:log_position_const", log_position_const);
86
87
88 INSERT_RADIATION_LENGTH = 0.89;
89 INSERT_CRITICAL_ENERGY = 0.00964;
90 INSERT_SHOWER_OFFSET = 1.0;
91
92 INSERT_PAR1=1.345;
93 INSERT_PAR2=0.04;
94 INSERT_PAR3=1.16;
95 INSERT_PAR4=2.;
96 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR1",INSERT_PAR1);
97 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR2",INSERT_PAR2);
98 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR3",INSERT_PAR3);
99 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR4",INSERT_PAR4);
100
101
102}
103
104//------------------
105// brun
106//------------------
107jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
108{
109
110 map<string, double> depth_correction_params;
111 if(loop->GetCalib("FCAL/depth_correction_params", depth_correction_params)) {
112 jerr << "Problem loading FCAL/depth_correction_params from CCDB!" << endl;
113 } else {
114 FCAL_RADIATION_LENGTH = depth_correction_params["radiation_length"];
115 FCAL_CRITICAL_ENERGY = depth_correction_params["critical_energy"];
116 FCAL_SHOWER_OFFSET = depth_correction_params["shower_offset"];
117 }
118
119
120 // Get calibration constants
121 map<string, double> fcal_parms;
122 loop->GetCalib("FCAL/fcal_parms", fcal_parms);
123 if (fcal_parms.find("FCAL_C_EFFECTIVE")!=fcal_parms.end()){
124 FCAL_C_EFFECTIVE = fcal_parms["FCAL_C_EFFECTIVE"];
125 if(debug_level>0)jout<<"FCAL_C_EFFECTIVE = "<<FCAL_C_EFFECTIVE<<endl;
126 } else {
127 jerr<<"Unable to get FCAL_C_EFFECTIVE from FCAL/fcal_parms in Calib database!"<<endl;
128 }
129
130 DApplication *dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
131 const DGeometry *geom = dapp->GetDGeometry(runnumber);
132
133 if (geom) {
134 geom->GetTargetZ(m_zTarget);
135 loop->GetSingle(fcalGeom);
136 m_FCALfront=fcalGeom->fcalFrontZ();
137 m_insertFront=fcalGeom->insertFrontZ();
138 }
139 else{
140
141 cerr << "No geometry accessible." << endl;
142 return RESOURCE_UNAVAILABLE;
143 }
144 // 29/03/2020 ijaegle@jlab.org add x,y
145 jana::JCalibration *jcalib = japp->GetJCalibration(runnumber);
146 std::map<string, float> beam_spot;
147 jcalib->Get("PHOTON_BEAM/beam_spot", beam_spot);
148
149 // by default, load non-linear shower corrections from the CCDB
150 // but allow these to be overridden by command line parameters
151 energy_dependence_correction_vs_ring.clear();
152 if(LOAD_NONLIN_CCDB > 0.1) {
153 map<string, double> shower_calib_piecewise;
154 loop->GetCalib("FCAL/shower_calib_piecewise", shower_calib_piecewise);
155 cutoff_energy = shower_calib_piecewise["cutoff_energy"];
156 linfit_slope = shower_calib_piecewise["linfit_slope"];
157 linfit_intercept = shower_calib_piecewise["linfit_intercept"];
158 expfit_param1 = shower_calib_piecewise["expfit_param1"];
159 expfit_param2 = shower_calib_piecewise["expfit_param2"];
160 expfit_param3 = shower_calib_piecewise["expfit_param3"];
161 m_beamSpotX = 0;
162 m_beamSpotY = 0;
163
164 if(debug_level>0) {
165 jout << "cutoff_energy = " << cutoff_energy << endl;
166 jout << "linfit_slope = " << linfit_slope << endl;
167 jout << "linfit_intercept = " << linfit_intercept << endl;
168 jout << "expfit_param1 = " << expfit_param1 << endl;
169 jout << "expfit_param2 = " << expfit_param2<< endl;
170 jout << "expfit_param3 = " << expfit_param3 << endl;
171 }
172 loop->GetCalib("FCAL/energy_dependence_correction_vs_ring", energy_dependence_correction_vs_ring);
173 if (energy_dependence_correction_vs_ring.size() > 0 && energy_dependence_correction_vs_ring[0][0] != 0) {
174 m_beamSpotX = beam_spot.at("x");
175 m_beamSpotY = beam_spot.at("y");
176 if (debug_level > 0) {
177 TString str_coef[] = {"A", "B", "C", "D", "E", "F"};
178 for (int i = 0; i < 24; i ++) {
179 //for (int j = 0; j < 6; j ++) {
180 for (int j = 0; j < 3; j ++) {
181 jout << "Ring # " << i << Form(" %s", str_coef[j].Data()) << energy_dependence_correction_vs_ring[i][j];
182 }
183 jout << endl;
184 }
185 }
186 }
187 }
188
189 if (LOAD_TIMING_CCDB > 0.1) {
190 // Get timing correction polynomial, J. Mirabelli 10/31/17
191 map<string,double> timing_correction;
192 loop->GetCalib("FCAL/shower_timing_correction", timing_correction);
193 timeConst0 = timing_correction["P0"];
194 timeConst1 = timing_correction["P1"];
195 timeConst2 = timing_correction["P2"];
196 timeConst3 = timing_correction["P3"];
197 timeConst4 = timing_correction["P4"];
198
199 if(debug_level>0) {
200 jout << "timeConst0 = " << timeConst0 << endl;
201 jout << "timeConst1 = " << timeConst1 << endl;
202 jout << "timeConst2 = " << timeConst2 << endl;
203 jout << "timeConst3 = " << timeConst3 << endl;
204 jout << "timeConst4 = " << timeConst4 << endl;
205 }
206 }
207
208 jerror_t result = LoadCovarianceLookupTables(eventLoop);
209 if (result!=NOERROR) return result;
210
211 INSERT_C_EFFECTIVE=FCAL_C_EFFECTIVE;
212
213 return NOERROR;
214}
215
216
217jerror_t DFCALShower_factory::erun(void) {
218 // delete lookup tables to prevent memory leak
219 for (int i=0; i<5; i++) {
220 for (int j=0; j<=i; j++) {
221 delete CovarianceLookupTable[i][j];
222 CovarianceLookupTable[i][j] = nullptr;
223 }
224 }
225 return NOERROR;
226}
227
228
229//------------------
230// evnt
231//------------------
232jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber)
233{
234 vector<const DFCALCluster*> fcalClusters;
235 eventLoop->Get(fcalClusters);
236 if(fcalClusters.size()<1)return NOERROR;
237
238 // Use the center of the target as an approximation for the vertex position
239 // 29/03/2020 ijaegle@jlab.org add beam center in x,y
240 DVector3 vertex(m_beamSpotX, m_beamSpotY, m_zTarget);
241
242 vector< const DTrackWireBased* > allWBTracks;
243 eventLoop->Get( allWBTracks );
244 vector< const DTrackWireBased* > wbTracks = filterWireBasedTracks( allWBTracks );
245
246 // Loop over list of DFCALCluster objects and calculate the "Non-linear" corrected
247 // energy and position for each. We'll use a logarithmic energy-weighting to
248 // find the final position and error.
249 for( vector< const DFCALCluster* >::const_iterator clItr = fcalClusters.begin();
250 clItr != fcalClusters.end(); ++clItr ){
251 const DFCALCluster* cluster=*clItr;
252
253 // energy weighted time provides better resolution:
254 double cTime = cluster->getTimeEWeight();
255
256 double zback=m_FCALfront + fcalGeom->blockLength();
257 double c_effective=FCAL_C_EFFECTIVE;
258
259 int channel = cluster->getChannelEmax();
260 DVector2 pos=fcalGeom->positionOnFace(channel);
261 // Check if the cluster is in the insert
262 bool in_insert=fcalGeom->inInsert(channel);
263 if (in_insert){
264 zback=m_insertFront + fcalGeom->insertBlockLength();
265 c_effective=INSERT_C_EFFECTIVE;
266 in_insert=true;
267 }
268
269 // Get corrected energy, position, and errZ
270 double Ecorrected;
271 DVector3 pos_corrected;
272 double errZ;
273 double radius = pos.Mod();
274 int ring_nb = (int) (radius / (5 * k_cm));
275 GetCorrectedEnergyAndPosition( cluster, ring_nb , Ecorrected, pos_corrected, errZ, &vertex,in_insert);
276
277 DVector3 pos_log;
278 GetLogWeightedPosition( cluster, pos_log, Ecorrected, &vertex );
279
280 if (Ecorrected>0.){
281 //up to this point, all times have been times at which light reaches
282 //the back of the detector. Here we correct for the time that it
283 //takes the Cherenkov light to reach the back of the detector
284 //so that the t reported is roughly the time of the shower at the
285 //position pos_corrected
286 cTime -= ( zback - pos_corrected.Z() )/c_effective;
287
288 //Apply time-walk correction/global timing offset
289 cTime += ( timeConst0 + timeConst1 * Ecorrected + timeConst2 * TMath::Power( Ecorrected, 2 ) +
290 timeConst3 * TMath::Power( Ecorrected, 3 ) + timeConst4 * TMath::Power( Ecorrected, 4 ) );
291
292 // Make the DFCALShower object
293 DFCALShower* shower = new DFCALShower;
294
295 shower->setEnergy( Ecorrected );
296 shower->setPosition( pos_corrected );
297 shower->setPosition_log( pos_log );
298 shower->setTime ( cTime );
299
300 if (in_insert==false){
301 FillCovarianceMatrix( shower );
302 }
303 else{
304 // Some guesses for insert resolution, currently hard-coded...
305 double sigx=0.1016/sqrt(Ecorrected)+0.2219;
306 shower->ExyztCovariance(1,1)=sigx*sigx;
307 shower->ExyztCovariance(2,2)=sigx*sigx;
308 shower->ExyztCovariance(0,0)=Ecorrected*Ecorrected*(0.01586/Ecorrected
309 +0.0002342/(Ecorrected*Ecorrected)
310 +1.695e-6);
311 for (unsigned int i=0;i<5;i++){
312 for(unsigned int j=0;j<5;j++){
313 if (i!=j) shower->ExyztCovariance(i,j)=0.;
314 }
315
316 }
317 }
318
319 if( VERBOSE > 2 ){
320 printf("FCAL shower: } E=%f x=%f y=%f z=%f t=%f\n",
321 shower->getEnergy(),shower->getPosition().X(),shower->getPosition().Y(),shower->getPosition().Z(),shower->getTime());
322 printf("FCAL shower: dE=%f dx=%f dy=%f dz=%f dt=%f\n",
323 shower->EErr(),shower->xErr(),shower->yErr(),shower->zErr(),shower->tErr());
324 printf("FCAL shower: Ex=%f Ey=%f Ez=%f Et=%f xy=%f\n",
325 shower->EXcorr(),shower->EYcorr(),shower->EZcorr(),shower->ETcorr(),shower->XYcorr());
326 printf("FCAL shower: xz=%f xt=%f yz=%f yt=%f zt=%f\n",
327 shower->XZcorr(),shower->XTcorr(),shower->YZcorr(),shower->YTcorr(),shower->ZTcorr());
328 }
329
330 // now fill information related to shower shape and nearby
331 // tracks -- useful for splitoff rejection later
332
333 double docaTr = 1E6;
334 double timeTr = 1E6;
335 double xTr = 0;
336 double yTr = 0;
337
338 double flightTime;
339 DVector3 projPos, projMom;
340
341 // find the closest track to the shower -- here we loop over the best FOM
342 // wire-based track for every track candidate not just the ones associated
343 // with the topology
344 for( size_t iTrk = 0; iTrk < wbTracks.size(); ++iTrk ){
345
346 if( !wbTracks[iTrk]->GetProjection( SYS_FCAL, projPos, &projMom, &flightTime ) ) continue;
347
348 // need to swim fcalPos to common z for DOCA calculation -- this really
349 // shouldn't be in the loop if the z-value of projPos doesn't change
350 // with each track
351
352 DVector3 fcalFacePos = ( shower->getPosition() - vertex );
353 fcalFacePos.SetMag( fcalFacePos.Mag() * projPos.Z() / fcalFacePos.Z() );
354
355 double distance = ( fcalFacePos - projPos ).Mag();
356
357 if( distance < docaTr ){
358
359 docaTr = distance;
360 // this is the time from the center of the target to the detector -- to compare with
361 // the FCAL time, one needs to have the t0RF at the center of the target. That
362 // comparison happens at a later stage in the analysis.
363 timeTr = ( wbTracks[iTrk]->position().Z() - vertex.Z() ) / SPEED_OF_LIGHT29.9792458 + flightTime;
364 xTr = projPos.X();
365 yTr = projPos.Y();
366 }
367 }
368
369 shower->setDocaTrack( docaTr );
370 shower->setTimeTrack( timeTr );
371
372 // now compute some variables at the hit level
373
374 vector< const DFCALHit* > fcalHits;
375 cluster->Get( fcalHits );
376 shower->setNumBlocks( fcalHits.size() );
377
378 double e9e25, e1e9;
379 getE1925FromHits( e1e9, e9e25, fcalHits, getMaxHit( fcalHits ) );
380 shower->setE1E9( e1e9 );
381 shower->setE9E25( e9e25 );
382
383 double sumU = 0;
384 double sumV = 0;
385 // if there is no nearest track, the defaults for xTr and yTr will result
386 // in using the beam axis as the directional axis
387 getUVFromHits( sumU, sumV, fcalHits,
388 DVector3( shower->getPosition().X(), shower->getPosition().Y(), 0 ),
389 DVector3( xTr, yTr, 0 ) );
390
391 shower->setSumU( sumU );
392 shower->setSumV( sumV );
393
394 shower->AddAssociatedObject( cluster );
395
396 _data.push_back(shower);
397 }
398 }
399
400 return NOERROR;
401}
402
403//--------------------------------
404// GetCorrectedEnergyAndPosition
405//
406// Non-linear and depth corrections should be fixed within DFCALShower member functions
407//--------------------------------
408 void DFCALShower_factory::GetCorrectedEnergyAndPosition(const DFCALCluster* cluster, int ring_nb, double &Ecorrected, DVector3 &pos_corrected, double &errZ, const DVector3 *vertex,bool in_insert)
409{
410 // Non-linear energy correction are done here
411 //int MAXITER = 1000;
412
413 DVector3 posInCal = cluster->getCentroid();
414 float x0 = posInCal.Px();
415 float y0 = posInCal.Py();
416 double Eclust = cluster->getEnergy();
417
418 double Ecutoff = 0;
419 double A = 0;
420 double B = 0;
421 double C = 0;
422 double D = 0;
423 double E = 0;
424 double Egamma = Eclust;
425 Ecorrected = 0;
426
427 // block properties
428 double radiation_length=FCAL_RADIATION_LENGTH;
429 double shower_offset=FCAL_SHOWER_OFFSET;
430 double critical_energy=FCAL_CRITICAL_ENERGY;
431 double zfront=m_FCALfront;
432
433 // Check for presence of insert
434 if (in_insert){
435 radiation_length=INSERT_RADIATION_LENGTH;
436 shower_offset=INSERT_SHOWER_OFFSET;
437 critical_energy=INSERT_CRITICAL_ENERGY;
438 zfront=m_insertFront;
439
440 A=INSERT_PAR1;
441 B=INSERT_PAR2;
442 C=INSERT_PAR3;
443 D=INSERT_PAR4;
444 if (Eclust<D){
445 Egamma=A*Eclust/(1.+B*Eclust);
446 }
447 else Egamma=A*D/(1.+D*B)+C*(Eclust-D);
448 }
449 else{
450 int ring_region = -1;
451 if (0 <= ring_nb && ring_nb <= 2)
452 ring_region = 0;
453 else if (3 <= ring_nb && ring_nb <= 4)
454 ring_region = 1;
455 else if (ring_nb == 5)
456 ring_region = 2;
457 else if (6 <= ring_nb && ring_nb <= 7)
458 ring_region = 3;
459 else if (8 <= ring_nb && ring_nb <= 9)
460 ring_region = 4;
461 else if (10 <= ring_nb && ring_nb <= 11)
462 ring_region = 5;
463 else if (12 <= ring_nb && ring_nb <= 17)
464 ring_region = 6;
465 else if (18 <= ring_nb && ring_nb <= 20)
466 ring_region = 7;
467 else if (21 <= ring_nb && ring_nb <= 23)
468 ring_region = 8;
469
470 // 06/04/2020 ijaegle@jlab.org allows two different energy dependence correction
471 if (LOAD_NONLIN_CCDB > 0.1) {
472
473 // Method I: IU way, one overall correction
474 Egamma = 0;
475 Ecutoff = cutoff_energy;
476 A = linfit_slope;
477 B = linfit_intercept;
478 C = expfit_param1;
479 D = expfit_param2;
480 E = expfit_param3;
481 // 06/02/2016 Shower Non-linearity Correction by Adesh.
482 // 29/03/2020 ijaegle@jlab.org the linear part correction is applied in (some) data/sim. backward comptability?
483 if ( Eclust <= Ecutoff ) {
484
485 Egamma = Eclust / (A * Eclust + B); // Linear part
486
487 } else
488 // 29/03/2020 ijaegle@jlab.org this correction is always applied if all C=2 & D=E=0 then Egamma = Eclust
489 if ( Eclust > Ecutoff ) {
490
491 Egamma = Eclust / (C - exp(-D * Eclust + E)); // Non-linear part
492
493 }
494 // 29/03/2020 ijaegle@jlab.org if all C=D=E=0 by mistake then Egamma = - Eclust
495 // End Correction method I
496
497 // Method II: PRIMEXD way, correction per ring
498 if (C == 2 && D == 0 && E == 0 && energy_dependence_correction_vs_ring.size() > 0 && ring_region != -1) {
499
500 Egamma = 0;
Value stored to 'Egamma' is never read
501 A = energy_dependence_correction_vs_ring[ring_region][0];
502 B = energy_dependence_correction_vs_ring[ring_region][1];
503 C = energy_dependence_correction_vs_ring[ring_region][2];
504 //D = energy_dependence_correction_vs_ring[ring_nb][3];
505 //E = energy_dependence_correction_vs_ring[ring_nb][4];
506 //F = energy_dependence_correction_vs_ring[ring_nb][5];
507 //Egamma = Eclust / (A + B * Eclust + C * pow(Eclust, 2) + D * pow(Eclust, 3) + E * pow(Eclust, 4) + F * pow(Eclust, 5));
508 //Egamma = Eclust / (A + B * Eclust + C * pow(Eclust, 2));
509 Egamma = Eclust / (A - exp(-B * Eclust + C));
510 }
511 // End Correction method II
512 }
513 //End energy dependence correction
514 }
515 if (Egamma <= 0 && Eclust > 0) Egamma = Eclust;
516
517 // then depth corrections
518 if ( Egamma > 0 ) {
519 float dxV = x0-vertex->X();
520 float dyV = y0-vertex->Y();
521 float zV = vertex->Z();
522
523 double z0 = zfront - zV;
524 double zMax = radiation_length*(shower_offset+log(Egamma/critical_energy));
525
526 double zed = z0;
527 double zed1 = z0 + zMax;
528
529 double r0 = sqrt(dxV*dxV + dyV*dyV );
530
531 int niter;
532 for ( niter=0; niter<100; niter++) {
533 double tt = r0/zed1;
534 zed = z0 + zMax/sqrt( 1 + tt*tt );
535 if ( fabs( (zed-zed1) ) < 0.001) {
536 break;
537 }
538 zed1 = zed;
539 }
540
541 posInCal.SetZ( zed + zV );
542 errZ = zed - zed1;
543 }
544
545 Ecorrected = Egamma;
546 pos_corrected = posInCal;
547
548}
549
550
551
552jerror_t
553DFCALShower_factory::FillCovarianceMatrix(DFCALShower *shower){
554 /// This function takes a FCALShower object and using the internal variables
555 /// overwrites any existing covaraince matrix using lookup tables.
556
557 // Get edges of lookup table histograms (assume that all histograms have the same limits.)
558 TAxis *xaxis = CovarianceLookupTable[0][0]->GetXaxis();
559 TAxis *yaxis = CovarianceLookupTable[0][0]->GetYaxis();
560 float minElookup = xaxis->GetBinLowEdge(1);
561 float maxElookup = xaxis->GetBinUpEdge(xaxis->GetNbins());
562 float minthlookup = yaxis->GetBinLowEdge(1);
563 float maxthlookup = yaxis->GetBinUpEdge(yaxis->GetNbins());
564
565 float shower_E = shower->getEnergy();
566 float shower_x = shower->getPosition().X();
567 float shower_y = shower->getPosition().Y();
568 float shower_z = shower->getPosition().Z();
569 float shower_r = sqrt(shower_x*shower_x + shower_y*shower_y);
570 float shower_theta = atan2(shower_r,shower_z);
571 float thlookup = shower_theta/3.14159265*180;
572 float Elookup = shower_E;
573
574 // Adjust values: in order to use Interpolate() must be within histogram range
575 if (Elookup<minElookup) Elookup=minElookup;
576 if (Elookup>maxElookup) Elookup=maxElookup-0.0001; // move below edge, on edge doesn't work.
577 if (thlookup<minthlookup) thlookup=minthlookup;
578 if (thlookup>maxthlookup) thlookup=maxthlookup-0.0001;
579 if (VERBOSE>3) printf("(%f,%F) limits (%f,%f) (%f,%f)\n",Elookup,thlookup,minElookup,maxElookup,minthlookup,maxthlookup);
580
581 DMatrixDSym ErphiztCovariance(5);
582 for (int i=0; i<5; i++) {
583 for (int j=0; j<=i; j++) {
584 float val = CovarianceLookupTable[i][j]->Interpolate(Elookup, thlookup);
585 if (i==0 && j==0) val *= shower_E; // E variance is divided by energy in CCDB
586 ErphiztCovariance(i,j) = ErphiztCovariance(j,i) = val;
587 }
588 }
589
590 float shower_phi = atan2(shower_y,shower_x);
591 float cosPhi = cos(shower_phi);
592 float sinPhi = sin(shower_phi);
593 DMatrix rotationmatrix(5,5);
594 rotationmatrix(0,0) = 1;
595 rotationmatrix(3,3) = 1;
596 rotationmatrix(4,4) = 1;
597 rotationmatrix(1,1) = cosPhi;
598 rotationmatrix(1,2) = -sinPhi;
599 rotationmatrix(2,1) = sinPhi;
600 rotationmatrix(2,2) = cosPhi;
601
602 if (VERBOSE>3) {printf("(E,r,phi,z,t) "); ErphiztCovariance.Print(); }
603 DMatrixDSym &D = ErphiztCovariance.Similarity(rotationmatrix);
604 for (int i=0; i<5; i++) {
605 for (int j=0; j<5; j++)
606 shower->ExyztCovariance(i, j) = D(i, j);
607 }
608 if (VERBOSE>2) {printf("(E,x,y,z,t) "); shower->ExyztCovariance.Print(); }
609
610 return NOERROR;
611}
612
613
614jerror_t
615DFCALShower_factory::LoadCovarianceLookupTables(JEventLoop *eventLoop){
616 std::thread::id this_id = std::this_thread::get_id();
617 stringstream idstring;
618 idstring << this_id;
619 if (VERBOSE>0) printf("DFCALShower_factory::LoadCovarianceLookupTables(): Thread %s\n",idstring.str().c_str());
620
621 bool USECCDB=0;
622 bool DUMMYTABLES=0;
623 // if filename specified try to use filename else get info from CCDB
624 if (COVARIANCEFILENAME == "") USECCDB=1;
625
626 map<string,string> covariance_data;
627 if (USECCDB) {
628 // load information for covariance matrix
629 if (eventLoop->GetJCalibration()->GetCalib("/FCAL/shower_covariance", covariance_data)) {
630 jerr << "Error loading /FCAL/shower_covariance !" << endl;
631 DUMMYTABLES=1;
632 }
633 if (covariance_data.size() == 15) { // there are 15 elements in the covariance matrix
634 // for example, print it all out
635 if (VERBOSE>0) {
636 for(auto element : covariance_data) {
637 cout << "\nTEST: " << element.first << " = " << element.second << endl;
638 }
639 }
640 } else {
641 jerr << "Wrong number of elements /FCAL/shower_covariance !" << endl;
642 DUMMYTABLES=1;
643 }
644 }
645
646 for (int i=0; i<5; i++) {
647 for (int j=0; j<=i; j++) {
648
649 japp->RootWriteLock();
650 // change directory to memory so that histograms are not saved to file
651 TDirectory *savedir = gDirectory(TDirectory::CurrentDirectory());
652
653 char histname[255];
654 sprintf(histname,"covariance_%i%i_thread%s",i,j,idstring.str().c_str());
655 // Read in string
656 ifstream ifs;
657 string line;
658 stringstream ss;
659 if (USECCDB) {
660 stringstream matrixname;
661 matrixname << "covmatrix_" << i << j;
662 if (VERBOSE>1) cout << "Using CCDB \"" << matrixname.str() << "\" " << covariance_data[matrixname.str()] << endl;
663 ss.str(covariance_data[matrixname.str()]);
664 } else {
665 char filename[255];
666 sprintf(filename,"%s%i%i.txt",COVARIANCEFILENAME.c_str(),i,j);
667 if (VERBOSE>0) cout << filename << std::endl;
668 ifs.open(filename);
669 if (! ifs.is_open()) {
670 jerr << " Error: Cannot open file! " << filename << std::endl;
671 DUMMYTABLES=1;
672 } else {
673 getline(ifs, line, '\n');
674 ss.str(line);
675 if (VERBOSE>1) cout << filename << " dump: " <<line<<endl;
676 }
677 }
678 if (DUMMYTABLES) {
679 // create dummy histogram since something went wrong
680 CovarianceLookupTable[i][j] = new TH2F(histname,"Covariance histogram",10,0,12,10,0,12);
681 CovarianceLookupTable[i][j]->SetDirectory(nullptr);
682 } else {
683 // Parse string
684 int nxbins, nybins;
685 ss>>nxbins;
686 ss>>nybins;
687 if (VERBOSE>1) printf("parsed dump: bins (%i,%i)\n",nxbins,nybins);
688 Float_t xbins[nxbins+1];
689 Float_t ybins[nybins+1];
690 for (int count=0; count<=nxbins; count++) {
691 ss>>xbins[count];
692 if (VERBOSE>1) printf("(%i,%f) ",count,xbins[count]);
693 }
694 if (VERBOSE>1) printf("\n");
695 for (int count=0; count<=nybins; count++) {
696 ss>>ybins[count];
697 if (VERBOSE>1) printf("(%i,%f) ",count,ybins[count]);
698 }
699 if (VERBOSE>1) printf("\n");
700 int xbin=1;
701 double cont;
702 int ybin=1;
703 // create histogram
704 CovarianceLookupTable[i][j] = new TH2F(histname,"Covariance histogram",nxbins,xbins,nybins,ybins);
705 CovarianceLookupTable[i][j]->SetDirectory(nullptr);
706 // fill histogram
707 while(ss>>cont){
708 if (VERBOSE>1) printf("(%i,%i) (%i,%i) %e ",i,j,xbin,ybin,cont);
709 CovarianceLookupTable[i][j]->SetBinContent(xbin,ybin,cont);
710 ybin++;
711 if (ybin>nybins) { xbin++; ybin=1; }
712 }
713 if (VERBOSE>1) printf("\n");
714 // Close file
715 ifs.close();
716 }
717 savedir->cd();
718 japp->RootUnLock();
719 }
720 }
721 return NOERROR;
722}
723
724unsigned int
725DFCALShower_factory::getMaxHit( const vector< const DFCALHit* >& hitVec ) const {
726
727 unsigned int maxIndex = 0;
728
729 double eMaxSh = 0;
730
731 for( vector< const DFCALHit* >::const_iterator hit = hitVec.begin();
732 hit != hitVec.end(); ++hit ){
733
734 if( (**hit).E > eMaxSh ){
735
736 eMaxSh = (**hit).E;
737 maxIndex = hit - hitVec.begin();
738 }
739 }
740
741 return maxIndex;
742}
743
744void
745DFCALShower_factory::getUVFromHits( double& sumUSh, double& sumVSh,
746 const vector< const DFCALHit* >& hits,
747 const DVector3& showerVec,
748 const DVector3& trackVec ) const {
749
750 // This method forms an axis pointing from the shower to nearest track
751 // and computes the energy-weighted second moment of the shower along
752 // and perpendicular to this axis. True photons are fairly symmetric
753 // and have similar values of sumU and sumV whereas splitoffs tend
754 // to be asymmetric in these variables.
755
756 DVector3 u = ( showerVec - trackVec ).Unit();
757 DVector3 z( 0, 0, 1 );
758 DVector3 v = u.Cross( z );
759
760 DVector3 hitLoc( 0, 0, 0 );
761
762 sumUSh = 0;
763 sumVSh = 0;
764
765 double sumE = 0;
766
767 for( vector< const DFCALHit* >::const_iterator hit = hits.begin();
768 hit != hits.end(); ++hit ){
769
770 hitLoc.SetX( (**hit).x - showerVec.X() );
771 hitLoc.SetY( (**hit).y - showerVec.Y() );
772
773 sumUSh += (**hit).E * pow( u.Dot( hitLoc ), 2 );
774 sumVSh += (**hit).E * pow( v.Dot( hitLoc ), 2 );
775
776 sumE += (**hit).E;
777 }
778
779 sumUSh /= sumE;
780 sumVSh /= sumE;
781}
782
783void
784DFCALShower_factory::getE1925FromHits( double& e1e9Sh, double& e9e25Sh,
785 const vector< const DFCALHit* >& hits,
786 unsigned int maxIndex ) const {
787
788 double E9 = 0;
789 double E25 = 0;
790
791 const DFCALHit* maxHit = hits[maxIndex];
792
793 for( vector< const DFCALHit* >::const_iterator hit = hits.begin();
794 hit != hits.end(); ++hit ){
795
796 if( fabs( (**hit).x - maxHit->x ) < 4.5 && fabs( (**hit).y - maxHit->y ) < 4.5 )
797 E9 += (**hit).E;
798
799 if( fabs( (**hit).x - maxHit->x ) < 8.5 && fabs( (**hit).y - maxHit->y ) < 8.5 )
800 E25 += (**hit).E;
801 }
802
803 e1e9Sh = maxHit->E/E9;
804 e9e25Sh = E9/E25;
805}
806
807
808vector< const DTrackWireBased* >
809DFCALShower_factory::filterWireBasedTracks( vector< const DTrackWireBased* >& wbTracks ) const {
810
811 vector< const DTrackWireBased* > finalTracks;
812 map< unsigned int, vector< const DTrackWireBased* > > sortedTracks;
813
814 // first sort the wire based tracks into lists with a common candidate id
815 // this means that they all come from the same track in the detector
816
817 for( unsigned int i = 0; i < wbTracks.size(); ++i ){
818
819 unsigned int id = wbTracks[i]->candidateid;
820
821 if( sortedTracks.find( id ) == sortedTracks.end() ){
822
823 sortedTracks[id] = vector< const DTrackWireBased* >();
824 }
825
826 sortedTracks[id].push_back( wbTracks[i] );
827 }
828
829 // now loop through that list of unique tracks and for each set
830 // of wire based tracks, choose the one with the highest FOM
831 // (this is choosing among different particle hypotheses)
832
833 for( map< unsigned int, vector< const DTrackWireBased* > >::const_iterator
834 anId = sortedTracks.begin();
835 anId != sortedTracks.end(); ++anId ){
836
837 double maxFOM = 0;
838 unsigned int bestIndex = 0;
839
840 for( unsigned int i = 0; i < anId->second.size(); ++i ){
841
842 if( anId->second[i]->Ndof < 15 ) continue;
843
844 if( anId->second[i]->FOM > maxFOM ){
845
846 maxFOM = anId->second[i]->FOM;
847 bestIndex = i;
848 }
849 }
850
851 finalTracks.push_back( anId->second[bestIndex] );
852 }
853
854 return finalTracks;
855}
856
857
858void DFCALShower_factory::GetLogWeightedPosition( const DFCALCluster* cluster, DVector3 &pos_log, double Egamma, const DVector3 *vertex )
859{
860
861 DVector3 posInCal = cluster->getCentroid();
862
863 const vector< DFCALCluster::DFCALClusterHit_t > locHitVector = cluster->GetHits();
864
865 int loc_nhits = (int)locHitVector.size();
866 if( loc_nhits < 1 ) {
867 pos_log = posInCal;
868 return;
869 }
870
871 //------ Loop over hits ------//
872
873 double sW = 0.0;
874 double xpos = 0.0;
875 double ypos = 0.0;
876 double W;
877
878 double ecluster = cluster->getEnergy();
879
880 for( int ih = 0; ih < loc_nhits; ih++ ) {
881
882 DFCALCluster::DFCALClusterHit_t locHit = locHitVector[ih];
883
884 double xcell = locHit.x;
885 double ycell = locHit.y;
886 double ecell = locHit.E;
887
888 W = log_position_const + log( ecell / ecluster );
889 if( W > 0. ) {
890 sW += W;
891 xpos += xcell * W;
892 ypos += ycell * W;
893 }
894
895 }
896
897 double x1, y1;
898 if( sW ) {
899 x1 = xpos / sW;
900 y1 = ypos / sW;
901 } else {
902 cout << "\nBad Cluster Logged in DFCALShower_factory::GetLogWeightedPosition" << endl;
903 x1 = 0.;
904 y1 = 0.;
905 }
906
907
908 // Shower Depth Corrections (copied from GetCorrectedEnergyAndPosition function)
909
910 if ( Egamma > 0 ) {
911 float dxV = x1 - vertex->X();
912 float dyV = y1 - vertex->Y();
913 float zV = vertex->Z();
914
915 double z0 = m_FCALfront - zV;
916 double zMax = FCAL_RADIATION_LENGTH*(FCAL_SHOWER_OFFSET
917 + log(Egamma/FCAL_CRITICAL_ENERGY));
918
919 double zed = z0;
920 double zed1 = z0 + zMax;
921
922 double r0 = sqrt( dxV*dxV + dyV*dyV );
923
924 int niter;
925 for ( niter=0; niter<100; niter++) {
926 double tt = r0/zed1;
927 zed = z0 + zMax/sqrt( 1 + tt*tt );
928 if ( fabs( (zed-zed1) ) < 0.001) {
929 break;
930 }
931 zed1 = zed;
932 }
933
934 posInCal.SetZ( zed + zV );
935 }
936
937 posInCal.SetX( x1 );
938 posInCal.SetY( y1 );
939
940
941 pos_log = posInCal;
942
943
944 return;
945
946}
947
948