Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding Hitmap and Average Nhits vs Eta map #113

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Adding Hitmap and Average Nhits vs Eta map #113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
241 changes: 241 additions & 0 deletions benchmarks/tracking_performances/NhitsvsEta_ePIC.C
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,241 @@
// Code to draw average number of hits vs eta at the generated level
// Shyam Kumar; [email protected]; [email protected]
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<Float_t> charge(myReader, "MCParticles.charge");
TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x");
TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y");
TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x");
TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y");
TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z");
TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus");
TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");

TTreeReaderArray<Double_t> *vtx_si_x, *vtx_si_y, *vtx_si_z;
TTreeReaderArray<Double_t> *barrel_si_x, *barrel_si_y, *barrel_si_z;
TTreeReaderArray<Double_t> *disks_si_x, *disks_si_y, *disks_si_z;
TTreeReaderArray<Double_t> *endcap_etof_x, *endcap_etof_y, *endcap_etof_z;
TTreeReaderArray<Double_t> *barrel_mm_x, *barrel_mm_y, *barrel_mm_z;
TTreeReaderArray<Double_t> *barrel_tof_x, *barrel_tof_y, *barrel_tof_z;
TTreeReaderArray<Double_t> *out_mm_x, *out_mm_y, *out_mm_z;
TTreeReaderArray<Double_t> *endcap_fmm_x, *endcap_fmm_y, *endcap_fmm_z;
TTreeReaderArray<Double_t> *endcap_bmm_x, *endcap_bmm_y, *endcap_bmm_z;

// check the quality flag for hits from primary tracks
TTreeReaderArray<Int_t> *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<Double_t>(myReader, "VertexBarrelHits.position.x");
vtx_si_y = new TTreeReaderArray<Double_t>(myReader, "VertexBarrelHits.position.y");
vtx_si_z = new TTreeReaderArray<Double_t>(myReader, "VertexBarrelHits.position.z");
vtx_si_quality = new TTreeReaderArray<Int_t>(myReader, "VertexBarrelHits.quality");

// SVT OB
barrel_si_x = new TTreeReaderArray<Double_t>(myReader, "SiBarrelHits.position.x");
barrel_si_y = new TTreeReaderArray<Double_t>(myReader, "SiBarrelHits.position.y");
barrel_si_z = new TTreeReaderArray<Double_t>(myReader, "SiBarrelHits.position.z");
barrel_si_quality = new TTreeReaderArray<Int_t>(myReader, "SiBarrelHits.quality");

// SVT Disks
disks_si_x = new TTreeReaderArray<Double_t>(myReader, "TrackerEndcapHits.position.x");
disks_si_y = new TTreeReaderArray<Double_t>(myReader, "TrackerEndcapHits.position.y");
disks_si_z = new TTreeReaderArray<Double_t>(myReader, "TrackerEndcapHits.position.z");
disks_si_quality = new TTreeReaderArray<Int_t>(myReader, "TrackerEndcapHits.quality");

// ETOF Hits
endcap_etof_x = new TTreeReaderArray<Double_t>(myReader, "TOFEndcapHits.position.x");
endcap_etof_y = new TTreeReaderArray<Double_t>(myReader, "TOFEndcapHits.position.y");
endcap_etof_z = new TTreeReaderArray<Double_t>(myReader, "TOFEndcapHits.position.z");
endcap_etof_quality = new TTreeReaderArray<Int_t>(myReader, "TOFEndcapHits.quality");

// Inner MPGD
barrel_mm_x = new TTreeReaderArray<Double_t>(myReader, "MPGDBarrelHits.position.x");
barrel_mm_y = new TTreeReaderArray<Double_t>(myReader, "MPGDBarrelHits.position.y");
barrel_mm_z = new TTreeReaderArray<Double_t>(myReader, "MPGDBarrelHits.position.z");
barrel_mm_quality = new TTreeReaderArray<Int_t>(myReader, "MPGDBarrelHits.quality");

// BarrelTOF
barrel_tof_x = new TTreeReaderArray<Double_t>(myReader, "TOFBarrelHits.position.x");
barrel_tof_y = new TTreeReaderArray<Double_t>(myReader, "TOFBarrelHits.position.y");
barrel_tof_z = new TTreeReaderArray<Double_t>(myReader, "TOFBarrelHits.position.z");
barrel_tof_quality = new TTreeReaderArray<Int_t>(myReader, "TOFBarrelHits.quality");

//Outer MPGD Hits
out_mm_x = new TTreeReaderArray<Double_t>(myReader, "OuterMPGDBarrelHits.position.x");
out_mm_y = new TTreeReaderArray<Double_t>(myReader, "OuterMPGDBarrelHits.position.y");
out_mm_z = new TTreeReaderArray<Double_t>(myReader, "OuterMPGDBarrelHits.position.z");
out_mm_quality = new TTreeReaderArray<Int_t>(myReader, "OuterMPGDBarrelHits.quality");

//Forward MPGD
endcap_fmm_x = new TTreeReaderArray<Double_t>(myReader, "ForwardMPGDEndcapHits.position.x");
endcap_fmm_y = new TTreeReaderArray<Double_t>(myReader, "ForwardMPGDEndcapHits.position.y");
endcap_fmm_z = new TTreeReaderArray<Double_t>(myReader, "ForwardMPGDEndcapHits.position.z");
endcap_fmm_quality = new TTreeReaderArray<Int_t>(myReader, "ForwardMPGDEndcapHits.quality");

//Backward MPGD
endcap_bmm_x = new TTreeReaderArray<Double_t>(myReader, "BackwardMPGDEndcapHits.position.x");
endcap_bmm_y = new TTreeReaderArray<Double_t>(myReader, "BackwardMPGDEndcapHits.position.y");
endcap_bmm_z = new TTreeReaderArray<Double_t>(myReader, "BackwardMPGDEndcapHits.position.z");
endcap_bmm_quality = new TTreeReaderArray<Int_t>(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<TVector3> 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));
}


Loading
Loading