9 #include <boost/format.hpp>
23 , energy_array(nullptr)
24 , theta_array(nullptr)
29 TFile*
f =
new TFile(fname.c_str());
32 std::cout <<
"BEmcProfile: Error when opening profile data file " << fname << std::endl;
38 TH1F* hen =
static_cast<TH1F*
>(f->Get(
"hen"));
41 std::cout <<
"BEmcProfile: Error when loading profile data: hen" << std::endl;
46 TH1F* hth =
static_cast<TH1F*
>(f->Get(
"hth"));
49 std::cout <<
"BEmcProfile: Error when loading profile data: hth" << std::endl;
54 nen = hen->GetNbinsX();
55 nth = hth->GetNbinsX();
58 for (
int i = 0;
i <
nen;
i++)
63 for (
int i = 0;
i <
nth;
i++)
70 std::cout <<
"BEmcProfile: Loading profile data from file " << fname << std::endl;
72 std::cout <<
" " << nen <<
" bins in energy: ";
73 for (
int i = 0;
i <
nen;
i++)
77 std::cout << std::endl;
78 std::cout <<
" " << nth <<
" bins in theta: ";
79 for (
int i = 0;
i <
nth;
i++)
83 std::cout << std::endl;
85 hmean =
new TH1F*[nen * nth *
NP];
87 hr4 =
new TH1F*[nen *
nth];
96 for (
int ie = 0; ie <
nen; ie++)
98 for (
int ip = 0; ip <
NP; ip++)
101 hh =
static_cast<TH1F*
>(f->Get(hname.c_str()));
104 std::cout <<
"BEmcProfile: Could not load histogram " << hname
105 <<
", Error when loading profile data for hmean it = "
106 << it <<
", ie = " << ie <<
"ip = " << ip << std::endl;
110 hmean[ii] =
static_cast<TH1F*
>(hh->Clone());
113 hh =
static_cast<TH1F*
>(f->Get(hname.c_str()));
116 std::cout <<
"BEmcProfile: Could not load histogram " << hname
117 <<
", Error when loading profile data for hsigma it = "
118 << it <<
", ie = " << ie <<
", ip = " << ip << std::endl;
122 hsigma[ii] =
static_cast<TH1F*
>(hh->Clone());
128 hh =
static_cast<TH1F*
>(f->Get(hname.c_str()));
132 std::cout <<
"BEmcProfile: Error when loading profile data for hr4 it = "
133 << it <<
", ie = " << ie << std::endl;
137 hr4[ii2] =
static_cast<TH1F*
>(hh->Clone());
178 int nn = plist->size();
189 int iy0 = -1, iz0 = -1;
191 for (
int i = 0;
i < nn;
i++)
193 ee = (*plist)[
i].amp;
197 ich = (*plist)[
i].ich;
210 for (
int i = 0;
i < nn;
i++)
212 ee = (*plist)[
i].amp;
213 ich = (*plist)[
i].ich;
216 if (ee >
thresh && abs(iz - iz0) <= 3 && abs(iy - iy0) <= 3)
223 float zcg = sz / etot;
224 float ycg = sy / etot;
225 int iz0cg = int(zcg + 0.5);
226 int iy0cg = int(ycg + 0.5);
227 float ddz = fabs(zcg - iz0cg);
228 float ddy = fabs(ycg - iy0cg);
244 float e1, e2, e3, e4;
266 float e1t = (e1 + e2 + e3 + e4) / etot;
267 float e2t = (e1 + e2 - e3 - e4) / etot;
268 float e3t = (e1 - e2 - e3 + e4) / etot;
269 float e4t = (e3) / etot;
275 for (
int ip = 0; ip <
NP; ip++)
277 PredictEnergy(ip, en, theta, phi, ddz, ddy, ep[ip], err[ip]);
284 err[ip] = sqrt(err[ip] * err[ip] + 4 * enoise * enoise / etot / etot);
288 err[ip] = sqrt(err[ip] * err[ip] + 1 * enoise * enoise / etot / etot);
293 chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
294 chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
295 chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
296 chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
299 float prob = TMath::Prob(chi2, ndf);
409 if (ip < 0 || ip >=
NP)
411 std::cout <<
"Error in BEmcProfile::PredictEnergy: profile index = "
412 << ip <<
" but should be from 0 to " <<
NP - 1 << std::endl;
418 std::cout <<
"Error in BEmcProfile::PredictEnergy: energy = "
419 << energy <<
" but should be >0" << std::endl;
424 std::cout <<
"Error in BEmcProfile::PredictEnergy: theta = "
425 << theta <<
" but should be >=0" << std::endl;
429 if (ddz < 0 || ddz > 0.5)
431 std::cout <<
"Error in BEmcProfile::PredictEnergy: ddz = "
432 << ddz <<
" but should be from 0 to 0.5" << std::endl;
435 if (ddy < 0 || ddy > 0.5)
437 std::cout <<
"Error in BEmcProfile::PredictEnergy: ddy = "
438 << ddy <<
" but should be from 0 to 0.5" << std::endl;
485 float rr = sqrt((0.5 - ddz) * (0.5 - ddz) + (0.5 - ddy) * (0.5 - ddy));
503 int ii11 = ip + ie1 *
NP + it1 *
nen *
NP;
504 int ii21 = ip + ie2 * NP + it1 *
nen *
NP;
505 int ii12 = ip + ie1 * NP + it2 *
nen *
NP;
506 int ii22 = ip + ie2 * NP + it2 *
nen *
NP;
512 float pr11 =
hmean[ii11]->GetBinContent(ibin);
513 float pr21 =
hmean[ii21]->GetBinContent(ibin);
514 float prt1 = pr11 + (pr21 - pr11) / (log(en2) - log(en1)) * (log(energy) - log(en1));
520 float er11 =
hsigma[ii11]->GetBinContent(ibin);
521 float er21 =
hsigma[ii21]->GetBinContent(ibin);
522 float ert1 = er11 + (er21 - er11) / (1. / sqrt(en2) - 1. / sqrt(en1)) * (1. / sqrt(energy) - 1. / sqrt(en1));
528 float pr12 =
hmean[ii12]->GetBinContent(ibin);
529 float pr22 =
hmean[ii22]->GetBinContent(ibin);
530 float prt2 = pr12 + (pr22 - pr12) / (log(en2) - log(en1)) * (log(energy) - log(en1));
536 float er12 =
hsigma[ii12]->GetBinContent(ibin);
537 float er22 =
hsigma[ii22]->GetBinContent(ibin);
538 float ert2 = er12 + (er22 - er12) / (1. / sqrt(en2) - 1. / sqrt(en1)) * (1. / sqrt(energy) - 1. / sqrt(en1));
546 float pr = prt1 + (prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
551 float er = ert1 + (ert2 - ert1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
565 if (ibin <
hmean[ii11]->GetNbinsX())
567 if (
hmean[ii11]->GetBinContent(ibin + 1) > 0)
572 float dd = (
hmean[ii11]->GetBinContent(ibin2) -
573 hmean[ii11]->GetBinContent(ibin1)) /
580 er = sqrt(er * er + dd * dd);
596 std::cout <<
"Error in BEmcProfile::PredictEnergyR: energy = " << energy
597 <<
" but should be >0" << std::endl;
602 std::cout <<
"Error in BEmcProfile::PredictEnergyR: theta = " << theta
603 <<
" but should be >=0" << std::endl;
609 std::cout <<
"Error in BEmcProfile::PredictEnergyR: rr = " << rr
610 <<
" but should be >1" << std::endl;
655 int ii11 = ie1 + it1 *
nen;
656 int ii21 = ie2 + it1 *
nen;
657 int ii12 = ie1 + it2 *
nen;
658 int ii22 = ie2 + it2 *
nen;
660 int ibin =
hr4[ii11]->FindBin(rr);
664 float pr11 =
hr4[ii11]->GetBinContent(ibin);
665 float pr21 =
hr4[ii21]->GetBinContent(ibin);
666 float prt1 = pr11 + (pr21 - pr11) / (log(en2) - log(en1)) * (log(energy) - log(en1));
672 float pr12 =
hr4[ii12]->GetBinContent(ibin);
673 float pr22 =
hr4[ii22]->GetBinContent(ibin);
674 float prt2 = pr12 + (pr22 - pr12) / (log(en2) - log(en1)) * (log(energy) - log(en1));
682 float pr = prt1 + (prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
693 int nn = plist->size();
699 for (
int i = 0;
i < nn;
i++)
701 int ich = (*plist)[
i].ich;
704 if (iy == iyt && iz == izt)
706 return (*plist)[
i].amp;