packages/CaloReco/RawClusterPositionCorrection.cc

0001 #include "RawClusterPositionCorrection.h"
0002
0003 #include <calobase/RawCluster.h>
0004 #include <calobase/RawClusterContainer.h>
0005 #include <calobase/RawClusterUtility.h>
0006 #include <calobase/RawTower.h>
0007 #include <calobase/RawTowerContainer.h>
0008 #include <calobase/RawTowerDefs.h>  // for decode_index1, decode_in...
0009 #include <calobase/RawTowerGeomContainer.h>
0010 #include <calobase/TowerInfo.h>
0011 #include <calobase/TowerInfoContainer.h>
0012
0013 #include <cdbobjects/CDBHistos.h>  // for CDBHistos
0014 #include <cdbobjects/CDBTTree.h>   // for CDBTTree
0015
0016 #include <ffamodules/CDBInterface.h>
0017
0018 #include <phparameter/PHParameters.h>
0019
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021 #include <fun4all/SubsysReco.h>
0022
0023 #include <phool/PHCompositeNode.h>
0024 #include <phool/PHIODataNode.h>
0025 #include <phool/PHNode.h>
0026 #include <phool/PHNodeIterator.h>
0027 #include <phool/PHObject.h>
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>
0030
0031 #include <TH1.h>
0032 #include <TH2.h>
0033 #include <TSystem.h>
0034
0035 #include <algorithm>  // for max
0036
0037 #include <cassert>
0038 #include <cmath>
0039 #include <fstream>
0040 #include <iostream>
0041 #include <map>
0042 #include <sstream>
0043 #include <stdexcept>
0044 #include <string>
0045 #include <utility>  // for pair
0046
0047 RawClusterPositionCorrection::RawClusterPositionCorrection(const std::string &name)
0048   : SubsysReco(std::string("RawClusterPositionCorrection_") + name)
0049   , _det_name(name)
0050   , bins_eta(384)
0051   , bins_phi(64)
0052   , iEvent(0)
0053 {
0054 }
0055 RawClusterPositionCorrection::~RawClusterPositionCorrection()
0056 {
0057   delete cdbHisto;
0058   delete cdbttree;
0059 }
0060
0061 int RawClusterPositionCorrection::InitRun(PHCompositeNode *topNode)
0062 {
0063   CreateNodeTree(topNode);
0064
0065   // access the cdb and get cdbtree
0066   std::string m_calibName = "cemc_PDC_NorthSouth_8x8_23instru";
0067   std::string calibdir = CDBInterface::instance()->getUrl(m_calibName);
0068   if (!calibdir.empty())
0069   {
0070     cdbttree = new CDBTTree(calibdir);
0071   }
0072   else
0073   {
0074     std::cout << std::endl
0075               << "did not find CDB tree" << std::endl;
0076     gSystem->Exit(1);
0077     exit(1);
0078   }
0079
0080   h2NorthSector = new TH2F("h2NorthSector", "Cluster; towerid #eta; towerid #phi", bins_eta, 47.5, 95.5, bins_phi, -0.5, 7.5);
0081   h2SouthSector = new TH2F("h2SouthSector", "Cluster; towerid #eta; towerid #phi", bins_eta, -0.5, 47.5, bins_phi, -0.5, 7.5);
0082
0083   /// north
0084   std::string m_fieldname = "cemc_PDC_NorthSector_8x8_clusE";
0085   std::string m_fieldname_ecore = "cemc_PDC_NorthSector_8x8_clusEcore";
0086   float calib_constant = 0;
0087
0088   // Read in the calibration factors and store in the array
0089   for (int i = 0; i < bins_phi; ++i)
0090   {
0091     std::vector<float> dumvec;
0092     std::vector<float> dumvec2;
0093     for (int j = 0; j < bins_eta; ++j)
0094     {
0095       int key = (i * bins_eta) + j;
0096       calib_constant = cdbttree->GetFloatValue(key, m_fieldname);
0097       dumvec.push_back(calib_constant);
0098       calib_constant = cdbttree->GetFloatValue(key, m_fieldname_ecore);
0099       dumvec2.push_back(calib_constant);
0100     }
0101     calib_constants_north.push_back(dumvec);
0102     calib_constants_north_ecore.push_back(dumvec2);
0103   }
0104   /// south
0105   m_fieldname = "cemc_PDC_SouthSector_8x8_clusE";
0106   m_fieldname_ecore = "cemc_PDC_SouthSector_8x8_clusEcore";
0107
0108   // Read in the calibration factors and store in the array
0109   for (int i = 0; i < bins_phi; ++i)
0110   {
0111     std::vector<float> dumvec;
0112     std::vector<float> dumvec2;
0113     for (int j = 0; j < bins_eta; ++j)
0114     {
0115       int key = (i * bins_eta) + j;
0116       calib_constant = cdbttree->GetFloatValue(key, m_fieldname);
0117       dumvec.push_back(calib_constant);
0118       calib_constant = cdbttree->GetFloatValue(key, m_fieldname_ecore);
0119       dumvec2.push_back(calib_constant);
0120     }
0121     calib_constants_south.push_back(dumvec);
0122     calib_constants_south_ecore.push_back(dumvec2);
0123   }
0124
0125   // Load PDC final stage correction
0126   calibdir = CDBInterface::instance()->getUrl("cemc_PDC_ResidualCorr");
0127   if (!calibdir.empty())
0128   {
0129     cdbHisto = new CDBHistos(calibdir);
0130     cdbHisto->LoadCalibrations();
0131     //pdcCorrFlat = cdbHisto->getHisto("h1_res_p");
0132     pdcCorrFlat = cdbHisto->getHisto("h_res_E_eta");
0133   }
0134   else
0135   {
0136     std::cout << std::endl
0137               << "did not find CDB histo" << std::endl;
0138     gSystem->Exit(1);
0139     exit(1);
0140   }
0141   return Fun4AllReturnCodes::EVENT_OK;
0142 }
0143
0144 int RawClusterPositionCorrection::process_event(PHCompositeNode *topNode)
0145 {
0146   // make sure new cluster container is empty
0147   _recalib_clusters->Reset();
0148
0149   if (Verbosity() >= Fun4AllBase::VERBOSITY_SOME)
0150   {
0151     if (iEvent % 100 == 0) { std::cout << "Progress: " << iEvent << std::endl;
0152 }
0153     ++iEvent;
0154   }
0155
0156   std::string rawClusNodeName = "CLUSTER_" + _det_name;
0157   if (m_UseTowerInfo)
0158   {
0159     rawClusNodeName = "CLUSTERINFO_" + _det_name;
0160   }
0161
0162   RawClusterContainer *rawclusters = findNode::getClass<RawClusterContainer>(topNode, rawClusNodeName);
0163   if (!rawclusters)
0164   {
0165     std::cout << "No " << _det_name << " Cluster Container found while in RawClusterPositionCorrection, can't proceed!!!" << std::endl;
0166     return Fun4AllReturnCodes::ABORTEVENT;
0167   }
0168
0169   RawTowerContainer *_towers = nullptr;
0170   TowerInfoContainer *_towerinfos = nullptr;
0171
0172   if (!m_UseTowerInfo)
0173   {
0174     _towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_" + _det_name);
0175     if (!_towers)
0176     {
0177       std::cout << "No calibrated " << _det_name << " tower info found while in RawClusterPositionCorrection, can't proceed!!!" << std::endl;
0178       return Fun4AllReturnCodes::ABORTEVENT;
0179     }
0180   }
0181   else
0182   {
0183     std::string towerinfoNodename = "TOWERINFO_CALIB_" + _det_name;
0184
0185     _towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerinfoNodename);
0186     if (!_towerinfos)
0187     {
0188       std::cerr << Name() << "::" << _det_name << "::" << __PRETTY_FUNCTION__
0189                 << " " << towerinfoNodename << " Node missing, doing bail out!"
0190                 << std::endl;
0191
0192       return Fun4AllReturnCodes::DISCARDEVENT;
0193     }
0194   }
0195
0196   std::string towergeomnodename = "TOWERGEOM_" + _det_name;
0197   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0198   if (!towergeom)
0199   {
0200     std::cout << PHWHERE << ": Could not find node " << towergeomnodename << std::endl;
0201     return Fun4AllReturnCodes::ABORTEVENT;
0202   }
0203   const int nphibin = towergeom->get_phibins();
0204
0205   // loop over the clusters
0206   RawClusterContainer::ConstRange begin_end = rawclusters->getClusters();
0207   RawClusterContainer::ConstIterator iter;
0208
0209   for (iter = begin_end.first; iter != begin_end.second; ++iter)
0210   {
0211     //    RawClusterDefs::keytype key = iter->first;
0212     RawCluster *cluster = iter->second;
0213
0214     float clus_energy = cluster->get_energy();
0215
0216     std::vector<float> toweretas;
0217     std::vector<float> towerphis;
0218     std::vector<float> towerenergies;
0219
0220     // loop over the towers in the cluster
0221     RawCluster::TowerConstRange towers = cluster->get_towers();
0222     RawCluster::TowerConstIterator toweriter;
0223     for (toweriter = towers.first;
0224          toweriter != towers.second;
0225          ++toweriter)
0226     {
0227       if (!m_UseTowerInfo)
0228       {
0229         RawTower *tower = _towers->getTower(toweriter->first);
0230
0231         int towereta = tower->get_bineta();
0232         int towerphi = tower->get_binphi();
0233         double towerenergy = tower->get_energy();
0234
0235         // put the etabin, phibin, and energy into the corresponding vectors
0236         toweretas.push_back(towereta);
0237         towerphis.push_back(towerphi);
0238         towerenergies.push_back(towerenergy);
0239       }
0240       else
0241       {
0242         // std::cout << "clus " << iter->first << " tow key " << toweriter->first << " decoded to " << towerindex << std::endl;
0243
0244         int iphi = RawTowerDefs::decode_index2(toweriter->first);  // index2 is phi in CYL
0245         int ieta = RawTowerDefs::decode_index1(toweriter->first);  // index1 is eta in CYL
0246         unsigned int towerkey = iphi + ((unsigned int) (ieta) << 16U);
0247
0248         assert(_towerinfos);
0249
0250         unsigned int towerindex = _towerinfos->decode_key(towerkey);
0251
0252         TowerInfo *towinfo = _towerinfos->get_tower_at_channel(towerindex);
0253
0254         double towerenergy = towinfo->get_energy();
0255
0256         // put the eta, phi, energy into corresponding vectors
0257         toweretas.push_back(ieta);
0258         towerphis.push_back(iphi);
0259         towerenergies.push_back(towerenergy);
0260       }
0261     }
0262
0263     int ntowers = toweretas.size();
0264     //    std::cout << "jf " <<  ntowers << std::endl;
0265     assert(ntowers >= 1);
0266
0267     // loop over the towers to determine the energy
0268     // weighted eta and phi position of the cluster
0269
0270     float etamult = 0;
0271     float etasum = 0;
0272     float phimult = 0;
0273     float phisum = 0;
0274
0275     for (int j = 0; j < ntowers; j++)
0276     {
0277       float energymult = towerenergies.at(j) * toweretas.at(j);
0278       etamult += energymult;
0279       etasum += towerenergies.at(j);
0280
0281       int phibin = towerphis.at(j);
0282
0283       if (phibin - towerphis.at(0) < -nphibin / 2.0)
0284       {
0285         phibin += nphibin;
0286       }
0287       else if (phibin - towerphis.at(0) > +nphibin / 2.0)
0288       {
0289         phibin -= nphibin;
0290       }
0291       assert(std::abs(phibin - towerphis.at(0)) <= nphibin / 2.0);
0292
0293       energymult = towerenergies.at(j) * phibin;
0294       phimult += energymult;
0295       phisum += towerenergies.at(j);
0296     }
0297
0298     float avgphi = phimult / phisum;
0299     if (std::isnan(avgphi)) { continue;
0300 }
0301
0302     float avgeta = etamult / etasum;
0303
0304     if (avgphi < 0)
0305     {
0306       avgphi += nphibin;
0307     }
0308
0309     avgphi = fmod(avgphi, nphibin);
0310
0311     if (avgphi >= 255.5) { avgphi -= bins_phi;
0312 }
0313
0314     avgphi = fmod(avgphi + 0.5, 8) - 0.5;  // wrapping [-0.5, 255.5] to [-0.5, 7.5]
0315     int etabin = -99;
0316     int phibin = -99;
0317
0318     // check if the cluster is in the north or south sector
0319     if (avgeta < 47.5)
0320     {
0321       etabin = h2SouthSector->GetXaxis()->FindBin(avgeta) - 1;
0322     }
0323     else
0324     {
0325       etabin = h2NorthSector->GetXaxis()->FindBin(avgeta) - 1;
0326     }
0327     phibin = h2NorthSector->GetYaxis()->FindBin(avgphi) - 1;  // can use either h2NorthSector or h2SouthSector since both have the same phi binning
0328
0329     if ((phibin < 0 || etabin < 0) && Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
0330     {
0331       std::cout << "couldn't recalibrate cluster, something went wrong??" << std::endl;
0332     }
0333
0334     float eclus_recalib_val = 1;
0335     float ecore_recalib_val = 1;
0336     if (phibin > -1 && etabin > -1)
0337     {
0338       if (avgeta < 47.5)
0339       {
0340         eclus_recalib_val = calib_constants_south[phibin][etabin];
0341         ecore_recalib_val = calib_constants_south_ecore[phibin][etabin];
0342       }
0343       else
0344       {
0345         eclus_recalib_val = calib_constants_north[phibin][etabin];
0346         ecore_recalib_val = calib_constants_north_ecore[phibin][etabin];
0347       }
0348     }
0349     RawCluster *recalibcluster = dynamic_cast<RawCluster *>(cluster->CloneMe());
0350
0351     assert(recalibcluster);
0352     //    if (m_UseTowerInfo)
0353     recalibcluster->set_energy(clus_energy / eclus_recalib_val);
0354     recalibcluster->set_ecore(cluster->get_ecore() / ecore_recalib_val);
0355
0356     CLHEP::Hep3Vector vertex(0,0,0);
0357     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recalibcluster, vertex);
0358     float clusEta = E_vec_cluster.pseudoRapidity();
0359
0360     if (cluster->get_ecore()    >= pdcCorrFlat->GetXaxis()->GetXmin()
0361         && cluster->get_ecore() <  pdcCorrFlat->GetXaxis()->GetXmax()
0362         && clusEta              >= pdcCorrFlat->GetYaxis()->GetXmin()
0363         && clusEta              <  pdcCorrFlat->GetYaxis()->GetXmax()
0364        )
0365     {
0366
0367       int ecoreBin = pdcCorrFlat->GetXaxis()->FindBin(recalibcluster->get_ecore());
0368       int etaBin = pdcCorrFlat -> GetYaxis() -> FindBin(clusEta);
0369       float pdcCalib =  pdcCorrFlat -> GetBinContent(ecoreBin, etaBin);
0370       //float pdcCalib = pdcCorrFlat->GetBinContent(ecoreBin);
0371       if (pdcCalib < 0.1) { pdcCalib = 1;
0372 }
0373
0374       recalibcluster->set_ecore(recalibcluster->get_ecore() / pdcCalib);
0375     }
0376
0377     _recalib_clusters->AddCluster(recalibcluster);
0378
0379     if (Verbosity() >= Fun4AllBase::VERBOSITY_EVEN_MORE && clus_energy > 1)
0380     {
0381       std::cout << "Input eclus cluster energy: " << clus_energy << std::endl;
0382       std::cout << "Recalib value: " << eclus_recalib_val << std::endl;
0383       std::cout << "phibin: " << phibin << ", etabin: " << etabin << std::endl;
0384       std::cout << "Recalibrated eclus cluster energy: "
0385                 << clus_energy / eclus_recalib_val << std::endl;
0386       std::cout << "Input ecore cluster energy: "
0387                 << cluster->get_ecore() << std::endl;
0388       std::cout << "Recalib value: " << ecore_recalib_val << std::endl;
0389       std::cout << "Recalibrated ecore cluster energy: "
0390                 << cluster->get_ecore() / ecore_recalib_val << std::endl;
0391     }
0392   }
0393
0394   return Fun4AllReturnCodes::EVENT_OK;
0395 }
0396
0397
0398 void RawClusterPositionCorrection::CreateNodeTree(PHCompositeNode *topNode)
0399 {
0400   PHNodeIterator iter(topNode);
0401
0402   // Get the DST Node
0403   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0404
0405   // Check that it is there
0406   if (!dstNode)
0407   {
0408     std::cout << "DST Node missing, quitting" << std::endl;
0409     throw std::runtime_error("failed to find DST node in RawClusterPositionCorrection::CreateNodeTree");
0410   }
0411
0412   // Get the _det_name subnode
0413   PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _det_name));
0414
0415   // Check that it is there
0416   if (!cemcNode)
0417   {
0418     cemcNode = new PHCompositeNode(_det_name);
0419     dstNode->addNode(cemcNode);
0420   }
0421
0422   // Check to see if the cluster recalib node is on the nodetree
0423   std::string ClusterCorrNodeName = "CLUSTER_POS_COR_" + _det_name;
0424   if (m_UseTowerInfo)
0425   {
0426     ClusterCorrNodeName = "CLUSTERINFO_POS_COR_" + _det_name;
0427   }
0428   _recalib_clusters = findNode::getClass<RawClusterContainer>(topNode, ClusterCorrNodeName);
0429   if (_recalib_clusters)
0430   {
0431     _recalib_clusters->Clear();
0432   }
0433   else
0434   {
0435     _recalib_clusters = new RawClusterContainer();
0436     PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_recalib_clusters, ClusterCorrNodeName, "PHObject");
0437     cemcNode->addNode(clusterNode);
0438   }
0439 }
0440 int RawClusterPositionCorrection::End(PHCompositeNode * /*topNode*/)
0441 {
0442   return Fun4AllReturnCodes::EVENT_OK;
0443 }