CalorimeterEmbedding/src/caloTowerEmbed.cc

0001 #include "caloTowerEmbed.h"
0002
0003 #include <calobase/TowerInfo.h>  // for TowerInfo
0004 #include <calobase/TowerInfoContainer.h>
0005 #include <calobase/TowerInfoContainerv1.h>
0006 #include <calobase/TowerInfoContainerv2.h>
0007 #include <calobase/TowerInfov1.h>
0008 #include <calobase/TowerInfoDefs.h>
0009 #include <calobase/RawTowerGeom.h>
0010 #include <calobase/RawTowerGeomContainer.h>
0011
0012 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0013
0014 #include <ffamodules/CDBInterface.h>
0015
0016 #include <ffaobjects/EventHeader.h>
0017
0018 #include <fun4all/Fun4AllServer.h>
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 #include <fun4all/SubsysReco.h>  // for SubsysReco
0021
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/PHIODataNode.h>    // for PHIODataNode
0024 #include <phool/PHNode.h>          // for PHNode
0025 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0026 #include <phool/PHObject.h>        // for PHObject
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>
0029 #include <phool/recoConsts.h>
0030
0031 #include <cstdlib>    // for exit
0032 #include <exception>  // for exception
0033 #include <iostream>   // for operator<<, basic_ostream
0034 #include <stdexcept>  // for runtime_error
0035
0036
0037 //____________________________________________________________________________..
0038 caloTowerEmbed::caloTowerEmbed(const std::string &name)
0039   : SubsysReco(name)
0040   , m_runNumber(-1)
0041   , m_eventNumber(-1)
0042 {
0043   std::cout << "caloTowerEmbed::caloTowerEmbed(const std::string &name) Calling ctor" << std::endl;
0044 }
0045
0046 //____________________________________________________________________________..
0047 caloTowerEmbed::~caloTowerEmbed()
0048 {
0049   std::cout << "caloTowerEmbed::~caloTowerEmbed() Calling dtor" << std::endl;
0050 }
0051
0052 //____________________________________________________________________________..
0053 int caloTowerEmbed::InitRun(PHCompositeNode *topNode)
0054 {
0055
0056   Fun4AllServer *se = Fun4AllServer::instance();
0057
0058   PHCompositeNode *simTopNode = se->topNode("TOPSim");
0059
0060   EventHeader *evtHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0061
0062   if (evtHeader)
0063   {
0064     m_runNumber = evtHeader->get_RunNumber();
0065   }
0066   else
0067   {
0068     m_runNumber = -1;
0069   }
0070   if(Verbosity()) std::cout << "at run" << m_runNumber << std::endl;
0071
0072   try
0073   {
0074     CreateNodeTree(topNode);
0075   }
0076   catch (std::exception &e)
0077   {
0078     std::cout << e.what() << std::endl;
0079     return Fun4AllReturnCodes::ABORTRUN;
0080   }
0081   topNode->print();
0082   simTopNode->print();
0083
0084   return Fun4AllReturnCodes::EVENT_OK;
0085 }
0086
0087 //____________________________________________________________________________..
0088 int caloTowerEmbed::process_event(PHCompositeNode* topNode)
0089 {
0090   ++m_eventNumber;
0091
0092   for(int caloIndex=0; caloIndex<3; caloIndex++){
0093
0094     if(Verbosity()) std::cout << "event " << m_eventNumber << " working on caloIndex " << caloIndex;
0095     if(caloIndex == 0) std::cout << " CEMC" << std::endl;
0096     else if(caloIndex == 1) std::cout << " HCALIN" << std::endl;
0097     else if(caloIndex == 2) std::cout << " HCALOUT" << std::endl;
0098     RawTowerDefs::keytype keyData = 0;
0099     RawTowerDefs::keytype keySim = 0;
0100
0101     unsigned int ntowers = _data_towers[caloIndex]->size();
0102     for (unsigned int channel = 0; channel < ntowers; channel++)
0103       {
0104     unsigned int data_key = _data_towers[caloIndex]->encode_key(channel);
0105     unsigned int sim_key = _sim_towers[caloIndex]->encode_key(channel);
0106
0107     int ieta_data = _data_towers[caloIndex]->getTowerEtaBin(data_key);
0108     int iphi_data = _data_towers[caloIndex]->getTowerPhiBin(data_key);
0109     int ieta_sim = _sim_towers[caloIndex]->getTowerEtaBin(sim_key);
0110     int iphi_sim = _sim_towers[caloIndex]->getTowerPhiBin(sim_key);
0111
0112     if(caloIndex == 0){
0113       if(!m_useRetower){
0114         keyData = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta_data, iphi_data);
0115         keySim = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta_sim, iphi_sim);
0116       }
0117       else{
0118         keyData = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta_data, iphi_data);
0119         keySim = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta_sim, iphi_sim);
0120       }
0121     }else if(caloIndex == 1){
0122       keyData = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta_data, iphi_data);
0123       keySim = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta_sim, iphi_sim);
0124     }else if(caloIndex == 2){
0125       keyData = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta_data, iphi_data);
0126       keySim = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta_sim, iphi_sim);
0127     }
0128
0129     TowerInfo *caloinfo_data = _data_towers[caloIndex]->get_tower_at_channel(channel);
0130     TowerInfo *caloinfo_sim = _sim_towers[caloIndex]->get_tower_at_channel(channel);
0131
0132     float data_E = caloinfo_data->get_energy();
0133     float sim_E = caloinfo_sim->get_energy();
0134     float embed_E = data_E + sim_E;
0135
0136     float data_phi = 0.0;
0137     float sim_phi = 0.0;
0138
0139     float data_eta = 0.0;
0140     float sim_eta = 0.0;
0141
0142     if(caloIndex == 0){
0143       if(!m_useRetower){
0144         data_phi = tower_geom->get_tower_geometry(keyData)->get_phi();
0145         data_eta = tower_geom->get_tower_geometry(keyData)->get_eta();
0146
0147         sim_phi = tower_geom->get_tower_geometry(keySim)->get_phi();
0148         sim_eta = tower_geom->get_tower_geometry(keySim)->get_eta();
0149       }else{
0150         data_phi = tower_geomIH->get_tower_geometry(keyData)->get_phi();
0151         data_eta = tower_geomIH->get_tower_geometry(keyData)->get_eta();
0152
0153         sim_phi = tower_geomIH->get_tower_geometry(keySim)->get_phi();
0154         sim_eta = tower_geomIH->get_tower_geometry(keySim)->get_eta();
0155       }
0156
0157       if(data_phi == sim_phi && data_eta == sim_eta){
0158         _data_towers[caloIndex]->get_tower_at_channel(channel)->set_energy(embed_E);
0159       }else{
0160         if(Verbosity()) std::cout << "eta and phi values in CEMC do not match between data and simulation, removing this event" << std::endl;
0161         return Fun4AllReturnCodes::ABORTEVENT;
0162       }
0163
0164     }else if(caloIndex == 1){
0165       data_phi = tower_geomIH->get_tower_geometry(keyData)->get_phi();
0166       data_eta = tower_geomIH->get_tower_geometry(keyData)->get_eta();
0167
0168       sim_phi = tower_geomIH->get_tower_geometry(keySim)->get_phi();
0169       sim_eta = tower_geomIH->get_tower_geometry(keySim)->get_eta();
0170
0171       if(data_phi == sim_phi && data_eta == sim_eta){
0172         _data_towers[caloIndex]->get_tower_at_channel(channel)->set_energy(embed_E);
0173       }else{
0174         if(Verbosity()) std::cout << "eta and phi values in HCALIN do not match between data and simulation, removing this event" << std::endl;
0175         return Fun4AllReturnCodes::ABORTEVENT;
0176       }
0177
0178     }else if(caloIndex == 2){
0179       data_phi = tower_geomOH->get_tower_geometry(keyData)->get_phi();
0180       data_eta = tower_geomOH->get_tower_geometry(keyData)->get_eta();
0181
0182       sim_phi = tower_geomOH->get_tower_geometry(keySim)->get_phi();
0183       sim_eta = tower_geomOH->get_tower_geometry(keySim)->get_eta();
0184
0185       if(data_phi == sim_phi && data_eta == sim_eta){
0186         _data_towers[caloIndex]->get_tower_at_channel(channel)->set_energy(embed_E);
0187       }else{
0188         if(Verbosity()) std::cout << "eta and phi values in HCALOUT do not match between data and simulation, removing this event" << std::endl;
0189         return Fun4AllReturnCodes::ABORTEVENT;
0190       }
0191
0192     }
0193       }//end loop over channels
0194   }//end loop over calorimeters
0195
0196   return Fun4AllReturnCodes::EVENT_OK;
0197 }
0198
0199 int caloTowerEmbed::End(PHCompositeNode *topNode)
0200 {
0201
0202   return Fun4AllReturnCodes::EVENT_OK;
0203
0204 }
0205
0206 void caloTowerEmbed::CreateNodeTree(PHCompositeNode *topNode)
0207 {
0208   std::cout << "creating node" << std::endl;
0209
0210   Fun4AllServer *se = Fun4AllServer::instance();
0211   PHCompositeNode *simTopNode = se->topNode("TOPSim");
0212
0213   tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0214   if(!tower_geom){
0215     std::cerr << Name() << "::" << __PRETTY_FUNCTION__
0216           << "tower geom CEMC missing, doing nothing." << std::endl;
0217     throw std::runtime_error(
0218                  "Failed to find TOWERGEOM_CEMC node");
0219   }
0220
0221   tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0222   if(!tower_geomIH){
0223     std::cerr << Name() << "::" << __PRETTY_FUNCTION__
0224           << "tower geom HCALIN missing, doing nothing." << std::endl;
0225     throw std::runtime_error(
0226                  "Failed to find TOWERGEOM_HCALIN node");
0227   }
0228
0229
0230   tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0231   if(!tower_geomOH){
0232     std::cerr << Name() << "::" << __PRETTY_FUNCTION__
0233           << "tower geom HCALOUT missing, doing nothing." << std::endl;
0234     throw std::runtime_error(
0235                  "Failed to find TOWERGEOM_HCALOUT node");
0236   }
0237
0238   PHNodeIterator dataIter(topNode);
0239   PHNodeIterator simIter(simTopNode);
0240
0241   //data top node first
0242
0243   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(dataIter.findFirst(
0244       "PHCompositeNode", "DST"));
0245   if (!dstNode)
0246   {
0247     std::cerr << Name() << "::" << __PRETTY_FUNCTION__
0248               << "DST Node missing, doing nothing." << std::endl;
0249     throw std::runtime_error(
0250         "Failed to find DST node in RawTowerCalibration::CreateNodes");
0251   }
0252
0253
0254   PHCompositeNode *dstNodeSim = dynamic_cast<PHCompositeNode *>(simIter.findFirst(
0255       "PHCompositeNode", "DST"));
0256   if (!dstNodeSim)
0257   {
0258     std::cerr << Name() << "::" << __PRETTY_FUNCTION__
0259               << "DSTsim Node missing, doing nothing." << std::endl;
0260     throw std::runtime_error(
0261         "Failed to find DSTsim node in RawTowerCalibration::CreateNodes");
0262   }
0263
0264   for(int i=0; i<3; i++){
0265
0266     std::string detector = "CEMC";
0267     if(i==0 && m_useRetower) detector = "CEMC_RETOWER";
0268     if(i==1) detector = "HCALIN";
0269     else if(i==2) detector = "HCALOUT";
0270
0271
0272     //data
0273     std::string DataTowerNodeName = "TOWERINFO_CALIB_" + detector;
0274     _data_towers[i] = findNode::getClass<TowerInfoContainer>(dstNode,
0275                                    DataTowerNodeName);
0276     if (!_data_towers[i])
0277       {
0278     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0279           << "Tower Calib Node missing, doing nothing." << std::endl;
0280     throw std::runtime_error(
0281                  "Failed to find Calib Tower node in caloTowerEmbed::CreateNodes");
0282       }
0283
0284
0285     //sim
0286     std::string SimTowerNodeName = "TOWERINFO_CALIB_" + detector;
0287     _sim_towers[i] = findNode::getClass<TowerInfoContainer>(dstNodeSim,
0288                                                            SimTowerNodeName);
0289     if (!_sim_towers[i])
0290       {
0291     std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0292           << "Tower Calib Sim Node missing, doing nothing." << std::endl;
0293     throw std::runtime_error(
0294                  "Failed to find Calib Tower Sim node in caloTowerEmbed::CreateNodes");
0295       }
0296
0297   }//end loop over calorimeter types
0298
0299   return;
0300 }