g4simulation/g4epd/PHG4EPDModuleReco.cc

0001 #include "PHG4EPDModuleReco.h"
0002
0003 #include <calobase/TowerInfo.h>
0004 #include <calobase/TowerInfoContainer.h>
0005 #include <calobase/TowerInfoContainerv1.h>
0006 #include <calobase/TowerInfoDefs.h>
0007
0008 #include <epd/EPDDefs.h>
0009 #include <epd/EpdGeomV2.h>
0010
0011 #include <g4main/PHG4Hit.h>
0012 #include <g4main/PHG4HitContainer.h>
0013 #include <g4main/PHG4HitDefs.h>  // for hit_idbits
0014
0015 #include <phparameter/PHParameterInterface.h>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <fun4all/SubsysReco.h>  // for SubsysReco
0019
0020 #include <phool/PHCompositeNode.h>
0021 #include <phool/PHIODataNode.h>
0022 #include <phool/PHNode.h>  // for PHNode
0023 #include <phool/PHNodeIterator.h>
0024 #include <phool/PHObject.h>  // for PHObject
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>  // for PHWHERE
0027
0028 #include <TSystem.h>
0029
0030 #include <cmath>
0031 #include <cstdlib>
0032 #include <iostream>
0033 #include <map>      // for _Rb_tree_const_iterator
0034 #include <utility>  // for pair
0035
0036 PHG4EPDModuleReco::PHG4EPDModuleReco(const std::string &name)
0037   : SubsysReco(name)
0038   , PHParameterInterface(name)
0039 {
0040   InitializeParameters();
0041   FillTilePhiArray();
0042   FillTilePhi0Array();
0043 }
0044
0045 int PHG4EPDModuleReco::InitRun(PHCompositeNode *topNode)
0046 {
0047   UpdateParametersWithMacro();
0048
0049   tmin = get_double_param("tmin");
0050   tmax = get_double_param("tmax");
0051   m_DeltaT = get_double_param("delta_t");
0052   m_EpdMpv = get_double_param("epdmpv");
0053
0054   CreateNodes(topNode);
0055   PHNodeIterator node_itr(topNode);
0056   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(node_itr.findFirst("PHCompositeNode", "RUN"));
0057   if (!runNode)
0058   {
0059     std::cout << PHWHERE << "RUN Node not found - that is fatal" << std::endl;
0060     gSystem->Exit(1);
0061     exit(1);
0062   }
0063
0064   PHNodeIterator runiter(runNode);
0065   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", m_Detector));
0066   if (!DetNode)
0067   {
0068     DetNode = new PHCompositeNode(m_Detector);
0069     runNode->addNode(DetNode);
0070   }
0071
0072   EpdGeom *epdGeom = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0073   if (!epdGeom)
0074   {
0075     epdGeom = new EpdGeomV2();
0076     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(epdGeom, "TOWERGEOM_EPD", "PHObject");
0077     DetNode->addNode(newNode);
0078   }
0079
0080   // fill epd geometry
0081   unsigned int epdchannels = 744;
0082   for (unsigned int ch = 0; ch < epdchannels; ch++)
0083   {
0084     unsigned int thiskey = TowerInfoDefs::encode_epd(ch);
0085     epdGeom->set_z(thiskey, GetTileZ(TowerInfoDefs::get_epd_arm(thiskey)));
0086     epdGeom->set_r(thiskey, GetTileR(TowerInfoDefs::get_epd_rbin(thiskey)));
0087     if (TowerInfoDefs::get_epd_rbin(thiskey) == 0)
0088     {
0089       epdGeom->set_phi0(thiskey, GetTilePhi0(TowerInfoDefs::get_epd_phibin(thiskey)));
0090     }
0091     else
0092     {
0093       epdGeom->set_phi(thiskey, GetTilePhi(TowerInfoDefs::get_epd_phibin(thiskey)));
0094     }
0095   }
0096
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 int PHG4EPDModuleReco::process_event(PHCompositeNode *topNode)
0101 {
0102   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_Hitnodename);
0103   if (!g4hit)
0104   {
0105     std::cout << "Could not locate g4 hit node " << m_Hitnodename << std::endl;
0106     exit(1);
0107   }
0108
0109   TowerInfoContainer *m_TowerInfoContainer = findNode::getClass<TowerInfoContainer>(topNode, m_TowerInfoNodeName);
0110   if (!m_TowerInfoContainer)
0111   {
0112     std::cout << PHWHERE << "Could not locate TowerInfoContainer node " << m_TowerInfoNodeName << std::endl;
0113     exit(1);
0114   }
0115
0116   TowerInfoContainer *m_TowerInfoContainer_calib = findNode::getClass<TowerInfoContainer>(topNode, m_TowerInfoNodeName_calib);
0117   if (!m_TowerInfoContainer_calib)
0118   {
0119     std::cout << PHWHERE << "Could not locate TowerInfoContainer node " << m_TowerInfoNodeName_calib << std::endl;
0120     exit(1);
0121   }
0122
0123   PHG4HitContainer::ConstIterator hiter;
0124   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0125   for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0126   {
0127     if (hiter->second->get_t(0) > tmax)
0128     {
0129       continue;
0130     }
0131     if (hiter->second->get_t(1) < tmin)
0132     {
0133       continue;
0134     }
0135     if (hiter->second->get_t(1) - hiter->second->get_t(0) > m_DeltaT)
0136     {
0137       continue;
0138     }
0139
0140     int sim_tileid = (hiter->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
0141     int this_tile = EPDDefs::get_tileid(sim_tileid);
0142     int this_sector = EPDDefs::get_sector(sim_tileid);
0143     int this_arm = EPDDefs::get_arm(sim_tileid);
0144
0145     m_EpdTile_e[this_arm][this_sector][this_tile] += hiter->second->get_light_yield();
0146     m_EpdTile_Calib_e[this_arm][this_sector][this_tile] += hiter->second->get_light_yield() / m_EpdMpv;
0147
0148   }  // end loop over g4hits
0149
0150   for (unsigned int k = 0; k < 2; k++)
0151   {
0152     for (int i = 0; i < 12; i++)
0153     {
0154       for (int j = 0; j < 31; j++)
0155       {
0156         unsigned int globalphi = Getphimap(j) + 2 * i;
0157         unsigned int r = Getrmap(j);
0158
0159         if (r == 0)
0160         {
0161           globalphi = i;
0162         }
0163
0164         unsigned int key = TowerInfoDefs::encode_epd(k, r, globalphi);
0165         unsigned int ch = m_TowerInfoContainer->decode_key(key);
0166         m_TowerInfoContainer->get_tower_at_channel(ch)->set_energy(m_EpdTile_e[k][i][j]);
0167         m_TowerInfoContainer_calib->get_tower_at_channel(ch)->set_energy(m_EpdTile_Calib_e[k][i][j]);
0168       }
0169     }
0170   }
0171
0172   return Fun4AllReturnCodes::EVENT_OK;
0173 }
0174
0175 void PHG4EPDModuleReco::SetDefaultParameters()
0176 {
0177   set_default_double_param("tmax", 60.0);
0178   set_default_double_param("tmin", -20.0);
0179   set_default_double_param("delta_t", 100.);
0180   set_default_double_param("epdmpv", 1.);
0181   return;
0182 }
0183
0184 void PHG4EPDModuleReco::set_timing_window(const double tmi, const double tma)
0185 {
0186   set_double_param("tmin", tmi);
0187   set_double_param("tmax", tma);
0188 }
0189
0190 int PHG4EPDModuleReco::Getrmap(int rindex)
0191 {
0192   static const int rmap[31] = {0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15};
0193
0194   return rmap[rindex];
0195 }
0196
0197 int PHG4EPDModuleReco::Getphimap(int phiindex)
0198 {
0199   static const int phimap[31] = {0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1};
0200
0201   return phimap[phiindex];
0202 }
0203
0204 float PHG4EPDModuleReco::GetTileR(int thisr)
0205 {
0206   static const float tileR[16] = {6.8, 11.2, 15.6, 20.565, 26.095, 31.625, 37.155, 42.685, 48.215, 53.745, 59.275, 64.805, 70.335, 75.865, 81.395, 86.925};
0207   return tileR[thisr];
0208 }
0209
0210 float PHG4EPDModuleReco::GetTileZ(int thisz)
0211 {
0212   static const float tileZ[2] = {-316.0, 316.0};
0213   return tileZ[thisz];
0214 }
0215
0216 void PHG4EPDModuleReco::CreateNodes(PHCompositeNode *topNode)
0217 {
0218   PHNodeIterator iter(topNode);
0219   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0220   if (!dstNode)
0221   {
0222     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0223     gSystem->Exit(1);
0224     exit(1);
0225   }
0226
0227   PHNodeIterator dstiter(dstNode);
0228   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
0229   if (!DetNode)
0230   {
0231     DetNode = new PHCompositeNode(m_Detector);
0232     dstNode->addNode(DetNode);
0233   }
0234
0235   m_TowerInfoNodeName = "TOWERINFO_" + m_EPDSimTowerNodePrefix + "_" + m_Detector;  // detector name and prefix are set by now
0236   TowerInfoContainer *m_TowerInfoContainer = findNode::getClass<TowerInfoContainer>(DetNode, m_TowerInfoNodeName);
0237   if (m_TowerInfoContainer == nullptr)
0238   {
0239     m_TowerInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::SEPD);
0240     PHIODataNode<PHObject> *TowerInfoNode = new PHIODataNode<PHObject>(m_TowerInfoContainer, m_TowerInfoNodeName, "PHObject");
0241     DetNode->addNode(TowerInfoNode);
0242   }
0243
0244   m_TowerInfoNodeName_calib = "TOWERINFO_" + m_EPDCalibTowerNodePrefix + "_" + m_Detector;  // detector name and prefix are set by now
0245   TowerInfoContainer *m_TowerInfoContainer_calib = findNode::getClass<TowerInfoContainer>(DetNode, m_TowerInfoNodeName_calib);
0246   if (m_TowerInfoContainer_calib == nullptr)
0247   {
0248     m_TowerInfoContainer_calib = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::SEPD);
0249     PHIODataNode<PHObject> *TowerInfoNodecalib = new PHIODataNode<PHObject>(m_TowerInfoContainer_calib, m_TowerInfoNodeName_calib, "PHObject");
0250     DetNode->addNode(TowerInfoNodecalib);
0251   }
0252
0253   return;
0254 }
0255 int PHG4EPDModuleReco::ResetEvent(PHCompositeNode * /*topNode*/)
0256 {
0257   // this only works for initializing to zero
0258   m_EpdTile_Calib_e = {};
0259   m_EpdTile_e = {};
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263 void PHG4EPDModuleReco::Detector(const std::string &detector)
0264 {
0265   m_Detector = detector;
0266   m_Hitnodename = "G4HIT_" + m_Detector;
0267 }
0268
0269 void PHG4EPDModuleReco::FillTilePhiArray()
0270 {
0271   size_t i = 0;
0272   for (float &iter : m_tilephi)
0273   {
0274     iter = ((2 * i + 1) * M_PI) / 24;
0275     i++;
0276   }
0277   return;
0278 }
0279
0280 void PHG4EPDModuleReco::FillTilePhi0Array()
0281 {
0282   size_t i = 0;
0283   for (float &iter : m_tilephi0)
0284   {
0285     iter = ((2 * i + 1) * M_PI) / 12;
0286     i++;
0287   }
0288   return;
0289 }
0290
0291 float PHG4EPDModuleReco::GetTilePhi(int thisphi)
0292 {
0293   return m_tilephi.at(thisphi);
0294 }
0295
0296 float PHG4EPDModuleReco::GetTilePhi0(int thisphi0)
0297 {
0298   return m_tilephi0.at(thisphi0);
0299 }