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 using namespace std;
0018
0019 G4EdepNtuple::G4EdepNtuple(const std::string &name, const std::string &filename)
0020   : SubsysReco(name)
0021   , nblocks(0)
0022   , hm(nullptr)
0023   , _filename(filename)
0024   , ntup(nullptr)
0025   , outfile(nullptr)
0026 {
0027 }
0028
0029 G4EdepNtuple::~G4EdepNtuple()
0030 {
0031   delete hm;
0032 }
0033
0034 int G4EdepNtuple::Init(PHCompositeNode * /*unused*/)
0035 {
0036   hm = new Fun4AllHistoManager(Name());
0037   outfile = new TFile(_filename.c_str(), "RECREATE");
0038   ntup = new TNtuple("edepntup", "G4Edeps", "detid:layer:edep");
0039   return 0;
0040 }
0041
0042 int G4EdepNtuple::process_event(PHCompositeNode *topNode)
0043 {
0044   ostringstream nodename;
0045   set<string>::const_iterator iter;
0046   map<int, double> layer_edep_map;
0047   map<int, double>::const_iterator edepiter;
0048   for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
0049   {
0050     layer_edep_map.clear();
0051     int detid = (_detid.find(*iter))->second;
0052     nodename.str("");
0053     nodename << "G4HIT_" << *iter;
0054     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
0055     if (hits)
0056     {
0057       double esum = 0;
0058       //          double numhits = hits->size();
0059       //          nhits[i]->Fill(numhits);
0060       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0061       for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0062
0063       {
0064         layer_edep_map[hit_iter->second->get_layer()] += hit_iter->second->get_edep();
0065         esum += hit_iter->second->get_edep();
0066       }
0067       for (edepiter = layer_edep_map.begin(); edepiter != layer_edep_map.end(); ++edepiter)
0068       {
0069         ntup->Fill(detid, edepiter->first, edepiter->second);
0070       }
0071       ntup->Fill(detid, -1, esum);  // fill sum over all layers for each detector
0072     }
0073   }
0074   return 0;
0075 }
0076
0077 int G4EdepNtuple::End(PHCompositeNode * /*topNode*/)
0078 {
0079   outfile->cd();
0080   ntup->Write();
0081   outfile->Write();
0082   outfile->Close();
0083   delete outfile;
0084   hm->dumpHistos(_filename, "UPDATE");
0085   return 0;
0086 }
0087
0088 void G4EdepNtuple::AddNode(const std::string &name, const int detid)
0089 {
0090   _node_postfix.insert(name);
0091   _detid[name] = detid;
0092   return;
0093 }