packages/CaloReco/BEmcProfile.cc

0001 #include "BEmcProfile.h"
0002 #include "BEmcCluster.h"
0003
0004 #include <TFile.h>
0005 #include <TH1.h>  // for TH1F
0006 #include <TMath.h>
0007 #include <TROOT.h>
0008
0009 #include <algorithm>
0010 #include <boost/format.hpp>
0011
0012 #include <cmath>    // for sqrt, log, pow, fabs
0013 #include <cstdlib>  // for abs
0014 #include <iostream>
0015 #include <memory>  // for allocator_traits<>::value_type
0016
0017 const int NP = 4;  // Number of profiles in a bin
0018
0019 BEmcProfile::BEmcProfile(const std::string& fname)
0020   : bloaded(false)
0021   , thresh(0.01)
0022   , nen(0)
0023   , nth(0)
0024   , energy_array(nullptr)
0025   , theta_array(nullptr)
0026   , hmean(nullptr)
0027   , hsigma(nullptr)
0028   , m_Verbosity(0)
0029 {
0030   TFile* f = new TFile(fname.c_str());
0031   if (f->IsZombie())
0032   {
0033     std::cout << "BEmcProfile: Error when opening profile data file " << fname << std::endl;
0034     return;
0035   }
0036
0037   gROOT->cd();
0038
0039   TH1F* hen = dynamic_cast<TH1F*>(f->Get("hen"));
0040   if (!hen)
0041   {
0042     std::cout << "BEmcProfile: Error when loading profile data: hen" << std::endl;
0043     f->Close();
0044     return;
0045   }
0046
0047   TH1F* hth = dynamic_cast<TH1F*>(f->Get("hth"));
0048   if (!hth)
0049   {
0050     std::cout << "BEmcProfile: Error when loading profile data: hth" << std::endl;
0051     f->Close();
0052     return;
0053   }
0054
0055   nen = hen->GetNbinsX();
0056   nth = hth->GetNbinsX();
0057
0058   energy_array = new float[nen];
0059   for (int i = 0; i < nen; i++)
0060   {
0061     energy_array[i] = hen->GetBinContent(i + 1);
0062   }
0063   theta_array = new float[nth];
0064   for (int i = 0; i < nth; i++)
0065   {
0066     theta_array[i] = hth->GetBinContent(i + 1);
0067   }
0068
0069   if (Verbosity())
0070   {
0071     std::cout << "BEmcProfile: Loading profile data from file " << fname << std::endl;
0072
0073     std::cout << " " << nen << " bins in energy: ";
0074     for (int i = 0; i < nen; i++)
0075     {
0076       std::cout << boost::format("%4.1f, ") % energy_array[i];
0077     }
0078     std::cout << std::endl;
0079     std::cout << "  " << nth << " bins in theta:  ";
0080     for (int i = 0; i < nth; i++)
0081     {
0082       std::cout << boost::format("%4.2f, ") % theta_array[i];
0083     }
0084     std::cout << std::endl;
0085   }
0086 // NOLINTNEXTLINE(bugprone-implicit-widening-of-multiplication-result)
0087   hmean = new TH1F*[nen * nth * NP];
0088 // NOLINTNEXTLINE(bugprone-implicit-widening-of-multiplication-result)
0089   hsigma = new TH1F*[nen * nth * NP];
0090 // NOLINTNEXTLINE(bugprone-implicit-widening-of-multiplication-result)
0091   hr4 = new TH1F*[nen * nth];
0092
0093   TH1F* hh;
0094   int ii = 0;
0095   int ii2 = 0;
0096
0097   std::string hname;
0098   for (int it = 0; it < nth; it++)
0099   {
0100     for (int ie = 0; ie < nen; ie++)
0101     {
0102       for (int ip = 0; ip < NP; ip++)
0103       {
0104         hname = boost::str(boost::format("hmean%d_en%d_t%d") % (ip + 1) % ie % it);
0105         hh = dynamic_cast<TH1F*>(f->Get(hname.c_str()));
0106         if (!hh)
0107         {
0108           std::cout << "BEmcProfile: Could not load histogram " << hname
0109                     << ", Error when loading profile data for hmean it = "
0110                     << it << ", ie = " << ie << "ip = " << ip << std::endl;
0111           f->Close();
0112           return;
0113         }
0114         hmean[ii] = dynamic_cast<TH1F*>(hh->Clone());
0115
0116         hname = boost::str(boost::format("hsigma%d_en%d_t%d") % (ip + 1) % ie % it);
0117         hh = dynamic_cast<TH1F*>(f->Get(hname.c_str()));
0118         if (!hh)
0119         {
0120           std::cout << "BEmcProfile: Could not load histogram " << hname
0121                     << ", Error when loading profile data for hsigma it = "
0122                     << it << ", ie = " << ie << ", ip = " << ip << std::endl;
0123           f->Close();
0124           return;
0125         }
0126         hsigma[ii] = dynamic_cast<TH1F*>(hh->Clone());
0127
0128         ii++;
0129       }
0130
0131       hname = boost::str(boost::format("hr4_en%d_t%d") % ie % it);
0132       hh = dynamic_cast<TH1F*>(f->Get(hname.c_str()));
0133
0134       if (!hh)
0135       {
0136         std::cout << "BEmcProfile: Error when loading profile data for hr4 it = "
0137                   << it << ", ie = " << ie << std::endl;
0138         f->Close();
0139         return;
0140       }
0141       hr4[ii2] = dynamic_cast<TH1F*>(hh->Clone());
0142       ii2++;
0143     }
0144   }
0145
0146   f->Close();
0147   bloaded = true;
0148 }
0149
0150 BEmcProfile::~BEmcProfile()
0151 {
0152   if (bloaded)
0153   {
0154     delete[] energy_array;
0155     delete[] theta_array;
0156     for (int i = 0; i < nth * nen * NP; i++)
0157     {
0158       delete hmean[i];
0159       delete hsigma[i];
0160     }
0161     delete[] hmean;
0162     delete[] hsigma;
0163
0164     for (int i = 0; i < nth * nen; i++)
0165     {
0166       delete hr4[i];
0167     }
0168     delete[] hr4;
0169   }
0170 }
0171
0172 float BEmcProfile::GetProb(std::vector<EmcModule>* plist, int NX, float en, float theta, float phi)
0173 {
0174   float enoise = 0.01;  // 10 MeV per tower
0175
0176   if (!bloaded)
0177   {
0178     //    std::cout << "Error in BEmcProfile::GetProb: profiles not loaded" << std::endl;
0179     return -1;
0180   }
0181
0182   int nn = plist->size();
0183   if (nn <= 0)
0184   {
0185     return -1;
0186   }
0187
0188   // z coordinate below means x coordinate
0189
0190   float ee;
0191   int ich;  // iy, iz;
0192
0193   int iy0 = -1;
0194   int iz0 = -1;
0195   float emax = 0;
0196   for (int i = 0; i < nn; i++)
0197   {
0198     ee = (*plist)[i].amp;
0199     if (ee > emax)
0200     {
0201       emax = ee;
0202       ich = (*plist)[i].ich;
0203       iy0 = ich / NX;
0204       iz0 = ich % NX;
0205     }
0206   }
0207   if (emax <= 0)
0208   {
0209     return -1;
0210   }
0211
0212   float etot = 0;
0213   float sz = 0;
0214   float sy = 0;
0215   for (int i = 0; i < nn; i++)
0216   {
0217     ee = (*plist)[i].amp;
0218     ich = (*plist)[i].ich;
0219     int iy = ich / NX;
0220     int iz = ich % NX;
0221     if (ee > thresh && abs(iz - iz0) <= 3 && abs(iy - iy0) <= 3)
0222     {
0223       etot += ee;
0224       sz += iz * ee;
0225       sy += iy * ee;
0226     }
0227   }
0228   float zcg = sz / etot;
0229   float ycg = sy / etot;
0230 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0231   int iz0cg = int(zcg + 0.5);
0232 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0233   int iy0cg = int(ycg + 0.5);
0234   float ddz = std::fabs(zcg - iz0cg);
0235   float ddy = std::fabs(ycg - iy0cg);
0236
0237   int isz = 1;
0238   if (zcg - iz0cg < 0)
0239   {
0240     isz = -1;
0241   }
0242   int isy = 1;
0243   if (ycg - iy0cg < 0)
0244   {
0245     isy = -1;
0246   }
0247
0248   // 4 central towers: 43
0249   //                   12
0250   // Tower 1 - central one
0251   float e1;
0252   float e2;
0253   float e3;
0254   float e4;
0255   e1 = GetTowerEnergy(iy0cg, iz0cg, plist, NX);
0256   e2 = GetTowerEnergy(iy0cg, iz0cg + isz, plist, NX);
0257   e3 = GetTowerEnergy(iy0cg + isy, iz0cg + isz, plist, NX);
0258   e4 = GetTowerEnergy(iy0cg + isy, iz0cg, plist, NX);
0259   if (e1 < thresh)
0260   {
0261     e1 = 0;
0262   }
0263   if (e2 < thresh)
0264   {
0265     e2 = 0;
0266   }
0267   if (e3 < thresh)
0268   {
0269     e3 = 0;
0270   }
0271   if (e4 < thresh)
0272   {
0273     e4 = 0;
0274   }
0275
0276   float e1t = (e1 + e2 + e3 + e4) / etot;
0277   float e2t = (e1 + e2 - e3 - e4) / etot;
0278   float e3t = (e1 - e2 - e3 + e4) / etot;
0279   float e4t = (e3) / etot;
0280   //  float rr = sqrt((0.5-ddz)*(0.5-ddz)+(0.5-ddy)*(0.5-ddy));
0281
0282   // Predicted values
0283   float ep[NP];
0284   float err[NP];
0285   for (int ip = 0; ip < NP; ip++)
0286   {
0287     PredictEnergy(ip, en, theta, phi, ddz, ddy, ep[ip], err[ip]);
0288     if (ep[ip] < 0)
0289     {
0290       return -1;
0291     }
0292     if (ip < 3)
0293     {
0294       err[ip] = std::sqrt((err[ip] * err[ip]) + (4 * enoise * enoise / etot / etot));
0295     }
0296     else
0297     {
0298       err[ip] = std::sqrt((err[ip] * err[ip]) + (1 * enoise * enoise / etot / etot));
0299     }
0300   }
0301
0302   float chi2 = 0.;
0303   chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
0304   chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
0305   chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
0306   chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
0307   int ndf = 4;
0308
0309   float prob = TMath::Prob(chi2, ndf);
0310
0311   return prob;
0312 }
0313
0314 /*
0315 float BEmcProfile::GetProbTest(std::vector<EmcModule>* plist, int NX, float en, float theta, float& test_rr, float& test_et, float& test_ep, float& test_err)
0316 // This is only for test purposes. Can be removed if not needed
0317 {
0318   int nn = plist->size();
0319   if( nn <= 0 ) return -1;
0320
0321   // z coordinate below means x coordinate
0322
0323   float ee;
0324   int ich, iy, iz;
0325
0326   int iy0=-1, iz0=-1;
0327   float emax=0;
0328   for( int i=0; i<nn; i++ ) {
0329     ee = (*plist)[i].amp;
0330     if( ee>emax ) {
0331       emax = ee;
0332       ich = (*plist)[i].ich;
0333       iy0 = ich/NX;
0334       iz0 = ich%NX;
0335     }
0336   }
0337   if( emax<=0 ) return -1;
0338
0339   float etot=0;
0340   float sz=0;
0341   float sy=0;
0342   for( int i=0; i<nn; i++ ) {
0343     ee = (*plist)[i].amp;
0344     ich = (*plist)[i].ich;
0345     iy = ich/NX;
0346     iz = ich%NX;
0347     if( ee>thresh && abs(iz-iz0)<=3 && abs(iy-iy0)<=3 ) {
0348       etot += ee;
0349       sz += iz*ee;
0350       sy += iy*ee;
0351     }
0352   }
0353   float zcg = sz/etot;
0354   float ycg = sy/etot;
0355   int iz0cg = int(zcg+0.5);
0356   int iy0cg = int(ycg+0.5);
0357   float ddz = fabs(zcg-iz0cg);
0358   float ddy = fabs(ycg-iy0cg);
0359
0360   int isz = 1;
0361   if( zcg-iz0cg < 0 ) isz = -1;
0362   int isy = 1;
0363   if( ycg-iy0cg < 0 ) isy = -1;
0364
0365   // 4 central towers: 43
0366   //                   12
0367   // Tower 1 - central one
0368   float e1, e2, e3, e4;
0369   e1 = GetTowerEnergy(iy0cg,    iz0cg,     plist, NX);
0370   e2 = GetTowerEnergy(iy0cg,    iz0cg+isz, plist, NX);
0371   e3 = GetTowerEnergy(iy0cg+isy,iz0cg+isz, plist, NX);
0372   e4 = GetTowerEnergy(iy0cg+isy,iz0cg,     plist, NX);
0373   if( e1<thresh ) e1 = 0;
0374   if( e2<thresh ) e2 = 0;
0375   if( e3<thresh ) e3 = 0;
0376   if( e4<thresh ) e4 = 0;
0377
0378   float e1t = (e1+e2+e3+e4)/etot;
0379   float e2t = (e1+e2-e3-e4)/etot;
0380   float e3t = (e1-e2-e3+e4)/etot;
0381   float e4t = (e3)/etot;
0382   float rr = sqrt((0.5-ddz)*(0.5-ddz)+(0.5-ddy)*(0.5-ddy));
0383
0384   // Predicted values
0385   float ep[NP];
0386   float err[NP];
0387   for( int ip=0; ip<NP; ip++ )
0388   {
0389     PredictEnergy(ip, en, theta, ddz, ddy, ep[ip], err[ip]);
0390   }
0391   float chi2 = 0.;
0392   chi2 += (ep[0]-e1t)*(ep[0]-e1t)/err[0]/err[0];
0393   chi2 += (ep[1]-e2t)*(ep[1]-e2t)/err[1]/err[1];
0394   chi2 += (ep[2]-e3t)*(ep[2]-e3t)/err[2]/err[2];
0395   chi2 += (ep[3]-e4t)*(ep[3]-e4t)/err[3]/err[3];
0396   int ndf = 4;
0397
0398   float prob = TMath::Prob(chi2, ndf);
0399
0400   test_rr = rr;
0401   test_et = e1t;
0402   test_ep = ep[0];
0403   test_err = err[0];
0404
0405   return prob;
0406 }
0407 */
0408
0409 void BEmcProfile::PredictEnergy(int ip, float energy, float theta, float /*phi*/, float ddz, float ddy, float& ep, float& err)
0410 // ip changes from 0 to NP-1, meaning the profile index 1,2,..,NP
0411 {
0412   ep = err = -1;
0413
0414   if (!bloaded)
0415   {
0416     // std::cout << "Error in BEmcProfile::PredictEnergy: profiles not loaded" << std::endl;
0417     return;
0418   }
0419
0420   if (ip < 0 || ip >= NP)
0421   {
0422     std::cout << "Error in BEmcProfile::PredictEnergy: profile index = "
0423               << ip << " but should be from 0 to " << NP - 1 << std::endl;
0424     return;
0425   }
0426
0427   if (energy <= 0)
0428   {
0429     std::cout << "Error in BEmcProfile::PredictEnergy: energy = "
0430               << energy << " but should be >0" << std::endl;
0431     return;
0432   }
0433   if (theta < 0)
0434   {
0435     std::cout << "Error in BEmcProfile::PredictEnergy: theta = "
0436               << theta << " but should be >=0" << std::endl;
0437     return;
0438   }
0439
0440   if (ddz < 0 || ddz > 0.5)
0441   {
0442     std::cout << "Error in BEmcProfile::PredictEnergy: ddz = "
0443               << ddz << " but should be from 0 to 0.5" << std::endl;
0444   }
0445
0446   if (ddy < 0 || ddy > 0.5)
0447   {
0448     std::cout << "Error in BEmcProfile::PredictEnergy: ddy = "
0449               << ddy << " but should be from 0 to 0.5" << std::endl;
0450   }
0451
0452   // Safety margin (slightly away from bin edge)
0453   //  if (ddz < 0.021) ddz = 0.021;
0454   //  if (ddz > 0.489) ddz = 0.489;
0455   //  if (ddy < 0.021) ddy = 0.021;
0456   //  if (ddy > 0.489) ddy = 0.489;
0457
0458   // Energy bin
0459   int ie2 = 0;
0460   while (ie2 < nen && energy > energy_array[ie2])
0461   {
0462     ie2++;
0463   }
0464   if (ie2 == 0)
0465   {
0466     ie2 = 1;
0467   }
0468   else if (ie2 >= nen)
0469   {
0470     ie2 = nen - 1;
0471   }
0472   int ie1 = ie2 - 1;
0473   //  int ie1 = ie2-2; // For a test()
0474
0475   // Theta bin
0476   int it2 = 0;
0477   while (it2 < nth && theta > theta_array[it2])
0478   {
0479     it2++;
0480   }
0481   if (it2 == 0)
0482   {
0483     it2 = 1;
0484   }
0485   else if (it2 >= nth)
0486   {
0487     it2 = nth - 1;
0488   }
0489   int it1 = it2 - 1;
0490   //  int it1 = it2-2; // For a test()
0491
0492   //  std::cout << "Energy bin= " << ie1 << " " << ie2 << " ("
0493   //       << energy << ")  Theta bin= " << it1 << " " << it2
0494   //       << " (" << theta << ")" << std::endl;
0495
0496   float rr = sqrt(((0.5 - ddz) * (0.5 - ddz)) + ((0.5 - ddy) * (0.5 - ddy)));
0497
0498   float xx = rr;
0499   if (ip == 1)
0500   {
0501     xx = ddy;
0502   }
0503   else if (ip == 2)
0504   {
0505     xx = ddz;
0506   }
0507
0508   float en1 = energy_array[ie1];
0509   float en2 = energy_array[ie2];
0510   float th1 = theta_array[it1];
0511   float th2 = theta_array[it2];
0512
0513   // 1st index - ie, second index - it
0514   int ii11 = ip + (ie1 * NP) + (it1 * nen * NP);
0515   int ii21 = ip + (ie2 * NP) + (it1 * nen * NP);
0516   int ii12 = ip + (ie1 * NP) + (it2 * nen * NP);
0517   int ii22 = ip + (ie2 * NP) + (it2 * nen * NP);
0518
0519   int ibin = hmean[ii11]->FindBin(xx);
0520
0521   // Log (1/sqrt) energy dependence of mean (sigma)
0522   //
0523   float pr11 = hmean[ii11]->GetBinContent(ibin);
0524   float pr21 = hmean[ii21]->GetBinContent(ibin);
0525   float prt1 = pr11 + ((pr21 - pr11) / (std::log(en2) - std::log(en1)) * (std::log(energy) - std::log(en1)));
0526   prt1 = std::max<float>(prt1, 0);
0527
0528   float er11 = hsigma[ii11]->GetBinContent(ibin);
0529   float er21 = hsigma[ii21]->GetBinContent(ibin);
0530   float ert1 = er11 + ((er21 - er11) / (1. / std::sqrt(en2) - 1. / std::sqrt(en1)) * (1. / std::sqrt(energy) - 1. / std::sqrt(en1)));
0531   ert1 = std::max<float>(ert1, 0);
0532
0533   float pr12 = hmean[ii12]->GetBinContent(ibin);
0534   float pr22 = hmean[ii22]->GetBinContent(ibin);
0535   float prt2 = pr12 + ((pr22 - pr12) / (std::log(en2) - std::log(en1)) * (std::log(energy) - std::log(en1)));
0536   prt2 = std::max<float>(prt2, 0);
0537
0538   float er12 = hsigma[ii12]->GetBinContent(ibin);
0539   float er22 = hsigma[ii22]->GetBinContent(ibin);
0540   float ert2 = er12 + ((er22 - er12) / (1. / std::sqrt(en2) - 1. / std::sqrt(en1)) * (1. / std::sqrt(energy) - 1. / std::sqrt(en1)));
0541   ert2 = std::max<float>(ert2, 0);
0542
0543   // Quadratic theta dependence of mean and sigma
0544   //
0545   float pr = prt1 + ((prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2)));
0546   pr = std::max<float>(pr, 0);
0547   float er = ert1 + ((ert2 - ert1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2)));
0548   er = std::max<float>(er, 0);
0549
0550   // Additional error due to binning in xx
0551   //
0552   int ibin1 = ibin;
0553   if (ibin > 1)
0554   {
0555     ibin1 = ibin - 1;
0556   }
0557   int ibin2 = ibin;
0558   if (ibin < hmean[ii11]->GetNbinsX())
0559   {
0560     if (hmean[ii11]->GetBinContent(ibin + 1) > 0)
0561     {
0562       ibin2 = ibin + 1;
0563     }
0564   }
0565   float dd = (hmean[ii11]->GetBinContent(ibin2) -
0566               hmean[ii11]->GetBinContent(ibin1)) /
0567              2.;
0568   //  if( fabs(dd)>er )
0569   // {
0570   //   std::cout << "ie = " << ie1 << ", it = " << it1 << ", bin = "
0571   //     << ibin << ": " << er << " " << dd << std::endl;
0572   // }
0573   er = std::sqrt((er * er) + (dd * dd));
0574
0575   ep = pr;
0576   err = er;
0577 }
0578
0579 float BEmcProfile::PredictEnergyR(float energy, float theta, float /*phi*/, float rr)
0580 {
0581   if (!bloaded)
0582   {
0583     //    std::cout << "Error in BEmcProfile::PredictEnergyR: profiles not loaded" << std::endl;
0584     return 0;
0585   }
0586
0587   if (energy <= 0)
0588   {
0589     std::cout << "Error in BEmcProfile::PredictEnergyR: energy = " << energy
0590               << " but should be >0" << std::endl;
0591     return 0;
0592   }
0593   if (theta < 0)
0594   {
0595     std::cout << "Error in BEmcProfile::PredictEnergyR: theta = " << theta
0596               << " but should be >=0" << std::endl;
0597     return 0;
0598   }
0599
0600   if (rr < 1)
0601   {
0602     std::cout << "Error in BEmcProfile::PredictEnergyR: rr = " << rr
0603               << " but should be >1" << std::endl;
0604   }
0605
0606   // Energy bin
0607   int ie2 = 0;
0608   while (ie2 < nen && energy > energy_array[ie2])
0609   {
0610     ie2++;
0611   }
0612   if (ie2 == 0)
0613   {
0614     ie2 = 1;
0615   }
0616   else if (ie2 >= nen)
0617   {
0618     ie2 = nen - 1;
0619   }
0620   int ie1 = ie2 - 1;
0621   //  int ie1 = ie2-2; // For a test()
0622
0623   // Theta bin
0624   int it2 = 0;
0625   while (it2 < nth && theta > theta_array[it2])
0626   {
0627     it2++;
0628   }
0629   if (it2 == 0)
0630   {
0631     it2 = 1;
0632   }
0633   else if (it2 >= nth)
0634   {
0635     it2 = nth - 1;
0636   }
0637   int it1 = it2 - 1;
0638   //  int it1 = it2-2; // For a test()
0639
0640   //  printf("Energy bin= %d %d (%f)  Theta bin= %d %d (%f)\n",ie1,ie2,energy,it1,it2,theta);
0641
0642   float en1 = energy_array[ie1];
0643   float en2 = energy_array[ie2];
0644   float th1 = theta_array[it1];
0645   float th2 = theta_array[it2];
0646
0647   // 1st index - ie, second index - it
0648   int ii11 = ie1 + (it1 * nen);
0649   int ii21 = ie2 + (it1 * nen);
0650   int ii12 = ie1 + (it2 * nen);
0651   int ii22 = ie2 + (it2 * nen);
0652
0653   int ibin = hr4[ii11]->FindBin(rr);
0654
0655   // Log (1/sqrt) energy dependence of mean (sigma)
0656   //
0657   float pr11 = hr4[ii11]->GetBinContent(ibin);
0658   float pr21 = hr4[ii21]->GetBinContent(ibin);
0659   float prt1 = pr11 + ((pr21 - pr11) / (std::log(en2) - std::log(en1)) * (std::log(energy) - std::log(en1)));
0660   prt1 = std::max<float>(prt1, 0);
0661
0662   float pr12 = hr4[ii12]->GetBinContent(ibin);
0663   float pr22 = hr4[ii22]->GetBinContent(ibin);
0664   float prt2 = pr12 + ((pr22 - pr12) / (std::log(en2) - std::log(en1)) * (std::log(energy) - std::log(en1)));
0665   prt2 = std::max<float>(prt2, 0);
0666
0667   // Quadratic theta dependence of mean and sigma
0668   //
0669   float pr = prt1 + ((prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2)));
0670   pr = std::max<float>(pr, 0);
0671
0672   return pr;
0673 }
0674
0675 float BEmcProfile::GetTowerEnergy(int iy, int iz, std::vector<EmcModule>* plist, int NX)
0676 {
0677   int nn = plist->size();
0678   if (nn <= 0)
0679   {
0680     return 0;
0681   }
0682
0683   for (int i = 0; i < nn; i++)
0684   {
0685     int ich = (*plist)[i].ich;
0686     int iyt = ich / NX;
0687     int izt = ich % NX;
0688     if (iy == iyt && iz == izt)
0689     {
0690       return (*plist)[i].amp;
0691     }
0692   }
0693   return 0;
0694 }