g4simulation/g4histos/G4EdepNtuple.cc

0001 #include "G4EdepNtuple.h"
0002
0003 #include <g4main/PHG4Hit.h>
0004 #include <g4main/PHG4HitContainer.h>
0005
0006 #include <fun4all/Fun4AllHistoManager.h>
0007 #include <fun4all/SubsysReco.h>  // for SubsysReco
0008
0009 #include <phool/getClass.h>
0010
0011 #include <TFile.h>
0012 #include <TNtuple.h>
0013
0014 #include <sstream>
0015 #include <utility>  // for pair
0016
0017 G4EdepNtuple::G4EdepNtuple(const std::string &name, const std::string &filename)
0018   : SubsysReco(name)
0019   , _filename(filename)
0020 {
0021 }
0022
0023 G4EdepNtuple::~G4EdepNtuple()
0024 {
0025   delete hm;
0026 }
0027
0028 int G4EdepNtuple::Init(PHCompositeNode * /*unused*/)
0029 {
0030   delete hm; // make cppcheck happy
0031   hm = new Fun4AllHistoManager(Name());
0032   outfile = new TFile(_filename.c_str(), "RECREATE");
0033   ntup = new TNtuple("edepntup", "G4Edeps", "detid:layer:edep");
0034   return 0;
0035 }
0036
0037 int G4EdepNtuple::process_event(PHCompositeNode *topNode)
0038 {
0039   std::string nodename;
0040   std::map<int, double> layer_edep_map;
0041   for (const auto &iter : _node_postfix)
0042   {
0043     layer_edep_map.clear();
0044     int detid = (_detid.find(iter))->second;
0045     nodename = "G4HIT_" + iter;
0046     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0047     if (hits)
0048     {
0049       double esum = 0;
0050       //          double numhits = hits->size();
0051       //          nhits[i]->Fill(numhits);
0052       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0053       for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0054
0055       {
0056         layer_edep_map[hit_iter->second->get_layer()] += hit_iter->second->get_edep();
0057         esum += hit_iter->second->get_edep();
0058       }
0059       for (auto edepiter : layer_edep_map)
0060       {
0061         ntup->Fill(detid, edepiter.first, edepiter.second);
0062       }
0063       ntup->Fill(detid, -1, esum);  // fill sum over all layers for each detector
0064     }
0065   }
0066   return 0;
0067 }
0068
0069 int G4EdepNtuple::End(PHCompositeNode * /*topNode*/)
0070 {
0071   outfile->cd();
0072   ntup->Write();
0073   outfile->Write();
0074   outfile->Close();
0075   delete outfile;
0076   hm->dumpHistos(_filename, "UPDATE");
0077   return 0;
0078 }
0079
0080 void G4EdepNtuple::AddNode(const std::string &name, const int detid)
0081 {
0082   _node_postfix.insert(name);
0083   _detid[name] = detid;
0084   return;
0085 }