17 #include <ecl/modules/eclShowerCorrection/ECLShowerCorrectorModule.h>
20 #include <framework/logging/Logger.h>
26 #include <ecl/dataobjects/ECLShower.h>
27 #include <ecl/dbobjects/ECLShowerCorrectorLeakageCorrection.h>
28 #include <ecl/dbobjects/ECLShowerEnergyCorrectionTemporary.h>
31 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
46 m_leakageCorrectionPtr_bgx0("ecl_shower_corrector_leakage_corrections"),
47 m_leakageCorrectionPtr_thetaGeo_phase2bgx1("ECLLeakageCorrection_thetaGeometry_phase2"),
48 m_leakageCorrectionPtr_thetaGeo_phase3bgx1("ECLLeakageCorrection_thetaGeometry_phase3"),
49 m_leakageCorrectionPtr_phiGeo_phase2bgx1("ECLLeakageCorrection_phiGeometry_phase2"),
50 m_leakageCorrectionPtr_phiGeo_phase3bgx1("ECLLeakageCorrection_phiGeometry_phase3"),
51 m_leakageCorrectionPtr_thetaEn_phase2bgx1("ECLLeakageCorrection_thetaEnergy_phase2"),
52 m_leakageCorrectionPtr_thetaEn_phase3bgx1("ECLLeakageCorrection_thetaEnergy_phase3"),
53 m_leakageCorrectionPtr_phiEn_phase2bgx1("ECLLeakageCorrection_phiEnergy_phase2"),
54 m_leakageCorrectionPtr_phiEn_phase3bgx1("ECLLeakageCorrection_phiEnergy_phase3"),
55 m_eclShowers(eclShowerArrayName()),
56 m_eventLevelClusteringInfo(eventLevelClusteringInfoName())
60 setDescription(
"ECLShowerCorrectorModule: Corrects for MC truth to reconstruction shower and highest crystal energy differences");
61 setPropertyFlags(c_ParallelProcessingCertified);
72 B2DEBUG(28,
"ECLShowerCorrectorModule::initialize()");
90 double backgroundLevel = 0.0;
101 const double energy = eclShower.getEnergy();
102 const double energyHighest = eclShower.getEnergyHighestCrystal();
103 const double theta = eclShower.getTheta() * TMath::RadToDeg();
104 const double phi = eclShower.getPhi() * TMath::RadToDeg();
107 double correctionFactor = 1.0;
109 if (backgroundLevel < 0.1) correctionFactor =
getLeakageCorrection(theta, phi, energy, backgroundLevel);
112 B2DEBUG(28,
"theta=" << theta <<
", phi=" << phi <<
", E=" << energy <<
", BG=" << backgroundLevel <<
" --> correction factor=" <<
115 if (correctionFactor < 1.e-5 or correctionFactor > 10.) {
116 B2ERROR(
"theta=" << theta <<
", phi=" << phi <<
", E=" << energy <<
" --> Correction factor=" << correctionFactor <<
117 " is very small/too large! Resetting to 1.0.");
118 correctionFactor = 1.0;
121 const double correctedEnergy = energy * correctionFactor;
122 const double correctedEnergyHighest = energyHighest * correctionFactor;
123 B2DEBUG(28,
"Correction factor=" << correctionFactor <<
", correctedEnergy=" << correctedEnergy <<
", correctedEnergyHighest=" <<
124 correctedEnergyHighest);
127 eclShower.setEnergy(correctedEnergy);
128 eclShower.setEnergyHighestCrystal(correctedEnergyHighest);
162 B2DEBUG(28,
"Preparing leakage corrections, thetaGeometry phase3, theta and energy boundaries : theta_min=" <<
182 B2DEBUG(28,
"Preparing leakage corrections, phiGeometry phase3, phi and energy boundaries : phi_min=" <<
202 B2DEBUG(28,
"Preparing leakage corrections, thetaEnergy phase3, theta and energy boundaries : theta_min=" <<
223 B2DEBUG(28,
"Preparing leakage corrections, phiEnergy phase3, phi and energy boundaries : phi_min=" <<
m_leakage_bgx1_limits[7][0]
294 const double background)
298 if (background > 0.75) add_type = 1;
303 const int type_theta_geo = 0 + add_type;
304 const int type_phi_geo = 2 + add_type;
305 const int type_theta_en = 4 + add_type ;
306 const int type_phi_en = 6 + add_type;
308 B2DEBUG(28,
"Index to select payloads: type_theta_geo = " << type_theta_geo <<
" , type_phi_geo = " << type_phi_geo
309 <<
" type_theta_en = " << type_theta_en <<
" , type_phi_en = " << type_phi_en);
311 double theta_clip = theta;
312 double phi_clip = phi;
313 double energy_clip = energy;
315 B2DEBUG(28,
"Cluster info to compute leakage corrections: theta=" << theta_clip <<
", phi=" << phi_clip <<
", E=" << energy_clip);
322 const double corr_theta_geo =
m_leakage_bgx1[type_theta_geo].Interpolate(theta_clip, energy_clip);
323 B2DEBUG(28,
"After corr_theta_geo computation : theta=" << theta_clip <<
", phi=" << phi_clip <<
", E=" << energy_clip <<
324 " , corr_theta_geo : " << corr_theta_geo);
325 theta_clip = theta; energy_clip = energy;
331 const double corr_phi_geo =
m_leakage_bgx1[type_phi_geo].Interpolate(phi_clip, energy_clip);
332 B2DEBUG(28,
"After corr_phi_geo computation : theta=" << theta_clip <<
", phi=" << phi_clip <<
", E=" << energy_clip <<
333 " , corr_phi_geo : " << corr_phi_geo);
334 phi_clip = phi; energy_clip = energy;
340 const double corr_theta_en =
m_leakage_bgx1[type_theta_en].Interpolate(theta_clip, energy_clip);
341 B2DEBUG(28,
"After corr_theta_en computation : theta=" << theta_clip <<
", phi=" << phi_clip <<
", E=" << energy_clip <<
342 " , corr_theta_en : " << corr_theta_en);
343 theta_clip = theta; energy_clip = energy;
349 const double corr_phi_en =
m_leakage_bgx1[type_phi_en].Interpolate(phi_clip, energy_clip);
350 B2DEBUG(28,
"After corr_phi_en computation : theta=" << theta_clip <<
", phi=" << phi_clip <<
", E=" << energy_clip <<
351 " , corr_phi_en : " << corr_phi_en);
358 const double corr = (1 + (corr_theta_geo + corr_phi_geo) / energy_clip) * corr_theta_en * corr_phi_en;
360 B2DEBUG(28,
"Geometrical correction factors: corr_theta_geo=" << corr_theta_geo <<
", corr_phi_geo=" << corr_phi_geo);
361 B2DEBUG(28,
"Energy correction factors: corr_theta_en=" << corr_theta_en <<
", corr_phi_en=" << corr_phi_en);
362 B2DEBUG(28,
"Final correction=" << corr);
373 const double background)
const
377 if (background > 0.1) bf = 1;
386 const double phiMod = phi + 180;
395 if (x <= phiBinWidth / 2.0) {
397 x1 = phiBinWidth / 2.0;
403 if (x >= (ii * phiBinWidth + phiBinWidth / 2.0) and x < ((ii + 1)*phiBinWidth + phiBinWidth / 2.0)) {
404 x0 = (ii * phiBinWidth + phiBinWidth / 2.0);
405 x1 = ((ii + 1) * phiBinWidth + phiBinWidth / 2.0);
418 B2DEBUG(28,
"m_numOfPhiBins=" <<
m_numOfPhiBins <<
" x0Bin=" << x0Bin <<
" x1Bin=" << x1Bin);
428 int reg1thetaBin = 0;
442 y0Bin = reg1thetaBin;
443 y1Bin = reg1thetaBin + 1;
457 int reg2thetaBin = 0;
471 y0Bin = reg2thetaBin;
472 y1Bin = reg2thetaBin + 1;
480 int reg3thetaBin = 0;
494 y0Bin = reg3thetaBin;
495 y1Bin = reg3thetaBin + 1;
509 B2DEBUG(28,
"No valid theta region found (theta=" << theta <<
" deg), return correction factor of 1.0.");
513 B2DEBUG(28,
"y0Bin=" << y0Bin <<
" y1Bin=" << y1Bin);
549 B2DEBUG(28,
"z0Bin=" << z0Bin <<
" z1Bin=" << z1Bin);
551 xd = (x - x0) / (x1 - x0);
552 yd = (y - y0) / (y1 - y0);
553 zd = (z - z0) / (z1 - z0);
556 double c00, c01, c10, c11, c0, c1;
562 c0 = c00 * (1 - yd) + c10 * yd;
563 c1 = c01 * (1 - yd) + c11 * yd;
565 result = c0 * (1 - zd) + c1 * zd;
567 double c00, c01, c10, c11, c0, c1;
573 c0 = c00 * (1 - yd) + c10 * yd;
574 c1 = c01 * (1 - yd) + c11 * yd;
576 result = c0 * (1 - zd) + c1 * zd;
578 float c00, c01, c10, c11, c0, c1;
584 c0 = c00 * (1 - yd) + c10 * yd;
585 c1 = c01 * (1 - yd) + c11 * yd;
587 result = c0 * (1 - zd) + c1 * zd;