packages/CaloReco/CaloRecoUtility.cc

0001 #include "CaloRecoUtility.h"
0002
0003 #include "BEmcCluster.h"
0004 #include "BEmcRec.h"  // for BEmcRec
0005 #include "BEmcRecCEMC.h"
0006
0007 #include <calobase/RawCluster.h>
0008 #include <calobase/RawTowerDefs.h>  // for decode_index1, decode_index2
0009
0010 #include <ffamodules/CDBInterface.h>
0011
0012 #include <phool/phool.h>  // for PHWHERE
0013
0014 #include <cmath>    // for atan2, log, cos, fabs, sin, sqrt
0015 #include <cstdlib>  // for exit, getenv
0016 #include <iostream>
0017 #include <map>  // for _Rb_tree_const_iterator, operat...
0018 #include <string>
0019 #include <utility>  // for pair
0020 #include <vector>   // for vector
0021
0022 void CaloRecoUtility::ShowerDepthCorrZVertex(RawCluster* clus, float vz)
0023 {
0024   //
0025
0026   float xA = clus->get_x();
0027   float yA = clus->get_y();
0028   float zA = clus->get_z();
0029
0030   float zC = zA;
0031
0032   //      bemcC.CorrectShowerDepth(clus->get_energy(),xt,yt,zt,xC,yC, zC);
0033   // Correction in z
0034   // Just tuned for sim data ... don't fully understand why it works like that
0035
0036   float logE = log(0.1);
0037   if (clus->get_energy() > 0.1)
0038   {
0039     logE = std::log(clus->get_energy());
0040   }
0041
0042   float rA = std::sqrt((xA * xA) + (yA * yA));
0043   //  float theta_twr = GetTowerTheta(xA,yA,zA);
0044   float theta_twr;
0045   if (std::fabs(zA) <= 15)
0046   {
0047     theta_twr = 0;
0048   }
0049   else if (zA > 15)
0050   {
0051     theta_twr = std::atan2(zA - 15, rA);
0052   }
0053   else
0054   {
0055     theta_twr = std::atan2(zA + 15, rA);
0056   }
0057
0058   //  fVz = 0;
0059
0060   float theta_tr = std::atan2(zA - vz, rA);
0061   float L = -1.3 + (0.7 * logE);  // Shower CG in long. direction
0062   float dz = L * std::sin(theta_tr - theta_twr) / std::cos(theta_twr);
0063
0064   dz -= vz * 0.10;
0065
0066   zC = zA - dz;
0067
0068   clus->set_z(zC);
0069 }
0070
0071 void CaloRecoUtility::ProbCorrsZVertex(RawCluster* clus, float vz)
0072 {
0073   if (!_profLoaded)
0074   {
0075     LoadProfile();
0076   }
0077
0078   float vvv[3];
0079   vvv[0] = 0;
0080   vvv[1] = 0;
0081   vvv[2] = vz;
0082
0083   _bemc->SetVertex(vvv);
0084   // construct hitlist, which is just the towers hit in the EmcModule data structure format
0085   //   (for this cluster)
0086
0087   EmcModule vhit;
0088   std::vector<EmcModule> HitList;
0089   HitList.erase(HitList.begin(), HitList.end());
0090   int ich;
0091
0092   RawCluster::TowerConstRange towers = clus->get_towers();
0093   RawCluster::TowerConstIterator toweriter;
0094
0095   for (toweriter = towers.first; toweriter != towers.second; ++toweriter)
0096   {
0097     int ix = RawTowerDefs::decode_index2(toweriter->first);  // index2 is phi in CYL
0098     int iy = RawTowerDefs::decode_index1(toweriter->first);  // index1 is eta in CYL
0099     /*
0100     ix -= BINX0;
0101     iy -= BINY0;
0102
0103     if (ix >= 0 && ix < NBINX && iy >= 0 && iy < NBINY)
0104       {
0105
0106     */
0107     //      ich = iy * NBINX+ ix;
0108     ich = iy * 256 + ix;
0109     // add key field to vhit
0110     vhit.ich = ich;
0111     vhit.amp = toweriter->second;
0112     // vhit.amp = tower->get_energy() * fEnergyNorm;  // !!! Global Calibration
0113     vhit.tof = 0.01;  // not used
0114
0115     HitList.push_back(vhit);
0116
0117     //} // if check NBINS
0118   }
0119
0120   float chi2 = 0;
0121   int ndf = 0;
0122   float prob = _bemc->GetProb(HitList, clus->get_energy(), clus->get_x(), clus->get_y(), clus->get_z(), chi2, ndf);
0123
0124   clus->set_prob(prob);
0125   if (ndf > 0)
0126   {
0127     clus->set_chi2(chi2 / ndf);
0128   }
0129   else
0130   {
0131     clus->set_chi2(0);
0132   }
0133 }
0134
0135 void CaloRecoUtility::LoadProfile()
0136 {
0137   const char* calibroot = getenv("CALIBRATIONROOT");
0138
0139   if (calibroot == nullptr)
0140   {
0141     std::cout << PHWHERE << "CALIBRATIONROOT env var not set" << std::endl;
0142     exit(1);
0143   }
0144   std::string fname_emc_prof = calibroot;
0145   fname_emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0146
0147   std::cout << "CaloRecoUtility:::loading emc_prof from " << fname_emc_prof << std::endl;
0148
0149   std::string url = CDBInterface::instance()->getUrl("EMCPROFILE", fname_emc_prof);
0150   _bemc->LoadProfile(url);
0151
0152   _profLoaded = true;
0153 }
0154
0155 CaloRecoUtility::CaloRecoUtility()
0156   : _bemc(new BEmcRecCEMC())
0157 {
0158
0159
0160   _bemc->SetDim(256, 96);
0161
0162   _bemc->SetTowerThreshold(0.030);
0163
0164   float fProbNoiseParam = 0.04;
0165   _bemc->SetProbNoiseParam(fProbNoiseParam);
0166 }
0167
0168 // this two stupid functions are  only here because of a
0169 // cppcheck warning that should have been suppressed
0170 //  "recommended to have a copy constructor/op= because of
0171 //  dynamic alloc resource"
0172 CaloRecoUtility::CaloRecoUtility(CaloRecoUtility& cru)
0173 {
0174
0175   _profLoaded = false;
0176
0177   if (cru._bemc == nullptr)
0178   {
0179     _bemc = nullptr;
0180     return;
0181   }
0182
0183   _bemc = new BEmcRecCEMC();
0184
0185   _bemc->SetDim(256, 96);
0186
0187   _bemc->SetTowerThreshold(0.030);
0188
0189   float fProbNoiseParam = 0.04;
0190   _bemc->SetProbNoiseParam(fProbNoiseParam);
0191 }
0192
0193 CaloRecoUtility& CaloRecoUtility::operator=(const CaloRecoUtility& cru)
0194 {
0195   if (this == &cru)
0196   {
0197     return *this;
0198   }
0199
0200   _profLoaded = false;
0201   if (cru._bemc == nullptr)
0202   {
0203     _bemc = nullptr;
0204     return *this;
0205   }
0206
0207   _bemc = new BEmcRecCEMC();
0208
0209   _bemc->SetDim(256, 96);
0210
0211   _bemc->SetTowerThreshold(0.030);
0212
0213   float fProbNoiseParam = 0.04;
0214   _bemc->SetProbNoiseParam(fProbNoiseParam);
0215
0216   return *this;
0217 }
0218
0219 CaloRecoUtility::~CaloRecoUtility()
0220 {
0221   delete _bemc;
0222 }