pid/scripts/mrich_analysis.cxx

0001 ////////////////////////////////////////
0002 // Read reconstruction ROOT output file
0003 // Plot variables
0004 ////////////////////////////////////////
0005
0006 #include "ROOT/RDataFrame.hxx"
0007 #include <iostream>
0008 #include <fmt/core.h>
0009
0010 #include "dd4pod/Geant4ParticleCollection.h"
0011 #include "dd4pod/PhotoMultiplierHitCollection.h"
0012
0013 #include "TCanvas.h"
0014 #include "TStyle.h"
0015 #include "TMath.h"
0016 #include "TH1.h"
0017 #include "TF1.h"
0018 #include "TH1D.h"
0019
0020 using ROOT::RDataFrame;
0021 using namespace ROOT::VecOps;
0022
0023 void mrich_analysis(const char* input_fname = "sim_output/sim_pid_backward_e-_5GeV.edm4hep.root", const char* input_pname = "e-")
0024 {
0025   // Setting for graphs
0026   gROOT->SetStyle("Plain");
0027   gStyle->SetOptFit(1);
0028   gStyle->SetLineWidth(2);
0029   gStyle->SetPadTickX(1);
0030   gStyle->SetPadTickY(1);
0031   gStyle->SetPadGridX(1);
0032   gStyle->SetPadGridY(1);
0033   gStyle->SetPadLeftMargin(0.14);
0034   gStyle->SetPadRightMargin(0.14);
0035
0036   ROOT::EnableImplicitMT();
0037   ROOT::RDataFrame d0("events", input_fname);
0038
0039   // Number of hits
0040   auto nhits = [] (const std::vector<dd4pod::PhotoMultiplierHitData>& evt) {return (int) evt.size(); };
0041
0042   // Define variables
0043   auto d1 = d0.Define("nhits", nhits, {"MRICHHits"});
0044
0045   // Define Histograms
0046   auto hNhits =
0047       d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
0048                   100, 0.0, 2000.0},
0049                  "nhits");
0050   auto hXYhits =
0051       d1.Histo2D({"hXYhits", "Hit positions for events; Horizontal position [mm]; Vertical position [mm]",
0052                   1000, -1000.0, +1000.0, 1000, -1000.0, +1000.0},
0053                  "MRICHHits.position.x", "MRICHHits.position.y");
0054
0055   // Event Counts
0056   auto nevents_thrown = d1.Count();
0057   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
0058
0059   // Draw Histograms
0060   {
0061     TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
0062     auto h = hXYhits->DrawCopy();
0063     c1->SaveAs(fmt::format("results/pid/backward/mrich/mrich_{}_hits_xy.png",input_pname).c_str());
0064     c1->SaveAs(fmt::format("results/pid/backward/mrich/mrich_{}_hits_xy.pdf",input_pname).c_str());
0065   }
0066   {
0067     TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
0068     c2->SetLogy(1);
0069     auto h = hNhits->DrawCopy();
0070     //h->GetYaxis()->SetTitleOffset(1.4);
0071     h->SetLineWidth(2);
0072     h->SetLineColor(kBlue);
0073     c2->SaveAs(fmt::format("results/pid/backward/mrich/mrich_{}_nhits.png",input_pname).c_str());
0074     c2->SaveAs(fmt::format("results/pid/backward/mrich/mrich_{}_nhits.pdf",input_pname).c_str());
0075   }
0076 }