9 #include <boost/format.hpp>
31 TFile*
f =
new TFile(fname.c_str());
34 cout <<
"BEmcProfile: Error when opening profile data file " << fname << endl;
40 TH1F* hen = (TH1F*) f->Get(
"hen");
43 cout <<
"BEmcProfile: Error when loading profile data: hen" << endl;
48 TH1F* hth = (TH1F*) f->Get(
"hth");
51 cout <<
"BEmcProfile: Error when loading profile data: hth" << endl;
56 nen = hen->GetNbinsX();
57 nth = hth->GetNbinsX();
60 for (
int i = 0; i <
nen; i++)
energy_array[i] = hen->GetBinContent(i + 1);
62 for (
int i = 0; i <
nth; i++)
theta_array[i] = hth->GetBinContent(i + 1);
66 cout <<
"BEmcProfile: Loading profile data from file " << fname << endl;
68 cout <<
" " << nen <<
" bins in energy: ";
69 for (
int i = 0; i <
nen; i++)
74 cout <<
" " << nth <<
" bins in theta: ";
75 for (
int i = 0; i <
nth; i++)
81 hmean =
new TH1F*[nen * nth *
NP];
83 hr4 =
new TH1F*[nen *
nth];
92 for (
int ie = 0; ie <
nen; ie++)
94 for (
int ip = 0; ip <
NP; ip++)
96 hname = boost::str(boost::format(
"hmean%d_en%d_t%d") % (ip + 1) % ie %
it);
97 hh = (TH1F*) f->Get(hname.c_str());
100 cout <<
"BEmcProfile: Could not load histogram " << hname
101 <<
", Error when loading profile data for hmean it = "
102 << it <<
", ie = " << ie <<
"ip = " << ip << endl;
106 hmean[ii] = (TH1F*) hh->Clone();
108 hname = boost::str(boost::format(
"hsigma%d_en%d_t%d") % (ip + 1) % ie % it);
109 hh = (TH1F*) f->Get(hname.c_str());
112 cout <<
"BEmcProfile: Could not load histogram " << hname
113 <<
", Error when loading profile data for hsigma it = "
114 << it <<
", ie = " << ie <<
", ip = " << ip << endl;
118 hsigma[ii] = (TH1F*) hh->Clone();
123 hname = boost::str(boost::format(
"hr4_en%d_t%d") % ie %
it);
124 hh = (TH1F*) f->Get(hname.c_str());
128 cout <<
"BEmcProfile: Error when loading profile data for hr4 it = "
129 << it <<
", ie = " << ie << endl;
133 hr4[ii2] = (TH1F*) hh->Clone();
149 for (
int i = 0; i <
nth *
nen *
NP; i++)
157 for (
int i = 0; i <
nth *
nen; i++)
175 int nn = plist->size();
176 if (nn <= 0)
return -1;
183 int iy0 = -1, iz0 = -1;
185 for (
int i = 0; i <
nn; i++)
187 ee = (*plist)[i].amp;
191 ich = (*plist)[i].ich;
196 if (emax <= 0)
return -1;
201 for (
int i = 0; i <
nn; i++)
203 ee = (*plist)[i].amp;
204 ich = (*plist)[i].ich;
207 if (ee >
thresh &&
abs(iz - iz0) <= 3 &&
abs(iy - iy0) <= 3)
214 float zcg = sz / etot;
215 float ycg = sy / etot;
216 int iz0cg =
int(zcg + 0.5);
217 int iy0cg =
int(ycg + 0.5);
218 float ddz = fabs(zcg - iz0cg);
219 float ddy = fabs(ycg - iy0cg);
222 if (zcg - iz0cg < 0) isz = -1;
224 if (ycg - iy0cg < 0) isy = -1;
239 float e1t = (e1 + e2 + e3 +
e4) / etot;
240 float e2t = (e1 + e2 - e3 -
e4) / etot;
241 float e3t = (e1 - e2 - e3 +
e4) / etot;
242 float e4t = (
e3) / etot;
248 for (
int ip = 0; ip <
NP; ip++)
250 PredictEnergy(ip, en, theta, phi, ddz, ddy, ep[ip], err[ip]);
251 if (ep[ip] < 0)
return -1;
253 err[ip] = sqrt(err[ip] * err[ip] + 4 * enoise * enoise / etot / etot);
255 err[ip] = sqrt(err[ip] * err[ip] + 1 * enoise * enoise / etot / etot);
259 chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
260 chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
261 chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
262 chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
265 float prob = TMath::Prob(chi2, ndf);
375 if (ip < 0 || ip >=
NP)
377 cout <<
"Error in BEmcProfile::PredictEnergy: profile index = "
378 << ip <<
" but should be from 0 to " <<
NP - 1 << endl;
384 cout <<
"Error in BEmcProfile::PredictEnergy: energy = "
385 << energy <<
" but should be >0" << endl;
390 cout <<
"Error in BEmcProfile::PredictEnergy: theta = "
391 << theta <<
" but should be >=0" << endl;
395 if (ddz < 0 || ddz > 0.5)
397 cout <<
"Error in BEmcProfile::PredictEnergy: ddz = "
398 << ddz <<
" but should be from 0 to 0.5" << endl;
401 if (ddy < 0 || ddy > 0.5)
403 cout <<
"Error in BEmcProfile::PredictEnergy: ddy = "
404 << ddy <<
" but should be from 0 to 0.5" << endl;
425 while (it2 < nth && theta >
theta_array[it2]) it2++;
437 float rr = sqrt((0.5 - ddz) * (0.5 - ddz) + (0.5 - ddy) * (0.5 - ddy));
451 int ii11 = ip + ie1 *
NP + it1 *
nen *
NP;
452 int ii21 = ip + ie2 * NP + it1 *
nen *
NP;
453 int ii12 = ip + ie1 * NP + it2 *
nen *
NP;
454 int ii22 = ip + ie2 * NP + it2 *
nen *
NP;
456 int ibin =
hmean[ii11]->FindBin(xx);
460 float pr11 =
hmean[ii11]->GetBinContent(ibin);
461 float pr21 =
hmean[ii21]->GetBinContent(ibin);
462 float prt1 = pr11 + (pr21 - pr11) / (log(en2) - log(en1)) * (log(energy) - log(en1));
463 if (prt1 < 0) prt1 = 0;
465 float er11 =
hsigma[ii11]->GetBinContent(ibin);
466 float er21 =
hsigma[ii21]->GetBinContent(ibin);
467 float ert1 = er11 + (er21 - er11) / (1. / sqrt(en2) - 1. / sqrt(en1)) * (1. / sqrt(energy) - 1. / sqrt(en1));
468 if (ert1 < 0) ert1 = 0;
470 float pr12 =
hmean[ii12]->GetBinContent(ibin);
471 float pr22 =
hmean[ii22]->GetBinContent(ibin);
472 float prt2 = pr12 + (pr22 - pr12) / (log(en2) - log(en1)) * (log(energy) - log(en1));
473 if (prt2 < 0) prt2 = 0;
475 float er12 =
hsigma[ii12]->GetBinContent(ibin);
476 float er22 =
hsigma[ii22]->GetBinContent(ibin);
477 float ert2 = er12 + (er22 - er12) / (1. / sqrt(en2) - 1. / sqrt(en1)) * (1. / sqrt(energy) - 1. / sqrt(en1));
478 if (ert2 < 0) ert2 = 0;
482 float pr = prt1 + (prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
484 float er = ert1 + (ert2 - ert1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
490 if (ibin > 1) ibin1 = ibin - 1;
492 if (ibin <
hmean[ii11]->GetNbinsX())
494 if (
hmean[ii11]->GetBinContent(ibin + 1) > 0)
499 float dd = (
hmean[ii11]->GetBinContent(ibin2) -
500 hmean[ii11]->GetBinContent(ibin1)) /
507 er = sqrt(er * er + dd * dd);
525 cout <<
"Error in BEmcProfile::PredictEnergyR: energy = " << energy
526 <<
" but should be >0" << endl;
531 cout <<
"Error in BEmcProfile::PredictEnergyR: theta = " << theta
532 <<
" but should be >=0" << endl;
538 cout <<
"Error in BEmcProfile::PredictEnergyR: rr = " << rr
539 <<
" but should be >1" << endl;
554 while (it2 < nth && theta >
theta_array[it2]) it2++;
570 int ii11 = ie1 + it1 *
nen;
571 int ii21 = ie2 + it1 *
nen;
572 int ii12 = ie1 + it2 *
nen;
573 int ii22 = ie2 + it2 *
nen;
575 int ibin =
hr4[ii11]->FindBin(rr);
579 float pr11 =
hr4[ii11]->GetBinContent(ibin);
580 float pr21 =
hr4[ii21]->GetBinContent(ibin);
581 float prt1 = pr11 + (pr21 - pr11) / (log(en2) - log(en1)) * (log(energy) - log(en1));
582 if (prt1 < 0) prt1 = 0;
584 float pr12 =
hr4[ii12]->GetBinContent(ibin);
585 float pr22 =
hr4[ii22]->GetBinContent(ibin);
586 float prt2 = pr12 + (pr22 - pr12) / (log(en2) - log(en1)) * (log(energy) - log(en1));
587 if (prt2 < 0) prt2 = 0;
591 float pr = prt1 + (prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
600 int nn = plist->size();
601 if (nn <= 0)
return 0;
603 for (
int i = 0; i <
nn; i++)
605 int ich = (*plist)[i].ich;
608 if (iy == iyt && iz == izt)
610 return (*plist)[i].amp;