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/EpdGeomV1.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 <cstdlib>
0031 #include <iostream>
0032 #include <map>      // for _Rb_tree_const_iterator
0033 #include <utility>  // for pair
0034
0035 PHG4EPDModuleReco::PHG4EPDModuleReco(const std::string &name)
0036   : SubsysReco(name)
0037   , PHParameterInterface(name)
0038 {
0039   InitializeParameters();
0040 }
0041
0042 int PHG4EPDModuleReco::InitRun(PHCompositeNode *topNode)
0043 {
0044   UpdateParametersWithMacro();
0045
0046   tmin = get_double_param("tmin");
0047   tmax = get_double_param("tmax");
0048   m_DeltaT = get_double_param("delta_t");
0049   m_EpdMpv = get_double_param("epdmpv");
0050
0051   CreateNodes(topNode);
0052   PHNodeIterator node_itr(topNode);
0053   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(node_itr.findFirst("PHCompositeNode", "RUN"));
0054   if (!runNode)
0055   {
0056     std::cout << PHWHERE << "RUN Node not found - that is fatal" << std::endl;
0057     gSystem->Exit(1);
0058     exit(1);
0059   }
0060
0061   PHNodeIterator runiter(runNode);
0062   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", m_Detector));
0063   if (!DetNode)
0064   {
0065     DetNode = new PHCompositeNode(m_Detector);
0066     runNode->addNode(DetNode);
0067   }
0068
0069   EpdGeom *epdGeom = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0070   if (!epdGeom)
0071   {
0072     epdGeom = new EpdGeomV1();
0073     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(epdGeom, "TOWERGEOM_EPD", "PHObject");
0074     DetNode->addNode(newNode);
0075   }
0076
0077   // fill epd geometry
0078   unsigned int epdchannels = 744;
0079   for (unsigned int ch = 0; ch < epdchannels; ch++)
0080   {
0081     unsigned int thiskey = TowerInfoDefs::encode_epd(ch);
0082     epdGeom->set_z(thiskey, GetTileZ(TowerInfoDefs::get_epd_arm(thiskey)));
0083     epdGeom->set_r(thiskey, GetTileR(TowerInfoDefs::get_epd_rbin(thiskey)));
0084     if (TowerInfoDefs::get_epd_rbin(thiskey) == 0)
0085     {
0086       epdGeom->set_phi0(thiskey, GetTilePhi0(TowerInfoDefs::get_epd_phibin(thiskey)));
0087     }
0088     else
0089     {
0090       epdGeom->set_phi(thiskey, GetTilePhi(TowerInfoDefs::get_epd_phibin(thiskey)));
0091     }
0092   }
0093
0094   return Fun4AllReturnCodes::EVENT_OK;
0095 }
0096
0097 int PHG4EPDModuleReco::process_event(PHCompositeNode *topNode)
0098 {
0099   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_Hitnodename);
0100   if (!g4hit)
0101   {
0102     std::cout << "Could not locate g4 hit node " << m_Hitnodename << std::endl;
0103     exit(1);
0104   }
0105
0106   TowerInfoContainer *m_TowerInfoContainer = findNode::getClass<TowerInfoContainer>(topNode, m_TowerInfoNodeName);
0107   if (!m_TowerInfoContainer)
0108   {
0109     std::cout << PHWHERE << "Could not locate TowerInfoContainer node " << m_TowerInfoNodeName << std::endl;
0110     exit(1);
0111   }
0112
0113   TowerInfoContainer *m_TowerInfoContainer_calib = findNode::getClass<TowerInfoContainer>(topNode, m_TowerInfoNodeName_calib);
0114   if (!m_TowerInfoContainer_calib)
0115   {
0116     std::cout << PHWHERE << "Could not locate TowerInfoContainer node " << m_TowerInfoNodeName_calib << std::endl;
0117     exit(1);
0118   }
0119
0120   PHG4HitContainer::ConstIterator hiter;
0121   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0122   for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0123   {
0124     if (hiter->second->get_t(0) > tmax)
0125     {
0126       continue;
0127     }
0128     if (hiter->second->get_t(1) < tmin)
0129     {
0130       continue;
0131     }
0132     if (hiter->second->get_t(1) - hiter->second->get_t(0) > m_DeltaT)
0133     {
0134       continue;
0135     }
0136
0137     int sim_tileid = (hiter->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
0138     int this_tile = EPDDefs::get_tileid(sim_tileid);
0139     int this_sector = EPDDefs::get_sector(sim_tileid);
0140     int this_arm = EPDDefs::get_arm(sim_tileid);
0141
0142     m_EpdTile_e[this_arm][this_sector][this_tile] += hiter->second->get_light_yield();
0143     m_EpdTile_Calib_e[this_arm][this_sector][this_tile] += hiter->second->get_light_yield() / m_EpdMpv;
0144
0145   }  // end loop over g4hits
0146
0147   for (unsigned int k = 0; k < 2; k++)
0148   {
0149     for (int i = 0; i < 12; i++)
0150     {
0151       for (int j = 0; j < 31; j++)
0152       {
0153         unsigned int globalphi = Getphimap(j) + 2 * i;
0154         unsigned int r = Getrmap(j);
0155
0156         if (r == 0)
0157         {
0158           globalphi = i;
0159         }
0160
0161         unsigned int key = TowerInfoDefs::encode_epd(k, r, globalphi);
0162         unsigned int ch = m_TowerInfoContainer->decode_key(key);
0163         m_TowerInfoContainer->get_tower_at_channel(ch)->set_energy(m_EpdTile_e[k][i][j]);
0164         m_TowerInfoContainer_calib->get_tower_at_channel(ch)->set_energy(m_EpdTile_Calib_e[k][i][j]);
0165       }
0166     }
0167   }
0168
0169   return Fun4AllReturnCodes::EVENT_OK;
0170 }
0171
0172 void PHG4EPDModuleReco::SetDefaultParameters()
0173 {
0174   set_default_double_param("tmax", 60.0);
0175   set_default_double_param("tmin", -20.0);
0176   set_default_double_param("delta_t", 100.);
0177   set_default_double_param("epdmpv", 1.);
0178   return;
0179 }
0180
0181 void PHG4EPDModuleReco::set_timing_window(const double tmi, const double tma)
0182 {
0183   set_double_param("tmin", tmi);
0184   set_double_param("tmax", tma);
0185 }
0186
0187 int PHG4EPDModuleReco::Getrmap(int rindex)
0188 {
0189   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};
0190
0191   return rmap[rindex];
0192 }
0193
0194 int PHG4EPDModuleReco::Getphimap(int phiindex)
0195 {
0196   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};
0197
0198   return phimap[phiindex];
0199 }
0200
0201 float PHG4EPDModuleReco::GetTilePhi(int thisphi)
0202 {
0203   static const float tilephi[24] = {0.13089969, 0.39269908, 0.65449847, 0.91629786, 1.17809725, 1.43989663, 1.70169602, 1.96349541, 2.2252948, 2.48709418, 2.74889357, 3.01069296, 3.27249235, 3.53429174, 3.79609112, 4.05789051, 4.3196899, 4.58148929, 4.84328867, 5.10508806, 5.36688745, 5.62868684, 5.89048623, 6.15228561};
0204   return tilephi[thisphi];
0205 }
0206
0207 float PHG4EPDModuleReco::GetTilePhi0(int thisphi0)
0208 {
0209   static const float tilephi0[12] = {0.26179939, 0.78539816, 1.30899694, 1.83259571, 2.35619449, 2.87979327, 3.40339204, 3.92699082, 4.45058959, 4.97418837, 5.49778714, 6.02138592};
0210   return tilephi0[thisphi0];
0211 }
0212
0213 float PHG4EPDModuleReco::GetTileR(int thisr)
0214 {
0215   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};
0216   return tileR[thisr];
0217 }
0218
0219 float PHG4EPDModuleReco::GetTileZ(int thisz)
0220 {
0221   static const float tileZ[2] = {-316.0, 316.0};
0222   return tileZ[thisz];
0223 }
0224
0225 void PHG4EPDModuleReco::CreateNodes(PHCompositeNode *topNode)
0226 {
0227   PHNodeIterator iter(topNode);
0228   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0229   if (!dstNode)
0230   {
0231     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0232     gSystem->Exit(1);
0233     exit(1);
0234   }
0235
0236   PHNodeIterator dstiter(dstNode);
0237   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
0238   if (!DetNode)
0239   {
0240     DetNode = new PHCompositeNode(m_Detector);
0241     dstNode->addNode(DetNode);
0242   }
0243
0244   m_TowerInfoNodeName = "TOWERINFO_" + m_EPDSimTowerNodePrefix + "_" + m_Detector;  // detector name and prefix are set by now
0245   TowerInfoContainer *m_TowerInfoContainer = findNode::getClass<TowerInfoContainer>(DetNode, m_TowerInfoNodeName);
0246   if (m_TowerInfoContainer == nullptr)
0247   {
0248     m_TowerInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::SEPD);
0249     PHIODataNode<PHObject> *TowerInfoNode = new PHIODataNode<PHObject>(m_TowerInfoContainer, m_TowerInfoNodeName, "PHObject");
0250     DetNode->addNode(TowerInfoNode);
0251   }
0252
0253   m_TowerInfoNodeName_calib = "TOWERINFO_" + m_EPDCalibTowerNodePrefix + "_" + m_Detector;  // detector name and prefix are set by now
0254   TowerInfoContainer *m_TowerInfoContainer_calib = findNode::getClass<TowerInfoContainer>(DetNode, m_TowerInfoNodeName_calib);
0255   if (m_TowerInfoContainer_calib == nullptr)
0256   {
0257     m_TowerInfoContainer_calib = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::SEPD);
0258     PHIODataNode<PHObject> *TowerInfoNodecalib = new PHIODataNode<PHObject>(m_TowerInfoContainer_calib, m_TowerInfoNodeName_calib, "PHObject");
0259     DetNode->addNode(TowerInfoNodecalib);
0260   }
0261
0262   return;
0263 }
0264 int PHG4EPDModuleReco::ResetEvent(PHCompositeNode * /*topNode*/)
0265 {
0266   // this only works for initializing to zero
0267   m_EpdTile_Calib_e = {};
0268   m_EpdTile_e = {};
0269   // if you ever want to initialize to a different value, do it this way:
0270   //   for (auto &entry : epd_tile_calib_e)
0271   //   {
0272   //     for (auto &entry1 : entry)
0273   //     {
0274   // entry1.fill(NAN);
0275   //     }
0276   //   }
0277   return Fun4AllReturnCodes::EVENT_OK;
0278 }
0279
0280 void PHG4EPDModuleReco::Detector(const std::string &detector)
0281 {
0282   m_Detector = detector;
0283   m_Hitnodename = "G4HIT_" + m_Detector;
0284 }