Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:13 

0001 #include <TFile.h>
0002 #include <TH1.h>
0003 #include <TH2.h>
0004 #include <TCanvas.h>
0005 #include <TTreeReader.h>
0006 #include <TTreeReaderValue.h>
0007 #include <TTreeReaderArray.h>
0008 #include <vector>
0009 #include <iostream>
0010 #include <TDatabasePDG.h>
0011 #include <TString.h>
0012 #include <set>
0013 #include <TChain.h>
0014 #include <TGraph.h>
0015 #include <TLegend.h>
0016 #include <TMath.h>
0017 
0018 //#include "sPhenixStyle.h"
0019 //#include "sPhenixStyle.C"
0020 
0021 #include "JetAnalysis.h"
0022 #include <ranges>
0023 
0024 
0025 
0026 //JetAnalysis::JetAnalysis(){
0027 //}
0028 
0029 std::tuple<TCanvas*, TPad**, TH1D**> PrepareCanvas4Pad(UInt_t xAxisBins, Double_t *xAxis);
0030 
0031 bool CheckValue(ROOT::Internal::TTreeReaderValueBase& value) {
0032    if (value.GetSetupStatus() < 0) {
0033       std::cerr << "Error " << value.GetSetupStatus()
0034                 << "setting up reader for " << value.GetBranchName() << '\n';
0035       return false;
0036    }
0037    return true;
0038 }
0039 
0040 /*void JetAnalysis::MatchedJetContainer::Reset()
0041 {
0042     vtxtype = GlobalVertex::VTXTYPE::UNDEFINED;
0043     reco_pt = NAN;
0044     reco_jet_nChConstituents = NAN;
0045     reco_jet_nConstituents = NAN;
0046     reco_constituents_Sdxy.clear();
0047     reco_jet_Sdxy_1N = NAN;
0048     reco_jet_Sdxy_2N = NAN;
0049     reco_jet_Sdxy_3N = NAN;
0050     reco_jet_Sdxy_4N = NAN;
0051     truth_jet_flavour = {};
0052     truth_pt = NAN;
0053     truth_jet_nChConstituents = NAN;
0054     truth_jet_nConstituents = NAN;
0055 }*/
0056 
0057 bool JetAnalysis::CheckValue(ROOT::Internal::TTreeReaderValueBase& value) {
0058    if (value.GetSetupStatus() < 0) {
0059       std::cerr << "Error " << value.GetSetupStatus()
0060                 << "setting up reader for " << value.GetBranchName() << '\n';
0061       return false;
0062    }
0063    return true;
0064 }
0065 
0066 // Analyze the tree "MyTree" in the file passed into the function.
0067 // Returns false in case of errors.
0068 bool JetAnalysis::processCondor(const std::string &dataFiles) {
0069    TChain *tc = new TChain("AntiKt_r04/Data");
0070    tc->Add(dataFiles.data());
0071    TTreeReader reader(tc);
0072    
0073    jet_container_ = new TTreeReaderValue<FullJetFinder::Container>(reader,"JetContainer");
0074    TString filename = work_dir_ + "outtree_v2.root";
0075    //TFile *f = new TFile(filename,"RECREATE");
0076    //TTree *ttree_out = new TTree("Data", "Data");
0077    //MatchedJetContainer matched_jets_out;
0078    //ttree_out->Branch( "matched_jets", &matched_jets_out );
0079 
0080 
0081    TH1D *h_reco_all = new TH1D("h_reco_all","reco_all",400,0,100);
0082    h_reco_all->GetXaxis()->SetTitle("pt_jet_reco");
0083 
0084    TH1D *h_truth_all = new TH1D("h_truth_all","truth_all",400,0,100);
0085    h_truth_all->GetXaxis()->SetTitle("pt_jet^truth");
0086 
0087    TH1D *h_reco_unmatch = new TH1D("h_reco_unmatch","reco_unmatch",400,0,100);
0088    h_reco_unmatch->GetXaxis()->SetTitle("pt_jet^reco");
0089 
0090    TH1D *h_reco_match = new TH1D("h_reco_match","reco_match",400,0,100);
0091    h_reco_match->GetXaxis()->SetTitle("pt_jet^reco");
0092 
0093    TH1D *h_truth_match = new TH1D("h_truth_match","truth_match",400,0,100);
0094    h_truth_match->GetXaxis()->SetTitle("pt_jet^truth");
0095 
0096    TH1D *h_phi_all = new TH1D("h_phi_all","phi_all",800,-4,4);
0097    h_phi_all->GetXaxis()->SetTitle("phi");
0098 
0099    TH1D *h_phi_match = new TH1D("h_phi_match","phi_match",800,-4,4);
0100    h_phi_match->GetXaxis()->SetTitle("phi_match");
0101 
0102    TH1D *h_eta_all = new TH1D("h_eta_all","eta_all",300,-1.5,1.5);
0103    h_eta_all->GetXaxis()->SetTitle("eta");
0104 
0105    TH1D *h_eta_match = new TH1D("h_eta_match","eta_match",300,-1.5,1.5);
0106    h_eta_match->GetXaxis()->SetTitle("eta_match");
0107 
0108    TH1D *h_phi_truth_all = new TH1D("h_phi_truth_all","phi_truth_all",800,-4,4);
0109    h_phi_truth_all->GetXaxis()->SetTitle("phi");
0110 
0111    TH1D *h_phi_truth_match = new TH1D("h_phi_truth_match","phi_truth_match",800,-4,4);
0112    h_phi_truth_match->GetXaxis()->SetTitle("phi_match");
0113 
0114    TH1D *h_eta_truth_all = new TH1D("h_eta_truth_all","eta_truth_all",300,-1.5,1.5);
0115    h_eta_truth_all->GetXaxis()->SetTitle("eta");
0116 
0117    TH1D *h_eta_truth_match = new TH1D("h_eta_truth_match","eta_truth_match",300,-1.5,1.5);
0118    h_eta_match->GetXaxis()->SetTitle("eta_match");
0119 
0120 
0121    TH1D *h_eta_track = new TH1D("h_eta_track","h_eta_track",300,-1.5,1.5);
0122    h_eta_track->GetXaxis()->SetTitle("eta_match");
0123 
0124    TH1D *h_phi_track = new TH1D("h_phi_track","h_phi_track",800,-4,4);
0125    h_phi_track->GetXaxis()->SetTitle("phi");
0126 
0127    TH1D *h_eta_track_match = new TH1D("h_eta_track_match","h_eta_track_match",300,-1.5,1.5);
0128    h_eta_track_match->GetXaxis()->SetTitle("eta_match");
0129 
0130    TH1D *h_phi_track_match = new TH1D("h_phi_track_match","h_phi_track_match",800,-4,4);
0131    h_phi_track_match->GetXaxis()->SetTitle("phi");
0132 
0133    TH2D *h_phi_pt_track = new TH2D("h_phi_pt_track","h_phi_pt_track",60,0,30,800,-4,4);
0134    TH2D *h_dca_xy = new TH2D("h_dca_xy","dca_xy",60,0,30,100,-2,2);
0135    TH2D *h_dca_xy_unc = new TH2D("dca_xy_unc","dca_xy_unc",60,0,30,100,0,0.02);
0136    TH2D *h_nmvtx = new TH2D("h_nmvtx","h_nmvtx_all",60,0,30,6,0,6);
0137    TH2D *h_nintt = new TH2D("h_nintt","h_nintt_all",60,0,30,5,0,5);
0138    TH2D *h_ntpc = new TH2D("h_ntpc","h_ntpc_all",60,0,30,50,0,50);
0139    TH2D *h_chi2ndof = new TH2D("h_chi2ndof","h_chi2ndof_all",60,0,30,80,0,20);
0140 
0141    h_phi_pt_track->GetXaxis()->SetTitle("Track p_T");
0142    h_dca_xy->GetXaxis()->SetTitle("Track p_T");
0143    h_dca_xy_unc->GetXaxis()->SetTitle("Track p_T");
0144    h_nmvtx->GetXaxis()->SetTitle("Track p_T");
0145    h_nintt->GetXaxis()->SetTitle("Track p_T");
0146    h_ntpc->GetXaxis()->SetTitle("Track p_T");
0147    h_chi2ndof->GetXaxis()->SetTitle("Track p_T");
0148 
0149    TH2D *h_jet_area = new TH2D("h_jet_area","h_jet_area",60,0,1.2,60,0,1.2);
0150    h_jet_area->GetXaxis()->SetTitle("area_jet_truth");
0151    h_jet_area->GetYaxis()->SetTitle("area_jet_reco");
0152    
0153    int num = 5;
0154    TH1D *h_sDCA[num][3]; //variables, flavour (0 - lf, 1 - C, 2 - B)
0155     std::string tagName[5] = {"Sdxy","1_N","2_N","3_N","4_N"};//"dxy","dxy_unc",
0156     for(int i = 0; i<num;i++){ 
0157       TString name = "matching_"+tagName[i];
0158       h_sDCA[i][0] = new TH1D(name+"_LF",name+"_LF",100,-40,40);
0159       h_sDCA[i][1] = new TH1D(name+"_C",name+"_C",100,-40,40);
0160       h_sDCA[i][2] = new TH1D(name+"_B",name+"_B",100,-40,40);
0161       for(int j = 0; j<3;j++){ 
0162          h_sDCA[i][j]->GetYaxis()->SetTitle("Probability distribution");
0163          h_sDCA[i][j]->GetXaxis()->SetTitle("2D Impact Parameter significance");
0164       }
0165     }
0166 
0167    TCanvas *c2 = new TCanvas("c2","c2",2400,800);
0168    c2->Divide(5,1);
0169 
0170     TH1I *h_matching[4];
0171     std::string matName[4] = {"total","LF","C","B"};
0172     for(int i = 0; i<4;i++){ 
0173       TString name = "matching_"+matName[i];
0174       h_matching[i]= new TH1I(name,name,5,0,5);
0175       h_matching[i]->GetXaxis()->SetBinLabel(1,"Uniquely matched");
0176       h_matching[i]->GetXaxis()->SetBinLabel(2,"Unmatched R");
0177       h_matching[i]->GetXaxis()->SetBinLabel(3,"Unmatched T");
0178       h_matching[i]->GetXaxis()->SetBinLabel(4,"2T -> 1R");
0179       h_matching[i]->GetXaxis()->SetBinLabel(5,"1T -> 2R");
0180     }
0181 
0182     TH1I *h_pflow[4];
0183     std::string matName2[4] = {"Uniquely_matched","Unmatched_R","2T_to_1R","1T_to_2R"};
0184     for(int i = 0; i<4;i++){ 
0185       TString name = "h_pflow_"+matName2[i];
0186       h_pflow[i]= new TH1I(name,name,6,-1.5,4.5);
0187       h_pflow[i]->GetXaxis()->SetBinLabel(1,"UNASSIGNED");
0188       h_pflow[i]->GetXaxis()->SetBinLabel(2,"MATCHED_CHARGED_HADRON");
0189       h_pflow[i]->GetXaxis()->SetBinLabel(3,"UNMATCHED_CHARGED_HADRON");
0190       h_pflow[i]->GetXaxis()->SetBinLabel(4,"UNMATCHED_EM_PARTICLE");
0191       h_pflow[i]->GetXaxis()->SetBinLabel(5,"UNMATCHED_NEUTRAL_HADRON");
0192       h_pflow[i]->GetXaxis()->SetBinLabel(6,"LEFTOVER_EM_PARTICLE");
0193     }
0194 
0195    TCanvas *c3 = new TCanvas("c3","c3",2400,800);
0196     c3->Divide(8,1);
0197 
0198    TH1I *h_primvtx[6];
0199    std::string matName3[6] = {"All","z<10cm","matched_jets_to_vtx","c-jets","b-jets","unmatched R"};
0200    for(int i = 0; i<6;i++){ 
0201       TString name = "h_primvtx_"+matName3[i];
0202       h_primvtx[i]= new TH1I(name,name,3,-0.5,2.5);
0203       h_primvtx[i]->GetXaxis()->SetBinLabel(1,"SVTX_MBD");
0204       h_primvtx[i]->GetXaxis()->SetBinLabel(2,"SVTX");
0205       h_primvtx[i]->GetXaxis()->SetBinLabel(3,"MBD");
0206    }
0207 
0208    TH1I *h_n_vtx = new TH1I("h_n_vtx","number of primary vertices in event",5,0,5);
0209    TH1I *h_n_vtx_jet = new TH1I("h_n_vtx_jet","number of primary vertices inside jet",5,0,5);
0210 
0211    // Read a TriggerInfo object from the tree entries:
0212    // Now iterate through the TTree entries and fill a histogram.
0213 
0214   // int nentries = reader.GetEntries();
0215    int nentries = tc->GetEntries();
0216    std::cout<<"events "<<nentries<<std::endl;
0217 
0218 
0219    // Create a new ROOT file
0220 TFile *file = new TFile("/sphenix/tg/tg01/hf/jkvapil/JET30_R22_npile_full/HF_matched_jets.root", "RECREATE");
0221 
0222 // Create a new TTree
0223 TTree *tree = new TTree("Data", "A tree with various branches");
0224 
0225 // Define variables for branches
0226 float id;
0227 float vertex_x, vertex_y, vertex_z;
0228 float vertex_x_unc, vertex_y_unc, vertex_z_unc;
0229 float vertex_chisq;
0230 int vertex_ndf;
0231 
0232 float reco_area;
0233 int reco_num_ChConstituents;
0234 float reco_px, reco_py, reco_pz, reco_pt, reco_eta, reco_phi, reco_m, reco_e;
0235 float truth_area;
0236 int truth_num_ChConstituents;
0237 float truth_px, truth_py, truth_pz, truth_pt, truth_eta, truth_phi, truth_m, truth_e;
0238 int truth_flavour;
0239 
0240 // Define vectors for leaves
0241 std::vector<float> track_e, track_eta, track_phi, track_px, track_py, track_pz, track_pt;
0242 std::vector<int> track_charge;
0243 std::vector<float> track_DCA_xy, track_DCA_xy_unc, track_sDCA_xy, track_DCA3d, track_sDCA3d;
0244 std::vector<int> track_n_mvtx, track_n_intt, track_n_tpc;
0245 std::vector<float> track_chisq;
0246 std::vector<int> track_ndf;
0247 
0248 // Create branches
0249 tree->Branch("vertex_x", &vertex_x, "vertex_x/F");
0250 tree->Branch("vertex_y", &vertex_y, "vertex_y/F");
0251 tree->Branch("vertex_z", &vertex_z, "vertex_z/F");
0252 tree->Branch("vertex_x_unc", &vertex_x_unc, "vertex_x_unc/F");
0253 tree->Branch("vertex_y_unc", &vertex_y_unc, "vertex_y_unc/F");
0254 tree->Branch("vertex_z_unc", &vertex_z_unc, "vertex_z_unc/F");
0255 tree->Branch("vertex_chisq", &vertex_chisq, "vertex_chisq/F");
0256 tree->Branch("vertex_ndf", &vertex_ndf, "vertex_ndf/I");
0257 tree->Branch("id", &id, "id/F");
0258 tree->Branch("reco_area", &reco_area, "reco_area/F");
0259 tree->Branch("reco_num_ChConstituents", &reco_num_ChConstituents, "reco_num_ChConstituents/I");
0260 tree->Branch("reco_px", &reco_px, "reco_px/F");
0261 tree->Branch("reco_py", &reco_py, "reco_py/F");
0262 tree->Branch("reco_pz", &reco_pz, "reco_pz/F");
0263 tree->Branch("reco_pt", &reco_pt, "reco_pt/F");
0264 tree->Branch("reco_eta", &reco_eta, "reco_eta/F");
0265 tree->Branch("reco_phi", &reco_phi, "reco_phi/F");
0266 tree->Branch("reco_m", &reco_m, "reco_m/F");
0267 tree->Branch("reco_e", &reco_e, "reco_e/F");
0268 tree->Branch("reco_sdxy_1N", &reco_sdxy_1N, "reco_sdxy_1N/F");
0269 tree->Branch("reco_sdxy_2N", &reco_sdxy_2N, "reco_sdxy_2N/F");
0270 tree->Branch("reco_sdxy_3N", &reco_sdxy_3N, "reco_sdxy_3N/F");
0271 tree->Branch("truth_area", &truth_area, "truth_area/F");
0272 tree->Branch("truth_num_ChConstituents", &truth_num_ChConstituents, "truth_num_ChConstituents/I");
0273 tree->Branch("truth_px", &truth_px, "truth_px/F");
0274 tree->Branch("truth_py", &truth_py, "truth_py/F");
0275 tree->Branch("truth_pz", &truth_pz, "truth_pz/F");
0276 tree->Branch("truth_pt", &truth_pt, "truth_pt/F");
0277 tree->Branch("truth_eta", &truth_eta, "truth_eta/F");
0278 tree->Branch("truth_phi", &truth_phi, "truth_phi/F");
0279 tree->Branch("truth_m", &truth_m, "truth_m/F");
0280 tree->Branch("truth_e", &truth_e, "truth_e/F");
0281 tree->Branch("truth_flavour", &truth_flavour, "truth_flavour/I");
0282 
0283 
0284 // Create leaves
0285 tree->Branch("track_e","std::vector<float>", &track_e);
0286 tree->Branch("track_eta", &track_eta);
0287 tree->Branch("track_phi", &track_phi);
0288 tree->Branch("track_px", &track_px);
0289 tree->Branch("track_py", &track_py);
0290 tree->Branch("track_pz", &track_pz);
0291 tree->Branch("track_pt", &track_pt);
0292 tree->Branch("track_charge", &track_charge);
0293 tree->Branch("track_DCA_xy", &track_DCA_xy);
0294 tree->Branch("track_DCA_xy_unc", &track_DCA_xy_unc);
0295 tree->Branch("track_sDCA_xy", &track_sDCA_xy);
0296 tree->Branch("track_DCA3d", &track_DCA3d);
0297 tree->Branch("track_sDCA3d", &track_sDCA3d);
0298 tree->Branch("track_n_mvtx", &track_n_mvtx);
0299 tree->Branch("track_n_intt", &track_n_intt);
0300 tree->Branch("track_n_tpc", &track_n_tpc);
0301 tree->Branch("track_chisq", &track_chisq);
0302 tree->Branch("track_ndf", &track_ndf);
0303 
0304    //loop over events
0305    while (reader.Next()) {
0306       if (!CheckValue(*jet_container_)) return false;
0307       printProgress(reader.GetCurrentEntry(),nentries);
0308 
0309       //vectors to hold jet ID for matching
0310       std::vector<int> matchedR;
0311       std::vector<int> matchedT;
0312       std::vector<int> JetIDR;
0313       std::vector<int> JetIDT;
0314 
0315       //if jet does not pass a selection it will not reach here and then can have single jet with id 1 while missing id 0
0316       //this is a map between jetid and vector position
0317       std::map<int,int> JetIDR_map;
0318       std::map<int,int> JetIDT_map;
0319 
0320       //std::cout<<"recojetid: ";
0321 
0322       //IDs of reco jets
0323       int local_index = -1;
0324       for (auto reco_jet : (**jet_container_).recojets) {
0325          local_index++;
0326          if(reco_jet.pt < 10.0 || reco_jet.pt > 120.0) continue;
0327          if(std::abs(reco_jet.eta) > 0.7) continue;
0328          if(reco_jet.area < 3.1415*0.4*0.4*0.6) continue;
0329          JetIDR.push_back(reco_jet.id);
0330          JetIDR_map.insert(std::pair<int, int>(reco_jet.id, local_index));  
0331 
0332          h_reco_all->Fill(reco_jet.pt);
0333          h_eta_all->Fill(reco_jet.eta);
0334          h_phi_all->Fill(reco_jet.phi);
0335          for (auto chtrk : reco_jet.chConstituents){
0336             h_phi_track->Fill(chtrk.phi);
0337             h_eta_track->Fill(chtrk.eta);
0338          }
0339          
0340         // std::cout<<reco_jet.id<<" ";
0341       }
0342      // std::cout<<std::endl;
0343      // std::cout<<"truthjetid: ";
0344 
0345       //IDs of truth jet
0346       local_index = -1;
0347       for (auto truth_jet : (**jet_container_).truthjets) {
0348          local_index++;
0349          if(truth_jet.pt < 30.0 || truth_jet.pt > 120.0) continue;
0350          if(std::abs(truth_jet.eta) > 0.7) continue;
0351          if(truth_jet.area < 3.1415*0.4*0.4*0.6) continue;
0352          JetIDT.push_back(truth_jet.id);
0353          JetIDT_map.insert(std::pair<int, int>(truth_jet.id, local_index));
0354          h_truth_all->Fill(truth_jet.pt);
0355          h_eta_truth_all->Fill(truth_jet.eta);
0356          h_phi_truth_all->Fill(truth_jet.phi);
0357          //std::cout<<truth_jet.id<<" ";
0358       }
0359 
0360      // std::cout<<std::endl;
0361 
0362      
0363             
0364 
0365 
0366 
0367 
0368       //main matching loop
0369       std::vector<std::pair<int,int>> matchedPair;
0370       for (auto reco_jet : (**jet_container_).recojets) {
0371          if(reco_jet.pt < 10.0 || reco_jet.pt > 120.0) continue;
0372          if(std::abs(reco_jet.eta) > 0.7) continue;
0373          if(reco_jet.area < 3.1415*0.4*0.4*0.6) continue;
0374          for (auto truth_jet : (**jet_container_).truthjets) {
0375             if(truth_jet.pt < 30.0 || truth_jet.pt > 120.0) continue;
0376             if(std::abs(truth_jet.eta) > 0.7) continue;
0377             if(truth_jet.area < 3.1415*0.4*0.4*0.6) continue;
0378          
0379             double etaR = reco_jet.eta;
0380             double etaT = truth_jet.eta;
0381             double phiR = reco_jet.phi;
0382             double phiT = truth_jet.phi;
0383             double dR = std::sqrt((etaR-etaT)*(etaR-etaT)+(phiR-phiT)*(phiR-phiT));
0384             if (dR > 0.3){
0385                // unmatched
0386                continue;
0387             } 
0388             else{
0389                matchedPair.push_back(std::pair(reco_jet.id,truth_jet.id));
0390                matchedR.push_back(reco_jet.id);
0391                matchedT.push_back(truth_jet.id);
0392             }
0393          }
0394       }
0395 
0396       //are there unmatched reco jets?
0397       std::vector<int> unmatchedR;
0398       std::set_difference(JetIDR.begin(), JetIDR.end(), matchedR.begin(), matchedR.end(), std::inserter(unmatchedR, unmatchedR.begin()));
0399       //unmatched truth
0400       std::vector<int> unmatchedT;
0401       std::set_difference(JetIDT.begin(), JetIDT.end(), matchedT.begin(), matchedT.end(), std::inserter(unmatchedT, unmatchedT.begin()));
0402       //are there doubly matched jets?
0403       //2 reco -> 1 truth
0404       std::set<int> doublymatchedR = findDuplicates(matchedR);
0405       //truth -> 2 reco
0406       std::set<int> doublymatchedT = findDuplicates(matchedT);
0407 
0408       /*for (auto truth_jet : (**jet_container_).truthjets){  
0409          Flavour flavour = getFlavour(truth_jet);
0410          if(std::find(unmatchedT.begin(), unmatchedT.end(), truth_jet.id)!=unmatchedT.end()){
0411             h_matching[0]->Fill(2.5);
0412             if(flavour.isLF) h_matching[1]->Fill(2.5);
0413             if(flavour.isC & !flavour.isB) h_matching[2]->Fill(2.5);
0414             if(flavour.isB) h_matching[3]->Fill(2.5);
0415          }
0416       }*/
0417 //std::cout<<"matched pair: "<<std::endl;
0418       //matching stats
0419       for(auto &id : matchedPair){
0420          if(std::find(doublymatchedR.begin(), doublymatchedR.end(),id.first) != doublymatchedR.end()) continue;
0421          if(std::find(doublymatchedT.begin(), doublymatchedT.end(),id.second) != doublymatchedT.end()) continue;
0422          //std::cout<<"<"<<id.first<<","<<id.second<<"> ";
0423          Flavour flavour = getFlavour((**jet_container_).truthjets.at(JetIDT_map[id.second]));  
0424          h_matching[0]->Fill(0.5);
0425          if(flavour.isLF) h_matching[1]->Fill(0.5);
0426          if(flavour.isC & !flavour.isB) h_matching[2]->Fill(0.5);
0427          if(flavour.isB) h_matching[3]->Fill(0.5);
0428          h_reco_match->Fill((**jet_container_).recojets.at(JetIDR_map[id.first]).pt);
0429          h_truth_match->Fill((**jet_container_).truthjets.at(JetIDT_map[id.second]).pt);
0430          h_phi_match->Fill((**jet_container_).recojets.at(JetIDR_map[id.first]).phi);
0431          h_eta_match->Fill((**jet_container_).recojets.at(JetIDR_map[id.first]).eta);
0432          h_phi_truth_match->Fill((**jet_container_).truthjets.at(JetIDT_map[id.second]).phi);
0433          h_eta_truth_match->Fill((**jet_container_).truthjets.at(JetIDT_map[id.second]).eta);
0434       }
0435 //std::cout<<std::endl<<"Unmatch R ID "<<std::endl;
0436       for (auto& id : unmatchedR){
0437          //std::cout<<id<<" ";
0438          h_matching[0]->Fill(1.5);
0439          h_reco_unmatch->Fill((**jet_container_).recojets.at(JetIDR_map[id]).pt);
0440          //if((**jet_container_).recojets.at(JetIDR_map[id]).chConstituents.size() > 0){
0441             /*if((**jet_container_).primaryVertex.at((**jet_container_).recojets.at(JetIDR_map[id]).chConstituents.at(0).vtx_id).vtxtype == GlobalVertex::VTXTYPE::SVTX_MBD){
0442                   h_primvtx[5]->Fill(0);
0443                } 
0444                if((**jet_container_).primaryVertex.at((**jet_container_).recojets.at(JetIDR_map[id]).chConstituents.at(0).vtx_id).vtxtype == GlobalVertex::VTXTYPE::SVTX){
0445                   h_primvtx[5]->Fill(1);
0446                } 
0447                if((**jet_container_).primaryVertex.at((**jet_container_).recojets.at(JetIDR_map[id]).chConstituents.at(0).vtx_id).vtxtype == GlobalVertex::VTXTYPE::MBD){
0448                   h_primvtx[5]->Fill(2);
0449             } */
0450          //}
0451          for(auto con : (**jet_container_).recojets.at(JetIDR_map[id]).chConstituents){
0452             h_pflow[1]->Fill(static_cast<int>(con.pflowtype));
0453             h_phi_track_match->Fill(con.phi);
0454             h_eta_track_match->Fill(con.eta);
0455             h_phi_pt_track->Fill(con.pt,con.phi);
0456          } 
0457          for(auto con : (**jet_container_).recojets.at(JetIDR_map[id]).neConstituents) h_pflow[1]->Fill(static_cast<int>(con.pflowtype));
0458       }
0459 //std::cout<<std::endl<<"Unmatch T ID "<<std::endl;
0460       for(auto& id : unmatchedT){
0461          //std::cout<<id<<" ";
0462          Flavour flavour = getFlavour((**jet_container_).truthjets.at(JetIDT_map[id]));
0463          h_matching[0]->Fill(2.5);
0464          if(flavour.isLF) h_matching[1]->Fill(2.5);
0465          if(flavour.isC & !flavour.isB) h_matching[2]->Fill(2.5);
0466          if(flavour.isB) h_matching[3]->Fill(2.5);
0467       }
0468 //std::cout<<std::endl<<"Doubly R ID "<<std::endl;   
0469       for (auto &id: doublymatchedR) {
0470          for(auto id2 :matchedPair){
0471             if(id2.first == id){
0472                //std::cout<<id2.second<<" ";
0473                Flavour flavour = getFlavour((**jet_container_).truthjets.at(JetIDT_map[id2.second]));
0474                h_matching[0]->Fill(3.5);
0475                if(flavour.isLF) h_matching[1]->Fill(3.5);
0476                if(flavour.isC & !flavour.isB) h_matching[2]->Fill(3.5);
0477                if(flavour.isB) h_matching[3]->Fill(3.5);          
0478             }
0479          }
0480       }
0481 //std::cout<<std::endl<<"Doubly T ID "<<std::endl;  
0482       for (auto &id: doublymatchedT) {
0483          for(auto id2 :matchedPair){
0484             if(id2.second == id){
0485                //std::cout<<id2.second<<" ";
0486                Flavour flavour = getFlavour((**jet_container_).truthjets.at(JetIDT_map[id2.second]));
0487                h_matching[0]->Fill(4.5);
0488                if(flavour.isLF) h_matching[1]->Fill(4.5);
0489                if(flavour.isC & !flavour.isB) h_matching[2]->Fill(4.5);
0490                if(flavour.isB) h_matching[3]->Fill(4.5);          
0491             }
0492          }
0493       }    
0494       //std::cout<<std::endl;
0495 /*
0496          for (auto reco_jet : (**jet_container_).recojets){            
0497             if(std::find(unmatchedR.begin(), unmatchedR.end(),reco_jet.id) != unmatchedR.end()){
0498                for(auto con : reco_jet.chConstituents){
0499                    pflow[1]->Fill(static_cast<int>(con.pflowtype));
0500                   if((**jet_container_).primaryVertex.at(con.vtx_id).vtxtype == GlobalVertex::VTXTYPE::SVTX_MBD) primvtx[5]->Fill(0);
0501                   if((**jet_container_).primaryVertex.at(con.vtx_id).vtxtype == GlobalVertex::VTXTYPE::SVTX) primvtx[5]->Fill(1);
0502                   if((**jet_container_).primaryVertex.at(con.vtx_id).vtxtype == GlobalVertex::VTXTYPE::MBD) primvtx[5]->Fill(2);
0503 
0504                }
0505                for(auto con : reco_jet.neConstituents) pflow[1]->Fill(static_cast<int>(con.pflowtype));
0506                
0507             } 
0508          }*/
0509 
0510 
0511 //std::cout<<"B"<<std::endl;
0512       
0513      
0514       
0515 
0516       /*
0517       for (int ijetR = 0; ijetR < **n_reco_jet_; ijetR++) {
0518          r->Fill((**reco_jet_pt_).at(ijetR));
0519          if(std::find(unmatchedR.begin(), unmatchedR.end(),(**reco_jet_id_).at(ijetR)) != unmatchedR.end()){
0520             ru->Fill((**reco_jet_pt_).at(ijetR));
0521          }
0522 
0523       }
0524          for (int ijetT = 0; ijetT < **n_truth_jet_; ijetT++) { 
0525             t->Fill((**truth_jet_pt_).at(ijetT));
0526          }
0527       */
0528 
0529       //get vertex info
0530       //int num_vtx = 0;
0531       /*for(auto vertex : (**jet_container_).primaryVertex){
0532          num_vtx++;
0533          if(vertex.vtxtype == GlobalVertex::VTXTYPE::SVTX_MBD) h_primvtx[0]->Fill(0);
0534          if(vertex.vtxtype == GlobalVertex::VTXTYPE::SVTX) h_primvtx[0]->Fill(1);
0535          if(vertex.vtxtype == GlobalVertex::VTXTYPE::MBD) h_primvtx[0]->Fill(2);
0536          if(abs(vertex.z) < 10){
0537             if(vertex.vtxtype == GlobalVertex::VTXTYPE::SVTX_MBD) h_primvtx[1]->Fill(0);
0538             if(vertex.vtxtype == GlobalVertex::VTXTYPE::SVTX) h_primvtx[1]->Fill(1);
0539             if(vertex.vtxtype == GlobalVertex::VTXTYPE::MBD) h_primvtx[1]->Fill(2);
0540          }
0541       }*/
0542       //h_n_vtx ->Fill(num_vtx);
0543 
0544 //std::cout<<"C"<<std::endl;
0545       //Main analysis loop
0546       int iid = 0;
0547       for (auto truth_jet : (**jet_container_).truthjets){
0548          Flavour flavour = getFlavour(truth_jet);  
0549          for (auto reco_jet : (**jet_container_).recojets){            
0550    //std::cout<<"D"<<std::endl;         
0551             //find matched combination
0552             if (std::find(matchedPair.begin(), matchedPair.end(), std::pair(reco_jet.id,truth_jet.id)) == matchedPair.end()) continue;
0553             //reject doubly matched reco jets  
0554             if(std::find(doublymatchedR.begin(), doublymatchedR.end(),reco_jet.id) != doublymatchedR.end()){
0555                //for(auto con : reco_jet.chConstituents) h_pflow[2]->Fill(static_cast<int>(con.pflowtype));
0556                //for(auto con : reco_jet.neConstituents) h_pflow[2]->Fill(static_cast<int>(con.pflowtype));
0557                continue;
0558             } 
0559             //reject doubly matched truth jets
0560             if(std::find(doublymatchedT.begin(), doublymatchedT.end(),truth_jet.id) != doublymatchedT.end()){
0561                //for(auto con : reco_jet.chConstituents) h_pflow[3]->Fill(static_cast<int>(con.pflowtype));
0562                //for(auto con : reco_jet.neConstituents) h_pflow[3]->Fill(static_cast<int>(con.pflowtype));
0563                continue;
0564             } 
0565 
0566             //bool good_vertex = false;
0567             //for( auto chtrk : reco_jet.chConstituents){
0568             //   if(abs((**jet_container_).primaryVertex.at(chtrk.vtx_id).z) < 10) good_vertex = true;
0569            // }
0570             //if(!good_vertex) continue;
0571 
0572             
0573 
0574        
0575 
0576           /*  int n_vtx_in_jet = 0;
0577             int vtx_1_id_tmp = -1;
0578             int vtx_2_id_tmp = -1;
0579 
0580             for(auto con : reco_jet.chConstituents){
0581                if(vtx_1_id_tmp == -1){
0582                   vtx_1_id_tmp = con.vtx_id;
0583                   n_vtx_in_jet++;
0584                }
0585                else if ((vtx_1_id_tmp != con.vtx_id) && vtx_2_id_tmp == -1){
0586                   vtx_2_id_tmp = con.vtx_id;
0587                   n_vtx_in_jet++;
0588                }
0589                else if((vtx_1_id_tmp != con.vtx_id) && (vtx_2_id_tmp != con.vtx_id)){
0590                   n_vtx_in_jet++;
0591                }
0592             }
0593 
0594             h_n_vtx_jet->Fill(n_vtx_in_jet);
0595 
0596             if(n_vtx_in_jet != 1) continue;*/
0597 
0598           /*  if((**jet_container_).primaryVertex.at(reco_jet.chConstituents.at(0).vtx_id).vtxtype == GlobalVertex::VTXTYPE::SVTX_MBD){
0599                h_primvtx[2]->Fill(0);
0600                if (flavour.isC & !flavour.isB) h_primvtx[3]->Fill(0);
0601                if (flavour.isB) h_primvtx[4]->Fill(0);
0602             } 
0603             if((**jet_container_).primaryVertex.at(reco_jet.chConstituents.at(0).vtx_id).vtxtype == GlobalVertex::VTXTYPE::SVTX){
0604                h_primvtx[2]->Fill(1);
0605                if (flavour.isC & !flavour.isB) h_primvtx[3]->Fill(1);
0606                if (flavour.isB) h_primvtx[4]->Fill(1);
0607             } 
0608             if((**jet_container_).primaryVertex.at(reco_jet.chConstituents.at(0).vtx_id).vtxtype == GlobalVertex::VTXTYPE::MBD){
0609                h_primvtx[2]->Fill(2);
0610                if (flavour.isC & !flavour.isB) h_primvtx[3]->Fill(2);
0611                if (flavour.isB) h_primvtx[4]->Fill(2);
0612             } */
0613 
0614             h_jet_area->Fill(truth_jet.area, reco_jet.area);
0615 
0616             //for(auto con : reco_jet.chConstituents) h_pflow[0]->Fill(static_cast<int>(con.pflowtype));
0617             //for(auto con : reco_jet.neConstituents) h_pflow[0]->Fill(static_cast<int>(con.pflowtype));
0618 //std::cout<<"E"<<std::endl;
0619             //reset output container
0620             //matched_jets_out.Reset();
0621 
0622 
0623             
0624             vertex_x = (**jet_container_).primaryVertex.x;
0625             vertex_y = (**jet_container_).primaryVertex.y;
0626             vertex_z = (**jet_container_).primaryVertex.z;
0627             vertex_x_unc = (**jet_container_).primaryVertex.x_unc;
0628             vertex_y_unc = (**jet_container_).primaryVertex.y_unc;
0629             vertex_z_unc = (**jet_container_).primaryVertex.z_unc;
0630             vertex_chisq = (**jet_container_).primaryVertex.chisq;
0631             vertex_ndf = (**jet_container_).primaryVertex.ndf;
0632             id = iid;
0633             iid++;
0634 
0635             reco_area = reco_jet.area;
0636             reco_num_ChConstituents = reco_jet.num_ChConstituents;
0637             reco_px = reco_jet.px; 
0638             reco_py = reco_jet.py; 
0639             reco_pz = reco_jet.pz; 
0640             reco_pt = reco_jet.pt; 
0641             reco_eta = reco_jet.eta; 
0642             reco_phi = reco_jet.phi; 
0643             reco_m = reco_jet.m; 
0644             reco_e = reco_jet.e;
0645             truth_area = truth_jet.area;
0646             truth_num_ChConstituents = truth_jet.num_ChConstituents;
0647             truth_px = truth_jet.px; 
0648             truth_py = truth_jet.py; 
0649             truth_pz = truth_jet.pz; 
0650             truth_pt = truth_jet.pt; 
0651             truth_eta = truth_jet.eta; 
0652             truth_phi = truth_jet.phi; 
0653             truth_m = truth_jet.m; 
0654             truth_e = truth_jet.e; 
0655 
0656             Flavour fl = getFlavour(truth_jet);
0657             if(fl.isB) truth_flavour = 2;
0658             else if(fl.isC) truth_flavour = 1;
0659             else truth_flavour = 0;
0660 
0661             track_e.clear();
0662                track_eta.clear();
0663                track_phi.clear();
0664                track_px.clear();
0665                track_py.clear();
0666                track_pz.clear();
0667                track_pt.clear();
0668                track_charge.clear();
0669                track_DCA_xy.clear();
0670                track_DCA_xy_unc.clear();
0671                track_sDCA_xy.clear();
0672                track_DCA3d.clear();
0673                track_sDCA3d.clear();
0674                track_n_mvtx.clear();
0675                track_n_intt.clear();
0676                track_n_tpc.clear();
0677                track_chisq.clear();
0678                track_ndf.clear();
0679 
0680                      
0681 
0682             std::vector<double> sDCA_cut;
0683             
0684             for( auto chtrk : reco_jet.chConstituents){
0685                h_dca_xy->Fill(chtrk.pt,chtrk.DCA_xy);
0686                h_dca_xy_unc->Fill(chtrk.pt,chtrk.DCA_xy_unc);
0687                h_nmvtx->Fill(chtrk.pt,chtrk.n_mvtx);
0688                h_nintt->Fill(chtrk.pt,chtrk.n_intt);
0689                h_ntpc->Fill(chtrk.pt,chtrk.n_tpc);
0690                h_chi2ndof->Fill(chtrk.pt,chtrk.chisq/chtrk.ndf);
0691 
0692 
0693                track_e.push_back(chtrk.e);
0694                track_eta.push_back(chtrk.eta);
0695                track_phi.push_back(chtrk.phi);
0696                track_px.push_back(chtrk.pt * TMath::Cos(chtrk.phi));
0697                track_py.push_back(chtrk.pt * TMath::Sin(chtrk.phi));
0698                track_pz.push_back(chtrk.pt * TMath::SinH(chtrk.eta));
0699                track_pt.push_back(chtrk.pt);
0700                track_charge.push_back(chtrk.charge);
0701                track_DCA_xy.push_back(chtrk.DCA_xy);
0702                track_DCA_xy_unc.push_back(chtrk.DCA_xy_unc);
0703                track_sDCA_xy.push_back(chtrk.sDCA_xy);
0704                track_DCA3d.push_back(chtrk.DCA3d);
0705                track_sDCA3d.push_back(chtrk.sDCA3d);
0706                track_n_mvtx.push_back(chtrk.n_mvtx);
0707                track_n_intt.push_back(chtrk.n_intt);
0708                track_n_tpc.push_back(chtrk.n_tpc);
0709                track_chisq.push_back(chtrk.chisq);
0710                track_ndf.push_back(chtrk.ndf);
0711 
0712                if(abs(chtrk.sDCA3d) < 40){
0713                   //track quality parameters
0714                   if (chtrk.n_mvtx < 3) continue;
0715                   if (chtrk.n_intt < 2) continue;
0716                   if (chtrk.n_tpc < 40) continue;
0717                   if (chtrk.chisq/chtrk.ndf > 3.0) continue;
0718                   if (chtrk.pt < 0.5) continue;
0719                   if (chtrk.DCA_xy_unc > 3.0) continue;
0720 
0721                   sDCA_cut.push_back(chtrk.sDCA3d);  
0722                   //matched_jets_out.reco_constituents_Sdxy.push_back(chtrk.sDCA3d);
0723                   if (flavour.isLF) h_sDCA[0][0]->Fill(chtrk.sDCA3d);
0724                   if (flavour.isC & !flavour.isB) h_sDCA[0][1]->Fill(chtrk.sDCA3d);
0725                   if (flavour.isB) h_sDCA[0][2]->Fill(chtrk.sDCA3d);         
0726                }  
0727             }
0728 
0729             reco_sdxy_1N = -999;
0730             reco_sdxy_2N = -999;
0731             reco_sdxy_3N = -999;
0732 
0733             //sort from greater to lowest
0734             std::nth_element(sDCA_cut.begin(), sDCA_cut.begin()+1, sDCA_cut.end(), std::greater<double>());
0735 
0736             //most displaced track
0737             if(sDCA_cut.size() >= 1){
0738                //matched_jets_out.reco_jet_Sdxy_1N = *(sDCA_cut.begin());
0739                reco_sdxy_1N = *(sDCA_cut.begin());
0740                if (flavour.isLF) h_sDCA[1][0]->Fill(*(sDCA_cut.begin()));
0741                if (flavour.isC & !flavour.isB) h_sDCA[1][1]->Fill(*(sDCA_cut.begin()));
0742                if (flavour.isB) h_sDCA[1][2]->Fill(*(sDCA_cut.begin()));   
0743             } 
0744             //secnd most displaced track
0745             if(sDCA_cut.size() >= 2){
0746                reco_sdxy_2N = *(sDCA_cut.begin()+1);
0747                //matched_jets_out.reco_jet_Sdxy_2N = *(sDCA_cut.begin()+1);
0748                if (flavour.isLF) h_sDCA[1][0]->Fill(*(sDCA_cut.begin()+1));
0749                if (flavour.isC & !flavour.isB) h_sDCA[1][1]->Fill(*(sDCA_cut.begin()+1));
0750                if (flavour.isB) h_sDCA[1][2]->Fill(*(sDCA_cut.begin()+1));   
0751             } 
0752             //third most displaced track
0753             if(sDCA_cut.size() >= 3){
0754                reco_sdxy_3N = *(sDCA_cut.begin()+2);
0755                //matched_jets_out.reco_jet_Sdxy_3N = *(sDCA_cut.begin()+2);
0756                if (flavour.isLF) h_sDCA[1][0]->Fill(*(sDCA_cut.begin()+2));
0757                if (flavour.isC & !flavour.isB) h_sDCA[1][1]->Fill(*(sDCA_cut.begin()+2));
0758                if (flavour.isB) h_sDCA[1][2]->Fill(*(sDCA_cut.begin()+2));   
0759             } 
0760             //fourth most dispalced track
0761             if(sDCA_cut.size() >= 4){
0762                //matched_jets_out.reco_jet_Sdxy_4N = *(sDCA_cut.begin()+3);
0763                if (flavour.isLF) h_sDCA[1][0]->Fill(*(sDCA_cut.begin()));
0764                if (flavour.isC & !flavour.isB) h_sDCA[1][1]->Fill(*(sDCA_cut.begin()+3));
0765                if (flavour.isB) h_sDCA[1][2]->Fill(*(sDCA_cut.begin()+3));   
0766             } 
0767 
0768            /* matched_jets_out.reco_pt = reco_jet.pt;
0769             matched_jets_out.reco_jet_nChConstituents = reco_jet.num_ChConstituents;
0770             matched_jets_out.reco_jet_nConstituents = reco_jet.num_Constituents;
0771             matched_jets_out.truth_jet_flavour = flavour;
0772             matched_jets_out.truth_pt = truth_jet.pt;
0773             matched_jets_out.truth_jet_nChConstituents = truth_jet.num_ChConstituents;
0774             matched_jets_out.truth_jet_nConstituents = truth_jet.num_Constituents;*/
0775 
0776             //ttree_out->Fill();
0777             tree->Fill();
0778          }
0779       }
0780    } // TTree entry / event loop
0781 
0782   /* std::string makeDirectory = "mkdir -p " + work_dir_ + "QA_histo";
0783    system(makeDirectory.c_str());
0784    TString qapath = work_dir_ + "QA_histo/";
0785 
0786    h_reco_all->SaveAs(qapath + h_reco_all->GetName()+".png");
0787    h_truth_all->SaveAs(qapath + h_truth_all->GetName()+".png");
0788    h_reco_unmatch->SaveAs(qapath + h_reco_unmatch->GetName()+".png");
0789    h_reco_match->SaveAs(qapath + h_reco_match->GetName()+".png");
0790    h_truth_match->SaveAs(qapath + h_truth_match->GetName()+".png");
0791    h_phi_all->SaveAs(qapath + h_phi_all->GetName()+".png");
0792    h_phi_match->SaveAs(qapath + h_phi_match->GetName()+".png");
0793    h_eta_all->SaveAs(qapath + h_eta_all->GetName()+".png");
0794    h_eta_match->SaveAs(qapath + h_eta_match->GetName()+".png");
0795    h_phi_truth_all->SaveAs(qapath + h_phi_truth_all->GetName()+".png");
0796    h_phi_truth_match->SaveAs(qapath + h_phi_truth_match->GetName()+".png");
0797    h_eta_truth_all->SaveAs(qapath + h_eta_truth_all->GetName()+".png");
0798    h_eta_truth_match->SaveAs(qapath + h_eta_truth_match->GetName()+".png");
0799    h_dca_xy->SaveAs(qapath + h_dca_xy->GetName()+".png");
0800    h_dca_xy_unc->SaveAs(qapath + h_dca_xy_unc->GetName()+".png");
0801    h_nmvtx->SaveAs(qapath + h_nmvtx->GetName()+".png");
0802    h_nintt->SaveAs(qapath + h_nintt->GetName()+".png");
0803    h_ntpc->SaveAs(qapath + h_ntpc->GetName()+".png");
0804    h_chi2ndof->SaveAs(qapath + h_chi2ndof->GetName()+".png");
0805    for(int i = 0; i<num;i++){ 
0806       for(int j = 0; j<3;j++){ 
0807          h_sDCA[i][j]->SaveAs(qapath + h_sDCA[i][j]->GetName()+".png");
0808       }
0809    }
0810    for(int i = 0; i<4;i++){ 
0811       h_matching[i]->SaveAs(qapath + h_matching[i]->GetName()+".png");
0812    }
0813    for(int i = 0; i<4;i++){ 
0814       h_pflow[i]->SaveAs(qapath + h_pflow[i]->GetName()+".png");
0815    }
0816    for(int i = 0; i<6;i++){ 
0817       h_primvtx[i]->SaveAs(qapath + h_primvtx[i]->GetName()+".png");
0818    }
0819    h_n_vtx->SaveAs(qapath + h_n_vtx->GetName()+".png");
0820    h_n_vtx_jet->SaveAs(qapath + h_n_vtx_jet->GetName()+".png");
0821    h_jet_area->SaveAs(qapath + h_jet_area->GetName()+".png");*/
0822 
0823 
0824    TDirectory *qadir = file->mkdir("QA_histo");
0825    qadir->cd();
0826    h_reco_all->Write();
0827    h_truth_all->Write();
0828    h_reco_unmatch->Write();
0829    h_reco_match->Write();
0830    h_truth_match->Write();
0831    h_phi_all->Write();
0832    h_phi_match->Write();
0833    h_eta_all->Write();
0834    h_eta_match->Write();
0835    h_phi_truth_all->Write();
0836    h_phi_truth_match->Write();
0837    h_eta_truth_all->Write();
0838    h_eta_truth_match->Write();
0839    h_phi_track->Write();
0840    h_eta_track->Write();
0841    h_phi_track_match->Write();
0842    h_eta_track_match->Write();
0843    h_phi_pt_track->Write();
0844    h_dca_xy->Write();
0845    h_dca_xy_unc->Write();
0846    h_nmvtx->Write();
0847    h_nintt->Write();
0848    h_ntpc->Write();
0849    h_chi2ndof->Write();
0850    for(int i = 0; i<num;i++){ 
0851       for(int j = 0; j<3;j++){ 
0852          h_sDCA[i][j]->Write();
0853       }
0854    }
0855    for(int i = 0; i<4;i++){ 
0856       h_matching[i]->Write();
0857    }
0858    for(int i = 0; i<4;i++){ 
0859       h_pflow[i]->Write();
0860    }
0861    for(int i = 0; i<6;i++){ 
0862       h_primvtx[i]->Write();
0863    }
0864    h_n_vtx->Write();
0865    h_n_vtx_jet->Write();
0866    h_jet_area->Write();
0867 
0868    file->cd();
0869    tree->Write("Data", TObject::kOverwrite);
0870    file->Close();
0871    return true;
0872 }
0873 
0874 JetAnalysis::Flavour JetAnalysis::getFlavour(FullJetFinder::TruthJets truth_jet){
0875    Flavour flavour = {false, false, false};
0876    if(!(truth_jet.constituents_PDG_ID.empty())){
0877       for(auto iPDGid : truth_jet.constituents_PDG_ID){
0878          if(TString((TDatabasePDG::Instance()->GetParticle(iPDGid))->ParticleClass()) == "CharmedMeson" || TString((TDatabasePDG::Instance()->GetParticle(iPDGid))->ParticleClass()) == "CharmedBaryon"){
0879             flavour.isC = true;    
0880          }
0881          else if(TString((TDatabasePDG::Instance()->GetParticle(iPDGid))->ParticleClass()) == "B-Meson" || TString((TDatabasePDG::Instance()->GetParticle(iPDGid))->ParticleClass()) == "B-Baryon"){
0882             flavour.isB = true;
0883          }
0884       }
0885    }
0886    else
0887       flavour.isLF = true;
0888 
0889    return flavour;
0890 }
0891 
0892 std::set<int> JetAnalysis::findDuplicates(std::vector<int> vec){
0893    std::set<int> duplicates;
0894    std::sort(vec.begin(), vec.end());
0895    std::set<int> distinct(vec.begin(), vec.end());
0896    std::set_difference(vec.begin(), vec.end(), distinct.begin(), distinct.end(),
0897    std::inserter(duplicates, duplicates.end()));
0898    return duplicates;
0899 }
0900 
0901 void JetAnalysis::printProgress(int cur, int total){  
0902    if((cur % (total/100))==0){
0903       float frac = float(cur)/float(total) + 0.01;
0904       int barWidth = 70;
0905       std::cout << "[";
0906       int pos = barWidth * frac;
0907       for (int i = 0; i < barWidth; ++i) {
0908          if (i < pos) std::cout << "=";
0909          else if (i == pos) std::cout << ">";
0910          else std::cout << " ";
0911       }
0912       std::cout << "] " << int(frac * 100.0) << " %\r";
0913       std::cout.flush();
0914       if(cur >= total*0.99) std::cout<<std::endl;
0915    }
0916 }
0917 
0918 /*
0919 
0920 bool JetAnalysis::scanPurityEfficiency() {
0921    //load TTree produced by processing processCondor
0922    TString filename = work_dir_ + "outtree_v2.root";
0923    TFile *f = new TFile(filename,"READ");
0924    if(!f){
0925       std::cout<<filename<<" does not exist. Run processCondor function first or verify path"<<std::endl;
0926       return false;
0927    }
0928    TTreeReader reader("Data", f);
0929    TTreeReaderValue<JetAnalysis::MatchedJetContainer> *jet_container3 = new TTreeReaderValue<JetAnalysis::MatchedJetContainer>(reader,"matched_jets");
0930 
0931    const int n_points = 250;
0932    //algo power: 1st, 2nd 3rd ; flavour: 0 -b, 1- c, 2- lf, 3- inclusive; number of points in graphs
0933    TH1D *purity[3][4][n_points];
0934    TH1D *efficiency[3][4][n_points];
0935    TH1D *efficiency_flavour[4]; //flavour: 0 -b, 1- c, 2- lf, 3- inclusive
0936    double bins[]{10,120};
0937    //TH2D *RM_all = new TH2D("RM_all","RM_all",200,10,110,200,5,105);
0938    //TH2D *RM_b = new TH2D("RM_b","RM_b",40,10,50,45,5,50);
0939    for (int ipower = 0; ipower < 3; ipower++){
0940       for (int iflavour = 0; iflavour < 4; iflavour++){
0941          for (int ipoint = 0; ipoint < n_points; ipoint++){
0942             TString name = "purity_power"+std::to_string(ipower)+"_flavor_"+std::to_string(iflavour)+"_sdxy_"+std::to_string(ipoint);
0943             purity[ipower][iflavour][ipoint] = new TH1D(name,name,1,bins);
0944          }
0945       }
0946    }
0947    for (int iflavour = 0; iflavour < 4; iflavour++){
0948       TString name = "efficiency_flavour_"+std::to_string(iflavour);
0949       efficiency_flavour[iflavour] = new TH1D(name,name,1,bins);   
0950    }
0951 
0952   //RM_all->GetXaxis()->SetTitle("pT truth");
0953   //RM_all->GetYaxis()->SetTitle("pT reco");
0954   //RM_b->GetXaxis()->SetTitle("pT truth");
0955   //RM_b->GetYaxis()->SetTitle("pT reco");
0956 
0957   
0958 
0959 
0960    //TH1D *resolution[10];
0961    //for (int k = 0; k < 1 ; k++){
0962    //   std::string num_text = std::to_string(bins[k]);
0963    //   std::string rounded = num_text.substr(0, num_text.find(".")+2);
0964    //   std::string num_text2 = std::to_string(bins[k+1]);
0965    //   std::string rounded2 = num_text2.substr(0, num_text2.find(".")+2);
0966    //   TString name = "resolution_pt_truth_"+rounded+"_"+rounded2;
0967    //   resolution[k] = new TH1D(name,name,40,-1,1);
0968    //   resolution[k]->GetXaxis()->SetTitle("dPT = (reco-truth)/truth");
0969    //   resolution[k]->GetYaxis()->SetTitle("probability");
0970    // }
0971   
0972       //  float tag[4]{0,2.5,5,10};
0973 
0974 
0975    TCanvas *canvas;
0976     TPad **Pad;
0977     TH1D **placeholder;
0978   UInt_t nbins = 80+2;
0979   Double_t axis[81+2];
0980   for (int i = 1; i < 82 ; i++){
0981      axis[i] = -40 + i;
0982   }
0983   axis[0] = -45;
0984   axis[82] = 45;
0985 
0986     UInt_t bins2 = 80;
0987   Double_t axis2[81];
0988   for (int i = 0; i < 81 ; i++){
0989      axis2[i] = -40 + i;
0990   }
0991 
0992   TH1D *Signi[3][4];
0993 
0994   for (int i = 0; i < 3 ; i++){
0995     for (int j = 0; j < 4 ; j++){
0996       TString name = "flavour"+std::to_string(i)+"_histo_"+std::to_string(j);
0997       std::cout<<"creating "<<name<<std::endl; 
0998       Signi[i][j] = new TH1D(name,name,bins2,axis2);
0999   }
1000   }
1001 
1002   int c_count = 0;
1003   int b_count = 0;
1004   int lf_count = 0;
1005   int cb_count = 0;
1006 
1007 
1008    while (reader.Next()) {
1009       if (!CheckValue(*jet_container3)) return false;
1010       auto jets = (**jet_container3);
1011 
1012       for (int ipoint = 0; ipoint < n_points ; ipoint++){
1013          if(jets.reco_jet_Sdxy_1N > 0.1*ipoint){
1014             if(jets.truth_jet_flavour.isB) purity[0][0][ipoint]->Fill(jets.reco_pt);
1015             if(jets.truth_jet_flavour.isC && !jets.truth_jet_flavour.isB) purity[0][1][ipoint]->Fill(jets.reco_pt);
1016             if(jets.truth_jet_flavour.isLF) purity[0][2][ipoint]->Fill(jets.reco_pt);
1017             purity[0][3][ipoint]->Fill(jets.reco_pt);
1018          }
1019          if(jets.reco_jet_Sdxy_2N > 0.1*ipoint){
1020             if(jets.truth_jet_flavour.isB) purity[1][0][ipoint]->Fill(jets.reco_pt);
1021             if(jets.truth_jet_flavour.isC && !jets.truth_jet_flavour.isB) purity[1][1][ipoint]->Fill(jets.reco_pt);
1022             if(jets.truth_jet_flavour.isLF) purity[1][2][ipoint]->Fill(jets.reco_pt);
1023             purity[1][3][ipoint]->Fill(jets.reco_pt);
1024             //if(*m_truth_jet_flavour == 2 && *m_reco_jet_Sdxy_1N > tag[1]){
1025             //RM_b->Fill(*m_truth_jet_pt,*m_reco_jet_pt);
1026             //for (int k = 0; k < 9 ; k++){
1027             //if(*m_reco_jet_pt > 5){
1028             //if(bins[k] < *m_truth_jet_pt && *m_truth_jet_pt < bins[k+1])
1029             //resolution[k]->Fill(((*m_reco_jet_pt)-(*m_truth_jet_pt))/(*m_truth_jet_pt));
1030             //}
1031            // }
1032 
1033             //}
1034          }
1035          if(jets.reco_jet_Sdxy_3N > 0.1*ipoint){
1036             if(jets.truth_jet_flavour.isB) purity[2][0][ipoint]->Fill(jets.reco_pt);
1037             if(jets.truth_jet_flavour.isC && !jets.truth_jet_flavour.isB) purity[2][1][ipoint]->Fill(jets.reco_pt);
1038             if(jets.truth_jet_flavour.isLF) purity[2][2][ipoint]->Fill(jets.reco_pt);
1039             purity[2][3][ipoint]->Fill(jets.reco_pt);
1040          }
1041       }
1042 
1043       if(jets.truth_jet_flavour.isB) efficiency_flavour[0]->Fill(jets.reco_pt);
1044       if(jets.truth_jet_flavour.isC && !jets.truth_jet_flavour.isB) efficiency_flavour[1]->Fill(jets.reco_pt);
1045       if(jets.truth_jet_flavour.isLF) efficiency_flavour[2]->Fill(jets.reco_pt);
1046       efficiency_flavour[3]->Fill(jets.reco_pt);
1047       //RM_all->Fill(jets.truth_pt,jets.reco_pt);
1048 
1049       if(jets.truth_jet_flavour.isLF) lf_count++;
1050       if(jets.truth_jet_flavour.isB) b_count++;
1051       if(jets.truth_jet_flavour.isC) c_count++;
1052       if(jets.truth_jet_flavour.isC && jets.truth_jet_flavour.isB) cb_count++;
1053 
1054 
1055       //lf
1056       if(jets.truth_jet_flavour.isLF){
1057         for( auto elem : jets.reco_constituents_Sdxy)Signi[0][0]->Fill(elem);
1058         Signi[0][1]->Fill(jets.reco_jet_Sdxy_1N);
1059         Signi[0][2]->Fill(jets.reco_jet_Sdxy_2N);
1060         Signi[0][3]->Fill(jets.reco_jet_Sdxy_3N);
1061       }
1062       if(jets.truth_jet_flavour.isC && !jets.truth_jet_flavour.isB){
1063         for( auto elem : jets.reco_constituents_Sdxy)Signi[1][0]->Fill(elem);
1064         Signi[1][1]->Fill(jets.reco_jet_Sdxy_1N);
1065         Signi[1][2]->Fill(jets.reco_jet_Sdxy_2N);
1066         Signi[1][3]->Fill(jets.reco_jet_Sdxy_3N);
1067       }
1068       if(jets.truth_jet_flavour.isB){
1069         for( auto elem : jets.reco_constituents_Sdxy)Signi[2][0]->Fill(elem);
1070         Signi[2][1]->Fill(jets.reco_jet_Sdxy_1N);
1071         Signi[2][2]->Fill(jets.reco_jet_Sdxy_2N);
1072         Signi[2][3]->Fill(jets.reco_jet_Sdxy_3N);
1073       }
1074 
1075 
1076    } // end event loop while (reader.Next()) {
1077 
1078       std::cout<<" C: "<<c_count<<" B: "<<b_count<<" LF: "<<lf_count<<" CB: "<<cb_count<<std::endl;
1079 
1080    
1081    for (int ipower = 0; ipower < 3; ipower++){
1082       for (int iflavour = 0; iflavour < 3; iflavour++){
1083          for (int ipoint = 0; ipoint < n_points; ipoint++){
1084             TString name = "eff_power"+std::to_string(ipower)+"_flavor_"+std::to_string(iflavour)+"_sdxy_"+std::to_string(ipoint);
1085             efficiency[ipower][iflavour][ipoint] = (TH1D*)(purity[ipower][iflavour][ipoint]->Clone(name));
1086             //if(!efficiency[ipower][iflavour][ipoint])std::cout<<name<<std::endl;
1087             efficiency[ipower][iflavour][ipoint]->Divide(efficiency[ipower][iflavour][ipoint],efficiency_flavour[iflavour],1,1,"B");
1088             purity[ipower][iflavour][ipoint]->Divide(purity[ipower][iflavour][ipoint],purity[ipower][3][ipoint],1.,1.,"B");
1089          }
1090       }
1091    }
1092 
1093    for (int ipower = 0; ipower < 3; ipower++){
1094       for (int iflavour = 0; iflavour < 3; iflavour++){
1095          for (int ipoint = 0; ipoint < n_points; ipoint++){
1096             if(ipower==0){
1097                purity[ipower][iflavour][ipoint]->SetMarkerColor(kBlue);
1098                purity[ipower][iflavour][ipoint]->SetLineColor(kBlue);
1099                efficiency[ipower][iflavour][ipoint]->SetMarkerColor(kBlue);
1100                efficiency[ipower][iflavour][ipoint]->SetLineColor(kBlue);
1101             }
1102             if(ipower==1){
1103                purity[ipower][iflavour][ipoint]->SetMarkerColor(kGreen);
1104                purity[ipower][iflavour][ipoint]->SetLineColor(kGreen);
1105                efficiency[ipower][iflavour][ipoint]->SetMarkerColor(kGreen);
1106                efficiency[ipower][iflavour][ipoint]->SetLineColor(kGreen);
1107             }
1108             if(ipower==2){
1109                purity[ipower][iflavour][ipoint]->SetMarkerColor(kRed);
1110                purity[ipower][iflavour][ipoint]->SetLineColor(kRed);
1111                efficiency[ipower][iflavour][ipoint]->SetMarkerColor(kRed);
1112                efficiency[ipower][iflavour][ipoint]->SetLineColor(kRed);
1113             }   
1114          }
1115       }
1116    }
1117 
1118 
1119 
1120    double eff_1_b[n_points], pur_1_b[n_points];
1121    double eff_2_b[n_points], pur_2_b[n_points];
1122    double eff_3_b[n_points], pur_3_b[n_points];
1123    double eff_1_c[n_points], pur_1_c[n_points];
1124    double eff_2_c[n_points], pur_2_c[n_points];
1125    double eff_3_c[n_points], pur_3_c[n_points];
1126    double eff_1_lf[n_points], eff_2_lf[n_points], eff_3_lf[n_points];
1127 
1128    for (int ipoint = 0; ipoint < n_points ; ipoint++){
1129       eff_1_b[ipoint] = efficiency[0][0][ipoint]->GetBinContent(1);
1130       eff_2_b[ipoint] = efficiency[1][0][ipoint]->GetBinContent(1);
1131       eff_3_b[ipoint] = efficiency[2][0][ipoint]->GetBinContent(1);
1132       eff_1_c[ipoint] = efficiency[0][1][ipoint]->GetBinContent(1);
1133       eff_2_c[ipoint] = efficiency[1][1][ipoint]->GetBinContent(1);
1134       eff_3_c[ipoint] = efficiency[2][1][ipoint]->GetBinContent(1);
1135       eff_1_lf[ipoint] = efficiency[0][2][ipoint]->GetBinContent(1);
1136       eff_2_lf[ipoint] = efficiency[1][2][ipoint]->GetBinContent(1);
1137       eff_3_lf[ipoint] = efficiency[2][2][ipoint]->GetBinContent(1);
1138 
1139       pur_1_b[ipoint] = purity[0][0][ipoint]->GetBinContent(1);
1140       pur_2_b[ipoint] = purity[1][0][ipoint]->GetBinContent(1);
1141       pur_3_b[ipoint] = purity[2][0][ipoint]->GetBinContent(1);
1142       pur_1_c[ipoint] = purity[0][1][ipoint]->GetBinContent(1);
1143       pur_2_c[ipoint] = purity[1][1][ipoint]->GetBinContent(1);
1144       pur_3_c[ipoint] = purity[2][1][ipoint]->GetBinContent(1);
1145    }
1146 
1147    auto g_1_eff_b_pur_b = new TGraph(n_points,eff_1_b,pur_1_b);
1148    g_1_eff_b_pur_b->SetTitle("Beauty Eff vs Purity;eff_b;pur_b");
1149    g_1_eff_b_pur_b->SetMarkerColor(kGreen);
1150    g_1_eff_b_pur_b->SetLineColor(kGreen);
1151    auto g_2_eff_b_pur_b = new TGraph(n_points,eff_2_b,pur_2_b);
1152    g_2_eff_b_pur_b->SetTitle("Eff vs Purity;eff_b;pur_b");
1153    g_2_eff_b_pur_b->SetMarkerColor(kRed);
1154    g_2_eff_b_pur_b->SetLineColor(kRed);
1155    auto g_3_eff_b_pur_b = new TGraph(n_points,eff_3_b,pur_3_b);
1156    g_3_eff_b_pur_b->SetTitle("Eff vs Purity;eff_b;pur_b");
1157    g_3_eff_b_pur_b->SetMarkerColor(kBlue);
1158    g_3_eff_b_pur_b->SetLineColor(kBlue);
1159 
1160    auto g_1_eff_c_pur_c = new TGraph(n_points,eff_1_c,pur_1_c);
1161    g_1_eff_c_pur_c->SetTitle("Charm Eff vs Purity;eff_c;pur_c");
1162    g_1_eff_c_pur_c->SetMarkerColor(kGreen);
1163    g_1_eff_c_pur_c->SetLineColor(kGreen);
1164    auto g_2_eff_c_pur_c = new TGraph(n_points,eff_2_c,pur_2_c);
1165    g_2_eff_c_pur_c->SetTitle("Eff vs Purity;eff_c;pur_c");
1166    g_2_eff_c_pur_c->SetMarkerColor(kRed);
1167    g_2_eff_c_pur_c->SetLineColor(kRed);
1168    auto g_3_eff_c_pur_c = new TGraph(n_points,eff_3_c,pur_3_c);
1169    g_3_eff_c_pur_c->SetTitle("Eff vs Purity;eff_c;pur_c");
1170    g_3_eff_c_pur_c->SetMarkerColor(kBlue);
1171    g_3_eff_c_pur_c->SetLineColor(kBlue);
1172 
1173    auto g_1_eff_b_eff_c = new TGraph(n_points,eff_1_b,eff_1_c);
1174    g_1_eff_b_eff_c->SetTitle("Eff charm vs Eff beauty;eff_b;eff_c");
1175    g_1_eff_b_eff_c->SetMarkerColor(kGreen);
1176    g_1_eff_b_eff_c->SetLineColor(kGreen);
1177    auto g_2_eff_b_eff_c = new TGraph(n_points,eff_2_b,eff_2_c);
1178    g_2_eff_b_eff_c->SetTitle("Eff vs Eff;eff_b;eff_c");
1179    g_2_eff_b_eff_c->SetMarkerColor(kRed);
1180    g_2_eff_b_eff_c->SetLineColor(kRed);
1181    auto g_3_eff_b_eff_c = new TGraph(n_points,eff_3_b,eff_3_c);
1182    g_3_eff_b_eff_c->SetTitle("Eff vs Eff;eff_b;eff_c");
1183    g_3_eff_b_eff_c->SetMarkerColor(kBlue);
1184    g_3_eff_b_eff_c->SetLineColor(kBlue);
1185 
1186    auto g_1_eff_b_eff_lf = new TGraph(n_points,eff_1_b,eff_1_lf);
1187    g_1_eff_b_eff_lf->SetTitle("Eff LF vs Eff beauty;eff_b;eff_lf");
1188    g_1_eff_b_eff_lf->SetMarkerColor(kGreen);
1189    g_1_eff_b_eff_lf->SetLineColor(kGreen);
1190    auto g_2_eff_b_eff_lf = new TGraph(n_points,eff_2_b,eff_2_lf);
1191    g_2_eff_b_eff_lf->SetTitle("Eff vs Eff;eff_b;eff_lf");
1192    g_2_eff_b_eff_lf->SetMarkerColor(kRed);
1193    g_2_eff_b_eff_lf->SetLineColor(kRed);
1194    auto g_3_eff_b_eff_lf = new TGraph(n_points,eff_3_b,eff_3_lf);
1195    g_3_eff_b_eff_lf->SetTitle("Eff vs Eff;eff_b;eff_lf");
1196    g_3_eff_b_eff_lf->SetMarkerColor(kBlue);
1197    g_3_eff_b_eff_lf->SetLineColor(kBlue);
1198 
1199    TLegend *leg = new TLegend(0.65,0.35,0.9,0.55);
1200    leg->SetHeader("sDCA_xy: 0-25 (0.1 step)","C");
1201    leg->AddEntry(g_1_eff_b_pur_b,"1st most displaced track","lp");
1202    leg->AddEntry(g_2_eff_b_pur_b,"2nd most displaced track","lp");
1203    leg->AddEntry(g_3_eff_b_pur_b,"3rd most displaced track","lp");
1204 
1205 
1206 
1207    //for (int k = 0; k < 9 ; k++){
1208    //   resolution[k]->Scale(1./resolution[k]->Integral());
1209    //   resolution[k]->Scale(1.,"width");
1210    //}
1211 
1212    //   TCanvas *c5 = new TCanvas("c5","c5",2000,1200);
1213    //c5->Divide(6,2);
1214    //c5->cd(1);
1215    //c5->cd(1)->SetLogz();
1216    //RM_all->Draw("colz");
1217    //c5->cd(2);
1218    //c5->cd(2)->SetLogz();
1219    //RM_b->Draw("colz");
1220    //for (int k = 0; k < 9 ; k++){
1221    //   c5->cd(k+3);
1222    //   c5->cd(k+3)->SetLogy();
1223    //   resolution[k]->Draw();
1224    //}
1225 
1226       //   TCanvas *c5 = new TCanvas("c5","c5",2000,1200);
1227    //c5->Divide(6,2);
1228   // c5->cd(1);
1229   // c5->cd(1)->SetLogz();
1230  //  RM_all->Draw("colz");
1231       
1232 std::cout<<"I AM HERE"<<std::endl;
1233 
1234 
1235    TCanvas *c_1 = new TCanvas("c_1","c_1",1600,1200);
1236    c_1->cd();
1237    g_1_eff_b_pur_b->GetXaxis()->SetRangeUser(0,1);
1238    g_1_eff_b_pur_b->GetYaxis()->SetRangeUser(0,1);
1239    g_1_eff_b_pur_b->Draw("AC*");
1240    g_2_eff_b_pur_b->Draw("C* same");
1241    g_3_eff_b_pur_b->Draw("C* same");
1242    leg->Draw();
1243 
1244 
1245    TCanvas *c_2 = new TCanvas("c_2","c_2",1600,1200);
1246    c_2->cd();
1247    g_1_eff_c_pur_c->GetXaxis()->SetRangeUser(0,1);
1248    g_1_eff_c_pur_c->GetYaxis()->SetRangeUser(0,1);
1249    g_1_eff_c_pur_c->Draw("AC*");
1250    g_2_eff_c_pur_c->Draw("C* same");
1251    g_3_eff_c_pur_c->Draw("C* same");
1252    leg->Draw();
1253 
1254    TCanvas *c_3 = new TCanvas("c_3","c_3",1600,1200);
1255    c_3->cd();
1256    c_3->SetLogy();
1257    g_1_eff_b_eff_c->GetXaxis()->SetRangeUser(0,1);
1258    g_1_eff_b_eff_c->GetYaxis()->SetRangeUser(1E-3,1);
1259    g_1_eff_b_eff_c->Draw("AC*");
1260    g_2_eff_b_eff_c->Draw("C* same");
1261    g_3_eff_b_eff_c->Draw("C* same");
1262    leg->Draw();
1263 
1264     TCanvas *c_4 = new TCanvas("c_4","c_4",1600,1200);
1265    c_4->cd();
1266    c_4->SetLogy();
1267    g_1_eff_b_eff_lf->GetXaxis()->SetRangeUser(0,1);
1268    g_1_eff_b_eff_lf->GetYaxis()->SetRangeUser(1E-5,1);
1269    g_1_eff_b_eff_lf->Draw("AC*");
1270    g_2_eff_b_eff_lf->Draw("C* same");
1271    g_3_eff_b_eff_lf->Draw("C* same");
1272    leg->Draw();
1273 
1274 
1275 for (int i = 0; i < 3 ; i++){
1276     for (int j = 0; j < 4 ; j++){
1277       Signi[i][j]->Scale(1./Signi[i][j]->Integral());
1278       Signi[i][j]->Scale(1,"width");
1279       Signi[i][j]->SetLineWidth(1.2);
1280       Signi[0][j]->SetMarkerColor(kBlue);
1281       Signi[1][j]->SetMarkerColor(kGreen);
1282       Signi[2][j]->SetMarkerColor(kRed);
1283       Signi[0][j]->SetLineColor(kBlue);
1284       Signi[1][j]->SetLineColor(kGreen);
1285       Signi[2][j]->SetLineColor(kRed);
1286     }
1287   }
1288 
1289   for (int j = 0; j < 4 ; j++){
1290       Signi[0][j]->SetMarkerStyle(25);
1291       Signi[1][j]->SetMarkerStyle(20);
1292       Signi[2][j]->SetMarkerStyle(21);
1293     }
1294 
1295 
1296 
1297 
1298   std::tie(canvas,Pad,placeholder) =PrepareCanvas4Pad(nbins,axis);
1299 
1300     std::cout<<"A"<<std::endl;
1301 
1302   TLegend *leg6 = new TLegend(.02,.7,.95,.95);
1303   leg6->SetFillStyle(0);
1304   leg6->AddEntry("","#it{#bf{sPHENIX}} Internal","");
1305   leg6->AddEntry("","PYTHIA 8 + GEANT 4","");
1306   leg6->AddEntry("","p+p #sqrt{s} = 200 GeV with no pileup","");
1307   Pad[0]->cd();
1308   leg6->Draw("same");
1309   Signi[0][0]->Draw("same");
1310   Signi[1][0]->Draw("same");
1311   Signi[2][0]->Draw("same");
1312 
1313   std::cout<<"B"<<std::endl;
1314 
1315   TLegend *leg2 = new TLegend(-0.1,.6,.95,.95);
1316   leg2->SetFillStyle(0);
1317   leg2->AddEntry("","Full jets, Anti-#it{k}_{T}","");
1318   leg2->AddEntry("","#it{R}=0.4, |#it{#eta}_{jet}| < 0.7","");
1319   leg2->AddEntry("","#it{p}_{T,jet}^{rec} > 10 GeV/#it{c}","");
1320   leg2->AddEntry("","#it{p}_{T,jet}^{truth} > 30 GeV/#it{c}","");
1321   Pad[1]->cd();
1322   leg2->Draw("same");
1323   Signi[0][1]->Draw("same");
1324   Signi[1][1]->Draw("same");
1325   Signi[2][1]->Draw("same");
1326   TLegend *leg4 = new TLegend(.57,.6,.97,.95);
1327   leg4->SetFillStyle(0);
1328   leg4->AddEntry(Signi[2][1],"b-jet","p");
1329    leg4->AddEntry(Signi[1][1],"c-jet","p");
1330   leg4->AddEntry(Signi[0][1],"light-flavour-jet","p");
1331    leg4->Draw("same");
1332 
1333   std::cout<<"C"<<std::endl;
1334 
1335   TLegend *leg3 = new TLegend(-0.12,.52,.95,.95);
1336   leg3->SetFillStyle(0);
1337   leg3->AddEntry("","Track Quality:","");
1338   leg3->AddEntry("","  #it{p}_{T} > 500 MeV/#it{c}","");
1339   leg3->AddEntry("","  #it{#chi^{2}}/#it{nDOF} < 3","");
1340   leg3->AddEntry("","  #it{Clusters}_{MVTX} > 3","");
1341   leg3->AddEntry("","  #it{Clusters}_{INTT} > 2","");
1342   leg3->AddEntry("","  #it{Clusters}_{TPC} > 40","");
1343   Pad[2]->cd();
1344   leg3->Draw("same");
1345   Signi[0][2]->Draw("same");
1346   Signi[1][2]->Draw("same");
1347   Signi[2][2]->Draw("same");
1348 
1349 std::cout<<"D"<<std::endl;
1350 
1351  
1352   
1353 
1354   Pad[3]->cd();
1355   Signi[0][3]->Draw("same");
1356   Signi[1][3]->Draw("same");
1357   Signi[2][3]->Draw("same");
1358 std::cout<<"E"<<std::endl;
1359   
1360   //TCanvas *c = new TCanvas("c", "c", 685, 630);
1361   //gPad->SetLogz();
1362   //gPad->SetRightMargin(0.15);
1363 
1364   //TLegend *leg = new TLegend(.12,.78,.4,.9);
1365   //leg->SetFillStyle(0);
1366   //leg->AddEntry("","#it{#bf{sPHENIX}} Preliminary","");
1367   //leg->AddEntry("","Au+Au #sqrt{s_{NN}} = 200 GeV","");
1368 
1369   //TH2D *h = new TH2D("", "", 100, 0, 1, 100, 0, 1);
1370   //h->GetYaxis()->SetNdivisions(405);
1371   //h->GetXaxis()->SetNdivisions(405);
1372   //h->SetXTitle("x-axis title [arb. units]");
1373   //h->SetYTitle("y-axis title [arb. units]");
1374   //h->Draw("col z");
1375   //leg->Draw("same");
1376   //TLatex l;
1377   //l.SetNDC();
1378   //l.SetTextFont(43);
1379   //l.SetTextSize(25);
1380   //l.DrawLatex(0.7, 0.965, "#it{7/21/2023}");
1381 
1382   canvas->Print("plot1.pdf");
1383   canvas->Print("plot1.ps");
1384   canvas->Print("plot1.png");
1385    //cc->Divide(10,10);
1386 
1387       //   for (int k = 0; k < 4 ; k++){
1388        //     cc->cd(4*i+k+1);
1389        //     //eff2[i][0][k]->Draw();
1390        //     //eff2[i][1][k]->Draw();
1391        //     eff2[i][2][k]->Draw("e0");
1392 
1393        //  }
1394 
1395 
1396 
1397 //   for (int k = 0; k < 100 ; k++){
1398 //      cc->cd(k+1);
1399 //         eff2[0][2][k]->GetYaxis()->SetRangeUser(0,1);
1400 //         std::string num_text = std::to_string(tag[k]);
1401 //         std::string rounded = num_text.substr(0, num_text.find(".")+2);
1402 //TString name = "efficiency, Sdxy > "+ rounded;//std::to_string(tag[k]);
1403 //         eff2[0][2][k]->SetTitle(name);
1404 //   eff2[0][2][k]->GetYaxis()->SetTitle("b-jet efficiency");
1405 //   eff2[0][2][k]->GetXaxis()->SetTitle("pt jet reco");
1406 //      eff2[0][2][k]->Draw("e0");
1407 //      eff2[1][2][k]->Draw("same e0");
1408 //      eff2[2][2][k]->Draw("same e0");
1409 //   cc->cd(k+5);
1410 //   name = "purity, Sdxy > "+ rounded;
1411 //         purity[0][2][k]->SetTitle(name);
1412 //   purity[0][2][k]->GetYaxis()->SetRangeUser(0,0.15);
1413 //   purity[0][2][k]->GetYaxis()->SetTitle("b-jet purity");
1414 //   purity[0][2][k]->GetXaxis()->SetTitle("pt jet reco");
1415 //   purity[0][2][k]->Draw("e0");
1416 //   purity[1][2][k]->Draw("same e0");
1417 //   purity[2][2][k]->Draw("same e0");
1418   // }
1419 
1420 //   cc->cd(1);
1421 //   purity[0][0]->GetYaxis()->SetRangeUser(0,1);
1422 //   purity[0][0]->Draw();
1423 //   purity[0][1]->Draw("same");
1424 //   purity[0][2]->Draw("same");
1425 //     cc->cd(2);
1426 //   purity[1][0]->GetYaxis()->SetRangeUser(0,1);
1427 //   purity[1][0]->Draw();
1428 //   purity[1][1]->Draw("same");
1429  //  purity[1][2]->Draw("same");
1430  //    cc->cd(2);
1431 //   purity[2][0]->GetYaxis()->SetRangeUser(0,1);
1432 //   purity[2][0]->Draw();
1433 //   purity[2][1]->Draw("same");
1434 //   purity[2][2]->Draw("same");
1435 
1436 //   for (unsigned int ii = 0; ii < 3; ii++){
1437 //     std::cout << ii << std::endl;
1438 //      cc->cd(ii+1);
1439 //      std::cout<<"plotting "<<ii<<" 0"<<std::endl;
1440 //       purity[ii][0]->GetYaxis()->SetRangeUser(0,1);
1441 //      purity[ii][0]->Draw();
1442 //       std::cout<<"plotting "<<ii<<" 1"<<std::endl;
1443 //      purity[ii][1]->Draw("same");
1444 //       std::cout<<"plotting "<<ii<<" 2"<<std::endl;
1445 //      purity[ii][2]->Draw("same");
1446      
1447 
1448 //   }
1449 
1450 
1451 return true;
1452 }
1453 */
1454 /*
1455 void JetAnalysis::printHiisto(TH1* h){
1456    TCanvas c1;
1457    c1.Print("hsimple.ps[");
1458    while ((key = (TKey*)keyList())) {
1459       TClass *cl = gROOT->GetClass(key->GetClassName());
1460       if (!cl->InheritsFrom("TH1")) continue;
1461       TH1 *h = (TH1*)key->ReadObj();
1462       h->Draw();
1463       c1.Print("hsimple.ps");
1464    }
1465    c1.Print("hsimple.ps]");
1466 
1467 }
1468 */
1469 
1470 
1471 std::tuple<TCanvas*, TPad**, TH1D**> PrepareCanvas4Pad(UInt_t xAxisBins, Double_t *xAxis){
1472     //SetsPhenixStyle();
1473     //prepare main canvas
1474     TCanvas *FinalSpectrum = new TCanvas("FinalSpectrum3", "FinalSpectrum3",0,45,(685+3*630),630);
1475     //gStyle->SetOptStat(0);
1476     //gStyle->SetOptTitle(0);
1477     FinalSpectrum->SetHighLightColor(2);
1478     FinalSpectrum->Range(0,0,1,1);
1479     FinalSpectrum->SetFillColor(0);
1480     FinalSpectrum->SetBorderMode(0);
1481     FinalSpectrum->SetBorderSize(2);
1482     FinalSpectrum->SetFrameBorderMode(0);
1483     TPad **pad = new TPad*[static_cast<ULong_t>(4)];
1484     //0  1   2   3
1485 
1486     Double_t marginLeft=0.05;
1487     Double_t marginRight=0.02;
1488     Double_t innerPadWidth=(1-marginLeft-marginRight)/4.;
1489     Double_t marginLeftForXAxis=0.02;
1490 
1491     Double_t marginTop=0.04;
1492     Double_t marginBottom=0.14;
1493     Double_t innerPadHeight=(1-marginTop-marginBottom)/1.;
1494     //Double_t marginBottomForYAxis=0.0;
1495 
1496     //pads for the plots
1497     auto xbaseL = [innerPadWidth,marginLeft,marginLeftForXAxis](int m) ->Double_t {return (m*innerPadWidth+marginLeft-marginLeftForXAxis);};
1498     auto xbaseR = [innerPadWidth,marginLeft](int m) ->Double_t {return (m*innerPadWidth+marginLeft);};
1499     //auto ybaseL = [innerPadHeight,marginBottom,marginBottomForYAxis](int m) ->Double_t {return (m*innerPadHeight+marginBottom-marginBottomForYAxis);};
1500     //auto ybaseR = [innerPadHeight,marginBottom](int m) ->Double_t {return (m*innerPadHeight+marginBottom);};
1501 
1502     //pads for the plots
1503     for(Int_t ipad = 0;ipad < 4;ipad++){
1504         std::cout<<"padID: "<<ipad<<std::endl;
1505         std::cout<<xbaseL(ipad)<<" "<<xbaseR(ipad)<<std::endl;
1506         if(ipad ==0) pad[ipad] = new TPad("pad0", "pad0",0,         0 ,xbaseR(1) ,1);
1507         if(ipad ==1) pad[ipad] = new TPad("pad1", "pad1",xbaseL(1), 0 ,xbaseR(2) ,1);
1508         if(ipad ==2) pad[ipad] = new TPad("pad2", "pad2",xbaseL(2), 0 ,xbaseR(3) ,1);
1509         if(ipad ==3) pad[ipad] = new TPad("pad3", "pad3",xbaseL(3), 0 ,1         ,1);
1510 
1511         pad[ipad]->Draw();
1512         pad[ipad]->cd();
1513         pad[ipad]->SetLogy();
1514 
1515         //left column
1516          if(ipad ==0){
1517             pad[ipad]->SetLeftMargin(static_cast<Float_t>(marginLeft/(marginLeft+innerPadWidth)));
1518             pad[ipad]->SetRightMargin(0);
1519             pad[ipad]->Modified();
1520             pad[ipad]->SetFillStyle(0);
1521         }
1522         //right column
1523         else if (ipad ==1 || ipad ==2){
1524             pad[ipad]->SetLeftMargin(static_cast<Float_t>(marginLeftForXAxis/(innerPadWidth+marginLeftForXAxis)));
1525             pad[ipad]->SetRightMargin(0);
1526             pad[ipad]->Modified();
1527             pad[ipad]->SetFillStyle(0);
1528         }
1529         else if (ipad ==3){
1530             pad[ipad]->SetLeftMargin(static_cast<Float_t>(marginLeftForXAxis/(innerPadWidth+marginLeftForXAxis+marginRight)));
1531             pad[ipad]->SetRightMargin(static_cast<Float_t>(marginRight/(innerPadWidth+marginRight)));
1532             pad[ipad]->Modified();
1533             pad[ipad]->SetFillStyle(0);
1534         }
1535 
1536         pad[ipad]->SetBottomMargin(marginBottom/(marginBottom+innerPadHeight));
1537         pad[ipad]->SetTopMargin(marginTop/(2*innerPadHeight+marginTop));
1538 
1539         pad[ipad]->SetFillStyle(0);
1540         pad[ipad]->Modified();
1541         pad[ipad]->Update();
1542         FinalSpectrum->cd();
1543     }
1544 
1545     //ratio Y-axis title placeholder
1546     TH1D **placeholder = new TH1D*[static_cast<ULong_t>(4)];
1547 
1548 
1549     //other axis placeholders
1550     for(Int_t ipad = 0;ipad <4;ipad++){
1551             pad[ipad]->cd();
1552             placeholder[ipad] = new TH1D(Form("placeholder2%d",ipad),"Central Values",static_cast<Int_t>(xAxisBins), xAxis);
1553             //set ranges
1554     
1555             placeholder[ipad]->SetMinimum(1E-5);
1556             placeholder[ipad]->SetMaximum(100);
1557 
1558             //if(ipad>3)placeholder[ipad]->GetXaxis()->SetTitle("#it{p}_{T,ch. jet} (GeV/#it{c})");
1559             //if(ipad<2)placeholder[ipad]->GetYaxis()->SetTitle("#frac{d^{2}#it{#sigma}^{#it{R}=0.2}}{d#it{p}_{T,ch. jet}d#it{#eta}_{ch. jet}}/#frac{d^{2}#it{#sigma}^{#it{R}=#it{X}}}{d#it{p}_{T,ch. jet}d#it{#eta}_{ch. jet}}");
1560             //else if(ipad==4)placeholder[ipad]->GetYaxis()->SetTitle("MC/Data");
1561 
1562             //TMathText l;
1563             //l.
1564             TString tt = "\\mathscr{S} DCA_{xy}";
1565             //tt + "\\mathscr{S}";
1566             //tt + "DCA_{xy}";
1567 
1568             if(ipad==0)placeholder[ipad]->GetYaxis()->SetTitle("Probability");
1569             if(ipad==0)placeholder[ipad]->GetXaxis()->SetTitle("#it{S}DCA_{xy}");
1570             if(ipad==1)placeholder[ipad]->GetXaxis()->SetTitle("1^{ st} largest #it{S}DCA_{xy}");
1571             if(ipad==2)placeholder[ipad]->GetXaxis()->SetTitle("2^{nd} largest #it{S}DCA_{xy}");
1572             if(ipad==3)placeholder[ipad]->GetXaxis()->SetTitle("3^{rd} largest #it{S}DCA_{xy}");
1573 
1574             /*placeholder[ipad]->GetYaxis()->SetTitleFont(43);
1575             placeholder[ipad]->GetXaxis()->SetTitleFont(43);
1576             placeholder[ipad]->GetYaxis()->SetLabelFont(43);
1577             placeholder[ipad]->GetXaxis()->SetLabelFont(43);
1578             placeholder[ipad]->GetYaxis()->SetTitleSize(0.05);
1579             placeholder[ipad]->GetXaxis()->SetTitleSize(0.05);
1580             placeholder[ipad]->GetYaxis()->SetLabelSize(0.05);
1581             placeholder[ipad]->GetXaxis()->SetLabelSize(0.05);*/
1582 
1583             placeholder[ipad]->GetYaxis()->SetTitleOffset(1.2f);//0.95*(gPad->GetHNDC())/scaleHeightPads/resizeTextFactor);
1584             placeholder[ipad]->GetXaxis()->SetTitleOffset(1.2f);// not number before, no number was optimized before font 43 was introduced
1585             if(ipad>0){
1586             placeholder[ipad]->GetYaxis()->SetTitleSize(0);
1587             placeholder[ipad]->GetYaxis()->SetLabelSize(0);
1588             }
1589 
1590             placeholder[ipad]->GetXaxis()->SetNdivisions(405);
1591             placeholder[ipad]->GetYaxis()->SetNdivisions(405);
1592             /*else placeholder[ipad]->GetYaxis()->SetDecimals();
1593             if(ipad>1)placeholder[ipad]->GetYaxis()->SetNdivisions(5, 5, 0, kTRUE);
1594             //because root is so advanced it decided to put different ticks lenght for 2 particular plots - lets just force it to look same
1595             if(ipad==4 || ipad==5) placeholder[ipad]->GetYaxis()->SetTickLength(0.04f);
1596             if(ipad>=2) placeholder[ipad]->GetXaxis()->SetTickLength(0.06);*/
1597 
1598             placeholder[ipad]->Draw();
1599         }
1600 
1601     return std::make_tuple(FinalSpectrum,pad, placeholder);
1602 
1603 }
[ 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!