g4simulation/g4caloprototype/Prototype2RawTowerBuilder.cc

0001 #include "Prototype2RawTowerBuilder.h"
0002 #include "PHG4PrototypeHcalDefs.h"
0003
0004 #include <calobase/RawTowerContainer.h>
0005 #include <calobase/RawTowerDefs.h>                      // for convert_name_...
0006 #include <calobase/RawTowerGeomContainer.h>             // for RawTowerGeomC...
0007 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0008 #include <calobase/RawTowerGeomv1.h>
0009 #include <calobase/RawTower.h>                          // for RawTower
0010 #include <calobase/RawTowerv1.h>
0011
0012 #include <g4detectors/PHG4ScintillatorSlat.h>
0013 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0014
0015 #include <phparameter/PHParameters.h>
0016 #include <phparameter/PHParameterInterface.h>           // for PHParameterIn...
0017
0018 #include <fun4all/Fun4AllReturnCodes.h>
0019 #include <fun4all/Fun4AllServer.h>
0020 #include <fun4all/SubsysReco.h>                         // for SubsysReco
0021
0022 #include <pdbcalbase/PdbParameterMapContainer.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHIODataNode.h>
0026 #include <phool/PHNode.h>                               // for PHNode
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/PHObject.h>                             // for PHObject
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>                                // for PHWHERE
0031
0032 #include <TSystem.h>
0033
0034 #include <cmath>                                       // for fabs, NAN
0035 #include <iostream>
0036 #include <map>
0037 #include <utility>                                      // for pair
0038
0039 using namespace std;
0040
0041 Prototype2RawTowerBuilder::Prototype2RawTowerBuilder(const std::string &name)
0042   : SubsysReco(name)
0043   , PHParameterInterface(name)
0044   , m_Detector("NONE")
0045   , m_Emin(NAN)
0046   , m_CheckEnergyConservationFlag(0)
0047   , m_TowerEnergySrc(kLightYield)
0048   , m_NumCellToTower(5)
0049
0050 {
0051   InitializeParameters();
0052 }
0053
0054 int Prototype2RawTowerBuilder::InitRun(PHCompositeNode *topNode)
0055 {
0056   PHNodeIterator iter(topNode);
0057   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0058   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0059   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0060   if (!runNode || !dstNode || !parNode)
0061   {
0062     cout << PHWHERE << "Top Nodes missing, quitting after printing node tree" << endl;
0063     Fun4AllServer *se = Fun4AllServer::instance();
0064     se->Print("NODETREE");
0065     gSystem->Exit(1);
0066   }
0067   string paramnodename = "G4TOWERPARAM_" + m_Detector;
0068   string geonodename = "G4TOWERGEO_" + m_Detector;
0069
0070   if (m_SimTowerNodePrefix.empty())
0071   {
0072     // no prefix, consistent with older convension
0073     m_TowerNodeName = "TOWER_" + m_Detector;
0074   }
0075   else
0076   {
0077     m_TowerNodeName = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
0078   }
0079   RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, m_TowerNodeName);
0080   if (!towers)
0081   {
0082     PHNodeIterator dstiter(dstNode);
0083     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
0084     if (!DetNode)
0085     {
0086       DetNode = new PHCompositeNode(m_Detector);
0087       dstNode->addNode(DetNode);
0088     }
0089     towers = new RawTowerContainer(RawTowerDefs::convert_name_to_caloid(m_Detector));
0090     PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(towers, m_TowerNodeName, "PHObject");
0091     DetNode->addNode(towerNode);
0092   }
0093   // order first default,
0094   // then parameter from g4detector on node tree
0095   ReadParamsFromNodeTree(topNode);
0096   // then macro setting
0097   UpdateParametersWithMacro();
0098   PHNodeIterator runIter(runNode);
0099   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_Detector));
0100   if (!RunDetNode)
0101   {
0102     RunDetNode = new PHCompositeNode(m_Detector);
0103     runNode->addNode(RunDetNode);
0104   }
0105   SaveToNodeTree(RunDetNode, paramnodename);
0106   // save this to the parNode for use
0107   PHNodeIterator parIter(parNode);
0108   PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", m_Detector));
0109   if (!ParDetNode)
0110   {
0111     ParDetNode = new PHCompositeNode(m_Detector);
0112     parNode->addNode(ParDetNode);
0113   }
0114   PutOnParNode(ParDetNode, geonodename);
0115
0116   m_TowerGeomNodeName = "TOWERGEOM_" + m_Detector;
0117   RawTowerGeomContainer *rawtowergeom = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
0118   if (!rawtowergeom)
0119   {
0120     rawtowergeom = new RawTowerGeomContainer_Cylinderv1(RawTowerDefs::convert_name_to_caloid(m_Detector));
0121     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(rawtowergeom, m_TowerGeomNodeName, "PHObject");
0122     runNode->addNode(newNode);
0123   }
0124   rawtowergeom->set_phibins(4);
0125   rawtowergeom->set_etabins(4);
0126   for (int irow = 0; irow < rawtowergeom->get_phibins(); irow++)
0127   {
0128     for (int icolumn = 0; icolumn < rawtowergeom->get_etabins(); icolumn++)
0129     {
0130       RawTowerGeomv1 *tg = new RawTowerGeomv1(RawTowerDefs::encode_towerid(RawTowerDefs::convert_name_to_caloid(m_Detector), icolumn, irow));
0131       tg->set_center_x(irow * 10 + icolumn);
0132       tg->set_center_y(irow * 10 + icolumn);
0133       tg->set_center_z(irow * 10 + icolumn);
0134
0135       rawtowergeom->add_tower_geometry(tg);
0136     }
0137   }
0138
0139   if (Verbosity() >= 1)
0140   {
0141     cout << "Prototype2RawTowerBuilder::InitRun :";
0142     if (m_TowerEnergySrc == kEnergyDeposition)
0143     {
0144       cout << "save Geant4 energy deposition as the weight of the cells"
0145            << endl;
0146     }
0147     else if (m_TowerEnergySrc == kLightYield)
0148     {
0149       cout << "save light yield as the weight of the cells" << endl;
0150     }
0151   }
0152   m_NumCellToTower = get_int_param(PHG4PrototypeHcalDefs::scipertwr);
0153   m_Emin = get_double_param("emin");
0154   return Fun4AllReturnCodes::EVENT_OK;
0155 }
0156
0157 int Prototype2RawTowerBuilder::process_event(PHCompositeNode *topNode)
0158 {
0159   if (Verbosity() > 3)
0160   {
0161     std::cout << PHWHERE << "Process event entered" << std::endl;
0162   }
0163
0164   // get cells
0165   string cellnodename = "G4CELL_" + m_Detector;
0166   PHG4ScintillatorSlatContainer *slats = findNode::getClass<PHG4ScintillatorSlatContainer>(topNode, cellnodename);
0167   if (!slats)
0168   {
0169     cout << PHWHERE << " " << cellnodename << " Node missing, doing nothing." << endl;
0170     return Fun4AllReturnCodes::ABORTEVENT;
0171   }
0172   RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, m_TowerNodeName);
0173   if (!towers)
0174   {
0175     cout << PHWHERE << " " << m_TowerNodeName << " Node missing, doing nothing." << endl;
0176     return Fun4AllReturnCodes::ABORTEVENT;
0177   }
0178
0179   // loop over all slats in an event
0180   PHG4ScintillatorSlatContainer::ConstIterator cell_iter;
0181   PHG4ScintillatorSlatContainer::ConstRange cell_range = slats->getScintillatorSlats();
0182   for (cell_iter = cell_range.first; cell_iter != cell_range.second;
0183        ++cell_iter)
0184   {
0185     PHG4ScintillatorSlat *cell = cell_iter->second;
0186
0187     if (Verbosity() > 2)
0188     {
0189       std::cout << PHWHERE << " print out the cell:" << std::endl;
0190       cell->identify();
0191     }
0192     short twrrow = get_tower_row(cell->get_row());
0193     // add the energy to the corresponding tower
0194     // towers are addressed column/row to make the mapping more intuitive
0195     RawTower *tower = towers->getTower(cell->get_column(), twrrow);
0196     if (!tower)
0197     {
0198       tower = new RawTowerv1();
0199       tower->set_energy(0);
0200       towers->AddTower(cell->get_column(), twrrow, tower);
0201     }
0202     double cell_weight = 0;
0203     if (m_TowerEnergySrc == kEnergyDeposition)
0204     {
0205       cell_weight = cell->get_edep();
0206     }
0207     else if (m_TowerEnergySrc == kLightYield)
0208     {
0209       cell_weight = cell->get_light_yield();
0210     }
0211     else if (m_TowerEnergySrc == kIonizationEnergy)
0212     {
0213       cell_weight = cell->get_eion();
0214     }
0215
0216     tower->add_ecell(cell->get_key(), cell_weight);
0217
0218     tower->set_energy(tower->get_energy() + cell_weight);
0219   }
0220   double towerE = 0;
0221   if (m_CheckEnergyConservationFlag)
0222   {
0223     double cellE = slats->getTotalEdep();
0224     towerE = towers->getTotalEdep();
0225     if (fabs(cellE - towerE) / cellE > 1e-5)
0226     {
0227       cout << "towerE: " << towerE << ", cellE: " << cellE << ", delta: "
0228            << cellE - towerE << endl;
0229     }
0230   }
0231   if (Verbosity())
0232   {
0233     towerE = towers->getTotalEdep();
0234   }
0235
0236   towers->compress(m_Emin);
0237   if (Verbosity())
0238   {
0239     cout << "Energy lost by dropping towers with less than " << m_Emin
0240          << " energy, lost energy: " << towerE - towers->getTotalEdep()
0241          << endl;
0242     towers->identify();
0243     RawTowerContainer::ConstRange begin_end = towers->getTowers();
0244     RawTowerContainer::ConstIterator iter;
0245     for (iter = begin_end.first; iter != begin_end.second; ++iter)
0246     {
0247       iter->second->identify();
0248     }
0249   }
0250
0251   return Fun4AllReturnCodes::EVENT_OK;
0252 }
0253
0254 short Prototype2RawTowerBuilder::get_tower_row(const short cellrow) const
0255 {
0256   short twrrow = cellrow / m_NumCellToTower;
0257   return twrrow;
0258 }
0259
0260 void Prototype2RawTowerBuilder::SetDefaultParameters()
0261 {
0262   set_default_int_param(PHG4PrototypeHcalDefs::scipertwr, 0);
0263   set_default_double_param("emin", 1.e-6);
0264   return;
0265 }
0266
0267 void Prototype2RawTowerBuilder::ReadParamsFromNodeTree(PHCompositeNode *topNode)
0268 {
0269   PHParameters *pars = new PHParameters("temp");
0270   // we need the number of scintillator plates per tower
0271   string geonodename = "G4GEOPARAM_" + m_Detector;
0272   PdbParameterMapContainer *saveparams = findNode::getClass<PdbParameterMapContainer>(topNode, geonodename);
0273   if (!saveparams)
0274   {
0275     cout << "could not find " << geonodename << endl;
0276     Fun4AllServer *se = Fun4AllServer::instance();
0277     se->Print("NODETREE");
0278     return;
0279   }
0280   pars->FillFrom(saveparams, 0);
0281   set_int_param(PHG4PrototypeHcalDefs::scipertwr, pars->get_int_param(PHG4PrototypeHcalDefs::scipertwr));
0282   delete pars;
0283   return;
0284 }
0285
0286 void Prototype2RawTowerBuilder::Print(const std::string &/*what*/) const
0287 {
0288   cout << "m_NumCellToTower: " << m_NumCellToTower << endl;
0289   return;
0290 }