Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-05-23 08:12:16 

0001 #pragma once 
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 
0004 #include <TTree.h>
0005 #include <TLegend.h>
0006 #include <TCanvas.h>
0007 #include <TH1.h>
0008 #include <TH2.h>
0009 #include <TExec.h>
0010 
0011 #include "Special_colors.h"
0012 //#include "PlottingConfigs.h"
0013 #include "sPhenixStyle.h"
0014 #include "sPhenixStyle.C"
0015 //#include "DrawingMacro.C"
0016  
0017 //#include <vandyclasses/EventInfo.h>
0018 //#include <vandyclasses/JetInfo.h>
0019 //#include <vandyclasses/Tower.h>
0020 
0021 #include <vector>
0022 #include <array>
0023 #include <math.h>
0024 #include <string>
0025 #include <numbers>
0026 #include <utility>
0027 
0028 R__LOAD_LIBRARY(libVandyClasses.so);
0029 
0030 //PlottingConfigs* conf {nullptr};
0031 Skaydis_colors* stylepoints {nullptr};
0032 
0033 TExec* les {nullptr};
0034 TExec* bi {nullptr};
0035 TExec* trans {nullptr};
0036 TExec* enby {nullptr};
0037 std::string tag{""};
0038 bool isData {false}; 
0039 
0040 TH1F* wEECtruthtowers   {nullptr};
0041 TH1F* wEECtruthparticles{nullptr};
0042 TH1F* wEECrecotowers    {nullptr};
0043 
0044 TH1F* pairtruthtowers   {nullptr};
0045 TH1F* pairtruthparticles{nullptr};
0046 TH1F* pairrecotowers    {nullptr};
0047 
0048 const std::vector<double> new_dPhiBins = {
0049 //    0.0000000, 0.041592654, 0.14159265,
0050     0.0000000, 0.14159265,
0051     0.34159265, 0.54159265, 0.74159265, 0.94159265,
0052     1.2415927, 1.5415927, 1.8415927, 2.1415927,
0053     2.3415927, 2.5415927, 2.7415927, 2.9415927,
0054     3.0415927, 3.1415927
0055 };
0056 std::vector<float> pairweights = {0.,5e-5};
0057 void makePairweightbins()
0058 {
0059     float upper = 2;
0060     float second_to_upper = 1;
0061     float loglower = std::log10(pairweights[1]);
0062     float logupper = std::log10(second_to_upper);
0063     int nlogbins = 17;
0064     float delta = (logupper-loglower)/((float) nlogbins-2);
0065     for(int i=1; i < nlogbins; i++){
0066            float binedge = std::pow(10, loglower + i*delta);
0067            if( binedge == pairweights[1] ) continue;
0068            else if (binedge >= second_to_upper) break;
0069            pairweights.push_back(binedge);
0070     }
0071     pairweights.push_back(second_to_upper);
0072     pairweights.push_back(upper);
0073     return;
0074 }
0075 void InitJetHistos(
0076         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* holding_array,
0077         std::vector<JetInfo>* jt,
0078         TH1F* pt, 
0079         int type,
0080         std::string nametag,
0081         std::string tag
0082         )
0083 {
0084     pt=new TH1F(std::format("h_{0}_jet_r0{1}", nametag, type+2).c_str(), 
0085             std::format("{} Jet R=0.{}; p_{{T}} [GeV]; #sigma", tag, type+2).c_str(),
0086                25, -0.5, 99.5 );
0087     std::pair<std::vector<JetInfo>*, TH1F*> ha = std::make_pair(jt, pt);
0088     holding_array->at(type)=&ha;
0089     return; 
0090 }
0091 void InitTowerHistos(
0092         TH1F* h_tower,
0093         std::string nametag, 
0094         std::string tag
0095         )
0096 {
0097     h_tower = new TH1F(std::format("h_{}_e", nametag).c_str(), 
0098             std::format("{} Energy; E [GeV]; N_{{towers}}", tag).c_str(),
0099             1000, 1e-4, 10);
0100     return; 
0101 }
0102 
0103 void InitwEECHistos()
0104 {
0105 
0106     wEECtruthtowers     = new TH1F(
0107         "h_truth_wEEC", 
0108         "Truth tower wEEC; #Delta #varphi; wEEC", 
0109         32, new_dPhiBins.data() );
0110 
0111     wEECtruthparticles  = new TH1F(
0112         "h_part_wEEC", 
0113         "Truth particle level wEEC; #Delta #varphi; wEEC", 
0114         32, new_dPhiBins.data() );
0115 
0116     wEECrecotowers      = new TH1F(
0117         "h_reco_wEEC", 
0118         "Reco tower wEEC; #Delta #varphi; wEEC", 
0119         32, new_dPhiBins.data() );
0120     return;
0121 }
0122 
0123 void InitPairHistos()
0124 {
0125     pairtruthtowers     = new TH1F(
0126         "h_truth_pair", 
0127         "Truth tower pair; #frac{E_{T,1}E_{T,2}}{p_{T, avg}^{2}}; N_{pair}", 
0128         pairweights.size(), pairweights.data() );
0129 
0130     pairtruthparticles  = new TH1F(
0131         "h_part_pair", 
0132         "Truth particle level pair; #frac{E_{T,1}E_{T,2}}{p_{T, avg}^{2}}; N_{pair}", 
0133         pairweights.size(), pairweights.data() );
0134 
0135     pairrecotowers      = new TH1F(
0136         "h_reco_pair", 
0137         "Reco tower pair; #frac{E_{T,1}E_{T,2}}{p_{T, avg}^{2}}; N_{pair}", 
0138         pairweights.size(), pairweights.data() );
0139     return;
0140 }
0141 
0142 void SetJetBranches(
0143         TTree* t, 
0144         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* truthjets, 
0145         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* towerjets, 
0146         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* truthtowerjets
0147         )
0148 {
0149     if(!isData)
0150     {
0151         for(int i=0; i<(int)truthjets->size(); i++)
0152         {
0153             if(!truthjets->at(i)) continue;
0154             t->SetBranchAddress(std::format("TruthJetInfo_r0{}", i+2).c_str(), &(truthjets->at(i)->first));
0155         }
0156         for(int i=0; i<(int)truthtowerjets->size(); i++)
0157         {
0158             if(!t->FindBranch(std::format("TruthTowerJetInfo_r0{}", i+2).c_str())) break;
0159             if(!truthtowerjets->at(i) ) continue;
0160             t->SetBranchAddress(std::format("TruthTowerJetInfo_r0{}", i+2).c_str(), &(truthtowerjets->at(i)->first));
0161         }
0162     }
0163     for(int i=0; i<(int)towerjets->size(); i++)
0164     {
0165         if(!towerjets->at(i)) continue;
0166         t->SetBranchAddress(std::format("JetInfo_r0{}", i+2).c_str(), &(towerjets->at(i)->first));
0167     }
0168     return; 
0169 }
0170 
0171 void SetTowerBranches(
0172         TTree* t,
0173         std::vector<Tower>* particle,
0174         std::vector<Tower>* truth_tower,
0175         std::vector<Tower>* reco_tower
0176         )
0177 {
0178     if(isData)
0179     {
0180         if(particle)    t->SetBranchAddress("TruthParticles", &particle);
0181         if(truth_tower) t->SetBranchAddress("TruthTowers", &truth_tower);
0182     }
0183     
0184     if(reco_tower)  t->SetBranchAddress("TowerInfo", &reco_tower);
0185     return;
0186 }
0187 
0188 void getwEEC(std::vector<Tower>* towers, TH1F* output_hist, float Q2)
0189 {
0190     for(int i=0; i<(int)towers->size()-1; i++)
0191     {
0192         float pt1 = std::pow(towers->at(i).px(), 2) + std::pow(towers->at(i).py(), 2);
0193         float phi1 = std::atan2(towers->at(i).py(), towers->at(i).px());
0194         pt1 = std::sqrt(pt1);
0195         for(int j=i+1; j<(int)towers->size(); j++)
0196         {
0197             float pt2 = std::pow(towers->at(j).px(), 2) + std::pow(towers->at(j).py(), 2);
0198             pt2 = std::sqrt(pt2);
0199             float phi2 = std::atan2(towers->at(j).py(), towers->at(j).px());
0200             float pw = pt1 * pt2/Q2;
0201             float dphi = phi1 - phi2;
0202             if(dphi < M_PI ) dphi+=2*M_PI;
0203             if(dphi > M_PI ) dphi-=2*M_PI;
0204             output_hist->Fill(dphi, pw);
0205         }
0206     }
0207     return;
0208 }
0209 
0210 void eventLoop(
0211         int event_n,    
0212         TTree* t, 
0213         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>*    truthjets, 
0214         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>*    towerjets, 
0215         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>*    truthtowerjets,
0216         std::pair<std::vector<Tower>*, TH1F*>*      truthparticles,
0217         std::pair<std::vector<Tower>*, TH1F*>*      truthtowers,
0218         std::pair<std::vector<Tower>*, TH1F*>*      recotowers,
0219         EventInfo*                  event
0220           )
0221 {
0222     std::cout<<"Event " <<event_n <<" out of " <<t->GetEntries() <<std::endl;
0223     t->GetEntry(event_n);
0224     std::cout<<__LINE__<<std::endl;
0225     bool isDijet = true;
0226     float Q2 = 1.;
0227     for(int i=0; i< 5/*(int)event->is_dijet_event->size()*/; i++)
0228     {
0229         if(!event->is_dijet_event(i)) isDijet=false;
0230     }
0231     if(!isDijet) return;
0232     std::cout<<__LINE__<<std::endl;
0233 
0234     if(!isData)
0235     {
0236         bool isTruthDijet = true;
0237         for(int i=0; i< 5/*(int)event->is_dijetTruth_event()->size()*/; i++)
0238         {
0239             if(!event->is_dijetTruth_event(i)) isTruthDijet=false;
0240         }
0241         if(!isTruthDijet) return;
0242         Q2 = std::pow(event->get_leadTruth_pT(2),2) +std::pow( event->get_subleadTruth_pT(2), 2);
0243         for(int i=0; i<(int)truthjets->size(); i++)
0244         {
0245             if(!truthjets->at(i) || 
0246                     !truthjets->at(i)->second
0247                     ) continue;
0248             for(int j=0; j<(int) truthjets->at(i)->first->size(); j++)
0249                 truthjets->at(i)->second->Fill(truthjets->at(i)->first->at(j).pt());
0250         }
0251         getwEEC(truthparticles->first, wEECtruthparticles, Q2);
0252     std::cout<<__LINE__<<std::endl;
0253 
0254         for(int i=0; i<(int)truthtowerjets->size(); i++)
0255         {
0256             if(!t->FindBranch(std::format("TruthTowerJetInfo_r0{}", i+2).c_str())) break;
0257             if(!truthtowerjets->at(i) || 
0258                     !truthtowerjets->at(i)->second
0259                     ) continue;
0260             for(int j=0; j<(int) truthtowerjets->at(i)->first->size(); j++)
0261                 truthtowerjets->at(i)->second->Fill(truthtowerjets->at(i)->first->at(j).pt());
0262         }
0263     std::cout<<__LINE__<<std::endl;
0264 
0265         for(int i=0; i<(int)truthtowers->first->size(); i++)
0266         {
0267             if(!truthtowers || 
0268                     !truthtowers->first ||
0269                     !truthtowers->second
0270               ) continue;
0271             truthtowers->second->Fill(truthtowers->first->at(i).e());
0272         }
0273     std::cout<<__LINE__<<std::endl;
0274 
0275         getwEEC(truthtowers->first, wEECtruthtowers, Q2);
0276 
0277     }
0278     Q2 = std::pow(event->get_lead_pT(2),2) + std::pow(event->get_sublead_pT(2),2);
0279     std::cout<<__LINE__<<std::endl;
0280 
0281     for(int i=0; i<(int)towerjets->size(); i++)
0282     {
0283         if(!towerjets->at(i) || 
0284                 !towerjets->at(i)->second
0285                 ) continue;
0286             for(int j=0; j<(int) towerjets->at(i)->first->size(); j++)
0287                 towerjets->at(i)->second->Fill(towerjets->at(i)->first->at(j).pt());
0288     }
0289     std::cout<<__LINE__<<std::endl;
0290     for(int i=0; i<(int)recotowers->first->size(); i++)
0291     {
0292         if(!recotowers || 
0293                 !recotowers->first ||
0294                 !recotowers->second
0295           ) continue;
0296         recotowers->second->Fill(recotowers->first->at(i).e());
0297         
0298         getwEEC(recotowers->first, wEECrecotowers, Q2);
0299     }
0300     std::cout<<__LINE__<<std::endl;
0301     return; 
0302 }
0303 
0304 void DrawPlots(
0305         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>*    truthjets, 
0306         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>*    towerjets, 
0307         std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5>*    truthtowerjets,
0308         std::pair<std::vector<Tower>*, TH1F*>*      truthparticles,
0309         std::pair<std::vector<Tower>*, TH1F*>*      truthtowers,
0310         std::pair<std::vector<Tower>*, TH1F*>*      recotowers,
0311         float weight,
0312         std::string output_file_name
0313     )
0314 {
0315     TFile* output = new TFile(output_file_name.c_str(), "RECREATE");
0316     return;
0317 }
0318 
0319 int TruthSkimmerQA(std::string input_file_name, std::string outputQA_name, bool isD=false, std::string gen="Pythia")
0320 {
0321     tag         = gen;
0322     stylepoints     = new Skaydis_colors();
0323     isData      = isD;
0324     //setting up the pallets
0325     trans=new TExec("trans_pallet", "gStyle->SetPalette(100, style_points->Trans_gradient_PT);");
0326     bi=new TExec("bi_pallet", "gStyle->SetPalette(100, style_points->Trans_gradient_PT);");
0327     les=new TExec("les_pallet", "gStyle->SetPalette(100, style_points->Trans_gradient_PT);");
0328     enby=new TExec("enby_pallet", "gStyle->SetPalette(100, style_points->Trans_gradient_PT);");
0329     TFile* input_file = new TFile(input_file_name.c_str(), "READ");
0330     input_file->cd();
0331     //make the branches
0332     
0333     //truth jets
0334     std::vector<JetInfo>* truthjet02 {nullptr};
0335     std::vector<JetInfo>* truthjet03 {nullptr};
0336     std::vector<JetInfo>* truthjet04 {nullptr};
0337     std::vector<JetInfo>* truthjet05 {nullptr};
0338     std::vector<JetInfo>* truthjet06 {nullptr};
0339     //accompaning histograms 
0340     TH1F* truthjet02pt {nullptr};   
0341     TH1F* truthjet03pt {nullptr};   
0342     TH1F* truthjet04pt {nullptr};   
0343     TH1F* truthjet05pt {nullptr};   
0344     TH1F* truthjet06pt {nullptr};   
0345     
0346     //tower jets
0347     std::vector<JetInfo>* towerjet02 {nullptr};
0348     std::vector<JetInfo>* towerjet03 {nullptr};
0349     std::vector<JetInfo>* towerjet04 {nullptr};
0350     std::vector<JetInfo>* towerjet05 {nullptr};
0351     std::vector<JetInfo>* towerjet06 {nullptr};
0352     //accompaning histograms 
0353     TH1F* towerjet02pt {nullptr};   
0354     TH1F* towerjet03pt {nullptr};   
0355     TH1F* towerjet04pt {nullptr};   
0356     TH1F* towerjet05pt {nullptr};   
0357     TH1F* towerjet06pt {nullptr};   
0358     
0359     //truthtower jets
0360     std::vector<JetInfo>* truthtowerjet02 {nullptr};
0361     std::vector<JetInfo>* truthtowerjet03 {nullptr};
0362     std::vector<JetInfo>* truthtowerjet04 {nullptr};
0363     std::vector<JetInfo>* truthtowerjet05 {nullptr};
0364     std::vector<JetInfo>* truthtowerjet06 {nullptr};
0365     //accompaning histograms 
0366     TH1F* truthtowerjet02pt {nullptr};  
0367     TH1F* truthtowerjet03pt {nullptr};  
0368     TH1F* truthtowerjet04pt {nullptr};  
0369     TH1F* truthtowerjet05pt {nullptr};  
0370     TH1F* truthtowerjet06pt {nullptr};  
0371 
0372     //to have and to hold
0373     std::array< std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* truthjets = new std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5> {};
0374     std::array< std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* towerjets = new std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5> {};
0375     std::array< std::pair<std::vector<JetInfo>*, TH1F*>*, 5>* truthtowerjets = new std::array<std::pair<std::vector<JetInfo>*, TH1F*>*, 5> {};
0376 
0377     InitJetHistos(truthjets, truthjet02, truthjet02pt, 0, "truth", "Truth Jets");
0378     InitJetHistos(truthjets, truthjet03, truthjet03pt, 1, "truth", "Truth Jets");
0379     InitJetHistos(truthjets, truthjet04, truthjet04pt, 2, "truth", "Truth Jets");
0380     InitJetHistos(truthjets, truthjet05, truthjet05pt, 3, "truth", "Truth Jets");
0381     InitJetHistos(truthjets, truthjet06, truthjet06pt, 4, "truth", "Truth Jets");
0382 
0383     InitJetHistos(towerjets, towerjet02, towerjet02pt, 0, "reco_tower", "Reco Tower Jets");
0384     InitJetHistos(towerjets, towerjet03, towerjet03pt, 1, "reco_tower", "Reco Tower Jets");
0385     InitJetHistos(towerjets, towerjet04, towerjet04pt, 2, "reco_tower", "Reco Tower Jets");
0386     InitJetHistos(towerjets, towerjet05, towerjet05pt, 3, "reco_tower", "Reco Tower Jets");
0387     InitJetHistos(towerjets, towerjet06, towerjet06pt, 4, "reco_tower", "Reco Tower Jets");
0388 
0389     InitJetHistos(truthtowerjets, truthtowerjet02, truthtowerjet02pt, 0, "truth_tower", "Truth Tower Jets");
0390     InitJetHistos(truthtowerjets, truthtowerjet03, truthtowerjet03pt, 1, "truth_tower", "Truth Tower Jets");
0391     InitJetHistos(truthtowerjets, truthtowerjet04, truthtowerjet04pt, 2, "truth_tower", "Truth Tower Jets");
0392     InitJetHistos(truthtowerjets, truthtowerjet05, truthtowerjet05pt, 3, "truth_tower", "Truth Tower Jets");
0393     InitJetHistos(truthtowerjets, truthtowerjet06, truthtowerjet06pt, 4, "truth_tower", "Truth Tower Jets");
0394 
0395     TTree* t = (TTree*) input_file->Get("T");
0396     SetJetBranches(t, truthjets, towerjets, truthtowerjets);
0397     std::vector<Tower>* truthparticles_V = new std::vector<Tower> {};
0398     std::vector<Tower>* truthtowers_V = new std::vector<Tower> {};
0399     std::vector<Tower>* recotowers_V = new std::vector<Tower> {};
0400 
0401     SetTowerBranches(t, truthparticles_V, truthtowers_V, recotowers_V);
0402     TH1F* truthparticles_H  {nullptr};
0403     TH1F* truthtowers_H {nullptr};
0404     TH1F* recotowers_H  {nullptr};
0405 
0406     InitTowerHistos(truthparticles_H, "t_p", "Truth Particles");
0407     InitTowerHistos(truthtowers_H, "t_t", "Truth Towers");
0408     InitTowerHistos(recotowers_H, "r_t", "Reco Towers");
0409     std::pair<std::vector<Tower>*, TH1F*>* truthparticles 
0410         = new std::pair<std::vector<Tower>*, TH1F*> {truthparticles_V, truthparticles_H};
0411     std::pair<std::vector<Tower>*, TH1F*>* truthtowers 
0412         = new std::pair<std::vector<Tower>*, TH1F*> {truthtowers_V, truthtowers_H};
0413     std::pair<std::vector<Tower>*, TH1F*>* recotowers 
0414         = new std::pair<std::vector<Tower>*, TH1F*> {recotowers_V, recotowers_H};
0415     EventInfo* event = new EventInfo();
0416     t->SetBranchAddress("EventInfo", &event);
0417     int n = t->GetEntries();
0418     for(int i=1; i<n/1000; i++) 
0419         eventLoop(i, t, 
0420             truthjets, towerjets, truthtowerjets, 
0421             truthparticles, truthtowers, recotowers, 
0422             event); 
0423     t->GetEntry(2);
0424     float sigma = event->get_cross_section();
0425     float weight = sigma /(float)n;
0426     DrawPlots(
0427             truthjets, towerjets, truthtowerjets, 
0428             truthparticles, truthtowers, recotowers,
0429             weight, outputQA_name 
0430             );
0431     input_file->Close();
0432     return 0;
0433 };
0434 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!