g4simulation/g4histos/G4HitNtuple.cc

0001 #include "G4HitNtuple.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 <TH1.h>
0013 #include <TNtuple.h>
0014
0015 #include <sstream>
0016 #include <utility>  // for pair
0017
0018 using namespace std;
0019
0020 G4HitNtuple::G4HitNtuple(const std::string &name, const std::string &filename)
0021   : SubsysReco(name)
0022   , nblocks(0)
0023   , hm(nullptr)
0024   , _filename(filename)
0025   , ntup(nullptr)
0026   , outfile(nullptr)
0027 {
0028 }
0029
0030 G4HitNtuple::~G4HitNtuple()
0031 {
0032   //  delete ntup;
0033   delete hm;
0034 }
0035
0036 int G4HitNtuple::Init(PHCompositeNode * /*unused*/)
0037 {
0038   hm = new Fun4AllHistoManager(Name());
0039   outfile = new TFile(_filename.c_str(), "RECREATE");
0040   ntup = new TNtuple("hitntup", "G4Hits", "detid:lyr:slat:x0:y0:z0:x1:y1:z1:edep");
0041   //  ntup->SetDirectory(0);
0042   TH1 *h1 = new TH1F("edep1GeV", "edep 0-1GeV", 1000, 0, 1);
0043   eloss.push_back(h1);
0044   h1 = new TH1F("edep100GeV", "edep 0-100GeV", 1000, 0, 100);
0045   eloss.push_back(h1);
0046   return 0;
0047 }
0048
0049 int G4HitNtuple::process_event(PHCompositeNode *topNode)
0050 {
0051   ostringstream nodename;
0052   set<string>::const_iterator iter;
0053   vector<TH1 *>::const_iterator eiter;
0054   for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
0055   {
0056     int detid = (_detid.find(*iter))->second;
0057     nodename.str("");
0058     nodename << "G4HIT_" << *iter;
0059     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
0060     if (hits)
0061     {
0062       double esum = 0;
0063       //          double numhits = hits->size();
0064       //          nhits[i]->Fill(numhits);
0065       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0066       for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0067
0068       {
0069         esum += hit_iter->second->get_edep();
0070         ntup->Fill(detid,
0071                    hit_iter->second->get_layer(),
0072                    hit_iter->second->get_scint_id(),
0073                    hit_iter->second->get_x(0),
0074                    hit_iter->second->get_y(0),
0075                    hit_iter->second->get_z(0),
0076                    hit_iter->second->get_x(1),
0077                    hit_iter->second->get_y(1),
0078                    hit_iter->second->get_z(1),
0079                    hit_iter->second->get_edep());
0080       }
0081       for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
0082       {
0083         (*eiter)->Fill(esum);
0084       }
0085     }
0086   }
0087   return 0;
0088 }
0089
0090 int G4HitNtuple::End(PHCompositeNode * /*topNode*/)
0091 {
0092   outfile->cd();
0093   ntup->Write();
0094   outfile->Write();
0095   outfile->Close();
0096   delete outfile;
0097   hm->dumpHistos(_filename, "UPDATE");
0098   return 0;
0099 }
0100
0101 void G4HitNtuple::AddNode(const std::string &name, const int detid)
0102 {
0103   _node_postfix.insert(name);
0104   _detid[name] = detid;
0105   return;
0106 }