zdc/scripts/analysis_zdc_particles.cxx

0001 ////////////////////////////////////////
0002 // Read reconstruction ROOT output file
0003 // Plot variables
0004 ////////////////////////////////////////
0005
0006 #include "ROOT/RDataFrame.hxx"
0007 #include <iostream>
0008
0009 #include "edm4hep/MCParticleCollection.h"
0010 #include "edm4hep/SimCalorimeterHitCollection.h"
0011 #include "edm4eic/CalorimeterHitCollection.h"
0012 #include "edm4eic/CalorimeterHitData.h"
0013
0014 #include "TCanvas.h"
0015 #include "TStyle.h"
0016 #include "TMath.h"
0017 #include "TH1.h"
0018 #include "TF1.h"
0019 #include "TH1D.h"
0020
0021 using ROOT::RDataFrame;
0022 using namespace ROOT::VecOps;
0023
0024 void analysis_zdc_particles(
0025   const std::string& input_fname = "sim_output/sim_zdc_uniform_neutrons.edm4hep.root",
0026   const std::string& results_path = "results/far_forward/zdc/"
0027 ) {
0028   // Setting for graphs
0029   gROOT->SetStyle("Plain");
0030   gStyle->SetOptFit(1);
0031   gStyle->SetLineWidth(2);
0032   gStyle->SetPadTickX(1);
0033   gStyle->SetPadTickY(1);
0034   gStyle->SetPadGridX(1);
0035   gStyle->SetPadGridY(1);
0036   gStyle->SetPadLeftMargin(0.14);
0037   gStyle->SetPadRightMargin(0.14);
0038
0039   ROOT::EnableImplicitMT();
0040   ROOT::RDataFrame d0("events", input_fname);
0041
0042   // Thrown Energy [GeV]
0043   auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0044     auto p = input[2];
0045     auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
0046     return energy;
0047   };
0048
0049   // Number of hits
0050   auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
0051
0052   // Energy deposition [GeV]
0053   auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0054     auto total_edep = 0.0;
0055     for (const auto& i: evt)
0056       total_edep += i.energy;
0057     return total_edep;
0058   };
0059
0060   // Sampling fraction = Esampling / Ethrown
0061   auto fsam = [](const double sampled, const double thrown) {
0062     return sampled / thrown;
0063   };
0064
0065   // Hit position X
0066   auto hit_x_position = [&] (const std::vector<edm4hep::SimCalorimeterHitData>& hits) {
0067     std::vector<double> result;
0068     for(const auto& h: hits)
0069         result.push_back(h.position.x); //mm
0070   return result;
0071   };
0072
0073   // Hit position Y
0074   auto hit_y_position = [&] (const std::vector<edm4hep::SimCalorimeterHitData>& hits) {
0075         std::vector<double> result;
0076         for(const auto& h: hits)
0077         result.push_back(h.position.y); //mm
0078   return result;
0079   };
0080
0081   // Hit position Z
0082   auto hit_z_position = [&] (const std::vector<edm4hep::SimCalorimeterHitData>& hits) {
0083         std::vector<double> result;
0084         for(const auto& h: hits)
0085         result.push_back(h.position.z); //mm
0086   return result;
0087   };
0088
0089   // Define variables
0090   auto d1 = d0
0091     .Define("Ethr", Ethr, {"MCParticles"})
0092     .Define("nhits_Ecal", nhits, {"EcalFarForwardZDCHits"})
0093     .Define("Esim_Ecal", Esim, {"EcalFarForwardZDCHits"})
0094     .Define("fsam_Ecal", fsam, {"Esim_Ecal", "Ethr"})
0095     .Define("nhits_Hcal", nhits, {"HcalFarForwardZDCHits"})
0096     .Define("Esim_Hcal", Esim, {"HcalFarForwardZDCHits"})
0097     .Define("fsam_Hcal", fsam, {"Esim_Hcal", "Ethr"})
0098     .Define("hit_x_position", hit_x_position, {"HcalFarForwardZDCHits"})
0099       .Define("hit_y_position", hit_y_position, {"HcalFarForwardZDCHits"})
0100   ;
0101
0102   // Define Histograms
0103   auto hEthr = d1.Histo1D(
0104       {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5},
0105       "Ethr");
0106   auto hNhits_Ecal =
0107       d1.Histo1D({"hNhits_Ecal", "Number of Ecal hits per events; Number of Ecal hits; Events",
0108                   100, 0.0, 2000.0},
0109                  "nhits_Ecal");
0110   auto hEsim_Ecal = d1.Histo1D(
0111       {"hEsim_Ecal", "Ecal Energy Deposit; Ecal Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
0112       "Esim_Ecal");
0113   auto hfsam_Ecal = d1.Histo1D(
0114       {"hfsam_Ecal", "Ecal Sampling Fraction; Ecal Sampling Fraction; Events", 100, 0.0, 0.1},
0115       "fsam_Ecal");
0116   auto hNhits_Hcal =
0117       d1.Histo1D({"hNhits_Hcal", "Number of Hcal hits per events; Number of Hcal hits; Events",
0118                   100, 0.0, 2000.0},
0119                  "nhits_Hcal");
0120   auto hEsim_Hcal = d1.Histo1D(
0121       {"hEsim_Hcal", "Hcal Energy Deposit; Hcal Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
0122       "Esim_Hcal");
0123   auto hfsam_Hcal = d1.Histo1D(
0124       {"hfsam_Hcal", "Hcal Sampling Fraction; Hcal Sampling Fraction; Events", 100, 0.0, 0.1},
0125       "fsam_Hcal");
0126
0127   auto hXY_HitMap_Hcal = d1.Histo2D({"hXY_HitMap_Hcal", "hit position Y vs. X histogram; hit position X [mm]; hit position Y [mm]", 8,-1300,-500, 8, -400, 400}, "hit_x_position", "hit_y_position");
0128
0129
0130   // Event Counts
0131   auto nevents_thrown = d1.Count();
0132   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
0133
0134   // Draw Histograms
0135   {
0136     TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
0137     c1->SetLogy(1);
0138     auto h = hEthr->DrawCopy();
0139     //h->GetYaxis()->SetTitleOffset(1.4);
0140     h->SetLineWidth(2);
0141     h->SetLineColor(kBlue);
0142     c1->SaveAs(TString(results_path + "/zdc_Ethr.png"));
0143     c1->SaveAs(TString(results_path + "/zdc_Ethr.pdf"));
0144   }
0145
0146   // Ecal
0147   {
0148     TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
0149     c2->SetLogy(1);
0150     auto h = hNhits_Ecal->DrawCopy();
0151     //h->GetYaxis()->SetTitleOffset(1.4);
0152     h->SetLineWidth(2);
0153     h->SetLineColor(kBlue);
0154     c2->SaveAs(TString(results_path + "/zdc_nhits_Ecal.png"));
0155     c2->SaveAs(TString(results_path + "/zdc_nhits_Ecal.pdf"));
0156   }
0157
0158   {
0159     TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
0160     c3->SetLogy(1);
0161     auto h = hEsim_Ecal->DrawCopy();
0162     //h->GetYaxis()->SetTitleOffset(1.4);
0163     h->SetLineWidth(2);
0164     h->SetLineColor(kBlue);
0165     c3->SaveAs(TString(results_path + "/zdc_Esim_Ecal.png"));
0166     c3->SaveAs(TString(results_path + "/zdc_Esim_Ecal.pdf"));
0167   }
0168
0169   std::cout << "derp4\n";
0170   {
0171     TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
0172     c4->SetLogy(1);
0173     auto h = hfsam_Ecal->DrawCopy();
0174     //h->GetYaxis()->SetTitleOffset(1.4);
0175     h->SetLineWidth(2);
0176     h->SetLineColor(kBlue);
0177     //h->Fit("gaus", "", "", 0.01, 0.1);
0178     //h->GetFunction("gaus")->SetLineWidth(2);
0179     //h->GetFunction("gaus")->SetLineColor(kRed);
0180     c4->SaveAs(TString(results_path + "/zdc_fsam_Ecal.png"));
0181     c4->SaveAs(TString(results_path + "/zdc_fsam_Ecal.pdf"));
0182   }
0183
0184   // Hcal
0185   {
0186     TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
0187     c2->SetLogy(1);
0188     auto h = hNhits_Hcal->DrawCopy();
0189     //h->GetYaxis()->SetTitleOffset(1.4);
0190     h->SetLineWidth(2);
0191     h->SetLineColor(kBlue);
0192     c2->SaveAs(TString(results_path + "/zdc_nhits_Hcal.png"));
0193     c2->SaveAs(TString(results_path + "/zdc_nhits_Hcal.pdf"));
0194   }
0195
0196   {
0197     TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
0198     c3->SetLogy(1);
0199     auto h = hEsim_Hcal->DrawCopy();
0200     //h->GetYaxis()->SetTitleOffset(1.4);
0201     h->SetLineWidth(2);
0202     h->SetLineColor(kBlue);
0203     c3->SaveAs(TString(results_path + "/zdc_Esim_Hcal.png"));
0204     c3->SaveAs(TString(results_path + "/zdc_Esim_Hcal.pdf"));
0205   }
0206
0207   std::cout << "derp4\n";
0208   {
0209     TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
0210     c4->SetLogy(1);
0211     auto h = hfsam_Hcal->DrawCopy();
0212     //h->GetYaxis()->SetTitleOffset(1.4);
0213     h->SetLineWidth(2);
0214     h->SetLineColor(kBlue);
0215     //h->Fit("gaus", "", "", 0.01, 0.1);
0216     //h->GetFunction("gaus")->SetLineWidth(2);
0217     //h->GetFunction("gaus")->SetLineColor(kRed);
0218     c4->SaveAs(TString(results_path + "/zdc_fsam_Hcal.png"));
0219     c4->SaveAs(TString(results_path + "/zdc_fsam_Hcal.pdf"));
0220   }
0221   {
0222     TCanvas* c5 = new TCanvas("c5", "c5", 600, 500);
0223     hXY_HitMap_Hcal->Draw("COLZ");
0224     c5->SaveAs(TString(results_path + "/zdc_xy_map_Hcal.png"));
0225     c5->SaveAs(TString(results_path + "/zdc_xy_map_Hcal.pdf"));
0226   }
0227 }