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