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   hm = new Fun4AllHistoManager(Name());
0037   outfile = new TFile(_filename.c_str(), "RECREATE");
0038   ntup = new TNtuple("snglntup", "G4Sngls", "phi:theta:e:detid:layer:x0:y0:z0:x1:y1:z1:edep");
0039   ntup_e = new TNtuple("sngl_e", "G4SnglEdeps", "phi:theta:e:detid:layer:edep");
0040   //  ntup->SetDirectory(0);
0041   TH1 *h1 = new TH1F("edep1GeV", "edep 0-1GeV", 1000, 0, 1);
0042   eloss.push_back(h1);
0043   h1 = new TH1F("edep100GeV", "edep 0-100GeV", 1000, 0, 100);
0044   eloss.push_back(h1);
0045   return 0;
0046 }
0047
0048 int G4SnglNtuple::process_event(PHCompositeNode *topNode)
0049 {
0050   // get the primary particle which did this to us....
0051   PHG4TruthInfoContainer *truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0052   double px = std::numeric_limits<double>::quiet_NaN();
0053   double py = std::numeric_limits<double>::quiet_NaN();
0054   double pz = std::numeric_limits<double>::quiet_NaN();
0055   double e = std::numeric_limits<double>::quiet_NaN();
0056   double pt = std::numeric_limits<double>::quiet_NaN();
0057   double phi = std::numeric_limits<double>::quiet_NaN();
0058   double theta = std::numeric_limits<double>::quiet_NaN();
0059   if (truthInfoList)
0060   {
0061     const PHG4TruthInfoContainer::Range primRange = truthInfoList->GetPrimaryParticleRange();
0062     px = primRange.first->second->get_px();
0063     py = primRange.first->second->get_py();
0064     pz = primRange.first->second->get_pz();
0065     e = primRange.first->second->get_e();
0066     pt = sqrt(px * px + py * py);
0067     phi = atan2(py, px);
0068     theta = atan2(pt, pz);
0069   }
0070   std::ostringstream nodename;
0071   std::set<std::string>::const_iterator iter;
0072   std::vector<TH1 *>::const_iterator eiter;
0073
0074   std::map<int, double> layer_edep_map;
0075   std::map<int, double>::const_iterator edepiter;
0076
0077   for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
0078   {
0079     layer_edep_map.clear();
0080     int detid = (_detid.find(*iter))->second;
0081     nodename.str("");
0082     nodename << "G4HIT_" << *iter;
0083     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
0084     if (hits)
0085     {
0086       double esum = 0;
0087       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0088       for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0089       {
0090         layer_edep_map[hit_iter->second->get_layer()] += hit_iter->second->get_edep();
0091         esum += hit_iter->second->get_edep();
0092         ntup->Fill(phi, theta, e, detid,
0093                    hit_iter->second->get_detid(),
0094                    hit_iter->second->get_x(0),
0095                    hit_iter->second->get_y(0),
0096                    hit_iter->second->get_z(0),
0097                    hit_iter->second->get_x(1),
0098                    hit_iter->second->get_y(1),
0099                    hit_iter->second->get_z(1),
0100                    hit_iter->second->get_edep());
0101       }
0102       for (edepiter = layer_edep_map.begin(); edepiter != layer_edep_map.end(); ++edepiter)
0103       {
0104         ntup_e->Fill(phi, theta, e, detid, edepiter->first, edepiter->second);
0105       }
0106       ntup_e->Fill(phi, theta, e, detid, -1, esum);  // fill sum over all layers for each detector
0107
0108       for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
0109       {
0110         (*eiter)->Fill(esum);
0111       }
0112     }
0113   }
0114   return 0;
0115 }
0116
0117 int G4SnglNtuple::End(PHCompositeNode * /*topNode*/)
0118 {
0119   outfile->cd();
0120   ntup->Write();
0121   outfile->Write();
0122   outfile->Close();
0123   delete outfile;
0124   hm->dumpHistos(_filename, "UPDATE");
0125   return 0;
0126 }
0127
0128 void G4SnglNtuple::AddNode(const std::string &name, const int detid)
0129 {
0130   _node_postfix.insert(name);
0131   _detid[name] = detid;
0132   return;
0133 }