Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2026-04-06 08:09:42

0001 #pragma once 
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 
0004 #include <TFile.h>
0005 #include <TDirectory.h>
0006 #include <TLegend.h>
0007 #include <TCanvas.h>
0008 #include <TH1.h>
0009 #include <TH2.h>
0010 #include "HerwigQAPlottingConfig.h"
0011 #include "HerwigQAPlottingConfig.C"
0012 #include "Special_colors.h"
0013 #include "sPhenixStyle.h"
0014 #include "sPhenixStyle.C"
0015 
0016 #include <string>
0017 #include <vector>
0018 #include <format>
0019 #include <map>
0020 #include <array>
0021 #include <utility>
0022 
0023 HerwigQAPlottingConfig* conf;
0024 TFile* storage;
0025 std::string tag {};
0026 Skaydis_colors* style_points;
0027 struct particle_colls
0028 {
0029     TH1F* pt;
0030     TH1I* n;
0031     TH2F* phi_eta;
0032 };
0033 void GetSampleSpectra(
0034         std::string Generator="Herwig",
0035             std::vector<TH1F*>* r04_all = new std::vector<TH1F*> {}, 
0036         std::vector<TH1F*>* r04_lead = new std::vector<TH1F*> {},
0037         std::vector<TFile*>* fs = new std::vector<TFile*> {}
0038         )
0039 {
0040     bool isHerwig = false;
0041     if(Generator.find("erwig") != std::string::npos) isHerwig=true;
0042     for(int i = 0; i<(int) fs->size(); i++)
0043     {
0044         fs->at(i)->cd();
0045         TDirectory* d_gen   = (TDirectory*) fs->at(i)->GetDirectory(std::format("{}_jets", Generator).c_str());
0046         
0047         TH1F* h_all     = (TH1F*)d_gen->Get("h_jet_r04_pt");    
0048         TH1F* h_lead    = (TH1F*)d_gen->Get("h_lead_jet_r04_pt");   
0049         if(isHerwig){
0050             h_all->SetMarkerStyle(20+i);
0051             if(i==4) h_all->SetMarkerStyle(29);
0052             else if(i==5) h_all->SetMarkerStyle(33);
0053             else if(i==6) h_all->SetMarkerStyle(34);
0054             else if (i==7) h_all->SetMarkerStyle(37);
0055 
0056             h_lead->SetMarkerStyle(20+i);
0057             if(i==4) h_lead->SetMarkerStyle(29);
0058             else if(i==5) h_lead->SetMarkerStyle(33);
0059             else if(i==6) h_lead->SetMarkerStyle(35);
0060             else if (i==7) h_lead->SetMarkerStyle(37);
0061         }
0062         else{   
0063             h_all->SetMarkerStyle(i+24);
0064             h_all->SetMarkerSize(2);
0065             if(i==6) h_all->SetMarkerStyle(30);
0066             if(i==7) h_all->SetMarkerStyle(38);
0067 
0068             h_lead->SetMarkerStyle(i+24);
0069             h_lead->SetMarkerSize(2);
0070             if(i==6) h_lead->SetMarkerStyle(30);
0071             if(i==7) h_lead->SetMarkerStyle(38);
0072         }   
0073         r04_all->push_back(h_all);
0074         r04_lead->push_back(h_lead);
0075     }
0076     return;
0077 }
0078 void SampleRatios(
0079         std::vector<TH1F*>* base, std::vector<TH1F*>* ratios, TH1F* total, 
0080         TCanvas* ratios_canvi, TLegend* head_leg, TLegend* data_leg, 
0081         HerwigQAPlottingConfig* co, 
0082         std::string Generator, std::vector<std::string> cut
0083 )
0084 {
0085     data_leg->SetNColumns(2);
0086     std::vector<TPad*>* pads = co->AddPads(ratios_canvi);
0087     for(int i = 1; i<(int)base->size(); i++)
0088     {
0089         if ( i > 0 ) 
0090         {
0091             pads->at(1)->cd();
0092             TH1F* last = (TH1F*)base->at(i-1)->Clone();
0093             TH1F* prev_ratio = co->GetRatioPlot(base->at(i), last);
0094             data_leg->AddEntry( prev_ratio, std::format("{} {} / {}", Generator, cut.at(i), cut.at(i-1)).c_str(), "pl");
0095             ratios->push_back(prev_ratio);
0096             prev_ratio->Draw("pmc plc same");
0097             prev_ratio->SetYTitle(" Sample / Previous Sample " );
0098             prev_ratio->GetYaxis()->SetRangeUser(0, 4);
0099             pads->at(0)->cd();
0100         }
0101             TH1F* total_ratio = co->GetRatioPlot(base->at(i), total);
0102             total_ratio->SetMarkerStyle(2*i+39);
0103             data_leg->AddEntry( total_ratio, std::format("{} {} / Combined Spectrum ", Generator, cut.at(i)).c_str(), "pl");
0104             total_ratio->Draw("pmc plc same");
0105             total_ratio->SetYTitle(" Sample / Combined " );
0106             total_ratio->GetYaxis()->SetRangeUser(0, 1);
0107     }
0108     pads->at(1)->cd();
0109     data_leg->Draw();
0110     head_leg->Draw();
0111     ratios_canvi->cd();
0112     pads->at(1)->Draw();
0113     pads->at(0)->Draw();
0114     return;
0115                     
0116 
0117 }
0118 void PlotThePlots(
0119         std::vector<TH1F*>* herwig_all, 
0120             std::vector<TH1F*>* pythia_all,  
0121         std::array<TCanvas*, 3>* AllCanvi,
0122             std::array<TLegend*, 3>* AllHeadLeg, std::array<TLegend*, 3>* AllDataLeg,
0123         HerwigQAPlottingConfig* co,
0124         std::vector<std::string> cuts
0125     )
0126 {
0127     
0128     std::vector<TPad*>* pall    = co->AddPads(AllCanvi->at(0));
0129     for(int i = 0 ; i<(int) herwig_all->size(); i++)
0130     {
0131         TH1F* hCo = (TH1F*) herwig_all->at(i)->Clone();
0132         TH1F* pCo = (TH1F*) pythia_all->at(i)->Clone();
0133         TH1F* rCo = co->GetRatioPlot(hCo, pCo);
0134         rCo -> SetMarkerStyle(2*i+39);
0135         pall->at(1)->cd();
0136         pall->at(1)->SetLogy();
0137         hCo->Draw("same plc pmc");
0138         if(i==0){
0139             float max_val=hCo->GetMaximum();
0140             hCo->GetYaxis()->SetRangeUser(0.001, max_val*3);
0141         }
0142         pCo->Draw("same plc pmc");
0143         AllDataLeg->at(0)->AddEntry(hCo, std::format("Herwig {}", cuts.at(i)).c_str(), "pl");
0144         AllDataLeg->at(0)->AddEntry(pCo, std::format("Pythia {}", cuts.at(i)).c_str(), "pl");
0145         AllDataLeg->at(0)->AddEntry(rCo, std::format("Herwig {} / Pythia {}", cuts.at(i), cuts.at(i)).c_str(), "pl");
0146         pall->at(0)->cd();
0147         rCo->Draw("plc pmc same");
0148         rCo->GetYaxis()->SetRangeUser(0, 2);
0149         rCo->SetYTitle("Herwig / Pythia");
0150         
0151         AllCanvi->at(1)->cd();
0152         herwig_all->at(i)->Draw("same plc pmc");
0153         if(i==0)
0154         {
0155             float max_val=herwig_all->at(i)->GetMaximum();
0156             herwig_all->at(i)->GetYaxis()->SetRangeUser(0.001, max_val*3);
0157         }
0158         AllDataLeg->at(1)->AddEntry(herwig_all->at(i), cuts.at(i).c_str(), "pl");
0159         
0160         AllCanvi->at(2)->cd();
0161         pythia_all->at(i)->Draw("same plc pmc");
0162         if(i==0)
0163         {
0164             float max_val=pythia_all->at(i)->GetMaximum();
0165             pythia_all->at(i)->GetYaxis()->SetRangeUser(0.001, max_val*3);
0166         }
0167         AllDataLeg->at(2)->AddEntry(pythia_all->at(i), cuts.at(i).c_str(), "pl");
0168     }
0169     pall->at(1)->cd();
0170     AllHeadLeg->at(0)->Draw();
0171     AllDataLeg->at(0)->Draw();
0172     AllCanvi->at(0)->cd();
0173     pall->at(0)->Draw();
0174     pall->at(1)->Draw();
0175 
0176     AllCanvi->at(1)->cd();
0177     AllHeadLeg->at(1)->Draw();
0178     AllDataLeg->at(1)->Draw();
0179     AllCanvi->at(1)->SetLogy();
0180 
0181     AllCanvi->at(2)->cd();
0182     AllHeadLeg->at(2)->Draw();
0183     AllDataLeg->at(2)->Draw();
0184     AllCanvi->at(2)->SetLogy();
0185 
0186     return;
0187 }
0188 
0189 void PlotBothThePlots(
0190         std::vector<TH1F*>* herwig_all, std::vector<TH1F*>* herwig_lead,
0191             std::vector<TH1F*>* pythia_all, std::vector<TH1F*>* pythia_lead, 
0192         std::array<TCanvas*, 3>* AllCanvi, std::array<TCanvas*, 3>* LeadCanvi,
0193         HerwigQAPlottingConfig* co,
0194         std::vector<std::string> cuts
0195     )
0196 {
0197     std::array<TLegend*, 3>* AllHeadLeg     = new std::array<TLegend*, 3> {};
0198         std::array<TLegend*, 3>* AllDataLeg     = new std::array<TLegend*, 3> {};
0199     std::array<TLegend*, 3>* LeadHeadLeg    = new std::array<TLegend*, 3> {};
0200     std::array<TLegend*, 3>* LeadDataLeg    = new std::array<TLegend*, 3> {};
0201     for(int i=0; i<3; i++)
0202     {
0203         TLegend* l_head = new TLegend(0.6, 0.7, 0.8, 1); 
0204         TLegend* l_data = new TLegend(0.6, 0.4, 0.8, 0.7);
0205         co->SetLegend(l_data);
0206         if(i == 0) co->SetsPhenixHeaderLegend(l_head, "none", "All Jets");
0207         else if(i == 1) co->SetsPhenixHeaderLegend(l_head, true, "All Jets");
0208         else if(i == 2) co->SetsPhenixHeaderLegend(l_head, false, "All Jets");
0209         l_head->SetTextSize(0.04f);
0210         AllHeadLeg->at(i)   = l_head;
0211         AllDataLeg->at(i)   = l_data;
0212     }
0213 
0214     for(int i=0; i<3; i++)
0215     {
0216         TLegend* l_head = new TLegend(0.6, 0.7, 0.8, 1); 
0217         TLegend* l_data = new TLegend(0.6, 0.4, 0.8, 0.7);
0218         co->SetLegend(l_data);
0219         if(i == 0)  co->SetsPhenixHeaderLegend(l_head, "none", "Lead Jets");
0220         else if(i == 1) co->SetsPhenixHeaderLegend(l_head, true, "Lead Jets");
0221         else if(i == 2) co->SetsPhenixHeaderLegend(l_head, false, "Lead Jets");
0222         LeadHeadLeg->at(i)  = l_head;
0223         LeadDataLeg->at(i)  = l_data;
0224     }
0225     AllDataLeg->at(0)->SetNColumns(3);
0226     LeadDataLeg->at(0)->SetNColumns(3);
0227     PlotThePlots(herwig_all, pythia_all, AllCanvi, AllHeadLeg, AllDataLeg, co, cuts);
0228     PlotThePlots(herwig_lead, pythia_lead, LeadCanvi, LeadHeadLeg, LeadDataLeg, co, cuts);
0229     return;
0230 }
0231 
0232 void PlotCombinedSpectrum(std::vector<TFile*>* fs, std::vector<std::string> cutnames)
0233 {
0234     SetsPhenixStyle();
0235     TCanvas* cAll       = new TCanvas("all_canvas", "all_Jets");
0236     TCanvas* cLead      = new TCanvas("lead_canvas", "Lead_jets");
0237     
0238     TCanvas* cAllPythia = new TCanvas("all_Pythia_canvas", "all_Pythia_Jets");
0239     TCanvas* cLeadPythia    = new TCanvas("lead_Pythia_canvas", "Lead_Pythia_jets");
0240     
0241     TCanvas* cAllHerwig = new TCanvas("all_Herwig_canvas", "all_Herwig_Jets");
0242     TCanvas* cLeadHerwig    = new TCanvas("lead_Herwig_canvas", "Lead_Herwig_jets");
0243     
0244 
0245     
0246     TCanvas* cAll_ratio_Pythia  = new TCanvas("all_ratio_Pythia_canvas", "all_ratio_Pythia_Jets");
0247     TCanvas* cLead_ratio_Pythia = new TCanvas("lead_ratio_Pythia_canvas", "Lead_ratio_Pythia_jets");
0248     
0249     TCanvas* cAll_ratio_Herwig  = new TCanvas("all_ratio_Herwig_canvas", "all_ratio_Herwig_Jets");
0250     TCanvas* cLead_ratio_Herwig = new TCanvas("lead_ratio_Herwig_canvas", "Lead_ratio_Herwig_jets");
0251     
0252 
0253     std::array<TCanvas*, 3> All     {cAll, cAllHerwig, cAllPythia}; 
0254     std::array<TCanvas*, 3> Lead    {cLead, cLeadHerwig, cLeadPythia};
0255     
0256     std::array<TCanvas*, 2> Herwig_ratios {cAll_ratio_Herwig, cLead_ratio_Herwig};
0257     std::array<TCanvas*, 2> Pythia_ratios {cAll_ratio_Pythia, cLead_ratio_Pythia};
0258     
0259     std::array< std::array< TCanvas*, 2 >*, 2 > RatioCanvi {&Herwig_ratios, &Pythia_ratios};    
0260     std::array< std::array< std::array<TLegend*, 2>*, 2 >*, 2 >* RatioLegs 
0261         = new std::array< std::array< std::array< TLegend*, 2>*, 2>*, 2>{}; 
0262     
0263     std::vector<TH1F*>* herwig_pt   = new std::vector<TH1F*>();
0264     std::vector<TH1F*>* herwig_l_pt = new std::vector<TH1F*>();
0265     
0266     std::vector<TH1F*>* pythia_pt   = new std::vector<TH1F*>();
0267     std::vector<TH1F*>* pythia_l_pt = new std::vector<TH1F*>();
0268     
0269     HerwigQAPlottingConfig* co = new HerwigQAPlottingConfig();
0270     Skaydis_colors* odd_colors = new Skaydis_colors();
0271     gStyle->SetPalette(100, odd_colors->Enby_gradient_PT);
0272 
0273 
0274     GetSampleSpectra("Herwig", herwig_pt, herwig_l_pt, fs);
0275     GetSampleSpectra("Pythia", pythia_pt, pythia_l_pt, fs);
0276     
0277     for(int i=0; i<2; i++)
0278     {
0279         bool isHerwig       = false;
0280         if( i==0 ) isHerwig = true;
0281         else isHerwig       = false;
0282         std::array< std::array<TLegend*, 2>*, 2 >* Legs = 
0283                 new std::array<std::array<TLegend*, 2>*, 2>{};
0284         for(int j=0; j<2; j++)
0285         {
0286             RatioCanvi.at(i)->at(j)->cd();
0287             std::string LJ = "All Jets";
0288             if(j==1) LJ="Lead Jets";
0289             TLegend* l_head = new TLegend(0.05, 0.8, 0.3, 1);
0290             TLegend* l_data = new TLegend(0.32, 0.7, 1, 0.95);
0291                     
0292             co->SetLegend(l_data);
0293             co->SetsPhenixHeaderLegend(l_head, isHerwig, LJ);
0294             std::array<TLegend*, 2>* l = new std::array<TLegend*, 2> {l_head, l_data};
0295             Legs->at(j)=l;
0296         }
0297         RatioLegs->at(i)=Legs;
0298     }
0299     
0300     TH1F* ha1=(TH1F*)herwig_pt->at(0)->Clone();
0301     TH1F* pa1=(TH1F*)pythia_pt->at(0)->Clone();
0302     for(int i = 1; i<(int)herwig_pt->size(); i++) ha1->Add(herwig_pt->at(i));
0303     for(int i = 1; i<(int)pythia_pt->size(); i++) pa1->Add(pythia_pt->at(i));
0304 
0305     TH1F* hj1=(TH1F*)herwig_l_pt->at(0)->Clone();
0306     TH1F* pj1=(TH1F*)pythia_l_pt->at(0)->Clone();
0307     for(int i = 1; i<(int)herwig_l_pt->size(); i++) hj1->Add(herwig_l_pt->at(i));
0308     for(int i = 1; i<(int)pythia_l_pt->size(); i++) pj1->Add(pythia_l_pt->at(i));
0309     TCanvas* cA=new TCanvas("Combined_all", "Combined_all");
0310     TCanvas* cL=new TCanvas("Combined_lead", "Combined_lead");
0311     std::vector<TPad*>* pA = co->AddPads(cA);
0312     std::vector<TPad*>* pL = co->AddPads(cL);
0313     pA->at(1)->cd();
0314     pA->at(1)->SetLogy();
0315     ha1->Draw("pmc plc");
0316     pa1->Draw("pmc plc same");
0317     ha1->SetMarkerStyle(20);
0318     pa1->SetMarkerStyle(24);
0319     TH1F* ra1 = co->GetRatioPlot(ha1, pa1);
0320 
0321     TLegend* l1 = new TLegend(0.7, 0.8, 1,1);
0322     TLegend* l2 = new TLegend(0.7, 0.6, 1, 0.8);
0323     co->SetsPhenixHeaderLegend(l1, "none");
0324     co->SetLegend(l2);
0325     
0326     l1->AddEntry("","Combined Samples all jets","");
0327     l2->AddEntry(ha1, "Herwig", "pl");
0328     l2->AddEntry(pa1, "Pythia", "pl");
0329     l2->AddEntry(ra1, "Herwig / Pythia", "pl");
0330     l1->Draw();
0331     l2->Draw();
0332     pA->at(0)->cd();
0333     ra1->Draw("pmc plc");
0334     ra1->SetYTitle("Herwig / Pythia");
0335     ra1->GetYaxis()->SetRangeUser(0,2);
0336     ra1->SetMarkerStyle(39);
0337     cA->cd();
0338     pA->at(1)->Draw();
0339     pA->at(0)->Draw();
0340 
0341     cL->cd();
0342     pL->at(1)->cd();
0343     pL->at(1)->SetLogy();
0344     hj1->Draw("pmc plc");
0345     pj1->Draw("pmc plc same");
0346     hj1->SetMarkerStyle(20);
0347     pj1->SetMarkerStyle(24);
0348     TH1F* rj1 = co->GetRatioPlot(hj1, pj1);
0349 
0350     TLegend* l3 = new TLegend(0.7, 0.8, 1,1);
0351     TLegend* l4 = new TLegend(0.7, 0.6, 1, 0.8);
0352     co->SetsPhenixHeaderLegend(l3, "none");
0353     co->SetLegend(l4);
0354     
0355     l3->AddEntry("","Combined Samples Lead jets","");
0356     l4->AddEntry(hj1, "Herwig", "pl");
0357     l4->AddEntry(pj1, "Pythia", "pl");
0358     l4->AddEntry(rj1, "Herwig / Pythia", "pl");
0359     l3->Draw();
0360     l4->Draw();
0361     pL->at(0)->cd();
0362     rj1->Draw("pmc plc");
0363     rj1->GetYaxis()->SetRangeUser(0,2);
0364     rj1->SetYTitle("Herwig / Pythia");
0365     rj1->SetMarkerStyle(39);
0366     cL->cd();
0367     pL->at(1)->Draw();
0368     pL->at(0)->Draw();
0369     
0370     std::vector<TH1F*> rHA  {};
0371     std::vector<TH1F*> rHL  {};
0372     std::vector<TH1F*> rPA  {};
0373     std::vector<TH1F*> rPL  {};
0374     SampleRatios(herwig_pt, &rHA ,ha1, 
0375             RatioCanvi.at(0)->at(0), RatioLegs->at(0)->at(0)->at(0) , RatioLegs->at(0)->at(0)->at(1), 
0376             co, "Herwig", cutnames);
0377     SampleRatios(herwig_l_pt, &rHL ,hj1, 
0378             RatioCanvi.at(0)->at(1), RatioLegs->at(0)->at(1)->at(0) , RatioLegs->at(0)->at(1)->at(1),
0379                 co, "Herwig", cutnames);
0380     SampleRatios(pythia_pt, &rPA ,pa1,
0381                 RatioCanvi.at(1)->at(0), RatioLegs->at(1)->at(0)->at(0) , RatioLegs->at(1)->at(0)->at(1),
0382                 co, "Pythia", cutnames);
0383     SampleRatios(pythia_l_pt, &rPL ,pj1, 
0384             RatioCanvi.at(1)->at(1), RatioLegs->at(1)->at(1)->at(0) , RatioLegs->at(0)->at(0)->at(1), 
0385             co, "Pythia", cutnames);
0386     PlotBothThePlots(
0387         herwig_pt, herwig_l_pt,
0388             pythia_pt, pythia_l_pt, 
0389         &All, &Lead,
0390         co, cutnames);
0391     std::vector<TCanvas*> Canvases_to_Print;
0392     for(int i = 0; i<(int) All.size();  i++) Canvases_to_Print.push_back(All.at(i));
0393     for(int i = 0; i<(int) Lead.size(); i++) Canvases_to_Print.push_back(Lead.at(i));
0394     for(int i = 0; i< 2; i++)
0395             for(int j = 0; j < 2; j++)
0396             Canvases_to_Print.push_back(RatioCanvi.at(i)->at(j));
0397     Canvases_to_Print.push_back(cA);
0398     Canvases_to_Print.push_back(cL);
0399 
0400     co->DumpHistos("~/Herwig_comb_pt_QA.pdf", &Canvases_to_Print);
0401     for(int i=0; i<(int)Canvases_to_Print.size(); i++)
0402             Canvases_to_Print.at(i)->SaveSource(std::format("{}.C", Canvases_to_Print.at(i)->GetName()).c_str());
0403 
0404     return;
0405 }
0406             
0407             
0408 void MakeParticleRatios( std::map<std::string, particle_colls*>* h,  std::map<std::string, particle_colls*>* p,  std::map<std::string, particle_colls*>* r)
0409 {
0410     for(auto m:*h) 
0411     {
0412         if(!m.second) continue;
0413         auto pt = conf->GetRatioPlot(h->at(m.first)->pt, p->at(m.first)->pt);
0414         auto n = conf->GetRatioPlot(h->at(m.first)->n, p->at(m.first)->n);
0415         auto phi_eta = conf->GetRatioPlot(h->at(m.first)->phi_eta, p->at(m.first)->phi_eta);
0416         particle_colls* part = new particle_colls;
0417             part->pt    = pt;
0418         part->n     = n;
0419             part->phi_eta   = phi_eta;
0420         r->insert(std::pair<std::string, particle_colls*> (m.first,part));
0421     }
0422     return;
0423 }   
0424 void CollectEventPlots(TDirectory* event, std::vector<TH1F*>* particle_level, std::vector<TH2F*>* particle_correlations,  std::map<std::string, particle_colls*>* particle_coll)
0425 {
0426     //Lots of 2D Plots here as well, so follow along with the PhotonJetApproach
0427     event->cd();
0428     
0429     //get the 1D Plots across all particle types
0430     TH1F* particle_eta  = (TH1F*) event->Get("h_particle_eta");
0431     TH1F* particle_phi  = (TH1F*) event->Get("h_particle_phi");
0432     TH1F* particle_e    = (TH1F*) event->Get("h_particle_e");
0433     TH1F* particle_et   = (TH1F*) event->Get("h_particle_et");
0434     TH1F* particle_pt   = (TH1F*) event->Get("h_particle_pt");
0435     TH1F* total_E       = (TH1F*) event->Get("h_total_E");
0436     particle_eta->Scale(1/(float)total_E->GetEntries());
0437     particle_phi->Scale(1/(float)total_E->GetEntries());
0438     particle_eta->SetYTitle(" N_{hits} / N_{evts}");
0439     particle_phi->SetYTitle(" N_{hits} / N_{evts}");
0440     particle_level->push_back(particle_eta);
0441     particle_level->push_back(particle_phi);
0442     particle_level->push_back(particle_e);
0443     particle_level->push_back(particle_et);
0444     particle_level->push_back(particle_pt);
0445     particle_level->push_back(total_E);
0446 
0447     //now get the 2D Plots  
0448     TH2F* h_particle_et_eta     = (TH2F*) event->Get("h_particle_et_eta");
0449     TH2F* h_particle_et_phi     = (TH2F*) event->Get("h_particle_et_phi");
0450     TH2F* h_particle_pt_eta     = (TH2F*) event->Get("h_particle_pt_eta");
0451     TH2F* h_particle_pt_phi     = (TH2F*) event->Get("h_particle_pt_phi");
0452     TH2F* h_particle_e_eta      = (TH2F*) event->Get("h_particle_e_eta");
0453     TH2F* h_particle_e_phi      = (TH2F*) event->Get("h_particle_e_phi");
0454     TH2F* h_particle_phi_eta    = (TH2F*) event->Get("h_particle_phi_eta");
0455 
0456     particle_correlations->push_back(h_particle_et_eta);
0457     particle_correlations->push_back(h_particle_et_phi);
0458     particle_correlations->push_back(h_particle_pt_eta);
0459     particle_correlations->push_back(h_particle_pt_phi);
0460     particle_correlations->push_back(h_particle_e_eta);
0461     particle_correlations->push_back(h_particle_e_phi);
0462     particle_correlations->push_back(h_particle_phi_eta);
0463 
0464     //Individual particle type comp
0465     TDirectory* particles = (TDirectory*) event->GetDirectory("ParticleCat");
0466     particles->cd();
0467     std::vector<std::string> particle_names {"electron", "proton", "neutron", "pion", "photon"};
0468     for(auto n:particle_names)
0469     {
0470         TH1F* h_p_pt    = (TH1F*) particles->Get(std::format("h_{}_pt", n).c_str());
0471         TH1I* h_p_n = (TH1I*) particles->Get(std::format("h_{}_n", n).c_str());
0472         TH2F* h_p_ph    = (TH2F*) particles->Get(std::format("h_{}_phi_eta", n).c_str());
0473         if(!h_p_pt) continue;
0474         particle_colls* part = new particle_colls;
0475             part->pt    = h_p_pt;
0476             part->n     = h_p_n;
0477             part->phi_eta   = h_p_ph;
0478         particle_coll->insert(std::pair<std::string, particle_colls*>(n, part));
0479     }
0480     return;
0481 }
0482 void PlotEventObs(TH1F* herwig_obs, TH1F* pythia_obs, TCanvas* c1)
0483 {
0484     //plotting the kinematic variables
0485     std::string var = herwig_obs->GetXaxis()->GetTitle();
0486     std::vector<TPad*>* pads = conf->AddPads(c1);
0487     pads->at(1)->cd();
0488     TLegend* l_data=new TLegend(0.7, 0.4, 1, 0.58);
0489     TLegend* l_header=new TLegend(0.5, 0.6, 0.7, 1);
0490     conf->SetsPhenixHeaderLegend(l_header, tag);
0491     l_header->AddEntry("", std::format("Event {}", var).c_str(), "");
0492     conf->SetLegend(l_data);
0493     if(var.find("[GeV]") != std::string::npos) pads->at(1)->SetLogy();
0494     auto ratio = conf->GetRatioPlot(herwig_obs, pythia_obs);
0495     pythia_obs->SetMarkerStyle(24);
0496     herwig_obs->SetMarkerStyle(20);
0497     THStack* s=new THStack("s", "s");
0498     s->Add(herwig_obs);
0499     s->Add(pythia_obs);
0500     s->Draw("axis");
0501     s->Draw("plc pmc nostack same e1");
0502     s->GetXaxis()->SetTitle(var.c_str());
0503     s->GetYaxis()->SetTitle(herwig_obs->GetYaxis()->GetTitle());
0504     s->GetYaxis()->SetRangeUser(herwig_obs->GetMinimum()*0.33, 3*(herwig_obs->GetMaximum()));
0505     l_data->AddEntry(herwig_obs);
0506     l_data->AddEntry(pythia_obs);
0507     l_header->Draw();
0508     l_data->Draw("same");
0509     pads->at(0)->cd();
0510     ratio->SetMarkerStyle(34);
0511     ratio->Draw("pmc plc e1");
0512     l_data->AddEntry(ratio, "Herwig / Pythia");
0513     c1->cd();
0514     pads->at(0)->Draw();
0515     pads->at(1)->Draw();
0516     return;
0517 }
0518 void PlotEventObs(TH2F* herwig_obs, TH2F* pythia_obs, TCanvas* c1)
0519 {
0520     //plotting the kinematic variables
0521     std::string var = herwig_obs->GetXaxis()->GetTitle();
0522     std::string var2 = herwig_obs->GetYaxis()->GetTitle();
0523     std::vector<TPad*>* pads = new std::vector<TPad*>();
0524     TPad* p1=new TPad("p1", "p1", 0, 0.52, 0.48, 1);
0525     TPad* p2=new TPad("p2", "p2", 0.52, 0.52, 1, 1);
0526     TPad* p3=new TPad("p3", "p3", 0, 0, 0.48, 0.48);
0527     TPad* p4=new TPad("p4", "p4", 0.52, 0, 1, 0.48);
0528     pads->push_back(p1);
0529     pads->push_back(p2);
0530     pads->push_back(p3);
0531     pads->push_back(p4);
0532     pads->at(0)->cd();
0533     TLegend* l_header=new TLegend(0.1, 0.1, 1, 1);
0534     conf->SetsPhenixHeaderLegend(l_header, tag);
0535     l_header->AddEntry("", std::format("Event {} - {}", var, var2).c_str(), "");
0536     if(var.find("[GeV]") != std::string::npos) pads->at(1)->SetLogx();
0537     if(var2.find("[GeV]") != std::string::npos) pads->at(1)->SetLogy();
0538     pads->at(1)->SetLogz();
0539     if(var.find("[GeV]") == std::string::npos &&  var2.find("[GeV]") == std::string::npos)
0540     {
0541         pads->at(1)->SetLogz(0);
0542         gStyle->SetPalette(100, style_points->Enby_log_gradient_PT);    
0543     }
0544     else{
0545         gStyle->SetPalette(100, style_points->Enby_gradient_PT);    
0546     }
0547     auto ratio = conf->GetRatioPlot(herwig_obs, pythia_obs);
0548     l_header->Draw();
0549     pads->at(1)->cd();
0550     herwig_obs->Draw("colz");
0551     TLegend* l_data=new TLegend(0.5, 0.8, 1, 1);
0552     conf->SetLegend(l_data);
0553     l_data->SetTextSize(0.05f);
0554     l_data->AddEntry(herwig_obs, std::format("Herwig {}", var).c_str(), "");
0555     l_data->Draw();
0556     pads->at(2)->cd();
0557     if(var.find("[GeV]") != std::string::npos) pads->at(2)->SetLogx();
0558     if(var2.find("[GeV]") != std::string::npos) pads->at(2)->SetLogy();
0559     pads->at(2)->SetLogz();
0560     if(var.find("[GeV]") == std::string::npos &&  var2.find("[GeV]") == std::string::npos) pads->at(2)->SetLogz(0);
0561     pythia_obs->Draw("colz");
0562     TLegend* l_data_p=new TLegend(0.5, 0.8, 1, 1);
0563     conf->SetLegend(l_data_p);
0564     l_data_p->SetTextSize(0.05f);
0565     l_data_p->AddEntry(pythia_obs, std::format("Pythia {}", var).c_str(), "");
0566     l_data_p->Draw();
0567     
0568     pads->at(3)->cd();
0569     if(var.find("[GeV]") != std::string::npos) pads->at(3)->SetLogx();
0570     if(var2.find("[GeV]") != std::string::npos) pads->at(3)->SetLogy();
0571     pads->at(3)->SetLogz();
0572     if(var.find("[GeV]") == std::string::npos &&  var2.find("[GeV]") == std::string::npos) pads->at(3)->SetLogz(0);
0573     ratio->SetMarkerStyle(34);
0574     ratio->Draw("colz");
0575     TLegend* l_data_r=new TLegend(0.5, 0.8, 1, 1);
0576     conf->SetLegend(l_data_r);
0577     l_data_r->SetTextSize(0.05f);
0578     l_data_r->AddEntry(pythia_obs, std::format(" Herwig / Pythia {}", var).c_str(), "");
0579     l_data_r->Draw();
0580     c1->cd();
0581     pads->at(0)->Draw();
0582     pads->at(1)->Draw();
0583     pads->at(2)->Draw();
0584     pads->at(3)->Draw();
0585     return;
0586 }
0587 void PlotEventObs(TH1F* herwig_particle, TH1F* pythia_particle, std::map<std::string, particle_colls*>* h_pc, std::map<std::string, particle_colls*>* p_pc, std::map<std::string, particle_colls*>* r_pc, TCanvas* c1)
0588 {
0589     //comparing the particle species pt
0590     std::vector<TPad*>* pads=conf->AddPads(c1);
0591     pads->at(1)->cd();
0592     pads->at(1)->SetLogy();
0593     TLegend* l_data=new TLegend(0.3, 0.4, 1, 0.58);
0594     TLegend* l_header=new TLegend(0.5, 0.6, 0.7, 1);
0595     conf->SetsPhenixHeaderLegend(l_header, tag);
0596     l_header->AddEntry("", "Particle Species p_{T}", "");
0597     conf->SetLegend(l_data);
0598     l_data->SetNColumns(3);
0599     THStack* herwig = new THStack("herwig","herwig");
0600     THStack* pythia = new THStack("pythia","pythia");
0601     THStack* ratio = new THStack("ratio","ratio");
0602     herwig_particle->SetMarkerStyle(20);
0603     pythia_particle->SetMarkerStyle(24);
0604     herwig->Add(herwig_particle);
0605     pythia->Add(pythia_particle);
0606     pads->at(0)->cd();
0607     TH1F* ratio_particle=conf->GetRatioPlot(herwig_particle, pythia_particle);
0608     ratio_particle->SetMarkerStyle(37);
0609     ratio->Add(ratio_particle);
0610     l_data->AddEntry(herwig_particle, "Herwig All Final State Particles");
0611     l_data->AddEntry(pythia_particle, "Pythia All Final State Particles");
0612     l_data->AddEntry(ratio_particle, "Herwig / Pythia Particles");
0613     int i=1;
0614     for(auto m: *h_pc)
0615     {
0616         h_pc->at(m.first)->pt->SetMarkerStyle(i+20);
0617         p_pc->at(m.first)->pt->SetMarkerStyle(i+24);
0618         r_pc->at(m.first)->pt->SetMarkerStyle(2*i+39);
0619         if(i >= 4){
0620              m.second->pt->SetMarkerStyle(2*(i-4)+34);
0621              r_pc->at(m.first)->pt->SetMarkerStyle(2*i+39);
0622         }
0623         herwig->Add(h_pc->at(m.first)->pt);
0624         pythia->Add(p_pc->at(m.first)->pt);
0625         ratio->Add(r_pc->at(m.first)->pt);
0626         l_data->AddEntry( m.second->pt, std::format("Herwig Final State {}", m.first).c_str());
0627         l_data->AddEntry( p_pc->at(m.first)->pt, std::format("Pythia Final State {}", m.first).c_str());
0628         l_data->AddEntry( r_pc->at(m.first)->pt, std::format("Herwig / Pythia {}", m.first).c_str());
0629     }
0630     c1->cd();
0631     pads->at(1)->cd();
0632     bool herwig_high = false;
0633     if (herwig->GetMaximum() > pythia->GetMaximum() ) herwig_high=true;
0634     if (herwig_high) herwig->Draw("axis");
0635     else pythia->Draw("axis");
0636     herwig->Draw("pmc plc nostack same e1");
0637     pythia->Draw("same pmc plc nostack e1");
0638     herwig->GetXaxis()->SetTitle(herwig_particle->GetXaxis()->GetTitle());
0639     herwig->GetYaxis()->SetTitle(herwig_particle->GetYaxis()->GetTitle());
0640     herwig->GetYaxis()->SetRangeUser(herwig_particle->GetMinimum()*0.33, herwig_particle->GetMaximum()*3);
0641     l_data->Draw();
0642     l_header->Draw();
0643     pads->at(0)->cd();
0644     ratio->Draw("pmc plc nostack e1");
0645     ratio->GetYaxis()->SetRangeUser(0.01, 2.);
0646     ratio->GetYaxis()->SetTitle(" Herwig / Pythia" ); 
0647     ratio->Draw("pmc plc nostack e1");
0648     c1->cd();
0649     pads->at(0)->Draw();
0650     pads->at(1)->Draw();
0651     return; 
0652 }
0653 void PlotEventObs(TH1I* herwig_particle, TH1I* pythia_particle, std::map<std::string, particle_colls*>* h_pc, std::map<std::string, particle_colls*>* p_pc, std::map<std::string, particle_colls*>* r_pc, TCanvas* c1)
0654 {
0655     //comparing the particle species n 
0656     std::vector<TPad*>* pads=conf->AddPads(c1);
0657     pads->at(1)->cd();
0658     pads->at(1)->SetLogy();
0659     TLegend* l_data=new TLegend(0.2, 0.4, 0.6, 0.58);
0660     TLegend* l_header=new TLegend(0.5, 0.6, 0.7, 1);
0661     conf->SetsPhenixHeaderLegend(l_header, tag);
0662     l_header->AddEntry("", "Particle Species Count", "");
0663     conf->SetLegend(l_data);
0664     l_data->SetNColumns(3);
0665     THStack* herwig = new THStack("herwig","herwig");
0666     THStack* pythia = new THStack("pythia","pythia");
0667     THStack* ratio = new THStack("ratio","ratio");
0668     herwig_particle->SetMarkerStyle(20);
0669     pythia_particle->SetMarkerStyle(24);
0670     herwig_particle->GetXaxis()->SetRangeUser(0,200);
0671     herwig->Add(herwig_particle);
0672     pythia->Add(pythia_particle);
0673     pads->at(0)->cd();
0674     TH1I* ratio_particle=conf->GetRatioPlot(herwig_particle, pythia_particle);
0675     ratio_particle->SetMarkerStyle(37);
0676     ratio->Add(ratio_particle);
0677     l_data->AddEntry(herwig_particle, "Herwig All Final State Particles");
0678     l_data->AddEntry(pythia_particle, "Pythia All Final State Particles");
0679     l_data->AddEntry(ratio_particle, "Herwig / Pythia Particles");
0680     int i=1;
0681     for(auto m: *h_pc)
0682     {
0683         h_pc->at(m.first)->n->SetMarkerStyle(i+20);
0684         p_pc->at(m.first)->n->SetMarkerStyle(i+24);
0685         r_pc->at(m.first)->n->SetMarkerStyle(2*i+39);
0686         if(i >= 4){
0687              h_pc->at(m.first)->n->SetMarkerStyle(2*(i-4)+34);
0688              r_pc->at(m.first)->n->SetMarkerStyle(2*i+39);
0689         }
0690         pads->at(1)->cd();
0691         pads->at(0)->cd();
0692         herwig->Add(h_pc->at(m.first)->n);
0693         pythia->Add(p_pc->at(m.first)->n);
0694         ratio->Add(r_pc->at(m.first)->n);
0695         l_data->AddEntry( h_pc->at(m.first)->n, std::format("Herwig Final State {}", m.first).c_str());
0696         l_data->AddEntry( p_pc->at(m.first)->n, std::format("Pythia Final State {}", m.first).c_str());
0697         l_data->AddEntry( r_pc->at(m.first)->n, std::format("Herwig / Pythia {}", m.first).c_str());
0698     }
0699     c1->cd();
0700     pads->at(1)->cd();
0701     bool herwig_high = false;
0702     if (herwig->GetMaximum() > pythia->GetMaximum() ) herwig_high=true;
0703     if (herwig_high) herwig->Draw("axis");
0704     else pythia->Draw("axis");
0705     herwig->Draw("pmc plc nostack same e1");
0706     pythia->Draw("same pmc plc nostack e1");
0707     herwig->GetXaxis()->SetRangeUser(0, 200);
0708     pythia->GetXaxis()->SetRangeUser(0, 200);
0709     herwig->GetXaxis()->SetTitle(herwig_particle->GetXaxis()->GetTitle());
0710     l_data->Draw();
0711     l_header->Draw();
0712     pads->at(0)->cd();
0713     ratio->Draw("pmc plc nostack e1");
0714     ratio->GetXaxis()->SetTitle(herwig_particle->GetXaxis()->GetTitle());
0715     ratio->GetYaxis()->SetTitle(" Herwig / Pythia" ); 
0716     ratio->GetXaxis()->SetRangeUser(0, 200);
0717     ratio->GetYaxis()->SetRangeUser(0.01, 2.);
0718     ratio->Draw("pmc plc nostack e1");  
0719     c1->cd();
0720     pads->at(0)->Draw();
0721     pads->at(1)->Draw();
0722     return; 
0723 }
0724 void PlotEventObs(TH2F* particle, std::map<std::string, particle_colls*>* pc, TCanvas* c1, std::string generator)
0725 {
0726     //comparing the particle species phi-eta hits
0727     std::vector<TPad*>* pads = conf->Canvas2DDivide(c1);
0728     pads->at(0)->cd();
0729     TLegend* l_header=new TLegend(0.5, 0.6, 1, 1);
0730     conf->SetsPhenixHeaderLegend(l_header, tag);
0731     l_header->SetTextSize(0.03f);
0732     l_header->AddEntry("", generator.c_str(), "");
0733     l_header->AddEntry("", "Final State #eta - #varphi hit distribution", "");
0734     particle->Draw("colz");
0735     l_header->Draw();
0736     int i=1;
0737     for(auto m:*pc)
0738     {
0739         pads->at(i)->cd();
0740         TLegend* l_data=new TLegend(0.5, 0.8, 1, 0.9);
0741         conf->SetLegend(l_data);
0742         l_data->AddEntry("", m.first.c_str(), "");
0743         m.second->phi_eta->Draw("colz");
0744         l_data->Draw("same");
0745         i++;
0746     }
0747     c1->cd();
0748     for(int i=0; i<(int)pads->size(); i++) pads->at(i)->Draw();
0749     return;
0750 }
0751 
0752 
0753 void PlotEventPlots(TDirectory* herwig_event, TDirectory* pythia_event, std::vector<TCanvas*>* Canvi)
0754 {
0755     std::vector<TH1F*>* herwig_particle_level   = new std::vector<TH1F*> ();
0756     std::vector<TH2F*>* herwig_particle_correlations= new std::vector<TH2F*> ();
0757     TH1I* herwig_particle_n             = (TH1I*) herwig_event->Get("h_particle_n");
0758     std::map<std::string, particle_colls*>* herwig_particle_coll = new std::map<std::string, particle_colls*> ();
0759     std::vector<TH1F*>* pythia_particle_level   = new std::vector<TH1F*> ();
0760     std::vector<TH2F*>* pythia_particle_correlations= new std::vector<TH2F*> ();
0761     TH1I* pythia_particle_n             = (TH1I*) pythia_event->Get("h_particle_n");
0762     std::map<std::string, particle_colls*>* pythia_particle_coll  = new std::map<std::string, particle_colls*> ();
0763     CollectEventPlots(herwig_event, herwig_particle_level, herwig_particle_correlations, herwig_particle_coll);
0764     CollectEventPlots(pythia_event, pythia_particle_level, pythia_particle_correlations, pythia_particle_coll);
0765     std::map<std::string, particle_colls*>* ratio_particle_coll = new std::map<std::string, particle_colls*> ();
0766     MakeParticleRatios( herwig_particle_coll, pythia_particle_coll, ratio_particle_coll);
0767     gStyle->SetPalette(100, style_points->Lesbian_gradient_PT); 
0768     for(int i=0; i<(int)herwig_particle_level->size(); i++){
0769         std::string var = herwig_particle_level->at(i)->GetXaxis()->GetTitle();
0770         TCanvas* c1=new TCanvas(std::format("Canv_{}", var).c_str(),std::format("Canv_{}", var).c_str()); 
0771         PlotEventObs(herwig_particle_level->at(i), pythia_particle_level->at(i), c1);
0772         Canvi->push_back(c1);
0773         if(var.find("p_{T}") != std::string::npos){
0774             TCanvas* c2 = new TCanvas("pt", "pt");
0775             PlotEventObs(herwig_particle_level->at(i), pythia_particle_level->at(i), herwig_particle_coll, pythia_particle_coll, ratio_particle_coll, c2);
0776             Canvi->push_back(c2);
0777         }
0778     }
0779     for(int i=0; i<(int)herwig_particle_correlations->size(); i++){
0780         std::string var = herwig_particle_correlations->at(i)->GetXaxis()->GetTitle();
0781         std::string var2 = herwig_particle_correlations->at(i)->GetYaxis()->GetTitle();
0782         TCanvas* c1=new TCanvas(std::format("Canv2d_{}_{}", var, var2).c_str(),std::format("Canv2d_{}_{}", var, var2).c_str()); 
0783         PlotEventObs(herwig_particle_correlations->at(i), pythia_particle_correlations->at(i), c1);
0784         Canvi->push_back(c1);
0785         if(var.find("eta") !=std::string::npos && var2.find("varphi") != std::string::npos) 
0786         {
0787             TCanvas* ch = new TCanvas("etph_h", "etph_h");
0788             TCanvas* cp = new TCanvas("etph_p", "etph_p");
0789             TCanvas* cr = new TCanvas("etph_r", "etph_r");
0790             TH2F* ratio = conf->GetRatioPlot(herwig_particle_correlations->at(i), pythia_particle_correlations->at(i));
0791             PlotEventObs(herwig_particle_correlations->at(i), herwig_particle_coll, ch, "Herwig7.2");
0792             PlotEventObs(pythia_particle_correlations->at(i), pythia_particle_coll, cp, "Pythia8");
0793             PlotEventObs(ratio, ratio_particle_coll, cr, "Herwig / Pythia");
0794         }
0795     }
0796     gStyle->SetPalette(100, style_points->Trans_gradient_PT);   
0797     TCanvas* c3 = new TCanvas("n", "n");
0798     PlotEventObs(herwig_particle_n, pythia_particle_n, herwig_particle_coll, pythia_particle_coll, ratio_particle_coll, c3);
0799     Canvi->push_back(c3);
0800     return;
0801 
0802 }
0803 void PlotSingleGenPhotonJetObs(std::vector<TH2F*>* obs, TCanvas* c1, std::vector<TPad*>* pads, std::string label)
0804 {
0805     c1->cd();
0806     pads->at(0)->cd();
0807     TLegend* l_label=new TLegend(0.1, 0.1, 1, 1);
0808     conf->SetsPhenixHeaderLegend(l_label, tag);
0809     l_label->AddEntry("", label.c_str(), "");
0810     l_label->AddEntry("", std::format("Jet {} Versus Photon {}", obs->at(0)->GetXaxis()->GetTitle(), obs->at(0)->GetYaxis()->GetTitle()).c_str(), "");
0811     for(int i=0; i<(int)obs->size(); i++)
0812     {
0813         pads->at(i+1)->cd();
0814         TLegend* l_data=new TLegend(0.7, 0.9, 0.9, 1);
0815         conf->SetLegend(l_data);
0816         l_data->AddEntry(obs->at(i), obs->at(i)->GetTitle(), "");
0817         obs->at(i)->Draw("colz");
0818         l_data->Draw("same");
0819         c1->cd();
0820         pads->at(i+1)->Draw();
0821     }
0822     return;
0823 }
0824 void PlotPhotonJetObs(std::vector<TH2F*>* herwig_obs, std::vector<TH2F*>* pythia_obs, TCanvas* herwig_Canvas, TCanvas* pythia_Canvas, TCanvas* ratio_Canvas)
0825 {
0826     std::vector<TPad*>* subdivided_herwig = conf->Canvas2DDivide(herwig_Canvas);
0827     std::vector<TPad*>* subdivided_pythia = conf->Canvas2DDivide(pythia_Canvas);
0828     std::vector<TPad*>* subdivided_ratio = conf->Canvas2DDivide(ratio_Canvas);
0829     std::vector<TH2F*>* ratio_obs = conf->GetRatioPlots(herwig_obs, pythia_obs);
0830     PlotSingleGenPhotonJetObs(herwig_obs, herwig_Canvas, subdivided_herwig, "Herwig7.2");
0831     PlotSingleGenPhotonJetObs(pythia_obs, pythia_Canvas, subdivided_pythia, "Pythia8");
0832     PlotSingleGenPhotonJetObs(ratio_obs, ratio_Canvas, subdivided_ratio, "Herwig / Pythia");
0833     return;
0834 }
0835 void CollectJets(TDirectory* jets, std::array<std::vector<TH1F*>*, 4>* all_jets, std::array<std::vector<TH1F*>*, 4>* lead_jets, std::vector<TH1I*>* all_jet_n, std::vector<TH1I*>* all_jet_comp, std::vector<TH1I*>* lead_jet_comp)
0836 {
0837     for(int i=2; i<7; i++)
0838     {
0839         jets->cd();
0840         TDirectory* r_dir=(TDirectory*)jets->GetDirectory(std::format("R0{}Jets", i).c_str());
0841         r_dir->cd();
0842         all_jets->at(0)->push_back((TH1F*)r_dir->Get(std::format("h_jet_r0{}_pt", i).c_str()));
0843         all_jets->at(1)->push_back((TH1F*)r_dir->Get(std::format("h_jet_r0{}_e", i).c_str()));
0844         all_jets->at(2)->push_back((TH1F*)r_dir->Get(std::format("h_jet_r0{}_phi", i).c_str()));
0845         all_jets->at(3)->push_back((TH1F*)r_dir->Get(std::format("h_jet_r0{}_eta", i).c_str()));
0846         all_jet_n->push_back((TH1I*)r_dir->Get(std::format("h_jet_r0{}_n",i).c_str()));
0847         all_jet_comp->push_back((TH1I*)r_dir->Get(std::format("h_jet_r0{}_comp",i).c_str()));
0848         TDirectory* leaddir=(TDirectory*)r_dir->GetDirectory("Lead_Jets");
0849         lead_jets->at(0)->push_back((TH1F*)leaddir->Get(std::format("h_lead_jet_r0{}_pt", i).c_str()));
0850         lead_jets->at(1)->push_back((TH1F*)leaddir->Get(std::format("h_lead_jet_r0{}_e", i).c_str()));
0851         lead_jets->at(2)->push_back((TH1F*)leaddir->Get(std::format("h_lead_jet_r0{}_phi", i).c_str()));
0852         lead_jets->at(3)->push_back((TH1F*)leaddir->Get(std::format("h_lead_jet_r0{}_eta", i).c_str()));
0853         lead_jet_comp->push_back((TH1I*)leaddir->Get(std::format("h_lead_jet_r0{}_comp",i).c_str()));
0854         r_dir->cd();
0855         float n_evts = (float)lead_jets->at(0)->back()->GetEntries();
0856         all_jets->at(2)->back()->Scale(1/n_evts);
0857         all_jets->at(2)->back()->SetYTitle(" N_{jets} / N_{evts} ");
0858         all_jets->at(3)->back()->Scale(1/n_evts);
0859         all_jets->at(3)->back()->SetYTitle(" N_{jets} / N_{evts} ");
0860         lead_jets->at(2)->back()->Scale(1/n_evts);
0861         lead_jets->at(2)->back()->SetYTitle(" N_{jets} / N_{evts} ");
0862         lead_jets->at(3)->back()->Scale(1/n_evts);
0863         lead_jets->at(3)->back()->SetYTitle(" N_{jets} / N_{evts} ");
0864     }   
0865     return;
0866 
0867 }
0868 void CollectPhotons(TDirectory* photons, std::array<TH1F*, 4>* all_photons, std::array<TH1F*, 4>* lead_photons, std::array<TH1F*, 4>* direct_photons, std::array<TH1F*, 4>* frag_photons, TH1I* all_photons_n, TH1I* direct_photons_n, TH1I* frag_photons_n)
0869 {
0870     photons->cd();
0871     all_photons->at(0)=(TH1F*)photons->Get("h_photons_pt");
0872     all_photons->at(1)=(TH1F*)photons->Get("h_photons_e");
0873     all_photons->at(2)=(TH1F*)photons->Get("h_photons_phi");
0874     all_photons->at(3)=(TH1F*)photons->Get("h_photons_eta");
0875     all_photons_n=(TH1I*)photons->Get("h_photon_n");
0876     TDirectory* leaddir=(TDirectory*)photons->GetDirectory("Lead_Photons");
0877     lead_photons->at(0)=(TH1F*)leaddir->Get("h_lead_photons_pt");
0878     lead_photons->at(1)=(TH1F*)leaddir->Get("h_lead_photons_e");
0879     lead_photons->at(2)=(TH1F*)leaddir->Get("h_lead_photons_phi");
0880     lead_photons->at(3)=(TH1F*)leaddir->Get("h_lead_photons_eta");
0881     photons->cd();
0882     TDirectory* directdir=(TDirectory*)photons->GetDirectory("Direct_Photons");
0883     direct_photons->at(0)=(TH1F*)directdir->Get("h_direct_photons_pt");
0884     direct_photons->at(1)=(TH1F*)directdir->Get("h_direct_photons_e");
0885     direct_photons->at(2)=(TH1F*)directdir->Get("h_direct_photons_phi");
0886     direct_photons->at(3)=(TH1F*)directdir->Get("h_direct_photons_eta");
0887     direct_photons_n=(TH1I*)directdir->Get("h_direct_photon_n");
0888     photons->cd();
0889     TDirectory* fragdir=(TDirectory*)photons->GetDirectory("Fragmentation_Photons");
0890     frag_photons->at(0)=(TH1F*)fragdir->Get("h_frag_photons_pt");
0891     frag_photons->at(1)=(TH1F*)fragdir->Get("h_frag_photons_e");
0892     frag_photons->at(2)=(TH1F*)fragdir->Get("h_frag_photons_phi");
0893     frag_photons->at(3)=(TH1F*)fragdir->Get("h_frag_photons_eta");
0894     frag_photons_n=(TH1I*)fragdir->Get("h_frag_photon_n");
0895     photons->cd();
0896     float n_evts = (float) lead_photons->at(0)->GetEntries();
0897     //make phi and eta dist an average across events
0898     all_photons->at(2)->Scale(1/n_evts);
0899     all_photons->at(2)->SetYTitle(" N_{#gamma} / N_{evts} ");
0900     all_photons->at(3)->Scale(1/n_evts);
0901     all_photons->at(3)->SetYTitle(" N_{#gamma} / N_{evts} ");
0902     lead_photons->at(2)->Scale(1/n_evts);
0903     lead_photons->at(2)->SetYTitle(" N_{#gamma} / N_{evts} ");
0904     lead_photons->at(3)->Scale(1/n_evts);
0905     lead_photons->at(3)->SetYTitle(" N_{#gamma} / N_{evts} ");
0906     direct_photons->at(2)->Scale(1/n_evts);
0907     direct_photons->at(2)->SetYTitle(" N_{#gamma} / N_{evts} ");
0908     direct_photons->at(3)->Scale(1/n_evts);
0909     direct_photons->at(3)->SetYTitle(" N_{#gamma} / N_{evts} ");
0910     frag_photons->at(2)->Scale(1/n_evts);
0911     frag_photons->at(2)->SetYTitle(" N_{#gamma} / N_{evts} ");
0912     frag_photons->at(3)->Scale(1/n_evts);
0913     frag_photons->at(3)->SetYTitle(" N_{#gamma} / N_{evts} ");
0914     return;
0915 
0916 }
0917 void PlotJetObs(std::vector<TH1F*>* herwig_obs, std::vector<TH1F*>* pythia_obs, TCanvas* obs_canv)
0918 {
0919     //just does the plotting of one variable at a time
0920     if(!herwig_obs || !pythia_obs) return;
0921     TLegend* l_data=new TLegend(0.7, 0.4, 1, 0.58);
0922     TLegend* l_header=new TLegend(0.5, 0.6, 0.7, 1);
0923     conf->SetsPhenixHeaderLegend(l_header, tag);
0924     l_header->AddEntry("", std::format("Jet {}", herwig_obs->at(0)->GetXaxis()->GetTitle()).c_str(), "");
0925     conf->SetLegend(l_data);
0926     l_data->SetNColumns(3);
0927     std::vector<TPad*>* pads=conf->AddPads(obs_canv);
0928     std::string var=herwig_obs->at(0)->GetXaxis()->GetTitle();
0929     bool ispt = false;
0930     if(var.find("p_{T}") !=std::string::npos) ispt=true;
0931     if(ispt)
0932     {
0933         conf->ScaleXS(herwig_obs, true);
0934         conf->ScaleXS(pythia_obs, false);
0935         pads->at(1)->SetLogy();
0936     }
0937     else if(var.find("[GeV]") != std::string::npos) pads->at(1)->SetLogy();
0938     pads->at(1)->cd();
0939     THStack* herwig = new THStack("herwig", "herwig");
0940     THStack* pythia = new THStack("herwig", "herwig");
0941     for(int i= 0 ; i <(int)herwig_obs->size()-1;  i++)
0942     {
0943         pads->at(1)->cd();
0944         herwig_obs->at(i)->SetMarkerStyle(i+20);
0945         herwig_obs->at(i)->SetMarkerSize(2);
0946         if(i==5) herwig_obs->at(i)->SetMarkerStyle(33);
0947         else if(i==6) herwig_obs->at(i)->SetMarkerStyle(34);
0948         herwig->Add(herwig_obs->at(i));
0949         pythia_obs->at(i)->SetLineColor(herwig_obs->at(i)->GetLineColor());
0950         pythia_obs->at(i)->SetMarkerColor(herwig_obs->at(i)->GetMarkerColor());
0951         pythia_obs->at(i)->SetMarkerStyle(i+24);
0952         pythia_obs->at(i)->SetMarkerSize(2);
0953         pythia->Add(pythia_obs->at(i));
0954         l_data->AddEntry(herwig_obs->at(i), std::format("Herwig {}", herwig_obs->at(i)->GetTitle()).c_str());
0955         l_data->AddEntry(pythia_obs->at(i), std::format("Pythia ", pythia_obs->at(i)->GetTitle()).c_str());
0956         pads->at(0)->cd();
0957         TH1F* h_ratio = conf->GetRatioPlot(herwig_obs->at(i), pythia_obs->at(i));
0958         h_ratio->SetMarkerStyle(2*i+37);
0959         h_ratio->Draw("plc pmc same e1");
0960         l_data->AddEntry(h_ratio);
0961     }
0962     pads->at(1)->cd();
0963     bool herwig_high = false;
0964     if (herwig->GetMaximum() > pythia->GetMaximum() ) herwig_high=true;
0965     if (herwig_high) herwig->Draw("axis nostack");
0966     else pythia->Draw("axis nostack");
0967     herwig->Draw("plc pmc nostack same e1");
0968     pythia->Draw("plc pmc same nostack e1");
0969     pythia->GetXaxis()->SetTitle(herwig_obs->at(0)->GetXaxis()->GetTitle());
0970     pythia->GetYaxis()->SetTitle(herwig_obs->at(0)->GetYaxis()->GetTitle());
0971     herwig->GetYaxis()->SetRangeUser(herwig_obs->at(0)->GetMinimum()*0.33, herwig_obs->at(0)->GetMaximum()*3);
0972     l_data->Draw();
0973     l_header->Draw();
0974     obs_canv->cd();
0975     pads->at(0)->Draw();
0976     pads->at(1)->Draw();
0977     return;
0978 }
0979 void PlotJetObs(std::vector<TH1I*>* herwig_obs, std::vector<TH1I*>* pythia_obs, TCanvas* obs_canv)
0980 {
0981     //just does the plotting of one variable at a time
0982     if(!herwig_obs || !pythia_obs) return;
0983     TLegend* l_data=new TLegend(0.6, 0.3, 1, 0.5);
0984     TLegend* l_header=new TLegend(0.8, 0.6, 1, 1);
0985     conf->SetsPhenixHeaderLegend(l_header, tag);
0986     conf->SetLegend(l_data);
0987     l_data->SetNColumns(3);
0988     std::vector<TPad*>* pads=conf->AddPads(obs_canv);
0989     conf->ScaleXS(herwig_obs, true);
0990     conf->ScaleXS(pythia_obs, false);
0991     THStack* herwig = new THStack("herwig", "herwig");
0992     THStack* pythia = new THStack("herwig", "herwig");
0993     for(int i=0; i<(int)herwig_obs->size(); i++)
0994     {
0995         pads->at(1)->cd();
0996         herwig_obs->at(i)->SetMarkerStyle(i+20);
0997         herwig_obs->at(i)->SetMarkerSize(2);
0998         if(i==5) herwig_obs->at(i)->SetMarkerStyle(33);
0999         else if(i==6) herwig_obs->at(i)->SetMarkerStyle(34);
1000         herwig->Add(herwig_obs->at(i));
1001         pythia_obs->at(i)->SetMarkerStyle(i+24);
1002         pythia_obs->at(i)->SetMarkerSize(2);
1003         pythia->Add(pythia_obs->at(i));
1004         l_data->AddEntry(herwig_obs->at(i), std::format("Herwig {}", herwig_obs->at(i)->GetTitle()).c_str());
1005         l_data->AddEntry(pythia_obs->at(i), std::format("Pythia {}", pythia_obs->at(i)->GetTitle()).c_str());
1006         pads->at(0)->cd();
1007         TH1I* h_ratio = conf->GetRatioPlot(herwig_obs->at(i), pythia_obs->at(i));
1008         h_ratio->SetMarkerStyle(2*(i)+37);
1009         h_ratio->Draw("same plc pmc");
1010         l_data->AddEntry(h_ratio);
1011     }
1012     pads->at(1)->cd();
1013     pads->at(1)->SetLogy();
1014     l_data->Draw();
1015     l_header->Draw();
1016     bool herwig_high = false;
1017     if (herwig->GetMaximum() > pythia->GetMaximum() ) herwig_high=true;
1018     if (herwig_high) herwig->Draw("axis");
1019     else pythia->Draw("axis");
1020     herwig->Draw("plc pmc nostack same e1");
1021     pythia->Draw("plc pmc same nostack e1");
1022     pythia->GetXaxis()->SetTitle(herwig_obs->at(0)->GetXaxis()->GetTitle());
1023     pythia->GetYaxis()->SetTitle(herwig_obs->at(0)->GetYaxis()->GetTitle());
1024     herwig->GetYaxis()->SetRangeUser(herwig_obs->at(0)->GetMinimum()*0.33, herwig_obs->at(0)->GetMaximum()*3);
1025     obs_canv->cd();
1026     pads->at(0)->Draw();
1027     pads->at(1)->Draw();
1028     return;
1029 }
1030 void PlotPhotonObs(TH1F* herwig_obs, TH1F* pythia_obs, TCanvas* obs_canv, int i)
1031 {
1032     //just does the plotting of one variable at a time
1033     if(!herwig_obs || !pythia_obs) return;
1034     TLegend* l_data=new TLegend(0.6, 0.3, 1, 0.5);
1035     TLegend* l_header=new TLegend(0.8, 0.6, 1, 1);
1036     conf->SetsPhenixHeaderLegend(l_header, tag);
1037     conf->SetLegend(l_data);
1038     l_data->SetNColumns(3);
1039     std::vector<TPad*>* pads=conf->AddPads(obs_canv);
1040     conf->ScaleXS(herwig_obs, true);
1041     conf->ScaleXS(pythia_obs, false);
1042     pads->at(1)->cd();
1043     herwig_obs->SetMarkerStyle(i+20);
1044     herwig_obs->SetMarkerSize(2);
1045     if(i==5) herwig_obs->SetMarkerStyle(33);
1046     else if(i==6) herwig_obs->SetMarkerStyle(34);
1047     THStack* st=new THStack("s1", "s1");
1048     st->Add(herwig_obs);
1049     pythia_obs->SetMarkerStyle(i+24);
1050     pythia_obs->SetMarkerSize(2);
1051     st->Add(pythia_obs);
1052     st->Draw("axis");
1053     st->Draw("plc pmc nostack same e1");
1054     l_data->AddEntry(herwig_obs, std::format("Herwig {}", herwig_obs->GetTitle()).c_str());
1055     l_data->AddEntry(pythia_obs, std::format("Pythia {}", pythia_obs->GetTitle()).c_str());
1056     pads->at(0)->cd();
1057     TH1F* h_ratio = conf->GetRatioPlot(herwig_obs, pythia_obs);
1058     h_ratio->SetMarkerStyle(2*(i)+37);
1059     h_ratio->Draw("same");
1060     l_data->AddEntry(h_ratio);
1061     pads->at(1)->cd();
1062     l_data->Draw();
1063     obs_canv->cd();
1064     pads->at(0)->Draw();
1065     pads->at(1)->Draw();
1066     return;
1067 }
1068 void PlotPhotonObs(TH1I* herwig_obs, TH1I* pythia_obs, TCanvas* obs_canv, int i)
1069 {
1070     //just does the plotting of one variable at a tie
1071     if(!herwig_obs || !pythia_obs) return;
1072     TLegend* l_data=new TLegend(0.6, 0.3, 1, 0.5);
1073     TLegend* l_header=new TLegend(0.8, 0.6, 1, 1);
1074     conf->SetsPhenixHeaderLegend(l_header, tag);
1075     conf->SetLegend(l_data);
1076     l_data->SetNColumns(3);
1077     std::vector<TPad*>* pads=conf->AddPads(obs_canv);
1078     conf->ScaleXS(herwig_obs, true);
1079     conf->ScaleXS(pythia_obs, false);
1080     pads->at(1)->cd();
1081     herwig_obs->SetLineColor((i)*22+1);
1082     herwig_obs->SetMarkerColor(i*22+1);
1083     herwig_obs->SetMarkerStyle(i+20);
1084     herwig_obs->SetMarkerSize(2);
1085     if(i==5) herwig_obs->SetMarkerStyle(33);
1086     else if(i==6) herwig_obs->SetMarkerStyle(34);
1087     herwig_obs->Draw("same e1");
1088     pythia_obs->SetLineColor((i)*22+1);
1089     pythia_obs->SetMarkerColor((i)*22+1);
1090     pythia_obs->SetMarkerStyle(i+24);
1091     pythia_obs->SetMarkerSize(2);
1092     pythia_obs->Draw("same e1");
1093     l_data->AddEntry(herwig_obs, std::format("Herwig {}", herwig_obs->GetTitle()).c_str());
1094     l_data->AddEntry(pythia_obs, std::format("Pythia {}", pythia_obs->GetTitle()).c_str());
1095     pads->at(0)->cd();
1096     TH1I* h_ratio = conf->GetRatioPlot(herwig_obs, pythia_obs);
1097     h_ratio->SetMarkerStyle(2*(i)+37);
1098     h_ratio->Draw("same");
1099     l_data->AddEntry(h_ratio);
1100     pads->at(1)->cd();
1101     l_data->Draw();
1102     obs_canv->cd();
1103     pads->at(0)->Draw();
1104     pads->at(1)->Draw();
1105     return;
1106 }
1107 void PlotJetPlots(TDirectory* herwig_jets, TDirectory* pythia_jets, std::vector<TCanvas*>*Canvi)
1108 {
1109     //all Jets
1110     std::vector<TH1F*>* herwig_pt       = new std::vector<TH1F*> ();
1111     std::vector<TH1F*>* herwig_e        = new std::vector<TH1F*> ();
1112     std::vector<TH1F*>* herwig_phi      = new std::vector<TH1F*> ();
1113     std::vector<TH1F*>* herwig_eta      = new std::vector<TH1F*> ();
1114     std::vector<TH1I*>* herwig_n_jets   = new std::vector<TH1I*> ();
1115     std::vector<TH1I*>* herwig_n_comp   = new std::vector<TH1I*> ();
1116     std::array<std::vector<TH1F*>*, 4>* herwig_jet_array =new std::array<std::vector<TH1F*>*, 4>();
1117     herwig_jet_array->at(0)=herwig_pt;
1118     herwig_jet_array->at(1)=herwig_e;
1119     herwig_jet_array->at(2)=herwig_phi;
1120     herwig_jet_array->at(3)=herwig_eta;
1121     //lead Jets
1122     std::vector<TH1F*>* herwig_lead_pt  = new std::vector<TH1F*> ();
1123     std::vector<TH1F*>* herwig_lead_e   = new std::vector<TH1F*> ();
1124     std::vector<TH1F*>* herwig_lead_phi = new std::vector<TH1F*> ();
1125     std::vector<TH1F*>* herwig_lead_eta = new std::vector<TH1F*> ();
1126     std::vector<TH1I*>* herwig_lead_n_comp  = new std::vector<TH1I*> ();
1127     std::array<std::vector<TH1F*>*, 4>* herwig_lead_jet_array =new std::array<std::vector<TH1F*>*, 4>();
1128     herwig_lead_jet_array->at(0)=herwig_lead_pt;
1129     herwig_lead_jet_array->at(1)=herwig_lead_e;
1130     herwig_lead_jet_array->at(2)=herwig_lead_phi;
1131     herwig_lead_jet_array->at(3)=herwig_lead_eta;
1132 
1133     CollectJets(herwig_jets, herwig_jet_array, herwig_lead_jet_array, herwig_n_jets, herwig_n_comp, herwig_lead_n_comp);
1134     //all Jets
1135     std::vector<TH1F*>* pythia_pt       = new std::vector<TH1F*> ();
1136     std::vector<TH1F*>* pythia_e        = new std::vector<TH1F*> ();
1137     std::vector<TH1F*>* pythia_phi      = new std::vector<TH1F*> ();
1138     std::vector<TH1F*>* pythia_eta      = new std::vector<TH1F*> ();
1139     std::vector<TH1I*>* pythia_n_jets   = new std::vector<TH1I*> ();
1140     std::vector<TH1I*>* pythia_n_comp   = new std::vector<TH1I*> ();
1141     std::array<std::vector<TH1F*>*, 4>* pythia_jet_array =new std::array<std::vector<TH1F*>*, 4>();
1142     pythia_jet_array->at(0)=pythia_pt;
1143     pythia_jet_array->at(1)=pythia_e;
1144     pythia_jet_array->at(2)=pythia_phi;
1145     pythia_jet_array->at(3)=pythia_eta;
1146     //lead Jets
1147     std::vector<TH1F*>* pythia_lead_pt  = new std::vector<TH1F*> ();
1148     std::vector<TH1F*>* pythia_lead_e   = new std::vector<TH1F*> ();
1149     std::vector<TH1F*>* pythia_lead_phi = new std::vector<TH1F*> ();
1150     std::vector<TH1F*>* pythia_lead_eta = new std::vector<TH1F*> ();
1151     std::vector<TH1I*>* pythia_lead_n_comp  = new std::vector<TH1I*> ();
1152     std::array<std::vector<TH1F*>*, 4>* pythia_lead_jet_array =new std::array<std::vector<TH1F*>*, 4>();
1153     pythia_lead_jet_array->at(0)=pythia_lead_pt;
1154     pythia_lead_jet_array->at(1)=pythia_lead_e;
1155     pythia_lead_jet_array->at(2)=pythia_lead_phi;
1156     pythia_lead_jet_array->at(3)=pythia_lead_eta;
1157     CollectJets(pythia_jets, pythia_jet_array, pythia_lead_jet_array, pythia_n_jets, pythia_n_comp, pythia_lead_n_comp);
1158     TCanvas* J_pt=new TCanvas("jet_pt", "jet_pt");
1159     PlotJetObs(herwig_jet_array->at(0), pythia_jet_array->at(0), J_pt);
1160     J_pt->SetLogy();
1161     TCanvas* J_l_pt=new TCanvas("jet_l_pt", "jet_l_pt");
1162     PlotJetObs(herwig_lead_jet_array->at(0), pythia_lead_jet_array->at(0), J_l_pt);
1163     J_l_pt->SetLogy();
1164     TCanvas* J_e=new TCanvas("jet_e", "jet_e");
1165     PlotJetObs(herwig_jet_array->at(1), pythia_jet_array->at(1), J_e);
1166     J_e->SetLogy();
1167     TCanvas* J_l_e=new TCanvas("jet_l_e", "jet_l_e");
1168     PlotJetObs(herwig_lead_jet_array->at(1), pythia_lead_jet_array->at(1), J_l_e);
1169     J_l_e->SetLogy();
1170     TCanvas* J_phi=new TCanvas("jet_phi", "jet_phi");
1171     PlotJetObs(herwig_jet_array->at(2), pythia_jet_array->at(2), J_phi);
1172     TCanvas* J_l_phi=new TCanvas("jet_l_phi", "jet_l_phi");
1173     PlotJetObs(herwig_lead_jet_array->at(2), pythia_lead_jet_array->at(2), J_l_phi);
1174     TCanvas* J_eta=new TCanvas("jet_eta", "jet_eta");
1175     PlotJetObs(herwig_jet_array->at(3), pythia_jet_array->at(3), J_eta);
1176     TCanvas* J_l_eta=new TCanvas("jet_l_eta", "jet_l_eta");
1177     PlotJetObs(herwig_lead_jet_array->at(3), pythia_lead_jet_array->at(3), J_l_eta);
1178     TCanvas* J_n_comp=new TCanvas("jet_n_comp", "jet_n_comp");
1179     PlotJetObs(herwig_n_comp, pythia_n_comp, J_n_comp);
1180     TCanvas* J_l_n_comp=new TCanvas("jet_l_n_comp", "jet_l_n_comp");
1181     PlotJetObs(herwig_lead_n_comp, pythia_lead_n_comp, J_l_n_comp);
1182     TCanvas* J_N=new TCanvas("jet_N", "jet_N");
1183     PlotJetObs(herwig_n_jets, pythia_n_jets, J_N);
1184 
1185     Canvi->push_back(J_pt);
1186     Canvi->push_back(J_l_pt);
1187     Canvi->push_back(J_e);
1188     Canvi->push_back(J_l_e);
1189     Canvi->push_back(J_phi);
1190     Canvi->push_back(J_l_phi);
1191     Canvi->push_back(J_eta);
1192     Canvi->push_back(J_l_eta);
1193     Canvi->push_back(J_n_comp);
1194     Canvi->push_back(J_l_n_comp);
1195     Canvi->push_back(J_N);
1196     
1197     storage->cd();
1198     TDirectory* hwj=(TDirectory*)storage->mkdir("Herwig_jets");
1199     hwj->cd();
1200     for(int i=0; i<(int)herwig_pt->size(); i++) herwig_pt->at(i)->Write();
1201     for(int i=0; i<(int)herwig_lead_pt->size(); i++) herwig_lead_pt->at(i)->Write();
1202     storage->cd();
1203     TDirectory* pj=(TDirectory*)storage->mkdir("Pythia_jets");
1204     pj->cd();
1205     for(int i=0; i<(int)pythia_pt->size(); i++) pythia_pt->at(i)->Write();
1206     for(int i=0; i<(int)pythia_lead_pt->size(); i++) pythia_lead_pt->at(i)->Write();
1207     storage->cd();
1208     
1209     return;
1210     
1211 }
1212 
1213 void PlotPhotonPlots(TDirectory* herwig_photons, TDirectory* pythia_photons, std::vector<TCanvas*>*Canvi)
1214 {
1215     //all Photons
1216     TH1F* herwig_pt     {};
1217     TH1F* herwig_e      {};
1218     TH1F* herwig_phi        {};
1219     TH1F* herwig_eta        {};
1220     TH1I* herwig_n_photons  {};
1221     std::array<TH1F*, 4>* herwig_photon_array =new std::array<TH1F*, 4>();
1222     herwig_photon_array->at(0)=herwig_pt;
1223     herwig_photon_array->at(1)=herwig_e;
1224     herwig_photon_array->at(2)=herwig_phi;
1225     herwig_photon_array->at(3)=herwig_eta;
1226     //lead Photons
1227     TH1F* herwig_lead_pt    {};
1228     TH1F* herwig_lead_e {};
1229     TH1F* herwig_lead_phi   {};
1230     TH1F* herwig_lead_eta   {};
1231     std::array<TH1F*, 4>* herwig_lead_photon_array =new std::array<TH1F*, 4>();
1232     herwig_lead_photon_array->at(0)=herwig_lead_pt;
1233     herwig_lead_photon_array->at(1)=herwig_lead_e;
1234     herwig_lead_photon_array->at(2)=herwig_lead_phi;
1235     herwig_lead_photon_array->at(3)=herwig_lead_eta;
1236     //direct Photons
1237     TH1F* herwig_direct_pt      {};
1238     TH1F* herwig_direct_e       {};
1239     TH1F* herwig_direct_phi     {};
1240     TH1F* herwig_direct_eta     {};
1241     TH1I* herwig_direct_n_photons   {};
1242     std::array<TH1F*, 4>* herwig_direct_photon_array =new std::array<TH1F*, 4>();
1243     herwig_direct_photon_array->at(0)=herwig_direct_pt;
1244     herwig_direct_photon_array->at(1)=herwig_direct_e;
1245     herwig_direct_photon_array->at(2)=herwig_direct_phi;
1246     herwig_direct_photon_array->at(3)=herwig_direct_eta;
1247     //all Photons
1248     TH1F* herwig_frag_pt        {};
1249     TH1F* herwig_frag_e     {};
1250     TH1F* herwig_frag_phi       {};
1251     TH1F* herwig_frag_eta       {};
1252     TH1I* herwig_frag_n_photons {};
1253     std::array<TH1F*, 4>* herwig_frag_photon_array =new std::array<TH1F*, 4>();
1254     herwig_frag_photon_array->at(0)=herwig_frag_pt;
1255     herwig_frag_photon_array->at(1)=herwig_frag_e;
1256     herwig_frag_photon_array->at(2)=herwig_frag_phi;
1257     herwig_frag_photon_array->at(3)=herwig_frag_eta;
1258     CollectPhotons(herwig_photons, herwig_photon_array, herwig_lead_photon_array, herwig_direct_photon_array, herwig_frag_photon_array, herwig_n_photons, herwig_direct_n_photons, herwig_frag_n_photons);
1259     
1260     //all Photons
1261     TH1F* pythia_pt     {};
1262     TH1F* pythia_e      {};
1263     TH1F* pythia_phi        {};
1264     TH1F* pythia_eta        {};
1265     TH1I* pythia_n_photons  {};
1266     std::array<TH1F*, 4>* pythia_photon_array =new std::array<TH1F*, 4>();
1267     pythia_photon_array->at(0)=pythia_pt;
1268     pythia_photon_array->at(1)=pythia_e;
1269     pythia_photon_array->at(2)=pythia_phi;
1270     pythia_photon_array->at(3)=pythia_eta;
1271     //lead Photons
1272     TH1F* pythia_lead_pt    {};
1273     TH1F* pythia_lead_e {};
1274     TH1F* pythia_lead_phi   {};
1275     TH1F* pythia_lead_eta   {};
1276     std::array<TH1F*, 4>* pythia_lead_photon_array =new std::array<TH1F*, 4>();
1277     pythia_lead_photon_array->at(0)=pythia_lead_pt;
1278     pythia_lead_photon_array->at(1)=pythia_lead_e;
1279     pythia_lead_photon_array->at(2)=pythia_lead_phi;
1280     pythia_lead_photon_array->at(3)=pythia_lead_eta;
1281     //direct Photons
1282     TH1F* pythia_direct_pt      {};
1283     TH1F* pythia_direct_e       {};
1284     TH1F* pythia_direct_phi     {};
1285     TH1F* pythia_direct_eta     {};
1286     TH1I* pythia_direct_n_photons   {};
1287     std::array<TH1F*, 4>* pythia_direct_photon_array =new std::array<TH1F*, 4>();
1288     pythia_direct_photon_array->at(0)=pythia_direct_pt;
1289     pythia_direct_photon_array->at(1)=pythia_direct_e;
1290     pythia_direct_photon_array->at(2)=pythia_direct_phi;
1291     pythia_direct_photon_array->at(3)=pythia_direct_eta;
1292     //all Photons
1293     TH1F* pythia_frag_pt        {};
1294     TH1F* pythia_frag_e     {};
1295     TH1F* pythia_frag_phi       {};
1296     TH1F* pythia_frag_eta       {};
1297     TH1I* pythia_frag_n_photons {};
1298     std::array<TH1F*, 4>* pythia_frag_photon_array =new std::array<TH1F*, 4>();
1299     pythia_frag_photon_array->at(0)=pythia_frag_pt;
1300     pythia_frag_photon_array->at(1)=pythia_frag_e;
1301     pythia_frag_photon_array->at(2)=pythia_frag_phi;
1302     pythia_frag_photon_array->at(3)=pythia_frag_eta;
1303     CollectPhotons(pythia_photons, pythia_photon_array, pythia_lead_photon_array, pythia_direct_photon_array, pythia_frag_photon_array, pythia_n_photons, pythia_direct_n_photons, pythia_frag_n_photons);
1304     
1305     TCanvas* J_pt=new TCanvas("photon_pt", "photon_pt");
1306     PlotPhotonObs(herwig_photon_array->at(0), pythia_photon_array->at(0), J_pt, 0);
1307     TCanvas* J_l_pt=new TCanvas("photon_l_pt", "photon_l_pt");
1308     PlotPhotonObs(herwig_lead_photon_array->at(0), pythia_lead_photon_array->at(0), J_l_pt, 1);
1309     TCanvas* J_e=new TCanvas("photon_e", "photon_e");
1310     PlotPhotonObs(herwig_photon_array->at(1), pythia_photon_array->at(1), J_e, 0);
1311     TCanvas* J_l_e=new TCanvas("photon_l_e", "photon_l_e");
1312     PlotPhotonObs(herwig_lead_photon_array->at(1), pythia_lead_photon_array->at(1), J_l_e, 1);
1313     TCanvas* J_phi=new TCanvas("photon_phi", "photon_phi");
1314     PlotPhotonObs(herwig_photon_array->at(2), pythia_photon_array->at(2), J_phi, 0);
1315     TCanvas* J_l_phi=new TCanvas("photon_l_pt", "photon_l_pt");
1316     PlotPhotonObs(herwig_lead_photon_array->at(2), pythia_lead_photon_array->at(2), J_l_phi, 1);
1317     TCanvas* J_eta=new TCanvas("photon_eta", "photon_eta");
1318     PlotPhotonObs(herwig_photon_array->at(3), pythia_photon_array->at(3), J_eta, 0);
1319     TCanvas* J_l_eta=new TCanvas("photon_l_eta", "photon_l_eta");
1320     PlotPhotonObs(herwig_lead_photon_array->at(3), pythia_lead_photon_array->at(3), J_l_eta, 1);
1321     TCanvas* J_N=new TCanvas("photon_N", "photon_N");
1322     PlotPhotonObs(herwig_n_photons, pythia_n_photons, J_N, 0);
1323 
1324     Canvi->push_back(J_pt);
1325     Canvi->push_back(J_l_pt);
1326     Canvi->push_back(J_e);
1327     Canvi->push_back(J_l_e);
1328     Canvi->push_back(J_phi);
1329     Canvi->push_back(J_l_phi);
1330     Canvi->push_back(J_eta);
1331     Canvi->push_back(J_l_eta);
1332     Canvi->push_back(J_N);
1333     
1334     TCanvas* J_d_pt=new TCanvas("direct_photon_pt", "direct_photon_pt");
1335     PlotPhotonObs(herwig_direct_photon_array->at(0), pythia_direct_photon_array->at(0), J_d_pt, 2);
1336     TCanvas* J_f_pt=new TCanvas("frag_photon_pt", "frag_photon_pt");
1337     PlotPhotonObs(herwig_frag_photon_array->at(0), pythia_frag_photon_array->at(0), J_f_pt, 3);
1338     TCanvas* J_d_e=new TCanvas("direct_photone", "direct_photone");
1339     PlotPhotonObs(herwig_direct_photon_array->at(1), pythia_direct_photon_array->at(1), J_d_e, 2);
1340     TCanvas* J_f_e=new TCanvas("frag_photon_e", "frag_photon_e");
1341     PlotPhotonObs(herwig_frag_photon_array->at(1), pythia_frag_photon_array->at(1), J_f_e, 3);
1342     TCanvas* J_d_phi=new TCanvas("direct_photonphi", "direct_photonphi");
1343     PlotPhotonObs(herwig_direct_photon_array->at(2), pythia_direct_photon_array->at(2), J_d_phi, 2);
1344     TCanvas* J_f_phi=new TCanvas("frag_photon_pt", "frag_photon_pt");
1345     PlotPhotonObs(herwig_frag_photon_array->at(2), pythia_frag_photon_array->at(2), J_f_phi, 3);
1346     TCanvas* J_d_eta=new TCanvas("direct_photoneta", "direct_photoneta");
1347     PlotPhotonObs(herwig_direct_photon_array->at(3), pythia_direct_photon_array->at(3), J_d_eta, 2);
1348     TCanvas* J_f_eta=new TCanvas("frag_photon_eta", "frag_photon_eta");
1349     PlotPhotonObs(herwig_frag_photon_array->at(3), pythia_frag_photon_array->at(3), J_f_eta, 3);
1350     TCanvas* J_d_N=new TCanvas("direct_photonN", "direct_photonN");
1351     PlotPhotonObs(herwig_direct_n_photons, pythia_direct_n_photons, J_d_N, 2);
1352     TCanvas* J_f_N=new TCanvas("frag_photonN", "frag_photonN");
1353     PlotPhotonObs(herwig_direct_n_photons, pythia_frag_n_photons, J_f_N, 3);
1354 
1355     Canvi->push_back(J_d_pt);
1356     Canvi->push_back(J_f_pt);
1357     Canvi->push_back(J_d_e);
1358     Canvi->push_back(J_f_e);
1359     Canvi->push_back(J_d_phi);
1360     Canvi->push_back(J_f_phi);
1361     Canvi->push_back(J_d_eta);
1362     Canvi->push_back(J_f_eta);
1363     Canvi->push_back(J_d_N);
1364     Canvi->push_back(J_f_N);
1365     return;
1366     
1367 }
1368 void PlotPhotonJetPlots(TDirectory* herwig_photons, TDirectory* pythia_photons, std::vector<TCanvas*>* Canvi)
1369 {
1370     //mainly 2D plots so I have to change the approach a bit here 
1371     TDirectory* H_pj = (TDirectory*)herwig_photons->GetDirectory("Photon-Jets");
1372     std::vector<TH2F*>* herwig_pt       = new std::vector<TH2F*> ();
1373     std::vector<TH2F*>* herwig_eta      = new std::vector<TH2F*> ();
1374     std::vector<TH2F*>* herwig_phi      = new std::vector<TH2F*> ();
1375     std::vector<TH2F*>* herwig_e        = new std::vector<TH2F*> ();
1376     std::vector<TH1F*>* herwig_dphi     = new std::vector<TH1F*> ();
1377     //CollectPhotonJets(H_pj, herwig_pt, herwig_eta, herwig_phi, herwig_e, herwig_dphi);
1378 
1379     TDirectory* P_pj = (TDirectory*)pythia_photons->GetDirectory("Photon-Jets");
1380     std::vector<TH2F*>* pythia_pt       = new std::vector<TH2F*> ();
1381     std::vector<TH2F*>* pythia_eta      = new std::vector<TH2F*> ();
1382     std::vector<TH2F*>* pythia_phi      = new std::vector<TH2F*> ();
1383     std::vector<TH2F*>* pythia_e        = new std::vector<TH2F*> ();
1384     std::vector<TH1F*>* pythia_dphi     = new std::vector<TH1F*> ();
1385     //CollectPhotonJets(P_pj, pythia_pt, pythia_eta, pythia_phi, pythia_e, pythia_dphi);
1386     
1387     TCanvas* J_pt=new TCanvas("photon_jet_pt");
1388     TCanvas* J_h_pt=new TCanvas("herwig_photon_jet_pt");
1389     TCanvas* J_p_pt=new TCanvas("pythia_photon_jet_pt");
1390     PlotPhotonJetObs(herwig_pt, pythia_pt, J_pt, J_h_pt, J_p_pt); 
1391     
1392     TCanvas* J_eta=new TCanvas("photon_jet_eta");
1393     TCanvas* J_h_eta=new TCanvas("photon_herwig_jet_eta");
1394     TCanvas* J_p_eta=new TCanvas("pythia_photon_jet_eta");
1395     PlotPhotonJetObs(herwig_eta, pythia_eta, J_eta, J_h_eta, J_p_eta); 
1396 
1397     TCanvas* J_phi=new TCanvas("photon_herwig_jet_phi");
1398     TCanvas* J_h_phi=new TCanvas("photon_herwig_jet_phi");
1399     TCanvas* J_p_phi=new TCanvas("pythia_photon_jet_phi");
1400     PlotPhotonJetObs(herwig_phi, pythia_phi, J_phi, J_h_phi, J_p_phi); 
1401 
1402     TCanvas* J_e=new TCanvas("photon_herwig_jet_e");
1403     TCanvas* J_h_e=new TCanvas("photon_herwig_jet_e");
1404     TCanvas* J_p_e=new TCanvas("pythia_photon_jet_e");
1405     PlotPhotonJetObs(herwig_e, pythia_e, J_e, J_h_e, J_p_e); 
1406 
1407         TCanvas* J_dphi=new TCanvas("photon_jet_dphi");
1408     PlotJetObs(herwig_dphi, pythia_dphi, J_dphi);
1409     Canvi->push_back(J_dphi);
1410     return;
1411 }
1412 void DoAllThePlotting(TFile* herwig_file, TFile* pythia_file, std::string trigger_tag="Jet30", float herwig_xs=1., float pythia_xs=1.)
1413 {
1414     //this does all the plotting of Pythia vrs herwig
1415     SetsPhenixStyle();  
1416     storage=new TFile(std::format("{}_scaled_pt.root", trigger_tag).c_str(), "RECREATE");
1417     conf =  new HerwigQAPlottingConfig(herwig_xs, pythia_xs);
1418     style_points=new Skaydis_colors();
1419     gStyle->SetPalette(100, style_points->Trans_gradient_PT);   
1420     conf->ExtractType(herwig_file);
1421     std::map<std::string, TDirectory*> top_dirs;
1422     std::vector<TCanvas*>* Canvi =new std::vector<TCanvas*> ();
1423     tag=trigger_tag;
1424     //find the relevant directories
1425     if(conf->isJet() || trigger_tag.find("MB") != std::string::npos){
1426         TDirectory* HJets=(TDirectory*)herwig_file->GetDirectory("Jets");
1427         TDirectory* PJets=(TDirectory*)pythia_file->GetDirectory("Jets");
1428         top_dirs["herwig jets"]=HJets;
1429         top_dirs["pythia jets"]=PJets;
1430         PlotJetPlots(HJets, PJets, Canvi);
1431     }
1432     if(conf->isPhoton() || trigger_tag.find("MB") != std::string::npos )
1433     {
1434         TDirectory* HPh=(TDirectory*)herwig_file->GetDirectory("Photons");
1435         TDirectory* PPh=(TDirectory*)pythia_file->GetDirectory("Photons");
1436         top_dirs["herwig photons"]=HPh;
1437         top_dirs["pythia photons"]=PPh;
1438         PlotPhotonPlots(HPh, PPh, Canvi);
1439     }
1440     if(conf->isPhoton() && conf->isJet() && trigger_tag.find("MB") == std::string::npos)
1441     {
1442         PlotJetPlots(top_dirs["herwig photons"], top_dirs["pythia photons"], Canvi);
1443     }
1444     TDirectory* HEvent=(TDirectory*)herwig_file->GetDirectory("Event");
1445     TDirectory* PEvent=(TDirectory*)pythia_file->GetDirectory("Event");
1446     top_dirs["herwig event"]=HEvent;
1447     top_dirs["pythia event"]=PEvent;    
1448     std::cout<<__LINE__<<std::endl;
1449     PlotEventPlots(HEvent, PEvent, Canvi);
1450     std::cout<<__LINE__<<std::endl;
1451     conf->DumpHistos(std::format("~/{}_output.pdf", trigger_tag), Canvi);
1452     storage->cd();
1453     storage->Write();
1454     storage->Close();
1455     return;
1456 }
1457 
1458 #endif