barrel_ecal/scripts/emcal_barrel_pi0_reader.cxx

0001 //////////////////////////
0002 // EMCAL Barrel detector
0003 // Pi0 dataset
0004 // M. Scott 05/2021
0005 //////////////////////////
0006 #include "HepMC3/GenEvent.h"
0007 #include "HepMC3/Print.h"
0008 #include "HepMC3/ReaderAscii.h"
0009 #include "HepMC3/WriterAscii.h"
0010
0011 #include <iostream>
0012
0013 #include "TROOT.h"
0014 #include "TH1F.h"
0015 #include "TH2F.h"
0016 #include "TH3F.h"
0017 #include "TMath.h"
0018 #include "TStyle.h"
0019 #include "TCanvas.h"
0020
0021 using namespace HepMC3;
0022
0023 void emcal_barrel_pi0_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_pi0.hepmc") {
0024   // Setting for graphs
0025   gROOT->SetStyle("Plain");
0026   gStyle->SetOptFit(1);
0027   gStyle->SetLineWidth(1);
0028   gStyle->SetPadTickX(1);
0029   gStyle->SetPadTickY(1);
0030   gStyle->SetPadGridX(1);
0031   gStyle->SetPadGridY(1);
0032   gStyle->SetPadLeftMargin(0.14);
0033   gStyle->SetPadRightMargin(0.17);
0034
0035   ReaderAscii hepmc_input(in_fname);
0036   int events_parsed = 0;
0037   GenEvent evt(Units::GEV, Units::MM);
0038
0039   // Histograms
0040   TH1F* h_pi0_energy = new TH1F("h_pi0_energy", "pi0 energy;E [GeV];Events",         100, -0.5, 30.5);
0041   TH1F* h_pi0_eta    = new TH1F("h_pi0_eta",    "pi0 #eta;#eta;Events",              100, -10.0, 10.0);
0042   TH1F* h_pi0_theta  = new TH1F("h_pi0_theta",  "pi0 #theta;#theta [degree];Events", 100, -0.5, 180.5);
0043   TH1F* h_pi0_phi    = new TH1F("h_pi0_phi",    "pi0 #phi;#phi [degree];Events",     100, -180.5, 180.5);
0044   TH2F* h_pi0_pzpt   = new TH2F("h_pi0_pzpt",   "pi0 pt vs pz;pt [GeV];pz [GeV]",    100, -0.5, 30.5, 100, -30.5, 30.5);
0045   TH2F* h_pi0_pxpy   = new TH2F("h_pi0_pxpy",   "pi0 px vs py;px [GeV];py [GeV]",    100, -30.5, 30.5, 100, -30.5, 30.5);
0046   TH3F* h_pi0_p      = new TH3F("h_pi0_p",      "pi0 p;px [GeV];py [GeV];pz [GeV]",  100, -30.5, 30.5, 100, -30.5, 30.5, 100, -30.5, 30.5);
0047
0048   while (!hepmc_input.failed()) {
0049     // Read event from input file
0050     hepmc_input.read_event(evt);
0051     // If reading failed - exit loop
0052     if (hepmc_input.failed())
0053       break;
0054
0055     for (const auto& v : evt.vertices()) {
0056       for (const auto& p : v->particles_out()) {
0057         if (p->pid() == 11) {
0058           h_pi0_energy->Fill(p->momentum().e());
0059           h_pi0_eta->Fill(p->momentum().eta());
0060           h_pi0_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
0061           h_pi0_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
0062           h_pi0_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
0063           h_pi0_pxpy->Fill(p->momentum().px(), p->momentum().py());
0064           h_pi0_p->Fill(p->momentum().px(), p->momentum().py(), p->momentum().pz());
0065         }
0066       }
0067     }
0068     evt.clear();
0069     events_parsed++;
0070   }
0071   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0072
0073   TCanvas* c = new TCanvas("c", "c", 500, 500);
0074   h_pi0_energy->GetYaxis()->SetTitleOffset(1.8);
0075   h_pi0_energy->SetLineWidth(2);
0076   h_pi0_energy->SetLineColor(kBlue);
0077   h_pi0_energy->DrawClone();
0078   c->SaveAs("results/input_emcal_barrel_pi0_energy.png");
0079   c->SaveAs("results/input_emcal_barrel_pi0_energy.pdf");
0080
0081   TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
0082   h_pi0_eta->GetYaxis()->SetTitleOffset(1.9);
0083   h_pi0_eta->SetLineWidth(2);
0084   h_pi0_eta->SetLineColor(kBlue);
0085   h_pi0_eta->DrawClone();
0086   c1->SaveAs("results/input_emcal_barrel_pi0_eta.png");
0087   c1->SaveAs("results/input_emcal_barrel_pi0_eta.pdf");
0088
0089   TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
0090   h_pi0_theta->GetYaxis()->SetTitleOffset(1.8);
0091   h_pi0_theta->SetLineWidth(2);
0092   h_pi0_theta->SetLineColor(kBlue);
0093   h_pi0_theta->DrawClone();
0094   c2->SaveAs("results/input_emcal_barrel_pi0_theta.png");
0095   c2->SaveAs("results/input_emcal_barrel_pi0_theta.pdf");
0096
0097   TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
0098   h_pi0_phi->GetYaxis()->SetTitleOffset(1.8);
0099   h_pi0_phi->SetLineWidth(2);
0100   h_pi0_phi->GetYaxis()->SetRangeUser(0.0, h_pi0_phi->GetMaximum() + 100.0);
0101   h_pi0_phi->SetLineColor(kBlue);
0102   h_pi0_phi->DrawClone();
0103   c3->SaveAs("results/input_emcal_barrel_pi0_phi.png");
0104   c3->SaveAs("results/input_emcal_barrel_pi0_phi.pdf");
0105
0106   TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
0107   h_pi0_pzpt->GetYaxis()->SetTitleOffset(1.4);
0108   h_pi0_pzpt->SetLineWidth(2);
0109   h_pi0_pzpt->SetLineColor(kBlue);
0110   h_pi0_pzpt->DrawClone("COLZ");
0111   c4->SaveAs("results/input_emcal_barrel_pi0_pzpt.png");
0112   c4->SaveAs("results/input_emcal_barrel_pi0_pzpt.pdf");
0113
0114   TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
0115   h_pi0_pxpy->GetYaxis()->SetTitleOffset(1.4);
0116   h_pi0_pxpy->SetLineWidth(2);
0117   h_pi0_pxpy->SetLineColor(kBlue);
0118   h_pi0_pxpy->DrawClone("COLZ");
0119   c5->SaveAs("results/input_emcal_barrel_pi0_pxpy.png");
0120   c5->SaveAs("results/input_emcal_barrel_pi0_pxpy.pdf");
0121
0122   TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
0123   h_pi0_p->GetYaxis()->SetTitleOffset(1.8);
0124   h_pi0_p->GetXaxis()->SetTitleOffset(1.6);
0125   h_pi0_p->GetZaxis()->SetTitleOffset(1.6);
0126   h_pi0_p->SetLineWidth(2);
0127   h_pi0_p->SetLineColor(kBlue);
0128   h_pi0_p->DrawClone();
0129   c6->SaveAs("results/input_emcal_barrel_pi0_p.png");
0130   c6->SaveAs("results/input_emcal_barrel_pi0_p.pdf");
0131 }