AntiSigma/src/EdepNtuple.cc

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