nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx

0001 #include <cmath>
0002 #include <fstream>
0003 #include <iostream>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "ROOT/RDataFrame.hxx"
0008 #include <TH1D.h>
0009 #include <TH2D.h>
0010 #include <TRandom3.h>
0011 #include <TFile.h>
0012 #include <TChain.h>
0013 #include <TTree.h>
0014 #include <TMath.h>
0015 #include <TVector3.h>
0016 #include <TCanvas.h>
0017 #include "TROOT.h"
0018 #include "TRandom.h"
0019 #include "TH3.h"
0020 #include "TLorentzVector.h"
0021 #include "TStyle.h"
0022
0023 #include "DD4hep/Detector.h"
0024 #include "DDRec/CellIDPositionConverter.h"
0025
0026 #include <podio/Frame.h>
0027 #include <podio/CollectionBase.h>
0028 #include "podio/ROOTReader.h"
0029 #include "podio/CollectionIDTable.h"
0030 #include "podio/ObjectID.h"
0031
0032 #include "edm4hep/MCParticleCollection.h"
0033 #include "edm4hep/MCParticle.h"
0034
0035 #include "edm4hep/SimCalorimeterHitCollection.h"
0036 #include "edm4hep/SimCalorimeterHit.h"
0037
0038 #include "edm4hep/CalorimeterHitCollection.h"
0039 #include "edm4hep/CalorimeterHit.h"
0040
0041 #include "edm4eic/ClusterCollection.h"
0042 #include "edm4eic/Cluster.h"
0043
0044 #include "edm4eic/CalorimeterHitCollection.h"
0045 #include "edm4eic/CalorimeterHit.h"
0046
0047 #include <edm4eic/vector_utils_legacy.h>
0048
0049 using namespace std;
0050 using namespace edm4hep;
0051
0052 dd4hep::Detector* det = NULL;
0053
0054 constexpr double ETA_MIN = -4.16, ETA_MAX = -1.16;
0055
0056 inline bool inNHCal(double eta) {return (eta >= ETA_MIN && eta <= ETA_MAX);}
0057
0058 inline string addPrefixAfterSlash(const string& path,
0059                                    const string& prefix) {
0060     const auto slash = path.find_last_of('/');
0061     const auto dot = path.find_last_of('.');
0062
0063     string extension = (dot != string::npos) ? path.substr(dot) : "";
0064
0065     if (slash == string::npos)
0066         return prefix + extension;
0067
0068     return path.substr(0, slash + 1) + prefix + extension;
0069 }
0070
0071 int sampling_fraction_analysis(const string &filename, string outname_pdf, string outname_png, TString compact_file)
0072 {
0073
0074     podio::ROOTReader *reader = new podio::ROOTReader();
0075     reader->openFile(filename);
0076     unsigned nEvents = reader->getEntries("events");
0077     cout << "Number of events: " << nEvents << endl;
0078
0079     det = &(dd4hep::Detector::getInstance());
0080     det->fromCompact(compact_file.Data());
0081     det->volumeManager();
0082     det->apply("DD4hepVolumeManager", 0, 0);
0083
0084     dd4hep::rec::CellIDPositionConverter cellid_converter(*det);
0085     auto idSpec = det->readout("HcalEndcapNHits").idSpec();
0086     auto decoder = idSpec.decoder();
0087     const int slice_index = decoder->index("slice");
0088     if (slice_index < 0) {
0089         cerr << "ERROR: 'slice' field not found in cell ID spec!" << endl;
0090         return 1;
0091     }
0092
0093     constexpr int NBINS = 80;
0094
0095     constexpr double E_MIN_GEV = 0.0;
0096     constexpr double E_MAX_GEV = 12.0;
0097     constexpr double SAMP_F_MIN = 0.0;
0098     constexpr double SAMP_F_MAX = 1.0;
0099     constexpr double SAMP_F_LOW = 0.2;
0100
0101     gStyle->SetTitleSize(0.04, "XYZ");
0102     gStyle->SetLabelSize(0.04, "XYZ");
0103     gStyle->SetPadLeftMargin(0.20);
0104     gStyle->SetPadRightMargin(0.15);
0105     gStyle->SetPadBottomMargin(0.15);
0106     gStyle->SetPadTopMargin(0.10);
0107
0108     TH2D *h_sampF_e = new TH2D("h_sampF_e", "nHCal sampling fraction vs. energy (e-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]; counts",
0109                                             NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_MAX);
0110
0111     TH2D *h_sampF_n = new TH2D("h_sampF_n", "nHCal sampling fraction vs. energy (neutron); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]; counts",
0112                                             NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_MAX);
0113
0114     TH2D *h_sampF_pi = new TH2D("h_sampF_pi", "nHCal sampling fraction vs. energy (#pi-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]; counts",
0115                                             NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_MAX);
0116
0117     TH2D *h_sampF_e_Ekin = new TH2D("h_sampF_e_Ekin", "nHCal sampling fraction vs. energy kin (e-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{kin}]; counts",
0118                                             NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_LOW);
0119
0120     TH2D *h_sampF_n_Ekin = new TH2D("h_sampF_n_Ekin", "nHCal sampling fraction vs. energy kin (neutron); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{kin}]; counts",
0121                                             NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_LOW);
0122
0123     TH2D *h_sampF_pi_Ekin = new TH2D("h_sampF_pi_Ekin", "nHCal sampling fraction vs. energy kin (#pi-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{kin}]; counts",
0124                                             NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_LOW);
0125
0126     for (unsigned ev = 0; ev < nEvents; ev++)
0127     {
0128         double hit_Esum = 0;
0129         double hit_scint_Esum = 0;
0130         double singlePart_Ekin = 0;
0131
0132         auto frameData = reader->readNextEntry(podio::Category::Event);
0133         if (!frameData)
0134         {
0135             cerr << "Invalid FrameData at event " << ev << endl;
0136             continue;
0137         }
0138
0139         podio::Frame frame(std::move(frameData));
0140
0141         const edm4hep::MCParticleCollection& MCParticles_coll  = frame.get<edm4hep::MCParticleCollection>("MCParticles");
0142         const edm4hep::SimCalorimeterHitCollection& SimCalorimeterHit_coll = frame.get<edm4hep::SimCalorimeterHitCollection>("HcalEndcapNHits");
0143
0144         if (!SimCalorimeterHit_coll.isValid())
0145         {
0146             cerr << "HcalEndcapNHits collection is not valid!" << endl;
0147         }
0148         if (!MCParticles_coll.isValid())
0149         {
0150             cerr << "MCParticles collection is not valid!" << endl;
0151         }
0152
0153         edm4hep::MCParticle mcpart =  MCParticles_coll.at(0);
0154         TLorentzVector v(mcpart.getMomentum().x, mcpart.getMomentum().y, mcpart.getMomentum().z, mcpart.getEnergy());
0155         if(!inNHCal(v.Eta())) continue;
0156         singlePart_Ekin = mcpart.getEnergy() - mcpart.getMass();
0157         int pdg = mcpart.getPDG();
0158
0159         for (const auto& hit : SimCalorimeterHit_coll)
0160         {
0161                 const int hit_slice = decoder->get(hit.getCellID(), slice_index);
0162                 hit_Esum += hit.getEnergy();
0163                 if(hit_slice == 3) hit_scint_Esum += hit.getEnergy();
0164         }
0165
0166         if (pdg == 11)               // e-
0167         {
0168             h_sampF_e->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
0169             h_sampF_e_Ekin->Fill(singlePart_Ekin, hit_scint_Esum/singlePart_Ekin);
0170         }
0171         else if (pdg == -211)        // pi-
0172         {
0173             h_sampF_pi->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
0174             h_sampF_pi_Ekin->Fill(singlePart_Ekin, hit_scint_Esum/singlePart_Ekin);
0175         }
0176         else if (pdg == 2112)        // neutron
0177         {
0178             h_sampF_n->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum);
0179             h_sampF_n_Ekin->Fill(singlePart_Ekin, hit_scint_Esum/singlePart_Ekin);
0180         }
0181     }
0182
0183     delete reader;
0184
0185     h_sampF_e->Sumw2();
0186     h_sampF_e_Ekin->Sumw2();
0187     h_sampF_pi->Sumw2();
0188     h_sampF_pi_Ekin->Sumw2();
0189     h_sampF_n->Sumw2();
0190     h_sampF_n_Ekin->Sumw2();
0191
0192     TProfile* p_sampF_e = h_sampF_e->ProfileX();
0193     p_sampF_e ->SetTitle("nHCal sampling fraction vs. energy (e-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]");
0194     TProfile* p_sampF_e_Ekin = h_sampF_e_Ekin->ProfileX();
0195     p_sampF_e_Ekin ->SetTitle("nHCal sampling fraction vs. energy kin (e-); E_{kin} [GeV]; f");
0196     TProfile* p_sampF_pi = h_sampF_pi->ProfileX();
0197     p_sampF_pi ->SetTitle("nHCal sampling fraction vs. energy (#pi-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]");
0198     TProfile* p_sampF_pi_Ekin = h_sampF_pi_Ekin->ProfileX();
0199     p_sampF_pi_Ekin ->SetTitle("nHCal sampling fraction vs. energy kin (#pi-); E_{kin} [GeV]; f");
0200     TProfile* p_sampF_n = h_sampF_n->ProfileX();
0201     p_sampF_n ->SetTitle("nHCal sampling fraction vs. energy (neutron); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]");
0202     TProfile* p_sampF_n_Ekin = h_sampF_n_Ekin->ProfileX();
0203     p_sampF_n_Ekin ->SetTitle("nHCal sampling fraction vs. energy kin (neutron); E_{kin} [GeV]; f");
0204
0205     TH1D *h_e = p_sampF_e->ProjectionX("h_e");
0206     TH1D *h_pi = p_sampF_pi->ProjectionX("h_pi");
0207
0208     TH1D *h_e_over_pi = (TH1D*)h_e->Clone("h_e_over_pi");
0209     h_e_over_pi->SetTitle("e/h ratio;E_{kin} [GeV];e/h");
0210     h_e_over_pi->Divide(h_e, h_pi, 1, 1);
0211
0212
0213     TCanvas *c_h = new TCanvas("canvas_h", "c_h", 1600, 800);
0214     c_h->Divide(2,2);
0215     c_h->cd(1);
0216     h_sampF_e->Draw("COLZ");
0217
0218     c_h->cd(2);
0219     h_sampF_pi->Draw("COLZ");
0220
0221     c_h->cd(3);
0222     h_sampF_n->Draw("COLZ");
0223
0224     c_h->SaveAs(outname_pdf.c_str());
0225     c_h->SaveAs(outname_png.c_str());
0226
0227     TCanvas *c_p = new TCanvas("canvas_p", "c_p", 1600, 800);
0228     c_p->Divide(2,2);
0229     c_p->cd(1);
0230     p_sampF_e->SetLineWidth(3); p_sampF_e->SetLineColor(kRed); p_sampF_e->SetMarkerColor(kRed);
0231     p_sampF_e->Draw();
0232     c_p->cd(2);
0233     p_sampF_pi->SetLineWidth(3); p_sampF_pi->SetLineColor(kRed); p_sampF_pi->SetMarkerColor(kRed);
0234     p_sampF_pi->Draw();
0235     c_p->cd(3);
0236     p_sampF_n->SetLineWidth(3); p_sampF_n->SetLineColor(kRed); p_sampF_n->SetMarkerColor(kRed);
0237     p_sampF_n->Draw();
0238     c_p->cd(4);
0239     h_e_over_pi->Draw();
0240
0241     c_p->SaveAs(addPrefixAfterSlash(outname_png, "prof_sampf_vs_Ehit").c_str());
0242     c_p->SaveAs(addPrefixAfterSlash(outname_pdf, "prof_sampf_vs_Ehit").c_str());
0243
0244     TH1D *h_e_Ekin = p_sampF_e_Ekin->ProjectionX("h_e_Ekin");
0245     TH1D *h_pi_Ekin = p_sampF_pi_Ekin->ProjectionX("h_pi_Ekin");
0246
0247     TH1D *h_e_over_pi_Ekin = (TH1D*)h_e_Ekin->Clone("h_e_over_pi_Ekin");
0248     h_e_over_pi_Ekin->SetTitle("e/#pi ratio;E_{kin} [GeV];e/#pi");
0249     h_e_over_pi_Ekin->Divide(h_e_Ekin, h_pi_Ekin, 1, 1);
0250
0251     TCanvas *c_h_Ekin = new TCanvas("canvas_h_Ekin", "c_h_Ekin", 1600, 800);
0252     c_h_Ekin->Divide(2,2);
0253     c_h_Ekin->cd(1);
0254     h_sampF_e_Ekin->Draw("COLZ");
0255
0256     c_h_Ekin->cd(2);
0257     h_sampF_pi_Ekin->Draw("COLZ");
0258
0259     c_h_Ekin->cd(3);
0260     h_sampF_n_Ekin->Draw("COLZ");
0261
0262     c_h_Ekin->SaveAs(addPrefixAfterSlash(outname_png, "hist_sampf_vs_Ekin").c_str());
0263     c_h_Ekin->SaveAs(addPrefixAfterSlash(outname_pdf, "hist_sampf_vs_Ekin").c_str());
0264
0265     TCanvas *c_p_Ekin = new TCanvas("canvas_p_Ekin", "c_p_Ekin", 1600, 800);
0266     c_p_Ekin->Divide(2,2);
0267     c_p_Ekin->cd(1);
0268     p_sampF_e_Ekin->SetLineWidth(3); p_sampF_e_Ekin->SetLineColor(kRed); p_sampF_e_Ekin->SetMarkerColor(kRed);
0269     p_sampF_e_Ekin->Draw();
0270     c_p_Ekin->cd(2);
0271     p_sampF_pi_Ekin->SetLineWidth(3); p_sampF_pi_Ekin->SetLineColor(kRed); p_sampF_pi_Ekin->SetMarkerColor(kRed);
0272     p_sampF_pi_Ekin->Draw();
0273     c_p_Ekin->cd(3);
0274     p_sampF_n_Ekin->SetLineWidth(3); p_sampF_n_Ekin->SetLineColor(kRed); p_sampF_n_Ekin->SetMarkerColor(kRed);
0275     p_sampF_n_Ekin->Draw();
0276     c_p_Ekin->cd(4);
0277     h_e_over_pi_Ekin->Draw();
0278
0279     c_p_Ekin->SaveAs(addPrefixAfterSlash(outname_png, "prof_sampf_vs_Ekin").c_str());
0280     c_p_Ekin->SaveAs(addPrefixAfterSlash(outname_pdf, "prof_sampf_vs_Ekin").c_str());
0281
0282     return 0;
0283 }
0284
0285 int main(int argc, char** argv) {
0286     if (argc < 5) {
0287         cerr << "Usage: " << argv[0] << " <input.root> <out.pdf> <out.png> <compact.xml>" << endl;
0288         return 1;
0289     }
0290     return sampling_fraction_analysis(argv[1], argv[2], argv[3], argv[4]);
0291 }