barrel_ecal/scripts/emcal_barrel_particles_reader.cxx

0001 //////////////////////////
0002 // EMCAL Barrel detector
0003 // Electron dataset
0004 // J.KIM 04/02/2021
0005 // M.Zurek 05/05/2021
0006 //////////////////////////
0007 #include "HepMC3/GenEvent.h"
0008 #include "HepMC3/Print.h"
0009 #include "HepMC3/ReaderAscii.h"
0010 #include "HepMC3/WriterAscii.h"
0011
0012 #include "TROOT.h"
0013 #include "TH1.h"
0014 #include "TH2.h"
0015 #include "TH3.h"
0016 #include "TStyle.h"
0017 #include "TCanvas.h"
0018 #include "TMath.h"
0019
0020 #include <iostream>
0021 #include <fmt/core.h>
0022
0023 #include "emcal_barrel_common_functions.h"
0024
0025 using namespace HepMC3;
0026
0027 void save_canvas(TCanvas* c, std::string label)
0028 {
0029   c->SaveAs(fmt::format("results/{}.png",label).c_str());
0030   c->SaveAs(fmt::format("results/{}.pdf",label).c_str());
0031 }
0032
0033 void save_canvas(TCanvas* c, std::string label, std::string particle_label)
0034 {
0035   std::string label_with_E = fmt::format("input_emcal_barrel_{}_{}", particle_label, label);
0036   save_canvas(c, label_with_E);
0037 }
0038
0039 void emcal_barrel_particles_reader(std::string particle_name = "electron") {
0040
0041   // Setting for graphs
0042   gROOT->SetStyle("Plain");
0043   gStyle->SetOptFit(1);
0044   gStyle->SetLineWidth(1);
0045   gStyle->SetPadTickX(1);
0046   gStyle->SetPadTickY(1);
0047   gStyle->SetPadGridX(1);
0048   gStyle->SetPadGridY(1);
0049   gStyle->SetPadLeftMargin(0.14);
0050   gStyle->SetPadRightMargin(0.17);
0051
0052   std::string in_fname = fmt::format("./data/emcal_barrel_{}.hepmc",particle_name);
0053   ReaderAscii hepmc_input(in_fname);
0054   int events_parsed = 0;
0055   GenEvent evt(Units::GEV, Units::MM);
0056
0057   // Histograms
0058   TH1F* h_energy = new TH1F(fmt::format("h_{}_energy",particle_name).c_str(), fmt::format("{} energy;E [GeV];Events",particle_name).c_str(),         100, -0.5, 30.5);
0059   TH1F* h_eta    = new TH1F(fmt::format("h_{}_eta",particle_name).c_str(),    fmt::format("{} #eta;#eta;Events",particle_name).c_str(),              100, -10.0, 10.0);
0060   TH1F* h_theta  = new TH1F(fmt::format("h_{}_theta",particle_name).c_str(),  fmt::format("{} #theta;#theta [degree];Events",particle_name).c_str(), 100, -0.5, 180.5);
0061   TH1F* h_phi    = new TH1F(fmt::format("h_{}_phi",particle_name).c_str(),    fmt::format("{} #phi;#phi [degree];Events",particle_name).c_str(),     100, -180.5, 180.5);
0062   TH2F* h_pzpt   = new TH2F(fmt::format("h_{}_pzpt",particle_name).c_str(),  fmt::format("{} pt vs pz;pt [GeV];pz [GeV]",particle_name).c_str(),    100, -0.5, 30.5, 100, -30.5, 30.5);
0063   TH2F* h_pxpy   = new TH2F(fmt::format("h_{}_pxpy",particle_name).c_str(),  fmt::format("{} px vs py;px [GeV];py [GeV]",particle_name).c_str(),    100, -30.5, 30.5, 100, -30.5, 30.5);
0064   TH3F* h_p      = new TH3F(fmt::format("h_{}_p",particle_name).c_str(),     fmt::format("{} p;px [GeV];py [GeV];pz [GeV]",particle_name).c_str(),  100, -30.5, 30.5, 100, -30.5, 30.5, 100, -30.5, 30.5);
0065
0066   while (!hepmc_input.failed()) {
0067     // Read event from input file
0068     hepmc_input.read_event(evt);
0069     // If reading failed - exit loop
0070     if (hepmc_input.failed())
0071       break;
0072
0073     auto [id, mass] = extract_particle_parameters(particle_name);
0074
0075     for (const auto& v : evt.vertices()) {
0076       for (const auto& p : v->particles_out()) {
0077         if (p->pid() == id) {
0078           h_energy->Fill(p->momentum().e());
0079           h_eta->Fill(p->momentum().eta());
0080           h_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
0081           h_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
0082           h_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
0083           h_pxpy->Fill(p->momentum().px(), p->momentum().py());
0084           h_p->Fill(p->momentum().px(), p->momentum().py(), p->momentum().pz());
0085         }
0086       }
0087     }
0088     evt.clear();
0089     events_parsed++;
0090   }
0091   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0092
0093   TCanvas* c = new TCanvas("c", "c", 500, 500);
0094   h_energy->GetYaxis()->SetTitleOffset(1.8);
0095   h_energy->SetLineWidth(2);
0096   h_energy->SetLineColor(kBlue);
0097   h_energy->DrawClone();
0098   save_canvas(c, "energy", particle_name);
0099
0100   TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
0101   h_eta->GetYaxis()->SetTitleOffset(1.9);
0102   h_eta->SetLineWidth(2);
0103   h_eta->SetLineColor(kBlue);
0104   h_eta->DrawClone();
0105   save_canvas(c1, "eta", particle_name);
0106
0107   TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
0108   h_theta->GetYaxis()->SetTitleOffset(1.8);
0109   h_theta->SetLineWidth(2);
0110   h_theta->SetLineColor(kBlue);
0111   h_theta->DrawClone();
0112   save_canvas(c2, "theta", particle_name);
0113
0114   TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
0115   h_phi->GetYaxis()->SetTitleOffset(1.8);
0116   h_phi->SetLineWidth(2);
0117   h_phi->GetYaxis()->SetRangeUser(0.0, h_phi->GetMaximum() + 100.0);
0118   h_phi->SetLineColor(kBlue);
0119   h_phi->DrawClone();
0120   save_canvas(c3, "phi", particle_name);
0121
0122   TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
0123   h_pzpt->GetYaxis()->SetTitleOffset(1.4);
0124   h_pzpt->SetLineWidth(2);
0125   h_pzpt->SetLineColor(kBlue);
0126   h_pzpt->DrawClone("COLZ");
0127   save_canvas(c4, "pzpt", particle_name);
0128
0129   TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
0130   h_pxpy->GetYaxis()->SetTitleOffset(1.4);
0131   h_pxpy->SetLineWidth(2);
0132   h_pxpy->SetLineColor(kBlue);
0133   h_pxpy->DrawClone("COLZ");
0134   save_canvas(c5, "pxpy", particle_name);
0135
0136   TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
0137   h_p->GetYaxis()->SetTitleOffset(1.8);
0138   h_p->GetXaxis()->SetTitleOffset(1.6);
0139   h_p->GetZaxis()->SetTitleOffset(1.6);
0140   h_p->SetLineWidth(2);
0141   h_p->SetLineColor(kBlue);
0142   h_p->DrawClone();
0143   save_canvas(c6, "p", particle_name);
0144 }
0145