diff --git a/benchmarks/tracking_performances/NhitsvsEta_ePIC.C b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C new file mode 100644 index 00000000..10eabf64 --- /dev/null +++ b/benchmarks/tracking_performances/NhitsvsEta_ePIC.C @@ -0,0 +1,242 @@ +// 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="", TString label="", TString output_prefix=".") + { + + 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("%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('/'); + + 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("%s;#eta_{mc};Nhits",label.Data())); + 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("%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 454abb33..77ef668a 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -45,10 +45,54 @@ 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_hit_maps: + input: + script_hitsmap="benchmarks/tracking_performances/draw_hits.C", + script_nhits_eta="benchmarks/tracking_performances/NhitsvsEta_ePIC.C", + 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.sim_hadd}", "'$OUTPUT_PREFIX'")' +root -l -b -q {input.script_nhits_eta}'("{input.sim_hadd}", "{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? sim=lambda wildcards: expand( diff --git a/benchmarks/tracking_performances/draw_hits.C b/benchmarks/tracking_performances/draw_hits.C new file mode 100644 index 00000000..d632c675 --- /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(TString filename="", TString output_prefix=".") +{ + +//==========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("%s",filename.Data())); + 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(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); + // 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(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(); + Double_t cputime = timer.CpuTime(); + printf("RealTime=%f seconds, CpuTime=%f seconds\n",realtime,cputime); + +}