From 3cde54b793619821d414392f0ed43941d7470e1c Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Tue, 20 Oct 2020 13:11:49 -0500
Subject: [PATCH] Included pi0 script in the benchmarks CI
Moved datasets scripts over benchmarks ecal/scripts directory
---
ecal/ecal_config.yml | 12 ++
ecal/emcal_pi0s.sh | 63 ++++++
ecal/scripts/emcal_pi0.cxx | 76 +++++++
ecal/scripts/emcal_pi0_reader.cxx | 126 ++++++++++++
ecal/scripts/makeplot_pi0.C | 315 ++++++++++++++++++++++++++++++
5 files changed, 592 insertions(+)
create mode 100644 ecal/emcal_pi0s.sh
create mode 100644 ecal/scripts/emcal_pi0.cxx
create mode 100644 ecal/scripts/emcal_pi0_reader.cxx
create mode 100644 ecal/scripts/makeplot_pi0.C
diff --git a/ecal/ecal_config.yml b/ecal/ecal_config.yml
index d445eb8c..69ffb803 100644
--- a/ecal/ecal_config.yml
+++ b/ecal/ecal_config.yml
@@ -12,3 +12,15 @@ ecal_1_emcal_electrons:
script:
- bash ecal/emcal_electrons.sh
+ecal_1_emcal_pi0s:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
+ tags:
+ - silicon
+ timeout: 12 hours 30 minutes
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ stage: run
+ script:
+ - bash ecal/emcal_pi0s.sh
diff --git a/ecal/emcal_pi0s.sh b/ecal/emcal_pi0s.sh
new file mode 100644
index 00000000..dd6e9caa
--- /dev/null
+++ b/ecal/emcal_pi0s.sh
@@ -0,0 +1,63 @@
+#!/bin/bash
+# Based on emcal_electrons.sh script
+
+# Detector
+if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then
+ export JUGGLER_DETECTOR="topside"
+fi
+# Number of events
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=100
+fi
+# File names
+export JUGGLER_FILE_NAME_TAG="emcal_uniform_pi0s"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+# Datasets
+#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+root -b -q "ecal/scripts/emcal_pi0.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "ecal/scripts/emcal_pi0_reader.cxx(\"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+# Detector TOPSiDE
+git clone https://eicweb.phy.anl.gov/EIC/detectors/${JUGGLER_DETECTOR}.git
+mkdir ${JUGGLER_DETECTOR}/build
+pushd ${JUGGLER_DETECTOR}/build
+cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
+popd
+
+echo "CHECK POINT FOR GEANT4 SIMULATION"
+pwd
+pushd ${JUGGLER_DETECTOR}
+pwd
+ls -l
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+echo "JUGGLER_FILE_NAME_TAG = ${JUGGLER_FILE_NAME_TAG}"
+echo "JUGGLER_SIM_FILE = ${JUGGLER_SIM_FILE}"
+# run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${JUGGLER_DETECTOR}.xml \
+ --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile ${JUGGLER_SIM_FILE}
+# Need to figure out how to pass file name to juggler from the commandline
+
+xenv -x /usr/local/Juggler.xenv gaudirun.py ../ecal/options/example_crystal.py
+ls -l
+
+popd
+ls -l
+
+pwd
+mkdir -p results
+root -b -q "ecal/scripts/makeplot_pi0.C(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+
+mkdir -p sim_output
+cp "${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}" sim_output/.
+cp "${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}" sim_output/.
+
diff --git a/ecal/scripts/emcal_pi0.cxx b/ecal/scripts/emcal_pi0.cxx
new file mode 100644
index 00000000..e07101d3
--- /dev/null
+++ b/ecal/scripts/emcal_pi0.cxx
@@ -0,0 +1,76 @@
+//////////////////////////////////////////////////////////////
+// Crystal EMCAL detector
+// Single Pion dataset
+// J.KIM 10/18/2020
+//////////////////////////////////////////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include <iostream>
+#include<random>
+#include<cmath>
+#include <math.h>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_100kEvt.hepmc")
+{
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom *r1 = new TRandom();
+
+ for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+ // FourVector(px,py,pz,e,pdgid,status)
+ // type 4 is beam
+ // pdgid 11 - electron
+ // pdgid 111 - pi0
+ // pdgid 2212 - proton
+ GenParticlePtr p1 =
+ std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+ GenParticlePtr p2 = std::make_shared<GenParticle>(
+ FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+ // Define momentum
+ Double_t p = r1->Uniform(0.0,30.0);
+ Double_t px;
+ Double_t py;
+ Double_t pz;
+ // Generates random vectors, uniformly distributed over the surface of a
+ // sphere of given radius, in this case momentum.
+ r1->Sphere(px, py, pz, p);
+
+ // type 1 is final state
+ // pdgid 111 - pi0 135 MeV/c^2
+ GenParticlePtr p3 = std::make_shared<GenParticle>(
+ FourVector(
+ px, py, pz,
+ sqrt(p*p + (0.134976*0.134976))),
+ 111, 1);
+
+ GenVertexPtr v1 = std::make_shared<GenVertex>();
+ v1->add_particle_in(p1);
+ v1->add_particle_in(p2);
+
+ v1->add_particle_out(p3);
+ evt.add_vertex(v1);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+
+ hepmc_output.write_event(evt);
+ if (events_parsed % 10000 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/ecal/scripts/emcal_pi0_reader.cxx b/ecal/scripts/emcal_pi0_reader.cxx
new file mode 100644
index 00000000..c43958ac
--- /dev/null
+++ b/ecal/scripts/emcal_pi0_reader.cxx
@@ -0,0 +1,126 @@
+//////////////////////////
+// Crystal EMCAL detector
+// Singe Pion dataset
+// J.KIM 10/18/2020
+//////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TH1F.h"
+#include <iostream>
+#include "TStyle.h"
+
+using namespace HepMC3;
+
+void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc"){
+
+ // Setting for graphs
+ gROOT->SetStyle("Plain");
+ gStyle->SetOptFit(1);
+ gStyle->SetLineWidth(2);
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.14);
+
+ ReaderAscii hepmc_input(in_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Histograms
+ TH1F* h_pi0_energy = new TH1F("h_pi0_energy", "pi0 energy;E [GeV];Events", 100,-0.5,30.5);
+ TH1F* h_pi0_eta = new TH1F("h_pi0_eta", "pi0 #eta;#eta;Events", 100,-10.0,10.0);
+ TH1F* h_pi0_theta = new TH1F("h_pi0_theta", "pi0 #theta;#theta [degree];Events", 100,-0.5,180.5);
+ TH1F* h_pi0_phi = new TH1F("h_pi0_phi", "pi0 #phi;#phi [degree];Events", 100,-180.5,180.5);
+ 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);
+ 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);
+ 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);
+
+ while(!hepmc_input.failed()) {
+ // Read event from input file
+ hepmc_input.read_event(evt);
+ // If reading failed - exit loop
+ if( hepmc_input.failed() ) break;
+
+ for(const auto& v : evt.vertices() ) {
+ for(const auto& p : v->particles_out() ) {
+ if(p->pid() == 111) {
+ h_pi0_energy->Fill(p->momentum().e());
+ h_pi0_eta->Fill(p->momentum().eta());
+ h_pi0_theta->Fill(p->momentum().theta()*TMath::RadToDeg());
+ h_pi0_phi->Fill(p->momentum().phi()*TMath::RadToDeg());
+ h_pi0_pzpt->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + p->momentum().py()*p->momentum().py()),p->momentum().pz());
+ h_pi0_pxpy->Fill(p->momentum().px(),p->momentum().py());
+ h_pi0_p->Fill(p->momentum().px(),p->momentum().py(),p->momentum().pz());
+ }
+ }
+ }
+ evt.clear();
+ events_parsed++;
+ }
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+
+ TCanvas* c = new TCanvas("c","c", 500, 500);
+ h_pi0_energy->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_energy->SetLineWidth(2);
+ h_pi0_energy->SetLineColor(kBlue);
+ h_pi0_energy->DrawClone();
+ c->SaveAs("results/input_pi0_energy_0GeVto30GeV.png");
+ c->SaveAs("results/input_pi0_energy_0GeVto30GeV.pdf");
+
+ TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
+ h_pi0_eta->GetYaxis()->SetTitleOffset(1.9);
+ h_pi0_eta->SetLineWidth(2);
+ h_pi0_eta->SetLineColor(kBlue);
+ h_pi0_eta->DrawClone();
+ c1->SaveAs("results/input_pi0_eta_0GeVto30GeV.png");
+ c1->SaveAs("results/input_pi0_eta_0GeVto30GeV.pdf");
+
+ TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
+ h_pi0_theta->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_theta->SetLineWidth(2);
+ h_pi0_theta->SetLineColor(kBlue);
+ h_pi0_theta->DrawClone();
+ c2->SaveAs("results/input_pi0_theta_0GeVto30GeV.png");
+ c2->SaveAs("results/input_pi0_theta_0GeVto30GeV.pdf");
+
+ TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
+ h_pi0_phi->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_phi->SetLineWidth(2);
+ h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+100.0);
+ h_pi0_phi->SetLineColor(kBlue);
+ h_pi0_phi->DrawClone();
+ c3->SaveAs("results/input_pi0_phi_0GeVto30GeV.png");
+ c3->SaveAs("results/input_pi0_phi_0GeVto30GeV.pdf");
+
+ TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
+ h_pi0_pzpt->GetYaxis()->SetTitleOffset(1.4);
+ h_pi0_pzpt->SetLineWidth(2);
+ h_pi0_pzpt->SetLineColor(kBlue);
+ h_pi0_pzpt->DrawClone("COLZ");
+ c4->SaveAs("results/input_pi0_pzpt_0GeVto30GeV.png");
+ c4->SaveAs("results/input_pi0_pzpt_0GeVto30GeV.pdf");
+
+ TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
+ h_pi0_pxpy->GetYaxis()->SetTitleOffset(1.4);
+ h_pi0_pxpy->SetLineWidth(2);
+ h_pi0_pxpy->SetLineColor(kBlue);
+ h_pi0_pxpy->DrawClone("COLZ");
+ c5->SaveAs("results/input_pi0_pxpy_0GeVto30GeV.png");
+ c5->SaveAs("results/input_pi0_pxpy_0GeVto30GeV.pdf");
+
+ TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
+ h_pi0_p->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_p->GetXaxis()->SetTitleOffset(1.6);
+ h_pi0_p->GetZaxis()->SetTitleOffset(1.6);
+ h_pi0_p->SetLineWidth(2);
+ h_pi0_p->SetLineColor(kBlue);
+ h_pi0_p->DrawClone();
+ c6->SaveAs("results/input_pi0_p_0GeVto30GeV.png");
+ c6->SaveAs("results/input_pi0_p_0GeVto30GeV.pdf");
+
+}
diff --git a/ecal/scripts/makeplot_pi0.C b/ecal/scripts/makeplot_pi0.C
new file mode 100644
index 00000000..94c5bf5b
--- /dev/null
+++ b/ecal/scripts/makeplot_pi0.C
@@ -0,0 +1,315 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+int makeplot_pi0(const char* input_fname = "../sim_output/sim_emcal_pi0s_output.root")
+{
+ // Setting figures
+ gROOT->SetStyle("Plain");
+ gStyle->SetLineWidth(3);
+ gStyle->SetOptStat("nem");
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.14);
+
+ // Input ROOT file
+ TFile *f = new TFile(input_fname,"READ");
+ TTree *t = (TTree *)f->Get("events");
+
+ // Set Branch status and addressed
+ t->SetMakeClass(1);
+ t->SetBranchStatus("*", 0);
+
+ Int_t EcalClusters_;
+ t->SetBranchStatus("EcalClusters", 1);
+ t->SetBranchAddress("EcalClusters", &EcalClusters_);
+
+ Int_t RecoEcalHits_;
+ t->SetBranchStatus("RecoEcalHits", 1);
+ t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
+
+ const Int_t kMaxEcalClusters = 100000;
+ Double_t cluster_x_pos[kMaxEcalClusters];
+ Double_t cluster_y_pos[kMaxEcalClusters];
+ Double_t cluster_z_pos[kMaxEcalClusters];
+ Float_t cluster_energy[kMaxEcalClusters];
+ t->SetBranchStatus("EcalClusters.position.x",1);
+ t->SetBranchStatus("EcalClusters.position.y",1);
+ t->SetBranchStatus("EcalClusters.position.z",1);
+ t->SetBranchStatus("EcalClusters.energy",1);
+ t->SetBranchAddress("EcalClusters.position.x",cluster_x_pos);
+ t->SetBranchAddress("EcalClusters.position.y",cluster_y_pos);
+ t->SetBranchAddress("EcalClusters.position.z",cluster_z_pos);
+ t->SetBranchAddress("EcalClusters.energy",cluster_energy);
+
+ const Int_t kMaxRecoEcalHits = 100000;
+ Double_t rec_x_pos[kMaxRecoEcalHits];
+ Double_t rec_y_pos[kMaxRecoEcalHits];
+ Double_t rec_energy[kMaxRecoEcalHits];
+ t->SetBranchStatus("RecoEcalHits.position.x",1);
+ t->SetBranchStatus("RecoEcalHits.position.y",1);
+ t->SetBranchStatus("RecoEcalHits.energy",1);
+ t->SetBranchAddress("RecoEcalHits.position.x",rec_x_pos);
+ t->SetBranchAddress("RecoEcalHits.position.y",rec_y_pos);
+ t->SetBranchAddress("RecoEcalHits.energy",rec_energy);
+
+ // Setting for Canvas
+ TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
+ TCanvas *c2 = new TCanvas("c2", "c2", 500, 500);
+ TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
+ TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
+ TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
+ TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
+ TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
+ TCanvas *c8 = new TCanvas("c8", "c8", 500, 500);
+ TCanvas *c9 = new TCanvas("c9", "c9", 500, 500);
+ TCanvas *c10 = new TCanvas("c10","c10", 500, 500);
+ TCanvas *c11 = new TCanvas("c11","c11", 500, 500);
+ TCanvas *c12 = new TCanvas("c12","c12", 500, 500);
+ TCanvas *c13 = new TCanvas("c13","c13", 500, 500);
+
+ // Declare histograms
+ TH1D *h1 = new TH1D("h1", "Scattering Angle(#theta)", 100,130.0,180.0);
+ TH1D *h2 = new TH1D("h2", "Pseudo-rapidity(#eta)", 100,-5.0,0.0);
+ TH2D *h3 = new TH2D("h3", "Cluster E vs Pseudo-rapidity", 100,-0.5,30.5,100,-5.0,0.0);
+ TH1D *h4 = new TH1D("h4", "Reconstructed energy per event", 100,-0.5,30.5);
+ TH1D *h5 = new TH1D("h5", "Number of Clusters per event", 5,-0.5,4.5);
+ TH1D *h6 = new TH1D("h6", "Scattering Angle(#theta) with CUT", 100,130.0,180.0);
+ TH1D *h7 = new TH1D("h7", "Pseudo-rapidity(#eta) with CUT", 100,-5.0,0.0);
+ TH2D *h8 = new TH2D("h8", "Cluster Hit Position", 62,-62.0,62.0,62,-62.0,62.0);
+ TH2D *h9 = new TH2D("h9", "All Hit Position", 62,-62.0,62.0,62,-62.0,62.0);
+ TH1D *h10 = new TH1D("h10","Invariant mass", 60,0.0,300.0);
+ TH1D *h11 = new TH1D("h11","E1", 100,-0.5,30.5);
+ TH1D *h12 = new TH1D("h12","E2", 100,-0.5,30.5);
+ TH1D *h13 = new TH1D("h13","angle", 100,0.0,90.0);
+
+
+ // Declare ellipse for boundary of crystal calorimeter
+ TEllipse *ell1 = new TEllipse(0.0, 0.0, 60.0, 60.0);
+ ell1->SetFillStyle(4000);
+ TEllipse *ell2 = new TEllipse(0.0, 0.0, 12.0, 12.0);
+ ell2->SetFillStyle(4000);
+
+ // Total number of entries
+ Int_t nentries = t->GetEntries();
+
+ // Variables are used in calculation
+ Double_t r; // Radius [cm]
+ Double_t phi; // Azimuth [degree]
+ Double_t theta; // Inclination [degree]
+ Double_t eta; // Pseudo-rapidity [unitless]
+ Float_t cluster_e; // Cluster energy [GeV]
+ Float_t total_cluster_e; // Add up clusters per event [GeV]
+ Double_t dot_product_pos_clusters; // dot product of positions of two photons
+ Double_t mag_pos2_cluster_1; // squared magnitude of position
+ Double_t mag_pos2_cluster_2; // squared magnitude of position
+ Double_t cosine_clusters; // cos(theta_photons)
+ Double_t theta_photons; // angle between two photons
+ Double_t invariant_mass; // M^2 = 2 * p_1 * p_2 * (1 - cos(theta_photons))
+
+ // Loop over event by event
+ for (int ievent = 0; ievent < nentries; ievent++)
+ {
+ t->GetEntry(ievent);
+
+ Int_t ncluster = EcalClusters_;
+ Int_t nreconhits = RecoEcalHits_;
+
+ total_cluster_e = 0.0;
+
+ h5->Fill(ncluster, 1.0);
+
+ // Loop over cluster by cluster
+ for (int icluster=0; icluster < ncluster; icluster++)
+ {
+ r = TMath::Sqrt((cluster_x_pos[icluster]*cluster_x_pos[icluster]) +
+ (cluster_y_pos[icluster]*cluster_y_pos[icluster]) +
+ (cluster_z_pos[icluster]*cluster_z_pos[icluster]));
+ phi = TMath::ATan(cluster_y_pos[icluster]/cluster_x_pos[icluster]) * TMath::RadToDeg();
+ theta = TMath::ACos(cluster_z_pos[icluster] / r) * TMath::RadToDeg();
+ eta = -1.0 * TMath::Log(TMath::Tan((theta*TMath::DegToRad())/2.0));
+ cluster_e = cluster_energy[icluster] / 1.e+3;
+ total_cluster_e += cluster_e;
+ }
+
+ // Select events with two clusters
+ // To calculate invariant mass
+ // M^2 = 2p1p2(1-cos(theta))
+ // p1 = E1
+ // p2 = E2
+ // theta: angle between two photons
+ if(ncluster == 2)
+ {
+ h1->Fill(theta, 1.0);
+ h2->Fill(eta, 1.0);
+ h3->Fill(cluster_e, eta, 1.0);
+ h4->Fill(total_cluster_e, 1.0);
+ h8->Fill(cluster_x_pos[0],cluster_y_pos[0], 1.0);
+
+ if(total_cluster_e > 0.5)
+ {
+ h6->Fill(theta, 1.0);
+ h7->Fill(eta, 1.0);
+ }
+
+ // Calculate invariant mass
+ dot_product_pos_clusters = cluster_x_pos[0]*cluster_x_pos[1] + cluster_y_pos[0]*cluster_y_pos[1] + cluster_z_pos[0]*cluster_z_pos[1];
+ mag_pos2_cluster_1 = (cluster_x_pos[0]*cluster_x_pos[0]) + (cluster_y_pos[0]*cluster_y_pos[0]) + (cluster_z_pos[0]*cluster_z_pos[0]);
+ mag_pos2_cluster_2 = (cluster_x_pos[1]*cluster_x_pos[1]) + (cluster_y_pos[1]*cluster_y_pos[1]) + (cluster_z_pos[1]*cluster_z_pos[1]);
+ cosine_clusters = (dot_product_pos_clusters/TMath::Sqrt(mag_pos2_cluster_1*mag_pos2_cluster_2));
+ theta_photons = TMath::ACos(cosine_clusters)*TMath::RadToDeg();
+
+ invariant_mass = TMath::Sqrt(2.0*cluster_energy[0]*cluster_energy[1]*(1.0 - cosine_clusters));
+
+ // Fill histograms
+ h10->Fill(invariant_mass, 1.0);
+ h11->Fill(cluster_energy[0], 1.0);
+ h12->Fill(cluster_energy[1], 1.0);
+ h13->Fill(theta_photons, 1.0);
+ }
+
+ // Loop over hit by hit
+ for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
+ h9->Fill(rec_x_pos[ireconhit],rec_y_pos[ireconhit], 1.0);
+ }
+
+ // Drawing and Saving figures
+ c1->cd();
+ h1->SetLineColor(kBlue);
+ h1->SetLineWidth(2);
+ h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+10.0);
+ h1->GetXaxis()->SetTitle("#theta [degree]");
+ h1->GetYaxis()->SetTitle("events");
+ h1->GetYaxis()->SetTitleOffset(1.4);
+ h1->DrawClone();
+ c1->SaveAs("results/pi0_theta_hist_0GeVto30GeV.png");
+ c1->SaveAs("results/pi0_theta_hist_0GeVto30GeV.pdf");
+
+ c2->cd();
+ h2->SetLineColor(kBlue);
+ h2->SetLineWidth(2);
+ h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+10.0);
+ h2->GetXaxis()->SetTitle("#eta");
+ h2->GetYaxis()->SetTitle("events");
+ h2->GetYaxis()->SetTitleOffset(1.4);
+ h2->DrawClone();
+ c2->SaveAs("results/pi0_eta_hist_0GeVto30GeV.png");
+ c2->SaveAs("results/pi0_eta_hist_0GeVto30GeV.pdf");
+
+ c3->cd();
+ h3->GetXaxis()->SetTitle("Cluster energy [GeV]");
+ h3->GetYaxis()->SetTitle("#eta");
+ h3->GetYaxis()->SetTitleOffset(1.4);
+ h3->DrawClone("COLZ");
+ c3->SaveAs("results/pi0_E_vs_eta_hist_0GeVto30GeV.png");
+ c3->SaveAs("results/pi0_E_vs_eta_hist_0GeVto30GeV.pdf");
+
+ c4->cd();
+ c4->SetLogy(1);
+ h4->SetLineColor(kBlue);
+ h4->SetLineWidth(2);
+ h4->GetXaxis()->SetTitle("reconstructed energy [GeV]");
+ h4->GetYaxis()->SetTitle("events");
+ h4->GetYaxis()->SetTitleOffset(1.4);
+ h4->DrawClone();
+ c4->SaveAs("results/pi0_Erec_hist_0GeVto30GeV.png");
+ c4->SaveAs("results/pi0_Erec_hist_0GeVto30GeV.pdf");
+
+ c5->cd();
+ c5->SetLogy(1);
+ h5->SetLineColor(kBlue);
+ h5->SetLineWidth(2);
+ h5->GetXaxis()->SetTitle("Number of Clusters");
+ h5->GetYaxis()->SetTitle("events");
+ h5->GetYaxis()->SetTitleOffset(1.4);
+ h5->DrawClone();
+ c5->SaveAs("results/pi0_ncluster_hist_0GeVto30GeV.png");
+ c5->SaveAs("results/pi0_ncluster_hist_0GeVto30GeV.pdf");
+
+ c6->cd();
+ h6->SetLineColor(kBlue);
+ h6->SetLineWidth(2);
+ h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+10.0);
+ h6->GetXaxis()->SetTitle("#theta [degree]");
+ h6->GetYaxis()->SetTitle("events");
+ h6->GetYaxis()->SetTitleOffset(1.4);
+ h6->DrawClone();
+ c6->SaveAs("results/pi0_theta_hist_CUT_0GeVto30GeV.png");
+ c6->SaveAs("results/pi0_theta_hist_CUT_0GeVto30GeV.pdf");
+
+ c7->cd();
+ h7->SetLineColor(kBlue);
+ h7->SetLineWidth(2);
+ h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+10.0);
+ h7->GetXaxis()->SetTitle("#eta");
+ h7->GetYaxis()->SetTitle("events");
+ h7->GetYaxis()->SetTitleOffset(1.4);
+ h7->DrawClone();
+ c7->SaveAs("results/pi0_eta_hist_CUT_0GeVto30GeV.png");
+ c7->SaveAs("results/pi0_eta_hist_CUT_0GeVto30GeV.pdf");
+
+ c8->cd();
+ h8->GetXaxis()->SetTitle("Hit position X [cm]");
+ h8->GetYaxis()->SetTitle("Hit position Y [cm]");
+ h8->GetYaxis()->SetTitleOffset(1.4);
+ h8->DrawClone("COLZ");
+ ell1->Draw("same");
+ ell2->Draw("same");
+ c8->SaveAs("results/pi0_hit_pos_cluster_0GeVto30GeV.png");
+ c8->SaveAs("results/pi0_hit_pos_cluster_0GeVto30GeV.pdf");
+
+ c9->cd();
+ h9->GetXaxis()->SetTitle("Hit position X [cm]");
+ h9->GetYaxis()->SetTitle("Hit position Y [cm]");
+ h9->GetYaxis()->SetTitleOffset(1.4);
+ h9->DrawClone("COLZ");
+ ell1->Draw("same");
+ ell2->Draw("same");
+ c9->SaveAs("results/pi0_hit_pos_all_0GeVto30GeV.png");
+ c9->SaveAs("results/pi0_hit_pos_all_0GeVto30GeV.pdf");
+
+ c10->cd();
+ h10->SetLineColor(kBlue);
+ h10->SetLineWidth(2);
+ h10->GetXaxis()->SetTitle("Invariant mass [MeV]");
+ h10->GetYaxis()->SetTitle("events");
+ h10->GetYaxis()->SetTitleOffset(1.4);
+ h10->DrawClone();
+ c10->SaveAs("results/pi0_invariant_mass_hist_0GeVto30GeV.png");
+ c10->SaveAs("results/pi0_invariant_mass_hist_0GeVto30GeV.pdf");
+
+ c11->cd();
+ h11->SetLineColor(kBlue);
+ h11->SetLineWidth(2);
+ h11->GetXaxis()->SetTitle("Cluster energy 1 [MeV]");
+ h11->GetYaxis()->SetTitle("events");
+ h11->GetYaxis()->SetTitleOffset(1.4);
+ h11->DrawClone();
+ c11->SaveAs("results/pi0_E1_hist_0GeVto30GeV.png");
+ c11->SaveAs("results/pi0_E1_hist_0GeVto30GeV.pdf");
+
+ c12->cd();
+ h12->SetLineColor(kBlue);
+ h12->SetLineWidth(2);
+ h12->GetXaxis()->SetTitle("Cluster energy 2 [MeV]");
+ h12->GetYaxis()->SetTitle("events");
+ h12->GetYaxis()->SetTitleOffset(1.4);
+ h12->DrawClone();
+ c12->SaveAs("results/pi0_E2_hist_0GeVto30GeV.png");
+ c12->SaveAs("results/pi0_E2_hist_0GeVto30GeV.pdf");
+
+ c13->cd();
+ h13->SetLineColor(kBlue);
+ h13->SetLineWidth(2);
+ h13->GetXaxis()->SetTitle("Angle between two photons [degree]");
+ h13->GetYaxis()->SetTitle("events");
+ h13->GetYaxis()->SetTitleOffset(1.4);
+ h13->DrawClone();
+ c13->SaveAs("results/pi0_angle_2photons_0GeVto30GeV.png");
+ c13->SaveAs("results/pi0_angle_2photons_0GeVto30GeV.pdf");
+
+ return 0;
+}
--
GitLab