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 <sstream>
0019 #include <utility>  // for pair, make_pair
0020
0021 using namespace std;
0022
0023 G4TowerNtuple::G4TowerNtuple(const std::string &name, const std::string &filename)
0024   : SubsysReco(name)
0025   , nblocks(0)
0026   , hm(nullptr)
0027   , _filename(filename)
0028   , ntup(nullptr)
0029   , outfile(nullptr)
0030 {
0031 }
0032
0033 G4TowerNtuple::~G4TowerNtuple()
0034 {
0035   //  delete ntup;
0036   delete hm;
0037 }
0038
0039 int G4TowerNtuple::Init(PHCompositeNode * /*unused*/)
0040 {
0041   hm = new Fun4AllHistoManager(Name());
0042   outfile = new TFile(_filename.c_str(), "RECREATE");
0043   ntup = new TNtuple("towerntup", "G4Towers", "detid:phi:eta:energy");
0044   //  ntup->SetDirectory(0);
0045   TH1 *h1 = new TH1F("energy1GeV", "energy 0-1GeV", 1000, 0, 1);
0046   eloss.push_back(h1);
0047   h1 = new TH1F("energy100GeV", "energy 0-100GeV", 1000, 0, 100);
0048   eloss.push_back(h1);
0049   return 0;
0050 }
0051
0052 int G4TowerNtuple::process_event(PHCompositeNode *topNode)
0053 {
0054   ostringstream nodename;
0055   ostringstream geonodename;
0056   set<string>::const_iterator iter;
0057   vector<TH1 *>::const_iterator eiter;
0058   for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
0059   {
0060     int detid = (_detid.find(*iter))->second;
0061     nodename.str("");
0062     nodename << "TOWERINFO_" << _tower_type[*iter];
0063     geonodename.str("");
0064     geonodename << "TOWERGEOM_" << *iter;
0065     RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, geonodename.str());
0066     if (!towergeom)
0067     {
0068       cout << "no geometry node " << geonodename.str() << " for " << *iter << endl;
0069       continue;
0070     }
0071
0072     TowerInfoContainer *towers = findNode::getClass<TowerInfoContainerv1>(topNode, nodename.str());
0073     if (towers)
0074     {
0075       double esum = 0;
0076       unsigned int nchannels = towers->size();
0077       for (unsigned int channel = 0; channel < nchannels; channel++)
0078       {
0079         TowerInfo *tower = towers->get_tower_at_channel(channel);
0080         double energy = tower->get_energy();
0081         if (!isfinite(energy))
0082         {
0083           cout << "invalid energy: " << energy << endl;
0084         }
0085         esum += energy;
0086         unsigned int towerkey = towers->encode_key(channel);
0087         int etabin = towers->getTowerEtaBin(towerkey);
0088         int phibin = towers->getTowerPhiBin(towerkey);
0089
0090         // to search the map fewer times, cache the geom object until the layer changes
0091         double phi = towergeom->get_phicenter(phibin);
0092         double eta = towergeom->get_etacenter(etabin);
0093         ntup->Fill(detid,
0094                    phi,
0095                    eta,
0096                    energy);
0097       }
0098       for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
0099       {
0100         (*eiter)->Fill(esum);
0101       }
0102     }
0103   }
0104   return 0;
0105 }
0106
0107 int G4TowerNtuple::End(PHCompositeNode * /*topNode*/)
0108 {
0109   outfile->cd();
0110   ntup->Write();
0111   outfile->Write();
0112   outfile->Close();
0113   delete outfile;
0114   hm->dumpHistos(_filename, "UPDATE");
0115   return 0;
0116 }
0117
0118 void G4TowerNtuple::AddNode(const std::string &name, const std::string &twrtype, const int detid)
0119 {
0120   ostringstream twrname;
0121   twrname << twrtype << "_" << name;
0122   _node_postfix.insert(name);
0123   _tower_type.insert(make_pair(name, twrname.str()));
0124   _detid[name] = detid;
0125   return;
0126 }