g4simulation/g4caloprototype/PHG4Prototype2HcalCellReco.cc

0001 #include "PHG4Prototype2HcalCellReco.h"
0002 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0003 #include <g4detectors/PHG4ScintillatorSlat.h>              // for PHG4ScintillatorSlat
0004 #include <g4detectors/PHG4ScintillatorSlatv1.h>
0005 #include <g4detectors/PHG4ScintillatorSlatDefs.h>          // for genkey, keytype
0006
0007 #include <g4main/PHG4Hit.h>
0008 #include <g4main/PHG4HitContainer.h>
0009
0010 #include <phparameter/PHParameterInterface.h>  // for PHParameterInterface
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/SubsysReco.h>                // for SubsysReco
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>
0018 #include <phool/PHNode.h>                      // for PHNode
0019 #include <phool/PHNodeIterator.h>
0020 #include <phool/PHObject.h>                    // for PHObject
0021 #include <phool/getClass.h>
0022 #include <phool/phool.h>                       // for PHWHERE
0023
0024 #include <TSystem.h>
0025
0026 #include <array>                               // for array, array<>::value_...
0027 #include <cmath>
0028 #include <cstdlib>
0029 #include <iostream>
0030 #include <map>                                 // for _Rb_tree_const_iterator
0031 #include <sstream>
0032 #include <utility>                             // for pair
0033
0034 using namespace std;
0035 // for hcal dimension
0036 // prototype outer hcal, 1st and 2nd prototype inner hcal, 3rd prototype inner hcal is 17
0037 #define ROWDIM 21
0038 // 12 scintillator tiles per row (index starting at 1)
0039 #define COLUMNDIM 13
0040 static array<array<PHG4ScintillatorSlat *, COLUMNDIM>, ROWDIM> slatarray = {{{nullptr}}};
0041
0042 PHG4Prototype2HcalCellReco::PHG4Prototype2HcalCellReco(const string &name)
0043   : SubsysReco(name)
0044   , PHParameterInterface(name)
0045   , m_CheckEnergyConservationFlag(0)
0046   , m_Tmin(NAN)
0047   , m_Tmax(NAN)
0048 {
0049   InitializeParameters();
0050 }
0051
0052 int PHG4Prototype2HcalCellReco::InitRun(PHCompositeNode *topNode)
0053 {
0054   PHNodeIterator iter(topNode);
0055
0056   // Looking for the DST node
0057   PHCompositeNode *dstNode;
0058   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0059   if (!dstNode)
0060   {
0061     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0062     gSystem->Exit(1);
0063     exit(1);
0064   }
0065   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0066   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0067
0068   string paramnodename = "G4CELLPARAM_" + m_Detector;
0069   string geonodename = "G4CELLGEO_" + m_Detector;
0070
0071   m_HitNodeName = "G4HIT_" + m_Detector;
0072   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0073   if (!g4hit)
0074   {
0075     cout << "Could not locate g4 hit node " << m_HitNodeName << endl;
0076     Fun4AllServer *se = Fun4AllServer::instance();
0077     se->Print("NODETREE");
0078     gSystem->Exit(1);
0079     exit(1);
0080   }
0081   m_CellNodeName = "G4CELL_" + m_Detector;
0082   PHG4ScintillatorSlatContainer *slats = findNode::getClass<PHG4ScintillatorSlatContainer>(topNode, m_CellNodeName);
0083   if (!slats)
0084   {
0085     PHNodeIterator dstiter(dstNode);
0086     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
0087     if (!DetNode)
0088     {
0089       DetNode = new PHCompositeNode(m_Detector);
0090       dstNode->addNode(DetNode);
0091     }
0092     slats = new PHG4ScintillatorSlatContainer();
0093     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(slats, m_CellNodeName, "PHObject");
0094     DetNode->addNode(newNode);
0095   }
0096   UpdateParametersWithMacro();
0097   PHNodeIterator runIter(runNode);
0098   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_Detector));
0099   if (!RunDetNode)
0100   {
0101     RunDetNode = new PHCompositeNode(m_Detector);
0102     runNode->addNode(RunDetNode);
0103   }
0104   SaveToNodeTree(RunDetNode, paramnodename);
0105   // save this to the parNode for use
0106   PHNodeIterator parIter(parNode);
0107   PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", m_Detector));
0108   if (!ParDetNode)
0109   {
0110     ParDetNode = new PHCompositeNode(m_Detector);
0111     parNode->addNode(ParDetNode);
0112   }
0113   PutOnParNode(ParDetNode, geonodename);
0114
0115   m_Tmin = get_double_param("tmin");
0116   m_Tmax = get_double_param("tmax");
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 int PHG4Prototype2HcalCellReco::process_event(PHCompositeNode *topNode)
0121 {
0122   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0123   if (!g4hit)
0124   {
0125     cout << "Could not locate g4 hit node " << m_HitNodeName << endl;
0126     exit(1);
0127   }
0128   PHG4ScintillatorSlatContainer *slats = findNode::getClass<PHG4ScintillatorSlatContainer>(topNode, m_CellNodeName);
0129   if (!slats)
0130   {
0131     cout << "could not locate cell node " << m_CellNodeName << endl;
0132     exit(1);
0133   }
0134
0135   PHG4HitContainer::LayerIter layer;
0136   pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
0137   for (layer = layer_begin_end.first; layer != layer_begin_end.second; ++layer)
0138   {
0139     PHG4HitContainer::ConstIterator hiter;
0140     PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
0141     for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0142     {
0143       if (hiter->second->get_t(0) > m_Tmax) continue;
0144       if (hiter->second->get_t(1) < m_Tmin) continue;
0145       short icolumn = hiter->second->get_scint_id();
0146       short irow = hiter->second->get_row();
0147       if (irow >= ROWDIM || irow < 0)
0148       {
0149         cout << "row " << irow
0150              << " exceed array size: " << ROWDIM
0151              << " adjust ROWDIM and recompile" << endl;
0152         gSystem->Exit(1);
0153       }
0154
0155       if (icolumn >= COLUMNDIM || icolumn < 0)
0156       {
0157         cout << "column: " << icolumn
0158              << " exceed array size: " << COLUMNDIM
0159              << " adjust COLUMNDIM and recompile" << endl;
0160         gSystem->Exit(1);
0161       }
0162
0163       if (!slatarray[irow][icolumn])
0164       {
0165         slatarray[irow][icolumn] = new PHG4ScintillatorSlatv1();
0166       }
0167       slatarray[irow][icolumn]->add_edep(hiter->second->get_edep(),
0168                                          hiter->second->get_eion(),
0169                                          hiter->second->get_light_yield());
0170       slatarray[irow][icolumn]->add_hit_key(hiter->first);
0171       // cout << "row: " << hiter->second->get_row()
0172       //       << ", column: " << hiter->second->get_scint_id() << endl;
0173       // checking ADC timing integration window cut
0174     }  // end loop over g4hits
0175     int nslathits = 0;
0176     for (int irow = 0; irow < ROWDIM; irow++)
0177     {
0178       for (int icolumn = 0; icolumn < COLUMNDIM; icolumn++)
0179       {
0180         if (slatarray[irow][icolumn])
0181         {
0182           PHG4ScintillatorSlatDefs::keytype key = PHG4ScintillatorSlatDefs::genkey(irow, icolumn);
0183           slats->AddScintillatorSlat(key, slatarray[irow][icolumn]);
0184           slatarray[irow][icolumn] = nullptr;
0185           nslathits++;
0186         }
0187       }
0188     }
0189     if (Verbosity() > 0)
0190     {
0191       cout << Name() << ": found " << nslathits << " slats with energy deposition" << endl;
0192     }
0193   }
0194
0195   if (m_CheckEnergyConservationFlag)
0196   {
0197     CheckEnergy(topNode);
0198   }
0199   return Fun4AllReturnCodes::EVENT_OK;
0200 }
0201
0202 int PHG4Prototype2HcalCellReco::CheckEnergy(PHCompositeNode *topNode)
0203 {
0204   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0205   PHG4ScintillatorSlatContainer *slats = findNode::getClass<PHG4ScintillatorSlatContainer>(topNode, m_CellNodeName);
0206   double sum_energy_g4hit = 0.;
0207   double sum_energy_cells = 0.;
0208   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0209   PHG4HitContainer::ConstIterator hiter;
0210   for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0211   {
0212     sum_energy_g4hit += hiter->second->get_edep();
0213   }
0214   PHG4ScintillatorSlatContainer::ConstRange cell_begin_end = slats->getScintillatorSlats();
0215   PHG4ScintillatorSlatContainer::ConstIterator citer;
0216   for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
0217   {
0218     sum_energy_cells += citer->second->get_edep();
0219   }
0220   // the fractional eloss for particles traversing eta bins leads to minute rounding errors
0221   if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
0222   {
0223     cout << "energy mismatch between cells: " << sum_energy_cells
0224          << " and hits: " << sum_energy_g4hit
0225          << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
0226          << endl;
0227     return -1;
0228   }
0229   else
0230   {
0231     if (Verbosity() > 0)
0232     {
0233       cout << Name() << ": total energy for this event: " << sum_energy_g4hit << " GeV" << endl;
0234     }
0235   }
0236   return 0;
0237 }
0238
0239 void PHG4Prototype2HcalCellReco::SetDefaultParameters()
0240 {
0241   set_default_double_param("tmax", 60.0);
0242   set_default_double_param("tmin", 0.0);
0243   return;
0244 }
0245
0246 void PHG4Prototype2HcalCellReco::set_timing_window(const double tmi, const double tma)
0247 {
0248   set_double_param("tmin", tmi);
0249   set_double_param("tmax", tma);
0250 }