Resolve "Add ZDC benchmark"

Closes #71 (closed)

• neutron example
• photon example

TODO: check crossing angle and make a separate issue on optimizing clustering algorithm

benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx

+////////////////////////////////////////
+// Read ROOT output file
+// Plot variables
+// Jihee Kim 09/2021
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <vector>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TGaxis.h"
+#include "TPaveStats.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TVector3.h"
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+#include "fmt/core.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+
+using ROOT::RDataFrame;
+
+void analysis_zdc_neutrons(const char* input_fname = "sim_zdc_uniform_neutron.root")
+{
+  // Setting for graphs
+  TStyle* kStyle = new TStyle("kStyle","Kim's Style");
+  kStyle->SetOptStat("emr");
+  kStyle->SetOptTitle(0);
+  kStyle->SetOptFit(1101);
+  kStyle->SetStatColor(0);
+  kStyle->SetStatW(0.15);
+  kStyle->SetStatH(0.10);
+  kStyle->SetStatX(0.9);
+  kStyle->SetStatY(0.9);
+  kStyle->SetStatBorderSize(1);
+  kStyle->SetLabelSize(0.045,"xyz");
+  kStyle->SetTitleSize(0.045,"xyz");
+  kStyle->SetTitleOffset(1.2,"y");
+  kStyle->SetLineWidth(2);
+  kStyle->SetTitleFont(42,"xyz");
+  kStyle->SetLabelFont(42,"xyz");
+  kStyle->SetCanvasDefW(700);
+  kStyle->SetCanvasDefH(500);
+  kStyle->SetCanvasColor(0);
+  kStyle->SetPadTickX(1);
+  kStyle->SetPadTickY(1);
+  kStyle->SetPadGridX(1);
+  kStyle->SetPadGridY(1);
+  kStyle->SetPadLeftMargin(0.1);
+  kStyle->SetPadRightMargin(0.1);
+  kStyle->SetPadTopMargin(0.1);
+  kStyle->SetPadBottomMargin(0.1);
+  TGaxis::SetMaxDigits(3);
+  gROOT->SetStyle("kStyle");
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Thrown Energy [GeV]
+  auto Ethr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.ps.x * p.ps.x + p.ps.y * p.ps.y + p.ps.z * p.ps.z + p.mass * p.mass);
+    return energy;
+  };
+
+  // Theta [mrad]
+  auto Thetathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto theta = p.ps.theta();
+    return theta*1000.0;
+  };
+
+  // Phi [rad]
+  auto Phithr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto phi = p.ps.phi();
+    return phi;
+  };
+  
+  // Eta
+  auto Etathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto eta = p.ps.eta();
+    return eta;
+  };  
+
+  // Define variables
+  auto d1 = d0.Define("Ethr",      Ethr,     {"mcparticles"})
+              .Define("Thetathr",  Thetathr, {"mcparticles"})
+              .Define("Phithr",    Phithr,   {"mcparticles"})
+	      .Define("Etathr",    Etathr,   {"mcparticles"})
+	      ;
+
+  // Define Histograms
+  auto hEthr        = d1.Histo1D({"hEthr",        "Thrown Energy; E_{thr} [GeV]; Events",                100, 0.0, 0.3},  "Ethr");
+  auto hThetathr    = d1.Histo1D({"hThetathr",    "Thrown Theta; #theta_{thr} [mrad]; Events",           100, 0.0, 5.0},  "Thetathr");
+  auto hPhithr      = d1.Histo1D({"hPhithr",      "Thrown Phi; #phi_{thr} [rad]; Events",                100, 0.0, 5.0},  "Phithr");
+  auto hEtathr      = d1.Histo1D({"hEtathr",      "Thrown Eta; #eta_{thr}; Events",                      100, 6.0, 15.0}, "Etathr");
+  auto hPhiThetathr = d1.Histo2D({"hPhiThetathr", "Thrown Phi vs Theta; #theta_{thr} [mrad]; #phi_{thr}", 20, 0.0, 5.0, 20, 0.0, 5.0}, 
+                  "Thetathr", "Phithr");
+
+  // Event Counts
+  auto nevents_thrown = d1.Count();
+  std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n"; 
+
+  // Draw Histograms
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1");
+    auto h = hEthr->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c1->SaveAs("results/far_forward/zdc/zdc_neutron_Ethr.png");
+    c1->SaveAs("results/far_forward/zdc/zdc_neutron_Ethr.pdf");
+  }
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2");
+    auto h = hThetathr->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c2->SaveAs("results/far_forward/zdc/zdc_neutron_Thetathr.png");
+    c2->SaveAs("results/far_forward/zdc/zdc_neutron_Thetathr.pdf");
+  }
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3");
+    auto h = hPhithr->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c3->SaveAs("results/far_forward/zdc/zdc_neutron_Phithr.png");
+    c3->SaveAs("results/far_forward/zdc/zdc_neutron_Phithr.pdf");
+  }
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4");
+    auto h = hEtathr->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c4->SaveAs("results/far_forward/zdc/zdc_neutron_Etathr.png");
+    c4->SaveAs("results/far_forward/zdc/zdc_neutron_Etathr.pdf");
+  }
+  {
+    TCanvas* c5 = new TCanvas("c5", "c5");
+    auto h = hPhiThetathr->DrawCopy("COLZ");
+    c5->SaveAs("results/far_forward/zdc/zdc_neutron_PhiThetathr.png");
+    c5->SaveAs("results/far_forward/zdc/zdc_neutron_PhiThetathr.pdf");
+  } 
+}
+