Skip to content
Snippets Groups Projects
Select Git revision
  • b08e9e16d1052e92c61dee6721c0ae42ef87d052
  • master default protected
  • update_imcal_ml
  • dRICH_residual_fixed_index
  • CherenkovPIDAnalysis_dRICH_103
  • sebouh137-master-patch-84798
  • zdc_sipmontile_ai
  • irt-algo
  • fix_direct_trigger
  • fix-include
  • irt-algo-sensor-normal
  • irt-algo-mod
  • wdconinc-master-patch-03839
  • update_imaging_ml_benchmarks
  • vgawas-new
  • eicrecon
  • vgawas-phy
  • ai_codesign
  • tracking-with-background-overlay
  • 86-ecal-benchmark-fails-but-job-does-not-fail
  • robin-ShaperBranch
21 results

emcal_electrons_analysis.cxx

Blame
  • jihee.kim's avatar
    Jihee Kim authored
    b08e9e16
    History
    emcal_electrons_analysis.cxx 6.59 KiB
    ////////////////////////////////////////
    // Read reconstruction ROOT output file
    // Plot variables
    ////////////////////////////////////////
    
    #include "ROOT/RDataFrame.hxx"
    #include <iostream>
    
    #include "dd4pod/Geant4ParticleCollection.h"
    #include "eicd/ClusterCollection.h"
    #include "eicd/ClusterData.h"
    
    #include "TCanvas.h"
    #include "TStyle.h"
    #include "TMath.h"
    
    // If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
    //  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
    //  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
    
    using ROOT::RDataFrame;
    using namespace ROOT::VecOps;
    
    void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.root")
    {
      // 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);
    
      ROOT::EnableImplicitMT();
      ROOT::RDataFrame d0("events", input_fname);
    
      // Number of Clusters
      auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
    
      // Cluster Energy [GeV]
      auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
      	std::vector<float> result;
      	for (const auto& i: evt)
        		result.push_back(i.energy);
      return result;
      };
    
      // Scattering Angle [degree]
      auto theta = [] (const std::vector<eic::ClusterData>& evt) {
    	std::vector<float> result;
    	for(const auto& i: evt) {
    		auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
    		result.push_back(TMath::ACos(i.position.z/r)*TMath::RadToDeg());
    	}
      return result;
      };
    
      // Pseudo-rapidity
      auto eta = [] (const std::vector<eic::ClusterData>& evt) {
    	std::vector<float> result;
    	for(const auto& i: evt) {
    		auto r  = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
    		auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg(); 
    		result.push_back(-1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0)));
    	}
      return result;
      };
    
      // Mass [GeV]
      auto mass2 = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
      	std::vector<float> result;
        	result.push_back(input[2].mass*input[2].mass);
      return result;
      };
    
      // Momentum [GeV]
      auto p_rec = [](const std::vector<float>& energy_term, const std::vector<float>& mass_term){return TMath::Sqrt((energy_term*energy_term)-(mass_term)); };
      //auto p_rec = [] (const std::vector<eic::ClusterData>& evt) {
      //      std::vector<float> result;
      //      for (const auto& i: evt)
      //              result.push_back(TMath::Sqrt(i.energy*i.energy - mass2));
      //return result;
      //};
    
      // Thrown Momentum [GeV]
      auto p_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
            std::vector<float> result;
            result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz));
      return result;
      };
    
      // Thrown Energy [GeV]
      auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
      	std::vector<float> result;
        	result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
      return result;
      };
    
      // Add Acceptance and Erec-Ethr for resolution
      auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
    	      .Define("clusterE", clusterE, {"EcalClusters"})
    	      .Define("theta",    theta,    {"EcalClusters"})
    	      .Define("eta",      eta,      {"EcalClusters"})
    	      .Define("mass2",    mass2,    {"mcparticles2"})
    	      .Define("p_rec",  p_rec,    {"clusterE","mass2"})
    	      .Define("p_thr",    p_thr,    {"mcparticles2"})
    	      .Define("E_thr",    E_thr,    {"mcparticles2"})
    	      ;
    
      // Define Histograms
      auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events",    5, -0.5, 4.5},      "ncluster");
      auto hClusterE = d1.Histo1D({"hClusterE", "Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 30.5},   "clusterE");
      auto hTheta    = d1.Histo1D({"hTheta",    "Scattering Angle; #theta [degree]; Events",    100, 130.0, 180.0}, "theta");
      auto hEta      = d1.Histo1D({"hEta",      "Pseudo-rapidity; #eta; Events",                100, -5.0, 0.0},    "eta");
      auto hPrec     = d1.Histo1D({"hPrec",     "Reconstructed Momentum; p_{rec}[GeV]; Events", 100, -0.5, 30.5},   "p_rec");
      auto hPthr     = d1.Histo1D({"hPthr",     "Thrown Momentum; p_{thr}[GeV]; Events",        100, -0.5, 30.5},   "p_thr");
      auto hEthr     = d1.Histo1D({"hEthr",     "Thrown Energy; E_{thr}[GeV]; Events",          100, -0.5, 30.5},   "E_thr");
    
      // Draw Histograms
      TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
      c1->SetLogy(1); 
      hNCluster->GetYaxis()->SetTitleOffset(1.6);
      hNCluster->SetLineWidth(2);
      hNCluster->SetLineColor(kBlue);
      hNCluster->DrawClone();
      c1->SaveAs("results/emcal_electrons_ncluster.png");
      c1->SaveAs("results/emcal_electrons_ncluster.pdf");
    
      TCanvas *c2 = new TCanvas("c2", "c2", 500, 500); 
      hClusterE->GetYaxis()->SetTitleOffset(1.4);
      hClusterE->SetLineWidth(2);
      hClusterE->SetLineColor(kBlue);
      hClusterE->DrawClone();
      c2->SaveAs("results/emcal_electrons_clusterE.png");
      c2->SaveAs("results/emcal_electrons_clusterE.pdf");
    
      TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
      hTheta->GetYaxis()->SetTitleOffset(1.4);
      hTheta->SetLineWidth(2);
      hTheta->SetLineColor(kBlue);
      hTheta->DrawClone();
      c3->SaveAs("results/emal_electrons_theta.png"); 
      c3->SaveAs("results/emal_electrons_theta.pdf");
    
      TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
      hEta->GetYaxis()->SetTitleOffset(1.4);
      hEta->SetLineWidth(2);
      hEta->SetLineColor(kBlue);
      hEta->DrawClone();
      c4->SaveAs("results/emcal_electrons_eta.png");
      c4->SaveAs("results/emcal_electrons_eta.pdf");
    
      TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
      hPrec->GetYaxis()->SetTitleOffset(1.4);
      hPrec->SetLineWidth(2);
      hPrec->SetLineColor(kBlue);
      hPrec->DrawClone();
      c5->SaveAs("results/emcal_electrons_Prec.png");
      c5->SaveAs("results/emcal_electrons_Prec.pdf");
    
      TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
      hPthr->GetYaxis()->SetTitleOffset(1.4);
      hPthr->SetLineWidth(2);
      hPthr->SetLineColor(kBlue);
      hPthr->DrawClone();
      c6->SaveAs("results/emcal_electrons_Pthr.png");
      c6->SaveAs("results/emcal_electrons_Pthr.pdf");
    
      TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
      hEthr->GetYaxis()->SetTitleOffset(1.4);
      hEthr->SetLineWidth(2);
      hEthr->SetLineColor(kBlue);
      hEthr->DrawClone();
      c7->SaveAs("results/emcal_electrons_Ethr.png");
      c7->SaveAs("results/emcal_electrons_Ethr.pdf");
    
    
      //h->Fit("gaus");
    }