Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:35 

0001 #include "RawTowerCalibration.h"
0002 
0003 #include <calobase/RawTower.h>
0004 #include <calobase/RawTowerContainer.h>
0005 #include <calobase/RawTowerDefs.h>
0006 #include <calobase/RawTowerGeom.h>
0007 #include <calobase/RawTowerGeomContainer.h>
0008 #include <calobase/RawTowerv2.h>
0009 #include <calobase/TowerInfo.h>
0010 #include <calobase/TowerInfoContainer.h>
0011 #include <calobase/TowerInfoContainerv1.h>
0012 #include <calobase/TowerInfoContainerv2.h>
0013 #include <calobase/TowerInfov1.h>
0014 #include <calobase/TowerInfov2.h>
0015 
0016 #include <phparameter/PHParameters.h>
0017 
0018 #include <cdbobjects/CDBTTree.h>
0019 
0020 #include <ffamodules/CDBInterface.h>
0021 
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/SubsysReco.h>
0024 
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/PHIODataNode.h>
0027 #include <phool/PHNode.h>
0028 #include <phool/PHNodeIterator.h>
0029 #include <phool/PHObject.h>
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>
0032 
0033 #include <TSystem.h>
0034 
0035 #include <cassert>
0036 #include <cstdlib>
0037 #include <exception>
0038 #include <fstream>
0039 #include <iostream>
0040 #include <map>
0041 #include <stdexcept>
0042 #include <string>
0043 #include <utility>
0044 
0045 RawTowerCalibration::RawTowerCalibration(const std::string &name)
0046   : SubsysReco(name)
0047   , _calib_algorithm(kNo_calibration)
0048   , m_Detector("NONE")
0049   , _calib_tower_node_prefix("CALIB")
0050   , _raw_tower_node_prefix("RAW")
0051   , _tower_calib_params(name)
0052 {
0053   //_tower_type = -1;
0054 }
0055 
0056 RawTowerCalibration::~RawTowerCalibration()
0057 {
0058   delete m_CDBTTree;
0059 }
0060 
0061 int RawTowerCalibration::InitRun(PHCompositeNode *topNode)
0062 {
0063   PHNodeIterator iter(topNode);
0064 
0065   // Looking for the DST node
0066   PHCompositeNode *dstNode;
0067   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode",
0068                                                            "DST"));
0069   if (!dstNode)
0070   {
0071     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0072               << "DST Node missing, doing nothing." << std::endl;
0073     exit(1);
0074   }
0075 
0076   try
0077   {
0078     CreateNodes(topNode);
0079   }
0080   catch (std::exception &e)
0081   {
0082     std::cout << e.what() << std::endl;
0083     return Fun4AllReturnCodes::ABORTRUN;
0084   }
0085 
0086   if (_calib_algorithm == kDbfile_tbt_gain_corr)
0087   {
0088     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0089               << "kDbfile_tbt_gain_corr  chosen but not implemented"
0090               << std::endl;
0091     return Fun4AllReturnCodes::ABORTRUN;
0092   }
0093 
0094   return Fun4AllReturnCodes::EVENT_OK;
0095 }
0096 
0097 int RawTowerCalibration::process_event(PHCompositeNode * /*topNode*/)
0098 {
0099   if (Verbosity())
0100   {
0101     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0102               << "Process event entered" << std::endl;
0103   }
0104 
0105   if (m_UseTowerInfo != 1)
0106   {
0107     RawTowerContainer::ConstRange begin_end = _raw_towers->getTowers();
0108     RawTowerContainer::ConstIterator rtiter;
0109 
0110     for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0111     {
0112       const RawTowerDefs::keytype key = rtiter->first;
0113 
0114       const RawTower *raw_tower = rtiter->second;
0115       assert(raw_tower);
0116 
0117       RawTowerGeom *raw_tower_geom = rawtowergeom->get_tower_geometry(
0118           raw_tower->get_id());
0119       assert(raw_tower_geom);
0120 
0121       if (_tower_type >= 0)
0122       {
0123         // Skip towers that don't match the type we are supposed to calibrate
0124         if (_tower_type != raw_tower_geom->get_tower_type())
0125         {
0126           continue;
0127         }
0128       }
0129 
0130       if (_calib_algorithm == kNo_calibration)
0131       {
0132         _calib_towers->AddTower(key, new RawTowerv2(*raw_tower));
0133       }
0134       else if (_calib_algorithm == kSimple_linear_calibration)
0135       {
0136         const double raw_energy = raw_tower->get_energy();
0137         const double calib_energy = (raw_energy - _pedstal_ADC) * _calib_const_GeV_ADC;
0138 
0139         RawTower *calib_tower = new RawTowerv2(*raw_tower);
0140         calib_tower->set_energy(calib_energy);
0141         _calib_towers->AddTower(key, calib_tower);
0142       }
0143       else if (_calib_algorithm == kTower_by_tower_calibration)
0144       {
0145         RawTowerDefs::CalorimeterId caloid = RawTowerDefs::decode_caloid(key);
0146         const int eta = raw_tower->get_bineta();
0147         const int phi = raw_tower->get_binphi();
0148 
0149         double tower_by_tower_calib = 1.;
0150         if (caloid == RawTowerDefs::LFHCAL)
0151         {
0152           const int l = raw_tower->get_binl();
0153           const std::string calib_const_name("calib_const_eta" + std::to_string(eta) + "_phi" + std::to_string(phi) + "_l" + std::to_string(l));
0154 
0155           tower_by_tower_calib = _tower_calib_params.get_double_param(calib_const_name);
0156 
0157           if (_pedestal_file == true)
0158           {
0159             const std::string pedstal_name("PedCentral_ADC_eta" + std::to_string(eta) + "_phi" + std::to_string(phi) + "_l" + std::to_string(l));
0160             _pedstal_ADC =
0161                 _tower_calib_params.get_double_param(pedstal_name);
0162           }
0163 
0164           if (_GeV_ADC_file == true)
0165           {
0166             const std::string GeVperADCname("GeVperADC_eta" + std::to_string(eta) + "_phi" + std::to_string(phi) + "_l" + std::to_string(l));
0167             _calib_const_GeV_ADC =
0168                 _tower_calib_params.get_double_param(GeVperADCname);
0169           }
0170         }
0171         else
0172         {
0173           const std::string calib_const_name("calib_const_eta" + std::to_string(eta) + "_phi" + std::to_string(phi));
0174 
0175           tower_by_tower_calib = _tower_calib_params.get_double_param(calib_const_name);
0176 
0177           if (_pedestal_file == true)
0178           {
0179             const std::string pedstal_name("PedCentral_ADC_eta" + std::to_string(eta) + "_phi" + std::to_string(phi));
0180             _pedstal_ADC =
0181                 _tower_calib_params.get_double_param(pedstal_name);
0182           }
0183 
0184           if (_GeV_ADC_file == true)
0185           {
0186             const std::string GeVperADCname("GeVperADC_eta" + std::to_string(eta) + "_phi" + std::to_string(phi));
0187             _calib_const_GeV_ADC =
0188                 _tower_calib_params.get_double_param(GeVperADCname);
0189           }
0190         }
0191         const double raw_energy = raw_tower->get_energy();
0192         const double calib_energy = (raw_energy - _pedstal_ADC) * _calib_const_GeV_ADC * tower_by_tower_calib;
0193 
0194         RawTower *calib_tower = new RawTowerv2(*raw_tower);
0195         calib_tower->set_energy(calib_energy);
0196         _calib_towers->AddTower(key, calib_tower);
0197       }
0198       // else if  // eventally this will be done exclusively of tow_by_tow
0199       else if (_calib_algorithm == kDbfile_tbt_gain_corr)
0200       {
0201         if (m_Detector.c_str()[0] == 'H')
0202         {
0203           std::string url = CDBInterface::instance()->getUrl("HCALTBYTCORR");
0204           if (url.empty())
0205           {
0206             std::cout << PHWHERE << " Could not get Hcal Calibration for domain HCALTBYTCORR" << std::endl;
0207             gSystem->Exit(1);
0208             exit(1);
0209           }
0210 
0211           m_CDBTTree = new CDBTTree(url);
0212         }
0213         else if (m_Detector.c_str()[0] == 'C')
0214         {
0215           std::string url = CDBInterface::instance()->getUrl("CEMCTBYTCORR");
0216           if (url.empty())
0217           {
0218             std::cout << PHWHERE << " Could not get Cemc Calibration for domain CEMCTBYTCORR" << std::endl;
0219             gSystem->Exit(1);
0220             exit(1);
0221           }
0222 
0223           m_CDBTTree = new CDBTTree(url);
0224         }
0225         if (!m_CDBTTree)
0226         {
0227           std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0228                     << "kDbfile_tbt_gain_corr  chosen but no file loaded" << std::endl;
0229           return Fun4AllReturnCodes::ABORTRUN;
0230         }
0231 
0232         float gain_factor = -888;
0233         //      gain_factor = _cal_dbfile->getCorr(key);
0234 
0235         const int eta = raw_tower->get_bineta();
0236         const int phi = raw_tower->get_binphi();
0237         unsigned int etaphikey = phi;
0238         etaphikey = (etaphikey << 16U) + eta;
0239 
0240         gain_factor = m_CDBTTree->GetFloatValue(etaphikey, "etaphi");
0241 
0242         const double raw_energy = raw_tower->get_energy();
0243         RawTower *calib_tower = new RawTowerv2(*raw_tower);
0244 
0245         // still include separate _calib_const_GeV_ADC factor
0246         // for global shifts.
0247 
0248         float corr_energy = raw_energy * gain_factor * _calib_const_GeV_ADC;
0249         calib_tower->set_energy(corr_energy);
0250         _calib_towers->AddTower(key, calib_tower);
0251       }
0252       else
0253       {
0254         std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0255                   << " invalid calibration algorithm #" << _calib_algorithm
0256                   << std::endl;
0257 
0258         return Fun4AllReturnCodes::ABORTRUN;
0259       }
0260     }
0261   }
0262 
0263   if (m_UseTowerInfo > 0)
0264   {
0265     unsigned int ntowers = _raw_towerinfos->size();
0266     for (unsigned int channel = 0; channel < ntowers; channel++)
0267     {
0268       // unsigned int key =rtiter->first;
0269       unsigned int key = _raw_towerinfos->encode_key(channel);
0270 
0271       TowerInfo *raw_tower = _raw_towerinfos->get_tower_at_channel(channel);
0272       assert(raw_tower);
0273 
0274       TowerInfo *calib_tower = new TowerInfov1();
0275 
0276       if (_calib_algorithm == kNo_calibration)
0277       {
0278         calib_tower->set_energy(0);
0279       }
0280 
0281       else if (_calib_algorithm == kSimple_linear_calibration)
0282       {
0283         const double raw_energy = raw_tower->get_energy();
0284         const double calib_energy = (raw_energy - _pedstal_ADC) * _calib_const_GeV_ADC;
0285         calib_tower->set_energy(calib_energy);
0286       }
0287       else if (_calib_algorithm == kTower_by_tower_calibration)
0288       {
0289         const int eta = _raw_towerinfos->getTowerEtaBin(key);
0290         const int phi = _raw_towerinfos->getTowerPhiBin(key);
0291         double tower_by_tower_calib = 1.;
0292         const std::string calib_const_name("calib_const_eta" + std::to_string(eta) + "_phi" + std::to_string(phi));
0293         tower_by_tower_calib = _tower_calib_params.get_double_param(calib_const_name);
0294         if (_pedestal_file == true)
0295         {
0296           const std::string pedstal_name("PedCentral_ADC_eta" + std::to_string(eta) + "_phi" + std::to_string(phi));
0297           _pedstal_ADC =
0298               _tower_calib_params.get_double_param(pedstal_name);
0299         }
0300         if (_GeV_ADC_file == true)
0301         {
0302           const std::string GeVperADCname("GeVperADC_eta" + std::to_string(eta) + "_phi" + std::to_string(phi));
0303           _calib_const_GeV_ADC =
0304               _tower_calib_params.get_double_param(GeVperADCname);
0305         }
0306         const double raw_energy = raw_tower->get_energy();
0307         const double calib_energy = (raw_energy - _pedstal_ADC) * _calib_const_GeV_ADC * tower_by_tower_calib;
0308         calib_tower->set_energy(calib_energy);
0309       }
0310       else if (_calib_algorithm == kDbfile_tbt_gain_corr)
0311       {
0312         if (!m_CDBTTree)
0313         {
0314           std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0315                     << "kDbfile_tbt_gain_corr  chosen but no file loaded" << std::endl;
0316           return Fun4AllReturnCodes::ABORTRUN;
0317         }
0318 
0319         float gain_factor = -888;
0320         const int eta = _raw_towerinfos->getTowerEtaBin(key);
0321         const int phi = _raw_towerinfos->getTowerPhiBin(key);
0322         unsigned int etaphikey = phi;
0323         etaphikey = (etaphikey << 16U) + eta;
0324 
0325         gain_factor = m_CDBTTree->GetFloatValue(etaphikey, "etaphi");
0326         const double raw_energy = raw_tower->get_energy();
0327         float corr_energy = raw_energy * gain_factor * _calib_const_GeV_ADC;
0328         calib_tower->set_energy(corr_energy);
0329       }
0330       else
0331       {
0332         std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0333                   << " invalid calibration algorithm #" << _calib_algorithm
0334                   << std::endl;
0335 
0336         return Fun4AllReturnCodes::ABORTRUN;
0337       }
0338 
0339       TowerInfo *calibrated_towerinfo = _calib_towerinfos->get_tower_at_channel(channel);
0340       calibrated_towerinfo->set_energy(calib_tower->get_energy());
0341       delete calib_tower;
0342     }
0343   }
0344 
0345   //  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0346 
0347   /*
0348   int towcount =0;
0349   RawTowerContainer::ConstRange begin_end2 = _calib_towers->getTowers();
0350   RawTowerContainer::ConstIterator rtiter2;
0351   std::cout << "Etowers  ---------------------================---------"
0352             << std::endl;
0353 
0354   for (rtiter2 = begin_end2.first; rtiter2 != begin_end2.second; ++rtiter2)
0355   {
0356     const RawTowerDefs::keytype key = rtiter2->first;
0357     const RawTower *aftcal_tower = rtiter2->second;
0358 
0359     if (towcount++ < 10)
0360       {
0361         std::cout << "E tow: " << key << "   " << aftcal_tower->get_energy()
0362                 << std::endl;
0363       }
0364 
0365 
0366   }
0367 
0368   */
0369 
0370   if (Verbosity())
0371   {
0372     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0373               << "input sum energy = " << _raw_towers->getTotalEdep()
0374               << ", output sum digitalized value = "
0375               << _calib_towers->getTotalEdep() << std::endl;
0376   }
0377   return Fun4AllReturnCodes::EVENT_OK;
0378 }
0379 
0380 int RawTowerCalibration::End(PHCompositeNode * /*topNode*/)
0381 {
0382   return Fun4AllReturnCodes::EVENT_OK;
0383 }
0384 
0385 void RawTowerCalibration::CreateNodes(PHCompositeNode *topNode)
0386 {
0387   PHNodeIterator iter(topNode);
0388   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
0389       "PHCompositeNode", "RUN"));
0390   if (!runNode)
0391   {
0392     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0393               << "Run Node missing, doing nothing." << std::endl;
0394     throw std::runtime_error(
0395         "Failed to find Run node in RawTowerCalibration::CreateNodes");
0396   }
0397 
0398   TowerGeomNodeName = "TOWERGEOM_" + m_Detector;
0399   rawtowergeom = findNode::getClass<RawTowerGeomContainer>(topNode,
0400                                                            TowerGeomNodeName);
0401   if (!rawtowergeom)
0402   {
0403     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0404               << " " << TowerGeomNodeName << " Node missing, doing bail out!"
0405               << std::endl;
0406     throw std::runtime_error(
0407         "Failed to find " + TowerGeomNodeName + " node in RawTowerCalibration::CreateNodes");
0408   }
0409 
0410   if (Verbosity() >= 1)
0411   {
0412     rawtowergeom->identify();
0413   }
0414 
0415   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
0416       "PHCompositeNode", "DST"));
0417   if (!dstNode)
0418   {
0419     std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0420               << "DST Node missing, doing nothing." << std::endl;
0421     throw std::runtime_error(
0422         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0423   }
0424 
0425   if (m_UseTowerInfo != 1)
0426   {
0427     RawTowerNodeName = "TOWER_" + _raw_tower_node_prefix + "_" + m_Detector;
0428     _raw_towers = findNode::getClass<RawTowerContainer>(dstNode,
0429                                                         RawTowerNodeName);
0430     if (!_raw_towers)
0431     {
0432       std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0433                 << " " << RawTowerNodeName << " Node missing, doing bail out!"
0434                 << std::endl;
0435       throw std::runtime_error(
0436           "Failed to find " + RawTowerNodeName + " node in RawTowerCalibration::CreateNodes");
0437     }
0438   }
0439   if (m_UseTowerInfo > 0)
0440   {
0441     RawTowerInfoNodeName = "TOWERINFO_" + _raw_tower_node_prefix + "_" + m_Detector;
0442     _raw_towerinfos = findNode::getClass<TowerInfoContainerv1>(dstNode, RawTowerInfoNodeName);
0443 
0444     if (!_raw_towerinfos)
0445     {
0446       std::cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
0447                 << " " << RawTowerInfoNodeName << " Node missing, doing bail out!"
0448                 << std::endl;
0449       throw std::runtime_error(
0450           "Failed to find " + RawTowerInfoNodeName + " node in RawTowerCalibration::CreateNodes");
0451     }
0452   }
0453 
0454   // Create the tower nodes on the tree
0455   PHNodeIterator dstiter(dstNode);
0456   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
0457       "PHCompositeNode", m_Detector));
0458   if (!DetNode)
0459   {
0460     DetNode = new PHCompositeNode(m_Detector);
0461     dstNode->addNode(DetNode);
0462   }
0463 
0464   // Be careful as a previous calibrator may have been registered for this detector
0465   if (m_UseTowerInfo != 1)
0466   {
0467     CaliTowerNodeName = "TOWER_" + _calib_tower_node_prefix + "_" + m_Detector;
0468     _calib_towers = findNode::getClass<RawTowerContainer>(DetNode,
0469                                                           CaliTowerNodeName);
0470     if (!_calib_towers)
0471     {
0472       _calib_towers = new RawTowerContainer(_raw_towers->getCalorimeterID());
0473       PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(_calib_towers, CaliTowerNodeName, "PHObject");
0474       DetNode->addNode(towerNode);
0475     }
0476   }
0477   if (m_UseTowerInfo > 0)
0478   {
0479     CaliTowerInfoNodeName = "TOWERINFO_" + _calib_tower_node_prefix + "_" + m_Detector;
0480     if (!m_UseTowerInfoV2)
0481     {
0482       _calib_towerinfos = findNode::getClass<TowerInfoContainerv1>(DetNode, CaliTowerInfoNodeName);
0483     }
0484     else
0485     {
0486       _calib_towerinfos = findNode::getClass<TowerInfoContainerv2>(DetNode, CaliTowerInfoNodeName);
0487     }
0488     if (!_calib_towerinfos)
0489     {
0490       TowerInfoContainerv1::DETECTOR detec;
0491       if (m_Detector == "CEMC")
0492       {
0493         detec = TowerInfoContainer::DETECTOR::EMCAL;
0494       }
0495       else if (m_Detector == "HCALIN" || m_Detector == "HCALOUT")
0496       {
0497         detec = TowerInfoContainer::DETECTOR::HCAL;
0498       }
0499       else
0500       {
0501         std::cout << PHWHERE << "Detector not implemented into the TowerInfoContainer object, defaulting to HCal implementation." << std::endl;
0502         detec = TowerInfoContainer::DETECTOR::HCAL;
0503       }
0504       if (!m_UseTowerInfoV2)
0505       {
0506         _calib_towerinfos = new TowerInfoContainerv1(detec);
0507       }
0508       else
0509       {
0510         _calib_towerinfos = new TowerInfoContainerv2(detec);
0511       }
0512 
0513       PHIODataNode<PHObject> *towerinfoNode = new PHIODataNode<PHObject>(_calib_towerinfos, CaliTowerInfoNodeName, "PHObject");
0514       DetNode->addNode(towerinfoNode);
0515     }
0516   }
0517 
0518   return;
0519 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!