analysis/Ana/PHAna.cc

0001 #include "PHAna.h"
0002
0003 #include <phool/phool.h>
0004 #include <phool/getClass.h>
0005 #include <g4main/PHG4HitContainer.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4VtxPoint.h>
0010 #include <fun4all/PHTFileServer.h>
0011 #include <fun4all/Fun4AllServer.h>
0012
0013 #include <TNtuple.h>
0014 #include <TH1.h>
0015
0016 #include <iostream>
0017
0018 using namespace std;
0019
0020 //____________________________________
0021 PHAna::PHAna(const string &name):
0022     SubsysReco(name),
0023     _flags( NONE )
0024 {
0025  //initialize
0026  _event = 0;
0027  _outfile = "phana.root";
0028  _truth = 0;
0029  _sf = 0;
0030  _hcalin_towers = 0;
0031 }
0032
0033 //___________________________________
0034 int PHAna::Init(PHCompositeNode *topNode)
0035 {
0036  cout << PHWHERE << " Openning file " << _outfile << endl;
0037  PHTFileServer::get().open( _outfile, "RECREATE");
0038
0039  if(_flags & TRUTH) _truth = new TNtuple("truth","truth","event:Energy:Rapidity:Pt:Phi:PID");
0040  if(_flags & HIST ) create_histos();
0041  if(_flags & SF ) _sf = new TNtuple("sf","sf","event:sfhcalin:sfhcalout:sfcemc");
0042  return 0;
0043 }
0044
0045 //__________________________________
0046 //Call user instructions for every event
0047 int PHAna::process_event(PHCompositeNode *topNode)
0048 {
0049  _event++;
0050  if(_event%1000==0) cout << PHWHERE << "Events processed: " << _event << endl;
0051  GetNodes(topNode);
0052  if(_flags & TRUTH) fill_truth_ntuple(topNode);
0053
0054  if(_flags & HIST) fill_histos(topNode);
0055
0056  if(_flags & SF) fill_sf_ntuple(topNode);
0057
0058  return 0;
0059 }
0060
0061 //__________________________________
0062 void PHAna::fill_histos(PHCompositeNode *topNode)
0063 {
0064  Fun4AllServer *se = Fun4AllServer::instance();
0065  Fun4AllHistoManager *hm = se->getHistoManager("HISTOS");
0066
0067  if(_truth_container)
0068  {
0069   PHG4VtxPoint *point = _truth_container->GetPrimaryVtx( _truth_container->GetPrimaryVertexIndex() );
0070   TH1F *h = (TH1F*) hm->getHisto("zvertex");
0071   if(h) h->Fill( point->get_z() );
0072   h = (TH1F*) hm->getHisto("xvertex");
0073   if(h) h->Fill( point->get_x() );
0074   h = (TH1F*) hm->getHisto("yvertex");
0075   if(h) h->Fill( point->get_y() );
0076  }
0077
0078 }
0079
0080 //_____________________________________
0081 void PHAna::fill_truth_ntuple(PHCompositeNode *topNode)
0082 {
0083  map<int, PHG4Particle*>::const_iterator particle_iter;
0084  float ntvars[7];
0085
0086  PHG4TruthInfoContainer::ConstRange primary_range =
0087      _truth_container->GetPrimaryParticleRange();
0088
0089  for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0090      particle_iter != primary_range.second; ++particle_iter)
0091   {
0092    PHG4Particle *particle = particle_iter->second;
0093    ntvars[0] = _event;
0094    ntvars[1] = particle->get_e();
0095    ntvars[2] = 0.5*log((particle->get_e()+particle->get_pz())/
0096                (particle->get_e()-particle->get_pz()));
0097    ntvars[3] = sqrt(Square(particle->get_px())+Square(particle->get_py()));
0098    ntvars[4] = atan2(particle->get_py(), particle->get_px());
0099    ntvars[5] = particle->get_pid();
0100   _truth->Fill( ntvars );
0101   }
0102
0103 }
0104
0105 //_____________________________________
0106 void PHAna::fill_sf_ntuple(PHCompositeNode *topNode)
0107 {
0108   //Total Outer HCal energy deposited and light yield
0109   float e_hcout = 0.;
0110   float ev_hcout = 0.;
0111   PHG4HitContainer::ConstRange hcalout_hit_range = _hcalout_hit_container->getHits();
0112    for(PHG4HitContainer::ConstIterator hit_iter = hcalout_hit_range.first;
0113        hit_iter != hcalout_hit_range.second; hit_iter++)
0114     {
0115      PHG4Hit *this_hit =  hit_iter->second;
0116      e_hcout += this_hit->get_edep();
0117      ev_hcout += this_hit->get_light_yield();
0118     }
0119
0120   //Total Outer HCal absorbed energy
0121   float ea_hcout = 0.;
0122   PHG4HitContainer::ConstRange hcalout_abs_hit_range = _hcalout_abs_hit_container->getHits();
0123   for (PHG4HitContainer::ConstIterator hit_iter = hcalout_abs_hit_range.first;
0124        hit_iter != hcalout_abs_hit_range.second; hit_iter++)
0125    {
0126     PHG4Hit *this_hit =  hit_iter->second;
0127     ea_hcout += this_hit->get_edep();
0128    }
0129
0130   //Total Inner HCal energy deposit and light yield
0131   float e_hcin =0.;
0132   float ev_hcin = 0.;
0133   PHG4HitContainer::ConstRange hcalin_hit_range = _hcalin_hit_container->getHits();
0134   for (PHG4HitContainer::ConstIterator hit_iter = hcalin_hit_range.first;
0135        hit_iter != hcalin_hit_range.second; hit_iter++)
0136    {
0137     PHG4Hit *this_hit =  hit_iter->second;
0138     e_hcin += this_hit->get_edep();
0139     ev_hcin += this_hit->get_light_yield();
0140    }
0141
0142   //Total Inner HCal absorbed energy
0143   float ea_hcin = 0.;
0144   PHG4HitContainer::ConstRange hcalin_abs_hit_range = _hcalin_abs_hit_container->getHits();
0145   for (PHG4HitContainer::ConstIterator hit_iter = hcalin_abs_hit_range.first;
0146        hit_iter != hcalin_abs_hit_range.second; hit_iter++)
0147    {
0148     PHG4Hit *this_hit =  hit_iter->second;
0149
0150     ea_hcin += this_hit->get_edep();
0151    }
0152
0153   PHG4HitContainer::ConstRange hcalin_spt_hit_range = _hcalin_spt_hit_container->getHits();
0154   for (PHG4HitContainer::ConstIterator hit_iter = hcalin_spt_hit_range.first;
0155        hit_iter != hcalin_spt_hit_range.second; hit_iter++)
0156    {
0157     PHG4Hit *this_hit =  hit_iter->second;
0158     ea_hcin += this_hit->get_edep();
0159    }
0160
0161   //Total EMCal deposited energy and light yield
0162   float e_cemc = 0.;
0163   float ev_cemc = 0.;
0164   PHG4HitContainer::ConstRange cemc_hit_range = _cemc_hit_container->getHits();
0165   for (PHG4HitContainer::ConstIterator hit_iter = cemc_hit_range.first;
0166        hit_iter != cemc_hit_range.second; hit_iter++)
0167    {
0168     PHG4Hit *this_hit =  hit_iter->second;
0169
0170     e_cemc += this_hit->get_edep();
0171     ev_cemc += this_hit->get_light_yield();
0172    }
0173
0174    //Total EMCal absorbed energy
0175    float ea_cemc = 0.;
0176    PHG4HitContainer::ConstRange cemc_abs_hit_range = _cemc_abs_hit_container->getHits();
0177   for (PHG4HitContainer::ConstIterator hit_iter = cemc_abs_hit_range.first;
0178        hit_iter != cemc_abs_hit_range.second; hit_iter++)
0179    {
0180     PHG4Hit *this_hit =  hit_iter->second;
0181     ea_cemc += this_hit->get_edep();
0182    }
0183
0184   PHG4HitContainer::ConstRange cemc_electronics_hit_range = _cemc_electronics_hit_container->getHits();
0185   for (PHG4HitContainer::ConstIterator hit_iter = cemc_electronics_hit_range.first;
0186        hit_iter != cemc_electronics_hit_range.second; hit_iter++)
0187    {
0188     PHG4Hit *this_hit =  hit_iter->second;
0189
0190     ea_cemc += this_hit->get_edep();
0191    }
0192
0193
0194   float ntvars[4];
0195   ntvars[0] = _event;
0196   ntvars[1] = ev_hcin/(e_hcin+ea_hcin);
0197   ntvars[2] = ev_hcout/(e_hcout+ea_hcout);
0198   ntvars[3] = ev_cemc/(e_cemc+ea_cemc);
0199
0200   _sf->Fill( ntvars );
0201 }
0202
0203 //___________________________________
0204 void PHAna::create_histos()
0205 {
0206  Fun4AllServer *se = Fun4AllServer::instance();
0207  Fun4AllHistoManager *hm = se->getHistoManager("HISTOS");
0208  if (!hm)
0209  {
0210   hm = new Fun4AllHistoManager("HISTOS");
0211   se->registerHistoManager(hm);
0212  }
0213  hm->registerHisto( new TH1F("zvertex","zvertex",120,-30,30) );
0214  hm->registerHisto( new TH1F("xvertex","xvertex",120,-5,5) );
0215  hm->registerHisto( new TH1F("yvertex","yvertex",120,-5,5) );
0216
0217 }
0218
0219 //___________________________________
0220 void PHAna::GetNodes(PHCompositeNode *topNode)
0221 {
0222  //DST objects
0223  //Truth container
0224   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0225   if(!_truth_container && _event<2) cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << endl;
0226
0227  //Outer HCal hit container
0228   _hcalout_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALOUT");
0229   if(!_hcalout_hit_container && _event<2) cout << PHWHERE << " G4HIT_HCALOUT node not found on node tree" << endl;
0230
0231   //Outer HCal absorber hit container
0232     _hcalout_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_HCALOUT");
0233    if(!_hcalout_abs_hit_container && _event<2) cout << PHWHERE << " G4HIT_ABSORBER_HCALOUT node not found on node tree" << endl;
0234
0235  //Inner HCal hit container
0236    _hcalin_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALIN");
0237     if(!_hcalin_hit_container && _event<2) cout << PHWHERE << " G4HIT_HCALIN node not found on node tree" << endl;
0238
0239  //Inner HCal absorber hit container
0240   _hcalin_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_HCALIN");
0241    if(!_hcalin_abs_hit_container && _event<2) cout << PHWHERE << " G4HIT_ABSORBER_HCALIN node not found on node tree" << endl;
0242
0243   //Inner HCal support structure hit container
0244      _hcalin_spt_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALIN_SPT");
0245    if(!_hcalin_spt_hit_container && _event<2) cout << PHWHERE << " G4HIT_HCALIN_SPT node not found on node tree" << endl;
0246
0247   //EMCAL hit container
0248     _cemc_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_CEMC");
0249    if(!_cemc_hit_container && _event<2) cout << PHWHERE << " G4HIT_CEMC node not found on node tree" << endl;
0250
0251   //EMCal absorber hit container
0252      _cemc_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_CEMC");
0253    if(!_cemc_abs_hit_container && _event<2) cout << PHWHERE << " G4HIT_ABSORBER_CEMC node not found on node tree" << endl;
0254
0255   //EMCal Electronics hit container
0256    _cemc_electronics_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_CEMC_ELECTRONICS");
0257    if(!_cemc_electronics_hit_container && _event<2) cout << PHWHERE << " G4HIT_CEMC_EMCELECTRONICS node not found on node tree" << endl;
0258
0259 }
0260
0261 //______________________________________
0262 int PHAna::End(PHCompositeNode *topNode)
0263 {
0264  PHTFileServer::get().cd( _outfile );
0265  if( _truth ) _truth->Write();
0266  if( _sf ) _sf->Write();
0267  if(_flags & HIST)
0268  {
0269  Fun4AllServer *se = Fun4AllServer::instance();
0270  Fun4AllHistoManager *hm = se->getHistoManager("HISTOS");
0271  for(unsigned int i=0; i<hm->nHistos(); i++) hm->getHisto(i)->Write();
0272  }
0273  return 0;
0274 }
0275