From a8c6d875240d46ab38944f12388119d0d7a6b029 Mon Sep 17 00:00:00 2001
From: Peter Steinberg <steinber@gmail.com>
Date: Thu, 10 Oct 2024 13:20:31 -0400
Subject: [PATCH] first commit of LFHCAL benchmarks (#48)
Co-authored-by: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>
---
.gitlab-ci.yml | 1 +
Snakefile | 1 +
benchmarks/lfhcal/LFHCAL_Performance.C | 156 ++++++++++++++++++
benchmarks/lfhcal/README.md | 11 ++
benchmarks/lfhcal/Snakefile | 145 +++++++++++++++++
benchmarks/lfhcal/config.yml | 48 ++++++
benchmarks/lfhcal/doCompare_widebins_mom.C | 176 +++++++++++++++++++++
7 files changed, 538 insertions(+)
create mode 100644 benchmarks/lfhcal/LFHCAL_Performance.C
create mode 100644 benchmarks/lfhcal/README.md
create mode 100644 benchmarks/lfhcal/Snakefile
create mode 100644 benchmarks/lfhcal/config.yml
create mode 100644 benchmarks/lfhcal/doCompare_widebins_mom.C
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index a56c79d2..f25c74e9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -132,6 +132,7 @@ include:
- local: 'benchmarks/tracking_performances_dis/config.yml'
- local: 'benchmarks/barrel_ecal/config.yml'
- local: 'benchmarks/barrel_hcal/config.yml'
+ - local: 'benchmarks/lfhcal/config.yml'
- local: 'benchmarks/zdc/config.yml'
- local: 'benchmarks/zdc_lyso/config.yml'
- local: 'benchmarks/zdc_photon/config.yml'
diff --git a/Snakefile b/Snakefile
index 9fbb633a..1390e867 100644
--- a/Snakefile
+++ b/Snakefile
@@ -7,6 +7,7 @@ include: "benchmarks/ecal_gaps/Snakefile"
include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
include: "benchmarks/tracking_performances_dis/Snakefile"
+include: "benchmarks/lfhcal/Snakefile"
include: "benchmarks/insert_muon/Snakefile"
include: "benchmarks/zdc_lambda/Snakefile"
include: "benchmarks/zdc_lyso/Snakefile"
diff --git a/benchmarks/lfhcal/LFHCAL_Performance.C b/benchmarks/lfhcal/LFHCAL_Performance.C
new file mode 100644
index 00000000..dcfb7595
--- /dev/null
+++ b/benchmarks/lfhcal/LFHCAL_Performance.C
@@ -0,0 +1,156 @@
+// Code to extract the Tracking Performances
+// Shyam Kumar; INFN Bari, Italy
+// shyam.kumar@ba.infn.it; shyam.kumar@cern.ch
+
+#include "TGraphErrors.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#include "TVector3.h"
+
+#define mpi 0.139 // 1.864 GeV/c^2
+
+void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.0, TString name = "")
+{
+
+ // style of the plot
+ gStyle->SetPalette(1);
+ gStyle->SetOptTitle(1);
+ gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
+ gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
+ gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
+ gStyle->SetHistLineWidth(2);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptStat(1);
+
+ TString dir = "";
+ TString dist_dir_mom = "mom_resol";
+
+ bool debug=true;
+ // Tree with reconstructed tracks
+ const int nbins_eta = 5;
+ int nfiles = 100;
+ double eta[nbins_eta+1]={1.2,1.5,2,2.5,3,3.5};
+ double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1};
+ TH1D *histp[nbins_eta];
+
+
+ for (int i=0; i<nbins_eta; i++){
+ histp[i] = new TH1D(Form("hist_etabin%d",i),Form("hist_etabin%d",i),600,-1,1);
+ histp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom));
+ histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
+ }
+
+ TFile* file = TFile::Open(filename.Data());
+ if (!file) {printf("file not found !!!"); return;}
+ TTreeReader myReader("events", file); // name of tree and file
+ if (debug) cout<<"Filename: "<<file->GetName()<<"\t NEvents: "<<myReader.GetEntries()<<endl;
+
+ // MC and Reco information
+ TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge");
+ TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x");
+ TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y");
+ TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
+ TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x");
+ TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y");
+ TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z");
+ TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus");
+ TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");
+
+ TTreeReaderArray<Float_t> pe_lc(myReader, "LFHCALClusters.energy");
+ TTreeReaderArray<Float_t> px_lc(myReader, "LFHCALClusters.position.x");
+ TTreeReaderArray<Float_t> py_lc(myReader, "LFHCALClusters.position.y");
+ TTreeReaderArray<Float_t> pz_lc(myReader, "LFHCALClusters.position.z");
+ TTreeReaderArray<Float_t> pe_ec(myReader, "EcalEndcapPClusters.energy");
+ TTreeReaderArray<Float_t> px_ec(myReader, "EcalEndcapPClusters.position.x");
+ TTreeReaderArray<Float_t> py_ec(myReader, "EcalEndcapPClusters.position.y");
+ TTreeReaderArray<Float_t> pz_ec(myReader, "EcalEndcapPClusters.position.z");
+
+ int count =0;
+ int matchId = 1; // Always matched track assigned the index 0
+ while (myReader.Next())
+ {
+ std::cout << "events = " << count++ << std::endl;
+ for (int j = 0; j < pdg.GetSize(); ++j)
+ {
+ if (status[j] !=1 && pdg.GetSize()!=1) continue;
+ Double_t pzmc = pz_mc[j];
+ Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]);
+ Double_t pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec
+ Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2));
+ Double_t phimc = TMath::ATan2(py_mc[j],px_mc[j]);
+ std::cout << "neutron p=" << pmc << " pt=" << ptmc << std::endl;
+
+ if (fabs(ptmc) < pTcut) continue;
+
+ float l_px_tot=0;
+ float l_py_tot=0;
+ float l_pz_tot=0;
+ float l_e_tot=0;
+
+ std::cout << "LFHCAL nclus=" << px_lc.GetSize() << " ECAL nclus=" << pe_ec.GetSize() << std::endl;
+ for (int jl = 0;jl<px_lc.GetSize();jl++)
+ {
+ float e = pe_lc[jl];
+ TVector3 v(px_lc[jl],py_lc[jl],pz_lc[jl]);
+ v.Print();
+ float eta = v.PseudoRapidity();
+ float phi = v.Phi();
+ float pt = e/cosh(eta);
+ std::cout << "LFHCAL clus: e=" << e << " eta=" << eta << " pt=" << pt << std::endl;
+ l_e_tot += e;
+ l_px_tot += pt*cos(phi);
+ l_py_tot += pt*sin(phi);
+ l_pz_tot += pt*sinh(eta);
+ }
+
+ float e_px_tot=0;
+ float e_py_tot=0;
+ float e_pz_tot=0;
+ float e_e_tot=0;
+
+ for (int je = 0;je<px_ec.GetSize();je++)
+ {
+ float e = pe_ec[je];
+ TVector3 v(px_ec[je],py_ec[je],pz_ec[je]);
+ float eta = v.PseudoRapidity();
+ float phi = v.Phi();
+ float pt = e/cosh(eta);
+ std::cout << "ECAL clus: e=" << e << " eta=" << eta << " pt=" << pt << std::endl;
+ e_e_tot += e;
+ e_px_tot += pt*cos(phi);
+ e_py_tot += pt*sin(phi);
+ e_pz_tot += pt*sinh(eta);
+ }
+
+ std::cout << "LFHCAL e=" <<l_e_tot << " ECAL e=" << e_e_tot << std::endl;
+ float px_tot = l_px_tot+e_px_tot;
+ float py_tot = l_py_tot+e_py_tot;
+ float pz_tot = l_pz_tot+e_pz_tot;
+ float e_tot = l_e_tot+e_e_tot;
+
+ float prec = sqrt(px_tot*px_tot+py_tot*py_tot+pz_tot*pz_tot);
+ float ptrec = sqrt(px_tot*px_tot+py_tot*py_tot);
+ float pzrec = pz_tot;
+
+ float p_resol = (e_tot-pmc)/pmc;
+ std::cout << "p_resol = " << p_resol << std::endl;
+ for (int ibin=0; ibin<nbins_eta; ++ibin){
+ if(etamc>eta[ibin] && etamc<eta[ibin+1]) histp[ibin]->Fill(p_resol);
+ }
+ } // Generated Tracks
+
+ }// event loop ends
+
+ TFile *fout_mom = new TFile(Form("%s/mom/lfhcal_mom_%1.1f_%s_%s.root",particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate");
+ fout_mom->cd();
+ for (int ibin=0; ibin<nbins_eta; ++ibin) histp[ibin]->Write();
+ fout_mom->Close();
+
+}
+
+
+
+
diff --git a/benchmarks/lfhcal/README.md b/benchmarks/lfhcal/README.md
new file mode 100644
index 00000000..5f56298b
--- /dev/null
+++ b/benchmarks/lfhcal/README.md
@@ -0,0 +1,11 @@
+Detector Benchmark for LFHCAL
+===============================
+
+## Overview
+This benchmark generates events with a range of particle species at fixed energies, all aimed at, or forward of, the LFHCAL, and produces basic plots of clustering performance, both in isolation and combined with the FEMCAL.
+
+
+## Contacts
+[Peter Steinberg](steinberg@bnl.gov)
+
+
diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile
new file mode 100644
index 00000000..459392c0
--- /dev/null
+++ b/benchmarks/lfhcal/Snakefile
@@ -0,0 +1,145 @@
+rule lfhcal_sim:
+ input:
+ steering_file="EPIC/EVGEN/SINGLE/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.steer",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+ output:
+ "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
+ log:
+ "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log",
+ wildcard_constraints:
+ ENERGY="[0-9]+[kMG]eV",
+ PARTICLE="(neutron|pi-|gamma)",
+ PHASE_SPACE="3to50deg",
+ INDEX="\d{4}",
+ params:
+ N_EVENTS=1000
+ shell:
+ """
+ddsim \
+ --runType batch \
+ --enableGun \
+ --steeringFile "{input.steering_file}" \
+ --random.seed 1{wildcards.INDEX} \
+ --filter.tracker edep0 \
+ -v INFO \
+ --numberOfEvents {params.N_EVENTS} \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --outputFile {output}
+"""
+
+
+rule lfhcal_recon:
+ input:
+ "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
+ output:
+ "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root",
+ log:
+ "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log",
+ wildcard_constraints:
+ INDEX="\d{4}",
+ shell: """
+env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
+ eicrecon {input} -Ppodio:output_file={output} \
+ -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits
+"""
+
+rule lfhcal_at_momentum:
+ input:
+ script="benchmarks/lfhcal/LFHCAL_Performance.C",
+ # TODO pass as a file list?
+ sim=lambda wildcards:
+ expand(
+ "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
+ DETECTOR_CONFIG="epic_craterlake",
+ ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV",
+ PHASE_SPACE=["3to50deg"],
+ INDEX=range(1),
+ )
+ if wildcards.CAMPAIGN == "local" else
+ ancient(expand(
+ "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
+ DETECTOR_CONFIG="epic_craterlake",
+ ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV",
+ PHASE_SPACE=["3to50deg"],
+ INDEX=range(1),
+ )),
+ output:
+ "{CAMPAIGN}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root",
+ combined_root=temp("{CAMPAIGN}/lfhcal_sim_{MOMENTUM}_{PARTICLE}.root"),
+ shell:
+ """
+hadd {output.combined_root} {input.sim}
+cd {wildcards.CAMPAIGN}
+root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15)'
+"""
+
+rule lfhcal_summary_at_eta:
+ input:
+ expand(
+ [
+ "{{CAMPAIGN}}/{{PARTICLE}}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{{PARTICLE}}.root",
+ ],
+ MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0],
+ ),
+ script="benchmarks/lfhcal/doCompare_widebins_mom.C",
+ output:
+ "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.png",
+ "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.root",
+ shell:
+ """
+if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then
+ set +e
+ EPIC_VERSION="${{DETECTOR_VERSION:-}}"
+ EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9\.]\+\).*/\\1/p')"
+ # Legacy detection
+ : ${{EPIC_VERSION:="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-/]\+\).*/\\1/p')"}}
+ set -e
+
+ echo "ePIC version: $EPIC_VERSION"
+ echo "EICrecon version: $EICRECON_VERSION"
+ EXTRA_LEGEND="ePIC $EPIC_VERSION / EICrecon $EICRECON_VERSION"
+else
+ EXTRA_LEGEND="ePIC Simulation {wildcards.CAMPAIGN}"
+fi
+cd {wildcards.CAMPAIGN}
+root -l -b -q ../{input.script}'("{wildcards.PARTICLE}", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")'
+"""
+
+LFHCAL_ETA_BINS = [1.2,1.5,2,2.5,3,3.5]
+
+rule lfhcal:
+ input:
+ lambda wildcards: expand(
+ [
+ "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.png",
+ "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.root",
+ ],
+ ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])],
+ PARTICLE=["neutron", "pi-", "gamma"] if wildcards.CAMPAIGN == "local" else ["pi-"],
+ )
+ output:
+ directory("results/lfhcal/{CAMPAIGN}/")
+ shell:
+ """
+mkdir {output}
+cp {input} {output}
+"""
+
+
+rule lfhcal_local:
+ input:
+ "results/lfhcal/local",
+
+
+rule lfhcal_campaigns:
+ input:
+ expand(
+ "results/lfhcal/{CAMPAIGN}",
+ CAMPAIGN=[
+ "23.10.0",
+ "23.11.0",
+ "23.12.0",
+ "24.03.1",
+ "24.04.0",
+ ],
+ )
diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml
new file mode 100644
index 00000000..f7b9ad3a
--- /dev/null
+++ b/benchmarks/lfhcal/config.yml
@@ -0,0 +1,48 @@
+sim:lfhcal:
+ extends: .det_benchmark
+ stage: simulate
+ parallel:
+ matrix:
+ - PARTICLE: ["neutron", "pi-"]
+ MOMENTUM: ["500MeV", "1GeV", "2GeV", "5GeV", "10GeV", "20GeV"]
+ script:
+ - |
+ snakemake --cache --cores 5 \
+ sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/3to50deg/${PARTICLE}_${MOMENTUM}_3to50deg.0000.eicrecon.tree.edm4eic.root
+
+bench:lfhcal:
+ extends: .det_benchmark
+ stage: benchmarks
+ needs:
+ - ["sim:lfhcal"]
+ script:
+ - snakemake --cores 3 lfhcal_local
+
+collect_results:lfhcal:
+ extends: .det_benchmark
+ stage: collect
+ needs:
+ - "bench:lfhcal"
+ script:
+ - ls -lrht
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output lfhcal_local
+ - mv results{_save,}/
+
+bench:lfhcal_campaigns:
+ extends: .det_benchmark
+ stage: benchmarks
+ #when: manual
+ script:
+ - snakemake --cores 1 lfhcal_campaigns
+
+collect_results:lfhcal_campaigns:
+ extends: .det_benchmark
+ stage: collect
+ needs:
+ - "bench:lfhcal_campaigns"
+ script:
+ - ls -lrht
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output lfhcal_campaigns
+ - mv results{_save,}/
diff --git a/benchmarks/lfhcal/doCompare_widebins_mom.C b/benchmarks/lfhcal/doCompare_widebins_mom.C
new file mode 100644
index 00000000..5258e5ef
--- /dev/null
+++ b/benchmarks/lfhcal/doCompare_widebins_mom.C
@@ -0,0 +1,176 @@
+// Code to compare the tracking performances: Truth seeding vs real seeding
+// Shyam Kumar; shyam.kumar@ba.infn.it; shyam055119@gmail.com
+
+#include "TGraphErrors.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#define mpi 0.139 // 1.864 GeV/c^2
+
+//void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.);
+void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, double range =0.3, Bool_t drawreq=1, TString extra_legend = "") // name = p, pt for getting p or pt dependence fitted results
+{
+
+ //=== style of the plot=========
+ gStyle->SetPalette(1);
+ gStyle->SetOptTitle(1);
+ gStyle->SetTitleOffset(1.0,"XY");
+ gStyle->SetTitleSize(.04,"XY");
+ gStyle->SetLabelSize(.04,"XY");
+ gStyle->SetHistLineWidth(2);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptStat(1);
+
+ const Int_t nfiles = 6;
+ double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,20.0};
+ std::vector<double> momV, momresolV, err_momresolV;
+ momV.clear(); momresolV.clear(); err_momresolV.clear();
+ TString symbolname = "";
+ if (particle == "pi-") symbolname = "#pi^{-}";
+ else symbolname = particle;
+ ofstream outfile;
+ outfile.open ("Mom_resol.txt",ios_base::app);
+
+ TF1 *f1=new TF1("f1","FitMomentumResolution",0.,30.0,2);
+ f1->SetParLimits(0,0.,0.1);
+ f1->SetParLimits(1,0.,5.0);
+
+ TCanvas *c_mom = new TCanvas("cmom","cmom",1400,1000);
+ c_mom->SetMargin(0.10, 0.05 ,0.1,0.05);
+ c_mom->SetGridy();
+
+ //Reading the root file
+ TFile *fmom[nfiles];
+ TGraphErrors *gr_mom;
+ TMultiGraph *mgMom;
+ TLegend *lmom;
+ mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %");
+
+ lmom = new TLegend(0.65,0.80,0.90,0.93);
+ lmom->SetTextSize(0.03);
+ lmom->SetBorderSize(0);
+ lmom->SetHeader(extra_legend.Data(), "C");
+ lmom->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.Data(), etamin, etamax), "C");
+
+ TF1 *func = new TF1("func","gaus",-0.5,0.5);
+
+ for (int i =0; i<nfiles; ++i){
+
+ TCanvas *cp = new TCanvas("cp","cp",1400,1000);
+ cp->SetMargin(0.10, 0.05 ,0.1,0.07);
+
+ //pi-/mom/lfhcal_mom_20.0_mom_resol_pi-.root
+ fmom[i] = TFile::Open(Form("./%s/mom/lfhcal_mom_%1.1f_mom_resol_%s.root",particle.Data(),mom[i],particle.Data()));
+
+ TH1D *hist = (TH1D*) fmom[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax));
+ hist->Rebin(2);
+ hist->SetName(Form("hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data()));
+ hist->SetTitle(Form("Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax));
+
+ double mu = hist->GetMean();
+ double sigma = hist->GetRMS();
+ double sigma_err = hist->GetRMSError();
+ /*
+ hist->GetXaxis()->SetRangeUser(-1.0*range,1.0*range);
+ func->SetRange(mu-2.0*sigma,mu+2.0*sigma); // fit with in 2 sigma range
+ hist->Fit(func,"NR+");
+ mu = func->GetParameter(1);
+ sigma = func->GetParameter(2);
+ func->SetRange(mu-2.0*sigma,mu+2.0*sigma);
+ hist->Fit(func,"R+");
+ float par2 = func->GetParameter(2)*100;
+ float par2_err = func->GetParError(2)*100;
+ */
+ momV.push_back(mom[i]);
+ momresolV.push_back(sigma);
+ err_momresolV.push_back(sigma_err);
+
+ cp->cd();
+ hist->Draw();
+ //cp->SaveAs(Form("Debug_Plots/%s/mom/mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax));
+ } // all files
+
+ const int size = momV.size();
+ double p[size], err_p[size], sigma_p[size], err_sigma_p[size];
+
+ for (int i=0; i<size; i++){
+ p[i] = momV.at(i);
+ sigma_p[i] = momresolV.at(i);
+ err_sigma_p[i] = err_momresolV.at(i);
+ err_p[i] = 0.;
+ }
+
+
+ TFile *fout = new TFile(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.root",particle.Data(),particle.Data(),etamin,etamax),"recreate");
+ TGraphErrors *gr1 = new TGraphErrors(size,p,sigma_p,err_p,err_sigma_p);
+ gr1->SetName("grseed");
+ gr1->SetMarkerStyle(25);
+ gr1->SetMarkerColor(kBlue);
+ gr1->SetMarkerSize(2.0);
+ gr1->SetTitle(";p (GeV/c);#sigmap/p");
+ gr1->GetXaxis()->CenterTitle();
+ gr1->GetYaxis()->CenterTitle();
+
+
+ mgMom->Add(gr1);
+ c_mom->cd();
+ mgMom->GetXaxis()->SetRangeUser(0.40,20.2);
+ mgMom->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr1->GetN(),gr1->GetY())); // 50% more of the maximum value on yaxis
+ mgMom->Draw("AP");
+ lmom->AddEntry(gr1,"Nominal");
+ lmom->Draw("same");
+ //draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax());
+ c_mom->SaveAs(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.png",particle.Data(),particle.Data(),etamin,etamax));
+
+ // Write the numbers in output file for comparisons
+ outfile << extra_legend << endl;
+ outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"p (GeV/c) \t"<<setw(20)<<"Resol #mum "<<endl;
+ for (Int_t i = 0; i<gr1->GetN(); ++i){
+ double x,y;
+ gr1->GetPoint(i,x,y);
+ outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<y<<endl;
+ }
+ outfile.close();
+
+ fout->cd();
+ mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax));
+ mgMom->Write();
+ fout->Close();
+}
+
+//===Fit Momentum Resolution
+float FitMomentumResolution(Double_t *x, Double_t *par)
+{
+ float func = sqrt(par[0]*par[0]*x[0]*x[0]+par[1]*par[1]);
+ return func;
+}
+
+//From Yellow report from section 11.2.2
+
+// 1.2,1.5,2,2.5,3,3.5
+
+/*
+void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.)
+{
+
+ TF1 *dd4hep_p;
+ if (etamin >= -3.5 && etamax <= -2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
+ else if (etamin >= -2.5 && etamax <= -1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax);
+ else if (etamin >= -1.0 && etamax <= 1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+0.5^2)",xmin,xmax);
+ else if (etamin >= 1.0 && etamax <= 2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax);
+ else if (etamin >= 2.5 && etamax <= 3.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
+ else return;
+ dd4hep_p->SetLineStyle(7);
+ dd4hep_p->SetLineColor(kMagenta);
+ dd4hep_p->SetLineWidth(3.0);
+ dd4hep_p->Draw("same");
+
+ TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
+ l->SetTextSize(0.03);
+ l->SetBorderSize(0);
+ l->AddEntry(dd4hep_p,"PWGReq","l");
+ l->Draw("same");
+ }
+*/
--
GitLab