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