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 hit_based_matching(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     //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;
0097     if(use_campaign) path = "./campaign_input/";
0098     else path = "./input/";
0099 
0100     TString beam_energies;
0101     if(energy_set == 0) beam_energies = "18x275";
0102     else if(energy_set == 1) beam_energies = "10x100";
0103 
0104     TString run_name;
0105     if(use_campaign) run_name.Form("list_%s_Q2_%d.txt",beam_energies.Data(),Q2_set);
0106     else run_name.Form("%s/Q2_%d/eicrecon_*root",beam_energies.Data(),Q2_set);
0107 
0108     //Open File

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

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

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

0152     TLorentzVector gen_vec;
0153     TVector3 gen_vertex;
0154 
0155     //Counters for truth-seeded associations

0156     int pure_track_counter_ts(0);
0157     int mixed_track_counter_ts(0);
0158 
0159     //Counter for real-seeded associations

0160     int pure_track_counter_rs(0);
0161     int mixed_track_counter_rs(0);
0162 
0163     TLorentzVector rec_vec;
0164     int counter(0);
0165 
0166     //Loop over events

0167     while (tr.Next()) {
0168 
0169         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0170         counter++;
0171 
0172         //Define vector of tuples containing (track index, mc particle status, association weight)

0173         std::vector<std::tuple<int, int, float>> vec;
0174 
0175         //Loop over (truth-seeded) hit-based associations

0176         for(int iassoc=0;iassoc<hit_assoc_weight_ts.GetSize();iassoc++){
0177 
0178             auto assoc_weight = hit_assoc_weight_ts[iassoc];
0179 
0180             h1a->Fill(assoc_weight);
0181 
0182             //Count pure and mixed tracks (avoid double-counting of mixed tracks)

0183             if(assoc_weight > 0.99)
0184                 pure_track_counter_ts++;
0185             else if(iassoc==0)
0186                 mixed_track_counter_ts++;
0187             else if(iassoc>0 && hit_assoc_track_ts[iassoc] != hit_assoc_track_ts[iassoc-1]) //I assume the associated track index is sorted

0188                 mixed_track_counter_ts++;
0189 
0190             //Get the associated track and MC particle

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

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

0205             int status = gen_status[mc_index];
0206 
0207             if(assoc_weight > 0.99){
0208                 h4a->Fill(status);
0209             }
0210             else{ //Fill vector of tuples for mixed tracks

0211                 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0212             }
0213 
0214         }
0215 
0216         //Sort vector by first tuple index (not needed since should already be in ascending order)

0217         std::sort(vec.begin(), vec.end());
0218 
0219         //Group based on the first tuple index (i.e. track index)

0220         size_t ivec = 0;
0221         while (ivec < vec.size()) {
0222 
0223             int current_first = std::get<0>(vec[ivec]);
0224             float current_max_weight = std::get<2>(vec[ivec]);
0225             size_t start = ivec;
0226             size_t index_max_weight = ivec;
0227 
0228             //Find all consecutive elements with the same first value

0229             while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0230                 float current_weight = std::get<2>(vec[ivec]);
0231                 if( current_weight > current_max_weight ){ 
0232                     current_max_weight = current_weight; 
0233                     index_max_weight = ivec;
0234                 }
0235                 ++ivec;
0236             }
0237 
0238             //Fill histograms with MC particle status values

0239             if (ivec - start > 1) {
0240                 for (size_t ii = start; ii < ivec; ++ii){
0241                     
0242                     if(ii == index_max_weight)
0243                         h5a->Fill( std::get<1>(vec[ii]) );
0244                     else
0245                         h6a->Fill( std::get<1>(vec[ii]) );
0246                 }
0247             }
0248         } //End while loop over vec

0249 
0250         //Clear vec and reset ivec for real-seeded tracking use

0251         vec.clear();
0252         ivec = 0;
0253         //End of truth-seeded tracking portion

0254 
0255         //Loop over (real-seeded) hit-based associations

0256         for(int iassoc=0;iassoc<hit_assoc_weight_rs.GetSize();iassoc++){
0257 
0258             auto assoc_weight = hit_assoc_weight_rs[iassoc];
0259 
0260             h1b->Fill(assoc_weight);
0261 
0262             //Count pure and mixed tracks (avoid double-counting of mixed tracks)

0263             if(assoc_weight > 0.99)
0264                 pure_track_counter_rs++;
0265             else if(iassoc==0)
0266                 mixed_track_counter_rs++;
0267             else if(iassoc>0 && hit_assoc_track_rs[iassoc] != hit_assoc_track_rs[iassoc-1]) //I assume the associated track index is sorted

0268                 mixed_track_counter_rs++;
0269  
0270             //Get the associated track and MC particle

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

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

0286             int status = gen_status[mc_index];
0287 
0288             if(assoc_weight > 0.99){
0289                 h4b->Fill(status);
0290             }
0291             else{ //Fill vector of tuples for mixed tracks

0292                 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0293             }
0294  
0295         }
0296 
0297         //Sort vector by first tuple index (not needed since should already be in ascending order)

0298         std::sort(vec.begin(), vec.end());
0299 
0300         //Group based on the first tuple index (i.e. track index)

0301         while (ivec < vec.size()) {
0302 
0303             int current_first = std::get<0>(vec[ivec]);
0304             float current_max_weight = std::get<2>(vec[ivec]);
0305             size_t start = ivec;
0306             size_t index_max_weight = ivec;
0307 
0308             //Find all consecutive elements with the same first value

0309             while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0310                 float current_weight = std::get<2>(vec[ivec]);
0311                 if( current_weight > current_max_weight ){ 
0312                     current_max_weight = current_weight; 
0313                     index_max_weight = ivec;
0314                 }
0315                 ++ivec;
0316             }
0317 
0318             //Fill histograms with MC particle status values

0319             if (ivec - start > 1) {
0320                 for (size_t ii = start; ii < ivec; ++ii){
0321                     
0322                     if(ii == index_max_weight)
0323                         h5b->Fill( std::get<1>(vec[ii]) );
0324                     else
0325                         h6b->Fill( std::get<1>(vec[ii]) );
0326                 }
0327             }
0328         } //End while loop over vec

0329         //End of real-seeded tracking portion

0330 
0331     } //End loop over events

0332 
0333     //Calculate fraction of pure tracks

0334     float frac_pure_ts = ((float) pure_track_counter_ts) / ((float) (pure_track_counter_ts + mixed_track_counter_ts));
0335     float frac_pure_rs = ((float) pure_track_counter_rs) / ((float) (pure_track_counter_rs + mixed_track_counter_rs));
0336 
0337     //Make plots

0338     //----

0339     //Association weights -- truth-seeding

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

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

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

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

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

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

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

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

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

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

0458     c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf[",beam_energies.Data(),Q2_set));
0459     c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0460     c1b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0461     c2a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0462     c2b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0463     c3a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0464     c3b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0465     c4a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0466     c4b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0467     c5a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0468     c6a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0469     c5b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0470     c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0471     c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf]",beam_energies.Data(),Q2_set));
0472 
0473 }