packages/jetbackground/RetowerCEMC.cc

0001 #include "RetowerCEMC.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/RawTowerv1.h>
0009
0010 #include <calobase/TowerInfo.h>
0011 #include <calobase/TowerInfoContainer.h>
0012 #include <calobase/TowerInfoDefs.h>
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/SubsysReco.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>
0019 #include <phool/PHNode.h>
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/PHObject.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>
0024
0025 // standard includes
0026 #include <cstdlib>
0027 #include <iostream>
0028 #include <utility>
0029
0030 RetowerCEMC::RetowerCEMC(const std::string &name)
0031   : SubsysReco(name)
0032 {
0033 }
0034
0035 int RetowerCEMC::InitRun(PHCompositeNode *topNode)
0036 {
0037   CreateNode(topNode);
0038   get_first_phi_index(topNode);
0039   if (_weighted_energy_distribution == 1)
0040   {
0041     get_weighted_fraction(topNode);
0042   }
0043   else
0044   {
0045     get_fraction(topNode);
0046   }
0047   return Fun4AllReturnCodes::EVENT_OK;
0048 }
0049
0050 int RetowerCEMC::process_event(PHCompositeNode *topNode)
0051 {
0052   if (Verbosity() > 0)
0053   {
0054     std::cout << "RetowerCEMC::process_event: entering" << std::endl;
0055   }
0056
0057   RawTowerContainer *towersEM3 = nullptr;
0058   TowerInfoContainer *towerinfosEM3 = nullptr;
0059   if (m_use_towerinfo)
0060   {
0061     EMTowerName = m_towerNodePrefix + "_CEMC";
0062     towerinfosEM3 = findNode::getClass<TowerInfoContainer>(topNode, EMTowerName);
0063     if (!towerinfosEM3)
0064     {
0065       std::cout << "RetowerCEMC::process_event: Cannot find node " << EMTowerName << std::endl;
0066       exit(1);
0067     }
0068   }
0069   else
0070   {
0071     towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0072     if (Verbosity() > 0)
0073     {
0074       std::cout << "RetowerCEMC::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC towers" << std::endl;
0075     }
0076   }
0077
0078   if (m_use_towerinfo)
0079   {
0080     unsigned int nchannels = towerinfosEM3->size();
0081     for (unsigned int channel = 0; channel < nchannels; channel++)
0082     {
0083       TowerInfo *tower = towerinfosEM3->get_tower_at_channel(channel);
0084       unsigned int channelkey = towerinfosEM3->encode_key(channel);
0085       int ieta = towerinfosEM3->getTowerEtaBin(channelkey);
0086       int iphi = towerinfosEM3->getTowerPhiBin(channelkey);
0087       rawtower_e[ieta][iphi] = tower->get_energy();
0088       rawtower_time[ieta][iphi] = tower->get_time_float();
0089       rawtower_status[ieta][iphi] = tower->get_isHot() || tower->get_isNoCalib() || tower->get_isNotInstr() || tower->get_isBadChi2();
0090     }
0091     EMRetowerName = m_towerNodePrefix + "_CEMC_RETOWER";
0092     TowerInfoContainer *emcal_retower = findNode::getClass<TowerInfoContainer>(topNode, EMRetowerName);
0093     if (Verbosity() > 0)
0094     {
0095       std::cout << "RetowerCEMC::process_event: filling " << EMRetowerName << " node" << std::endl;
0096     }
0097     for (int ieta_ihcal = 0; ieta_ihcal < neta_ihcal; ++ieta_ihcal)
0098     {
0099       for (int iphi_ihcal = 0; iphi_ihcal < nphi_ihcal; ++iphi_ihcal)
0100       {
0101         double retower_e_temp = 0;
0102         double retower_time_temp = 0;
0103         double retower_badarea = 0;
0104         for (int ieta_emcal = retower_lowerbound_originaltower_ieta[ieta_ihcal]; ieta_emcal <= retower_upperbound_originaltower_ieta[ieta_ihcal]; ++ieta_emcal)
0105         {
0106           for (int iphi_emcal = retower_first_lowerbound_originaltower_iphi + (iphi_ihcal * 4); iphi_emcal < retower_first_lowerbound_originaltower_iphi + iphi_ihcal * 4 + 4; ++iphi_emcal)
0107           {
0108             int iphi_emcal_wrap = iphi_emcal;
0109             if (iphi_emcal > nphi_emcal - 1)
0110             {
0111               iphi_emcal_wrap -= nphi_emcal;
0112             }
0113             double fraction_temp;
0114             if (ieta_emcal == retower_lowerbound_originaltower_ieta[ieta_ihcal])
0115             {
0116               fraction_temp = retower_lowerbound_originaltower_fraction[ieta_ihcal];
0117             }
0118             else if (ieta_emcal == retower_upperbound_originaltower_ieta[ieta_ihcal])
0119             {
0120               fraction_temp = retower_upperbound_originaltower_fraction[ieta_ihcal];
0121             }
0122             else
0123             {
0124               fraction_temp = 1;
0125             }
0126             if (rawtower_status[ieta_emcal][iphi_emcal_wrap])
0127             {
0128               retower_badarea += fraction_temp;
0129             }
0130             else
0131             {
0132               retower_e_temp += rawtower_e[ieta_emcal][iphi_emcal_wrap] * fraction_temp;
0133               retower_time_temp += rawtower_time[ieta_emcal][iphi_emcal_wrap] * rawtower_e[ieta_emcal][iphi_emcal_wrap] * fraction_temp;
0134             }
0135           }
0136         }
0137         unsigned int towerkey = TowerInfoDefs::encode_hcal(ieta_ihcal, iphi_ihcal);
0138         unsigned int towerindex = emcal_retower->decode_key(towerkey);
0139         TowerInfo *towerinfo = emcal_retower->get_tower_at_channel(towerindex);
0140         double scalefactor = retower_badarea / retower_totalarea[ieta_ihcal];
0141         if (scalefactor > _frac_cut)
0142         {
0143           towerinfo->set_energy(0);
0144           towerinfo->set_isHot(true);
0145         }
0146         else
0147         {
0148           towerinfo->set_energy(retower_e_temp / (double) (1 - scalefactor));
0149           if (retower_e_temp == 0)
0150           {
0151             towerinfo->set_time_float(0);
0152           }
0153           else
0154           {
0155             towerinfo->set_time_float((retower_time_temp / retower_e_temp));
0156           }
0157           towerinfo->set_chi2(scalefactor);
0158         }
0159       }
0160     }
0161   }
0162   else
0163   {
0164     for (int ieta_ihcal = 0; ieta_ihcal < neta_ihcal; ++ieta_ihcal)
0165     {
0166       for (int iphi_ihcal = 0; iphi_ihcal < nphi_ihcal; ++iphi_ihcal)
0167       {
0168         double retower_e_temp = 0;
0169         for (int ieta_emcal = retower_lowerbound_originaltower_ieta[ieta_ihcal]; ieta_emcal <= retower_upperbound_originaltower_ieta[ieta_ihcal]; ++ieta_emcal)
0170         {
0171           for (int iphi_emcal = retower_first_lowerbound_originaltower_iphi; iphi_emcal < retower_first_lowerbound_originaltower_iphi + iphi_ihcal * 4; ++iphi_emcal)
0172           {
0173             int iphi_emcal_wrap = iphi_emcal;
0174             if (iphi_emcal > nphi_emcal - 1)
0175             {
0176               iphi_emcal_wrap -= nphi_emcal;
0177             }
0178             RawTower *tower = towersEM3->getTower(ieta_emcal, iphi_emcal_wrap);
0179             double energy = tower->get_energy();
0180             double fraction_temp;
0181             if (ieta_emcal == retower_lowerbound_originaltower_ieta[ieta_ihcal])
0182             {
0183               fraction_temp = retower_lowerbound_originaltower_fraction[ieta_ihcal];
0184             }
0185             else if (ieta_emcal == retower_upperbound_originaltower_ieta[ieta_ihcal])
0186             {
0187               fraction_temp = retower_upperbound_originaltower_fraction[ieta_ihcal];
0188             }
0189             else
0190             {
0191               fraction_temp = 1;
0192             }
0193             retower_e_temp += energy * fraction_temp;
0194           }
0195         }
0196         RawTowerContainer *emcal_retower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0197         if (Verbosity() > 0)
0198         {
0199           std::cout << "RetowerCEMC::process_event: filling TOWER_CALIB_CEMC_RETOWER node, with initial size = " << emcal_retower->size() << std::endl;
0200         }
0201         RawTower *new_tower = new RawTowerv1();
0202         new_tower->set_energy(retower_e_temp);
0203         emcal_retower->AddTower(ieta_ihcal, iphi_ihcal, new_tower);
0204       }
0205     }
0206   }
0207   if (Verbosity() > 0)
0208   {
0209     std::cout << "RetowerCEMC::process_event: exiting" << std::endl;
0210   }
0211   return Fun4AllReturnCodes::EVENT_OK;
0212 }
0213
0214 int RetowerCEMC::CreateNode(PHCompositeNode *topNode)
0215 {
0216   PHNodeIterator iter(topNode);
0217   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0218   if (!dstNode)
0219   {
0220     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0221     return Fun4AllReturnCodes::ABORTRUN;
0222   }
0223   PHCompositeNode *emcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
0224   if (!emcalNode)
0225   {
0226     std::cout << PHWHERE << "EMCal Node note found, doing nothing." << std::endl;
0227   }
0228
0229   if (m_use_towerinfo)
0230   {
0231     EMRetowerName = m_towerNodePrefix + "_CEMC_RETOWER";
0232     IHTowerName = m_towerNodePrefix + "_HCALIN";
0233     TowerInfoContainer *test_emcal_retower = findNode::getClass<TowerInfoContainer>(topNode, EMRetowerName);
0234     TowerInfoContainer *hcal_towers = findNode::getClass<TowerInfoContainer>(topNode, IHTowerName);
0235     if (!test_emcal_retower)
0236     {
0237       if (Verbosity() > 0)
0238       {
0239         std::cout << "RetowerCEMC::CreateNode : creating " << EMRetowerName << " node " << std::endl;
0240       }
0241       if (!hcal_towers)
0242       {
0243         std::cout << PHWHERE << " Could not find input HCAL tower node: " << IHTowerName << std::endl;
0244         exit(1);
0245       }
0246       TowerInfoContainer *emcal_retower = dynamic_cast<TowerInfoContainer *>(hcal_towers->CloneMe());
0247       PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_retower, EMRetowerName, "PHObject");
0248       emcalNode->addNode(emcalTowerNode);
0249     }
0250     else
0251     {
0252       if (Verbosity() > 0)
0253       {
0254         std::cout << "RetowerCEMC::CreateNode : " << EMRetowerName << " already exists! " << std::endl;
0255       }
0256     }
0257   }
0258   else
0259   {
0260     RawTowerContainer *test_emcal_retower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0261     if (!test_emcal_retower)
0262     {
0263       if (Verbosity() > 0)
0264       {
0265         std::cout << "RetowerCEMC::CreateNode : creating TOWER_CALIB_CEMC_RETOWER node " << std::endl;
0266       }
0267       RawTowerContainer *emcal_retower = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
0268       PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_retower, "TOWER_CALIB_CEMC_RETOWER", "PHObject");
0269       emcalNode->addNode(emcalTowerNode);
0270     }
0271     else
0272     {
0273       if (Verbosity() > 0)
0274       {
0275         std::cout << "RetowerCEMC::CreateNode : TOWER_CALIB_CEMC_RETOWER already exists! " << std::endl;
0276       }
0277     }
0278   }
0279   return Fun4AllReturnCodes::EVENT_OK;
0280 }
0281
0282 void RetowerCEMC::get_first_phi_index(PHCompositeNode *topNode)
0283 {
0284   RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0285   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0286   if (!geomEM || !geomIH)
0287   {
0288     std::cout << "RetowerCEMC::get_first_phi_index: Could not locate geometry nodes" << std::endl;
0289     exit(1);
0290   }
0291
0292   bool find_first_lowerbound = false;
0293   int iphi_emcal = 0;
0294   while (iphi_emcal < nphi_emcal)
0295   {
0296     const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, 0, iphi_emcal);
0297     RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0298     int this_IHphibin = geomIH->get_phibin(tower_geom->get_phi());
0299     if (this_IHphibin == 0)
0300     {
0301       find_first_lowerbound = true;
0302       break;
0303     }
0304     iphi_emcal++;
0305   }
0306
0307   if (find_first_lowerbound && iphi_emcal == 0)
0308   {
0309     bool outofrange = false;
0310     int iphi_emcal_temp = nphi_emcal - 1;
0311     while (iphi_emcal_temp > iphi_emcal)
0312     {
0313       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, 0, iphi_emcal_temp);
0314       RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0315       int this_IHphibin = geomIH->get_phibin(tower_geom->get_phi());
0316       if (this_IHphibin == nphi_ihcal - 1)
0317       {
0318         outofrange = true;
0319         break;
0320       }
0321       iphi_emcal_temp--;
0322     }
0323     if (!outofrange)
0324     {
0325       std::cout << "RetowerCEMC::get_first_phi_index: could not find matching EMCal towers for iphi_ihcal = " << nphi_ihcal - 1 << std::endl;
0326       exit(1);
0327     }
0328     if (iphi_emcal_temp + 1 == nphi_emcal)
0329     {
0330       retower_first_lowerbound_originaltower_iphi = 0;
0331     }
0332     else
0333     {
0334       retower_first_lowerbound_originaltower_iphi = iphi_emcal_temp + 1;
0335     }
0336   }
0337   else if (!find_first_lowerbound)
0338   {
0339     std::cout << "RetowerCEMC::get_first_phi_index: could not find matching EMCal towers for iphi_ihcal = 0" << std::endl;
0340     exit(1);
0341   }
0342   else
0343   {
0344     retower_first_lowerbound_originaltower_iphi = iphi_emcal;
0345   }
0346 }
0347
0348 void RetowerCEMC::get_fraction(PHCompositeNode *topNode)
0349 {
0350   RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0351   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0352   if (!geomEM || !geomIH)
0353   {
0354     std::cout << "RetowerCEMC::get_fraction: Could not locate geometry nodes" << std::endl;
0355     exit(1);
0356   }
0357
0358   int ieta_emcal = 0;
0359   int ieta_ihcal = 0;
0360   while (ieta_emcal < neta_emcal && ieta_ihcal < neta_ihcal)
0361   {
0362     const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta_emcal, 0);
0363     RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0364     int this_IHetabin = geomIH->get_etabin(tower_geom->get_eta());
0365     if (this_IHetabin == ieta_ihcal)
0366     {
0367       retower_lowerbound_originaltower_ieta[ieta_ihcal] = ieta_emcal;
0368       ieta_ihcal++;
0369     }
0370     ieta_emcal++;
0371   }
0372   for (int ieta = 0; ieta < neta_ihcal; ++ieta)
0373   {
0374     if (ieta < neta_ihcal - 1)
0375     {
0376       retower_upperbound_originaltower_ieta[ieta] = retower_lowerbound_originaltower_ieta[ieta + 1] - 1;
0377     }
0378     else
0379     {
0380       retower_upperbound_originaltower_ieta[ieta] = neta_emcal - 1;
0381     }
0382     retower_lowerbound_originaltower_fraction[ieta] = 1;
0383     retower_upperbound_originaltower_fraction[ieta] = 1;
0384   }
0385 }
0386
0387 void RetowerCEMC::get_weighted_fraction(PHCompositeNode *topNode)
0388 {
0389   RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0390   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0391   if (!geomEM || !geomIH)
0392   {
0393     std::cout << "RetowerCEMC::get_weighted_fraction: Could not locate geometry nodes" << std::endl;
0394     exit(1);
0395   }
0396
0397   int ieta_emcal = 0;
0398   for (int ieta_ihcal = 0; ieta_ihcal < neta_ihcal; ++ieta_ihcal)
0399   {
0400     std::pair<double, double> range_ihcal = geomIH->get_etabounds(ieta_ihcal);
0401     double ihcal_lowerbound = range_ihcal.first;
0402     double ihcal_upperbound = range_ihcal.second;
0403
0404     bool found_lowerbound = false;
0405     bool found_upperbound = false;
0406     while ((!found_lowerbound || !found_upperbound) && ieta_emcal < neta_emcal)
0407     {
0408       std::pair<double, double> range_emcal = geomEM->get_etabounds(ieta_emcal);
0409       double emcal_lowerbound = range_emcal.first;
0410       double emcal_upperbound = range_emcal.second;
0411       if (!found_lowerbound)
0412       {
0413         if (emcal_upperbound > ihcal_lowerbound && emcal_lowerbound <= ihcal_lowerbound)
0414         {
0415           retower_lowerbound_originaltower_ieta[ieta_ihcal] = ieta_emcal;
0416           retower_lowerbound_originaltower_fraction[ieta_ihcal] = (emcal_upperbound - ihcal_lowerbound) / (emcal_upperbound - emcal_lowerbound);
0417           found_lowerbound = true;
0418         }
0419         if (emcal_upperbound > ihcal_lowerbound && emcal_lowerbound > ihcal_lowerbound)
0420         {
0421           retower_lowerbound_originaltower_ieta[ieta_ihcal] = ieta_emcal;
0422           retower_lowerbound_originaltower_fraction[ieta_ihcal] = 1;
0423           found_lowerbound = true;
0424         }
0425       }
0426       else
0427       {
0428         if (emcal_upperbound >= ihcal_upperbound && emcal_lowerbound < ihcal_upperbound)
0429         {
0430           retower_upperbound_originaltower_ieta[ieta_ihcal] = ieta_emcal;
0431           retower_upperbound_originaltower_fraction[ieta_ihcal] = (ihcal_upperbound - emcal_lowerbound) / (emcal_upperbound - emcal_lowerbound);
0432           found_upperbound = true;
0433         }
0434         if (emcal_upperbound > ihcal_upperbound && emcal_lowerbound > ihcal_upperbound)
0435         {
0436           ieta_emcal--;
0437           retower_upperbound_originaltower_ieta[ieta_ihcal] = ieta_emcal;
0438           retower_upperbound_originaltower_fraction[ieta_ihcal] = 1;
0439           found_upperbound = true;
0440         }
0441       }
0442       if (found_lowerbound && found_upperbound)
0443       {
0444         retower_totalarea[ieta_ihcal] = 4 * (retower_lowerbound_originaltower_fraction[ieta_ihcal] + retower_upperbound_originaltower_fraction[ieta_ihcal] + (retower_upperbound_originaltower_ieta[ieta_ihcal] - retower_lowerbound_originaltower_ieta[ieta_ihcal] - 1));
0445       }
0446       else
0447       {
0448         ieta_emcal++;
0449       }
0450     }
0451     if (!found_lowerbound)
0452     {
0453       std::cout << "RetowerCEMC::get_weighted_fraction: could not find lower bound EMCal towers for ieta_ihcal = " << ieta_ihcal << std::endl;
0454       exit(1);
0455     }
0456     else if (!found_upperbound)
0457     {
0458       std::cout << "RetowerCEMC::get_weighted_fraction: could not find upper bound EMCal towers for ieta_ihcal = " << ieta_ihcal << std::endl;
0459       exit(1);
0460     }
0461   }
0462 }