File indexing completed on 2024-09-27 07:02:37
0001
0002
0003
0004
0005
0006 #include "ROOT/RDataFrame.hxx"
0007 #include <iostream>
0008 #include <algorithm>
0009 #include <string>
0010
0011 #include "edm4hep/MCParticleCollection.h"
0012 #include "edm4hep/SimCalorimeterHitCollection.h"
0013
0014 #include "common_bench/benchmark.h"
0015 #include "common_bench/mt.h"
0016 #include "common_bench/util.h"
0017
0018 #include <boost/range/combine.hpp>
0019
0020 #include "DD4hep/IDDescriptor.h"
0021 #include "DD4hep/Readout.h"
0022 #include "DD4hep/Segmentations.h"
0023
0024 #include "TCanvas.h"
0025 #include "TStyle.h"
0026 #include "TMath.h"
0027 #include "TH1.h"
0028 #include "TF1.h"
0029 #include "TH1D.h"
0030 #include "TFitResult.h"
0031 #include "TLegend.h"
0032 #include "TString.h"
0033 #include "TGraph.h"
0034 #include "TGraph2D.h"
0035 #include "TGraphErrors.h"
0036 #include "TLine.h"
0037 #include "TError.h"
0038
0039 R__LOAD_LIBRARY(libfmt.so)
0040 #include "fmt/core.h"
0041 #include "DD4hep/Detector.h"
0042 #include "DDG4/Geant4Data.h"
0043 #include "DDRec/CellIDPositionConverter.h"
0044 #include "emcal_barrel_common_functions.h"
0045 #include <nlohmann/json.hpp>
0046
0047 using ROOT::RDataFrame;
0048 using namespace ROOT::VecOps;
0049
0050 void emcal_barrel_pion_rejection_analysis(
0051 const char* input_fname1 = "sim_output/sim_emcal_barrel_piRej_electron.edm4hep.root",
0052 const char* input_fname2 = "sim_output/sim_emcal_barrel_piRej_piminus.edm4hep.root"
0053 )
0054 {
0055
0056 gErrorIgnoreLevel = kFatal;
0057
0058
0059 gROOT->SetStyle("Plain");
0060 gStyle->SetOptFit(1);
0061 gStyle->SetLineWidth(2);
0062 gStyle->SetPadTickX(1);
0063 gStyle->SetPadTickY(1);
0064 gStyle->SetPadGridX(1);
0065 gStyle->SetPadGridY(1);
0066 gStyle->SetPadLeftMargin(0.14);
0067 gStyle->SetPadRightMargin(0.14);
0068
0069 ROOT::EnableImplicitMT();
0070 ROOT::RDataFrame d0("events", {input_fname1, input_fname2});
0071
0072
0073 if (! d0.HasColumn("EcalBarrelScFiHits")) {
0074 std::cout << "EcalBarrelScFiHits is required" << std::endl;
0075 return;
0076 }
0077
0078
0079 std::string detector_path = "";
0080 std::string detector_name = "athena";
0081 if(std::getenv("DETECTOR_PATH")) {
0082 detector_path = std::getenv("DETECTOR_PATH");
0083 }
0084 if(std::getenv("DETECTOR_CONFIG")) {
0085 detector_name = std::getenv("DETECTOR_CONFIG");
0086 }
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099 int layerNum;
0100 int dep_min = 1;
0101 int dep_max = 6;
0102
0103
0104 dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0105 detector.fromCompact(fmt::format("{}/{}.xml", detector_path, detector_name));
0106
0107 auto decoder = detector.readout("EcalBarrelImagingHits").idSpec().decoder();
0108 auto decoderScFi = detector.readout("EcalBarrelScFiHits").idSpec().decoder();
0109 auto layer_index = decoder->index("layer");
0110 auto layer_indexScFi = decoderScFi->index("layer");
0111
0112
0113 auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0114 return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z + input[2].mass*input[2].mass);
0115 };
0116
0117
0118 auto Pthr = [](std::vector<edm4hep::MCParticleData> const& input) {
0119 return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z);
0120 };
0121
0122
0123 auto Eta = [](std::vector<edm4hep::MCParticleData> const& input) {
0124 double E = TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z + input[2].mass*input[2].mass);
0125 return 0.5*TMath::Log((E + input[2].momentum.z) / (E - input[2].momentum.z));
0126 };
0127
0128
0129 auto pT = [](std::vector<edm4hep::MCParticleData> const& input) {
0130 return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y);
0131 };
0132
0133
0134 auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
0135
0136
0137 auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0138 double total_edep = 0.0;
0139 for (const auto& i: evt){
0140 total_edep += i.energy;
0141 }
0142 return total_edep;
0143 };
0144
0145
0146 auto Esim_dep2 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0147 auto total_edep = 0.0;
0148 for (const auto& i: evt) {
0149 if( decoder->get(i.cellID, layer_index) < 3 ){
0150 total_edep += i.energy;
0151 }
0152 }
0153 return total_edep;
0154 };
0155
0156
0157
0158 auto Esim_dep3 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0159 auto total_edep = 0.0;
0160 for (const auto& i: evt) {
0161 if( decoder->get(i.cellID, layer_index) < 4 ){
0162 total_edep += i.energy;
0163 }
0164 }
0165 return total_edep;
0166 };
0167
0168
0169 auto Esim_dep4 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0170 auto total_edep = 0.0;
0171 for (const auto& i: evt) {
0172 if( decoder->get(i.cellID, layer_index) < 5 ){
0173 total_edep += i.energy;
0174 }
0175 }
0176 return total_edep;
0177 };
0178
0179
0180 auto Esim_dep5 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0181 auto total_edep = 0.0;
0182 for (const auto& i: evt) {
0183 if( decoder->get(i.cellID, layer_index) < 6 ){
0184 total_edep += i.energy;
0185 }
0186 }
0187 return total_edep;
0188 };
0189
0190
0191 auto Esim_dep6 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0192 auto total_edep = 0.0;
0193 for (const auto& i: evt) {
0194 if( decoder->get(i.cellID, layer_index) < 7 ){
0195 total_edep += i.energy;
0196 }
0197 }
0198 return total_edep;
0199 };
0200
0201
0202 auto Esim_dep6_ScFi = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0203 auto total_edep = 0.0;
0204 for (const auto& i: evt) {
0205 if( decoderScFi->get(i.cellID, layer_indexScFi) < 7 ){
0206 total_edep += i.energy;
0207 }
0208 }
0209 return total_edep;
0210 };
0211
0212
0213 auto Esim_dep7 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0214 auto total_edep = 0.0;
0215 for (const auto& i: evt) {
0216 if( decoder->get(i.cellID, layer_index) < 8 ){
0217 total_edep += i.energy;
0218 }
0219 }
0220 return total_edep;
0221 };
0222
0223 auto Esim_dep8 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0224 auto total_edep = 0.0;
0225 for (const auto& i: evt) {
0226 if( decoder->get(i.cellID, layer_index) < 9 ){
0227 total_edep += i.energy;
0228 }
0229 }
0230 return total_edep;
0231 };
0232
0233
0234 auto Esim_dep9 = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0235 auto total_edep = 0.0;
0236 for (const auto& i: evt) {
0237 if( decoder->get(i.cellID, layer_index) < 10 ){
0238 total_edep += i.energy;
0239 }
0240 }
0241 return total_edep;
0242 };
0243
0244
0245
0246 auto Esim_dep = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0247 std::vector<double>res(20);
0248 for (const auto& i: evt) {
0249 res[decoder->get(i.cellID, layer_index)] += i.energy;
0250 }
0251 return res;
0252 };
0253
0254
0255 auto Esim_dep_sum = [&dep_min, &dep_max](const std::vector<double>& dep) {
0256 double res = 0;
0257 for (int i = dep_min; i < dep_max + 1; i++) {
0258 if (i >= dep.size()) continue;
0259 res += dep[i];
0260 }
0261 return res;
0262 };
0263
0264
0265 auto Esim_depN = [&layerNum](const std::vector<double>& dep) {
0266 return (layerNum < dep.size() ? dep[layerNum] : 0.0);
0267 };
0268
0269
0270 auto fsam = [](const double& sampled, const double& thrown) {
0271 return sampled / thrown;
0272 };
0273
0274
0275 auto fEp = [](const double& E_front, const double& mom) {
0276 return E_front / mom;
0277 };
0278
0279
0280 auto getpid = [](std::vector<edm4hep::MCParticleData> const& input) {
0281 return input[2].PDG;
0282 };
0283
0284
0285 auto getdau = [](std::vector<edm4hep::MCParticleData> const& input){
0286 return input[2].daughters_begin;
0287 };
0288
0289
0290 auto is_electron = [](std::vector<edm4hep::MCParticleData> const& input){
0291 return (input[2].PDG == 11 ? true : false);
0292 };
0293
0294
0295 auto is_piMinus = [](std::vector<edm4hep::MCParticleData> const& input){
0296 return (input[2].PDG == -211 ? true : false);
0297 };
0298
0299
0300 auto EDep_bool = [](std::vector<double> const& input){
0301 return (input[0] > input[1]);
0302 };
0303
0304 auto Diff = [](const double &esim, const double &edep){
0305 return esim-edep;
0306 };
0307
0308
0309
0310 auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"})
0311 .Define("Pthr", Pthr, {"MCParticles"})
0312 .Define("nhits", nhits, {"EcalBarrelImagingHits"})
0313 .Define("Esim", Esim, {"EcalBarrelImagingHits"})
0314 .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
0315 .Define("EsimOverP", fEp, {"Esim", "Pthr"})
0316 .Define("EsimScFiOverP", fEp, {"EsimScFi", "Pthr"})
0317 .Define("EsimTot", "EsimScFi+Esim")
0318 .Define("EsimTotOverP", fEp, {"EsimTot", "Pthr"})
0319 .Define("fsam", fsam, {"Esim","Ethr"})
0320 .Define("pid", getpid, {"MCParticles"})
0321 .Define("EDep", Esim_dep, {"EcalBarrelImagingHits"})
0322 .Define("EDepSum", Esim_dep_sum, {"EDep"})
0323 .Define("EDepN", Esim_depN, {"EDep"})
0324 .Define("EDep2", Esim_dep2, {"EcalBarrelImagingHits"})
0325 .Define("EDep3", Esim_dep3, {"EcalBarrelImagingHits"})
0326 .Define("EDep4", Esim_dep4, {"EcalBarrelImagingHits"})
0327 .Define("EDep5", Esim_dep5, {"EcalBarrelImagingHits"})
0328 .Define("EDep6", Esim_dep6, {"EcalBarrelImagingHits"})
0329 .Define("EDep6OverP", fEp, {"EDep6", "Pthr"})
0330 .Define("EOverP", fEp, {"EDep3", "Pthr"})
0331 .Define("Eta", Eta, {"MCParticles"})
0332 .Define("pT", pT, {"MCParticles"})
0333 .Define("EDepOverP", fEp, {"EDepN", "Pthr"})
0334 .Define("EDepOverPT", fEp, {"EDepN", "pT"})
0335 .Define("EDepSumOverP", fEp, {"EDepSum", "Pthr"})
0336 .Define("EDepSumOverPT", fEp, {"EDepSum", "pT"})
0337 .Define("EDepFrac", fEp, {"EDepSum", "Esim"})
0338 ;
0339
0340
0341 dep_min = 1;
0342 dep_max = 6;
0343 auto d_ele = d1.Filter(is_electron, {"MCParticles"});
0344 auto d_pim = d1.Filter(is_piMinus, {"MCParticles"});
0345
0346
0347 std::string currentCut = "(EDep6OverP>2.5e-3)&&(EDep6>5e-3)";
0348
0349
0350
0351 std::vector<std::string> var = {"Esim [GeV];", "EsimTot [GeV];", "EDep6 [GeV];", "EDep6/p;", "pT [GeV];", "#eta;", "EsimScFi [GeV]", "EsimScFi/p"};
0352 std::vector<std::string> var_save = {"Esim", "EsimTot", "EDep6", "EDep6OverP", "pT", "eta", "EsimScFi", "EsimScFiOverP"};
0353 std::vector<std::string> col = {"Esim", "EsimTot", "EDep6", "EDep6OverP", "pT", "Eta", "EsimScFi", "EsimScFiOverP"};
0354 std::vector<std::vector<double>> h1Ranges = {{0,0.2}, {0, 0.2}, {0,0.25}, {0, 0.02}, {0, 18}, {-1, 1}, {0,0.2}, {0,0.2}};
0355 for (int i = 0; i < var.size(); i++){
0356 std::string title = "#pi^{-}, e^{-};" + var[i] + " Events";
0357 auto he = d_ele.Histo1D({"he", title.c_str(), 100, h1Ranges[i][0], h1Ranges[i][1]}, col[i]);
0358 auto hp = d_pim.Histo1D({"hp", title.c_str(), 100, h1Ranges[i][0], h1Ranges[i][1]}, col[i]);
0359
0360 hp->GetYaxis()->SetTitleOffset(1.4);
0361 he->SetLineWidth(2);
0362 he->SetLineColor(kRed);
0363 hp->SetLineWidth(2);
0364 hp->SetLineColor(kBlue);
0365 auto c = new TCanvas("c", "c", 700, 500);
0366 auto leng = new TLegend(0.7, 0.7, 0.9, 0.9);
0367 if (var[i] != "EsimScFi/p"){
0368 hp->DrawClone();
0369 he->DrawClone("same");
0370 }
0371 else {
0372 he->DrawClone();
0373 hp->DrawClone("same");
0374 }
0375 c->Update();
0376
0377 leng->AddEntry(he.GetPtr(),"e^{-}","l");
0378 leng->AddEntry(hp.GetPtr(),"#pi^{-}","l");
0379 leng->Draw();
0380 c->SaveAs(("results/emcal_barrel_pion_rej_uncut_comb_" + var_save[i] + ".png").c_str());
0381 }
0382
0383
0384
0385
0386
0387 std::string cutEEta;
0388 std::vector<std::vector<double>> EBins = {{0,2}, {2, 4}, {4, 6}, {6, 9}, {9, 12}, {12, 19}};
0389 for (int i = 0; i < EBins.size(); i++){
0390 std::string minCut = "Pthr>="+std::to_string(EBins[i][0]);
0391 std::string maxCut = "Pthr<"+std::to_string(EBins[i][1]);
0392 cutEEta += "(" + minCut + "&&" + maxCut + "&&";
0393
0394 for (int j = -1; j < 1; j++){
0395 std::string title = "#pi^{-}, e^{-}";
0396 title += fmt::format(" : {} < E < {}", EBins[i][0], EBins[i][1]);
0397 title += fmt::format(" & {} < #eta < {};EDep6/p; Events", j, j+1);
0398 std::string etaCutMin = fmt::format("Eta>={}", j);
0399 std::string etaCutMax = fmt::format("Eta<{}",j+1);
0400 cutEEta += "(" + etaCutMin + "&&" + etaCutMax + "&&";
0401 auto he = d_ele.Filter(minCut).Filter(maxCut).Filter(etaCutMin).Filter(etaCutMax).Histo1D({"he", title.c_str(), 50, 0, 0.02}, "EDep6OverP");
0402 auto hp = d_pim.Filter(minCut).Filter(maxCut).Filter(etaCutMin).Filter(etaCutMax).Histo1D({"hp", title.c_str(), 50, 0, 0.02}, "EDep6OverP");
0403 auto hecopy = he->DrawCopy();
0404 hecopy->Fit("gaus", "", "", 0, hecopy->GetMaximum());
0405 double* res = hecopy->GetFunction("gaus")->GetParameters();
0406 cutEEta += fmt::format("EDep6OverP>={}", res[1] - 2.0*res[2]);
0407 cutEEta += fmt::format("&&EDep6OverP<{})||",res[1] + 3.0*res[2]);
0408
0409 hp->GetYaxis()->SetTitleOffset(1.4);
0410 he->SetLineWidth(2);
0411 he->SetLineColor(kRed);
0412 hp->SetLineWidth(2);
0413 hp->SetLineColor(kBlue);
0414 auto c = new TCanvas("c", "c", 700, 500);
0415 auto leng = new TLegend(0.7, 0.7, 0.9, 0.9);
0416 hp->DrawClone();
0417 he->DrawClone("same");
0418 c->Update();
0419
0420 leng->AddEntry(he.GetPtr(),"e^{-}","l");
0421 leng->AddEntry(hp.GetPtr(),"#pi^{-}","l");
0422 leng->Draw();
0423 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_uncut_comb_E{}Eta{}.png", i, j+1)).c_str());
0424 }
0425 cutEEta.pop_back();
0426 cutEEta.pop_back();
0427 cutEEta += ")||";
0428 }
0429 cutEEta.pop_back();
0430 cutEEta.pop_back();
0431 currentCut = cutEEta;
0432
0433
0434 auto d_pim_cut = d_pim.Filter(currentCut.c_str());
0435 auto d_ele_cut = d_ele.Filter(currentCut.c_str());
0436
0437
0438 std::vector<double> E = {5, 10, 18};
0439 std::vector<double> ledges5 = {2.8, 0.4, 0.3, 0.5};
0440 std::vector<double> ledges10 = {1.4, 0.5, 0.6, 1.0};
0441 std::vector<double> ledges18 = {0.9, 0.9, 1.0, 1.8};
0442 std::vector<double> maxRate5 = {0.1, 100, 500, 1000};
0443 std::vector<double> maxRate10 = {10, 400, 800, 1000};
0444 std::vector<double> maxRate18 = {200, 800, 1000, 100};
0445 std::vector<vector<double>> maxRate = {maxRate5, maxRate10, maxRate18};
0446 std::vector<vector<double>> lowEdges = {ledges5, ledges10, ledges18};
0447 double suppression = 1e-4;
0448 std::vector<vector<double>> rejRatios = lowEdges;
0449 std::vector<vector<double>> effEle = lowEdges;
0450 std::vector<vector<double>> effPim = lowEdges;
0451 std::vector<std::string> etaBin = {"Eta >= -3.5 && Eta < -2.0", "Eta >= -2.0 && Eta < -1.0", "Eta >= -1.0 && Eta < 0", "Eta >= 0 && Eta < 1.0"};
0452 std::vector<std::string> etaTitle = {"-3.5 < #eta < -2.0", "-2.0 < #eta < -1.0", "-1.0 < #eta < 0", "0 < #eta < 1.0"};
0453
0454
0455 std::vector<double> pBins = {0.1,0.2,0.3,0.4,0.5,1,2,3,4,5,10,12,14,16,18};
0456 EBins = {{0,2}, {2, 4}, {4, 6}, {6, 9}, {9, 12}, {12, 19}};
0457 auto tg = new TGraphErrors();
0458
0459 for (int i = 0; i < EBins.size(); i++){
0460 std::string filter = currentCut;
0461 filter += "&&(Pthr>=" + std::to_string(EBins[i][0]) + "&&Pthr<" + std::to_string(EBins[i][1]) + ")";
0462 double numer = (double)*d_ele_cut.Filter(filter.c_str()).Count();
0463 double denom = (double)*d_pim_cut.Filter(filter.c_str()).Count();
0464 double error = std::sqrt(std::pow(numer / denom, 2.0)*(1.0/numer + 1.0/denom));
0465 double ratio = numer / denom;
0466 if (denom == 0){ratio = 1; error = 1;}
0467 tg->SetPoint(i, 0.5*(EBins[i][0] + EBins[i][1]), ratio);
0468 tg->SetPointError(i, 0.5*(EBins[i][1] - EBins[i][0]), error);
0469 }
0470 double e_eff = (double)*d_ele_cut.Count() / (double)*d_ele.Count();
0471 tg->SetTitle(("#pi Rejection with #varepsilon_{e} = "+ std::to_string(e_eff)).c_str());
0472 tg->GetXaxis()->SetTitle("p [GeV]");
0473 tg->GetYaxis()->SetTitle("R_{e/#pi}");
0474 tg->SetMarkerColor(kBlue);
0475 tg->SetMarkerStyle(20);
0476
0477 auto cp = new TCanvas("cp", "cp");
0478 cp->SetLogy();
0479 cp->SetLogx();
0480 tg->DrawClone("ap");
0481 cp->SaveAs("results/emcal_barrel_pion_rej_RatioRej.png");
0482
0483
0484
0485
0486 dep_min = 1;
0487 dep_max = 6;
0488 for (int i = 0; i < 3; i++){
0489 for (int j = 2; j < 4; j++){
0490
0491
0492 std::string pCut = "Pthr>=" + std::to_string(lowEdges[i][j]) + "&&Pthr<" + std::to_string(E[i]);
0493 auto e_eta = d_ele.Filter(etaBin[j]).Filter(pCut);
0494 auto p_eta = d_pim.Filter(etaBin[j]).Filter(pCut);
0495
0496
0497 std::string title = "e^{-} (E = " + std::to_string((int)E[i]) + " GeV) : " + etaTitle[j] + "; p [GeV]; Events";
0498 auto he = e_eta.Histo1D({"he", title.c_str(), 100, lowEdges[i][j], E[i]}, "Pthr");
0499 he->SetLineColor(kBlue);
0500 auto he_cut = e_eta.Filter(currentCut).Histo1D({"he_cut", title.c_str(), 100, lowEdges[i][j], E[i]}, "Pthr");
0501 he_cut->GetYaxis()->SetTitleOffset(1.4);
0502 he_cut->SetLineWidth(2);
0503 he_cut->SetLineColor(kRed);
0504
0505 title = "#pi^{-} (E = " + std::to_string((int)E[i]) + " GeV) : " + etaTitle[j] + "; p [GeV]; Events";
0506 auto hp = p_eta.Histo1D({"hp", title.c_str(), 100, lowEdges[i][j], E[i]}, "Pthr");
0507 hp->SetLineColor(kBlue);
0508 auto hp_cut = p_eta.Filter(currentCut).Histo1D({"hp", title.c_str(), 100, lowEdges[i][j], E[i]}, "Pthr");
0509 hp_cut->GetYaxis()->SetTitleOffset(1.4);
0510 hp_cut->SetLineWidth(2);
0511 hp_cut->SetLineColor(kRed);
0512
0513 auto c = new TCanvas("c", "c", 700, 500);
0514 c->SetLogy(1);
0515 he->DrawClone();
0516 he_cut->DrawClone("same");
0517 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_mom_ele_E{}_eta{}.png", (int)E[i], j)).c_str());
0518 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_mom_ele_E{}_eta{}.pdf", (int)E[i], j)).c_str());
0519 c->Clear();
0520
0521 hp->DrawClone();
0522 hp_cut->DrawClone("same");
0523 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_mom_pim_E{}_eta{}.png", (int)E[i], j)).c_str());
0524 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_mom_pim_E{}_eta{}.pdf", (int)E[i], j)).c_str());
0525 c->Clear();
0526
0527
0528
0529 rejRatios[i][j] = (double)hp_cut->Integral() / (double)he_cut->Integral();
0530 effPim[i][j] = (double)hp_cut->Integral() / (double)hp->Integral();
0531 effEle[i][j] = (double)he_cut->Integral() / (double)he->Integral();
0532
0533 hp_cut->Divide(he.GetPtr());
0534 title = "#pi^{-}/e^{-} (E = " + std::to_string((int)E[i]) + " GeV) : " + etaTitle[j];
0535 hp_cut->SetTitle(title.c_str());
0536 hp_cut->SetLineColor(kBlack);
0537 hp_cut->DrawClone();
0538 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_ratio_pim_E{}_eta{}.png", (int)E[i], j)).c_str());
0539 c->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_ratio_pim_E{}_eta{}.pdf", (int)E[i], j)).c_str());
0540 c->Clear();
0541
0542
0543 std::vector<std::string> endStr = {";pT [GeV]; Events", ";EDep6/p; Events"};
0544 std::vector<std::string> var_save_loc = {"pT", "EDep6OverP"};
0545 std::vector<std::string> col_loc = {"pT", "EDep6OverP"};
0546 std::vector<std::vector<double>> h1Ranges = {{0, E[i]}, {0, 0.02}};
0547 for (int k = 0; k < 2; k++){
0548 auto cl = new TCanvas("cl", "cl", 700, 500);
0549 title = "e^{-} (E = " + std::to_string((int)E[i]) + " GeV) : " + etaTitle[j] + endStr[k];
0550 auto he1 = e_eta.Histo1D({"he", title.c_str(), 100, h1Ranges[k][0], h1Ranges[k][1]}, col_loc[k]);
0551 he1->SetLineColor(kBlue);
0552 auto he1_cut = e_eta.Filter(currentCut).Histo1D({"he_cut", title.c_str(), 100, h1Ranges[k][0], h1Ranges[k][1]}, col_loc[k]);
0553 he1_cut->GetYaxis()->SetTitleOffset(1.4);
0554 he1_cut->SetLineWidth(2);
0555 he1_cut->SetLineColor(kRed);
0556
0557 title = "#pi^{-} (E = " + std::to_string((int)E[i]) + " GeV) : " + etaTitle[j] + endStr[k];
0558 auto hp1 = p_eta.Histo1D({"hp", title.c_str(), 100, h1Ranges[k][0], h1Ranges[k][1]}, col_loc[k]);
0559 hp1->SetLineColor(kBlue);
0560 auto hp1_cut = p_eta.Filter(currentCut).Histo1D({"hp_cut", title.c_str(), 100, h1Ranges[k][0], h1Ranges[k][1]}, col_loc[k]);
0561 hp1_cut->GetYaxis()->SetTitleOffset(1.4);
0562 hp1_cut->SetLineWidth(2);
0563 hp1_cut->SetLineColor(kRed);
0564
0565 cl->SetLogy(1);
0566 he1->DrawClone();
0567 he1_cut->DrawClone("same");
0568 cl->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_{}_ele_E{}_eta{}.png", col_loc[k], (int)E[i], j)).c_str());
0569 cl->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_{}_ele_E{}_eta{}.pdf", col_loc[k], (int)E[i], j)).c_str());
0570 cl->Clear();
0571
0572 hp1->DrawClone();
0573 hp1_cut->DrawClone("same");
0574 cl->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_{}_pim_E{}_eta{}.png", col_loc[k], (int)E[i], j)).c_str());
0575 cl->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_{}_pim_E{}_eta{}.pdf", col_loc[k], (int)E[i], j)).c_str());
0576 cl->Clear();
0577
0578
0579 title = "#pi^{-}, e^{-} (E = " + std::to_string((int)E[i]) + " GeV) : " + etaTitle[j] + endStr[k];
0580 hp1_cut->SetLineColor(kBlue);
0581 he1_cut->SetLineColor(kRed);
0582 he1_cut->SetTitle(title.c_str());
0583
0584 auto leng = new TLegend(0.7, 0.7, 0.9, 0.9);
0585 he1_cut->DrawClone();
0586 hp1_cut->DrawClone("same");
0587 cl->Update();
0588
0589 leng->AddEntry(he1_cut.GetPtr(),"e^{-}","l");
0590 leng->AddEntry(hp1_cut.GetPtr(),"#pi^{-}","l");
0591 leng->Draw();
0592 cl->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_{}_comb_E{}_eta{}.png", col_loc[k], (int)E[i], j)).c_str());
0593 cl->SaveAs((fmt::format("results/emcal_barrel_pion_rej_cut_{}_comb_E{}_eta{}.pdf", col_loc[k], (int)E[i], j)).c_str());
0594
0595 }
0596 }
0597 }
0598
0599
0600
0601 std::string test_tag = "Barrel_emcal_pion_rejection";
0602
0603 std:string detectorEle = "Barrel_emcal";
0604
0605 for (int i = 0; i < etaTitle.size(); i++){
0606 etaTitle[i].erase(std::remove(etaTitle[i].begin(), etaTitle[i].end(), '#'), etaTitle[i].end());
0607 std::replace(etaTitle[i].begin(), etaTitle[i].end(), 'e', 'E');
0608 }
0609
0610
0611 common_bench::Test pion_rejection_E18_Eta2{
0612 {{"name", fmt::format("{}_E{}_EtaBin{}", test_tag, (int)E[0], 2)},
0613 {"title", "Pion Rejection1"},
0614 {"description", fmt::format("Pion rejection with E = {}, and {}", (int)E[0], etaTitle[2])},
0615 {"quantity", "pi-/e-"},
0616 {"cut", currentCut},
0617 {"e- efficiency", std::to_string(effEle[0][2])},
0618 {"pi- efficiency", std::to_string(effPim[0][2])},
0619 {"target", std::to_string(suppression * maxRate[0][2])}
0620 }
0621 };
0622 suppression * maxRate[0][2] >= rejRatios[0][2] ? pion_rejection_E18_Eta2.pass(rejRatios[0][2]) : pion_rejection_E18_Eta2.fail(rejRatios[0][2]);
0623
0624
0625 common_bench::Test pion_rejection_E18_Eta3{
0626 {{"name", fmt::format("{}_E{}_EtaBin{}", test_tag, (int)E[0], 3)},
0627 {"title", "Pion Rejection"},
0628 {"description", fmt::format("Pion rejection with E = {}, and {}", (int)E[0], etaTitle[3])},
0629 {"quantity", "pi-/e-"},
0630 {"cut", currentCut},
0631 {"e- efficiency", std::to_string(effEle[0][3])},
0632 {"pi- efficiency", std::to_string(effPim[0][3])},
0633 {"target", std::to_string(suppression * maxRate[0][3])}
0634 }
0635 };
0636 suppression * maxRate[0][3] >= rejRatios[0][3] ? pion_rejection_E18_Eta3.pass(rejRatios[0][3]) : pion_rejection_E18_Eta3.fail(rejRatios[0][3]);
0637
0638
0639 common_bench::Test pion_rejection_E10_Eta2{
0640 {{"name", fmt::format("{}_E{}_EtaBin{}", test_tag, (int)E[1], 2)},
0641 {"title", "Pion Rejection"},
0642 {"description", fmt::format("Pion rejection with E = {}, and {}", (int)E[1], etaTitle[2])},
0643 {"quantity", "pi-/e-"},
0644 {"cut", currentCut},
0645 {"e- efficiency", std::to_string(effEle[1][2])},
0646 {"pi- efficiency", std::to_string(effPim[1][2])},
0647 {"target", std::to_string(suppression * maxRate[1][2])}
0648 }
0649 };
0650 suppression * maxRate[1][2] >= rejRatios[1][2] ? pion_rejection_E10_Eta2.pass(rejRatios[1][2]) : pion_rejection_E10_Eta2.fail(rejRatios[1][2]);
0651
0652
0653 common_bench::Test pion_rejection_E10_Eta3{
0654 {{"name", fmt::format("{}_E{}_EtaBin{}", test_tag, (int)E[1], 3)},
0655 {"title", "Pion Rejection"},
0656 {"description", fmt::format("Pion rejection with E = {}, and {}", (int)E[1], etaTitle[3])},
0657 {"quantity", "pi-/e-"},
0658 {"e- efficiency", std::to_string(effEle[1][3])},
0659 {"pi- efficiency", std::to_string(effPim[1][3])},
0660 {"target", std::to_string(suppression * maxRate[1][3])}
0661 }
0662 };
0663 suppression * maxRate[1][3] >= rejRatios[1][3] ? pion_rejection_E10_Eta3.pass(rejRatios[1][3]) : pion_rejection_E10_Eta3.fail(rejRatios[1][3]);
0664
0665
0666 common_bench::Test pion_rejection_E5_Eta2{
0667 {{"name", fmt::format("{}_E{}_EtaBin{}", test_tag, (int)E[2], 2)},
0668 {"title", "Pion Rejection"},
0669 {"description", fmt::format("Pion rejection with E = {}, and {}", (int)E[2], etaTitle[2])},
0670 {"quantity", "pi-/e-"},
0671 {"cut", currentCut},
0672 {"e- efficiency", std::to_string(effEle[2][2])},
0673 {"pi- efficiency", std::to_string(effPim[2][2])},
0674 {"target", std::to_string(suppression * maxRate[2][2])}
0675 }
0676 };
0677 suppression * maxRate[2][2] >= rejRatios[2][2] ? pion_rejection_E5_Eta2.pass(rejRatios[2][2]) : pion_rejection_E5_Eta2.fail(rejRatios[2][2]);
0678
0679
0680 common_bench::Test pion_rejection_E5_Eta3{
0681 {{"name", fmt::format("{}_E{}_EtaBin{}", test_tag, (int)E[2], 3)},
0682 {"title", "Pion Rejection"},
0683 {"description", fmt::format("Pion rejection with E = {}, and {}", (int)E[2], etaTitle[3])},
0684 {"quantity", "pi-/e-"},
0685 {"cut", currentCut},
0686 {"e- efficiency", std::to_string(effEle[2][3])},
0687 {"pi- efficiency", std::to_string(effPim[2][3])},
0688 {"target", std::to_string(suppression * maxRate[2][3])}
0689 }
0690 };
0691 suppression * maxRate[2][3] >= rejRatios[2][3] ? pion_rejection_E5_Eta3.pass(rejRatios[2][3]) : pion_rejection_E5_Eta3.fail(rejRatios[2][3]);
0692
0693
0694 common_bench::write_test({pion_rejection_E18_Eta2,
0695 pion_rejection_E18_Eta3,
0696 pion_rejection_E10_Eta2,
0697 pion_rejection_E10_Eta3,
0698 pion_rejection_E5_Eta2,
0699 pion_rejection_E5_Eta3},
0700 fmt::format("results/{}_pion_rej.json", detectorEle));
0701 }