CaloAna/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
0017 // Cluster includes
0018 #include <calobase/RawCluster.h>
0019 #include <calobase/RawClusterContainer.h>
0020
0021 #include <fun4all/Fun4AllHistoManager.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023
0024 #include <phool/getClass.h>
0025
0026 #include <TFile.h>
0027 #include <TNtuple.h>
0028
0029 #include <cassert>
0030 #include <sstream>
0031 #include <string>
0032
0033 using namespace std;
0034
0035 CaloAna::CaloAna(const std::string& name, const std::string& filename)
0036   : SubsysReco(name)
0037   , detector("HCALIN")
0038   , outfilename(filename)
0039 {
0040 }
0041
0042 CaloAna::~CaloAna()
0043 {
0044   delete hm;
0045   delete g4hitntuple;
0046   delete g4cellntuple;
0047   delete towerntuple;
0048   delete clusterntuple;
0049 }
0050
0051 int CaloAna::Init(PHCompositeNode*)
0052 {
0053   hm = new Fun4AllHistoManager(Name());
0054   // create and register your histos (all types) here
0055   // TH1 *h1 = new TH1F("h1",....)
0056   // hm->registerHisto(h1);
0057   outfile = new TFile(outfilename.c_str(), "RECREATE");
0058   g4hitntuple = new TNtuple("hitntup", "G4Hits", "x0:y0:z0:x1:y1:z1:edep");
0059   g4cellntuple = new TNtuple("cellntup", "G4Cells", "phi:eta:edep");
0060   towerntuple = new TNtuple("towerntup", "Towers", "phi:eta:energy");
0061   clusterntuple = new TNtuple("clusterntup", "Clusters", "phi:z:energy:towers");
0062   return 0;
0063 }
0064
0065 int CaloAna::process_event(PHCompositeNode* topNode)
0066 {
0067   // For the calorimeters we have the following node name convention
0068   // where detector is the calorimeter name (CEMC, HCALIN, HCALOUT)
0069   // this is the order in which they are reconstructed
0070   //  G4HIT_<detector>: G4 Hits
0071   //  G4CELL_<detector>: Cells (combined hits inside a cell - scintillator, eta/phi bin)
0072   //  TOWER_SIM_<detector>: simulated tower (before pedestal and noise)
0073   //  TOWER_RAW_<detector>: Raw Tower (adc/tdc values - from sims or real data)
0074   //  TOWER_CALIB_<detector>: Calibrated towers
0075   //  CLUSTER_<detector>: clusters
0076   process_g4hits(topNode);
0077   process_g4cells(topNode);
0078   process_towers(topNode);
0079   process_clusters(topNode);
0080   return Fun4AllReturnCodes::EVENT_OK;
0081 }
0082
0083 int CaloAna::process_g4hits(PHCompositeNode* topNode)
0084 {
0085   ostringstream nodename;
0086
0087   // loop over the G4Hits
0088   nodename.str("");
0089   nodename << "G4HIT_" << detector;
0090   PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
0091   if (hits)
0092   {
0093     // this returns an iterator to the beginning and the end of our G4Hits
0094     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0095     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0096
0097     {
0098       // the pointer to the G4Hit is hit_iter->second
0099       g4hitntuple->Fill(hit_iter->second->get_x(0),
0100                         hit_iter->second->get_y(0),
0101                         hit_iter->second->get_z(0),
0102                         hit_iter->second->get_x(1),
0103                         hit_iter->second->get_y(1),
0104                         hit_iter->second->get_z(1),
0105                         hit_iter->second->get_edep());
0106     }
0107   }
0108   return Fun4AllReturnCodes::EVENT_OK;
0109 }
0110
0111 int CaloAna::process_g4cells(PHCompositeNode* topNode)
0112 {
0113   ostringstream nodename;
0114
0115   // loop over the G4Hits
0116   nodename.str("");
0117   nodename << "G4CELL_" << detector;
0118   PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, nodename.str());
0119   if (cells)
0120   {
0121     PHG4CellContainer::ConstRange cell_range = cells->getCells();
0122     int phibin = -999;
0123     int etabin = -999;
0124     for (PHG4CellContainer::ConstIterator cell_iter = cell_range.first; cell_iter != cell_range.second; cell_iter++)
0125     {
0126       if (cell_iter->second->has_binning(PHG4CellDefs::scintillatorslatbinning))
0127       {
0128         phibin = PHG4CellDefs::ScintillatorSlatBinning::get_row(cell_iter->second->get_cellid());
0129         etabin = PHG4CellDefs::ScintillatorSlatBinning::get_column(cell_iter->second->get_cellid());
0130       }
0131       else if (cell_iter->second->has_binning(PHG4CellDefs::sizebinning))
0132       {
0133         phibin = PHG4CellDefs::SizeBinning::get_phibin(cell_iter->second->get_cellid());
0134         etabin = PHG4CellDefs::SizeBinning::get_zbin(cell_iter->second->get_cellid());
0135       }
0136       else if (cell_iter->second->has_binning(PHG4CellDefs::spacalbinning))
0137       {
0138         phibin = PHG4CellDefs::SpacalBinning::get_phibin(cell_iter->second->get_cellid());
0139         etabin = PHG4CellDefs::SpacalBinning::get_etabin(cell_iter->second->get_cellid());
0140       }
0141       else
0142       {
0143         cout << "unknown cell binning, implement 0x" << hex << PHG4CellDefs::get_binning(cell_iter->second->get_cellid()) << dec << endl;
0144       }
0145       g4cellntuple->Fill(
0146           phibin,
0147           etabin,
0148           cell_iter->second->get_edep());
0149     }
0150   }
0151   return Fun4AllReturnCodes::EVENT_OK;
0152 }
0153
0154 int CaloAna::process_towers(PHCompositeNode* topNode)
0155 {
0156   ostringstream nodename;
0157   ostringstream geonodename;
0158
0159   // loop over the G4Hits
0160   nodename.str("");
0161   nodename << "TOWER_CALIB_" << detector;
0162   geonodename.str("");
0163   geonodename << "TOWERGEOM_" << detector;
0164   RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, geonodename.str().c_str());
0165   if (!towergeom)
0166   {
0167     return Fun4AllReturnCodes::EVENT_OK;
0168   }
0169   RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, nodename.str().c_str());
0170   if (towers)
0171   {
0172     // again pair of iterators to begin and end of tower map
0173     RawTowerContainer::ConstRange tower_range = towers->getTowers();
0174     for (RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++)
0175
0176     {
0177       int phibin = tower_iter->second->get_binphi();
0178       int etabin = tower_iter->second->get_bineta();
0179       double phi = towergeom->get_phicenter(phibin);
0180       double eta = towergeom->get_etacenter(etabin);
0181       towerntuple->Fill(phi,
0182                         eta,
0183                         tower_iter->second->get_energy());
0184     }
0185   }
0186   return Fun4AllReturnCodes::EVENT_OK;
0187 }
0188
0189 int CaloAna::process_clusters(PHCompositeNode* topNode)
0190 {
0191   ostringstream nodename;
0192
0193   // loop over the G4Hits
0194   nodename.str("");
0195   nodename << "CLUSTER_" << detector;
0196   RawClusterContainer* clusters = findNode::getClass<RawClusterContainer>(topNode, nodename.str().c_str());
0197   if (clusters)
0198   {
0199     RawClusterContainer::ConstRange cluster_range = clusters->getClusters();
0200     for (RawClusterContainer::ConstIterator cluster_iter = cluster_range.first; cluster_iter != cluster_range.second; cluster_iter++)
0201     {
0202       clusterntuple->Fill(cluster_iter->second->get_phi(),
0203                           cluster_iter->second->get_z(),
0204                           cluster_iter->second->get_energy(),
0205                           cluster_iter->second->getNTowers());
0206     }
0207   }
0208   return Fun4AllReturnCodes::EVENT_OK;
0209 }
0210
0211 int CaloAna::End(PHCompositeNode* /*topNode*/)
0212 {
0213   outfile->cd();
0214   g4hitntuple->Write();
0215   g4cellntuple->Write();
0216   towerntuple->Write();
0217   clusterntuple->Write();
0218   outfile->Write();
0219   outfile->Close();
0220   delete outfile;
0221   hm->dumpHistos(outfilename, "UPDATE");
0222   return 0;
0223 }