Bug Summary

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