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