1 | |
2 | |
3 | |
4 | |
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 | |
25 | |
26 | DFCALShower_factory::DFCALShower_factory() |
27 | { |
28 | |
29 | |
30 | LOAD_NONLIN_CCDB = true; |
31 | LOAD_TIMING_CCDB = true; |
32 | |
33 | gPARMS->SetDefaultParameter("FCAL:LOAD_NONLIN_CCDB", LOAD_NONLIN_CCDB); |
34 | gPARMS->SetDefaultParameter("FCAL:LOAD_TIMING_CCDB", LOAD_TIMING_CCDB); |
35 | |
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 | |
43 | |
44 | |
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 | |
73 | |
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; |
83 | COVARIANCEFILENAME = ""; |
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 | |
109 | |
110 | jerror_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 | |
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 | |
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 | |
153 | |
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 | |
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 | |
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 | |
220 | jerror_t DFCALShower_factory::erun(void) { |
221 | |
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 | |
234 | |
235 | jerror_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 | |
242 | |
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 | |
250 | |
251 | |
252 | for( vector< const DFCALCluster* >::const_iterator clItr = fcalClusters.begin(); |
253 | clItr != fcalClusters.end(); ++clItr ){ |
254 | const DFCALCluster* cluster=*clItr; |
255 | |
256 | |
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 | |
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 | |
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 | |
285 | |
286 | |
287 | |
288 | |
289 | cTime -= ( zback - pos_corrected.Z() )/c_effective; |
290 | |
291 | |
292 | cTime += ( timeConst0 + timeConst1 * Ecorrected + timeConst2 * TMath::Power( Ecorrected, 2 ) + |
293 | timeConst3 * TMath::Power( Ecorrected, 3 ) + timeConst4 * TMath::Power( Ecorrected, 4 ) ); |
294 | |
295 | |
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 | |
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 | |
334 | |
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 | |
345 | |
346 | |
347 | for( size_t iTrk = 0; iTrk < wbTracks.size(); ++iTrk ){ |
348 | |
349 | if( !wbTracks[iTrk]->GetProjection( SYS_FCAL, projPos, &projMom, &flightTime ) ) continue; |
350 | |
351 | |
352 | |
353 | |
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 | |
364 | |
365 | |
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 | |
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 | |
390 | |
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 | |
409 | |
410 | |
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 | |
415 | |
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 | |
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 | |
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 | |
456 | if (USE_RING_E_CORRECTION && energy_dependence_correction_vs_ring.size()>0){ |
457 | |
458 | Egamma=Eclust; |
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 | |
484 | |
485 | |
486 | |
487 | |
488 | Egamma = Eclust / (A - exp(-B * Eclust + C)); |
489 | } |
490 | |
491 | } else { |
492 | |
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 | |
501 | |
502 | if ( Eclust <= Ecutoff ) { |
503 | |
504 | Egamma = Eclust / (A * Eclust + B); |
505 | |
506 | } else { |
507 | |
508 | |
509 | Egamma = Eclust / (C - exp(-D * Eclust + E)); |
510 | } |
511 | } |
512 | } |
513 | |
514 | |
515 | if (Egamma <= 0 && Eclust > 0) Egamma = Eclust; |
516 | |
517 | |
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 | |
555 | |
556 | |
557 | |
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 | |
575 | if (Elookup<minElookup) Elookup=minElookup; |
576 | if (Elookup>maxElookup) Elookup=maxElookup-0.0001; |
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; |
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 | |
624 | if (COVARIANCEFILENAME == "") USECCDB=1; |
625 | |
626 | map<string,string> covariance_data; |
627 | if (USECCDB) { |
628 | |
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) { |
634 | |
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 | |
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 | |
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 | |
680 | CovarianceLookupTable[i][j] = new TH2F(histname,"Covariance histogram",10,0,12,10,0,12); |
681 | CovarianceLookupTable[i][j]->SetDirectory(nullptr); |
682 | } else { |
683 | |
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 | |
704 | CovarianceLookupTable[i][j] = new TH2F(histname,"Covariance histogram",nxbins,xbins,nybins,ybins); |
705 | CovarianceLookupTable[i][j]->SetDirectory(nullptr); |
706 | |
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 | |
715 | ifs.close(); |
716 | } |
717 | savedir->cd(); |
718 | japp->RootUnLock(); |
719 | } |
720 | } |
721 | return NOERROR; |
722 | } |
723 | |
724 | unsigned int |
725 | DFCALShower_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 | |
740 | unsigned int |
741 | DFCALShower_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 | |
760 | void |
761 | DFCALShower_factory::getUVFromHits( double& sumUSh, double& sumVSh, |
762 | const vector< const DFCALHit* >& hits, |
763 | const DVector3& showerVec, |
764 | const DVector3& trackVec ) const { |
765 | |
766 | |
767 | |
768 | |
769 | |
770 | |
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 | |
799 | void |
800 | DFCALShower_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 | |
824 | vector< const DTrackWireBased* > |
825 | DFCALShower_factory::filterWireBasedTracks( vector< const DTrackWireBased* >& wbTracks ) const { |
826 | |
827 | vector< const DTrackWireBased* > finalTracks; |
828 | map< unsigned int, vector< const DTrackWireBased* > > sortedTracks; |
829 | |
830 | |
831 | |
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 | |
846 | |
847 | |
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 | |
874 | void 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 | |
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 | |
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 | |