File indexing completed on 2024-11-15 08:59:22
0001
0002
0003
0004
0005
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 <fstream>
0022 #include <cstdlib>
0023 #include <fmt/core.h>
0024
0025 #include "hcal_barrel_common_functions.h"
0026
0027 using namespace HepMC3;
0028
0029 void save_canvas(TCanvas* c, std::string label)
0030 {
0031 c->SaveAs(fmt::format("results/{}.png",label).c_str());
0032 c->SaveAs(fmt::format("results/{}.pdf",label).c_str());
0033 }
0034
0035 void save_canvas(TCanvas* c, std::string label, std::string particle_label)
0036 {
0037 std::string label_with_E = fmt::format("input_hcal_barrel_{}_{}", particle_label, label);
0038 save_canvas(c, label_with_E);
0039 }
0040
0041 void hcal_barrel_particles_reader(std::string particle_name = "electron")
0042 {
0043
0044 gROOT->SetStyle("Plain");
0045 gStyle->SetOptFit(1);
0046 gStyle->SetLineWidth(1);
0047 gStyle->SetPadTickX(1);
0048 gStyle->SetPadTickY(1);
0049 gStyle->SetPadGridX(1);
0050 gStyle->SetPadGridY(1);
0051 gStyle->SetPadLeftMargin(0.14);
0052 gStyle->SetPadRightMargin(0.17);
0053
0054 std::string in_fname;
0055 auto env_fname = getenv("JUGGLER_GEN_FILE");
0056 if (env_fname != nullptr)
0057 in_fname = env_fname;
0058 else
0059 in_fname = fmt::format("./data/hcal_barrel_{}.hepmc", particle_name);
0060 ReaderAscii hepmc_input(in_fname);
0061 int events_parsed = 0;
0062 GenEvent evt(Units::GEV, Units::MM);
0063
0064
0065 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);
0066 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);
0067 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);
0068 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);
0069 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);
0070 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);
0071 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);
0072
0073 while (!hepmc_input.failed()) {
0074
0075 hepmc_input.read_event(evt);
0076
0077 if (hepmc_input.failed())
0078 break;
0079
0080 auto [id, mass] = extract_particle_parameters(particle_name);
0081
0082 for (const auto& v : evt.vertices()) {
0083 for (const auto& p : v->particles_out()) {
0084 if (p->pid() == id) {
0085 h_energy->Fill(p->momentum().e());
0086 h_eta->Fill(p->momentum().eta());
0087 h_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
0088 h_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
0089 h_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
0090 h_pxpy->Fill(p->momentum().px(), p->momentum().py());
0091 h_p->Fill(p->momentum().px(), p->momentum().py(), p->momentum().pz());
0092 }
0093 }
0094 }
0095 evt.clear();
0096 events_parsed++;
0097 }
0098 std::cout << "Events parsed and written: " << events_parsed << std::endl;
0099
0100 TCanvas* c = new TCanvas("c", "c", 500, 500);
0101 h_energy->GetYaxis()->SetTitleOffset(1.8);
0102 h_energy->SetLineWidth(2);
0103 h_energy->SetLineColor(kBlue);
0104 h_energy->DrawClone();
0105 save_canvas(c, "energy", particle_name);
0106
0107 TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
0108 h_eta->GetYaxis()->SetTitleOffset(1.9);
0109 h_eta->SetLineWidth(2);
0110 h_eta->SetLineColor(kBlue);
0111 h_eta->DrawClone();
0112 save_canvas(c1, "eta", particle_name);
0113
0114 TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
0115 h_theta->GetYaxis()->SetTitleOffset(1.8);
0116 h_theta->SetLineWidth(2);
0117 h_theta->SetLineColor(kBlue);
0118 h_theta->DrawClone();
0119 save_canvas(c2, "theta", particle_name);
0120
0121 TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
0122 h_phi->GetYaxis()->SetTitleOffset(1.8);
0123 h_phi->SetLineWidth(2);
0124 h_phi->GetYaxis()->SetRangeUser(0.0, h_phi->GetMaximum() + 100.0);
0125 h_phi->SetLineColor(kBlue);
0126 h_phi->DrawClone();
0127 save_canvas(c3, "phi", particle_name);
0128
0129 TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
0130 h_pzpt->GetYaxis()->SetTitleOffset(1.4);
0131 h_pzpt->SetLineWidth(2);
0132 h_pzpt->SetLineColor(kBlue);
0133 h_pzpt->DrawClone("COLZ");
0134 save_canvas(c4, "pzpt", particle_name);
0135
0136 TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
0137 h_pxpy->GetYaxis()->SetTitleOffset(1.4);
0138 h_pxpy->SetLineWidth(2);
0139 h_pxpy->SetLineColor(kBlue);
0140 h_pxpy->DrawClone("COLZ");
0141 save_canvas(c5, "pxpy", particle_name);
0142
0143 TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
0144 h_p->GetYaxis()->SetTitleOffset(1.8);
0145 h_p->GetXaxis()->SetTitleOffset(1.6);
0146 h_p->GetZaxis()->SetTitleOffset(1.6);
0147 h_p->SetLineWidth(2);
0148 h_p->SetLineColor(kBlue);
0149 h_p->DrawClone();
0150 save_canvas(c6, "p", particle_name);
0151 }
0152
