File indexing completed on 2026-05-15 07:41:26
0001 #include <cmath>
0002 #include <fstream>
0003 #include <iostream>
0004 #include <string>
0005 #include <vector>
0006 #include <utility>
0007 #include <algorithm>
0008 #include <limits>
0009 #include "TLorentzVector.h"
0010
0011 #include "ROOT/RDataFrame.hxx"
0012 #include <TH1D.h>
0013 #include <TH2D.h>
0014 #include <TRandom3.h>
0015 #include <TFile.h>
0016 #include <TChain.h>
0017 #include <TTree.h>
0018 #include <TMath.h>
0019 #include <TVector3.h>
0020 #include <TCanvas.h>
0021 #include <unordered_map>
0022 #include <unordered_set>
0023 #include <TLorentzVector.h>
0024 #include <TStyle.h>
0025 #include <TLine.h>
0026 #include <TLegend.h>
0027 #include <TGraphErrors.h>
0028 #include <TMultiGraph.h>
0029
0030 #include "TROOT.h"
0031 #include "TRandom.h"
0032 #include "TH3.h"
0033
0034 #include "DD4hep/Detector.h"
0035 #include "DDRec/CellIDPositionConverter.h"
0036 #include "DD4hep/Volumes.h"
0037 #include "DD4hep/Objects.h"
0038 #include "DD4hep/Shapes.h"
0039 #include <DD4hep/Detector.h>
0040 #include <DD4hep/VolumeManager.h>
0041 #include <DD4hep/Printout.h>
0042 #include <TGeoBBox.h>
0043 #include <TGeoMatrix.h>
0044
0045 #include "podio/Frame.h"
0046 #include "podio/CollectionBase.h"
0047 #include "podio/ROOTReader.h"
0048 #include "podio/CollectionIDTable.h"
0049 #include "podio/ObjectID.h"
0050
0051 #include "edm4hep/MCParticleCollection.h"
0052 #include "edm4hep/MCParticle.h"
0053
0054 #include "edm4hep/SimCalorimeterHitCollection.h"
0055 #include "edm4hep/SimCalorimeterHit.h"
0056
0057 #include "edm4hep/CalorimeterHit.h"
0058 #include "edm4hep/CalorimeterHitCollection.h"
0059
0060 #include "edm4eic/ClusterCollection.h"
0061 #include "edm4eic/Cluster.h"
0062 #include "edm4eic/ClusterData.h"
0063
0064 #include "edm4eic/CalorimeterHit.h"
0065 #include "edm4eic/CalorimeterHitCollection.h"
0066
0067 #include "edm4eic/InclusiveKinematicsCollection.h"
0068 #include "edm4eic/InclusiveKinematics.h"
0069 #include "edm4hep/utils/kinematics.h"
0070 #include "edm4hep/utils/vector_utils.h"
0071 #include "edm4eic/vector_utils_legacy.h"
0072
0073 #include "edm4eic/Track.h"
0074 #include "edm4eic/TrackSegment.h"
0075 #include "edm4eic/TrackPoint.h"
0076 #include "edm4eic/TrackParameters.h"
0077 #include "edm4eic/TrackSegmentCollection.h"
0078
0079 #include "edm4eic/ReconstructedParticleCollection.h"
0080 #include "edm4eic/ReconstructedParticle.h"
0081
0082 #include "edm4eic/MCRecoCalorimeterHitAssociationCollection.h"
0083 #include "edm4eic/MCRecoCalorimeterHitAssociation.h"
0084 #include "edm4eic/MCRecoParticleAssociationCollection.h"
0085 #include "edm4eic/MCRecoParticleAssociation.h"
0086
0087 using namespace std;
0088 using namespace ROOT;
0089 using namespace TMath;
0090 using namespace edm4hep;
0091
0092 dd4hep::Detector* det = NULL;
0093
0094 constexpr double ETA_MIN = -4.14, ETA_MAX = -1.16;
0095
0096 inline bool inNHCal(double eta) {return (eta >= ETA_MIN && eta <= ETA_MAX);}
0097
0098 inline string addPrefixAfterSlash(const string& path,
0099 const string& prefix) {
0100 const auto slash = path.find_last_of('/');
0101 if (slash == string::npos)
0102 return prefix + path;
0103 return path.substr(0, slash + 1) + prefix + path.substr(slash + 1);
0104 }
0105
0106 void MakeLogBins(double *array, int nbins, double binLo, double binHi)
0107 {
0108 double logMin = log10(binLo);
0109 double logMax = log10(binHi);
0110 double binWidth = (logMax - logMin) / nbins;
0111
0112 for (int i = 0; i <= nbins; ++i) {
0113 array[i] = pow(10, logMin + i * binWidth);
0114 }
0115 }
0116
0117 TLine* makeLine(double x1, double y1, double x2, double y2) {
0118 TLine* l = new TLine(x1, y1, x2, y2);
0119 l->SetLineColor(kRed);
0120 l->SetLineStyle(2);
0121 l->SetLineWidth(2);
0122 l->Draw("same");
0123 return l;
0124 }
0125
0126 static TGraph* makeEffGraph_B(TH1D* h_sel, TH1D* h_all,
0127 const char* name,
0128 int minAll = 2,
0129 bool logx = false)
0130 {
0131 if (!h_sel || !h_all) return nullptr;
0132
0133 std::unique_ptr<TH1D> h_eff((TH1D*)h_sel->Clone(Form("tmp_%s", name ? name : "eff")));
0134 h_eff->SetDirectory(nullptr);
0135 h_eff->Sumw2();
0136 h_eff->Divide(h_sel, h_all, 1.0, 1.0, "B");
0137
0138 const int nb = h_all->GetNbinsX();
0139
0140 std::unique_ptr<TGraphAsymmErrors> g_log;
0141 std::unique_ptr<TGraphErrors> g_lin;
0142
0143 if (logx) g_log = std::make_unique<TGraphAsymmErrors>();
0144 else g_lin = std::make_unique<TGraphErrors>();
0145
0146 int ip = 0;
0147 for (int b = 1; b <= nb; ++b) {
0148 const double all = h_all->GetBinContent(b);
0149 if (all < minAll) continue;
0150
0151 const double xl = h_all->GetXaxis()->GetBinLowEdge(b);
0152 const double xh = h_all->GetXaxis()->GetBinUpEdge(b);
0153
0154 const double y = 100.0 * h_eff->GetBinContent(b);
0155 const double ey = 100.0 * h_eff->GetBinError(b);
0156
0157 if (logx) {
0158 if (xl <= 0.0 || xh <= 0.0) continue;
0159 const double x = sqrt(xl * xh);
0160 const double exl = x - xl;
0161 const double exh = xh - x;
0162
0163 g_log->SetPoint(ip, x, y);
0164 g_log->SetPointError(ip, exl, exh, ey, ey);
0165 } else {
0166 const double x = h_all->GetBinCenter(b);
0167 const double ex = 0.5 * (xh - xl);
0168
0169 g_lin->SetPoint(ip, x, y);
0170 g_lin->SetPointError(ip, ex, ey);
0171 }
0172 ++ip;
0173 }
0174
0175 TGraph* g = logx ? (TGraph*)g_log.release() : (TGraph*)g_lin.release();
0176 g->SetName(Form("g_%s", name ? name : "eff"));
0177 return g;
0178 }
0179
0180 inline void drawEffPanel(TH1D* h_all,
0181 TH1D* const h_in[3],
0182 const char* title,
0183 const char* xlabel,
0184 bool logx)
0185 {
0186 gPad->SetGrid();
0187 if (logx) gPad->SetLogx();
0188
0189 auto* mg = new TMultiGraph();
0190 auto* leg = new TLegend(0.63, 0.7, 0.88, 0.88);
0191 leg->SetBorderSize(0);
0192 leg->SetFillStyle(0);
0193 leg->SetTextSize(0.04);
0194
0195 Color_t colors[3] = {kBlue, kRed, kBlack};
0196 Style_t markers[3] = {20, 21, 22};
0197
0198 int added = 0;
0199
0200 for (int i = 0; i < 3; ++i) {
0201 if (!h_all || !h_in[i]) continue;
0202 auto* g = makeEffGraph_B(h_in[i], h_all, Form("eff_%d", i), 2, logx);
0203 if (!g || g->GetN() == 0) continue;
0204
0205 g->SetMarkerColor(colors[i]);
0206 g->SetMarkerStyle(markers[i]);
0207 g->SetMarkerSize(1.0);
0208 g->SetLineColor(kBlack);
0209
0210 mg->Add(g, "PE");
0211 leg->AddEntry(g, Form("%d mu-in nHCAL", i), "p");
0212 ++added;
0213 }
0214
0215 if (h_all && h_in[0] && h_in[1] && h_in[2]) {
0216 std::unique_ptr<TH1D> h_sum((TH1D*)h_in[0]->Clone("h_sum"));
0217 h_sum->SetDirectory(nullptr);
0218 h_sum->Add(h_in[1]);
0219 h_sum->Add(h_in[2]);
0220
0221 auto* gsum = makeEffGraph_B(h_sum.get(), h_all, "eff_sum", 2, logx);
0222 if (gsum && gsum->GetN() > 0) {
0223 gsum->SetMarkerStyle(25);
0224 gsum->SetMarkerColor(kMagenta + 1);
0225 gsum->SetFillColor(kMagenta + 1);
0226 gsum->SetLineColor(kBlack);
0227 gsum->SetMarkerSize(1.0);
0228
0229 mg->Add(gsum, "PE");
0230 leg->AddEntry(gsum, "All eligible", "p");
0231 ++added;
0232 }
0233 }
0234
0235 mg->SetTitle(Form("%s;%s;geom. acc. [%%]", title ? title : "", xlabel ? xlabel : ""));
0236 mg->Draw("A");
0237
0238 gPad->Update();
0239 if (mg->GetHistogram()) {
0240 double ymax = mg->GetHistogram()->GetMaximum();
0241 mg->GetHistogram()->SetMaximum(ymax * 1.35);
0242 }
0243 gPad->Modified();
0244 gPad->Update();
0245
0246 if (logx) {
0247 auto* ax = mg->GetXaxis();
0248 if (ax) {
0249 double xmin = ax->GetXmin(), xmax = ax->GetXmax();
0250 if (xmin <= 0) xmin = 1e-12;
0251 mg->GetXaxis()->SetLimits(xmin, xmax);
0252 }
0253 }
0254
0255 leg->Draw();
0256 }
0257
0258 inline double hypot2(double a, double b){ return sqrt(a*a + b*b); }
0259 inline int sgnD(double v){ return (v>0) - (v<0); }
0260
0261 inline pair<double,double>
0262 trackXYatZ_fromTwoStates(double x1,double y1,double z1,
0263 double x2,double y2,double z2,
0264 double zTarget)
0265 {
0266 constexpr double EPSZ = 1e-12;
0267
0268 const double dz = z2 - z1;
0269 if (abs(dz) < EPSZ) {
0270 const double d1 = abs(zTarget - z1);
0271 const double d2 = abs(zTarget - z2);
0272 if (d1 < d2) return {x1, y1};
0273 if (d2 < d1) return {x2, y2};
0274 return {0.5*(x1+x2), 0.5*(y1+y2)};
0275 }
0276
0277 const double t = (zTarget - z1) / dz;
0278 const double x = x1 + t * (x2 - x1);
0279 const double y = y1 + t * (y2 - y1);
0280 return {x, y};
0281 }
0282
0283 inline pair<double,double>
0284 trackXYatZ(const TVector3& A, const TVector3& B, double zTarget){
0285 return trackXYatZ_fromTwoStates(
0286 A.x(), A.y(), A.z(), B.x(), B.y(), B.z(),
0287 zTarget
0288 );
0289 }
0290
0291
0292 int pion_rejection_analysis(const string& filename, string outname_pdf, string outname_png, TString compact_file) {
0293
0294 gStyle->SetOptStat(0);
0295 podio::ROOTReader reader;
0296 reader.openFile(filename);
0297 unsigned nEvents = reader.getEntries("events");
0298 cout << "Number of events: " << nEvents << endl;
0299
0300 det = &(dd4hep::Detector::getInstance());
0301 det->fromCompact(compact_file.Data());
0302
0303 double thickness_plastic_cm = 2.4;
0304 try {
0305 thickness_plastic_cm = det->constantAsDouble("HcalEndcapNPolystyreneThickness");
0306 } catch (...) {
0307 cerr << "[WARN] Using default plastic thickness = " << thickness_plastic_cm << " cm\n";
0308 }
0309
0310 constexpr int NBINS = 60;
0311
0312 constexpr double Q2_MIN = 1e-2;
0313 constexpr double Q2_MAX = 1.0;
0314 constexpr double X_MIN = 1e-6;
0315 constexpr double X_MAX = 1e-3;
0316 constexpr double Y_MIN = 0.0;
0317 constexpr double Y_MAX = 1.0;
0318 constexpr double W_MIN = 0.0;
0319 constexpr double W_MAX = 160.0;
0320
0321 constexpr int Z_NBINS = 100;
0322 constexpr double Z_MIN_MM = -4500.0;
0323 constexpr double Z_MAX_MM = -3600.0;
0324
0325 constexpr int DXY_NBINS = 120;
0326 constexpr double DXY_MIN_MM = -1000.0;
0327 constexpr double DXY_MAX_MM = 1000.0;
0328
0329 constexpr int DR_NBINS = 120;
0330 constexpr double DR_MIN_MM = 0.0;
0331 constexpr double DR_MAX_MM = 400.0;
0332
0333 constexpr double P_MIN_GEV = 0.0;
0334 constexpr double P_MAX_GEV = 25.0;
0335 constexpr int P_NBINS = 50;
0336
0337 constexpr double E_MIN_GEV = 2e-2;
0338 constexpr double E_MAX_GEV = 10.0;
0339
0340 constexpr int SIZE = 3;
0341 constexpr double MIP_ENERGY_GEV = 0.002;
0342 constexpr double E_CUTS[SIZE] = {1.5, 1.7, 2.0};
0343 constexpr double E_THRESH = E_CUTS[2];
0344 constexpr double LAYER_PROC[SIZE] = {0.6, 0.7, 0.8};
0345
0346 double DR_CUTS_CM[SIZE] = {7.0, 10.0, 13.0};
0347 double DR_THRESH_MM = 10*DR_CUTS_CM[2];
0348 int LAYER_MAX = 10;
0349 int LAYER_CUTS[SIZE] = {5, 6, 7};
0350 int LAYER_THRESH = LAYER_CUTS[0];
0351
0352 vector<double> Xbins(NBINS+1), Q2bins(NBINS+1), Ebins(NBINS+1);
0353 MakeLogBins(Q2bins.data(), NBINS, Q2_MIN, Q2_MAX);
0354 MakeLogBins(Xbins.data(), NBINS, X_MIN, X_MAX);
0355 MakeLogBins(Ebins.data(), NBINS, E_MIN_GEV, E_MAX_GEV);
0356
0357 gStyle->SetTitleSize(0.06, "XYZ");
0358 gStyle->SetLabelSize(0.06, "XYZ");
0359 gStyle->SetPadLeftMargin(0.15);
0360 gStyle->SetPadRightMargin(0.15);
0361 gStyle->SetPadBottomMargin(0.15);
0362 gStyle->SetPadTopMargin(0.10);
0363 gROOT->ForceStyle();
0364
0365 TH2D* hEtaPt = new TH2D("hEtaPt", "Pion #eta vs p_{T};#eta;p_{T} [GeV]", 100, -6., 6., 100, 0., 7.);
0366
0367 TH1D* hZ_proj = new TH1D("hZ_proj", "Pion track projections in nHCal; z [mm]; N_{proj}", NBINS, Z_MIN_MM, Z_MAX_MM);
0368
0369 TH1D* hZ_hits = new TH1D("hZ_hits", "Reconstructed hit z in nHCal; z [mm]; N_{hits}", NBINS, Z_MIN_MM, Z_MAX_MM);
0370
0371 TH3D* hDxDyZ_layer = new TH3D("hDxDyZ_layer","3D residuals (rec - proj): dx, dy vs z; dx [mm]; dy [mm]; z [mm]", NBINS, DXY_MIN_MM, DXY_MAX_MM, NBINS, DXY_MIN_MM, DXY_MAX_MM, NBINS, Z_MIN_MM, Z_MAX_MM);
0372
0373 TH2D* hDxDy_all = new TH2D("hDxDy_all",
0374 "Residuals (rec - proj): dx vs dy (all layers); dx [mm]; dy [mm]",
0375 DXY_NBINS, DXY_MIN_MM, DXY_MAX_MM, DXY_NBINS, DXY_MIN_MM, DXY_MAX_MM);
0376
0377 TH2D* hDrZ_layer = new TH2D("hDrZ_layer","Radial residual (rec - proj) vs z; dr [mm]; z [mm]", NBINS, DR_MIN_MM, DR_MAX_MM, NBINS, Z_MIN_MM, Z_MAX_MM);
0378
0379 TH1D* hDr_all = new TH1D("hDr_all",
0380 "Radial residual dr = #sqrt{dx^{2}+dy^{2}} (all layers); dr [mm]; N",
0381 DR_NBINS, DR_MIN_MM, DR_MAX_MM);
0382
0383 TH2D* hE_z = new TH2D("hE_z", "Reconstructed hit energy vs layer z; z_{layer} [mm]; E [GeV]", NBINS, Z_MIN_MM, Z_MAX_MM, NBINS, Ebins.data());
0384
0385 TH2D* hEsum_z = new TH2D("hEsum_z", "Layer energy sum vs z (reconstructed); z_{layer} [mm]; E_{sum} [GeV]",
0386 NBINS, Z_MIN_MM, Z_MAX_MM, NBINS, Ebins.data());
0387
0388 TH1D* hE = new TH1D("hE", "Reconstructed hit energy; E [GeV]; N_{hits}", 10000, 0.0, 1.0);
0389
0390 TH1D* hEsum = new TH1D("hEsum", "Layer energy sum (reconstructed); E_{sum} [GeV]; N_{energy}", 10000, 0.0, 1.0);
0391
0392 TH1D* hP_all_pion = new TH1D("hP_all_pion", "MC pion momentum (pions in nHCal acceptance); p_{MC} [GeV]; N_{pions}", P_NBINS, P_MIN_GEV, P_MAX_GEV);
0393
0394 TH1D* hP_pass_dr[3] = {
0395 new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[0]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0396 new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[1]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0397 new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[2]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0398 };
0399
0400 TH1D* hP_pass_ECut[3] = {
0401 new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[0]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0402 new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[1]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0403 new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[2]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0404 };
0405
0406 TH1D* hP_pass_LayerCut[3] = {
0407 new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[0]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0408 new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[1]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0409 new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[2]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV),
0410 };
0411
0412 for (unsigned ev = 0; ev < nEvents; ev++) {
0413 auto frameData = reader.readNextEntry(podio::Category::Event);
0414 if (!frameData) continue;
0415 podio::Frame frame(std::move(frameData));
0416
0417 auto getCol = [&](const std::string& name) -> const podio::CollectionBase* {
0418 return frame.get(name);
0419 };
0420
0421 const auto* mcColPtr = dynamic_cast<const edm4hep::MCParticleCollection*>(getCol("MCParticles"));
0422 const auto* recPartsPtr = dynamic_cast<const edm4eic::ReconstructedParticleCollection*>(getCol("ReconstructedParticles"));
0423 const auto* projSegsPtr = dynamic_cast<const edm4eic::TrackSegmentCollection*>(getCol("CalorimeterTrackProjections"));
0424 const auto* hcalRecPtr = dynamic_cast<const edm4eic::CalorimeterHitCollection*>(getCol("HcalEndcapNRecHits"));
0425 const auto* assocColPtr = dynamic_cast<const edm4eic::MCRecoParticleAssociationCollection*>(getCol("ReconstructedParticleAssociations"));
0426
0427 if (!mcColPtr || !recPartsPtr || !projSegsPtr || !hcalRecPtr || !assocColPtr) continue;
0428
0429 const auto& mcCol = *mcColPtr;
0430 const auto& recParts = *recPartsPtr;
0431 const auto& projSegs = *projSegsPtr;
0432 const auto& hcalRec = *hcalRecPtr;
0433 const auto& assocCol = *assocColPtr;
0434
0435 vector<edm4hep::MCParticle> vPions;
0436 vector<TLorentzVector> vLorentzPions;
0437 for (const auto& p : mcCol) {
0438 if (p.getPDG() != -211 || p.getGeneratorStatus() == 0) continue;
0439 vPions.push_back(p);
0440 TLorentzVector v(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
0441 vLorentzPions.push_back(v);
0442 hEtaPt->Fill(v.Eta(), v.Pt());
0443 if(inNHCal(v.Eta())) {hP_all_pion->Fill(v.P());}
0444 }
0445
0446 vector<edm4eic::ReconstructedParticle> matchedRecos;
0447 auto find_associated_reco = [&](const edm4hep::MCParticle& mc)->void {
0448 try {
0449 if (!assocCol.isValid() || assocCol.empty()) return;
0450
0451 const uint32_t mc_idx = mc.getObjectID().index;
0452
0453 for (const auto& assoc : assocCol) {
0454 if (assoc.getSimID() == mc_idx) {
0455 uint32_t ridx = assoc.getRecID();
0456 if (!recParts.isValid() || ridx >= recParts.size()) continue;
0457 auto reco = recParts.at(ridx);
0458 if (reco.isAvailable()) {
0459 matchedRecos.push_back(reco);
0460 }
0461 }
0462 }
0463 } catch (...) {}
0464 };
0465
0466 if (!assocCol.isValid() || assocCol.empty()) continue;
0467 for(const auto& p: vPions) if (p.getPDG() == -211) find_associated_reco(p);
0468
0469 if (!recParts.isValid() || !projSegs.isValid() || !hcalRec.isValid() || recParts.empty() || projSegs.empty() || hcalRec.empty()) continue;
0470
0471 set<int> uniqueCentersZ10;
0472 map<double, double> layerData;
0473 for (const auto& hit : hcalRec) {
0474 double z = hit.getPosition().z;
0475 double zBin = round(z);
0476
0477 int z10 = lround(z * 10.0);
0478 uniqueCentersZ10.insert(z10);
0479
0480 layerData[zBin] += hit.getEnergy();
0481 hZ_hits->Fill(z);
0482 hE_z->Fill(z, hit.getEnergy());
0483 hE->Fill(hit.getEnergy());
0484 }
0485
0486 vector<double> layerCentersZ;
0487 layerCentersZ.reserve(uniqueCentersZ10.size());
0488 for (int z10 : uniqueCentersZ10) layerCentersZ.push_back(z10 / 10.0);
0489 if(layerCentersZ.size() > LAYER_MAX) LAYER_MAX = layerCentersZ.size();
0490
0491 for(size_t n = 0; n < SIZE; n++) LAYER_CUTS[n] = static_cast<int>(LAYER_MAX*LAYER_PROC[n]);
0492 LAYER_THRESH = LAYER_CUTS[0];
0493
0494 for (const auto& [zValue, sumEnergy] : layerData) {hEsum_z->Fill(zValue, sumEnergy); hEsum->Fill(sumEnergy);}
0495
0496 vector<edm4eic::Track> allTracks;
0497 for (const auto& R : matchedRecos) {
0498 for (const auto& tr : R.getTracks()) {
0499 if (tr.isAvailable()) allTracks.push_back(tr);
0500 }
0501 }
0502 vector<edm4eic::TrackSegment> segsTagged;
0503 for (const auto& seg : projSegs) {
0504 auto linkedTr = seg.getTrack();
0505 if (!linkedTr.isAvailable()) continue;
0506 for (const auto& TT : allTracks) {
0507 if (linkedTr.getObjectID() == TT.getObjectID()) {
0508 segsTagged.push_back(seg);
0509 break;
0510 }
0511 }
0512 }
0513
0514 for (size_t s = 0; s < segsTagged.size(); ++s) {
0515
0516 vector<double> segMinDistance(LAYER_MAX, std::numeric_limits<double>::infinity());
0517 vector<double> segHitEnergy(LAYER_MAX, std::numeric_limits<double>::quiet_NaN());
0518 vector<int> count_DrCuts(SIZE, 0);
0519 vector<int> count_ECuts(SIZE, 0);
0520
0521 TVector3 A{}, B{};
0522 bool haveA = false, haveB = false;
0523
0524 auto points = segsTagged[s].getPoints();
0525
0526 for (const auto& pt : points) {
0527 if (pt.system != 113) continue;
0528 hZ_proj->Fill(pt.position.z);
0529 if (!haveA) {
0530 A.SetXYZ(pt.position.x, pt.position.y, pt.position.z);
0531 haveA = true;
0532 } else if (!haveB) {
0533 B.SetXYZ(pt.position.x, pt.position.y, pt.position.z);
0534 haveB = true;
0535 break;
0536 }
0537 }
0538
0539 if (!haveA || !haveB) {continue;}
0540
0541 for (size_t i = 0; i < LAYER_MAX; ++i) {
0542 double best_dr_in_layer = 1e10;
0543 double best_E_in_layer;
0544 double partLayerEnergySum = 0;
0545 double ratio_HitPartLayerEnergy = 0;
0546 dd4hep::DDSegmentation::CellID best_cid_in_layer;
0547
0548 double zc = layerCentersZ[i];
0549 auto [X, Y] = trackXYatZ(A, B, zc);
0550
0551 for (const auto& hit : hcalRec) {
0552 const auto& hp = hit.getPosition();
0553
0554 if (fabs(hp.z - zc) > 5.0 ) continue;
0555
0556 const double dx = X - hp.x;
0557 const double dy = Y - hp.y;
0558 const double dr = sqrt(dx*dx + dy*dy);
0559 if(dr < 30*10 ) partLayerEnergySum += hit.getEnergy();
0560
0561 hDxDyZ_layer->Fill(dx, dy, hp.z);
0562 hDxDy_all->Fill(dx, dy);
0563 hDrZ_layer->Fill(dr, hp.z);
0564 hDr_all->Fill(dr);
0565
0566 if (dr < best_dr_in_layer) {
0567 best_dr_in_layer = dr;
0568 best_E_in_layer = hit.getEnergy();
0569 }
0570 }
0571
0572 if (best_dr_in_layer < DR_THRESH_MM) {
0573 segMinDistance[i] = best_dr_in_layer;
0574 segHitEnergy[i] = best_E_in_layer;
0575 }
0576 ratio_HitPartLayerEnergy += segHitEnergy[i]/partLayerEnergySum;
0577 if (std::isnan(ratio_HitPartLayerEnergy) || fabs(ratio_HitPartLayerEnergy - 1) >= 0.2) continue;
0578 for(size_t j = 0; j < SIZE; ++j){
0579 if (segMinDistance[i] < DR_CUTS_CM[j] * 10 ){++count_DrCuts[j];}
0580 if (segHitEnergy[i] < (E_CUTS[j] * MIP_ENERGY_GEV * thickness_plastic_cm * 100)){++count_ECuts[j];}
0581 }
0582 }
0583 for (int j = 0; j < SIZE; ++j){
0584 if (count_DrCuts[j] > LAYER_THRESH) hP_pass_dr[j]->Fill(vLorentzPions[s].P());
0585 if (count_ECuts[j] > LAYER_THRESH) hP_pass_ECut[j]->Fill(vLorentzPions[s].P());
0586 if (count_DrCuts[SIZE-1] > LAYER_CUTS[j]) hP_pass_LayerCut[j]->Fill(vLorentzPions[s].P());
0587 }
0588 }
0589 }
0590
0591 TCanvas* c = new TCanvas("c", "Pion analysis", 1600, 1000);
0592 c->Divide(2,2);
0593
0594 c->cd(1); hEtaPt->Draw("COLZ");
0595
0596 c->SaveAs(outname_pdf.c_str());
0597 c->SaveAs(outname_png.c_str());
0598
0599 TCanvas* c_hitTrack = new TCanvas("c_hitTrack", "Pion hit and track analysis", 1600, 1000);
0600 c_hitTrack->Divide(2,2);
0601 c_hitTrack->cd(1); gPad->SetGrid(); hZ_proj->Draw();
0602 c_hitTrack->cd(2); gPad->SetGrid(); hZ_hits->Draw();
0603 c_hitTrack->SaveAs(addPrefixAfterSlash(outname_png, "hitTrack_").c_str());
0604 c_hitTrack->SaveAs(addPrefixAfterSlash(outname_pdf, "hitTrack_").c_str());
0605
0606 TCanvas* c_hitTrackDistance = new TCanvas("c_hitTrackDistance", "Pion hit-track distance analysis", 1600, 1000);
0607 c_hitTrackDistance->Divide(2,2);
0608 c_hitTrackDistance->cd(1); gPad->SetGrid(); hDxDyZ_layer->Draw("COLZ");
0609 c_hitTrackDistance->cd(2); gPad->SetGrid(); hDxDy_all->Draw("COLZ");
0610 c_hitTrackDistance->cd(3); gPad->SetGrid(); hDrZ_layer->Draw("COLZ");
0611 c_hitTrackDistance->cd(4); gPad->SetGrid(); hDr_all->Draw("COLZ");
0612 c_hitTrackDistance->SaveAs(addPrefixAfterSlash(outname_png, "hitTrackDistance_").c_str());
0613 c_hitTrackDistance->SaveAs(addPrefixAfterSlash(outname_pdf, "hitTrackDistance_").c_str());
0614
0615 TCanvas* c_Eff = new TCanvas("c_Eff", "Pion efficiency analysis", 1600, 1000);
0616 c_Eff->Divide(2,2);
0617 c_Eff->cd(1); gPad->SetGrid(); hP_all_pion->SetLineColor(kBlack); hP_all_pion->Draw();
0618 c_Eff->cd(2); gPad->SetGrid();
0619 TH1D* hEff_dr[3];
0620 for (int idr=0; idr<3; ++idr) {
0621 hEff_dr[idr] = (TH1D*)hP_pass_dr[idr]->Clone(
0622 Form("hEff_dr%.1fcm", DR_CUTS_CM[idr])
0623 );
0624 hEff_dr[idr]->SetDirectory(nullptr);
0625 hEff_dr[idr]->SetTitle(Form("Efficiency vs p_{MC} (Layers > %d); p_{MC} [GeV]; efficiency", LAYER_THRESH));
0626 hEff_dr[idr]->Divide(hP_pass_dr[idr], hP_all_pion, 1, 1, "B");
0627 }
0628
0629 hEff_dr[0]->SetLineColor(kBlue); hEff_dr[0]->SetMinimum(0.0); hEff_dr[0]->SetMaximum(0.3);
0630 hEff_dr[1]->SetLineColor(kRed);
0631 hEff_dr[2]->SetLineColor(kGreen+2);
0632
0633 hEff_dr[0]->Draw("HIST E");
0634 hEff_dr[1]->Draw("HIST E SAME");
0635 hEff_dr[2]->Draw("HIST E SAME");
0636 { auto leg_dr = new TLegend(0.60,0.70,0.88,0.88);
0637 for(int idr=0; idr<3; ++idr)
0638 {leg_dr->AddEntry(hEff_dr[idr],Form("dr < %.1f cm", DR_CUTS_CM[idr]),"l");}
0639 leg_dr->Draw();
0640 }
0641
0642 c_Eff->cd(3); gPad->SetGrid();
0643 TH1D* hEff_E[3];
0644 for (int ie=0; ie<3; ++ie) {
0645 hEff_E[ie] = (TH1D*)hP_pass_ECut[ie]->Clone(
0646 Form("hEff_E%.1fcm", E_CUTS[ie])
0647 );
0648 hEff_E[ie]->SetDirectory(nullptr);
0649 hEff_E[ie]->SetTitle(Form("Efficiency vs p_{MC} (Layers > %d, dr < %.1f cm); p_{MC} [GeV]; efficiency",LAYER_THRESH, DR_THRESH_MM/10));
0650 hEff_E[ie]->Divide(hP_pass_ECut[ie], hP_all_pion, 1, 1, "B");
0651 }
0652
0653 hEff_E[0]->SetLineColor(kBlue); hEff_E[0]->SetMinimum(0.0); hEff_E[0]->SetMaximum(0.3);
0654 hEff_E[1]->SetLineColor(kRed);
0655 hEff_E[2]->SetLineColor(kGreen+2);
0656
0657 hEff_E[0]->Draw("HIST E");
0658 hEff_E[1]->Draw("HIST E SAME");
0659 hEff_E[2]->Draw("HIST E SAME");
0660 { auto leg_E = new TLegend(0.60,0.70,0.88,0.88);
0661 for(int ie=0; ie<3; ++ie)
0662 {leg_E->AddEntry(hEff_E[ie],Form("E < %.1f MIP", E_CUTS[ie]),"l");}
0663 leg_E->Draw();
0664 }
0665
0666 c_Eff->cd(4); gPad->SetGrid();
0667 TH1D* hEff_Layer[3];
0668 for (int il=0; il<3; ++il) {
0669 hEff_Layer[il] = (TH1D*)hP_pass_LayerCut[il]->Clone(
0670 Form("hEff_Layer%d", LAYER_CUTS[il])
0671 );
0672 hEff_Layer[il]->SetDirectory(nullptr);
0673 hEff_Layer[il]->SetTitle(Form("Efficiency vs p_{MC} (dr < %.1f cm); p_{MC} [GeV]; efficiency", DR_THRESH_MM/10));
0674 hEff_Layer[il]->Divide(hP_pass_LayerCut[il], hP_all_pion, 1, 1, "B");
0675 }
0676
0677 hEff_Layer[0]->SetLineColor(kBlue); hEff_Layer[0]->SetMinimum(0.0); hEff_Layer[0]->SetMaximum(0.3);
0678 hEff_Layer[1]->SetLineColor(kRed);
0679 hEff_Layer[2]->SetLineColor(kGreen+2);
0680
0681 hEff_Layer[0]->Draw("HIST E");
0682 hEff_Layer[1]->Draw("HIST E SAME");
0683 hEff_Layer[2]->Draw("HIST E SAME");
0684 { auto leg_Layer = new TLegend(0.60,0.70,0.88,0.88);
0685 for(int il=0; il<3; ++il)
0686 {leg_Layer->AddEntry(hEff_Layer[il],Form("L > %d Layers", LAYER_CUTS[il]),"l");}
0687 leg_Layer->Draw();
0688 }
0689 c_Eff->SaveAs(addPrefixAfterSlash(outname_png, "matching_efficiency_").c_str());
0690 c_Eff->SaveAs(addPrefixAfterSlash(outname_pdf, "matching_efficiency_").c_str());
0691
0692 TCanvas* c_hEnergy = new TCanvas("c_hEnergy", "Pion hit energy analysis", 1600, 1000);
0693 c_hEnergy->Divide(2,2);
0694
0695 c_hEnergy->cd(1); gPad->SetGrid(); gPad->SetLogy(); hE_z->GetYaxis()->SetMoreLogLabels(); hE_z->Draw("COLZ");
0696 TProfile* pE_z = hE_z->ProfileX("pE_z"); pE_z->SetLineWidth(3); pE_z->SetLineColor(kRed); pE_z->SetMarkerColor(kRed);
0697 pE_z->SetDirectory(nullptr); pE_z->Draw("SAME");
0698
0699 const double Ecut_GeV = MIP_ENERGY_GEV * thickness_plastic_cm * 100;
0700 TLine* cut = new TLine(hE_z->GetXaxis()->GetXmin(), Ecut_GeV, hE_z->GetXaxis()->GetXmax(), Ecut_GeV);
0701 cut->SetLineColor(kRed+1); cut->SetLineStyle(2); cut->SetLineWidth(2);
0702 cut->Draw("SAME");
0703 auto* leg = new TLegend(0.58, 0.75, 0.84, 0.88); leg->SetTextSize(0.04);leg->AddEntry(cut, Form("E_{cut} = %.3f GeV", Ecut_GeV), "l");
0704 leg->Draw();
0705
0706 c_hEnergy->cd(2); gPad->SetGrid(); gPad->SetLogy(); hEsum_z->GetYaxis()->SetMoreLogLabels(); hEsum_z->Draw("COLZ");
0707 TProfile* pEsum_z = hEsum_z->ProfileX("pEsum_z"); pEsum_z->SetLineWidth(3); pEsum_z->SetLineColor(kRed); pEsum_z->SetMarkerColor(kRed);
0708 pEsum_z->SetDirectory(nullptr); pEsum_z->Draw("SAME");
0709 c_hEnergy->cd(3); gPad->SetGrid(); hE->Draw();
0710 c_hEnergy->cd(4); gPad->SetGrid(); hEsum->Draw();
0711
0712 c_hEnergy->SaveAs(addPrefixAfterSlash(outname_png, "hit_energy_").c_str());
0713 c_hEnergy->SaveAs(addPrefixAfterSlash(outname_pdf, "hit_energy_").c_str());
0714
0715 delete c;
0716 delete c_hitTrack;
0717 delete c_hitTrackDistance;
0718 delete c_Eff;
0719 delete c_hEnergy;
0720
0721 return 0;
0722 }
