Adding benchmarks for Backward hcal (#147)

Co-authored-by: lkosarz <42405623+lkosarz@users.noreply.github.com>
Co-authored-by: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>

.gitlab-ci.yml

@@ -149,6 +149,8 @@ include:
   - local: 'benchmarks/femc_electron/config.yml'
   - local: 'benchmarks/femc_photon/config.yml'
   - local: 'benchmarks/femc_pi0/config.yml'
+  - local: 'benchmarks/nhcal_acceptance/config.yml'
+  - local: 'benchmarks/nhcal_basic_distribution/config.yml'
   
 deploy_results:
   allow_failure: true
@@ -180,6 +182,8 @@ deploy_results:
     - "collect_results:femc_electron"
     - "collect_results:femc_photon"
     - "collect_results:femc_pi0"
+    - "collect_results:nhcal_acceptance"
+    - "collect_results:nhcal_basic_distribution"
   script:
     - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json
     - find results -print | sort | tee summary.txt

Snakefile

@@ -52,6 +52,8 @@ include: "benchmarks/insert_tau/Snakefile"
 include: "benchmarks/femc_electron/Snakefile"
 include: "benchmarks/femc_photon/Snakefile"
 include: "benchmarks/femc_pi0/Snakefile"
+include: "benchmarks/nhcal_acceptance/Snakefile"
+include: "benchmarks/nhcal_basic_distribution/Snakefile"
 
 use_s3 = config["remote_provider"].lower() == "s3"
 use_xrootd = config["remote_provider"].lower() == "xrootd"

benchmarks/nhcal_acceptance/Snakefile

+rule nhcal_acceptance_simulate:
+    input:
+        warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+    output:
+        "sim_output/nhcal_acceptance/E{ENERGY}GeV/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root",
+    params:
+        N_EVENTS=10000,
+        SEED=lambda wildcards: "1" + wildcards.INDEX,
+        DETECTOR_PATH=os.environ["DETECTOR_PATH"],
+        DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
+        ENERGY=lambda wildcards: wildcards.ENERGY,
+        DD4HEP_HASH=get_spack_package_hash("dd4hep"),
+        NPSIM_HASH=get_spack_package_hash("npsim"),
+    cache: True
+    shell:
+        """
+exec ddsim \
+    --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
+    --numberOfEvents {params.N_EVENTS} \
+    --random.seed {params.SEED} \
+    --enableGun \
+    -v WARNING \
+    --gun.particle pi- \
+    --gun.thetaMin 120*degree \
+    --gun.thetaMax 180*degree \
+    --gun.distribution uniform \
+    --gun.energy "{params.ENERGY}*GeV" \
+    --outputFile {output}
+"""
+
+
+rule nhcal_acceptance_combine:
+    input:
+        lambda wildcards: expand(
+            "sim_output/nhcal_acceptance/E{ENERGY}GeV/sim_{DETECTOR_CONFIG}.{INDEX:02d}.edm4hep.root", 
+            DETECTOR_CONFIG=wildcards.DETECTOR_CONFIG,
+            ENERGY=wildcards.ENERGY,
+            INDEX=range(int(wildcards.N)),
+        ),
+    wildcard_constraints:
+        N=r"\d+",
+        ENERGY=r"\d+"
+    output:
+        temp("sim_output/nhcal_acceptance/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root"),
+    shell:
+        """
+hadd -f {output} {input} 
+"""
+
+rule nhcal_acceptance_analysis:
+    input:
+        combined="sim_output/nhcal_acceptance/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root",
+        script="benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx",
+    output:
+        png=f"results/nhcal_acceptance/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.png",
+        pdf=f"results/nhcal_acceptance/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.pdf",
+    shell:
+        """
+    root -l -b -q '{input.script}("{input.combined}","{output.pdf}","{output.png}")'
+"""

benchmarks/nhcal_acceptance/config.yml

+sim:nhcal_acceptance:
+  extends: .det_benchmark
+  stage: simulate
+  parallel:
+    matrix:
+      - ENERGY: ["1GeV", "5GeV", "10GeV"]
+        INDEX_RANGE: ["0 4","5 9"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 5 $(for INDEX in $(seq -f '%02.0f' $INDEX_RANGE); do echo sim_output/nhcal_acceptance/E${ENERGY}/sim_epic_backward_hcal_only.${INDEX}.edm4hep.root; done)
+
+bench:nhcal_acceptance_analysis:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - "sim:nhcal_acceptance"
+  parallel:
+    matrix:
+      - ENERGY: ["1GeV", "5GeV", "10GeV"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/nhcal_acceptance/analysis_epic_backward_hcal_only_E${ENERGY}_combined_10files.pdf
+
+collect_results:nhcal_acceptance:
+  extends: .det_benchmark
+  stage: collect
+  needs: 
+    - "bench:nhcal_acceptance_analysis"
+  parallel:
+    matrix:
+      - ENERGY: ["1GeV", "5GeV", "10GeV"]
+  script:
+    - ls -lrht
+    - mv results{,_save}/ 
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_acceptance/analysis_epic_backward_hcal_only_E${ENERGY}_combined_10files.pdf
+    - mv results{_save,}/

benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx

+#include <TChain.h>
+#include <TTreeReader.h>
+#include <TTreeReaderArray.h>
+#include <TTreeReaderValue.h>
+#include <TH2D.h>
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <set>
+#include <iostream>
+#include <TString.h>
+
+using namespace std;
+
+int acceptance_analysis(TString filename, string outname_pdf, string outname_png) 
+{
+    TChain *chain = new TChain("events");
+    chain->Add(filename);
+    
+    TTreeReader reader(chain);
+    
+    TTreeReaderArray<int> mc_pdg(reader, "MCParticles.PDG");
+    TTreeReaderArray<int> mc_genStatus(reader, "MCParticles.generatorStatus");
+    TTreeReaderArray<double> mc_px(reader, "MCParticles.momentum.x");
+    TTreeReaderArray<double> mc_py(reader, "MCParticles.momentum.y");
+    TTreeReaderArray<double> mc_pz(reader, "MCParticles.momentum.z");
+
+    TTreeReaderArray<int> contrib_particle_idx(reader, "_HcalEndcapNHitsContributions_particle.index");
+    TTreeReaderArray<unsigned int> contrib_particle_cid(reader, "_HcalEndcapNHitsContributions_particle.collectionID");
+
+    int nEtaBins = 100;
+    int nPhiBins = 100;
+    double etaMin = -5, etaMax = 0;
+
+    TH2D* hEtaPhiAll = new TH2D("hEtaPhiAll", "All #pi- (status==1); #eta[1]; #phi[rad]",
+                                nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi());
+
+    TH2D* hEtaPhiDetected = new TH2D("hEtaPhiDetected", "#pi- detected in nHCal; #eta[1]; #phi[rad]",
+                                     nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi());
+
+    while (reader.Next()) 
+    {
+        map<int, pair<double, double>> pi_minus_eta_phi; 
+        set<int> detected;
+        for (size_t i = 0; i < mc_pdg.GetSize(); ++i) 
+        {
+            if (mc_pdg[i] == -211 && mc_genStatus[i] == 1) 
+            {
+                float px = mc_px[i];
+                float py = mc_py[i];
+                float pz = mc_pz[i];
+
+                float p = sqrt(px * px + py * py + pz * pz);
+                float eta = 0.5 * log((p + pz) / (p - pz + 1e-8));
+                float phi = atan2(py, px);
+
+                hEtaPhiAll->Fill(eta, phi);
+                pi_minus_eta_phi[i] = make_pair(eta, phi);
+            }
+        }
+
+        for (size_t i = 0; i < contrib_particle_idx.GetSize(); i++) {
+            int idx = contrib_particle_idx[i];
+            if (pi_minus_eta_phi.count(idx)) {
+                detected.insert(idx);
+            }
+        }
+
+        for (auto idx : detected) {
+            auto [eta, phi] = pi_minus_eta_phi[idx];
+            hEtaPhiDetected->Fill(eta, phi);
+        }
+    }
+
+    TH2D* hAcceptance = (TH2D*)hEtaPhiAll->Clone("hAcceptance");
+    hAcceptance->Divide(hEtaPhiDetected);
+    hAcceptance->SetTitle("#pi- detected/All");
+    hAcceptance->SetMinimum(0);
+    hAcceptance->SetMaximum(1);
+
+    TCanvas *canvas = new TCanvas("canvas", "pi- All", 1600, 600);
+    canvas->Divide(3,1);
+    canvas->cd(1);
+    hEtaPhiAll->Draw("COLZ");
+    canvas->cd(2);
+    hEtaPhiDetected->Draw("COLZ");
+    canvas->cd(3);
+    hAcceptance->Draw("COLZ");
+
+    canvas->SaveAs(outname_pdf.c_str());
+    canvas->SaveAs(outname_png.c_str());
+
+    return 0;
+}
+
+

benchmarks/nhcal_basic_distribution/Snakefile

+rule nhcal_basic_distribution_simulate:
+    input:
+        warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+    output:
+        "sim_output/nhcal_basic_distribution/E{ENERGY}GeV/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root",
+    params:
+        N_EVENTS=10000,
+        SEED=lambda wildcards: "1" + wildcards.INDEX,
+        DETECTOR_PATH=os.environ["DETECTOR_PATH"],
+        DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
+        ENERGY=lambda wildcards: wildcards.ENERGY,
+        DD4HEP_HASH=get_spack_package_hash("dd4hep"),
+        NPSIM_HASH=get_spack_package_hash("npsim"),
+    cache: True
+    shell:
+        """
+exec ddsim \
+    --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
+    --numberOfEvents {params.N_EVENTS} \
+    --random.seed {params.SEED} \
+    --enableGun \
+    -v WARNING \
+    --gun.particle neutron \
+    --gun.thetaMin 120*degree \
+    --gun.thetaMax 180*degree \
+    --gun.distribution uniform \
+    --gun.energy "{params.ENERGY}*GeV" \
+    --outputFile {output}
+"""
+
+
+rule nhcal_basic_distribution_combine:
+    input:
+        lambda wildcards: expand(
+            "sim_output/nhcal_basic_distribution/E{ENERGY:.1f}GeV/sim_{DETECTOR_CONFIG}.{INDEX:02d}.edm4hep.root", 
+            DETECTOR_CONFIG=wildcards.DETECTOR_CONFIG,
+            ENERGY=float(wildcards.ENERGY),
+            INDEX=range(int(wildcards.N)),
+        ),
+    wildcard_constraints:
+        N=r"\d+",
+        ENERGY=r"\d+(\.\d+)?"
+    output:
+        temp("sim_output/nhcal_basic_distribution/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root"),
+    shell:
+        """
+hadd -f {output} {input} 
+"""
+
+rule nhcal_basic_distribution_analysis:
+    input:
+        combined="sim_output/nhcal_basic_distribution/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root",
+        script="benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx",
+    output:
+        pdf=f"results/nhcal_basic_distribution/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.pdf",
+        png=f"results/nhcal_basic_distribution/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.png",
+    shell:
+        """
+    root -l -b -q '{input.script}("{input.combined}","{output.pdf}","{output.png}")'
+"""

benchmarks/nhcal_basic_distribution/config.yml

+sim:nhcal_basic_distribution:
+  extends: .det_benchmark
+  stage: simulate
+  parallel:
+    matrix:
+      - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"]
+        INDEX_RANGE: ["0 4","5 9"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 5 $(for INDEX in $(seq -f '%02.0f' $INDEX_RANGE); do echo sim_output/nhcal_basic_distribution/E${ENERGY}/sim_epic_backward_hcal_only.${INDEX}.edm4hep.root; done)
+
+sim:nhcal_basic_distribution_full:
+  extends: .det_benchmark
+  stage: simulate
+  parallel:
+    matrix:
+      - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"]
+        INDEX_RANGE: ["0 4","5 9"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 5 $(for INDEX in $(seq -f '%02.0f' $INDEX_RANGE); do echo sim_output/nhcal_basic_distribution/E${ENERGY}/sim_epic_full.${INDEX}.edm4hep.root; done)
+
+bench:nhcal_basic_distribution_analysis:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - "sim:nhcal_basic_distribution"
+  parallel:
+    matrix:
+      - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/nhcal_basic_distribution/analysis_epic_backward_hcal_only_E${ENERGY}_combined_10files.pdf
+
+bench:nhcal_basic_distribution_analysis_full:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - "sim:nhcal_basic_distribution_full"
+  parallel:
+    matrix:
+      - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"]
+  script:
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/nhcal_basic_distribution/analysis_epic_full_E${ENERGY}_combined_10files.pdf
+
+
+collect_results:nhcal_basic_distribution:
+  extends: .det_benchmark
+  stage: collect
+  needs: 
+    - "bench:nhcal_basic_distribution_analysis"
+    - "bench:nhcal_basic_distribution_analysis_full"
+  parallel:
+    matrix:
+      - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"]
+        DETECTOR_CONFIG: ["epic_backward_hcal_only", "epic_full"]
+  script:
+    - ls -lrht
+    - mv results{,_save}/ 
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_basic_distribution/analysis_${DETECTOR_CONFIG}_E${ENERGY}_combined_10files.pdf
+    - mv results{_save,}/
+        

benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx

+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "ROOT/RDataFrame.hxx"
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TRandom3.h>
+#include <TFile.h>
+#include <TChain.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TCanvas.h>
+
+#include "TROOT.h"
+#include "TRandom.h"
+#include "TH3.h"
+
+
+#include "DD4hep/Detector.h"
+#include "DDRec/CellIDPositionConverter.h"
+
+#include <podio/Frame.h>
+#include <podio/CollectionBase.h>
+#include "podio/ROOTReader.h"
+//#include <podio/ROOTFrameReader.h>
+#include "podio/CollectionIDTable.h"
+#include "podio/ObjectID.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/MCParticleCollectionData.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/MCParticleData.h"
+
+#include "edm4hep/SimCalorimeterHitCollectionData.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/SimCalorimeterHitData.h"
+#include "edm4hep/SimCalorimeterHit.h"
+
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/CalorimeterHitCollectionData.h"
+#include "edm4hep/CalorimeterHitCollection.h"
+#include "edm4hep/CalorimeterHitData.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/CalorimeterHitObj.h"
+
+#include "edm4eic/ClusterCollection.h"
+#include "edm4eic/Cluster.h"
+#include "edm4eic/ClusterData.h"
+
+#include "edm4eic/CalorimeterHit.h"
+#include "edm4eic/CalorimeterHitCollectionData.h"
+#include "edm4eic/CalorimeterHitCollection.h"
+#include "edm4eic/CalorimeterHitData.h"
+#include "edm4eic/CalorimeterHit.h"
+#include "edm4eic/CalorimeterHitObj.h"
+
+
+#include <edm4eic/vector_utils_legacy.h>
+#include <edm4hep/Vector3f.h>
+
+// #include "FileList.h"
+// #include "EICutil.h"
+// #include "BasicUtil.h"
+
+// #pragma link C++ class vector<edm4hep::MCParticleData>+;
+// #pragma link C++ class vector<podio::ObjectID>+;
+// #pragma link C++ class vector<edm4hep::SimCalorimeterHitData>+;
+
+using namespace std;
+using namespace ROOT;
+using namespace TMath;
+using namespace edm4hep;
+
+int basic_distribution_analysis(const string &filename, string outname_pdf, string outname_png) 
+{
+    podio::ROOTReader *reader = new podio::ROOTReader();
+    reader->openFile(filename);
+    unsigned nEvents = reader->getEntries("events");
+    cout << "Number of events: " << nEvents << endl;
+
+    vector<double> xPosition;
+    vector<double> yPosition;
+    vector<double> zPosition;
+    vector<double> rPosition;
+    vector<double> energy;
+    vector<double> totalEventEnergy;
+    vector<double> totalEventHits;
+    vector<double> zLayerValues;
+    vector<int> nHitsPerLayer;
+    vector<double> nEnergyPerLayer;
+    
+    for (unsigned ev = 0; ev < nEvents; ev++) 
+    {
+        auto frameData = reader->readNextEntry(podio::Category::Event);
+        if (!frameData) 
+        {
+            cerr << "Invalid FrameData at event " << ev << endl;
+            continue;
+        }
+
+        podio::Frame frame(std::move(frameData));
+
+        auto& hits = frame.get<edm4hep::SimCalorimeterHitCollection>("HcalEndcapNHits");
+        if (!hits.isValid())    
+        {
+            cerr << "HcalEndcapNHits collection is not valid!" << endl;
+        }
+
+        map<double, pair<int, double>> layerData;
+        double totalEnergy = 0;
+        int totalHits = 0;
+
+        for (const auto& hit : hits) 
+        {
+            totalHits++;
+            energy.push_back(hit.getEnergy());
+            totalEnergy += hit.getEnergy();
+
+            auto pos = hit.getPosition();
+            xPosition.push_back(pos.x);
+            yPosition.push_back(pos.y);
+            zPosition.push_back(pos.z);
+            rPosition.push_back(sqrt(pos.x * pos.x + pos.y * pos.y));
+
+            double zBin = round(pos.z);
+            layerData[zBin].first++;         
+            layerData[zBin].second += hit.getEnergy();
+        }
+
+        totalEventEnergy.push_back(totalEnergy);
+        totalEventHits.push_back(totalHits);
+
+        for (const auto& [zValue, stats] : layerData)
+        {
+            zLayerValues.push_back(zValue);
+            nHitsPerLayer.push_back(stats.first);
+            nEnergyPerLayer.push_back(stats.second);
+        }
+
+    }
+
+    auto minmax_xPosition        = minmax_element(xPosition.begin(), xPosition.end());
+    auto minmax_yPosition        = minmax_element(yPosition.begin(), yPosition.end());
+    auto minmax_zPosition        = minmax_element(zPosition.begin(), zPosition.end());
+    auto minmax_rPosition        = minmax_element(rPosition.begin(), rPosition.end());
+    auto minmax_totalEventEnergy = minmax_element(totalEventEnergy.begin(), totalEventEnergy.end());
+    auto minmax_totalEventHits   = minmax_element(totalEventHits.begin(), totalEventHits.end());
+    auto minmax_energy           = minmax_element(energy.begin(), energy.end());
+    auto minmax_zLayerValues     = minmax_element(zLayerValues.begin(), zLayerValues.end());
+    auto minmax_nHitsPerLayer    = minmax_element(nHitsPerLayer.begin(), nHitsPerLayer.end());
+    auto minmax_nEnergyPerLayer  = minmax_element(nEnergyPerLayer.begin(), nEnergyPerLayer.end());
+
+    int nbins = nEvents/1000;
+
+    TH1D *h_energyTotal = new TH1D("h_energyTotal", "Total Energy per Event; Energy per Event", 
+                                    nbins/2, *minmax_totalEventEnergy.first, *minmax_totalEventEnergy.second);
+    TH2D *h_layerEnergy = new TH2D("h_layerEnergy", "Energy in Layers (Z); Z; Energy", 
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second, 
+                                    nbins/4, *minmax_nEnergyPerLayer.first, *minmax_nEnergyPerLayer.second);
+    TProfile *p_layerEnergy = new TProfile("p_layerEnergy", "Energy in Layers (Z); Z; Mean Energy",
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second);
+    TH1D *h_hitCount    = new TH1D("h_hitCount", "Number of Hits per Event; Hits per Event", 
+                                    nbins/2, *minmax_totalEventHits.first, *minmax_totalEventHits.second);
+    TH2D *h_layerHits   = new TH2D("h_layerHits", "Number of Hits in Layers (Z); Layer(Z); Hits", 
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second, 
+                                    *minmax_nHitsPerLayer.second, *minmax_nHitsPerLayer.first, *minmax_nHitsPerLayer.second);
+    TProfile *p_layerHits = new TProfile("p_layerHits", "Number of Hits in Layers (Z); Z; Mean Hits",
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second);
+    TH2D *h_XYPos       = new TH2D("h_XYPos", "Hits position X,Y;X [mm];Y [mm]", 
+                                    nbins, *minmax_xPosition.first, *minmax_xPosition.second, 
+                                    nbins, *minmax_yPosition.first, *minmax_yPosition.second);
+    TH2D *h_ZRPos       = new TH2D("h_ZRPos", "Hits position Z,R;Z [mm];R [mm]", 
+                                    nbins/5, *minmax_zPosition.first, *minmax_zPosition.second, 
+                                    nbins, *minmax_rPosition.first, *minmax_rPosition.second);
+    TH2D *h_XYEnergy    = new TH2D("h_XYEnergy", "Hits energy X,Y;X [mm];Y [mm]", 
+                                    nbins, *minmax_xPosition.first, *minmax_xPosition.second, 
+                                    nbins, *minmax_yPosition.first, *minmax_yPosition.second);
+
+    for(int i = 0; i < xPosition.size(); i++)
+    {
+        
+        h_XYPos->Fill(xPosition[i], yPosition[i]);
+        h_ZRPos->Fill(zPosition[i], rPosition[i]);
+        h_XYEnergy->Fill(xPosition[i], yPosition[i], energy[i]);
+    }
+
+    for (int i = 0; i < zLayerValues.size(); i++)
+    {
+        h_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]);
+        p_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]);
+
+        h_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]);
+        p_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]);
+    }
+
+    for(int i = 0; i < nEvents; i++)
+    {
+        h_energyTotal->Fill(totalEventEnergy[i]);
+        h_hitCount->Fill(totalEventHits[i]);
+    }
+
+    gStyle->SetPadLeftMargin(0.13);
+
+    TCanvas *canvas = new TCanvas("canvas", "canvas", 1600, 800);
+    canvas->Divide(4,2);
+    canvas->cd(1);
+    h_energyTotal->Draw();
+    canvas->cd(2);
+    h_layerEnergy->Draw("COLZ");
+    p_layerEnergy->Draw("SAME");
+    canvas->cd(3);
+    h_hitCount->Draw();
+    canvas->cd(4);
+    h_layerHits->Draw("COLZ");
+    p_layerHits->Draw("SAME");
+    canvas->cd(5);
+    h_XYPos->Draw("COLZ");
+    canvas->cd(6);
+    h_ZRPos->Draw("COLZ");
+    canvas->cd(7);
+    h_XYEnergy->Draw("COLZ");
+
+    canvas->SaveAs(outname_pdf.c_str());
+    canvas->SaveAs(outname_png.c_str());
+
+    return 0;
+}