Tracking/HitCountEval/HitCountEval.cc

0001 /*!
0002  *  \file       HitCountEval.cc
0003  *  \brief      Evaluation module for PHG4TrackFastSim output
0004  *  \details     input: PHG4TruthInfoContainer G4TruthInfo
0005  *              PHG4HitContainer G4HIT_EGEM_0
0006  *              PHG4HitContainer G4HIT_EGEM_1
0007  *              PHG4HitContainer G4HIT_EGEM_3
0008  *          PHG4HitContainer G4HIT_FGEM_0
0009  *          PHG4HitContainer G4HIT_FGEM_1
0010  *          PHG4HitContainer G4HIT_FGEM_2
0011  *              PHG4HitContainer G4HIT_FGEM_3
0012  *          PHG4HitContainer G4HIT_FGEM_4
0013  *          PHG4HitContainer G4HIT_MAPS
0014  *          PHG4HitContainer G4HIT_SVTX
0015  *  \author     Giorgian Borca-Tasciuc <giorgian.borca-tasciuc@stonybrook.edu>
0016  */
0017 #include "HitCountEval.h"
0018
0019 #include <phool/getClass.h>
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021 #include <g4main/PHG4HitContainer.h>
0022 #include <g4main/PHG4TruthInfoContainer.h>
0023 #include <g4main/PHG4Particle.h>
0024 #include <g4main/PHG4Hit.h>
0025 #include <fun4all/PHTFileServer.h>
0026
0027 #include <TTree.h>
0028 #include <TH1F.h>
0029 #include <iostream>
0030
0031 const char *const HitCountEval::hit_containers[10] {"G4HIT_EGEM_0", "G4HIT_EGEM_1", "G4HIT_EGEM_3",
0032     "G4HIT_FGEM_0", "G4HIT_FGEM_1", "G4HIT_FGEM_2", "G4HIT_FGEM_3","G4HIT_FGEM_4", "G4HIT_MAPS",
0033     "G4HIT_SVTX"};
0034
0035
0036 #define NELEMS(a) (sizeof(a)/sizeof(a[0]))
0037
0038 HitCountEval::HitCountEval(const std::string &name, const std::string &file_name):
0039     SubsysReco {name},
0040     output_file_name {file_name},
0041     event {0}
0042 {
0043 }
0044
0045
0046 //----------------------------------------------------------------------------//
0047 //-- Init():
0048 //--   Intialize all histograms, trees, and ntuples
0049 //----------------------------------------------------------------------------//
0050 int HitCountEval::Init(PHCompositeNode *topNode) {
0051     std::cout << PHWHERE << " Openning file " << output_file_name << std::endl;
0052     PHTFileServer::get().open(output_file_name, "RECREATE");
0053
0054     // create TTree
0055     for (size_t i {0}; i < NELEMS(hit_containers); ++i) {
0056         hit_count[i] = new TH1F {(std::string(hit_containers[i]) + "-hits").c_str(),
0057             "Hit Count", 100, -4.5, 4.5};
0058         eta_count[i] = new TH1F {(std::string(hit_containers[i]) + "-etas").c_str(),
0059             "Eta Count", 100, -4.5, 4.5};
0060
0061         hits[i] = new TTree {hit_containers[i], ""};
0062         hits[i]->Branch("hit_id", &hit_id, "hit_id/l");
0063         hits[i]->Branch("x", &hit_x, "x/F");
0064         hits[i]->Branch("y", &hit_y, "y/F");
0065         hits[i]->Branch("z", &hit_z, "z/F");
0066         hits[i]->Branch("px", &hit_px, "px/F");
0067         hits[i]->Branch("py", &hit_py, "py/F");
0068         hits[i]->Branch("pz", &hit_pz, "pz/F");
0069         hits[i]->Branch("particle_initial_px", &particle_initial_px, "particle_initial_px/F");
0070         hits[i]->Branch("particle_initial_py", &particle_initial_py, "particle_initial_py/F");
0071         hits[i]->Branch("particle_initial_pz", &particle_initial_pz, "particle_initial_pz/F");
0072         hits[i]->Branch("layer", &hit_layer, "hit_layer/i");
0073         hits[i]->Branch("event", &hit_event, "event/i");
0074
0075         normalized_hits[i] = new TTree {(std::string(hit_containers[i]) + "_normalized").c_str(), ""};
0076         normalized_hits[i]->Branch("eta", &normalized_hit_eta, "eta/D");
0077         normalized_hits[i]->Branch("hit_count", &normalized_hit_count, "eta/D");
0078     }
0079
0080     return Fun4AllReturnCodes::EVENT_OK;
0081 }
0082
0083 //----------------------------------------------------------------------------//
0084 //-- process_event():
0085 //--   Call user instructions for every event.
0086 //--   This function contains the analysis structure.
0087 //----------------------------------------------------------------------------//
0088 int HitCountEval::process_event(PHCompositeNode *const topNode) {
0089     event++;
0090     if (verbosity >= 2 and event % 1000 == 0)
0091         std::cout << PHWHERE << "Events processed: " << event << std::endl;
0092
0093     const auto truth = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0094     const auto particle_range = truth->GetPrimaryParticleRange();
0095     for (auto it = particle_range.first; it != particle_range.second; ++it) {
0096         PHG4Particle const* g4particle{it->second};
0097         if(!g4particle) {
0098             std::cerr << "Missing PHG4Particle\n";
0099             continue;
0100         }
0101
0102         particle_initial_px =  g4particle->get_px();
0103         particle_initial_py =  g4particle->get_py();
0104         particle_initial_pz =  g4particle->get_pz();
0105
0106         const double mom {std::sqrt(particle_initial_px * particle_initial_px
0107                 + particle_initial_py * particle_initial_py
0108                 + particle_initial_pz * particle_initial_pz)};
0109         const double eta {0.5 *std::log((mom + particle_initial_pz) / (mom - particle_initial_pz))};
0110         for (size_t i {0}; i < NELEMS(hit_containers); ++i) {
0111             eta_count[i]->Fill(eta);
0112             const auto h = findNode::getClass<PHG4HitContainer>(topNode, hit_containers[i]);
0113             if (!h) {
0114                 std::cerr << "Couldn't find " << hit_containers[i] << '\n';
0115                 continue;
0116             }
0117
0118             const auto hit_range = h->getHits();
0119             for (auto hp = hit_range.first; hp != hit_range.second; ++hp) {
0120                 PHG4Hit *const hit {hp->second};
0121                 for (int j {0}; j < 2; ++j) {
0122                     hit_count[i]->Fill(eta);
0123                     hit_id = hit->get_hit_id();
0124                     hit_x = hit->get_x(j);
0125                     hit_y = hit->get_y(j);
0126                     hit_z = hit->get_z(j);
0127                     hit_px = hit->get_px(j);
0128                     hit_py = hit->get_py(j);
0129                     hit_pz = hit->get_pz(j);
0130                     hits[i]->Fill();
0131                 }
0132             }
0133         }
0134     }
0135
0136     return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138
0139 int HitCountEval::End(PHCompositeNode *topNode) {
0140     PHTFileServer::get().cd(output_file_name);
0141
0142     for (size_t i {0}; i < NELEMS(hit_containers); ++i) {
0143         const Long64_t nbins {hit_count[i]->GetSize() - 2}; /* -2 for underflow and overflow */
0144         for (Long64_t j {0}; j < nbins; ++j) {
0145             if (eta_count[i]->GetBinContent(j + 1) != 0) {
0146                 const Double_t n {hit_count[i]->GetBinContent(j + 1)};
0147                 const Double_t N {eta_count[i]->GetBinContent(j + 1)};
0148
0149                 normalized_hit_eta = hit_count[i]->GetBinCenter(j + 1);
0150                 normalized_hit_count = n / N;
0151                 normalized_hits[i]->Fill();
0152
0153             }
0154             normalized_hits[i]->Write();
0155             hits[i]->Write();
0156         }
0157     }
0158
0159
0160     return Fun4AllReturnCodes::EVENT_OK;
0161 }
0162
0163 void HitCountEval::set_filename(const char *const file) {
0164     if (file)
0165         output_file_name = file;
0166 }