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