Tracking/BackgroundStudy/TrackingBG.C

0001 void TrackingBG(TString type="Forced"){
0002
0003     //Define Style

0004     gStyle->SetOptStat(0);
0005     gStyle->SetPadBorderMode(0);
0006     gStyle->SetFrameBorderMode(0);
0007     gStyle->SetFrameLineWidth(2);
0008     gStyle->SetLabelSize(0.035,"X");
0009     gStyle->SetLabelSize(0.035,"Y");
0010     //gStyle->SetLabelOffset(0.01,"X");

0011     //gStyle->SetLabelOffset(0.01,"Y");

0012     gStyle->SetTitleXSize(0.04);
0013     gStyle->SetTitleXOffset(0.9);
0014     gStyle->SetTitleYSize(0.04);
0015     gStyle->SetTitleYOffset(0.9);
0016
0017     TFile *fout = new TFile("TrackingBG_10x275_Forced.root","RECREATE");
0018
0019     //Define histograms

0020     fout->cd();
0021     //--

0022     //Eta distributions for generated particles

0023     TH1 * hMC1 = new TH1D("hMC1","",100,-15,10); //DIS

0024     hMC1->SetLineColor(kBlack); hMC1->SetLineWidth(2);
0025     hMC1->GetXaxis()->SetTitle("Particle #eta");hMC1->GetXaxis()->CenterTitle();
0026
0027     TH1 * hMC2 = new TH1D("hMC2","",100,-15,10); //SR

0028     hMC2->SetLineColor(kBlue); hMC2->SetLineWidth(2);
0029     hMC2->GetXaxis()->SetTitle("Particle #eta");hMC2->GetXaxis()->CenterTitle();
0030
0031     TH1 * hMC3 = new TH1D("hMC3","",100,-15,10); //Bremstrahlung

0032     hMC3->SetLineColor(kOrange); hMC3->SetLineWidth(2);
0033     hMC3->GetXaxis()->SetTitle("Particle #eta");hMC3->GetXaxis()->CenterTitle();
0034
0035     TH1 * hMC4 = new TH1D("hMC4","",100,-15,10); //Coulomb

0036     hMC4->SetLineColor(kGreen); hMC4->SetLineWidth(2);
0037     hMC4->GetXaxis()->SetTitle("Particle #eta");hMC4->GetXaxis()->CenterTitle();
0038
0039     TH1 * hMC5 = new TH1D("hMC5","",100,-15,10); //Touschek

0040     hMC5->SetLineColor(kRed); hMC5->SetLineWidth(2);
0041     hMC5->GetXaxis()->SetTitle("Particle #eta");hMC5->GetXaxis()->CenterTitle();
0042
0043     TH1 * hMC6 = new TH1D("hMC6","",100,-15,10); //Proton beam gas

0044     hMC6->SetLineColor(kMagenta); hMC6->SetLineWidth(2);
0045     hMC6->GetXaxis()->SetTitle("Particle #eta");hMC6->GetXaxis()->CenterTitle();
0046     //--

0047     //Digitized hit rate on SVT Disks

0048     const int N_disks = 10;
0049     double disk_min[N_disks] = {-1040.0, -860.0, -660.0, -460.0, -260.0, 240.0, 440.0, 690.0, 990.0 , 1340.0};
0050     double disk_max[N_disks] = {-1000.0, -840.0, -640.0, -440.0, -240.0, 260.0, 460.0, 710.0, 1010.0, 1360.0};
0051     TString disk_name[N_disks] = {"E-Si Disk 4","E-Si Disk 3","E-Si Disk 2","E-Si Disk 1","E-Si Disk 0",
0052                                  "H-Si Disk 0","H-Si Disk 1","H-Si Disk 2","H-Si Disk 3","H-Si Disk 4"};
0053
0054     TH2 *hSVTDiskRate[N_disks];
0055
0056     //Loop over momentum points

0057     for(int idisk = 0; idisk < N_disks; idisk++){
0058
0059         TString hname = TString::Format("hSVTDiskRate_%d", idisk);
0060         TString htitle = TString::Format("Digitized hit Rate per RSU per 1 ms: %s",disk_name[idisk].Data());
0061
0062         hSVTDiskRate[idisk] = new TH2D(hname,htitle,50,-500,500,50,-500,500);
0063         hSVTDiskRate[idisk]->GetXaxis()->SetTitle("X [mm]");hSVTDiskRate[idisk]->GetXaxis()->CenterTitle();
0064         hSVTDiskRate[idisk]->GetYaxis()->SetTitle("Y [mm]");hSVTDiskRate[idisk]->GetYaxis()->CenterTitle();
0065
0066     }
0067     //--

0068     //Digitized hit rate on SVT Vertexing layers

0069     const int N_vtx = 3;
0070     double vtx_min[N_vtx] = {30.,46.,115.};
0071     double vtx_max[N_vtx] = {42.,60.,130.};
0072
0073     TString vtx_name[N_vtx] = {"SVT L0","SVT L1","SVT L2"};
0074
0075     TH2 *hSVTVtxRate[N_vtx];
0076
0077     //Loop over momentum points

0078     for(int ivtx = 0; ivtx < N_vtx; ivtx++){
0079
0080         TString hname = TString::Format("hSVTVtxRate_%d", ivtx);
0081         TString htitle = TString::Format("Digitized hit Rate per RSU per 1 ms: %s",vtx_name[ivtx].Data());
0082
0083         hSVTVtxRate[ivtx] = new TH2D(hname,htitle,20,-200,200,40,-400,400);
0084         hSVTVtxRate[ivtx]->GetXaxis()->SetTitle("Z [mm]");hSVTVtxRate[ivtx]->GetXaxis()->CenterTitle();
0085         hSVTVtxRate[ivtx]->GetYaxis()->SetTitle("r #phi [mm]");hSVTVtxRate[ivtx]->GetYaxis()->CenterTitle();
0086
0087     }
0088     //--

0089     //Digitized hit rate on SVT Barrel layers

0090     const int N_barrel = 2;
0091     double barrel_min[N_barrel] = {260.,410.};
0092     double barrel_max[N_barrel] = {280.,450.};
0093
0094     TString barrel_name[N_barrel] = {"SVT L3","SVT L4"};
0095
0096     TH2 *hSVTbarrelRate[N_barrel];
0097
0098     //Loop over momentum points

0099     for(int ibarrel = 0; ibarrel < N_barrel; ibarrel++){
0100
0101         TString hname = TString::Format("hSVTbarrelRate_%d", ibarrel);
0102         TString htitle = TString::Format("Digitized hit Rate per RSU per 1 ms: %s", barrel_name[ibarrel].Data());
0103
0104         hSVTbarrelRate[ibarrel] = new TH2D(hname,htitle,40,-400,400,150,-1500,1500);
0105         hSVTbarrelRate[ibarrel]->GetXaxis()->SetTitle("Z [mm]");hSVTbarrelRate[ibarrel]->GetXaxis()->CenterTitle();
0106         hSVTbarrelRate[ibarrel]->GetYaxis()->SetTitle("r #phi [mm]");hSVTbarrelRate[ibarrel]->GetYaxis()->CenterTitle();
0107
0108     }
0109     //--

0110
0111     //Open file

0112     TString input = "10x275_" + type + "/eicrecon_*.root";
0113
0114     TChain *tree = new TChain("events");
0115     tree->Add(input.Data()); //Events with backgrounds present

0116
0117     cout<<"Running analysis on "<<input<<"!"<<endl;
0118     cout<<"Analyzing "<<tree->GetEntries()<<" events!"<<endl;
0119
0120     //Create Array Reader

0121     TTreeReader tr(tree);
0122
0123     TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0124     TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0125     TTreeReaderArray<double> gen_px(tr, "MCParticles.momentum.x");
0126     TTreeReaderArray<double> gen_py(tr, "MCParticles.momentum.y");
0127     TTreeReaderArray<double> gen_pz(tr, "MCParticles.momentum.z");
0128     TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0129     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0130     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0131     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0132     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0133     TTreeReaderArray<float> SiEC_RecHit_posx(tr,"SiEndcapTrackerRecHits.position.x");
0134     TTreeReaderArray<float> SiEC_RecHit_posy(tr,"SiEndcapTrackerRecHits.position.y");
0135     TTreeReaderArray<float> SiEC_RecHit_posz(tr,"SiEndcapTrackerRecHits.position.z");
0136     TTreeReaderArray<float> SiVtx_RecHit_posx(tr,"SiBarrelVertexRecHits.position.x");
0137     TTreeReaderArray<float> SiVtx_RecHit_posy(tr,"SiBarrelVertexRecHits.position.y");
0138     TTreeReaderArray<float> SiVtx_RecHit_posz(tr,"SiBarrelVertexRecHits.position.z");
0139     TTreeReaderArray<float> Sibarrel_RecHit_posx(tr,"SiBarrelTrackerRecHits.position.x");
0140     TTreeReaderArray<float> Sibarrel_RecHit_posy(tr,"SiBarrelTrackerRecHits.position.y");
0141     TTreeReaderArray<float> Sibarrel_RecHit_posz(tr,"SiBarrelTrackerRecHits.position.z");
0142
0143     //Define other variables

0144     int counter(0);
0145     TLorentzVector gen_vec;
0146
0147     //Loop over events

0148     while (tr.Next()) {
0149
0150         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0151         counter++;
0152
0153         //Loop over generated particles

0154         for(int igen=0;igen<gen_status.GetSize();igen++){
0155
0156             auto status = gen_status[igen];
0157             gen_vec.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0158             auto eta = gen_vec.Eta();
0159
0160             if( status==6001 )
0161                 hMC6->Fill(eta);
0162             else if( status==5001 )
0163                 hMC5->Fill(eta);
0164             else if( status==4001 )
0165                 hMC4->Fill(eta);
0166             else if( status==3001 )
0167                 hMC3->Fill(eta);
0168             else if( status==2001 )
0169                 hMC2->Fill(eta);
0170             else if( status==1 )
0171                 hMC1->Fill(eta);
0172
0173         } //End loop over generated particles

0174
0175         //Loop over Digitized hits for SVT disks

0176         for(int ihit=0;ihit < SiEC_RecHit_posx.GetSize();ihit++){
0177
0178             for(int idisk = 0; idisk < N_disks; idisk++){
0179                 if( SiEC_RecHit_posz[ihit] > disk_min[idisk] && SiEC_RecHit_posz[ihit] < disk_max[idisk] )
0180                     hSVTDiskRate[idisk]->Fill(SiEC_RecHit_posx[ihit],SiEC_RecHit_posy[ihit]);
0181             }
0182
0183         } // End Loop over Digitized hits for SVT disks

0184
0185         //Loop over Digitized hits for SVT Vtx layers

0186         for(int ihit=0;ihit < SiVtx_RecHit_posx.GetSize();ihit++){
0187
0188             for(int ivtx = 0; ivtx < N_vtx; ivtx++){
0189
0190                 auto hit_r = std::hypot(SiVtx_RecHit_posx[ihit],SiVtx_RecHit_posy[ihit]);
0191                 auto hit_phi = std::atan2(SiVtx_RecHit_posy[ihit],SiVtx_RecHit_posx[ihit]);
0192                 auto hit_z = SiVtx_RecHit_posz[ihit];
0193
0194                 if( hit_r > vtx_min[ivtx] && hit_r < vtx_max[ivtx] )
0195                     hSVTVtxRate[ivtx]->Fill(hit_z,hit_r*hit_phi);
0196             }
0197
0198         } // End Loop over Digitized hits for SVT Vtx layers

0199
0200         //Loop over Digitized hits for SVT Barrel layers

0201         for(int ihit=0;ihit < Sibarrel_RecHit_posx.GetSize();ihit++){
0202
0203             for(int ibarrel = 0; ibarrel < N_barrel; ibarrel++){
0204
0205                 auto hit_r = std::hypot(Sibarrel_RecHit_posx[ihit],Sibarrel_RecHit_posy[ihit]);
0206                 auto hit_phi = std::atan2(Sibarrel_RecHit_posy[ihit],Sibarrel_RecHit_posx[ihit]);
0207                 auto hit_z = Sibarrel_RecHit_posz[ihit];
0208
0209                 if( hit_r > barrel_min[ibarrel] && hit_r < barrel_max[ibarrel] )
0210                     hSVTbarrelRate[ibarrel]->Fill(hit_z,hit_r*hit_phi);
0211             }
0212
0213         } // End Loop over Digitized hits for SVT Barrel layers

0214
0215     } //End loop over events

0216
0217     //Check SR rate

0218     cout <<"------"<<endl;
0219     cout <<"Total number of SR photons generated: "<< hMC2->GetEntries() << endl;
0220     cout <<"SR photon rate [MHz]: " << 1.0e-6 * ( hMC2->GetEntries() / (2.0e-6 * tree->GetEntries()) ) << endl;
0221     cout <<"------"<<endl;
0222
0223     //Scale hit rate histograms to 1 ms.

0224     //Each 'event' is 2us.

0225     auto SFac = 1.0e-3 / (2.0e-6 * tree->GetEntries());
0226
0227     for(int idisk = 0; idisk < N_disks; idisk++){
0228         hSVTDiskRate[idisk]->Scale(SFac);
0229     }
0230
0231     for(int ivtx = 0; ivtx < N_vtx; ivtx++){
0232         hSVTVtxRate[ivtx]->Scale(SFac);
0233     }
0234
0235     for(int ibarrel = 0; ibarrel < N_barrel; ibarrel++){
0236         hSVTbarrelRate[ibarrel]->Scale(SFac);
0237     }
0238
0239     //Write rate histograms to output ROOT file

0240     fout->cd();
0241     for(int idisk = 0; idisk < N_disks; idisk++){
0242         hSVTDiskRate[idisk]->Write();
0243     }
0244
0245     for(int ivtx = 0; ivtx < N_vtx; ivtx++){
0246         hSVTVtxRate[ivtx]->Write();
0247     }
0248
0249     for(int ibarrel = 0; ibarrel < N_barrel; ibarrel++){
0250         hSVTbarrelRate[ibarrel]->Write();
0251     }
0252
0253     // Make plots

0254     TCanvas *c1 = new TCanvas("c1");
0255
0256     auto frame1 = c1->DrawFrame(-15,0.1,10,1e8);
0257     frame1->SetTitle("10x275 GeV: Forced DIS Configuration");
0258     frame1->GetXaxis()->SetTitle("Particle #eta");frame1->GetXaxis()->CenterTitle();
0259     gPad->SetLogy();
0260
0261     hMC1->Draw("same");
0262     hMC2->Draw("same");
0263     hMC3->Draw("same");
0264     hMC4->Draw("same");
0265     hMC5->Draw("same");
0266     hMC6->Draw("same");
0267
0268     TLegend *leg1 = new TLegend(0.175, 0.65, 0.375, 0.85, "MCParticles.generatorStatus");
0269     leg1->AddEntry(hMC1, "1: DIS", "l");
0270     leg1->AddEntry(hMC2, "2001: SR", "l");
0271     leg1->AddEntry(hMC3, "3001: Bremstrahlung", "l");
0272     leg1->AddEntry(hMC4, "4001: Coulomb", "l");
0273     leg1->AddEntry(hMC5, "5001: Touschek", "l");
0274     leg1->AddEntry(hMC6, "6001: Proton beam gas", "l");
0275     leg1->Draw();
0276
0277     c1->Print("plots/TrackingBG_10x275_Forced.pdf[");
0278     c1->Print("plots/TrackingBG_10x275_Forced.pdf");
0279
0280     TCanvas *c2[N_disks];
0281     for(int idisk = 0; idisk < N_disks; idisk++){
0282         c2[idisk] = new TCanvas(Form("c2[%d]",idisk));
0283         hSVTDiskRate[idisk]->Draw("colz");
0284         c2[idisk]->Print("plots/TrackingBG_10x275_Forced.pdf");
0285
0286     }
0287
0288     TCanvas *c3[N_vtx];
0289     for(int ivtx = 0; ivtx < N_vtx; ivtx++){
0290         c3[ivtx] = new TCanvas(Form("c3[%d]",ivtx));
0291         hSVTVtxRate[ivtx]->Draw("colz");
0292         c3[ivtx]->Print("plots/TrackingBG_10x275_Forced.pdf");
0293
0294     }
0295
0296     TCanvas *c4[N_barrel];
0297     for(int ibarrel = 0; ibarrel < N_barrel; ibarrel++){
0298         c4[ibarrel] = new TCanvas(Form("c4[%d]",ibarrel));
0299         hSVTbarrelRate[ibarrel]->Draw("colz");
0300         c4[ibarrel]->Print("plots/TrackingBG_10x275_Forced.pdf");
0301
0302     }
0303
0304     c4[N_barrel-1]->Print("plots/TrackingBG_10x275_Forced.pdf]");
0305
0306     // Close output ROOT file

0307     fout->Close();
0308
0309 }