Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-01 08:59:06

0001 //------------------

0002 // energy_set 0 is 18x275 GeV

0003 // energy_set 1 is 10x100 GeV

0004 void DIS_reconstruction(int energy_set = 0, int Q2_set = 1, bool use_campaign = 0){
0005 
0006     //Define Style

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

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

0015     gStyle->SetTitleXSize(0.04);
0016     gStyle->SetTitleXOffset(0.9);
0017     gStyle->SetTitleYSize(0.04);
0018     gStyle->SetTitleYOffset(0.9);
0019 
0020     //Define histograms

0021  
0022     //---Eta distributions---

0023     TH1 *h1a = new TH1D("h1a","Generated Charged Particles",100,-4,4);
0024     h1a->GetXaxis()->SetTitle("#eta_{gen.}");h1a->GetXaxis()->CenterTitle();
0025     h1a->SetLineColor(kTeal);h1a->SetLineWidth(2);
0026 
0027     TH1 *h1a1 = new TH1D("h1a1","Generated Charged Particles",100,-4,4);  //Minimum momentum cut of Pt > 200 MeV/c

0028     h1a1->GetXaxis()->SetTitle("#eta_{gen.}");h1a1->GetXaxis()->CenterTitle();
0029     h1a1->SetLineColor(kRed);h1a1->SetLineWidth(2);
0030 
0031     TH1 *h1a2 = new TH1D("h1a2","Generated Charged Particles",100,-4,4);  //Minimum momentum cut of Pt > 500 MeV/c

0032     h1a2->GetXaxis()->SetTitle("#eta_{gen.}");h1a2->GetXaxis()->CenterTitle();
0033     h1a2->SetLineColor(kBlack);h1a2->SetLineWidth(2);
0034     h1a2->SetFillColor(kBlack);h1a2->SetFillStyle(3244);
0035 
0036     TH1 *h1b = new TH1D("h1b","Reconstructed Real-seeded tracks",100,-4,4);
0037     h1b->GetXaxis()->SetTitle("#eta_{rec.}");h1b->GetXaxis()->CenterTitle();
0038     h1b->SetLineColor(kGreen);h1b->SetLineWidth(2);
0039 
0040     TH1 *h1b1 = new TH1D("h1b1","Reconstructed Real-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 200 MeV/c

0041     h1b1->GetXaxis()->SetTitle("#eta_{rec.}");h1b1->GetXaxis()->CenterTitle();
0042     h1b1->SetLineColor(kBlue);h1b1->SetLineWidth(2);
0043 
0044     TH1 *h1b2 = new TH1D("h1b2","Reconstructed Real-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 500 MeV/c

0045     h1b2->GetXaxis()->SetTitle("#eta_{rec.}");h1b2->GetXaxis()->CenterTitle();
0046     h1b2->SetLineColor(kRed);h1b2->SetLineWidth(2);
0047     h1b2->SetMarkerColor(kRed);h1b2->SetMarkerStyle(kFullCrossX);//h1b2->SetMarkerSize(0.75);

0048 
0049     TH1 *h1c = new TH1D("h1c","Reconstructed Truth-seeded tracks",100,-4,4);
0050     h1c->GetXaxis()->SetTitle("#eta_{rec.}");h1c->GetXaxis()->CenterTitle();
0051     h1c->SetLineColor(kRed);h1c->SetLineWidth(2);
0052 
0053     TH1 *h1c1 = new TH1D("h1c1","Reconstructed Truth-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 200 MeV/c

0054     h1c1->GetXaxis()->SetTitle("#eta_{rec.}");h1c1->GetXaxis()->CenterTitle();
0055     h1c1->SetLineColor(kOrange);h1c1->SetLineWidth(2);
0056 
0057     TH1 *h1c2 = new TH1D("h1c2","Reconstructed Truth-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 500 MeV/c

0058     h1c2->GetXaxis()->SetTitle("#eta_{rec.}");h1c2->GetXaxis()->CenterTitle();
0059     h1c2->SetLineColor(kMagenta);h1c2->SetLineWidth(2);
0060     h1c2->SetMarkerColor(kMagenta);h1c2->SetMarkerStyle(kFullCrossX);//h1c2->SetMarkerSize(0.75);

0061 
0062     TH1 *h1rb1 = new TH1D("h1rb1","",100,-4,4); //Real-seeded tracks (Pt > 200 MeV/c cut)

0063     TH1 *h1rc1 = new TH1D("h1rc1","",100,-4,4); //Truth-seeded tracks (Pt > 200 MeV/c cut)

0064     TH1 *h1rb2 = new TH1D("h1rb2","",100,-4,4); //Real-seeded tracks (Pt > 500 MeV/c cut)

0065     TH1 *h1rc2 = new TH1D("h1rc2","",100,-4,4); //Truth-seeded tracks (Pt > 500 MeV/c cut)

0066 
0067     //Transverse Momentum distributions

0068     //For generated charged particles and tracks, place |eta| < 3.5 cut

0069     TH1 *h2a = new TH1D("h2a","Generated Charged Particles",100,0,15);
0070     h2a->GetXaxis()->SetTitle("P_{T} [GeV/c]");h2a->GetXaxis()->CenterTitle();
0071     h2a->SetLineColor(kBlack);h2a->SetLineWidth(2);
0072     h2a->SetFillColor(kBlack);h2a->SetFillStyle(3244);
0073 
0074     TH1 *h2b = new TH1D("h2b","Reconstructed Real-seeded tracks",100,0,15);
0075     h2b->GetXaxis()->SetTitle("P_{T} [GeV/c]");h2b->GetXaxis()->CenterTitle();
0076     h2b->SetLineColor(kRed);h2b->SetLineWidth(2);
0077     h2b->SetMarkerColor(kRed);h2b->SetMarkerStyle(kFullCrossX);
0078 
0079     TH1 *h2c = new TH1D("h2c","Reconstructed Truth-seeded tracks",100,0,15);
0080     h2c->GetXaxis()->SetTitle("P_{T} [GeV/c]");h2c->GetXaxis()->CenterTitle();
0081     h2c->SetLineColor(kMagenta);h2c->SetLineWidth(2);
0082     h2c->SetMarkerColor(kMagenta);h2c->SetMarkerStyle(kFullCrossX);
0083 
0084     //Transverse Momentum distributions for generated charged particles

0085     //with matching track. Only consider particles w/ |eta| < 3.5

0086     TH1 *h3a = new TH1D("h3a","",100,0,15);
0087     h3a->GetXaxis()->SetTitle("P_{T} [GeV/c]");h3a->GetXaxis()->CenterTitle();
0088     h3a->SetLineColor(kRed);h3a->SetLineWidth(2);
0089     h3a->SetMarkerColor(kRed);h3a->SetMarkerStyle(kFullCircle);    
0090 
0091     TH1 *h3b = new TH1D("h3b","",100,0,15);
0092     h3b->GetXaxis()->SetTitle("P_{T} [GeV/c]");h3b->GetXaxis()->CenterTitle();
0093     h3b->SetLineColor(kMagenta);h3b->SetLineWidth(2);
0094     h3b->SetMarkerColor(kMagenta);h3b->SetMarkerStyle(kOpenCircle);
0095 
0096     TH1 *h3ra = new TH1D("h3ra","",100,0,15);
0097     TH1 *h3rb = new TH1D("h3rb","",100,0,15);
0098 
0099     //File information

0100     TString path;
0101     if(use_campaign) path = "./campaign_input/";
0102     else path = "./input/";
0103 
0104     TString beam_energies;
0105     if(energy_set == 0) beam_energies = "18x275";
0106     else if(energy_set == 1) beam_energies = "10x100";
0107 
0108     TString run_name;
0109     if(use_campaign) run_name.Form("list_%s_Q2_%d.txt",beam_energies.Data(),Q2_set);
0110     else run_name.Form("%s/Q2_%d/eicrecon_*root",beam_energies.Data(),Q2_set);
0111     
0112     //Open File

0113     TString input = path + run_name;
0114     TChain *tree = new TChain("events");
0115 
0116     if(use_campaign){
0117     std::ifstream in(input);
0118     std::string file("");
0119     while (in >> file)
0120         tree->Add(file.data());
0121         in.close();
0122     }
0123     else{
0124         tree->Add(input.Data());
0125    }
0126 
0127     cout<<"Running file "<<input<<"!"<<endl;
0128     cout<<"Analyzing "<<tree->GetEntries()<<" events!"<<endl;
0129 
0130     //Create Array Reader

0131     TTreeReader tr(tree);
0132 
0133     TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0134     TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0135     TTreeReaderArray<double> gen_px(tr, "MCParticles.momentum.x");
0136     TTreeReaderArray<double> gen_py(tr, "MCParticles.momentum.y");
0137     TTreeReaderArray<double> gen_pz(tr, "MCParticles.momentum.z");
0138     TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass"); //Not important here

0139     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0140     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0141     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0142     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0143 
0144     TTreeReaderArray<float> rec_px(tr, "ReconstructedChargedParticles.momentum.x");
0145     TTreeReaderArray<float> rec_py(tr, "ReconstructedChargedParticles.momentum.y");
0146     TTreeReaderArray<float> rec_pz(tr, "ReconstructedChargedParticles.momentum.z");
0147     TTreeReaderArray<float> rec_mass(tr, "ReconstructedChargedParticles.mass");
0148     TTreeReaderArray<int> rec_type(tr, "ReconstructedChargedParticles.type"); //Type 0: Match to generated particle

0149     TTreeReaderArray<int> rec_pdg(tr, "ReconstructedChargedParticles.PDG"); //Uses PID lookup tables

0150 
0151     TTreeReaderArray<float> rec_ts_px(tr, "ReconstructedTruthSeededChargedParticles.momentum.x");
0152     TTreeReaderArray<float> rec_ts_py(tr, "ReconstructedTruthSeededChargedParticles.momentum.y");
0153     TTreeReaderArray<float> rec_ts_pz(tr, "ReconstructedTruthSeededChargedParticles.momentum.z");
0154     TTreeReaderArray<float> rec_ts_mass(tr, "ReconstructedTruthSeededChargedParticles.mass");
0155     TTreeReaderArray<int> rec_ts_type(tr, "ReconstructedTruthSeededChargedParticles.type"); //Match to generated particle

0156     TTreeReaderArray<int> rec_ts_pdg(tr, "ReconstructedTruthSeededChargedParticles.PDG"); //Uses PID lookup tables

0157 
0158     TTreeReaderArray<unsigned int> assoc_rec_id(tr, "ReconstructedChargedParticleAssociations.recID");
0159     TTreeReaderArray<unsigned int> assoc_sim_id(tr, "ReconstructedChargedParticleAssociations.simID");
0160     TTreeReaderArray<float> assoc_weight(tr, "ReconstructedChargedParticleAssociations.weight");
0161 
0162     TTreeReaderArray<unsigned int> assoc_ts_rec_id(tr, "ReconstructedTruthSeededChargedParticleAssociations.recID");
0163     TTreeReaderArray<unsigned int> assoc_ts_sim_id(tr, "ReconstructedTruthSeededChargedParticleAssociations.simID");
0164     TTreeReaderArray<float> assoc_ts_weight(tr, "ReconstructedTruthSeededChargedParticleAssociations.weight");
0165 
0166     //Other variables

0167     TLorentzVector gen_vec;
0168     TVector3 gen_vertex;
0169 
0170     TLorentzVector rec_vec;
0171     TVector3 track_vec; //Reconstructed track momentum vector

0172     int counter(0);
0173 
0174     //Loop over events

0175     while (tr.Next()) {
0176 
0177         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0178         counter++;
0179 
0180         //Loop over generated particles

0181         for(int igen=0;igen<gen_status.GetSize();igen++){
0182             
0183             auto charge = gen_charge[igen];
0184             auto status = gen_status[igen];
0185 
0186             //Require final-state, charged particle (no secondaries)

0187             if( status==1 && fabs(charge) > 0.01 ){
0188                 gen_vec.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0189                 gen_vertex.SetXYZ(gen_vx[igen],gen_vy[igen],gen_vz[igen]);
0190                 
0191         auto eta = gen_vec.Eta();
0192         auto pt = gen_vec.Pt();
0193 
0194                 //Fill eta histograms

0195                 h1a->Fill(eta);
0196                 if( pt > 0.2 ) h1a1->Fill(eta);
0197                 if( pt > 0.5 ) h1a2->Fill(eta);
0198 
0199         //Fill Pt histogram

0200         if( fabs(eta) < 3.5 ) h2a->Fill(pt);        
0201 
0202         //For |eta| < 3.5, find if there is an associated reconstructed track with weight > 0.8

0203         
0204         //Real-seeded tracks

0205         for(int iassoc = 0; iassoc < assoc_weight.GetSize(); iassoc++){
0206 
0207             if(fabs(eta) < 3.5 && assoc_sim_id[iassoc] == igen && assoc_weight[iassoc]>0.8){
0208                 //For a given track,, only one MC Particle can have a weight > 0.8

0209                 //---

0210                 //It is possible that a given MC Particle is associated to multiple

0211                 //tracks with weight > 0.8, although this is unlikely due to shared hit

0212                 //cut in ambiguity solver

0213                 h3a->Fill(pt);  
0214             }
0215         }
0216 
0217         //Truth-seeded tracks

0218         for(int iassoc = 0; iassoc < assoc_ts_weight.GetSize(); iassoc++){
0219 
0220                         if(fabs(eta) < 3.5 && assoc_ts_sim_id[iassoc] == igen && assoc_ts_weight[iassoc]>0.8){
0221                                 //For a given track,, only one MC Particle can have a weight > 0.8

0222                                 //---

0223                                 //It is possible that a given MC Particle is associated to multiple

0224                                 //tracks with weight > 0.8, although this is unlikely due to shared hit

0225                                 //cut in ambiguity solver

0226                                 h3b->Fill(pt);
0227                         }
0228                 }
0229 
0230             } //if generated final-state charged particle

0231         } //End loop over generated particles

0232         
0233         //Loop over reconstructed real-seeded charged particles (copy of tracks with PID info)

0234         int rec_mult = rec_type.GetSize();
0235 
0236         for(int irec=0;irec<rec_mult;irec++){
0237 
0238             rec_vec.SetXYZM(rec_px[irec],rec_py[irec],rec_pz[irec],rec_mass[irec]);
0239 
0240         auto eta = rec_vec.Eta();
0241             auto pt = rec_vec.Pt();            
0242 
0243             //Fill eta histograms

0244             h1b->Fill(eta);
0245             if( pt > 0.2 ) h1b1->Fill(eta);
0246             if( pt > 0.5 ) h1b2->Fill(eta);
0247             //if( rec_vec.Pt() > 0.2 && rec_type[irec] == 0 ) h1b2->Fill(rec_vec.Eta());

0248         
0249         //Fill pt histogram

0250         if( fabs(eta) < 3.5 ) h2b->Fill(pt);
0251 
0252         } //End loop over reconstructed particles

0253 
0254         //Loop over reconstructed truth-seeded charged particles (copy of tracks with PID info)

0255         int rec_ts_mult = rec_ts_type.GetSize();
0256 
0257         for(int irec=0;irec<rec_ts_mult;irec++){
0258 
0259             rec_vec.SetXYZM(rec_ts_px[irec],rec_ts_py[irec],rec_ts_pz[irec],rec_ts_mass[irec]);
0260             
0261         auto eta = rec_vec.Eta();
0262             auto pt = rec_vec.Pt();
0263     
0264             //Fill eta histograms

0265             h1c->Fill(rec_vec.Eta());
0266             if( rec_vec.Pt() > 0.2 ) h1c1->Fill(rec_vec.Eta());
0267             if( rec_vec.Pt() > 0.5 ) h1c2->Fill(rec_vec.Eta());
0268             //if( rec_vec.Pt() > 0.2 && rec_ts_type[irec] == 0 ) h1c2->Fill(rec_vec.Eta());

0269 
0270         //Fill pt histogram

0271         if( fabs(eta) < 3.5 ) h2c->Fill(pt);
0272 
0273         } //End loop over reconstructed particles

0274 
0275     } //End loop over events

0276 
0277     //Make ratio histograms

0278     h1rb1 = (TH1*) h1b1->Clone("h1rb1");
0279     h1rb1->Divide(h1a1);
0280     h1rb1->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 200 MeV/c");
0281     h1rb1->GetXaxis()->SetTitle("#eta");h1rb1->GetXaxis()->CenterTitle();
0282     h1rb1->GetYaxis()->SetTitle("Ratio");h1rb1->GetYaxis()->CenterTitle();
0283 
0284     h1rc1 = (TH1*) h1c1->Clone("h1rc1");
0285     h1rc1->Divide(h1a1);
0286     h1rc1->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 200 MeV/c");
0287     h1rc1->GetXaxis()->SetTitle("#eta");h1rc1->GetXaxis()->CenterTitle();
0288     h1rc1->GetYaxis()->SetTitle("Ratio");h1rc1->GetYaxis()->CenterTitle();
0289 
0290     h1rb2 = (TH1*) h1b2->Clone("h1rb2");
0291     h1rb2->Divide(h1a2);
0292     h1rb2->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 500 MeV/c");
0293     h1rb2->GetXaxis()->SetTitle("#eta");h1rb2->GetXaxis()->CenterTitle();
0294     h1rb2->GetYaxis()->SetTitle("Ratio");h1rb2->GetYaxis()->CenterTitle();
0295 
0296     h1rc2 = (TH1*) h1c2->Clone("h1rc2");
0297     h1rc2->Divide(h1a2);
0298     h1rc2->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 500 MeV/c");
0299     h1rc2->GetXaxis()->SetTitle("#eta");h1rc2->GetXaxis()->CenterTitle();
0300     h1rc2->GetYaxis()->SetTitle("Ratio");h1rc2->GetYaxis()->CenterTitle();
0301 
0302     h3ra = (TH1*) h3a->Clone("h3ra");
0303     h3ra->Divide(h2a);
0304 
0305     h3rb = (TH1*) h3b->Clone("h3rb");
0306     h3rb->Divide(h2a);
0307 
0308     //Make plots

0309     //Generated charged particles

0310     TCanvas *c1a = new TCanvas("c1a");
0311     h1a->Draw();
0312     h1a1->Draw("same");
0313     h1a2->Draw("same");
0314 
0315     TLegend *leg1a = new TLegend(0.125,0.7,0.375,0.875);
0316     leg1a->SetBorderSize(0);leg1a->SetFillStyle(0);
0317     leg1a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0318     leg1a->AddEntry("h1a","All generated charged particles","l");
0319     leg1a->AddEntry("h1a1","+ P_{T} > 200 MeV/c","l");
0320     leg1a->AddEntry("h1a2","+ P_{T} > 500 MeV/c","l");
0321     leg1a->Draw();
0322 
0323     //Real-seeded tracks

0324     TCanvas *c1b = new TCanvas("c1b");
0325     h1b->Draw();
0326     h1b1->Draw("same");
0327     h1b2->Draw("same");
0328 
0329     TLegend *leg1b = new TLegend(0.25,0.7,0.5,0.875);
0330     leg1b->SetBorderSize(0);leg1b->SetFillStyle(0);
0331     leg1b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0332     leg1b->AddEntry("h1b","All real-seeded tracks","l");
0333     leg1b->AddEntry("h1b1","+ P_{T} > 200 MeV/c","l");
0334     leg1b->AddEntry("h1b2","+ P_{T} > 500 MeV/c","l");
0335     leg1b->Draw();
0336 
0337     //Truth-seeded tracks

0338     TCanvas *c1c = new TCanvas("c1c");
0339     h1c->Draw();
0340     h1c1->Draw("same");
0341     h1c2->Draw("same");
0342 
0343     TLegend *leg1c = new TLegend(0.25,0.7,0.5,0.875);
0344     leg1c->SetBorderSize(0);leg1c->SetFillStyle(0);
0345     leg1c->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0346     leg1c->AddEntry("h1c","All truth-seeded tracks","l");
0347     leg1c->AddEntry("h1c1","+ P_{T} > 200 MeV/c","l");
0348     leg1c->AddEntry("h1c2","+ P_{T} > 500 MeV/c","l");
0349     leg1c->Draw();
0350 
0351     //Comparison 1

0352     TCanvas *c1d = new TCanvas("c1d");
0353     auto frame_d1 = c1d->DrawFrame(-4,0,4,1.2*h1a1->GetMaximum());
0354     frame_d1->GetXaxis()->SetTitle("#eta_{gen} or #eta_{rec}");frame_d1->GetXaxis()->CenterTitle();
0355     h1a1->Draw("same");
0356     h1b1->Draw("same");
0357     h1c1->Draw("same");
0358 
0359     TLegend *leg1d = new TLegend(0.125,0.675,0.575,0.875);
0360     leg1d->SetBorderSize(0);leg1d->SetFillStyle(0);
0361     leg1d->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0362     leg1d->AddEntry("h1a1","Generated charged particles w/ P_{T} > 200 MeV/c","l");
0363     leg1d->AddEntry("h1b1","Real-seeded tracks w/ P_{T} > 200 MeV/c","l");
0364     leg1d->AddEntry("h1c1","Truth-seeded tracks w/ P_{T} > 200 MeV/c","l");
0365     leg1d->Draw();
0366 
0367     //Comparison 2a

0368     TCanvas *c1e = new TCanvas("c1e");
0369     auto frame_e1 = c1e->DrawFrame(-4,0,4,1.2*h1a1->GetMaximum());
0370     frame_e1->GetXaxis()->SetTitle("#eta_{gen} or #eta_{rec}");frame_e1->GetXaxis()->CenterTitle();
0371     h1a2->Draw("same");
0372     h1b2->Draw("P same");
0373 
0374     TLegend *leg1e = new TLegend(0.125,0.675,0.575,0.875);
0375     leg1e->SetBorderSize(0);leg1e->SetFillStyle(0);
0376     leg1e->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0377     leg1e->AddEntry("h1a2","Generated charged particles w/ P_{T} > 500 MeV/c","fl");
0378     leg1e->AddEntry("h1b2","Real-seeded tracks w/ P_{T} > 500 MeV/c","p");
0379     leg1e->Draw();
0380 
0381     //Comparison 2b

0382     TCanvas *c1e1 = new TCanvas("c1e1");
0383     frame_e1->Draw();
0384     h1a2->Draw("same");
0385     h1c2->Draw("P same");
0386 
0387     TLegend *leg1e1 = new TLegend(0.125,0.675,0.575,0.875);
0388     leg1e1->SetBorderSize(0);leg1e1->SetFillStyle(0);
0389     leg1e1->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0390     leg1e1->AddEntry("h1a2","Generated charged particles w/ P_{T} > 500 MeV/c","fl");
0391     leg1e1->AddEntry("h1c2","Truth-seeded tracks w/ P_{T} > 500 MeV/c","p");
0392     leg1e1->Draw();
0393 
0394     //Comparison 1 ratio

0395     TCanvas *c1f = new TCanvas("c1f");
0396     h1rb1->Draw();
0397     h1rc1->Draw("same");
0398 
0399     TLegend *leg1f = new TLegend(0.575,0.25,0.875,0.45);
0400     leg1f->SetBorderSize(0);leg1f->SetFillStyle(0);
0401     leg1f->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0402     leg1f->AddEntry("h1rb1","Real-seeded tracking","l");
0403     leg1f->AddEntry("h1rc1","Truth-seeded tracking","l");
0404     leg1f->Draw();
0405 
0406     //Comparison 2 ratio

0407     TCanvas *c1g = new TCanvas("c1g");
0408     h1rb2->Draw();
0409     h1rc2->Draw("same");
0410 
0411     TLegend *leg1g = new TLegend(0.575,0.25,0.875,0.45);
0412     leg1g->SetBorderSize(0);leg1g->SetFillStyle(0);
0413     leg1g->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0414     leg1g->AddEntry("h1rb2","Real-seeded tracking","l");
0415     leg1g->AddEntry("h1rc2","Truth-seeded tracking","l");
0416     leg1g->Draw();
0417 
0418     //PT spectra -- Generated charged particles and Reconstructed truth-seeded tracks

0419     TCanvas *c2a = new TCanvas("c2a");
0420     auto frame_2a = c2a->DrawFrame(0,0.1,15,1.2*h2a->GetMaximum());
0421     frame_2a->GetXaxis()->SetTitle("P_{T,gen} or P_{T,rec} [GeV/c]");frame_2a->GetXaxis()->CenterTitle();
0422     gPad->SetLogy();
0423     h2a->Draw("same");
0424     h2b->Draw("P same");
0425 
0426     TLegend *leg2a = new TLegend(0.375,0.675,0.875,0.875);
0427     leg2a->SetBorderSize(0);leg2a->SetFillStyle(0);
0428     leg2a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0429     leg2a->AddEntry("h2a","Generated charged particles w/ -3.5 < #eta < 3.5","fl");
0430     leg2a->AddEntry("h2b","Real-seeded tracks w/ -3.5 < #eta < 3.5","p");
0431     leg2a->Draw();
0432 
0433     //PT spectra -- Generated charged particles and Reconstructed truth-seeded tracks

0434     TCanvas *c2b = new TCanvas("c2b");
0435     frame_2a->Draw();
0436     gPad->SetLogy();
0437     h2a->Draw("same");
0438     h2c->Draw("P same");
0439 
0440     TLegend *leg2b = new TLegend(0.375,0.675,0.875,0.875);
0441     leg2b->SetBorderSize(0);leg2b->SetFillStyle(0);
0442     leg2b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0443     leg2b->AddEntry("h2a","Generated charged particles w/ -3.5 < #eta < 3.5","fl");
0444     leg2b->AddEntry("h2c","Truth-seeded tracks w/ -3.5 < #eta < 3.5","p");
0445     leg2b->Draw();
0446 
0447     //Efficiency as a function of PT

0448     TCanvas *c3a = new TCanvas("c3a");
0449     auto frame_3a = c3a->DrawFrame(0,0,15,1.1);
0450     frame_3a->GetXaxis()->SetTitle("P_{T} [GeV/c]");frame_3a->GetXaxis()->CenterTitle();
0451     frame_3a->GetYaxis()->SetTitle("Efficiency");frame_3a->GetYaxis()->CenterTitle();
0452 
0453     h3ra->Draw("P same");
0454     h3rb->Draw("P same");
0455 
0456     TLegend *leg3a = new TLegend(0.4,0.3,0.8,0.5);
0457     leg3a->SetBorderSize(0);leg3a->SetFillStyle(0);
0458     leg3a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0459     leg3a->AddEntry("h3ra","Real-Seeded tracks","p");
0460     leg3a->AddEntry("h3rb","Truth-seeded tracks","p");
0461     leg3a->Draw();
0462 
0463     //Print plots to file

0464     c1a->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf[",beam_energies.Data(),Q2_set));
0465     c1a->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0466     c1b->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0467     c1c->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0468     c1d->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0469     c1e->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0470     c1e1->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0471     c1f->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0472     c1g->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0473     c2a->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0474     c2b->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0475     c3a->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0476     c3a->Print(Form("plots/DIS_reconstruction_E_%s_Q2_%d.pdf]",beam_energies.Data(),Q2_set));
0477 
0478 }