g4simulation/g4histos/G4SnglNtuple.cc

0001 #include "G4SnglNtuple.h"
0002
0003 #include <g4main/PHG4Hit.h>
0004 #include <g4main/PHG4HitContainer.h>
0005 #include <g4main/PHG4Particle.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <fun4all/SubsysReco.h>  // for SubsysReco
0010
0011 #include <phool/getClass.h>
0012
0013 #include <TFile.h>
0014 #include <TH1.h>
0015 #include <TNtuple.h>
0016
0017 #include <cmath>  // for atan2, sqrt
0018 #include <limits>
0019 #include <sstream>
0020 #include <utility>  // for pair
0021
0022 G4SnglNtuple::G4SnglNtuple(const std::string &name, const std::string &filename)
0023   : SubsysReco(name)
0024   , _filename(filename)
0025 {
0026 }
0027
0028 G4SnglNtuple::~G4SnglNtuple()
0029 {
0030   //  delete ntup;
0031   delete hm;
0032 }
0033
0034 int G4SnglNtuple::Init(PHCompositeNode * /*unused*/)
0035 {
0036   delete hm; // make cppcheck happy
0037   hm = new Fun4AllHistoManager(Name());
0038   outfile = new TFile(_filename.c_str(), "RECREATE");
0039   ntup = new TNtuple("snglntup", "G4Sngls", "phi:theta:e:detid:layer:x0:y0:z0:x1:y1:z1:edep");
0040   ntup_e = new TNtuple("sngl_e", "G4SnglEdeps", "phi:theta:e:detid:layer: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 G4SnglNtuple::process_event(PHCompositeNode *topNode)
0050 {
0051   // get the primary particle which did this to us....
0052   PHG4TruthInfoContainer *truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0053   double px = std::numeric_limits<double>::quiet_NaN();
0054   double py = std::numeric_limits<double>::quiet_NaN();
0055   double pz = std::numeric_limits<double>::quiet_NaN();
0056   double e = std::numeric_limits<double>::quiet_NaN();
0057   double pt = std::numeric_limits<double>::quiet_NaN();
0058   double phi = std::numeric_limits<double>::quiet_NaN();
0059   double theta = std::numeric_limits<double>::quiet_NaN();
0060   if (truthInfoList)
0061   {
0062     const PHG4TruthInfoContainer::Range primRange = truthInfoList->GetPrimaryParticleRange();
0063     px = primRange.first->second->get_px();
0064     py = primRange.first->second->get_py();
0065     pz = primRange.first->second->get_pz();
0066     e = primRange.first->second->get_e();
0067     pt = sqrt(px * px + py * py);
0068     phi = atan2(py, px);
0069     theta = atan2(pt, pz);
0070   }
0071   std::ostringstream nodename;
0072   std::set<std::string>::const_iterator iter;
0073   std::vector<TH1 *>::const_iterator eiter;
0074
0075   std::map<int, double> layer_edep_map;
0076   std::map<int, double>::const_iterator edepiter;
0077
0078   for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
0079   {
0080     layer_edep_map.clear();
0081     int detid = (_detid.find(*iter))->second;
0082     nodename.str("");
0083     nodename << "G4HIT_" << *iter;
0084     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
0085     if (hits)
0086     {
0087       double esum = 0;
0088       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0089       for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0090       {
0091         layer_edep_map[hit_iter->second->get_layer()] += hit_iter->second->get_edep();
0092         esum += hit_iter->second->get_edep();
0093         ntup->Fill(phi, theta, e, detid,
0094                    hit_iter->second->get_detid(),
0095                    hit_iter->second->get_x(0),
0096                    hit_iter->second->get_y(0),
0097                    hit_iter->second->get_z(0),
0098                    hit_iter->second->get_x(1),
0099                    hit_iter->second->get_y(1),
0100                    hit_iter->second->get_z(1),
0101                    hit_iter->second->get_edep());
0102       }
0103       for (edepiter = layer_edep_map.begin(); edepiter != layer_edep_map.end(); ++edepiter)
0104       {
0105         ntup_e->Fill(phi, theta, e, detid, edepiter->first, edepiter->second);
0106       }
0107       ntup_e->Fill(phi, theta, e, detid, -1, esum);  // fill sum over all layers for each detector
0108
0109       for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
0110       {
0111         (*eiter)->Fill(esum);
0112       }
0113     }
0114   }
0115   return 0;
0116 }
0117
0118 int G4SnglNtuple::End(PHCompositeNode * /*topNode*/)
0119 {
0120   outfile->cd();
0121   ntup->Write();
0122   outfile->Write();
0123   outfile->Close();
0124   delete outfile;
0125   hm->dumpHistos(_filename, "UPDATE");
0126   return 0;
0127 }
0128
0129 void G4SnglNtuple::AddNode(const std::string &name, const int detid)
0130 {
0131   _node_postfix.insert(name);
0132   _detid[name] = detid;
0133   return;
0134 }