Skip to content
Snippets Groups Projects
Commit 3c24d6b6 authored by Maria Zurek's avatar Maria Zurek
Browse files

Restructure dataframe

parent fa067d4a
No related branches found
No related tags found
1 merge request!153Draft: Resolve "Add energy scan for Barrel Ecal"
Hide whitespace changes
Inline Side-by-side
benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis_pi_e.cxx
+ 159
− 110
View file @ 3c24d6b6
...@@ -33,6 +33,25 @@ ...@@ -33,6 +33,25 @@
using ROOT::RDataFrame; using ROOT::RDataFrame;
using namespace ROOT::VecOps; using namespace ROOT::VecOps;
struct infoHist1D {
string name; // histo name
string column; // column name
string titles; // histo and axis titles
int nbins; // Nb of bins
double x0; // Lower hist limit
double x1; // Upper hist limit
};
struct infoHist2D {
string name; // Plot variable name
int nbinsx; // Nb of bins in x
double x0; // Lower hist limit in x
double x1; // Upper hist limit in x
int nbinsy; // Nb of bins in y
double y0; // Lower hist limit in y
double y1; // Upper hist limit in y
};
void save_canvas(TCanvas* c, std::string label) void save_canvas(TCanvas* c, std::string label)
{ {
c->SaveAs(fmt::format("results/{}.png",label).c_str()); c->SaveAs(fmt::format("results/{}.png",label).c_str());
...@@ -64,12 +83,19 @@ void set_histo_range(TH1D* h) ...@@ -64,12 +83,19 @@ void set_histo_range(TH1D* h)
std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_label) std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_label)
{ {
std::string input_fname_electron = fmt::format("sim_output/rec_emcal_barrel_{}_{}.root", "electron", E_label); std::string input_fname_electron = fmt::format("sim_output/rec_emcal_barrel_{}_{}.root", "electron", E_label);
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname_electron);
std::string input_fname_pionm = fmt::format("sim_output/rec_emcal_barrel_{}_{}.root", "pion-", E_label); std::string input_fname_pionm = fmt::format("sim_output/rec_emcal_barrel_{}_{}.root", "pion-", E_label);
ROOT::EnableImplicitMT(); ROOT::EnableImplicitMT();
ROOT::RDataFrame d2("events", input_fname_pionm); ROOT::RDataFrame d0("events", {input_fname_electron, input_fname_pionm});
auto isElectron = [](std::vector<dd4pod::Geant4ParticleData> const& input){
return (input[2].pdgID == 11 ? true : false);
};
auto isPion = [](std::vector<dd4pod::Geant4ParticleData> const& input){
return (input[2].pdgID == -211 ? true : false);
};
// Thrown Energy [GeV] // Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
...@@ -118,126 +144,149 @@ std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_l ...@@ -118,126 +144,149 @@ std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_l
}; };
// Define variables // Define variables
auto delectron = d0.Define("Ethr", Ethr, {"mcparticles2"}) auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"})
.Define("mom", momentum, {"mcparticles2"}) .Define("Mom", momentum, {"mcparticles2"})
.Define("mom_smeared", momentum_smeared, {"mom"}) .Define("MomSmeared", momentum_smeared, {"mom"})
.Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"}) .Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"})
.Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"}) .Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"})
.Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"}) .Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"})
.Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"}) .Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"})
.Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"}) .Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"})
.Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"}) .Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"})
.Define("E_over_p_ClusterScFi", E_over_mom, {"EClusterScFi", "mom_smeared"}) .Define("EoverPClusterScFi", E_over_mom, {"EClusterScFi", "MomSmeared"})
.Define("E_over_p_RecoScFi", E_over_mom, {"ErecScFi", "mom_smeared"}); .Define("EoverPRecoScFi", E_over_mom, {"ErecScFi", "MomSmeared"});
auto dpion = d2.Define("Ethr", Ethr, {"mcparticles2"}) auto dpions = d1.Filter(isPion, {"mcparticles2"});
.Define("mom", momentum, {"mcparticles2"}) auto delectrons = d1.Filter(isElectron, {"mcparticles2"});
.Define("mom_smeared", momentum_smeared, {"mom"})
.Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"}) enum typeColumns {Ethr, Mom, MomSmeared, EClusterScFi, NClusterScFi, EoverPScFi, EoverPRecScFi, nColumnsToPlot};
.Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"})
.Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"}) infoHist1D infoHists[nColumnsToPlot];
.Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"}) infoHists[Ethr] = {"hEthr", "Ethr", "Thrown Energy; Thrown Energy [GeV]; Events", 500, 0.0, 11.0};
.Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"}) infoHists[Mom] = {"hMom", "Mom", "Thrown Momentum; Thrown Momentum [GeV/c]; Events", 500, 0.0, 11.0};
.Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"}) infoHists[MomSmeared] = {"hMomSmeared", "MomSmeared", "Thrown Momentum Smeared; Thrown Momentum Smeared [GeV/c]; Events", 500, 0.0, 11.0};
.Define("E_over_p_ClusterScFi", E_over_mom, {"EClusterScFi", "mom_smeared"}) infoHists[EClusterScFi] = {"hEClusterScFi", "EClusterScFi", "Cluster Energy; Cluster Energy [GeV]; Events", 500, 0.0, 11.0};
.Define("E_over_p_RecoScFi", E_over_mom, {"ErecScFi", "mom_smeared"}); infoHists[NClusterScFi] = {"hNClusterScFi", "NClusterScFi", "Number of Clusters; # of Clusters; Events", 10, 0, 10};
infoHists[EoverPScFi] = {"hEoverPScFi", "EoverPClusterScFi", "Cluster Energy/Momentum; Cluster Energy/Momentum; Events", 480, 0., 1.2};
infoHists[EoverPRecScFi] = {"hEoverPScFi", "EoverPRecoScFi", "Reco Hits Energy/Momentum; Reco Hits Energy/Momentum; Events", 500, 0.0, 0.25};
// Define and draw Histograms
for(int col=0; col<nColumnsToPlot; col++){
auto infoHist = infoHists[col];
auto helectron = delectrons.Histo1D({fmt::format("{}_electron",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
auto hpion = dpions.Histo1D({fmt::format("{}_pion",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
auto c = new TCanvas(fmt::format("c{}",infoHist.comlumn).c_str(), fmt::format("c{}",infoHist.comlumn).c_str(), 700, 500);
helectron->GetYaxis()->SetTitleOffset(1.4);
helectron->SetLineWidth(2);
helectron->SetLineColor(kBlue);
helectron->DrawClone();
hpion->SetLineColor(kRed);
hpion->DrawClone("same");
save_canvas(c, infoHist.column.c_str(), E_label, "electron_pion");
}
// Define Histograms auto nlayers = 20;
enum typeLayerColumns {ErecScFi, ErecImg, nhitsScFi, nhitsImg, nLayerColumnsToPlot};
auto hEthr_el = delectron.Histo1D({"hEthr_el", "Thrown Energy; Thrown Energy [GeV]; Events", 1000, 0.0, 11.0},"Ethr"); std::string cname[nLayerColumnsToPlot] = {"ErecScFi_layer", "ErecImg_layer", "nhitsScFi_layer", "nhitsImg_layer"};
auto hMom_el = delectron.Histo1D({"hMom_el", "Thrown Momentum; Thrown Momentum [GeV/c]; Events", 1000, 0.0, 11.0},"mom"); TCanvas *clayer[nLayerColumnsToPlot];
auto hMomSmeared_el = delectron.Histo1D({"hMom_el", "Thrown Momentum Smeared; Thrown Momentum Smeared [GeV/c]; Events", 1000, 0.0, 11.0},"mom_smeared"); for(int col=0; col<nLayerColumnsToPlot; col++){
clayer[col] = new TCanvas(fmt::format("clayer{}",col).c_str(),fmt::format("clayer{}",col).c_str(),1250, 1000);
// auto hEClusterImg_el = delectron.Histo1D({"hEClusterImg", "Cluster Energy; Cluster Energy [GeV]; Events", 200, 0.0, 25.0},"EClusterImg"); clayer[col]->Divide(4,5);
// auto hNClusterImg_el = delectron.Histo1D({"hNClusterImg", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "NClusterImg"); }
// auto hEoverPImg_el = delectron.Histo1D({"hfsamImg", "Cluster Energy/E true; Cluster Energy/E true; Events", 100, 0.0, 1.5},"fsamClusterImg");
// auto hEoverPRecImg_el = delectron.Histo1D({"hfsamRecImg", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.1},"fsamRecImg");
auto hEClusterScFi_el = delectron.Histo1D({"hEClusterScFi_el", "Cluster Energy; Cluster Energy [GeV]; Events", 500, 0.0, 11.0},"EClusterScFi");
auto hNClusterScFi_el = delectron.Histo1D({"hNClusterScFi_el", "Number of Clusters; # of Clusters; Events", 10, 0, 10}, "NClusterScFi");
auto hEoverPScFi_el = delectron.Histo1D({"hEoverPScFi_el", "Cluster Energy/E true; Cluster Energy/E true; Events", 480, 0., 1.2},"E_over_p_ClusterScFi");
auto hEoverPRecScFi_el = delectron.Histo1D({"hEoverPRecScFi_el", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 100, 0.0, 0.25},"E_over_p_RecoScFi");
auto hEthr_pion = dpion.Histo1D({"hEthr_pion", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},"Ethr");
auto hMom_pion = dpion.Histo1D({"hMom_pion", "Thrown Momentum; Thrown Momentum [GeV/c]; Events", 100, 0.0, 25.0},"mom");
auto hMomSmeared_pion = dpion.Histo1D({"hMomSmeared_pion", "Thrown Momentum Smeared; Thrown Momentum Smeared [GeV/c]; Events", 100, 0.0, 25.0},"mom_smeared");
// auto hEClusterImg_pion = dpion.Histo1D({"hEClusterImg", "Cluster Energy; Cluster Energy [GeV]; Events", 200, 0.0, 25.0},"EClusterImg");
// auto hNClusterImg_pion = dpion.Histo1D({"hNClusterImg", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "NClusterImg");
// auto hEoverPImg_pion = dpion.Histo1D({"hfsamImg", "Cluster Energy/E true; Cluster Energy/E true; Events", 100, 0.0, 1.5},"fsamClusterImg");
// auto hEoverPRecImg_pion = dpion.Histo1D({"hfsamRecImg", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.1},"fsamRecImg");
auto hEClusterScFi_pion = dpion.Histo1D({"hEClusterScFi_pion", "Cluster Energy; Cluster Energy [GeV]; Events", 500, 0.0, 11.0},"EClusterScFi"); for(int layer = 1; layer < nlayers + 1; layer++) {
auto hNClusterScFi_pion = dpion.Histo1D({"hNClusterScFi_pion", "Number of Clusters; # of Clusters; Events", 10, 0, 10}, "NClusterScFi");
auto hEoverPScFi_pion = dpion.Histo1D({"hEoverPScFi_pion", "Cluster Energy/E true; Cluster Energy/E true; Events", 480, 0., 1.2},"E_over_p_ClusterScFi");
auto hEoverPRecScFi_pion = dpion.Histo1D({"hEoverPRecScFi_pion", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.25},"E_over_p_RecoScFi");
std::cout << "Layer to process: " << layer << std::endl;
// Energy reconstructed per layer
auto ERecInLayer = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
auto layer_edep = 0.0;
for (const auto& i: evt) {
if (i.layer == layer) {
layer_edep += i.energyDeposit;
}
}
return layer_edep;
};
// Hits reconstructed per layer
auto nhitsInLayer = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
auto nhits = 0;
for (const auto& i: evt) {
if (i.layer == layer) {
nhits++;
}
}
return nhits;
};
// dE up to the layer
auto ERecUpToLayer = [=](const std::vector<eic::CalorimeterHitData>& evt) {
double edep = 0.0;
for (const auto& i: evt){
if (i.layer < layer+1){edep += i.energy;}
}
return edep;
};
auto d2 = d0.Define(fmt::format("ErecImg_layer_{}",layer).c_str(), ERecInLayer, {"RecoEcalBarrelImagingHits"})
.Define(fmt::format("ErecScFi_layer_{}",layer).c_str(), ERecInLayer, {"EcalBarrelScFiHitsReco"})
.Define(fmt::format("Erec_layer_{}",layer).c_str(), fmt::format("EsimImg_layer_{}+EsimScFi_layer_{}",layer,layer).c_str())
.Define(fmt::format("nhitsImg_layer_{}",layer).c_str(), nhitsInLayer, {"RecoEcalBarrelImagingHits"})
.Define(fmt::format("nhitsScFi_layer_{}",layer).c_str(), nhitsInLayer, {"EcalBarrelScFiHitsReco"})
.Define(fmt::format("ErecSumImg_layer_{}",layer).c_str(), ERecUpToLayer, {"RecoEcalBarrelImagingHits"})
.Define(fmt::format("ErecSumScFi_layer_{}",layer).c_str(), ERecUpToLayer, {"EcalBarrelScFiHitsReco"});
auto dpionslayer = d2.Filter(isPion,{"mcparticles2"});
auto delectronslayer = d2.Filter(isElectron,{"mcparticles2"});
infoHist1D infoLayerHists[nLayerColumnsToPlot];
infoLayerHists[ErecImg] = {fmt::format("hErecImg_layer_{}",layer), fmt::format("ErecImg_layer_{}",layer), fmt::format("Energy Deposit in layer {}; Energy Deposit [Gev]; Events",layer), 200, 0.0, 0.04};
infoLayerHists[ErecScFi] = {fmt::format("hErecScFi_layer_{}",layer), fmt::format("ErecScFi_layer_{}",layer), fmt::format("Energy Deposit in layer {}; Energy Deposit [Gev]; Events",layer), 200, 0.0, 0.1};
infoLayerHists[nhitsImg] = {fmt::format("hnhitsImg_layer_{}_el",layer), fmt::format("nhitsImg_layer_{}",layer), fmt::format("Number of Reconstructed Hits in layer (); Number of hits; Events",layer), 500, 0, 500};
infoLayerHists[nhitsScFi] = {fmt::format("hnhitsScFi_layer_{}_el",layer), fmt::format("nhitsScFi_layer_{}",layer), fmt::format("Number of Reconstructed Hits in layer (); Number of hits; Events",layer), 500, 0, 500};
for(int col=0; col<nLayerColumnsToPlot; col++){
auto infoHist = infoLayerHists[col];
auto helectron = delectronslayer.Histo1D({fmt::format("{}_electron",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
auto hpion = dpionslayer.Histo1D({fmt::format("{}_pion",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
clayer[col]->cd(layer);
gPad->SetLogy();
helectron->GetYaxis()->SetTitleOffset(1.4);
helectron->SetLineWidth(2);
helectron->SetLineColor(kBlue);
helectron->DrawClone();
hpion->SetLineColor(kRed);
hpion->DrawClone("same");
// auto h = hEsim_layer->DrawCopy();
// auto up_range = h->GetMean() + 3*h->GetStdDev();
// auto down_range = h->GetMean() - 3*h->GetStdDev();
// h->SetLineWidth(2);
// h->SetLineColor(kBlue);
// h->GetXaxis()->SetRange(h->GetXaxis()->GetBinUpEdge(1), up_range); // skip 0th bin
// auto mean_layer = h->GetMean();
// auto rms_layer = h->GetStdDev();
// h->GetXaxis()->SetRange(); // reset the range
// h->GetXaxis()->SetRangeUser(0.,up_range);
// auto no_edep = (h->GetBinContent(1)/h->GetEntries())*100;
// gr_no_edep.SetPoint(gr_no_edep.GetN(),layer,no_edep);
// gr_edep_mean.SetPoint(gr_edep_mean.GetN(),layer, mean_layer);
// gr_edep_mean.SetPointError(gr_edep_mean.GetN()-1,0, rms_layer);
}
for(int col=0; col<nLayerColumnsToPlot; col++) save_canvas(clayer[col], cname[col] , E_label, "electron_pion");
// Event Counts // Event Counts
auto nevents_thrown = delectron.Count(); auto nevents_thrown = delectron.Count();
std::cout << "Number of Thrown Electron Events: " << (*nevents_thrown) << "\n"; std::cout << "Number of Thrown Electron Events: " << (*nevents_thrown) << "\n";
// Draw Histograms
{
TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
c1->Divide(2,2);
c1->SetLogy(1);
c1->cd(1);
auto h = hEthr_el->DrawCopy();
h->GetYaxis()->SetTitleOffset(1.4);
h->SetLineWidth(2);
h->SetLineColor(kBlue);
auto h_pion = hEthr_pion->DrawCopy("same");
h_pion->SetLineColor(kRed);
c1->cd(2);
hMom_el->GetYaxis()->SetTitleOffset(1.4);
hMom_el->SetLineWidth(2);
hMom_el->SetLineColor(kBlue);
hMom_el->DrawClone();
auto up_fit = hMom_el->GetMean() + 5*hMom_el->GetStdDev();
auto down_fit = hMom_el->GetMean() - 5*hMom_el->GetStdDev();
hMom_el->GetXaxis()->SetRangeUser(down_fit,up_fit);
hMom_pion->SetLineColor(kRed);
hMom_pion->DrawClone("same");
c1->cd(3);
hMomSmeared_el->GetYaxis()->SetTitleOffset(1.4);
hMomSmeared_el->SetLineWidth(2);
hMomSmeared_el->SetLineColor(kBlue);
hMomSmeared_el->DrawClone();
hMomSmeared_el->GetXaxis()->SetRangeUser(down_fit,up_fit);
hMomSmeared_pion->SetLineColor(kRed);
hMomSmeared_pion->DrawClone("same");
save_canvas(c1, "gen", E_label, "electron_pion");
}
TCanvas* c2 = new TCanvas("c2", "c2", 1400, 1000);
hNClusterScFi_el->GetYaxis()->SetTitleOffset(1.4);
hNClusterScFi_el->SetLineWidth(2);
hNClusterScFi_el->SetLineColor(kBlue);
hNClusterScFi_el->DrawClone();
hNClusterScFi_pion->SetLineColor(kRed);
hNClusterScFi_pion->DrawClone("same");
save_canvas(c2, "NCluster", E_label, "electron_pion");
TCanvas* c3 = new TCanvas("c3", "c3", 1400, 1000);
hEoverPScFi_el->GetYaxis()->SetTitleOffset(1.4);
hEoverPScFi_el->SetLineWidth(2);
hEoverPScFi_el->SetLineColor(kBlue);
hEoverPScFi_el->DrawClone();
hEoverPScFi_pion->SetLineColor(kRed);
hEoverPScFi_pion->DrawClone("same");
save_canvas(c3, "EoverPClusterScFi", E_label, "electron_pion");
TCanvas* c4 = new TCanvas("c4", "c4", 1400, 1000);
hEoverPRecScFi_el->GetYaxis()->SetTitleOffset(1.4);
hEoverPRecScFi_el->SetLineWidth(2);
hEoverPRecScFi_el->SetLineColor(kBlue);
hEoverPRecScFi_el->DrawClone();
hEoverPRecScFi_pion->SetLineColor(kRed);
hEoverPRecScFi_pion->DrawClone("same");
save_canvas(c4, "EoverPRecScFi", E_label, "electron_pion");
return std::make_tuple(0., 0.); return std::make_tuple(0., 0.);
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment