50 setDescription(
"Store quantities from single photon MC used to calculated ECL energy leakage corrections");
54 addParam(
"energies_forward",
m_energies_forward,
"generated photon energies, forward", std::vector<double> {0.030, 0.050, 0.100, 0.200, 0.483, 1.166, 2.816, 6.800});
55 addParam(
"energies_barrel",
m_energies_barrel,
"generated photon energies, barrel", std::vector<double> {0.030, 0.050, 0.100, 0.200, 0.458, 1.049, 2.402, 5.500});
56 addParam(
"energies_backward",
m_energies_backward,
"generated photon energies, backward", std::vector<double> {0.030, 0.050, 0.100, 0.200, 0.428, 0.917, 1.962, 4.200});
181 double firstBin = 2.01;
203 if (nMC != 1) {
return;}
211 double minAngle = 999.;
213 for (
int is = 0; is < nShower; is++) {
219 ROOT::Math::XYZVector position;
220 VectorUtil::setMagThetaPhi(
225 ROOT::Math::XYZVector direction = ROOT::Math::XYZVector(position) - vertex;
227 double angle = ROOT::Math::VectorUtil::Angle(direction, mcp3);
228 if (angle < minAngle) {
234 if (minShower == -1) {
return;}
255 const float genEnergyMeV = 1000.*mcLabE;
256 const float tolerance = std::max(0.002 * genEnergyMeV, 1.0);
281 double e3x3 = e3x3Orig;
284 const double eHighestCorr =
m_eclShowerArray[minShower]->getEnergyHighestCrystal();
286 const double eHighestRaw = eHighestCorr / correction;
287 const double e3x3Alt = eHighestRaw /
m_eclShowerArray[minShower]->getE1oE9();
288 if (e3x3 < 0.001) {e3x3 = e3x3Alt;}
306 std::vector<double> digitEnergies;
308 unsigned int nRelatedDigits = showerDigitRelations.size();
309 for (
unsigned int ir = 0; ir < nRelatedDigits; ir++) {
310 const auto calDigit = showerDigitRelations.object(ir);
311 const auto weight = showerDigitRelations.weight(ir);
312 digitEnergies.push_back(calDigit->getEnergy() * weight);
316 std::sort(digitEnergies.begin(), digitEnergies.end());
319 double eSumOfN = 1.e-10;
320 for (
int isum = 0; isum <
t_nCrys; isum++) {
321 const int i = (int)nRelatedDigits - 1 - isum;
322 if (i >= 0) {eSumOfN += digitEnergies[i];}
335 const int iy = iEnergy + 1;
337 const int ixNom = 3 * ig + 1;
338 const int ixLowerE = ixNom + 1;
339 const int ixHigherE = ixNom + 2;
343 const double logEHigher =
m_logPeakEnergy.GetBinContent(ixHigherE, iy);
345 const double biasNom =
m_bias.GetBinContent(ixNom, iy);
346 const double biasLower =
m_bias.GetBinContent(ixLowerE, iy);
347 const double biasHigher =
m_bias.GetBinContent(ixHigherE, iy);
354 const double logESumN = log(eSumOfN);
356 double logEOther = logELower;
357 double biasOther = biasLower;
358 double peakOther = peakLower;
359 if (logESumN > logENom) {
360 logEOther = logEHigher;
361 biasOther = biasHigher;
362 peakOther = peakHigher;
366 double bias = biasNom;
367 double peak = peakNom;
368 if (std::abs(logEOther - logENom) > 0.0001) {
369 bias = biasNom + (biasOther - biasNom) * (logESumN - logENom) / (logEOther - logENom);
370 peak = peakNom + (peakOther - peakNom) * (logESumN - logENom) / (logEOther - logENom);
379 ROOT::Math::XYZVector measuredLocation;
380 VectorUtil::setMagThetaPhi(
381 measuredLocation, radius, thetaLocation, phiLocation);
382 ROOT::Math::XYZVector measuredDirection = ROOT::Math::XYZVector(measuredLocation) - vertex;
383 t_locationError = radius * ROOT::Math::VectorUtil::Angle(measuredDirection, mcp3);
392 B2INFO(
" e3x3Orig " << e3x3Orig <<
" e3x3Alt " << e3x3Alt <<
" Eraw " << eRaw <<
" ESum " << eSumOfN <<
" eFrac " <<
394 B2INFO(
" 3 log E " << logELower <<
" " << logENom <<
" " << logEHigher <<
" log E " << logESumN);
395 B2INFO(
" 3 biases " << biasLower <<
" " << biasNom <<
" " << biasHigher <<
" bias " << bias);
396 B2INFO(
" 3 peaks " << peakLower <<
" " << peakNom <<
" " << peakHigher <<
" peak " << peak);