sEPDChannelAnalysis/src/CaloAna.cc

0001 #include "CaloAna.h"
0002
0003 // G4Hits includes
0004 #include <g4main/PHG4Hit.h>
0005 #include <g4main/PHG4HitContainer.h>
0006
0007 // G4Cells includes
0008 #include <g4detectors/PHG4Cell.h>
0009 #include <g4detectors/PHG4CellContainer.h>
0010
0011 // Tower includes
0012 #include <calobase/RawTower.h>
0013 #include <calobase/RawTowerContainer.h>
0014 #include <calobase/RawTowerGeom.h>
0015 #include <calobase/RawTowerGeomContainer.h>
0016 #include <calobase/TowerInfoContainerv4.h>
0017 #include <calobase/TowerInfov4.h>
0018
0019 // Cluster includes
0020 #include <calobase/RawCluster.h>
0021 #include <calobase/RawClusterContainer.h>
0022
0023 #include <fun4all/Fun4AllHistoManager.h>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025
0026 #include <phool/getClass.h>
0027
0028 #include <TFile.h>
0029 #include <TNtuple.h>
0030 #include <TTree.h>
0031 #include <TH2.h>
0032
0033 #include <Event/Event.h>
0034 #include <Event/packet.h>
0035 #include <cassert>
0036 #include <sstream>
0037 #include <string>
0038
0039 using namespace std;
0040
0041 CaloAna::CaloAna(const std::string& name, const std::string& filename,float zero_sup)
0042   : SubsysReco(name)
0043   , detector("sEPD")
0044   , m_outputFileName(filename)
0045   , min_energy(zero_sup)
0046   , num_channels(768)
0047   , m_event(-1)
0048 {
0049 }
0050
0051 CaloAna::~CaloAna()
0052 {
0053   delete hm;
0054   delete g4hitntuple;
0055   delete g4cellntuple;
0056   delete towerntuple;
0057   delete clusterntuple;
0058 }
0059
0060 int CaloAna::Init(PHCompositeNode*)
0061 {
0062   hm = new Fun4AllHistoManager(Name());
0063   // create and register your histos (all types) here
0064   // TH1 *h1 = new TH1F("h1",....)
0065   // hm->registerHisto(h1);
0066     m_histoFileName = m_outputFileName + "_histos.root";
0067         m_outputFileName = m_outputFileName + ".root";
0068
0069     hm->setOutfileName(m_histoFileName.c_str());
0070
0071   outfile = new TFile(m_outputFileName.c_str(), "RECREATE");
0072
0073
0074   towerntuple = new TTree("towerntup", "Towers");
0075   towerntuple->Branch("energy",&m_energy);
0076   towerntuple->Branch("etabin",&m_etabin);
0077   towerntuple->Branch("phibin",&m_phibin);
0078
0079
0080
0081
0082     for(int i = 0; i < num_channels; i++) {
0083         std::string histo_name = "channel" + std::to_string(i);
0084
0085         auto histo = new TH1F(histo_name.c_str(),histo_name.c_str(),1200,0,1200);
0086
0087         hm->registerHisto(histo);
0088     }
0089
0090   // for (int i = 0; i < 31;i++)
0091   //   {
0092   //     towerntuple->Branch(Form("waveform_%i",i),&m_waveform[i]);
0093   //   }
0094
0095   return 0;
0096 }
0097
0098 int CaloAna::process_event(PHCompositeNode* topNode)
0099 {
0100     ++m_event;
0101   process_towers(topNode);
0102
0103     if(m_event % 100 == 0) {
0104         std::cout << "Event " << m_event << std::endl;
0105     }
0106
0107   return Fun4AllReturnCodes::EVENT_OK;
0108 }
0109
0110
0111 int CaloAna::process_towers(PHCompositeNode* topNode)
0112 {
0113   ostringstream nodename;
0114   ostringstream geonodename;
0115
0116   nodename.str("");
0117   nodename << "TOWERS_SEPD_" << detector;
0118   geonodename.str("");
0119
0120   TowerInfoContainer* offlinetowers = findNode::getClass<TowerInfoContainerv4>(topNode, "TOWERS_SEPD");
0121
0122
0123     if(offlinetowers == NULL){
0124         std::cout << "No towers found!!!" << std::endl;
0125         return -1;
0126     }
0127
0128   int size = offlinetowers->size(); //online towers should be the same!
0129
0130
0131   for (int channel = 0; channel < size;channel++)
0132     {
0133       TowerInfo* offlinetower = offlinetowers->get_tower_at_channel(channel);
0134       float offlineenergy = offlinetower->get_energy();
0135       if(offlineenergy < min_energy){
0136         continue;
0137     }
0138     unsigned int towerkey = offlinetowers->encode_key(channel);
0139       int ieta = offlinetowers->getTowerEtaBin(towerkey);
0140       int iphi = offlinetowers->getTowerPhiBin(towerkey);
0141
0142
0143
0144       m_energy.push_back(offlineenergy);
0145       m_etabin.push_back(ieta);
0146       m_phibin.push_back(iphi);
0147
0148     std::string histo_name = "channel" + std::to_string(channel);
0149
0150     auto h_energy = dynamic_cast<TH1*>(hm->getHisto(histo_name.c_str()));
0151     h_energy->Fill(offlineenergy);
0152
0153       // towerntuple->Fill(offlineenergy,ieta,iphi);
0154     }
0155
0156   // cout << 2 << endl;
0157
0158   // cout << 5 << endl;
0159   towerntuple->Fill();
0160
0161   m_etabin.clear();
0162   m_phibin.clear();
0163   m_energy.clear();
0164
0165   return Fun4AllReturnCodes::EVENT_OK;
0166 }
0167
0168 int CaloAna::End(PHCompositeNode* /*topNode*/)
0169 {
0170     hm->dumpHistos(m_histoFileName.c_str(),"RECREATE");
0171
0172     outfile->cd();
0173
0174   towerntuple->Write();
0175
0176   outfile->Write();
0177     outfile->Close();
0178     delete outfile;
0179
0180   return 0;
0181 }
0182