g4simulation/g4histos/G4TowerNtuple.cc

0001 #include "G4TowerNtuple.h"
0002
0003 #include <calobase/RawTowerGeomContainer.h>
0004 #include <calobase/TowerInfo.h>
0005 #include <calobase/TowerInfoContainerv1.h>
0006
0007 #include <fun4all/Fun4AllHistoManager.h>
0008 #include <fun4all/SubsysReco.h>  // for SubsysReco
0009
0010 #include <phool/getClass.h>
0011
0012 #include <TFile.h>
0013 #include <TH1.h>
0014 #include <TNtuple.h>
0015
0016 #include <cmath>     // for isfinite
0017 #include <iostream>  // for operator<<, basic_ostream
0018 #include <utility>   // for pair, make_pair
0019
0020 G4TowerNtuple::G4TowerNtuple(const std::string &name, const std::string &filename)
0021   : SubsysReco(name)
0022   , m_filename(filename)
0023 {
0024 }
0025
0026 G4TowerNtuple::~G4TowerNtuple()
0027 {
0028   delete hm;
0029 }
0030
0031 int G4TowerNtuple::Init(PHCompositeNode * /*unused*/)
0032 {
0033   delete hm; // make cppcheck happy
0034   hm = new Fun4AllHistoManager(Name());
0035   outfile = new TFile(m_filename.c_str(), "RECREATE");
0036   ntup = new TNtuple("towerntup", "G4Towers", "detid:phi:eta:energy");
0037   //  ntup->SetDirectory(0);
0038   TH1 *h1 = new TH1F("energy1GeV", "energy 0-1GeV", 1000, 0, 1);
0039   eloss.push_back(h1);
0040   h1 = new TH1F("energy100GeV", "energy 0-100GeV", 1000, 0, 100);
0041   eloss.push_back(h1);
0042   return 0;
0043 }
0044
0045 int G4TowerNtuple::process_event(PHCompositeNode *topNode)
0046 {
0047   std::string nodename;
0048   std::string geonodename;
0049   for (const auto &iter : m_node_postfix)
0050   {
0051     int detid = (m_detid.find(iter))->second;
0052     nodename = "TOWERINFO_" + m_tower_type[iter];
0053     geonodename = "TOWERGEOM_" + iter;
0054     RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, geonodename);
0055     if (!towergeom)
0056     {
0057       std::cout << "no geometry node " << geonodename << " for " << iter << std::endl;
0058       continue;
0059     }
0060
0061     TowerInfoContainer *towers = findNode::getClass<TowerInfoContainer>(topNode, nodename);
0062     if (towers)
0063     {
0064       double esum = 0;
0065       unsigned int nchannels = towers->size();
0066       for (unsigned int channel = 0; channel < nchannels; channel++)
0067       {
0068         TowerInfo *tower = towers->get_tower_at_channel(channel);
0069         double energy = tower->get_energy();
0070         if (!std::isfinite(energy))
0071         {
0072           std::cout << "invalid energy: " << energy << std::endl;
0073         }
0074         esum += energy;
0075         unsigned int towerkey = towers->encode_key(channel);
0076         int etabin = towers->getTowerEtaBin(towerkey);
0077         int phibin = towers->getTowerPhiBin(towerkey);
0078
0079         // to search the map fewer times, cache the geom object until the layer changes
0080         double phi = towergeom->get_phicenter(phibin);
0081         double eta = towergeom->get_etacenter(etabin);
0082         ntup->Fill(detid,
0083                    phi,
0084                    eta,
0085                    energy);
0086       }
0087       for (auto *eiter : eloss)
0088       {
0089         eiter->Fill(esum);
0090       }
0091     }
0092   }
0093   return 0;
0094 }
0095
0096 int G4TowerNtuple::End(PHCompositeNode * /*topNode*/)
0097 {
0098   outfile->cd();
0099   ntup->Write();
0100   outfile->Write();
0101   outfile->Close();
0102   delete outfile;
0103   hm->dumpHistos(m_filename, "UPDATE");
0104   return 0;
0105 }
0106
0107 void G4TowerNtuple::AddNode(const std::string &name, const std::string &twrtype, const int detid)
0108 {
0109   std::string twrname;
0110   twrname = twrtype + "_" + name;
0111   m_node_postfix.insert(name);
0112   m_tower_type.insert(make_pair(name, twrname));
0113   m_detid[name] = detid;
0114   return;
0115 }