AntiSigma/src/TimingNtuple.cc

0001 #include "TimingNtuple.h"
0002
0003 #include <g4main/PHG4HitContainer.h>
0004 #include <g4main/PHG4Hit.h>
0005
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4VtxPoint.h>
0009
0010 #include <fun4all/Fun4AllHistoManager.h>
0011
0012 #include <phool/getClass.h>
0013
0014 #include <TFile.h>
0015 #include <TH1.h>
0016 #include <TH2.h>
0017 #include <TNtuple.h>
0018
0019 #include <boost/foreach.hpp>
0020
0021 #include<sstream>
0022
0023 using namespace std;
0024
0025 TimingNtuple::TimingNtuple(const std::string &name, const std::string &filename):
0026   SubsysReco( name ),
0027   nblocks(0),
0028   hm(NULL),
0029   _filename(filename),
0030   ntup(NULL),
0031   outfile(NULL)
0032 {}
0033
0034 TimingNtuple::~TimingNtuple()
0035 {
0036   //  delete ntup;
0037   delete hm;
0038 }
0039
0040
0041 int
0042 TimingNtuple::Init( PHCompositeNode* )
0043 {
0044   ostringstream hname, htit;
0045   hm = new  Fun4AllHistoManager(Name());
0046   outfile = new TFile(_filename.c_str(), "RECREATE");
0047   //ntup = new TNtuple("hitntup", "G4Hits", "detid:layer:x0:y0:z0:x1:y1:z1:edep");
0048   ntupt = new TNtuple("tnt", "G4Hits", "ID:t:dtotal");
0049 //  ntupavet = new TNtuple("time", "G4Hits", "ID:avt");
0050  //  ntup->SetDirectory(0);
0051   TH1 *h1 = new TH1F("edep1GeV","edep 0-1GeV",1000,0,1);
0052   eloss.push_back(h1);
0053   h1 = new TH1F("edep100GeV","edep 0-100GeV",1000,0,100);
0054   eloss.push_back(h1);
0055   return 0;
0056 }
0057
0058 int
0059 TimingNtuple::process_event( PHCompositeNode* topNode )
0060 {
0061   ostringstream nodename;
0062
0063   map<int, PHG4Particle*>::const_iterator particle_iter;
0064   PHG4TruthInfoContainer *_truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0065
0066   PHG4TruthInfoContainer::ConstRange primary_range =
0067     _truth_container->GetPrimaryParticleRange();
0068   float ntvars[5] = {0};
0069   for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0070        particle_iter != primary_range.second; ++particle_iter)
0071   {
0072     PHG4Particle *particle = particle_iter->second;
0073     ntvars[0] = atan2(particle->get_py(), particle->get_px());
0074     ntvars[1] = atan(sqrt(particle->get_py()*particle->get_py() +
0075               particle->get_px()*particle->get_px()) /
0076              particle->get_pz());
0077     if (ntvars[1] < 0)
0078       {
0079     ntvars[1]+=M_PI;
0080       }
0081     ntvars[2] = 0.5*log((particle->get_e()+particle->get_pz())/
0082             (particle->get_e()-particle->get_pz()));
0083     ntvars[3] = particle->get_e();
0084     ntvars[4] = particle->get_pz();
0085   }
0086
0087   set<string>::const_iterator iter;
0088   vector<TH1 *>::const_iterator eiter;
0089
0090   PHG4HitContainer *hits = nullptr;
0091
0092   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_CEMC");
0093
0094   if (hits)
0095     {
0096       //          double numhits = hits->size();
0097       //          nhits[i]->Fill(numhits);
0098       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0099       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0100
0101     {
0102       double dtotal = get_dtotal(hit_iter->second, ntvars[0],ntvars[1]);
0103       ntupt->Fill(0,hit_iter->second->get_avg_t(),dtotal);
0104     }
0105     }
0106
0107
0108   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_HCALIN");
0109
0110   if (hits)
0111     {
0112       //          double numhits = hits->size();
0113       //          nhits[i]->Fill(numhits);
0114       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0115       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0116
0117     {
0118       double dtotal = get_dtotal(hit_iter->second, ntvars[0],ntvars[1]);
0119       ntupt->Fill(1,hit_iter->second->get_avg_t(),dtotal);
0120     }
0121     }
0122
0123   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_HCALOUT");
0124
0125   if (hits)
0126     {
0127       //          double numhits = hits->size();
0128       //          nhits[i]->Fill(numhits);
0129       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0130       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0131
0132     {
0133       double dtotal = get_dtotal(hit_iter->second, ntvars[0],ntvars[1]);
0134       ntupt->Fill(2,hit_iter->second->get_avg_t(),dtotal);
0135     }
0136     }
0137
0138   return 0;
0139 }
0140
0141 int
0142 TimingNtuple::End(PHCompositeNode * topNode)
0143 {
0144   outfile->cd();
0145   //  ntup->Write();
0146   ntupt->Write();
0147   outfile->Write();
0148   outfile->Close();
0149   delete outfile;
0150   hm->dumpHistos(_filename, "UPDATE");
0151   return 0;
0152 }
0153
0154 void
0155 TimingNtuple::AddNode(const std::string &name, const int detid)
0156 {
0157  _node_postfix.insert(name);
0158  _detid[name] = detid;
0159  return;
0160 }
0161
0162 double
0163 TimingNtuple::get_dtotal(const PHG4Hit *hit, const double phi, const double theta)
0164 {
0165   double phi_hit =  atan2(hit->get_avg_y(),hit->get_avg_x());
0166   double theta_hit = atan(sqrt(hit->get_avg_x()*hit->get_avg_x()+
0167                    hit->get_avg_y()*hit->get_avg_y()) /
0168               hit->get_avg_z());
0169   double diffphi =  phi_hit - phi;
0170   if (diffphi > M_PI)
0171   {
0172     diffphi -= 2*M_PI;
0173   }
0174   else if (diffphi < - M_PI)
0175   {
0176     diffphi += 2*M_PI;
0177   }
0178   if (theta_hit < 0)
0179   {
0180     theta_hit += M_PI;
0181   }
0182   double difftheta = theta_hit-theta;
0183   double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0184   return deltasqrt;
0185 }