Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-01 07:56:16 

0001 // Code to draw average number of hits vs eta at the generated level
0002 // Shyam Kumar; shyam055119@gmail.com; shyam.kumar@ba.infn.it
0003 void NhitsvsEta_ePIC(TString filePath="", TString label="", TString output_prefix=".")
0004   {
0005   
0006    gStyle->SetPalette(1);
0007    gStyle->SetOptTitle(1);
0008    gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(1.0,"Y");
0009    gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
0010    gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
0011    gStyle->SetHistLineWidth(2);
0012    gStyle->SetTitleAlign(23);     
0013    gStyle->SetOptFit(1);
0014    gStyle->SetOptStat(0);
0015    
0016      // MC Track Properties
0017     TFile* file = new TFile(Form("%s",filePath.Data())); // Tree with tracks and hits
0018     TTreeReader myReader("events", file); // name of tree and file
0019     // Find the last occurrence of '/'
0020     Int_t lastSlashPos = filePath.Last('/');
0021 
0022    TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge"); 
0023    TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x"); 
0024    TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y"); 
0025    TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z"); 
0026    TTreeReaderArray<Double_t> px_mc(myReader, "MCParticles.momentum.x"); 
0027    TTreeReaderArray<Double_t> py_mc(myReader, "MCParticles.momentum.y"); 
0028    TTreeReaderArray<Double_t> pz_mc(myReader, "MCParticles.momentum.z"); 
0029    TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus"); 
0030    TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG"); 
0031 
0032   TTreeReaderArray<Double_t> *vtx_si_x, *vtx_si_y, *vtx_si_z;
0033   TTreeReaderArray<Double_t> *barrel_si_x, *barrel_si_y, *barrel_si_z; 
0034   TTreeReaderArray<Double_t> *disks_si_x, *disks_si_y, *disks_si_z; 
0035   TTreeReaderArray<Double_t> *endcap_etof_x, *endcap_etof_y, *endcap_etof_z;
0036   TTreeReaderArray<Double_t> *barrel_mm_x, *barrel_mm_y, *barrel_mm_z;
0037   TTreeReaderArray<Double_t> *barrel_tof_x, *barrel_tof_y, *barrel_tof_z;
0038   TTreeReaderArray<Double_t> *out_mm_x, *out_mm_y, *out_mm_z; 
0039   TTreeReaderArray<Double_t> *endcap_fmm_x, *endcap_fmm_y, *endcap_fmm_z; 
0040   TTreeReaderArray<Double_t> *endcap_bmm_x, *endcap_bmm_y, *endcap_bmm_z; 
0041 
0042   // check the quality flag for hits from primary tracks
0043   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;    
0044 
0045 
0046    // Hits on detectors
0047    // SVT IB
0048    vtx_si_x = new TTreeReaderArray<Double_t>(myReader, "VertexBarrelHits.position.x"); 
0049    vtx_si_y = new TTreeReaderArray<Double_t>(myReader, "VertexBarrelHits.position.y"); 
0050    vtx_si_z = new TTreeReaderArray<Double_t>(myReader, "VertexBarrelHits.position.z"); 
0051    vtx_si_quality = new TTreeReaderArray<Int_t>(myReader, "VertexBarrelHits.quality");
0052    
0053    // SVT OB
0054    barrel_si_x = new TTreeReaderArray<Double_t>(myReader, "SiBarrelHits.position.x"); 
0055    barrel_si_y = new TTreeReaderArray<Double_t>(myReader, "SiBarrelHits.position.y"); 
0056    barrel_si_z = new TTreeReaderArray<Double_t>(myReader, "SiBarrelHits.position.z");
0057    barrel_si_quality = new TTreeReaderArray<Int_t>(myReader, "SiBarrelHits.quality");
0058    
0059    // SVT Disks
0060    disks_si_x = new TTreeReaderArray<Double_t>(myReader, "TrackerEndcapHits.position.x"); 
0061    disks_si_y = new TTreeReaderArray<Double_t>(myReader, "TrackerEndcapHits.position.y"); 
0062    disks_si_z = new TTreeReaderArray<Double_t>(myReader, "TrackerEndcapHits.position.z");
0063    disks_si_quality = new TTreeReaderArray<Int_t>(myReader, "TrackerEndcapHits.quality");
0064    
0065    // ETOF Hits
0066    endcap_etof_x = new TTreeReaderArray<Double_t>(myReader, "TOFEndcapHits.position.x"); 
0067    endcap_etof_y = new TTreeReaderArray<Double_t>(myReader, "TOFEndcapHits.position.y"); 
0068    endcap_etof_z = new TTreeReaderArray<Double_t>(myReader, "TOFEndcapHits.position.z");
0069    endcap_etof_quality = new TTreeReaderArray<Int_t>(myReader, "TOFEndcapHits.quality"); 
0070     
0071    // Inner MPGD
0072    barrel_mm_x = new TTreeReaderArray<Double_t>(myReader, "MPGDBarrelHits.position.x"); 
0073    barrel_mm_y = new TTreeReaderArray<Double_t>(myReader, "MPGDBarrelHits.position.y"); 
0074    barrel_mm_z = new TTreeReaderArray<Double_t>(myReader, "MPGDBarrelHits.position.z");
0075    barrel_mm_quality = new TTreeReaderArray<Int_t>(myReader, "MPGDBarrelHits.quality");  
0076    
0077    // BarrelTOF
0078    barrel_tof_x = new TTreeReaderArray<Double_t>(myReader, "TOFBarrelHits.position.x"); 
0079    barrel_tof_y = new TTreeReaderArray<Double_t>(myReader, "TOFBarrelHits.position.y"); 
0080    barrel_tof_z = new TTreeReaderArray<Double_t>(myReader, "TOFBarrelHits.position.z");
0081    barrel_tof_quality = new TTreeReaderArray<Int_t>(myReader, "TOFBarrelHits.quality");  
0082    
0083     //Outer MPGD Hits
0084    out_mm_x = new TTreeReaderArray<Double_t>(myReader, "OuterMPGDBarrelHits.position.x"); 
0085    out_mm_y = new TTreeReaderArray<Double_t>(myReader, "OuterMPGDBarrelHits.position.y"); 
0086    out_mm_z = new TTreeReaderArray<Double_t>(myReader, "OuterMPGDBarrelHits.position.z");
0087    out_mm_quality = new TTreeReaderArray<Int_t>(myReader, "OuterMPGDBarrelHits.quality");
0088    
0089     //Forward MPGD 
0090    endcap_fmm_x = new TTreeReaderArray<Double_t>(myReader, "ForwardMPGDEndcapHits.position.x"); 
0091    endcap_fmm_y = new TTreeReaderArray<Double_t>(myReader, "ForwardMPGDEndcapHits.position.y"); 
0092    endcap_fmm_z = new TTreeReaderArray<Double_t>(myReader, "ForwardMPGDEndcapHits.position.z");
0093    endcap_fmm_quality = new TTreeReaderArray<Int_t>(myReader, "ForwardMPGDEndcapHits.quality");    
0094    
0095    //Backward MPGD 
0096    endcap_bmm_x = new TTreeReaderArray<Double_t>(myReader, "BackwardMPGDEndcapHits.position.x"); 
0097    endcap_bmm_y = new TTreeReaderArray<Double_t>(myReader, "BackwardMPGDEndcapHits.position.y"); 
0098    endcap_bmm_z = new TTreeReaderArray<Double_t>(myReader, "BackwardMPGDEndcapHits.position.z");
0099    endcap_bmm_quality = new TTreeReaderArray<Int_t>(myReader, "BackwardMPGDEndcapHits.quality"); 
0100         
0101 
0102    Double_t etamc = 0.; 
0103    int iEvent=-1; 
0104    
0105    TCanvas * c1 = new TCanvas("c1","coutput",1400,1000);
0106    c1->SetMargin(0.10, 0.05 ,0.1,0.08);
0107    c1->SetGridx();
0108    c1->SetGridy();
0109    
0110     TProfile* hits = new TProfile("hits","Nhits (#theta)",70,-3.5,3.5);
0111     hits->GetXaxis()->CenterTitle();
0112     hits->GetYaxis()->CenterTitle();
0113     hits->SetMinimum(0.); 
0114   
0115     std::vector<TVector3> hitPos; // The ePIC tracker
0116     double epsilon = 1.0e-5;
0117     Double_t pmc = 0.;
0118     bool debug = false;
0119     int nhits_SVTIB, nhits_SVTOB, nhits_InMPGD, nhits_BTOF, nhits_OutMPGD, nhits_SVTDisks, nhits_FwdMPGDDisks, nhits_BwdMPGDDisks, nhits_ETOF;
0120     
0121      while (myReader.Next()) 
0122      {
0123       hitPos.clear(); debug = false;
0124       if (iEvent<10) debug = true;
0125       iEvent++;
0126       if (debug) printf("<--------------------Event No. %d---------------------> \n",iEvent);      
0127       // Generated primary track
0128       if (charge.GetSize()>1) continue; // skip event for larger than 1 tracks
0129       Double_t eta_Track = 100.; // set it ouside from -3.5 to 3.5      
0130       for (int j = 0; j < charge.GetSize(); ++j){
0131       
0132       if (status[j] !=1) continue;
0133       pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]);   
0134       Double_t pzmc = pz_mc[j];
0135       Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/pmc))/2));
0136       eta_Track = etamc;
0137       Double_t particle = pdg[j];
0138       }
0139       if (fabs(eta_Track)>3.5) continue; //outside tracker acceptance
0140       // Associated hits with the primary track of momentum (mom)
0141     //  if (fabs(pmc-mom)> epsilon) continue; 
0142     //  if (eta_Track<3.4) continue; // Debug for the hits in a given eta range
0143       if (debug) printf("Eta of the generated track: %f, Momentum (GeV/c): %f \n",eta_Track,pmc);
0144      // ePIC SVT IB Tracker
0145      for (int j = 0; j < vtx_si_x->GetSize(); ++j){
0146        Double_t xhit = vtx_si_x->At(j); Double_t yhit = vtx_si_y->At(j); Double_t zhit = vtx_si_z->At(j);
0147        Int_t quality = vtx_si_quality->At(j);
0148        if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0149      }
0150      nhits_SVTIB =  hitPos.size(); hitPos.clear();
0151      if (debug)  printf("SVT IB Associated hits: %d \n",nhits_SVTIB);
0152      
0153      // ePIC SVT OB Tracker
0154      for (int j = 0; j < barrel_si_x->GetSize(); ++j){
0155        Double_t xhit = barrel_si_x->At(j); Double_t yhit = barrel_si_y->At(j); Double_t zhit = barrel_si_z->At(j);
0156        Int_t quality = barrel_si_quality->At(j); 
0157        if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0158      }
0159      nhits_SVTOB =  hitPos.size(); hitPos.clear();
0160      if (debug)  printf("SVT OB Associated hits: %d \n",nhits_SVTOB);  
0161     
0162     // Inner MPGD Tracker
0163     for (int j = 0; j < barrel_mm_x->GetSize(); ++j){
0164       Double_t xhit = barrel_mm_x->At(j); Double_t yhit = barrel_mm_y->At(j); Double_t zhit = barrel_mm_z->At(j);
0165       Int_t quality = barrel_mm_quality->At(j); 
0166       if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0167     }
0168      nhits_InMPGD =  hitPos.size(); hitPos.clear();
0169      if (debug)  printf("Inner MPGD Associated hits: %d \n",nhits_InMPGD);
0170       
0171     // Outer MPGD Tracker
0172     for (int j = 0; j < out_mm_x->GetSize(); ++j){
0173       Double_t xhit = out_mm_x->At(j); Double_t yhit = out_mm_y->At(j); Double_t zhit = out_mm_z->At(j);
0174       Int_t quality = out_mm_quality->At(j);    
0175       if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); 
0176     }
0177      nhits_OutMPGD =  hitPos.size(); hitPos.clear();
0178      if (debug)  printf("Outer MPGD Associated hits: %d \n",nhits_OutMPGD);
0179          
0180     // BTOF Tracker
0181     for (int j = 0; j < barrel_tof_x->GetSize(); ++j){
0182       Double_t xhit = barrel_tof_x->At(j); Double_t yhit = barrel_tof_y->At(j); Double_t zhit = barrel_tof_z->At(j);
0183       Int_t quality = barrel_tof_quality->At(j); 
0184       if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit)); 
0185     }
0186     
0187      nhits_BTOF =  hitPos.size(); hitPos.clear();
0188      if (debug)  printf("BTOF Associated hits: %d \n",nhits_BTOF);
0189    
0190     // ePIC SVT Disks Hits
0191     for (int j = 0; j < disks_si_x->GetSize(); ++j){
0192       Int_t quality = disks_si_quality->At(j); 
0193       Double_t xhit = disks_si_x->At(j); Double_t yhit = disks_si_y->At(j); Double_t zhit = disks_si_z->At(j);
0194       if (quality==0 && zhit>0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0195       else if (quality==0 && zhit<0) hitPos.push_back(TVector3(xhit,yhit,zhit));      
0196     }
0197     
0198      nhits_SVTDisks =  hitPos.size(); hitPos.clear();
0199      if (debug)  printf("SVT Disks Associated hits: %d \n",nhits_SVTDisks);
0200     
0201      // ETOF Tracker (Forward)
0202     for (int j = 0; j < endcap_etof_x->GetSize(); ++j){
0203     Double_t xhit = endcap_etof_x->At(j); Double_t yhit = endcap_etof_y->At(j); Double_t zhit = endcap_etof_z->At(j);
0204     Int_t quality = endcap_etof_quality->At(j);    
0205     if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0206     }  
0207     
0208      nhits_ETOF =  hitPos.size(); hitPos.clear();
0209      if (debug)  printf("ETOF Associated hits: %d \n",nhits_ETOF);
0210      
0211      // Forward MPGD
0212     for (int j = 0; j < endcap_fmm_x->GetSize(); ++j){
0213     Double_t xhit = endcap_fmm_x->At(j); Double_t yhit = endcap_fmm_y->At(j); Double_t zhit = endcap_fmm_z->At(j);
0214     Int_t quality = endcap_fmm_quality->At(j);    
0215     if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0216     } 
0217      nhits_FwdMPGDDisks =  hitPos.size(); hitPos.clear();
0218      if (debug)  printf("Forward MPGD Associated hits: %d \n",nhits_FwdMPGDDisks);
0219     
0220     // Backward MPGD
0221     for (int j = 0; j < endcap_bmm_x->GetSize(); ++j){
0222     Double_t xhit = endcap_bmm_x->At(j); Double_t yhit = endcap_bmm_y->At(j); Double_t zhit = endcap_bmm_z->At(j);
0223     Int_t quality = endcap_bmm_quality->At(j);    
0224     if (quality==0) hitPos.push_back(TVector3(xhit,yhit,zhit));
0225     } 
0226      
0227      nhits_BwdMPGDDisks =  hitPos.size(); hitPos.clear();
0228      if (debug)  printf("Backward MPGD Associated hits: %d \n",nhits_BwdMPGDDisks);
0229    
0230     int nhits = nhits_SVTIB + nhits_SVTOB + nhits_InMPGD + nhits_BTOF + nhits_OutMPGD + nhits_SVTDisks + nhits_FwdMPGDDisks + nhits_BwdMPGDDisks + nhits_ETOF;
0231     if (nhits>0) hits->Fill(eta_Track,nhits); 
0232 
0233  }
0234  
0235   c1->cd();
0236   gPad->SetTicks(1,1);
0237   hits->SetTitle(";#eta_{mc};Nhits");
0238   hits->SetLineWidth(2);
0239   hits->Draw("hist");
0240   TPaveText *pt = new TPaveText(0.1, 0.95, 0.9, 1.0, "NDC");
0241   pt->AddText(Form("p_{mc} = %1.1f GeV/c",pmc));
0242   pt->SetBorderSize(0);     
0243   pt->SetFillStyle(0);      
0244   pt->SetTextAlign(23); 
0245   pt->Draw();
0246 
0247   c1->SaveAs(Form("%s/Nhits_vs_eta.png", output_prefix.Data()));
0248   c1->SaveAs(Form("%s/Nhits_vs_eta.root", output_prefix.Data()));
0249 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!