g4simulation/g4histos/G4CellNtuple.cc

0001 #include "G4CellNtuple.h"
0002
0003 #include <g4detectors/PHG4Cell.h>
0004 #include <g4detectors/PHG4CellContainer.h>
0005 #include <g4detectors/PHG4CellDefs.h>  // for get_phibin
0006 #include <g4detectors/PHG4TpcCylinderGeom.h>
0007 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0008
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <fun4all/SubsysReco.h>  // for SubsysReco
0011
0012 #include <phool/getClass.h>
0013
0014 #include <TFile.h>
0015 #include <TH1.h>
0016 #include <TNtuple.h>
0017
0018 #include <cassert>
0019 #include <cmath>     // for isfinite
0020 #include <iostream>  // for operator<<
0021 #include <utility>   // for pair
0022
0023 G4CellNtuple::G4CellNtuple(const std::string &name, const std::string &filename)
0024   : SubsysReco(name)
0025   , _filename(filename)
0026 {
0027 }
0028
0029 G4CellNtuple::~G4CellNtuple()
0030 {
0031   delete hm;
0032 }
0033
0034 int G4CellNtuple::Init(PHCompositeNode * /*unused*/)
0035 {
0036   delete hm; // make cppcheck happy
0037   hm = new Fun4AllHistoManager(Name());
0038   outfile = new TFile(_filename.c_str(), "RECREATE");
0039   ntup = new TNtuple("cellntup", "G4Cells", "detid:layer:phi:eta:edep");
0040   //  ntup->SetDirectory(0);
0041   TH1 *h1 = new TH1F("edep1GeV", "edep 0-1GeV", 1000, 0, 1);
0042   eloss.push_back(h1);
0043   h1 = new TH1F("edep100GeV", "edep 0-100GeV", 1000, 0, 100);
0044   eloss.push_back(h1);
0045   return 0;
0046 }
0047
0048 int G4CellNtuple::process_event(PHCompositeNode *topNode)
0049 {
0050   std::string nodename;
0051   std::string geonodename;
0052   for (const auto &iter : _node_postfix)
0053   {
0054     int detid = (_detid.find(iter))->second;
0055     nodename = "G4CELL_" + iter;
0056     geonodename = "CYLINDERCELLGEOM_" + iter;
0057     PHG4TpcCylinderGeomContainer *cellgeos = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, geonodename);
0058     if (!cellgeos)
0059     {
0060       std::cout << "no geometry node " << geonodename << " for " << iter << std::endl;
0061       continue;
0062     }
0063     PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, nodename);
0064     if (cells)
0065     {
0066       int previouslayer = -1;
0067       PHG4TpcCylinderGeom *cellgeom = nullptr;
0068       double esum = 0;
0069       //          double numcells = cells->size();
0070       //          ncells[i]->Fill(numcells);
0071       //      std::cout << "number of cells: " << cells->size() << std::endl;
0072       PHG4CellContainer::ConstRange cell_range = cells->getCells();
0073       for (PHG4CellContainer::ConstIterator cell_iter = cell_range.first; cell_iter != cell_range.second; cell_iter++)
0074
0075       {
0076         double edep = cell_iter->second->get_edep();
0077         if (!std::isfinite(edep))
0078         {
0079           std::cout << "invalid edep: " << edep << std::endl;
0080         }
0081         esum += cell_iter->second->get_edep();
0082         int phibin = std::numeric_limits<int>::min();
0083         int etabin = std::numeric_limits<int>::min();
0084         double phi = std::numeric_limits<double>::quiet_NaN();
0085         double eta = std::numeric_limits<double>::quiet_NaN();
0086         int layer = cell_iter->second->get_layer();
0087         // to search the map fewer times, cache the geom object until the layer changes
0088         if (layer != previouslayer)
0089         {
0090           cellgeom = cellgeos->GetLayerCellGeom(layer);
0091           previouslayer = layer;
0092         }
0093         if (cell_iter->second->has_binning(PHG4CellDefs::etaphibinning))
0094         {
0095           phibin = PHG4CellDefs::EtaPhiBinning::get_phibin(cell_iter->second->get_cellid());
0096           etabin = PHG4CellDefs::EtaPhiBinning::get_etabin(cell_iter->second->get_cellid());
0097           phi = cellgeom->get_phicenter(phibin);
0098           eta = cellgeom->get_etacenter(etabin);
0099         }
0100         else if (cell_iter->second->has_binning(PHG4CellDefs::sizebinning))
0101         {
0102           phibin = PHG4CellDefs::SizeBinning::get_phibin(cell_iter->second->get_cellid());
0103           etabin = PHG4CellDefs::SizeBinning::get_zbin(cell_iter->second->get_cellid());
0104           phi = cellgeom->get_phicenter(phibin);
0105           eta = cellgeom->get_zcenter(etabin);
0106         }
0107         assert(cellgeom != nullptr);
0108         ntup->Fill(detid,
0109                    cell_iter->second->get_layer(),
0110                    phi,
0111                    eta,
0112                    cell_iter->second->get_edep());
0113       }
0114       for (auto *eiter : eloss)
0115       {
0116         eiter->Fill(esum);
0117       }
0118     }
0119   }
0120   return 0;
0121 }
0122
0123 int G4CellNtuple::End(PHCompositeNode * /*topNode*/)
0124 {
0125   outfile->cd();
0126   ntup->Write();
0127   outfile->Write();
0128   outfile->Close();
0129   delete outfile;
0130   hm->dumpHistos(_filename, "UPDATE");
0131   return 0;
0132 }
0133
0134 void G4CellNtuple::AddNode(const std::string &name, const int detid)
0135 {
0136   _node_postfix.insert(name);
0137   _detid[name] = detid;
0138   return;
0139 }