Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:18:32

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

0002 // energy_set 0 is 18x275 GeV

0003 // energy_set 1 is 10x100 GeV

0004 void hit_based_matching(int energy_set = 0, int Q2_set = 1){
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     //Track purity

0023     TH1 *h1a = new TH1D("h1a","Truth-seeded tracks",100,0,1.1);
0024     h1a->GetXaxis()->SetTitle("Fraction of track measurements from a given MC Particle");h1a->GetXaxis()->CenterTitle();
0025     h1a->GetYaxis()->SetTitle("Number of Tracks"); h1a->GetYaxis()->CenterTitle();
0026     h1a->SetLineColor(kBlack);h1a->SetLineWidth(2);
0027 
0028     TH1 *h1b = new TH1D("h1b","Real-seeded tracks",100,0,1.1);
0029     h1b->GetXaxis()->SetTitle("Fraction of track measurements from a given MC Particle");h1b->GetXaxis()->CenterTitle();
0030     h1b->GetYaxis()->SetTitle("Number of Tracks"); h1b->GetYaxis()->CenterTitle();
0031     h1b->SetLineColor(kBlack);h1b->SetLineWidth(2);
0032 
0033     //Momentum reconstruction

0034     TH2 *h2a = new TH2D("h2a","Truth-seeded tracks",100,0,50,100,0,50);
0035     h2a->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h2a->GetXaxis()->CenterTitle();
0036     h2a->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h2a->GetYaxis()->CenterTitle();
0037 
0038     TH2 *h2b = new TH2D("h2b","Real-seeded tracks",100,0,50,100,0,50);
0039     h2b->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h2b->GetXaxis()->CenterTitle();
0040     h2b->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h2b->GetYaxis()->CenterTitle();
0041 
0042     //Momentum reconstruction -- -3 < eta_gen < +3

0043     TH2 *h3a = new TH2D("h3a","Truth-seeded tracks",100,0,50,100,0,50);
0044     h3a->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h3a->GetXaxis()->CenterTitle();
0045     h3a->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h3a->GetYaxis()->CenterTitle();
0046 
0047     TH2 *h3b = new TH2D("h3b","Real-seeded tracks",100,0,50,100,0,50);
0048     h3b->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h3b->GetXaxis()->CenterTitle();
0049     h3b->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h3b->GetYaxis()->CenterTitle();
0050 
0051     //Status of associated particles for entirely pure tracks

0052     TH1 *h4a = new TH1D("h4a","Truth-seeded tracks: Entirely pure tracks",2,0,2);
0053     h4a->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0054     h4a->GetXaxis()->SetBinLabel(2,"Primary Particle");
0055     h4a->GetXaxis()->SetTitle("Status of associated MC particle");h4a->GetXaxis()->CenterTitle();
0056     h4a->GetYaxis()->SetTitle("Number of Tracks"); h4a->GetYaxis()->CenterTitle();
0057     h4a->SetLineColor(kBlack);h4a->SetLineWidth(2);
0058 
0059     TH1 *h4b = new TH1D("h4b","Real-seeded tracks: Entirely pure tracks",2,0,2);
0060     h4b->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0061     h4b->GetXaxis()->SetBinLabel(2,"Primary Particle");
0062     h4b->GetXaxis()->SetTitle("Status of associated MC particle");h4b->GetXaxis()->CenterTitle();
0063     h4b->GetYaxis()->SetTitle("Number of Tracks"); h4b->GetYaxis()->CenterTitle();
0064     h4b->SetLineColor(kBlack);h4b->SetLineWidth(2);
0065 
0066     //Status of associated particle for mixed tracks

0067     TH1 *h5a = new TH1D("h5a","Truth-seeded tracks: Status of most-associated particle for mixed tracks",2,0,2);
0068     h5a->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0069     h5a->GetXaxis()->SetBinLabel(2,"Primary Particle");
0070     h5a->GetXaxis()->SetTitle("Status of associated MC particle");h5a->GetXaxis()->CenterTitle();
0071     h5a->GetYaxis()->SetTitle("Number of Tracks"); h5a->GetYaxis()->CenterTitle();
0072     h5a->SetLineColor(kBlack);h5a->SetLineWidth(2);
0073 
0074     TH1 *h6a = new TH1D("h6a","Truth-seeded tracks: Status of less-associated particles for mixed tracks",2,0,2);
0075     h6a->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0076     h6a->GetXaxis()->SetBinLabel(2,"Primary Particle");
0077     h6a->GetXaxis()->SetTitle("Status of associated MC particle");h6a->GetXaxis()->CenterTitle();
0078     h6a->GetYaxis()->SetTitle("Number of Tracks"); h6a->GetYaxis()->CenterTitle();
0079     h6a->SetLineColor(kBlack);h6a->SetLineWidth(2);
0080 
0081     TH1 *h5b = new TH1D("h5b","Real-seeded tracks: Status of most-associated particle for mixed tracks",2,0,2);
0082     h5b->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0083     h5b->GetXaxis()->SetBinLabel(2,"Primary Particle");
0084     h5b->GetXaxis()->SetTitle("Status of associated MC particle");h5b->GetXaxis()->CenterTitle();
0085     h5b->GetYaxis()->SetTitle("Number of Tracks"); h5b->GetYaxis()->CenterTitle();
0086     h5b->SetLineColor(kBlack);h5b->SetLineWidth(2);
0087 
0088     TH1 *h6b = new TH1D("h6b","Real-seeded tracks: Status of less-associated particles for mixed tracks",2,0,2);
0089     h6b->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0090     h6b->GetXaxis()->SetBinLabel(2,"Primary Particle");
0091     h6b->GetXaxis()->SetTitle("Status of associated MC particle");h6b->GetXaxis()->CenterTitle();
0092     h6b->GetYaxis()->SetTitle("Number of Tracks"); h6b->GetYaxis()->CenterTitle();
0093     h6b->SetLineColor(kBlack);h6b->SetLineWidth(2);
0094 
0095     //File information

0096     TString path = "./input/";
0097     TString run_name;
0098     TString beam_energies;
0099 
0100     if(energy_set == 0 && Q2_set == 1){
0101         beam_energies = "18x275";
0102         run_name = "Q2_1/eicrecon_*.root"; //Pythia8 DIS events with 18x275 GeV and Q2>1 GeV2

0103     }
0104     else if(energy_set == 0 && Q2_set == 10){
0105         beam_energies = "18x275";
0106         run_name = "eicrecon_out_E_18_275_Q2_10.root"; //Pythia8 DIS events with 18x275 GeV and Q2>10 GeV2

0107     }
0108     else if(energy_set == 0 && Q2_set == 100){
0109         beam_energies = "18x275";
0110         run_name = "Q2_100/eicrecon_*.root"; //Pythia8 DIS events with 18x275 GeV and Q2>100 GeV2

0111     }
0112     else if(energy_set == 1 && Q2_set == 1){
0113         beam_energies = "10x100";
0114         run_name = "eicrecon_out_E_10_100_Q2_1.root"; //Pythia8 DIS events with 10x100 GeV and Q2>1 GeV2

0115     }
0116     else if(energy_set == 1 && Q2_set == 10){
0117         beam_energies = "10x100";
0118         run_name = "eicrecon_out_E_10_100_Q2_10.root"; //Pythia8 DIS events with 10x100 GeV and Q2>10 GeV2

0119     }
0120     else if(energy_set == 1 && Q2_set == 100){
0121         beam_energies = "10x100";
0122         run_name = "eicrecon_out_E_10_100_Q2_100.root"; //Pythia8 DIS events with 10x100 GeV and Q2>100 GeV2

0123     }
0124 
0125     //Open File

0126     TString input = path + run_name;
0127     TChain *tree = new TChain("events");
0128     tree->Add(input.Data());
0129 
0130     cout<<"Analyzing "<<tree->GetEntries()<<" events!"<<endl;
0131 
0132     //Create Array Reader

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

0141     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0142     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0143     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0144     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0145 
0146     TTreeReaderArray<float> hit_assoc_weight_ts(tr,"CentralCKFTruthSeededTrackAssociations.weight");
0147     TTreeReaderArray<int>   hit_assoc_mc_ts(tr, "_CentralCKFTruthSeededTrackAssociations_sim.index");
0148     TTreeReaderArray<int>   hit_assoc_track_ts(tr, "_CentralCKFTruthSeededTrackAssociations_rec.index");
0149 
0150     TTreeReaderArray<float> hit_assoc_weight_rs(tr,"CentralCKFTrackAssociations.weight");
0151     TTreeReaderArray<int>   hit_assoc_mc_rs(tr, "_CentralCKFTrackAssociations_sim.index");
0152     TTreeReaderArray<int>   hit_assoc_track_rs(tr, "_CentralCKFTrackAssociations_rec.index");
0153 
0154     TTreeReaderArray<float> track_qoverp_ts(tr, "CentralCKFTruthSeededTrackParameters.qOverP");
0155     TTreeReaderArray<float> track_qoverp_rs(tr, "CentralCKFTrackParameters.qOverP");
0156     
0157     //Other variables

0158     TLorentzVector gen_vec;
0159     TVector3 gen_vertex;
0160 
0161     //Counters for truth-seeded associations

0162     int pure_track_counter_ts(0);
0163     int mixed_track_counter_ts(0);
0164 
0165     //Counter for real-seeded associations

0166     int pure_track_counter_rs(0);
0167     int mixed_track_counter_rs(0);
0168 
0169     TLorentzVector rec_vec;
0170     int counter(0);
0171 
0172     //Loop over events

0173     while (tr.Next()) {
0174 
0175         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0176         counter++;
0177 
0178         //Define vector of tuples containing (track index, mc particle status, association weight)

0179         std::vector<std::tuple<int, int, float>> vec;
0180 
0181         //Loop over (truth-seeded) hit-based associations

0182         for(int iassoc=0;iassoc<hit_assoc_weight_ts.GetSize();iassoc++){
0183 
0184             auto assoc_weight = hit_assoc_weight_ts[iassoc];
0185 
0186             h1a->Fill(assoc_weight);
0187 
0188             //Count pure and mixed tracks (avoid double-counting of mixed tracks)

0189             if(assoc_weight > 0.99)
0190                 pure_track_counter_ts++;
0191             else if(iassoc==0)
0192                 mixed_track_counter_ts++;
0193             else if(iassoc>0 && hit_assoc_track_ts[iassoc] != hit_assoc_track_ts[iassoc-1]) //I assume the associated track index is sorted

0194                 mixed_track_counter_ts++;
0195 
0196             //Get the associated track and MC particle

0197             //N.B. We use the fact that (right now) the track collection index is the same as the track parameter collection index

0198             int mc_index = hit_assoc_mc_ts[iassoc];
0199             int trk_index = hit_assoc_track_ts[iassoc];
0200 
0201             gen_vec.SetXYZM(gen_px[mc_index],gen_py[mc_index],gen_pz[mc_index],gen_mass[mc_index]);
0202             auto gen_mom = gen_vec.P();
0203 
0204             auto trk_mom = fabs(1./track_qoverp_ts[trk_index]);
0205 
0206             h2a->Fill(gen_mom,trk_mom);
0207             if( fabs(gen_vec.Eta()) < 3 )
0208                 h3a->Fill(gen_mom,trk_mom);
0209 
0210             //Check status (i.e. primary or secondary) of associated particle

0211             int status = gen_status[mc_index];
0212 
0213             if(assoc_weight > 0.99){
0214                 h4a->Fill(status);
0215             }
0216             else{ //Fill vector of tuples for mixed tracks

0217                 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0218             }
0219 
0220         }
0221 
0222         //Sort vector by first tuple index (not needed since should already be in ascending order)

0223         std::sort(vec.begin(), vec.end());
0224 
0225         //Group based on the first tuple index (i.e. track index)

0226         size_t ivec = 0;
0227         while (ivec < vec.size()) {
0228 
0229             int current_first = std::get<0>(vec[ivec]);
0230             float current_max_weight = std::get<2>(vec[ivec]);
0231             size_t start = ivec;
0232             size_t index_max_weight = ivec;
0233 
0234             //Find all consecutive elements with the same first value

0235             while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0236                 float current_weight = std::get<2>(vec[ivec]);
0237                 if( current_weight > current_max_weight ){ 
0238                     current_max_weight = current_weight; 
0239                     index_max_weight = ivec;
0240                 }
0241                 ++ivec;
0242             }
0243 
0244             //Fill histograms with MC particle status values

0245             if (ivec - start > 1) {
0246                 for (size_t ii = start; ii < ivec; ++ii){
0247                     
0248                     if(ii == index_max_weight)
0249                         h5a->Fill( std::get<1>(vec[ii]) );
0250                     else
0251                         h6a->Fill( std::get<1>(vec[ii]) );
0252                 }
0253             }
0254         } //End while loop over vec

0255 
0256         //Clear vec and reset ivec for real-seeded tracking use

0257         vec.clear();
0258         ivec = 0;
0259         //End of truth-seeded tracking portion

0260 
0261         //Loop over (real-seeded) hit-based associations

0262         for(int iassoc=0;iassoc<hit_assoc_weight_rs.GetSize();iassoc++){
0263 
0264             auto assoc_weight = hit_assoc_weight_rs[iassoc];
0265 
0266             h1b->Fill(assoc_weight);
0267 
0268             //Count pure and mixed tracks (avoid double-counting of mixed tracks)

0269             if(assoc_weight > 0.99)
0270                 pure_track_counter_rs++;
0271             else if(iassoc==0)
0272                 mixed_track_counter_rs++;
0273             else if(iassoc>0 && hit_assoc_track_rs[iassoc] != hit_assoc_track_rs[iassoc-1]) //I assume the associated track index is sorted

0274                 mixed_track_counter_rs++;
0275  
0276             //Get the associated track and MC particle

0277             //N.B. We use the fact that (right now) the track collection index is the same as the track parameter collection index

0278             int mc_index = hit_assoc_mc_rs[iassoc];
0279             int trk_index = hit_assoc_track_rs[iassoc];
0280 
0281             gen_vec.SetXYZM(gen_px[mc_index],gen_py[mc_index],gen_pz[mc_index],gen_mass[mc_index]);
0282             auto gen_mom = gen_vec.P();
0283 
0284             auto trk_mom = fabs(1./track_qoverp_rs[trk_index]);
0285 
0286             h2b->Fill(gen_mom,trk_mom);
0287 
0288             if( fabs(gen_vec.Eta()) < 3 )
0289                 h3b->Fill(gen_mom,trk_mom);
0290 
0291             //Check status (i.e. primary or secondary) of associated particle

0292             int status = gen_status[mc_index];
0293 
0294             if(assoc_weight > 0.99){
0295                 h4b->Fill(status);
0296             }
0297             else{ //Fill vector of tuples for mixed tracks

0298                 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0299             }
0300  
0301         }
0302 
0303         //Sort vector by first tuple index (not needed since should already be in ascending order)

0304         std::sort(vec.begin(), vec.end());
0305 
0306         //Group based on the first tuple index (i.e. track index)

0307         while (ivec < vec.size()) {
0308 
0309             int current_first = std::get<0>(vec[ivec]);
0310             float current_max_weight = std::get<2>(vec[ivec]);
0311             size_t start = ivec;
0312             size_t index_max_weight = ivec;
0313 
0314             //Find all consecutive elements with the same first value

0315             while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0316                 float current_weight = std::get<2>(vec[ivec]);
0317                 if( current_weight > current_max_weight ){ 
0318                     current_max_weight = current_weight; 
0319                     index_max_weight = ivec;
0320                 }
0321                 ++ivec;
0322             }
0323 
0324             //Fill histograms with MC particle status values

0325             if (ivec - start > 1) {
0326                 for (size_t ii = start; ii < ivec; ++ii){
0327                     
0328                     if(ii == index_max_weight)
0329                         h5b->Fill( std::get<1>(vec[ii]) );
0330                     else
0331                         h6b->Fill( std::get<1>(vec[ii]) );
0332                 }
0333             }
0334         } //End while loop over vec

0335         //End of real-seeded tracking portion

0336 
0337     } //End loop over events

0338 
0339     //Calculate fraction of pure tracks

0340     float frac_pure_ts = ((float) pure_track_counter_ts) / ((float) (pure_track_counter_ts + mixed_track_counter_ts));
0341     float frac_pure_rs = ((float) pure_track_counter_rs) / ((float) (pure_track_counter_rs + mixed_track_counter_rs));
0342 
0343     //Make plots

0344     //----

0345     //Association weights -- truth-seeding

0346     TCanvas *c1a = new TCanvas("c1a");
0347     c1a->SetLogy();
0348     h1a->Draw();
0349 
0350     TLegend *leg1a = new TLegend(0.125,0.7,0.475,0.875);
0351     leg1a->SetBorderSize(0);leg1a->SetFillStyle(0);
0352     leg1a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0353     leg1a->AddEntry(h1a, Form("Fraction of entirely pure tracks = %.3f",frac_pure_ts),"l");
0354     leg1a->Draw();
0355 
0356     //Association weights -- real-seeding

0357     TCanvas *c1b = new TCanvas("c1b");
0358     c1b->SetLogy();
0359     h1b->Draw();
0360 
0361     TLegend *leg1b = new TLegend(0.125,0.7,0.475,0.875);
0362     leg1b->SetBorderSize(0);leg1b->SetFillStyle(0);
0363     leg1b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0364     leg1b->AddEntry(h1b, Form("Fraction of entirely pure tracks = %.3f",frac_pure_rs),"l");
0365     leg1b->Draw();
0366 
0367     //Track vs. Associated MC particle Momentum -- truth-seeding

0368     TCanvas *c2a = new TCanvas("c2a");
0369     c2a->SetLogz();
0370     h2a->Draw("colz");
0371 
0372     TLatex *tex2a = new TLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0373     tex2a->SetTextSize(0.03);
0374     tex2a->Draw();
0375 
0376     //Track vs. Associated MC particle Momentum -- real-seeding

0377     TCanvas *c2b = new TCanvas("c2b");
0378     c2b->SetLogz();
0379     h2b->Draw("colz");
0380 
0381     TLatex *tex2b = new TLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0382     tex2b->SetTextSize(0.03);
0383     tex2b->Draw();
0384 
0385     //Track vs. Associated MC particle Momentum -- truth-seeding, eta_gen restriction

0386     TCanvas *c3a = new TCanvas("c3a");
0387     c3a->SetLogz();
0388     h3a->Draw("colz");
0389 
0390     TLatex *tex3a = new TLatex();
0391     tex3a->SetTextSize(0.03);
0392     tex3a->DrawLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0393     tex3a->DrawLatex(25,5," -3 < #eta_{gen} < +3");
0394 
0395     //Track vs. Associated MC particle Momentum -- real-seeding, eta_gen restriction

0396     TCanvas *c3b = new TCanvas("c3b");
0397     c3b->SetLogz();
0398     h3b->Draw("colz");
0399 
0400     TLatex *tex3b = new TLatex();
0401     tex3b->SetTextSize(0.03);
0402     tex3b->DrawLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0403     tex3b->DrawLatex(25,5," -3 < #eta_{gen} < +3");
0404 
0405     //Status of associated particles for entirely pure tracks -- truth seeding

0406     TCanvas *c4a = new TCanvas("c4a");
0407     c4a->SetLogy();
0408     h4a->Draw();
0409 
0410     TLegend *leg4a = new TLegend(0.125,0.7,0.475,0.875);
0411     leg4a->SetBorderSize(0);leg4a->SetFillStyle(0);
0412     leg4a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0413     leg4a->Draw();
0414 
0415     //Status of associated particles for entirely pure tracks -- real seeding

0416     TCanvas *c4b = new TCanvas("c4b");
0417     c4b->SetLogy();
0418     h4b->Draw();
0419 
0420     TLegend *leg4b = new TLegend(0.125,0.7,0.475,0.875);
0421     leg4b->SetBorderSize(0);leg4b->SetFillStyle(0);
0422     leg4b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0423     leg4b->Draw();
0424 
0425     //Status of associated particles for mixed tracks -- truth seeding

0426     TCanvas *c5a = new TCanvas("c5a");
0427     c5a->SetLogy();
0428     h5a->Draw();
0429 
0430     TLegend *leg5a = new TLegend(0.125,0.7,0.475,0.875);
0431     leg5a->SetBorderSize(0);leg5a->SetFillStyle(0);
0432     leg5a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0433     leg5a->Draw();
0434 
0435     TCanvas *c6a = new TCanvas("c6a");
0436     c6a->SetLogy();
0437     h6a->Draw();
0438 
0439     TLegend *leg6a = new TLegend(0.525,0.7,0.875,0.875);
0440     leg6a->SetBorderSize(0);leg6a->SetFillStyle(0);
0441     leg6a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0442     leg6a->Draw();
0443 
0444     //Status of associated particles for mixed tracks -- real seeding

0445     TCanvas *c5b = new TCanvas("c5b");
0446     c5b->SetLogy();
0447     h5b->Draw();
0448 
0449     TLegend *leg5b = new TLegend(0.125,0.7,0.475,0.875);
0450     leg5b->SetBorderSize(0);leg5b->SetFillStyle(0);
0451     leg5b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0452     leg5b->Draw();
0453 
0454     TCanvas *c6b = new TCanvas("c6b");
0455     c6b->SetLogy();
0456     h6b->Draw();
0457 
0458     TLegend *leg6b = new TLegend(0.525,0.7,0.875,0.875);
0459     leg6b->SetBorderSize(0);leg6b->SetFillStyle(0);
0460     leg6b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0461     leg6b->Draw();
0462 
0463     //Print plots to file

0464     c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf[",beam_energies.Data(),Q2_set));
0465     c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0466     c1b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0467     c2a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0468     c2b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0469     c3a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0470     c3b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0471     c4a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0472     c4b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0473     c5a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0474     c6a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0475     c5b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0476     c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0477     c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf]",beam_energies.Data(),Q2_set));
0478 
0479 }