From cb926e3daa03b721d796a740d1e4bb5d6e584d52 Mon Sep 17 00:00:00 2001 From: Shyam Date: Sat, 14 Dec 2024 16:51:15 +0100 Subject: [PATCH 1/8] adding HitMap and Nhits vs Eta map --- .../tracking_performances/NhitsvsEta_ePIC.C | 241 ++++++++++++++ benchmarks/tracking_performances/draw_hits.C | 313 ++++++++++++++++++ 2 files changed, 554 insertions(+) create mode 100644 benchmarks/tracking_performances/NhitsvsEta_ePIC.C create mode 100644 benchmarks/tracking_performances/draw_hits.C diff --git a/benchmarks/tracking_performances/NhitsvsEta_ePIC.C b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C new file mode 100644 index 00000000..6eaa3395 --- /dev/null +++ b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C @@ -0,0 +1,241 @@ +// Code to draw average number of hits vs eta at the generated level +// Shyam Kumar; shyam055119@gmail.com; shyam.kumar@ba.infn.it +void NhitsvsEta_ePIC(Double_t mom=0.5) + { + + gStyle->SetPalette(1); + gStyle->SetOptTitle(1); + gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(1.0,"Y"); + gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y"); + gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y"); + gStyle->SetHistLineWidth(2); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + + // MC Track Properties + TFile* file = new TFile(Form("./sim%1.1f.edm4hep.root",mom)); // Tree with tracks and hits + TTreeReader myReader("events", file); // name of tree and file + + TTreeReaderArray charge(myReader, "MCParticles.charge"); + TTreeReaderArray vx_mc(myReader, "MCParticles.vertex.x"); + TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); + TTreeReaderArray vz_mc(myReader, "MCParticles.vertex.z"); + TTreeReaderArray px_mc(myReader, "MCParticles.momentum.x"); + TTreeReaderArray py_mc(myReader, "MCParticles.momentum.y"); + TTreeReaderArray pz_mc(myReader, "MCParticles.momentum.z"); + TTreeReaderArray status(myReader, "MCParticles.generatorStatus"); + TTreeReaderArray pdg(myReader, "MCParticles.PDG"); + + TTreeReaderArray *vtx_si_x, *vtx_si_y, *vtx_si_z; + TTreeReaderArray *barrel_si_x, *barrel_si_y, *barrel_si_z; + TTreeReaderArray *disks_si_x, *disks_si_y, *disks_si_z; + TTreeReaderArray *endcap_etof_x, *endcap_etof_y, *endcap_etof_z; + TTreeReaderArray *barrel_mm_x, *barrel_mm_y, *barrel_mm_z; + TTreeReaderArray *barrel_tof_x, *barrel_tof_y, *barrel_tof_z; + TTreeReaderArray *out_mm_x, *out_mm_y, *out_mm_z; + TTreeReaderArray *endcap_fmm_x, *endcap_fmm_y, *endcap_fmm_z; + TTreeReaderArray *endcap_bmm_x, *endcap_bmm_y, *endcap_bmm_z; + + // check the quality flag for hits from primary tracks + TTreeReaderArray *vtx_si_quality, *barrel_si_quality, *disks_si_quality, *endcap_etof_quality, *barrel_mm_quality, *barrel_tof_quality, *out_mm_quality, *endcap_fmm_quality, *endcap_bmm_quality; + + + // Hits on detectors + // SVT IB + vtx_si_x = new TTreeReaderArray(myReader, "VertexBarrelHits.position.x"); + vtx_si_y = new TTreeReaderArray(myReader, "VertexBarrelHits.position.y"); + vtx_si_z = new TTreeReaderArray(myReader, "VertexBarrelHits.position.z"); + vtx_si_quality = new TTreeReaderArray(myReader, "VertexBarrelHits.quality"); + + // SVT OB + barrel_si_x = new TTreeReaderArray(myReader, "SiBarrelHits.position.x"); + barrel_si_y = new TTreeReaderArray(myReader, "SiBarrelHits.position.y"); + barrel_si_z = new TTreeReaderArray(myReader, "SiBarrelHits.position.z"); + barrel_si_quality = new TTreeReaderArray(myReader, "SiBarrelHits.quality"); + + // SVT Disks + disks_si_x = new TTreeReaderArray(myReader, "TrackerEndcapHits.position.x"); + disks_si_y = new TTreeReaderArray(myReader, "TrackerEndcapHits.position.y"); + disks_si_z = new TTreeReaderArray(myReader, "TrackerEndcapHits.position.z"); + disks_si_quality = new TTreeReaderArray(myReader, "TrackerEndcapHits.quality"); + + // ETOF Hits + endcap_etof_x = new TTreeReaderArray(myReader, "TOFEndcapHits.position.x"); + endcap_etof_y = new TTreeReaderArray(myReader, "TOFEndcapHits.position.y"); + endcap_etof_z = new TTreeReaderArray(myReader, "TOFEndcapHits.position.z"); + endcap_etof_quality = new TTreeReaderArray(myReader, "TOFEndcapHits.quality"); + + // Inner MPGD + barrel_mm_x = new TTreeReaderArray(myReader, "MPGDBarrelHits.position.x"); + barrel_mm_y = new TTreeReaderArray(myReader, "MPGDBarrelHits.position.y"); + barrel_mm_z = new TTreeReaderArray(myReader, "MPGDBarrelHits.position.z"); + barrel_mm_quality = new TTreeReaderArray(myReader, "MPGDBarrelHits.quality"); + + // BarrelTOF + barrel_tof_x = new TTreeReaderArray(myReader, "TOFBarrelHits.position.x"); + barrel_tof_y = new TTreeReaderArray(myReader, "TOFBarrelHits.position.y"); + barrel_tof_z = new TTreeReaderArray(myReader, "TOFBarrelHits.position.z"); + barrel_tof_quality = new TTreeReaderArray(myReader, "TOFBarrelHits.quality"); + + //Outer MPGD Hits + out_mm_x = new TTreeReaderArray(myReader, "OuterMPGDBarrelHits.position.x"); + out_mm_y = new TTreeReaderArray(myReader, "OuterMPGDBarrelHits.position.y"); + out_mm_z = new TTreeReaderArray(myReader, "OuterMPGDBarrelHits.position.z"); + out_mm_quality = new TTreeReaderArray(myReader, "OuterMPGDBarrelHits.quality"); + + //Forward MPGD + endcap_fmm_x = new TTreeReaderArray(myReader, "ForwardMPGDEndcapHits.position.x"); + endcap_fmm_y = new TTreeReaderArray(myReader, "ForwardMPGDEndcapHits.position.y"); + endcap_fmm_z = new TTreeReaderArray(myReader, "ForwardMPGDEndcapHits.position.z"); + endcap_fmm_quality = new TTreeReaderArray(myReader, "ForwardMPGDEndcapHits.quality"); + + //Backward MPGD + endcap_bmm_x = new TTreeReaderArray(myReader, "BackwardMPGDEndcapHits.position.x"); + endcap_bmm_y = new TTreeReaderArray(myReader, "BackwardMPGDEndcapHits.position.y"); + endcap_bmm_z = new TTreeReaderArray(myReader, "BackwardMPGDEndcapHits.position.z"); + endcap_bmm_quality = new TTreeReaderArray(myReader, "BackwardMPGDEndcapHits.quality"); + + + Double_t etamc = 0.; + int iEvent=-1; + + TCanvas * c1 = new TCanvas("c1","coutput",1400,1000); + c1->SetMargin(0.10, 0.05 ,0.1,0.07); + c1->SetGridx(); + c1->SetGridy(); + + TProfile* hits = new TProfile("hits","Nhits (#theta)",70,-3.5,3.5); + hits->SetTitle(Form("p = %1.1f GeV/c;#eta_{mc};Nhits",mom)); + hits->GetXaxis()->CenterTitle(); + hits->GetYaxis()->CenterTitle(); + hits->SetMinimum(0.); + + std::vector hitPos; // The ePIC tracker + double epsilon = 1.0e-5; + bool debug = true; + int nhits_SVTIB, nhits_SVTOB, nhits_InMPGD, nhits_BTOF, nhits_OutMPGD, nhits_SVTDisks, nhits_FwdMPGDDisks, nhits_BwdMPGDDisks, nhits_ETOF; + + while (myReader.Next()) + { + hitPos.clear(); + iEvent++; + printf("<--------------------Event No. %d---------------------> \n",iEvent); + // Generated primary track + if (charge.GetSize()>1) continue; // skip event for larger than 1 tracks + Double_t eta_Track = 100.; // set it ouside from -3.5 to 3.5 + Double_t pmc = 0.; + for (int j = 0; j < charge.GetSize(); ++j){ + + if (status[j] !=1) continue; + pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); + Double_t pzmc = pz_mc[j]; + Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/pmc))/2)); + eta_Track = etamc; + Double_t particle = pdg[j]; + } + if (fabs(eta_Track)>3.5) continue; //outside tracker acceptance + // Associated hits with the primary track of momentum (mom) + if (fabs(pmc-mom)> epsilon) continue; + // if (eta_Track<3.4) continue; // Debug for the hits in a given eta range + printf("Eta of the generated track: %f, Momentum (GeV/c): %f \n",eta_Track,pmc); + // ePIC SVT IB Tracker + for (int j = 0; j < vtx_si_x->GetSize(); ++j){ + Double_t xhit = vtx_si_x->At(j); Double_t yhit = vtx_si_y->At(j); Double_t zhit = vtx_si_z->At(j); + Int_t quality = vtx_si_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + nhits_SVTIB = hitPos.size(); hitPos.clear(); + if (debug) printf("SVT IB Associated hits: %d \n",nhits_SVTIB); + + // ePIC SVT OB Tracker + for (int j = 0; j < barrel_si_x->GetSize(); ++j){ + Double_t xhit = barrel_si_x->At(j); Double_t yhit = barrel_si_y->At(j); Double_t zhit = barrel_si_z->At(j); + Int_t quality = barrel_si_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + nhits_SVTOB = hitPos.size(); hitPos.clear(); + if (debug) printf("SVT OB Associated hits: %d \n",nhits_SVTOB); + + // Inner MPGD Tracker + for (int j = 0; j < barrel_mm_x->GetSize(); ++j){ + Double_t xhit = barrel_mm_x->At(j); Double_t yhit = barrel_mm_y->At(j); Double_t zhit = barrel_mm_z->At(j); + Int_t quality = barrel_mm_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + nhits_InMPGD = hitPos.size(); hitPos.clear(); + if (debug) printf("Inner MPGD Associated hits: %d \n",nhits_InMPGD); + + // Outer MPGD Tracker + for (int j = 0; j < out_mm_x->GetSize(); ++j){ + Double_t xhit = out_mm_x->At(j); Double_t yhit = out_mm_y->At(j); Double_t zhit = out_mm_z->At(j); + Int_t quality = out_mm_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + nhits_OutMPGD = hitPos.size(); hitPos.clear(); + if (debug) printf("Outer MPGD Associated hits: %d \n",nhits_OutMPGD); + + // BTOF Tracker + for (int j = 0; j < barrel_tof_x->GetSize(); ++j){ + Double_t xhit = barrel_tof_x->At(j); Double_t yhit = barrel_tof_y->At(j); Double_t zhit = barrel_tof_z->At(j); + Int_t quality = barrel_tof_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + + nhits_BTOF = hitPos.size(); hitPos.clear(); + if (debug) printf("BTOF Associated hits: %d \n",nhits_BTOF); + + // ePIC SVT Disks Hits + for (int j = 0; j < disks_si_x->GetSize(); ++j){ + Int_t quality = disks_si_quality->At(j); + Double_t xhit = disks_si_x->At(j); Double_t yhit = disks_si_y->At(j); Double_t zhit = disks_si_z->At(j); + if (quality==0 && zhit>0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + else if (quality==0 && zhit<0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + + nhits_SVTDisks = hitPos.size(); hitPos.clear(); + if (debug) printf("SVT Disks Associated hits: %d \n",nhits_SVTDisks); + + // ETOF Tracker (Forward) + for (int j = 0; j < endcap_etof_x->GetSize(); ++j){ + Double_t xhit = endcap_etof_x->At(j); Double_t yhit = endcap_etof_y->At(j); Double_t zhit = endcap_etof_z->At(j); + Int_t quality = endcap_etof_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + + nhits_ETOF = hitPos.size(); hitPos.clear(); + if (debug) printf("ETOF Associated hits: %d \n",nhits_ETOF); + + // Forward MPGD + for (int j = 0; j < endcap_fmm_x->GetSize(); ++j){ + Double_t xhit = endcap_fmm_x->At(j); Double_t yhit = endcap_fmm_y->At(j); Double_t zhit = endcap_fmm_z->At(j); + Int_t quality = endcap_fmm_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + nhits_FwdMPGDDisks = hitPos.size(); hitPos.clear(); + if (debug) printf("Forward MPGD Associated hits: %d \n",nhits_FwdMPGDDisks); + + // Backward MPGD + for (int j = 0; j < endcap_bmm_x->GetSize(); ++j){ + Double_t xhit = endcap_bmm_x->At(j); Double_t yhit = endcap_bmm_y->At(j); Double_t zhit = endcap_bmm_z->At(j); + Int_t quality = endcap_bmm_quality->At(j); + if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); + } + + nhits_BwdMPGDDisks = hitPos.size(); hitPos.clear(); + if (debug) printf("Backward MPGD Associated hits: %d \n",nhits_BwdMPGDDisks); + + int nhits = nhits_SVTIB + nhits_SVTOB + nhits_InMPGD + nhits_BTOF + nhits_OutMPGD + nhits_SVTDisks + nhits_FwdMPGDDisks + nhits_BwdMPGDDisks + nhits_ETOF; + if (nhits>0) hits->Fill(eta_Track,nhits); + printf("Total Associated hits: %d \n",nhits); + + } + + c1->cd(); + gPad->SetTicks(1,1); + hits->SetLineWidth(2); + hits->Draw("hist"); + c1->SaveAs(Form("./Output/Nhitsvsmom%1.1f.png",mom)); + c1->SaveAs(Form("./Output/Nhitsvsmom%1.1f.root",mom)); +} + + diff --git a/benchmarks/tracking_performances/draw_hits.C b/benchmarks/tracking_performances/draw_hits.C new file mode 100644 index 00000000..bb2fabc6 --- /dev/null +++ b/benchmarks/tracking_performances/draw_hits.C @@ -0,0 +1,313 @@ +// Hits distribution of detectors +// Shyam Kumar:INFN Bari, shyam.kumar@ba.infn.it; shyam055119@gmail.com + +#include +#include +#include +#include +#include +#include + +void draw_hits(double mom=5.0) +{ + +//==========Style of the plot============ + gStyle->SetPalette(1); + gStyle->SetOptTitle(1); + gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y"); + gStyle->SetTitleSize(.04,"X");gStyle->SetTitleSize(.04,"Y"); + gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y"); + gStyle->SetHistLineWidth(2); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + +//=======Reading the root file DD4HEP=========== + TFile *f = TFile::Open(Form("sim%1.1f.edm4hep.root",mom)); + TTree *sim = (TTree*)f->Get("events"); + + // Timer Start + TStopwatch timer; + timer.Start(); + + TCanvas *c1 = new TCanvas("c1","c1",1200,1000); + c1->SetMargin(0.09, 0.03 ,0.1,0.06); + + // X-Y Hits + Int_t nbins = 320; + Double_t x= 100., y = 100.; + TH2D *hitsxy_vtx_si = new TH2D("hitsxy_vtx_si","hitsxy_vtx_si",nbins,-x,x,nbins,-y,y); + TH2D *hitsxy_barrel_si = new TH2D("hitsxy_barrel_si","hitsxy_barrel_si",nbins,-x,x,nbins,-y,y); + TH2D *hitsxy_barrel_mm_in = new TH2D("hitsxy_barrel_mm_in","hitsxy_barrel_mm_in",nbins,-x,x,nbins,-y,y); + TH2D *hitsxy_barrel_tof = new TH2D("hitsxy_barrel_tof","hitsxy_barrel_tof",nbins,-x,x,nbins,-y,y); + TH2D *hitsxy_barrel_mm_out = new TH2D("hitsxy_barrel_mm_out","hitsxy_barrel_mm_out",nbins,-x,x,nbins,-y,y); + + TString si_vtx_hitsXY ="VertexBarrelHits.position.y*0.1:VertexBarrelHits.position.x*0.1>>hitsxy_vtx_si"; + TString si_barrel_hitsXY ="SiBarrelHits.position.y*0.1:SiBarrelHits.position.x*0.1>>hitsxy_barrel_si"; + TString barrel_mpgd_in_hitsXY ="MPGDBarrelHits.position.y*0.1:MPGDBarrelHits.position.x*0.1>>hitsxy_barrel_mm_in"; + TString tof_barrel_hitsXY ="TOFBarrelHits.position.y*0.1:TOFBarrelHits.position.x*0.1>>hitsxy_barrel_tof"; + TString barrel_mpgd_out_hitsXY ="OuterMPGDBarrelHits.position.y*0.1:OuterMPGDBarrelHits.position.x*0.1>>hitsxy_barrel_mm_out"; + + sim->Draw(si_vtx_hitsXY.Data(),"",""); // Multiply by 0.1 for cm + sim->Draw(si_barrel_hitsXY.Data(),"",""); + sim->Draw(barrel_mpgd_in_hitsXY.Data(),"",""); + sim->Draw(tof_barrel_hitsXY.Data(),"",""); + sim->Draw(barrel_mpgd_out_hitsXY.Data(),"",""); + + hitsxy_vtx_si->SetMarkerStyle(31); + hitsxy_vtx_si->SetTitle("Hit Points (XY)"); + hitsxy_vtx_si->SetMarkerColor(kBlack); + hitsxy_vtx_si->SetMarkerSize(0.1); + hitsxy_vtx_si->SetLineColor(kBlack); + hitsxy_vtx_si->GetXaxis()->SetTitle("X [cm]"); + hitsxy_vtx_si->GetYaxis()->SetTitle("Y [cm]"); + hitsxy_vtx_si->GetXaxis()->CenterTitle(); + hitsxy_vtx_si->GetYaxis()->CenterTitle(); + + hitsxy_barrel_si->SetMarkerStyle(20); + hitsxy_barrel_si->SetMarkerSize(0.1); + hitsxy_barrel_si->SetMarkerColor(kMagenta); + hitsxy_barrel_si->SetLineColor(kMagenta); + + hitsxy_barrel_mm_in->SetMarkerStyle(20); + hitsxy_barrel_mm_in->SetMarkerSize(0.1); + hitsxy_barrel_mm_in->SetMarkerColor(kBlue); + hitsxy_barrel_mm_in->SetLineColor(kBlue); + + hitsxy_barrel_tof->SetMarkerStyle(20); + hitsxy_barrel_tof->SetMarkerSize(0.1); + hitsxy_barrel_tof->SetMarkerColor(kGreen); + hitsxy_barrel_tof->SetLineColor(kGreen);hitsxy_barrel_mm_out->SetMarkerStyle(20); + hitsxy_barrel_mm_out->SetMarkerSize(0.1); + hitsxy_barrel_mm_out->SetMarkerColor(kBlue-7); + hitsxy_barrel_mm_out->SetLineColor(kBlue-7); + + c1->cd(); + hitsxy_vtx_si->Draw(); + hitsxy_barrel_si->Draw("same"); + hitsxy_barrel_mm_in->Draw("same"); + hitsxy_barrel_tof->Draw("same"); + hitsxy_barrel_mm_out->Draw("same"); + + TLegend *l= new TLegend(0.65,0.85,0.90,1.0); + l->SetTextSize(0.025); + l->SetBorderSize(0); + l->AddEntry(hitsxy_vtx_si,"VertexBarrelHits"); + l->AddEntry(hitsxy_barrel_si,"SiBarrelHits"); + l->AddEntry(hitsxy_barrel_mm_in,"MPGDBarrelHits"); + l->AddEntry(hitsxy_barrel_tof,"TOFBarrelHits"); + l->AddEntry(hitsxy_barrel_mm_out,"OuterMPGDBarrelHits"); + l->Draw(); + c1->SaveAs("hitsxy_dd4hep.png"); + c1->SaveAs("hitsxy_dd4hep.eps"); + c1->SaveAs("hitsxy_dd4hep.root"); + + TCanvas *c2 = new TCanvas("c2","c2",1200,1000); + c2->SetMargin(0.09, 0.03 ,0.1,0.06); + // Y-Z Hits + Int_t nbinsx = 400, nbinsy=1800.; + x= 200.; y = 90.; + double xmin = -200.; + + TH2D *hitsrz_vtx_si = new TH2D("hitsrz_vtx_si","hitsrz_vtx_si",nbinsx,xmin,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_si = new TH2D("hitsrz_barrel_si","hitsrz_barrel_si",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_mm_in = new TH2D("hitsrz_barrel_mm_in","hitsrz_barrel_mm_in",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_mm_out = new TH2D("hitsrz_barrel_mm_out","hitsrz_barrel_mm_out",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_tof = new TH2D("hitsrz_barrel_tof","hitsrz_barrel_tof",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_disks_si = new TH2D("hitsrz_disks_si","hitsrz_disks_si",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_endcap_tof = new TH2D("hitsrz_endcap_tof","hitsrz_endcap_tof",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_fwd_mpgd = new TH2D("hitsrz_fwd_mpgd","hitsrz_fwd_mpgd",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_bwd_mpgd = new TH2D("hitsrz_bwd_mpgd","hitsrz_bwd_mpgd",nbinsx,-x,x,nbinsy,-y,y); + + TH2D *hitsrz_vtx_si_1 = new TH2D("hitsrz_vtx_si_1","hitsrz_vtx_si_1",nbinsx,xmin,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_si_1 = new TH2D("hitsrz_barrel_si_1","hitsrz_barrel_si_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_mm_in_1 = new TH2D("hitsrz_barrel_mm_in_1","hitsrz_barrel_mm_in_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_mm_out_1 = new TH2D("hitsrz_barrel_mm_out_1","hitsrz_barrel_mm_out_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_barrel_tof_1 = new TH2D("hitsrz_barrel_tof_1","hitsrz_barrel_tof_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_disks_si_1 = new TH2D("hitsrz_disks_si_1","hitsrz_disks_si_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_endcap_tof_1 = new TH2D("hitsrz_endcap_tof_1","hitsrz_endcap_tof_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_fwd_mpgd_1 = new TH2D("hitsrz_fwd_mpgd_1","hitsrz_fwd_mpgd_1",nbinsx,-x,x,nbinsy,-y,y); + TH2D *hitsrz_bwd_mpgd_1 = new TH2D("hitsrz_bwd_mpgd_1","hitsrz_bwd_mpgd_1",nbinsx,-x,x,nbinsy,-y,y); + + TString si_vtx_hitsrz_posR ="sqrt(VertexBarrelHits.position.x*VertexBarrelHits.position.x+VertexBarrelHits.position.y*VertexBarrelHits.position.y)*0.1:VertexBarrelHits.position.z*0.1>>hitsrz_vtx_si"; + + TString si_barrel_hitsrz_posR ="sqrt(SiBarrelHits.position.x*SiBarrelHits.position.x+SiBarrelHits.position.y*SiBarrelHits.position.y)*0.1:SiBarrelHits.position.z*0.1>>hitsrz_barrel_si"; + + TString barrel_mpgd_in_hitsrz_posR= "sqrt(MPGDBarrelHits.position.x*MPGDBarrelHits.position.x+MPGDBarrelHits.position.y*MPGDBarrelHits.position.y)*0.1:MPGDBarrelHits.position.z*0.1>>hitsrz_barrel_mm_in"; + + TString tof_barrel_hitsrz_posR ="sqrt(TOFBarrelHits.position.x*TOFBarrelHits.position.x+TOFBarrelHits.position.y*TOFBarrelHits.position.y)*0.1:TOFBarrelHits.position.z*0.1>>hitsrz_barrel_tof"; + + TString barrel_mpgd_out_hitsrz_posR ="sqrt(OuterMPGDBarrelHits.position.x*OuterMPGDBarrelHits.position.x+OuterMPGDBarrelHits.position.y*OuterMPGDBarrelHits.position.y)*0.1:OuterMPGDBarrelHits.position.z*0.1>>hitsrz_barrel_mm_out"; + + TString disks_si_hitsrz_posR ="sqrt(TrackerEndcapHits.position.x*TrackerEndcapHits.position.x+TrackerEndcapHits.position.y*TrackerEndcapHits.position.y)*0.1:TrackerEndcapHits.position.z*0.1>>hitsrz_disks_si"; + + TString endcap_tof_hitsrz_posR ="sqrt(TOFEndcapHits.position.x*TOFEndcapHits.position.x+TOFEndcapHits.position.y*TOFEndcapHits.position.y)*0.1:TOFEndcapHits.position.z*0.1>>hitsrz_endcap_tof"; + + TString fwd_mpgd_hitsrz_posR ="sqrt(ForwardMPGDEndcapHits.position.x*ForwardMPGDEndcapHits.position.x+ForwardMPGDEndcapHits.position.y*ForwardMPGDEndcapHits.position.y)*0.1:ForwardMPGDEndcapHits.position.z*0.1>>hitsrz_fwd_mpgd"; + + TString bwd_mpgd_hitsrz_posR ="sqrt(BackwardMPGDEndcapHits.position.x*BackwardMPGDEndcapHits.position.x+BackwardMPGDEndcapHits.position.y*BackwardMPGDEndcapHits.position.y)*0.1:BackwardMPGDEndcapHits.position.z*0.1>>hitsrz_bwd_mpgd"; + + + sim->Draw(si_vtx_hitsrz_posR.Data(),"VertexBarrelHits.position.y>0",""); // Multiply by 0.1 for cm + sim->Draw(si_barrel_hitsrz_posR.Data(),"SiBarrelHits.position.y>0",""); + sim->Draw(barrel_mpgd_in_hitsrz_posR.Data(),"MPGDBarrelHits.position.y>0",""); + sim->Draw(tof_barrel_hitsrz_posR.Data(),"TOFBarrelHits.position.y>0",""); + sim->Draw(disks_si_hitsrz_posR.Data(),"TrackerEndcapHits.position.y>0",""); + sim->Draw(endcap_tof_hitsrz_posR.Data(),"TOFEndcapHits.position.y>0",""); + sim->Draw(barrel_mpgd_out_hitsrz_posR.Data(),"OuterMPGDBarrelHits.position.y>0",""); + sim->Draw(fwd_mpgd_hitsrz_posR.Data(),"ForwardMPGDEndcapHits.position.y>0",""); + sim->Draw(bwd_mpgd_hitsrz_posR.Data(),"BackwardMPGDEndcapHits.position.y>0",""); + + TString si_vtx_hitsrz_negR ="-1.0*sqrt(VertexBarrelHits.position.x*VertexBarrelHits.position.x+VertexBarrelHits.position.y*VertexBarrelHits.position.y)*0.1:VertexBarrelHits.position.z*0.1>>hitsrz_vtx_si_1"; + + TString si_barrel_hitsrz_negR ="-1.0*sqrt(SiBarrelHits.position.x*SiBarrelHits.position.x+SiBarrelHits.position.y*SiBarrelHits.position.y)*0.1:SiBarrelHits.position.z*0.1>>hitsrz_barrel_si_1"; + + TString barrel_mpgd_in_hitsrz_negR= "-1.0*sqrt(MPGDBarrelHits.position.x*MPGDBarrelHits.position.x+MPGDBarrelHits.position.y*MPGDBarrelHits.position.y)*0.1:MPGDBarrelHits.position.z*0.1>>hitsrz_barrel_mm_in_1"; + + TString tof_barrel_hitsrz_negR ="-1.0*sqrt(TOFBarrelHits.position.x*TOFBarrelHits.position.x+TOFBarrelHits.position.y*TOFBarrelHits.position.y)*0.1:TOFBarrelHits.position.z*0.1>>hitsrz_barrel_tof_1"; + + TString barrel_mpgd_out_hitsrz_negR ="-1.0*sqrt(OuterMPGDBarrelHits.position.x*OuterMPGDBarrelHits.position.x+OuterMPGDBarrelHits.position.y*OuterMPGDBarrelHits.position.y)*0.1:OuterMPGDBarrelHits.position.z*0.1>>hitsrz_barrel_mm_out_1"; + + TString disks_si_hitsrz_negR ="-1.0*sqrt(TrackerEndcapHits.position.x*TrackerEndcapHits.position.x+TrackerEndcapHits.position.y*TrackerEndcapHits.position.y)*0.1:TrackerEndcapHits.position.z*0.1>>hitsrz_disks_si_1"; + + TString endcap_tof_hitsrz_negR ="-1.0*sqrt(TOFEndcapHits.position.x*TOFEndcapHits.position.x+TOFEndcapHits.position.y*TOFEndcapHits.position.y)*0.1:TOFEndcapHits.position.z*0.1>>hitsrz_endcap_tof_1"; + + TString fwd_mpgd_hitsrz_negR ="-1.0*sqrt(ForwardMPGDEndcapHits.position.x*ForwardMPGDEndcapHits.position.x+ForwardMPGDEndcapHits.position.y*ForwardMPGDEndcapHits.position.y)*0.1:ForwardMPGDEndcapHits.position.z*0.1>>hitsrz_fwd_mpgd_1"; + + TString bwd_mpgd_hitsrz_negR ="-1.0*sqrt(BackwardMPGDEndcapHits.position.x*BackwardMPGDEndcapHits.position.x+BackwardMPGDEndcapHits.position.y*BackwardMPGDEndcapHits.position.y)*0.1:BackwardMPGDEndcapHits.position.z*0.1>>hitsrz_bwd_mpgd_1"; + + sim->Draw(si_vtx_hitsrz_negR.Data(),"VertexBarrelHits.position.y<0",""); // Multiply by 0.1 for cm + sim->Draw(si_barrel_hitsrz_negR.Data(),"SiBarrelHits.position.y<0",""); + sim->Draw(barrel_mpgd_in_hitsrz_negR.Data(),"MPGDBarrelHits.position.y<0",""); + sim->Draw(tof_barrel_hitsrz_negR.Data(),"TOFBarrelHits.position.y<0",""); + sim->Draw(disks_si_hitsrz_negR.Data(),"TrackerEndcapHits.position.y<0",""); + sim->Draw(endcap_tof_hitsrz_negR.Data(),"TOFEndcapHits.position.y<0",""); + sim->Draw(barrel_mpgd_out_hitsrz_negR.Data(),"OuterMPGDBarrelHits.position.y<0",""); + sim->Draw(fwd_mpgd_hitsrz_negR.Data(),"ForwardMPGDEndcapHits.position.y<0",""); + sim->Draw(bwd_mpgd_hitsrz_negR.Data(),"BackwardMPGDEndcapHits.position.y<0",""); + + + hitsrz_vtx_si->SetMarkerStyle(31); + hitsrz_vtx_si->SetTitle(""); + hitsrz_vtx_si->SetMarkerSize(0.1); + hitsrz_vtx_si->SetLineColor(kBlack); + hitsrz_vtx_si->SetMarkerColor(kBlack); + hitsrz_vtx_si->GetXaxis()->SetTitle("Z [cm]"); + hitsrz_vtx_si->GetYaxis()->SetTitle("R [cm]"); + hitsrz_vtx_si->GetXaxis()->CenterTitle(); + hitsrz_vtx_si->GetYaxis()->CenterTitle(); + + hitsrz_vtx_si_1->SetMarkerSize(0.1); + hitsrz_vtx_si_1->SetLineColor(kBlack); + hitsrz_vtx_si_1->SetMarkerColor(kBlack); + + hitsrz_barrel_si->SetMarkerStyle(20); + hitsrz_barrel_si->SetMarkerSize(0.1); + hitsrz_barrel_si->SetMarkerColor(kMagenta); + hitsrz_barrel_si->SetLineColor(kMagenta); + + hitsrz_barrel_si_1->SetMarkerSize(0.1); + hitsrz_barrel_si_1->SetLineColor(kMagenta); + hitsrz_barrel_si_1->SetMarkerColor(kMagenta); + + hitsrz_barrel_mm_in->SetMarkerStyle(20); + hitsrz_barrel_mm_in->SetMarkerSize(0.1); + hitsrz_barrel_mm_in->SetLineColor(kBlue); + hitsrz_barrel_mm_in->SetMarkerColor(kBlue); + hitsrz_barrel_mm_in_1->SetMarkerSize(0.1); + hitsrz_barrel_mm_in_1->SetLineColor(kBlue); + hitsrz_barrel_mm_in_1->SetMarkerColor(kBlue); + + hitsrz_barrel_tof->SetMarkerStyle(20); + hitsrz_barrel_tof->SetMarkerSize(0.1); + hitsrz_barrel_tof->SetLineColor(kGreen); + hitsrz_barrel_tof->SetMarkerColor(kGreen); + hitsrz_barrel_tof_1->SetMarkerSize(0.1); + hitsrz_barrel_tof_1->SetLineColor(kGreen); + hitsrz_barrel_tof_1->SetMarkerColor(kGreen); + + + hitsrz_barrel_mm_out->SetMarkerStyle(20); + hitsrz_barrel_mm_out->SetMarkerSize(0.1); + hitsrz_barrel_mm_out->SetMarkerColor(kBlue-7); + hitsrz_barrel_mm_out->SetLineColor(kBlue-7); + + hitsrz_barrel_mm_out_1->SetMarkerSize(0.1); + hitsrz_barrel_mm_out_1->SetLineColor(kBlue-7); + hitsrz_barrel_mm_out_1->SetMarkerColor(kBlue-7); + hitsrz_endcap_tof->SetMarkerStyle(20); + hitsrz_endcap_tof->SetMarkerSize(0.1); + hitsrz_endcap_tof->SetMarkerColor(kCyan); + hitsrz_endcap_tof->SetLineColor(kCyan); + + hitsrz_endcap_tof_1->SetMarkerSize(0.1); + hitsrz_endcap_tof_1->SetLineColor(kCyan); + hitsrz_endcap_tof_1->SetMarkerColor(kCyan); + + hitsrz_disks_si->SetMarkerStyle(20); + hitsrz_disks_si->SetMarkerSize(0.1); + hitsrz_disks_si->SetMarkerColor(kRed); + hitsrz_disks_si->SetLineColor(kRed); + + hitsrz_disks_si_1->SetMarkerSize(0.1); + hitsrz_disks_si_1->SetLineColor(kRed); + hitsrz_disks_si_1->SetMarkerColor(kRed); + + hitsrz_fwd_mpgd->SetMarkerSize(0.1); + hitsrz_fwd_mpgd->SetLineColor(kRed-7); + hitsrz_fwd_mpgd->SetMarkerColor(kRed-7); + hitsrz_fwd_mpgd_1->SetMarkerSize(0.1); + hitsrz_fwd_mpgd_1->SetLineColor(kRed-7); + hitsrz_fwd_mpgd_1->SetMarkerColor(kRed-7); + hitsrz_bwd_mpgd->SetMarkerSize(0.1); + hitsrz_bwd_mpgd->SetLineColor(kOrange); + hitsrz_bwd_mpgd->SetMarkerColor(kOrange); + + hitsrz_bwd_mpgd_1->SetMarkerSize(0.1); + hitsrz_bwd_mpgd_1->SetLineColor(kOrange); + hitsrz_bwd_mpgd_1->SetMarkerColor(kOrange); + + c2->cd(); + hitsrz_vtx_si->Draw(""); + hitsrz_vtx_si_1->Draw("same"); + hitsrz_barrel_si->Draw("same"); + hitsrz_barrel_si_1->Draw("same"); + hitsrz_barrel_mm_in->Draw("same"); + hitsrz_barrel_mm_in_1->Draw("same"); + hitsrz_barrel_tof->Draw("same"); + hitsrz_barrel_tof_1->Draw("same"); + hitsrz_barrel_mm_out->Draw("same"); + hitsrz_barrel_mm_out_1->Draw("same"); + hitsrz_endcap_tof->Draw("same"); + hitsrz_endcap_tof_1->Draw("same"); + hitsrz_disks_si->Draw("same"); + hitsrz_disks_si_1->Draw("same"); + hitsrz_fwd_mpgd->Draw("same"); + hitsrz_fwd_mpgd_1->Draw("same"); + hitsrz_bwd_mpgd->Draw("same"); + hitsrz_bwd_mpgd_1->Draw("same"); + + TLegend *l1= new TLegend(0.11,0.88,0.95,0.99); + l1->SetNColumns(3); + l1->SetTextSize(0.025); + l1->SetBorderSize(0); + l1->AddEntry(hitsrz_vtx_si,"VertexBarrelHits"); + l1->AddEntry(hitsrz_barrel_si,"SiBarrelHits"); + l1->AddEntry(hitsrz_barrel_mm_in,"MPGDBarrelHits"); + l1->AddEntry(hitsrz_barrel_tof,"TOFBarrelHits"); + l1->AddEntry(hitsrz_barrel_mm_out,"OuterMPGDBarrelHits"); + l1->AddEntry(hitsrz_disks_si,"TrackerEndcapHits"); + l1->AddEntry(hitsrz_endcap_tof,"TOFEndcapHits"); + l1->AddEntry(hitsrz_fwd_mpgd,"ForwardMPGDEndcapHits"); + l1->AddEntry(hitsrz_bwd_mpgd,"BackwardMPGDEndcapHits"); + l1->Draw(); + + c2->SaveAs("hitsrz_dd4hep.png"); + c2->SaveAs("hitsrz_dd4hep.eps"); + c2->SaveAs("hitsrz_dd4hep.root"); + // Timer Stop + timer.Stop(); + Double_t realtime = timer.RealTime(); + Double_t cputime = timer.CpuTime(); + printf("RealTime=%f seconds, CpuTime=%f seconds\n",realtime,cputime); + +} From f16accf4b80f69e1b45f9d3615ac33117d00ea08 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Thu, 26 Dec 2024 12:30:54 +0100 Subject: [PATCH 2/8] update code with snakefile --- benchmarks/tracking_performances/draw_hits.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/tracking_performances/draw_hits.C b/benchmarks/tracking_performances/draw_hits.C index bb2fabc6..3f86cf21 100644 --- a/benchmarks/tracking_performances/draw_hits.C +++ b/benchmarks/tracking_performances/draw_hits.C @@ -8,7 +8,7 @@ #include #include -void draw_hits(double mom=5.0) +void draw_hits(TString filename="") { //==========Style of the plot============ @@ -22,7 +22,7 @@ void draw_hits(double mom=5.0) gStyle->SetOptStat(0); //=======Reading the root file DD4HEP=========== - TFile *f = TFile::Open(Form("sim%1.1f.edm4hep.root",mom)); + TFile *f = TFile::Open(Form("%s",filename.Data())); TTree *sim = (TTree*)f->Get("events"); // Timer Start From df57e14d530c3fb4bbc8b8fda9ca353e5e02c393 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Thu, 26 Dec 2024 12:31:38 +0100 Subject: [PATCH 3/8] update code --- .../tracking_performances/NhitsvsEta_ePIC.C | 24 ++++++++++++------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/benchmarks/tracking_performances/NhitsvsEta_ePIC.C b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C index 6eaa3395..ce4be54b 100644 --- a/benchmarks/tracking_performances/NhitsvsEta_ePIC.C +++ b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C @@ -1,6 +1,6 @@ // Code to draw average number of hits vs eta at the generated level // Shyam Kumar; shyam055119@gmail.com; shyam.kumar@ba.infn.it -void NhitsvsEta_ePIC(Double_t mom=0.5) +void NhitsvsEta_ePIC(TString filePath="") { gStyle->SetPalette(1); @@ -13,9 +13,14 @@ void NhitsvsEta_ePIC(Double_t mom=0.5) gStyle->SetOptStat(0); // MC Track Properties - TFile* file = new TFile(Form("./sim%1.1f.edm4hep.root",mom)); // Tree with tracks and hits + TFile* file = new TFile(Form("%s",filePath.Data())); // Tree with tracks and hits TTreeReader myReader("events", file); // name of tree and file - + // Find the last occurrence of '/' + Int_t lastSlashPos = filePath.Last('/'); + + // Extract the file name + TString filename = (lastSlashPos != kNPOS) ? filePath(lastSlashPos + 1, filePath.Length()) : filePath; + TTreeReaderArray charge(myReader, "MCParticles.charge"); TTreeReaderArray vx_mc(myReader, "MCParticles.vertex.x"); TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); @@ -103,9 +108,11 @@ void NhitsvsEta_ePIC(Double_t mom=0.5) c1->SetMargin(0.10, 0.05 ,0.1,0.07); c1->SetGridx(); c1->SetGridy(); - + + filename.Resize(filename.Sizeof()-14); // keep filename up to energy + TProfile* hits = new TProfile("hits","Nhits (#theta)",70,-3.5,3.5); - hits->SetTitle(Form("p = %1.1f GeV/c;#eta_{mc};Nhits",mom)); + hits->SetTitle(Form("%s;#eta_{mc};Nhits",filename.Data())); hits->GetXaxis()->CenterTitle(); hits->GetYaxis()->CenterTitle(); hits->SetMinimum(0.); @@ -135,7 +142,7 @@ void NhitsvsEta_ePIC(Double_t mom=0.5) } if (fabs(eta_Track)>3.5) continue; //outside tracker acceptance // Associated hits with the primary track of momentum (mom) - if (fabs(pmc-mom)> epsilon) continue; + // if (fabs(pmc-mom)> epsilon) continue; // if (eta_Track<3.4) continue; // Debug for the hits in a given eta range printf("Eta of the generated track: %f, Momentum (GeV/c): %f \n",eta_Track,pmc); // ePIC SVT IB Tracker @@ -234,8 +241,7 @@ void NhitsvsEta_ePIC(Double_t mom=0.5) gPad->SetTicks(1,1); hits->SetLineWidth(2); hits->Draw("hist"); - c1->SaveAs(Form("./Output/Nhitsvsmom%1.1f.png",mom)); - c1->SaveAs(Form("./Output/Nhitsvsmom%1.1f.root",mom)); + c1->SaveAs(Form("%s.png",filename.Data())); + c1->SaveAs(Form("Nhitsvsmom%s.root",filename.Data())); } - From d743aa9ee6f0eaeda7f8fe93b9d245b7fd292a10 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Thu, 26 Dec 2024 12:33:18 +0100 Subject: [PATCH 4/8] fixed indentation --- benchmarks/tracking_performances/draw_hits.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/tracking_performances/draw_hits.C b/benchmarks/tracking_performances/draw_hits.C index 3f86cf21..abf5c1b5 100644 --- a/benchmarks/tracking_performances/draw_hits.C +++ b/benchmarks/tracking_performances/draw_hits.C @@ -153,8 +153,8 @@ void draw_hits(TString filename="") sim->Draw(tof_barrel_hitsrz_posR.Data(),"TOFBarrelHits.position.y>0",""); sim->Draw(disks_si_hitsrz_posR.Data(),"TrackerEndcapHits.position.y>0",""); sim->Draw(endcap_tof_hitsrz_posR.Data(),"TOFEndcapHits.position.y>0",""); - sim->Draw(barrel_mpgd_out_hitsrz_posR.Data(),"OuterMPGDBarrelHits.position.y>0",""); - sim->Draw(fwd_mpgd_hitsrz_posR.Data(),"ForwardMPGDEndcapHits.position.y>0",""); + sim->Draw(barrel_mpgd_out_hitsrz_posR.Data(),"OuterMPGDBarrelHits.position.y>0",""); + sim->Draw(fwd_mpgd_hitsrz_posR.Data(),"ForwardMPGDEndcapHits.position.y>0",""); sim->Draw(bwd_mpgd_hitsrz_posR.Data(),"BackwardMPGDEndcapHits.position.y>0",""); TString si_vtx_hitsrz_negR ="-1.0*sqrt(VertexBarrelHits.position.x*VertexBarrelHits.position.x+VertexBarrelHits.position.y*VertexBarrelHits.position.y)*0.1:VertexBarrelHits.position.z*0.1>>hitsrz_vtx_si_1"; From 300cfb8bb514a695ed4369ce664f980bc0a04f35 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Thu, 26 Dec 2024 12:35:30 +0100 Subject: [PATCH 5/8] Updated Snakefile with simulation performances --- benchmarks/tracking_performances/Snakefile | 27 ++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index 454abb33..a810dd6d 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -45,10 +45,35 @@ exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ -Ppodio:output_collections=MCParticles,CentralCKFTrajectories,CentralCKFTrackParameters,CentralCKFSeededTrackParameters,CentralCKFTruthSeededTrackParameters,CentralTrackVertices """ +rule tracking_performance_sim_hadd: + input: + simoutput=lambda wildcards: + expand( + "sim_output/tracking_performance/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.edm4hep.root", + DETECTOR_CONFIG="epic_craterlake_tracking_only", + PARTICLE=wildcards.PARTICLE, + ENERGY=wildcards.ENERGY, + PHASE_SPACE=["3to50deg", "45to135deg", "130to177deg"], + INDEX=range(1), + ) + output: + "sim_output/tracking_performance/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PARTICLE}.{ENERGY}.edm4hep.root", + shell: + """ +hadd {output} {input.simoutput} +""" rule tracking_performance_at_momentum: input: script="benchmarks/tracking_performances/Tracking_Performances.C", + script_hitsmap="benchmarks/tracking_performances/draw_hits.C", + script_nhits_eta="benchmarks/tracking_performances/NhitsvsEta_ePIC.C", + outsim=lambda wildcards: + expand( + "sim_output/tracking_performance/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PARTICLE}.{ENERGY}.edm4hep.root", + DETECTOR_CONFIG="epic_craterlake_tracking_only", PARTICLE=wildcards.PARTICLE, + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV") + if wildcards.CAMPAIGN == "local" else None, # TODO pass as a file list? sim=lambda wildcards: expand( @@ -80,6 +105,8 @@ fi hadd {output.combined_root} {input.sim} cd {wildcards.CAMPAIGN} root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15, '$TRUTH_SEEDING')' +root -l -b -q ../{input.script_hitsmap}'("../{input.outsim}")' +root -l -b -q ../{input.script_nhits_eta}'("../{input.outsim}")' """ From 8cce96071d609fa8b2f7721951a2a42ecd3a0b92 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 27 Dec 2024 00:03:45 -0500 Subject: [PATCH 6/8] fix indentation, fix eval --- benchmarks/tracking_performances/Snakefile | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index a810dd6d..ff214186 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -68,12 +68,14 @@ rule tracking_performance_at_momentum: script="benchmarks/tracking_performances/Tracking_Performances.C", script_hitsmap="benchmarks/tracking_performances/draw_hits.C", script_nhits_eta="benchmarks/tracking_performances/NhitsvsEta_ePIC.C", - outsim=lambda wildcards: + outsim=lambda wildcards: expand( "sim_output/tracking_performance/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PARTICLE}.{ENERGY}.edm4hep.root", DETECTOR_CONFIG="epic_craterlake_tracking_only", PARTICLE=wildcards.PARTICLE, - ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV") - if wildcards.CAMPAIGN == "local" else None, + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + ) + if wildcards.CAMPAIGN == "local" else + [], # TODO pass as a file list? sim=lambda wildcards: expand( From 21711f9f74e3df75e43e63b2f054f9b4ad18c4f3 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 27 Dec 2024 00:43:08 -0500 Subject: [PATCH 7/8] move to a separate rule, avoid overwriting files --- .../tracking_performances/NhitsvsEta_ePIC.C | 13 +++---- benchmarks/tracking_performances/Snakefile | 35 +++++++++++++------ benchmarks/tracking_performances/draw_hits.C | 14 ++++---- 3 files changed, 36 insertions(+), 26 deletions(-) diff --git a/benchmarks/tracking_performances/NhitsvsEta_ePIC.C b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C index ce4be54b..10eabf64 100644 --- a/benchmarks/tracking_performances/NhitsvsEta_ePIC.C +++ b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C @@ -1,6 +1,6 @@ // Code to draw average number of hits vs eta at the generated level // Shyam Kumar; shyam055119@gmail.com; shyam.kumar@ba.infn.it -void NhitsvsEta_ePIC(TString filePath="") +void NhitsvsEta_ePIC(TString filePath="", TString label="", TString output_prefix=".") { gStyle->SetPalette(1); @@ -18,9 +18,6 @@ void NhitsvsEta_ePIC(TString filePath="") // Find the last occurrence of '/' Int_t lastSlashPos = filePath.Last('/'); - // Extract the file name - TString filename = (lastSlashPos != kNPOS) ? filePath(lastSlashPos + 1, filePath.Length()) : filePath; - TTreeReaderArray charge(myReader, "MCParticles.charge"); TTreeReaderArray vx_mc(myReader, "MCParticles.vertex.x"); TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); @@ -109,10 +106,8 @@ void NhitsvsEta_ePIC(TString filePath="") c1->SetGridx(); c1->SetGridy(); - filename.Resize(filename.Sizeof()-14); // keep filename up to energy - TProfile* hits = new TProfile("hits","Nhits (#theta)",70,-3.5,3.5); - hits->SetTitle(Form("%s;#eta_{mc};Nhits",filename.Data())); + hits->SetTitle(Form("%s;#eta_{mc};Nhits",label.Data())); hits->GetXaxis()->CenterTitle(); hits->GetYaxis()->CenterTitle(); hits->SetMinimum(0.); @@ -241,7 +236,7 @@ void NhitsvsEta_ePIC(TString filePath="") gPad->SetTicks(1,1); hits->SetLineWidth(2); hits->Draw("hist"); - c1->SaveAs(Form("%s.png",filename.Data())); - c1->SaveAs(Form("Nhitsvsmom%s.root",filename.Data())); + c1->SaveAs(Form("%s/Nhits_vs_eta.png", output_prefix.Data())); + c1->SaveAs(Form("%s/Nhits_vs_eta.root", output_prefix.Data())); } diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index ff214186..d9bb20c9 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -63,17 +63,34 @@ rule tracking_performance_sim_hadd: hadd {output} {input.simoutput} """ -rule tracking_performance_at_momentum: +rule tracking_performance_hit_maps: input: - script="benchmarks/tracking_performances/Tracking_Performances.C", script_hitsmap="benchmarks/tracking_performances/draw_hits.C", script_nhits_eta="benchmarks/tracking_performances/NhitsvsEta_ePIC.C", - outsim=lambda wildcards: - expand( - "sim_output/tracking_performance/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PARTICLE}.{ENERGY}.edm4hep.root", - DETECTOR_CONFIG="epic_craterlake_tracking_only", PARTICLE=wildcards.PARTICLE, - ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", - ) + sim_hadd=expand( + "sim_output/tracking_performance/{DETECTOR_CONFIG}/{{PARTICLE}}/{{MOMENTUM}}/{{PARTICLE}}.{{MOMENTUM}}.edm4hep.root", + DETECTOR_CONFIG="epic_craterlake_tracking_only", + ) + output: + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/hitsxy_dd4hep.png", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/hitsxy_dd4hep.eps", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/hitsxy_dd4hep.root", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/hitsrz_dd4hep.png", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/hitsrz_dd4hep.eps", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/hitsrz_dd4hep.root", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/Nhits_vs_eta.png", + "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/Nhits_vs_eta.root", + shell: """ +OUTPUT_PREFIX="$(dirname {output[0]})" +root -l -b -q {input.script_hitsmap}'("{input.outsim}", "'$OUTPUT_PREFIX'")' +root -l -b -q {input.script_nhits_eta}'("{input.outsim}", "{wildcards.MOMENTUM}", "'$OUTPUT_PREFIX'")' +""" + +rule tracking_performance_at_momentum: + input: + script="benchmarks/tracking_performances/Tracking_Performances.C", + hit_maps=lambda wildcards: + [ "{CAMPAIGN}/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/Nhits_vs_eta.png".format(**wildcards) ] if wildcards.CAMPAIGN == "local" else [], # TODO pass as a file list? @@ -107,8 +124,6 @@ fi hadd {output.combined_root} {input.sim} cd {wildcards.CAMPAIGN} root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15, '$TRUTH_SEEDING')' -root -l -b -q ../{input.script_hitsmap}'("../{input.outsim}")' -root -l -b -q ../{input.script_nhits_eta}'("../{input.outsim}")' """ diff --git a/benchmarks/tracking_performances/draw_hits.C b/benchmarks/tracking_performances/draw_hits.C index abf5c1b5..d632c675 100644 --- a/benchmarks/tracking_performances/draw_hits.C +++ b/benchmarks/tracking_performances/draw_hits.C @@ -8,7 +8,7 @@ #include #include -void draw_hits(TString filename="") +void draw_hits(TString filename="", TString output_prefix=".") { //==========Style of the plot============ @@ -97,9 +97,9 @@ void draw_hits(TString filename="") l->AddEntry(hitsxy_barrel_tof,"TOFBarrelHits"); l->AddEntry(hitsxy_barrel_mm_out,"OuterMPGDBarrelHits"); l->Draw(); - c1->SaveAs("hitsxy_dd4hep.png"); - c1->SaveAs("hitsxy_dd4hep.eps"); - c1->SaveAs("hitsxy_dd4hep.root"); + c1->SaveAs(Form("%s/hitsxy_dd4hep.png", output_prefix.Data())); + c1->SaveAs(Form("%s/hitsxy_dd4hep.eps", output_prefix.Data())); + c1->SaveAs(Form("%s/hitsxy_dd4hep.root", output_prefix.Data())); TCanvas *c2 = new TCanvas("c2","c2",1200,1000); c2->SetMargin(0.09, 0.03 ,0.1,0.06); @@ -301,9 +301,9 @@ void draw_hits(TString filename="") l1->AddEntry(hitsrz_bwd_mpgd,"BackwardMPGDEndcapHits"); l1->Draw(); - c2->SaveAs("hitsrz_dd4hep.png"); - c2->SaveAs("hitsrz_dd4hep.eps"); - c2->SaveAs("hitsrz_dd4hep.root"); + c2->SaveAs(Form("%s/hitsrz_dd4hep.png", output_prefix.Data())); + c2->SaveAs(Form("%s/hitsrz_dd4hep.eps", output_prefix.Data())); + c2->SaveAs(Form("%s/hitsrz_dd4hep.root", output_prefix.Data())); // Timer Stop timer.Stop(); Double_t realtime = timer.RealTime(); From daad459cd5902941184b33a6c400bb5860189af9 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Fri, 27 Dec 2024 08:24:57 +0100 Subject: [PATCH 8/8] Correcting input outsim-->sim_hadd --- benchmarks/tracking_performances/Snakefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index d9bb20c9..77ef668a 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -82,8 +82,8 @@ rule tracking_performance_hit_maps: "local/{SEEDING}/pi-/{MOMENTUM}/{SEEDING_IGNORE}/{PARTICLE}/Nhits_vs_eta.root", shell: """ OUTPUT_PREFIX="$(dirname {output[0]})" -root -l -b -q {input.script_hitsmap}'("{input.outsim}", "'$OUTPUT_PREFIX'")' -root -l -b -q {input.script_nhits_eta}'("{input.outsim}", "{wildcards.MOMENTUM}", "'$OUTPUT_PREFIX'")' +root -l -b -q {input.script_hitsmap}'("{input.sim_hadd}", "'$OUTPUT_PREFIX'")' +root -l -b -q {input.script_nhits_eta}'("{input.sim_hadd}", "{wildcards.MOMENTUM}", "'$OUTPUT_PREFIX'")' """ rule tracking_performance_at_momentum: