Bug Summary

File:libraries/FCAL/DFCALShower_factory.cc
Location:line 508, column 7
Description:Value stored to 'Egamma' is never read

Annotated Source Code

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