Bug Summary

File:libraries/FCAL/DFCALShower_factory.cc
Location:line 429, column 10
Description:Value stored to 'Egamma' during its initialization 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=1.365;
97 INSERT_PAR2=0.04;
98 INSERT_PAR3=1.185;
99 INSERT_PAR4=2.;
100 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR1",INSERT_PAR1);
101 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR2",INSERT_PAR2);
102 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR3",INSERT_PAR3);
103 gPARMS->SetDefaultParameter("FCAL:INSERT_PAR4",INSERT_PAR4);
104
105}
106
107//------------------
108// brun
109//------------------
110jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
111{
112
113 map<string, double> depth_correction_params;
114 if(loop->GetCalib("FCAL/depth_correction_params", depth_correction_params)) {
115 jerr << "Problem loading FCAL/depth_correction_params from CCDB!" << endl;
116 } else {
117 FCAL_RADIATION_LENGTH = depth_correction_params["radiation_length"];
118 FCAL_CRITICAL_ENERGY = depth_correction_params["critical_energy"];
119 FCAL_SHOWER_OFFSET = depth_correction_params["shower_offset"];
120 }
121
122
123 // Get calibration constants
124 map<string, double> fcal_parms;
125 loop->GetCalib("FCAL/fcal_parms", fcal_parms);
126 if (fcal_parms.find("FCAL_C_EFFECTIVE")!=fcal_parms.end()){
127 FCAL_C_EFFECTIVE = fcal_parms["FCAL_C_EFFECTIVE"];
128 if(debug_level>0)jout<<"FCAL_C_EFFECTIVE = "<<FCAL_C_EFFECTIVE<<endl;
129 } else {
130 jerr<<"Unable to get FCAL_C_EFFECTIVE from FCAL/fcal_parms in Calib database!"<<endl;
131 }
132
133 DApplication *dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
134 const DGeometry *geom = dapp->GetDGeometry(runnumber);
135
136 if (geom) {
137 geom->GetTargetZ(m_zTarget);
138 loop->GetSingle(fcalGeom);
139 m_FCALfront=fcalGeom->fcalFrontZ();
140 m_insertFront=fcalGeom->insertFrontZ();
141 }
142 else{
143
144 cerr << "No geometry accessible." << endl;
145 return RESOURCE_UNAVAILABLE;
146 }
147 // 29/03/2020 ijaegle@jlab.org add x,y
148 jana::JCalibration *jcalib = japp->GetJCalibration(runnumber);
149 std::map<string, float> beam_spot;
150 jcalib->Get("PHOTON_BEAM/beam_spot", beam_spot);
151
152 // by default, load non-linear shower corrections from the CCDB
153 // but allow these to be overridden by command line parameters
154 energy_dependence_correction_vs_ring.clear();
155 if(LOAD_NONLIN_CCDB) {
156 map<string, double> shower_calib_piecewise;
157 loop->GetCalib("FCAL/shower_calib_piecewise", shower_calib_piecewise);
158 cutoff_energy = shower_calib_piecewise["cutoff_energy"];
159 linfit_slope = shower_calib_piecewise["linfit_slope"];
160 linfit_intercept = shower_calib_piecewise["linfit_intercept"];
161 expfit_param1 = shower_calib_piecewise["expfit_param1"];
162 expfit_param2 = shower_calib_piecewise["expfit_param2"];
163 expfit_param3 = shower_calib_piecewise["expfit_param3"];
164 m_beamSpotX = 0;
165 m_beamSpotY = 0;
166
167 if(debug_level>0) {
168 jout << "cutoff_energy = " << cutoff_energy << endl;
169 jout << "linfit_slope = " << linfit_slope << endl;
170 jout << "linfit_intercept = " << linfit_intercept << endl;
171 jout << "expfit_param1 = " << expfit_param1 << endl;
172 jout << "expfit_param2 = " << expfit_param2<< endl;
173 jout << "expfit_param3 = " << expfit_param3 << endl;
174 }
175 loop->GetCalib("FCAL/energy_dependence_correction_vs_ring", energy_dependence_correction_vs_ring);
176 if (energy_dependence_correction_vs_ring.size() > 0 && energy_dependence_correction_vs_ring[0][0] != 0) {
177 m_beamSpotX = beam_spot.at("x");
178 m_beamSpotY = beam_spot.at("y");
179 if (debug_level > 0) {
180 TString str_coef[] = {"A", "B", "C", "D", "E", "F"};
181 for (int i = 0; i < 24; i ++) {
182 //for (int j = 0; j < 6; j ++) {
183 for (int j = 0; j < 3; j ++) {
184 jout << "Ring # " << i << Form(" %s", str_coef[j].Data()) << energy_dependence_correction_vs_ring[i][j];
185 }
186 jout << endl;
187 }
188 }
189 }
190 }
191
192 if (LOAD_TIMING_CCDB) {
193 // Get timing correction polynomial, J. Mirabelli 10/31/17
194 map<string,double> timing_correction;
195 loop->GetCalib("FCAL/shower_timing_correction", timing_correction);
196 timeConst0 = timing_correction["P0"];
197 timeConst1 = timing_correction["P1"];
198 timeConst2 = timing_correction["P2"];
199 timeConst3 = timing_correction["P3"];
200 timeConst4 = timing_correction["P4"];
201
202 if(debug_level>0) {
203 jout << "timeConst0 = " << timeConst0 << endl;
204 jout << "timeConst1 = " << timeConst1 << endl;
205 jout << "timeConst2 = " << timeConst2 << endl;
206 jout << "timeConst3 = " << timeConst3 << endl;
207 jout << "timeConst4 = " << timeConst4 << endl;
208 }
209 }
210
211 jerror_t result = LoadCovarianceLookupTables(eventLoop);
212 if (result!=NOERROR) return result;
213
214 INSERT_C_EFFECTIVE=FCAL_C_EFFECTIVE;
215
216 return NOERROR;
217}
218
219
220jerror_t DFCALShower_factory::erun(void) {
221 // delete lookup tables to prevent memory leak
222 for (int i=0; i<5; i++) {
223 for (int j=0; j<=i; j++) {
224 delete CovarianceLookupTable[i][j];
225 CovarianceLookupTable[i][j] = nullptr;
226 }
227 }
228 return NOERROR;
229}
230
231
232//------------------
233// evnt
234//------------------
235jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber)
236{
237 vector<const DFCALCluster*> fcalClusters;
238 eventLoop->Get(fcalClusters);
239 if(fcalClusters.size()<1)return NOERROR;
240
241 // Use the center of the target as an approximation for the vertex position
242 // 29/03/2020 ijaegle@jlab.org add beam center in x,y
243 DVector3 vertex(m_beamSpotX, m_beamSpotY, m_zTarget);
244
245 vector< const DTrackWireBased* > allWBTracks;
246 eventLoop->Get( allWBTracks );
247 vector< const DTrackWireBased* > wbTracks = filterWireBasedTracks( allWBTracks );
248
249 // Loop over list of DFCALCluster objects and calculate the "Non-linear" corrected
250 // energy and position for each. We'll use a logarithmic energy-weighting to
251 // find the final position and error.
252 for( vector< const DFCALCluster* >::const_iterator clItr = fcalClusters.begin();
253 clItr != fcalClusters.end(); ++clItr ){
254 const DFCALCluster* cluster=*clItr;
255
256 // energy weighted time provides better resolution:
257 double cTime = cluster->getTimeEWeight();
258
259 double zback=m_FCALfront + fcalGeom->blockLength();
260 double c_effective=FCAL_C_EFFECTIVE;
261
262 int channel = cluster->getChannelEmax();
263 DVector2 pos=fcalGeom->positionOnFace(channel);
264 // Check if the cluster is in the insert
265 bool in_insert=fcalGeom->inInsert(channel);
266 if (in_insert){
267 zback=m_insertFront + fcalGeom->insertBlockLength();
268 c_effective=INSERT_C_EFFECTIVE;
269 in_insert=true;
270 }
271
272 // Get corrected energy, position, and errZ
273 double Ecorrected;
274 DVector3 pos_corrected;
275 double errZ;
276 double radius = pos.Mod();
277 int ring_nb = (int) (radius / (5 * k_cm));
278 GetCorrectedEnergyAndPosition( cluster, ring_nb , Ecorrected, pos_corrected, errZ, &vertex,in_insert);
279
280 DVector3 pos_log;
281 GetLogWeightedPosition( cluster, pos_log, Ecorrected, &vertex );
282
283 if (Ecorrected>0.){
284 //up to this point, all times have been times at which light reaches
285 //the back of the detector. Here we correct for the time that it
286 //takes the Cherenkov light to reach the back of the detector
287 //so that the t reported is roughly the time of the shower at the
288 //position pos_corrected
289 cTime -= ( zback - pos_corrected.Z() )/c_effective;
290
291 //Apply time-walk correction/global timing offset
292 cTime += ( timeConst0 + timeConst1 * Ecorrected + timeConst2 * TMath::Power( Ecorrected, 2 ) +
293 timeConst3 * TMath::Power( Ecorrected, 3 ) + timeConst4 * TMath::Power( Ecorrected, 4 ) );
294
295 // Make the DFCALShower object
296 DFCALShower* shower = new DFCALShower;
297
298 shower->setEnergy( Ecorrected );
299 shower->setPosition( pos_corrected );
300 shower->setPosition_log( pos_log );
301 shower->setTime ( cTime );
302
303 if (in_insert==false){
304 FillCovarianceMatrix( shower );
305 }
306 else{
307 // Some guesses for insert resolution, currently hard-coded...
308 double sigx=0.1016/sqrt(Ecorrected)+0.2219;
309 shower->ExyztCovariance(1,1)=sigx*sigx;
310 shower->ExyztCovariance(2,2)=sigx*sigx;
311 shower->ExyztCovariance(0,0)=Ecorrected*Ecorrected*(0.01586/Ecorrected
312 +0.0002342/(Ecorrected*Ecorrected)
313 +1.695e-6);
314 for (unsigned int i=0;i<5;i++){
315 for(unsigned int j=0;j<5;j++){
316 if (i!=j) shower->ExyztCovariance(i,j)=0.;
317 }
318
319 }
320 }
321
322 if( VERBOSE > 2 ){
323 printf("FCAL shower: } E=%f x=%f y=%f z=%f t=%f\n",
324 shower->getEnergy(),shower->getPosition().X(),shower->getPosition().Y(),shower->getPosition().Z(),shower->getTime());
325 printf("FCAL shower: dE=%f dx=%f dy=%f dz=%f dt=%f\n",
326 shower->EErr(),shower->xErr(),shower->yErr(),shower->zErr(),shower->tErr());
327 printf("FCAL shower: Ex=%f Ey=%f Ez=%f Et=%f xy=%f\n",
328 shower->EXcorr(),shower->EYcorr(),shower->EZcorr(),shower->ETcorr(),shower->XYcorr());
329 printf("FCAL shower: xz=%f xt=%f yz=%f yt=%f zt=%f\n",
330 shower->XZcorr(),shower->XTcorr(),shower->YZcorr(),shower->YTcorr(),shower->ZTcorr());
331 }
332
333 // now fill information related to shower shape and nearby
334 // tracks -- useful for splitoff rejection later
335
336 double docaTr = 1E6;
337 double timeTr = 1E6;
338 double xTr = 0;
339 double yTr = 0;
340
341 double flightTime;
342 DVector3 projPos, projMom;
343
344 // find the closest track to the shower -- here we loop over the best FOM
345 // wire-based track for every track candidate not just the ones associated
346 // with the topology
347 for( size_t iTrk = 0; iTrk < wbTracks.size(); ++iTrk ){
348
349 if( !wbTracks[iTrk]->GetProjection( SYS_FCAL, projPos, &projMom, &flightTime ) ) continue;
350
351 // need to swim fcalPos to common z for DOCA calculation -- this really
352 // shouldn't be in the loop if the z-value of projPos doesn't change
353 // with each track
354
355 DVector3 fcalFacePos = ( shower->getPosition() - vertex );
356 fcalFacePos.SetMag( fcalFacePos.Mag() * projPos.Z() / fcalFacePos.Z() );
357
358 double distance = ( fcalFacePos - projPos ).Mag();
359
360 if( distance < docaTr ){
361
362 docaTr = distance;
363 // this is the time from the center of the target to the detector -- to compare with
364 // the FCAL time, one needs to have the t0RF at the center of the target. That
365 // comparison happens at a later stage in the analysis.
366 timeTr = ( wbTracks[iTrk]->position().Z() - vertex.Z() ) / SPEED_OF_LIGHT29.9792458 + flightTime;
367 xTr = projPos.X();
368 yTr = projPos.Y();
369 }
370 }
371
372 shower->setDocaTrack( docaTr );
373 shower->setTimeTrack( timeTr );
374
375 // now compute some variables at the hit level
376
377 vector< const DFCALHit* > fcalHits;
378 cluster->Get( fcalHits );
379 shower->setNumBlocks( fcalHits.size() );
380
381 double e9e25, e1e9;
382 getE1925FromHits( e1e9, e9e25, fcalHits,
383 getMaxHit(cluster->getChannelEmax(),fcalHits) );
384 shower->setE1E9( e1e9 );
385 shower->setE9E25( e9e25 );
386
387 double sumU = 0;
388 double sumV = 0;
389 // if there is no nearest track, the defaults for xTr and yTr will result
390 // in using the beam axis as the directional axis
391 getUVFromHits( sumU, sumV, fcalHits,
392 DVector3( shower->getPosition().X(), shower->getPosition().Y(), 0 ),
393 DVector3( xTr, yTr, 0 ) );
394
395 shower->setSumU( sumU );
396 shower->setSumV( sumV );
397
398 shower->AddAssociatedObject( cluster );
399
400 _data.push_back(shower);
401 }
402 }
403
404 return NOERROR;
405}
406
407//--------------------------------
408// GetCorrectedEnergyAndPosition
409//
410// Non-linear and depth corrections should be fixed within DFCALShower member functions
411//--------------------------------
412 void DFCALShower_factory::GetCorrectedEnergyAndPosition(const DFCALCluster* cluster, int ring_nb, double &Ecorrected, DVector3 &pos_corrected, double &errZ, const DVector3 *vertex,bool in_insert)
413{
414 // Non-linear energy correction are done here
415 //int MAXITER = 1000;
416
417 DVector3 posInCal = cluster->getCentroid();
418
419 float x0 = posInCal.Px();
420 float y0 = posInCal.Py();
421 double Eclust = cluster->getEnergy();
422
423 double Ecutoff = 0;
424 double A = 0;
425 double B = 0;
426 double C = 0;
427 double D = 0;
428 double E = 0;
429 double Egamma = Eclust;
Value stored to 'Egamma' during its initialization is never read
430 Ecorrected = 0;
431
432 // block properties
433 double radiation_length=FCAL_RADIATION_LENGTH;
434 double shower_offset=FCAL_SHOWER_OFFSET;
435 double critical_energy=FCAL_CRITICAL_ENERGY;
436 double zfront=m_FCALfront;
437
438 // Check for presence of insert
439 if (in_insert){
440 radiation_length=INSERT_RADIATION_LENGTH;
441 shower_offset=INSERT_SHOWER_OFFSET;
442 critical_energy=INSERT_CRITICAL_ENERGY;
443 zfront=m_insertFront;
444
445 A=INSERT_PAR1;
446 B=INSERT_PAR2;
447 C=INSERT_PAR3;
448 D=INSERT_PAR4;
449 if (Eclust<D){
450 Egamma=A*Eclust/(1.+B*Eclust);
451 }
452 else Egamma=A*D/(1.+D*B)+C*(Eclust-D);
453 }
454 else{
455 // 06/04/2020 ijaegle@jlab.org allows two different energy dependence correction
456 if (USE_RING_E_CORRECTION && energy_dependence_correction_vs_ring.size()>0){
457 // Method II: PRIMEXD way, correction per ring
458 Egamma=Eclust; // Initialize, before correction
459 int ring_region = -1;
460 if (0 <= ring_nb && ring_nb <= 2)
461 ring_region = 0;
462 else if (3 <= ring_nb && ring_nb <= 4)
463 ring_region = 1;
464 else if (ring_nb == 5)
465 ring_region = 2;
466 else if (6 <= ring_nb && ring_nb <= 7)
467 ring_region = 3;
468 else if (8 <= ring_nb && ring_nb <= 9)
469 ring_region = 4;
470 else if (10 <= ring_nb && ring_nb <= 11)
471 ring_region = 5;
472 else if (12 <= ring_nb && ring_nb <= 17)
473 ring_region = 6;
474 else if (18 <= ring_nb && ring_nb <= 20)
475 ring_region = 7;
476 else if (21 <= ring_nb && ring_nb <= 23)
477 ring_region = 8;
478 if (ring_region != -1) {
479 Egamma = 0;
480 A = energy_dependence_correction_vs_ring[ring_region][0];
481 B = energy_dependence_correction_vs_ring[ring_region][1];
482 C = energy_dependence_correction_vs_ring[ring_region][2];
483 //D = energy_dependence_correction_vs_ring[ring_nb][3];
484 //E = energy_dependence_correction_vs_ring[ring_nb][4];
485 //F = energy_dependence_correction_vs_ring[ring_nb][5];
486 //Egamma = Eclust / (A + B * Eclust + C * pow(Eclust, 2) + D * pow(Eclust, 3) + E * pow(Eclust, 4) + F * pow(Eclust, 5));
487 //Egamma = Eclust / (A + B * Eclust + C * pow(Eclust, 2));
488 Egamma = Eclust / (A - exp(-B * Eclust + C));
489 }
490 // End Correction method II
491 } else {
492 // Method I: IU way, one overall correction
493 Egamma = 0;
494 Ecutoff = cutoff_energy;
495 A = linfit_slope;
496 B = linfit_intercept;
497 C = expfit_param1;
498 D = expfit_param2;
499 E = expfit_param3;
500 // 06/02/2016 Shower Non-linearity Correction by Adesh.
501 // 29/03/2020 ijaegle@jlab.org the linear part correction is applied in (some) data/sim. backward comptability?
502 if ( Eclust <= Ecutoff ) {
503
504 Egamma = Eclust / (A * Eclust + B); // Linear part
505
506 } else {
507 // 29/03/2020 ijaegle@jlab.org this correction is always applied if all C=2 & D=E=0 then Egamma = Eclust
508 // if all C=D=E=0 by mistake then Egamma = - Eclust
509 Egamma = Eclust / (C - exp(-D * Eclust + E)); // Non-linear part
510 }
511 } // End Correction method I
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(int chan_Emax, const vector< const DFCALHit* >& hitVec ) const {
726
727 unsigned int maxIndex = 0;
728
729 for( vector< const DFCALHit* >::const_iterator hit = hitVec.begin();
730 hit != hitVec.end(); ++hit ){
731 if (fcalGeom->channel((**hit).row,(**hit).column)==chan_Emax){
732 maxIndex = hit - hitVec.begin();
733 break;
734 }
735 }
736 return maxIndex;
737}
738
739
740unsigned int
741DFCALShower_factory::getMaxHit( const vector< const DFCALHit* >& hitVec ) const {
742
743 unsigned int maxIndex = 0;
744
745 double eMaxSh = 0;
746
747 for( vector< const DFCALHit* >::const_iterator hit = hitVec.begin();
748 hit != hitVec.end(); ++hit ){
749
750 if( (**hit).E > eMaxSh ){
751
752 eMaxSh = (**hit).E;
753 maxIndex = hit - hitVec.begin();
754 }
755 }
756
757 return maxIndex;
758}
759
760void
761DFCALShower_factory::getUVFromHits( double& sumUSh, double& sumVSh,
762 const vector< const DFCALHit* >& hits,
763 const DVector3& showerVec,
764 const DVector3& trackVec ) const {
765
766 // This method forms an axis pointing from the shower to nearest track
767 // and computes the energy-weighted second moment of the shower along
768 // and perpendicular to this axis. True photons are fairly symmetric
769 // and have similar values of sumU and sumV whereas splitoffs tend
770 // to be asymmetric in these variables.
771
772 DVector3 u = ( showerVec - trackVec ).Unit();
773 DVector3 z( 0, 0, 1 );
774 DVector3 v = u.Cross( z );
775
776 DVector3 hitLoc( 0, 0, 0 );
777
778 sumUSh = 0;
779 sumVSh = 0;
780
781 double sumE = 0;
782
783 for( vector< const DFCALHit* >::const_iterator hit = hits.begin();
784 hit != hits.end(); ++hit ){
785
786 hitLoc.SetX( (**hit).x - showerVec.X() );
787 hitLoc.SetY( (**hit).y - showerVec.Y() );
788
789 sumUSh += (**hit).E * pow( u.Dot( hitLoc ), 2 );
790 sumVSh += (**hit).E * pow( v.Dot( hitLoc ), 2 );
791
792 sumE += (**hit).E;
793 }
794
795 sumUSh /= sumE;
796 sumVSh /= sumE;
797}
798
799void
800DFCALShower_factory::getE1925FromHits( double& e1e9Sh, double& e9e25Sh,
801 const vector< const DFCALHit* >& hits,
802 unsigned int maxIndex ) const {
803
804 double E9 = 0;
805 double E25 = 0;
806
807 const DFCALHit* maxHit = hits[maxIndex];
808
809 for( vector< const DFCALHit* >::const_iterator hit = hits.begin();
810 hit != hits.end(); ++hit ){
811
812 if( fabs( (**hit).x - maxHit->x ) < 4.5 && fabs( (**hit).y - maxHit->y ) < 4.5 )
813 E9 += (**hit).E;
814
815 if( fabs( (**hit).x - maxHit->x ) < 8.5 && fabs( (**hit).y - maxHit->y ) < 8.5 )
816 E25 += (**hit).E;
817 }
818
819 e1e9Sh = maxHit->E/E9;
820 e9e25Sh = E9/E25;
821}
822
823
824vector< const DTrackWireBased* >
825DFCALShower_factory::filterWireBasedTracks( vector< const DTrackWireBased* >& wbTracks ) const {
826
827 vector< const DTrackWireBased* > finalTracks;
828 map< unsigned int, vector< const DTrackWireBased* > > sortedTracks;
829
830 // first sort the wire based tracks into lists with a common candidate id
831 // this means that they all come from the same track in the detector
832
833 for( unsigned int i = 0; i < wbTracks.size(); ++i ){
834
835 unsigned int id = wbTracks[i]->candidateid;
836
837 if( sortedTracks.find( id ) == sortedTracks.end() ){
838
839 sortedTracks[id] = vector< const DTrackWireBased* >();
840 }
841
842 sortedTracks[id].push_back( wbTracks[i] );
843 }
844
845 // now loop through that list of unique tracks and for each set
846 // of wire based tracks, choose the one with the highest FOM
847 // (this is choosing among different particle hypotheses)
848
849 for( map< unsigned int, vector< const DTrackWireBased* > >::const_iterator
850 anId = sortedTracks.begin();
851 anId != sortedTracks.end(); ++anId ){
852
853 double maxFOM = 0;
854 unsigned int bestIndex = 0;
855
856 for( unsigned int i = 0; i < anId->second.size(); ++i ){
857
858 if( anId->second[i]->Ndof < 15 ) continue;
859
860 if( anId->second[i]->FOM > maxFOM ){
861
862 maxFOM = anId->second[i]->FOM;
863 bestIndex = i;
864 }
865 }
866
867 finalTracks.push_back( anId->second[bestIndex] );
868 }
869
870 return finalTracks;
871}
872
873
874void DFCALShower_factory::GetLogWeightedPosition( const DFCALCluster* cluster, DVector3 &pos_log, double Egamma, const DVector3 *vertex )
875{
876
877 DVector3 posInCal = cluster->getCentroid();
878
879 vector<const DFCALHit*> locHitVector;
880 cluster->Get(locHitVector);
881
882 int loc_nhits = (int)locHitVector.size();
883 if( loc_nhits < 1 ) {
884 pos_log = posInCal;
885 return;
886 }
887
888 //------ Loop over hits ------//
889
890 double sW = 0.0;
891 double xpos = 0.0;
892 double ypos = 0.0;
893 double W;
894
895 double ecluster = cluster->getEnergy();
896
897 for( int ih = 0; ih < loc_nhits; ih++ ) {
898
899 const DFCALHit *locHit = locHitVector[ih];
900
901 double xcell = locHit->x;
902 double ycell = locHit->y;
903 double ecell = locHit->E;
904
905 W = log_position_const + log( ecell / ecluster );
906 if( W > 0. ) {
907 sW += W;
908 xpos += xcell * W;
909 ypos += ycell * W;
910 }
911
912 }
913
914 double x1, y1;
915 if( sW ) {
916 x1 = xpos / sW;
917 y1 = ypos / sW;
918 } else {
919 cout << "\nBad Cluster Logged in DFCALShower_factory::GetLogWeightedPosition" << endl;
920 x1 = 0.;
921 y1 = 0.;
922 }
923
924
925 // Shower Depth Corrections (copied from GetCorrectedEnergyAndPosition function)
926
927 if ( Egamma > 0 ) {
928 float dxV = x1 - vertex->X();
929 float dyV = y1 - vertex->Y();
930 float zV = vertex->Z();
931
932 double z0 = m_FCALfront - zV;
933 double zMax = FCAL_RADIATION_LENGTH*(FCAL_SHOWER_OFFSET
934 + log(Egamma/FCAL_CRITICAL_ENERGY));
935
936 double zed = z0;
937 double zed1 = z0 + zMax;
938
939 double r0 = sqrt( dxV*dxV + dyV*dyV );
940
941 int niter;
942 for ( niter=0; niter<100; niter++) {
943 double tt = r0/zed1;
944 zed = z0 + zMax/sqrt( 1 + tt*tt );
945 if ( fabs( (zed-zed1) ) < 0.001) {
946 break;
947 }
948 zed1 = zed;
949 }
950
951 posInCal.SetZ( zed + zV );
952 }
953
954 posInCal.SetX( x1 );
955 posInCal.SetY( y1 );
956
957
958 pos_log = posInCal;
959
960
961 return;
962
963}
964
965