diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index a56c79d..f25c74e 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 9fbb633..1390e86 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 0000000..dcfb759
--- /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 0000000..5f56298
--- /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 0000000..459392c
--- /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 0000000..f7b9ad3
--- /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 0000000..5258e5e
--- /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");
+ }
+*/