g4simulation/g4histos/G4HitTTree.cc

0001 #include "G4HitTTree.h"
0002
0003 #include "G4RootHitContainer.h"
0004
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <fun4all/SubsysReco.h>  // for SubsysReco
0010
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/PHIODataNode.h>    // for PHIODataNode
0013 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0014 #include <phool/PHObject.h>        // for PHObject
0015 #include <phool/getClass.h>
0016
0017 #include <TH1.h>
0018 #include <TH2.h>
0019 #include <TSystem.h>
0020
0021 #include <iostream>  // for operator<<, endl, basic_ost...
0022 #include <map>       // for _Rb_tree_const_iterator
0023 #include <utility>   // for pair
0024
0025 using namespace std;
0026
0027 G4HitTTree::G4HitTTree(const std::string &name)
0028   : SubsysReco(name)
0029   , savehits(1)
0030   , evtno(0)
0031   , hm(nullptr)
0032   , etot_hist(nullptr)
0033   , eion_etot_hist(nullptr)
0034 {
0035   BlackHoleName("BH_1");  // initialize this to what we have in our common sims
0036 }
0037
0038 int G4HitTTree::Init(PHCompositeNode *topNode)
0039 {
0040   if (!_detector.size())
0041   {
0042     cout << "Detector not set via Detector(<name>) method" << endl;
0043     cout << "(it is the name appended to the G4HIT_<name> nodename)" << endl;
0044     cout << "you do not want to run like this, exiting now" << endl;
0045     gSystem->Exit(1);
0046   }
0047   hm = new Fun4AllHistoManager("HITHIST");
0048   etot_hist = new TH1F("etot", "total deposited energy", 200, 0, 20);
0049   hm->registerHisto(etot_hist);
0050   eion_etot_hist = new TH2F("eion_etot", "ionization energy vs total energy", 200, 0, 20, 200, 0, 20);
0051   hm->registerHisto(eion_etot_hist);
0052
0053   PHNodeIterator iter(topNode);
0054   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0055   G4RootHitContainer *hits = new G4RootHitContainer();
0056   PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(hits, _outnodename, "PHObject");
0057   dstNode->addNode(node);
0058   evtno = 0;
0059   return 0;
0060 }
0061
0062 int G4HitTTree::process_event(PHCompositeNode *topNode)
0063 {
0064   evtno++;
0065   G4RootHitContainer *hits = findNode::getClass<G4RootHitContainer>(topNode, _outnodename);
0066   PHG4HitContainer *g4hits = findNode::getClass<PHG4HitContainer>(topNode, _hitnodename);
0067   double etot = 0;
0068   double eion = 0;
0069   if (g4hits)
0070   {
0071     PHG4HitContainer::ConstRange hit_range = g4hits->getHits();
0072     // shower_z->Reset();
0073     //   cout << "Number of Hits: " << g4hits->size() << endl;
0074     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0075     {
0076       PHG4Hit *inhit = hit_iter->second;
0077       if (savehits)
0078       {
0079         hits->AddHit(inhit);
0080       }
0081       etot += inhit->get_edep();
0082       eion += inhit->get_eion();
0083     }
0084     etot_hist->Fill(etot);
0085     eion_etot_hist->Fill(etot, eion);
0086   }
0087   g4hits = findNode::getClass<PHG4HitContainer>(topNode, _absorbernodename);
0088   if (g4hits)
0089   {
0090     PHG4HitContainer::ConstRange hit_range = g4hits->getHits();
0091     // shower_z->Reset();
0092     //   cout << "Number of Hits: " << g4hits->size() << endl;
0093     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0094     {
0095       PHG4Hit *inhit = hit_iter->second;
0096       if (savehits)
0097       {
0098         hits->AddHit(inhit);
0099       }
0100     }
0101   }
0102   double eleak = 0;
0103   g4hits = findNode::getClass<PHG4HitContainer>(topNode, _blackholenodename);
0104   if (g4hits)
0105   {
0106     PHG4HitContainer::ConstRange hit_range = g4hits->getHits();
0107     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0108     {
0109       if (savehits)
0110       {
0111         //        PHG4Hit *g4h = hits->AddHit( *(hit_iter->second));
0112       }
0113       eleak += hit_iter->second->get_edep();
0114     }
0115   }
0116   hits->set_etotal(etot);
0117   hits->set_eion(eion);
0118   hits->set_leakage(eleak);
0119   hits->set_event(evtno);
0120   return 0;
0121 }
0122
0123 int G4HitTTree::End(PHCompositeNode * /*topNode*/)
0124 {
0125   hm->dumpHistos("HitHistos.root");
0126   delete hm;
0127   return 0;
0128 }
0129
0130 void G4HitTTree::Detector(const std::string &det)
0131 {
0132   _detector = det;
0133   _outnodename = "G4RootHit_" + det;
0134   _hitnodename = "G4HIT_" + det;
0135   _absorbernodename = "G4HIT_ABSORBER_" + det;
0136 }
0137
0138 void G4HitTTree::BlackHoleName(const std::string &bh)
0139 {
0140   _blackholenodename = "G4HIT_" + bh;
0141 }