Bug Summary

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