g4simulation/g4detectors/PHG4HcalCellReco.cc

0001 #include "PHG4HcalCellReco.h"
0002
0003 #include "PHG4Cell.h"  // for PHG4Cell
0004 #include "PHG4CellContainer.h"
0005 #include "PHG4CellDefs.h"  // for genkey, keytype
0006 #include "PHG4Cellv1.h"
0007
0008 #include <phparameter/PHParameterInterface.h>  // for PHParameterInterface
0009
0010 #include <g4main/PHG4Hit.h>
0011 #include <g4main/PHG4HitContainer.h>
0012 #include <g4main/PHG4HitDefs.h>  // for hit_idbits
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/SubsysReco.h>  // for SubsysReco
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>
0019 #include <phool/PHNode.h>  // for PHNode
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/PHObject.h>  // for PHObject
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>  // for PHWHERE
0024
0025 #include <TSystem.h>
0026
0027 #include <array>  // for array, array<>::value_...
0028 #include <cmath>
0029 #include <cstdlib>
0030 #include <iostream>
0031 #include <map>  // for _Rb_tree_const_iterator
0032 #include <sstream>
0033 #include <utility>  // for pair
0034
0035 // for hcal dimension
0036 #define ROWDIM 320
0037 #define COLUMNDIM 27
0038
0039 static std::array<std::array<PHG4Cell *, COLUMNDIM>, ROWDIM> slatarray = {{{nullptr}}};
0040
0041 PHG4HcalCellReco::PHG4HcalCellReco(const std::string &name)
0042   : SubsysReco(name)
0043   , PHParameterInterface(name)
0044 {
0045   InitializeParameters();
0046 }
0047
0048 int PHG4HcalCellReco::InitRun(PHCompositeNode *topNode)
0049 {
0050   PHNodeIterator iter(topNode);
0051
0052   // Looking for the DST node
0053   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0054   if (!dstNode)
0055   {
0056     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0057     exit(1);
0058   }
0059   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0060   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0061
0062   std::string paramnodename = "G4CELLPARAM_" + detector;
0063   std::string geonodename = "G4CELLGEO_" + detector;
0064   hitnodename = "G4HIT_" + detector;
0065   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0066   if (!g4hit)
0067   {
0068     std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
0069     topNode->print();
0070     gSystem->Exit(1);
0071     exit(1);
0072   }
0073   cellnodename = "G4CELL_" + detector;
0074   PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0075   if (!slats)
0076   {
0077     PHNodeIterator dstiter(dstNode);
0078     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
0079     if (!DetNode)
0080     {
0081       DetNode = new PHCompositeNode(detector);
0082       dstNode->addNode(DetNode);
0083     }
0084     slats = new PHG4CellContainer();
0085     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(slats, cellnodename, "PHObject");
0086     DetNode->addNode(newNode);
0087   }
0088   UpdateParametersWithMacro();
0089   // save this to the run wise tree to store on DST
0090   PHNodeIterator runIter(runNode);
0091   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
0092   if (!RunDetNode)
0093   {
0094     RunDetNode = new PHCompositeNode(detector);
0095     runNode->addNode(RunDetNode);
0096   }
0097   SaveToNodeTree(RunDetNode, paramnodename);
0098   // save this to the parNode for use
0099   PHNodeIterator parIter(parNode);
0100   PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
0101   if (!ParDetNode)
0102   {
0103     ParDetNode = new PHCompositeNode(detector);
0104     parNode->addNode(ParDetNode);
0105   }
0106   PutOnParNode(ParDetNode, geonodename);
0107   tmin = get_double_param("tmin");
0108   tmax = get_double_param("tmax");
0109   m_DeltaT = get_double_param("delta_t");
0110   return Fun4AllReturnCodes::EVENT_OK;
0111 }
0112
0113 int PHG4HcalCellReco::process_event(PHCompositeNode *topNode)
0114 {
0115   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0116   if (!g4hit)
0117   {
0118     std::cout << "Could not locate g4 hit node " << hitnodename << std::endl;
0119     exit(1);
0120   }
0121   PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0122   if (!slats)
0123   {
0124     std::cout << "could not locate cell node " << cellnodename << std::endl;
0125     exit(1);
0126   }
0127   if (std::isfinite(m_FixedEnergy))
0128   {
0129     int maxcolumn = 24;
0130     int maxrow = 320;
0131     if (detector == "HCALIN")
0132     {
0133       maxrow = 256;
0134     }
0135     for (int icolumn = 0; icolumn < maxcolumn; icolumn++)
0136     {
0137       for (int irow = 0; irow < maxrow; irow++)
0138       {
0139         PHG4CellDefs::keytype key = PHG4CellDefs::ScintillatorSlatBinning::genkey(0, icolumn, irow);
0140         PHG4Cell *cell = new PHG4Cellv1(key);
0141         cell->add_edep(m_FixedEnergy);
0142         cell->add_eion(m_FixedEnergy);
0143         cell->add_light_yield(m_FixedEnergy);
0144         cell->add_raw_light_yield(m_FixedEnergy);
0145         slats->AddCell(cell);
0146       }
0147     }
0148     return Fun4AllReturnCodes::EVENT_OK;
0149   }
0150   PHG4HitContainer::ConstIterator hiter;
0151   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0152   for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0153   {
0154     if (hiter->second->get_t(0) > tmax)
0155     {
0156       continue;
0157     }
0158     if (hiter->second->get_t(1) < tmin)
0159     {
0160       continue;
0161     }
0162     if (hiter->second->get_t(1) - hiter->second->get_t(0) > m_DeltaT)
0163     {
0164       continue;
0165     }
0166
0167     short icolumn = hiter->second->get_scint_id();
0168     int introw = (hiter->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
0169     if (introw >= ROWDIM || introw < 0)
0170     {
0171       std::cout << __PRETTY_FUNCTION__ << " row " << introw
0172                 << " exceed array size: " << ROWDIM
0173                 << " adjust ROWDIM and recompile" << std::endl;
0174       exit(1);
0175     }
0176     // after checking for size of introw so we do not run into
0177     // overflow issues, put this into the short we want later
0178     short irow = introw;
0179     if (icolumn >= COLUMNDIM || icolumn < 0)
0180     {
0181       std::cout << __PRETTY_FUNCTION__ << " column: " << icolumn
0182                 << " exceed array size: " << COLUMNDIM
0183                 << " adjust COLUMNDIM and recompile" << std::endl;
0184       exit(1);
0185     }
0186
0187     if (!slatarray[irow][icolumn])
0188     {
0189       // hcal has no layers so far, I do not want to make an expensive
0190       // call to the g4hits to find that out use 0 as layer number
0191       PHG4CellDefs::keytype key = PHG4CellDefs::ScintillatorSlatBinning::genkey(0, icolumn, irow);
0192       slatarray[irow][icolumn] = new PHG4Cellv1(key);
0193     }
0194     slatarray[irow][icolumn]->add_edep(hiter->second->get_edep());
0195     slatarray[irow][icolumn]->add_eion(hiter->second->get_eion());
0196     slatarray[irow][icolumn]->add_light_yield(hiter->second->get_light_yield());
0197     float raw_light = hiter->second->get_raw_light_yield();
0198     if (std::isfinite(raw_light))
0199     {
0200       slatarray[irow][icolumn]->add_raw_light_yield(raw_light);
0201     }
0202     slatarray[irow][icolumn]->add_edep(hiter->first, hiter->second->get_edep());
0203     slatarray[irow][icolumn]->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep());
0204   }  // end loop over g4hits
0205   int nslathits = 0;
0206   for (int irow = 0; irow < ROWDIM; irow++)
0207   {
0208     for (int icolumn = 0; icolumn < COLUMNDIM; icolumn++)
0209     {
0210       if (slatarray[irow][icolumn])
0211       {
0212         slats->AddCell(slatarray[irow][icolumn]);
0213         slatarray[irow][icolumn] = nullptr;
0214         nslathits++;
0215       }
0216     }
0217   }
0218   if (Verbosity() > 0)
0219   {
0220     std::cout << Name() << ": found " << nslathits << " slats with energy deposition" << std::endl;
0221   }
0222
0223   if (chkenergyconservation)
0224   {
0225     CheckEnergy(topNode);
0226   }
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229
0230 int PHG4HcalCellReco::CheckEnergy(PHCompositeNode *topNode)
0231 {
0232   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0233   PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0234   double sum_energy_g4hit = 0.;
0235   double sum_energy_cells = 0.;
0236   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0237   PHG4HitContainer::ConstIterator hiter;
0238   for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0239   {
0240     sum_energy_g4hit += hiter->second->get_edep();
0241   }
0242   PHG4CellContainer::ConstRange cell_begin_end = slats->getCells();
0243   PHG4CellContainer::ConstIterator citer;
0244   for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
0245   {
0246     sum_energy_cells += citer->second->get_edep();
0247   }
0248   // the fractional eloss for particles traversing eta bins leads to minute rounding errors
0249   if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
0250   {
0251     std::cout << "hint: timing cuts tmin/tmax will do this to you" << std::endl;
0252     std::cout << "energy mismatch between cells: " << sum_energy_cells
0253               << " and hits: " << sum_energy_g4hit
0254               << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
0255               << std::endl;
0256     return -1;
0257   }
0258   else
0259   {
0260     if (Verbosity() > 0)
0261     {
0262       std::cout << Name() << ": total energy for this event: " << sum_energy_g4hit << " GeV" << std::endl;
0263     }
0264   }
0265   return 0;
0266 }
0267
0268 void PHG4HcalCellReco::SetDefaultParameters()
0269 {
0270   set_default_double_param("tmax", 60.0);
0271   set_default_double_param("tmin", -20.0);  // collision has a timing spread around the triggered event. Accepting negative time too.
0272   set_default_double_param("delta_t", 100.);
0273   return;
0274 }
0275
0276 void PHG4HcalCellReco::set_timing_window(const double tmi, const double tma)
0277 {
0278   set_double_param("tmin", tmi);
0279   set_double_param("tmax", tma);
0280 }