FullJetFinder/macro/readJetTree.C

0001 /**
0002  * @file readJetTree.C
0003  *
0004  * @brief An example macro how to proccess TTrees from FullJetFinder.cc
0005  *
0006  * @ingroup FullJetFinder
0007  *
0008  * @author Jakub Kvapil
0009  * Contact: jakub.kvapil@cern.ch
0010  *
0011  */
0012
0013 #include <TFile.h>
0014 #include <TH1.h>
0015 #include <TH2.h>
0016 #include <TCanvas.h>
0017 #include <TTreeReader.h>
0018 #include <TTreeReaderValue.h>
0019 #include <TTreeReaderArray.h>
0020 #include <vector>
0021 #include <iostream>
0022 #include <TDatabasePDG.h>
0023 #include <TString.h>
0024 #include <set>
0025 #include <TChain.h>
0026 #include <TGraph.h>
0027 #include <string>
0028 #include <vector>
0029
0030 #pragma GCC diagnostic push
0031 #pragma GCC diagnostic ignored "-Wundefined-internal"
0032 #include <FullJetFinder.h>
0033 #pragma GCC diagnostic pop
0034
0035 R__LOAD_LIBRARY(libFullJetFinder.so)
0036
0037 bool CheckValue(ROOT::Internal::TTreeReaderValueBase& value) {
0038    if (value.GetSetupStatus() < 0) {
0039       std::cerr << "Error " << value.GetSetupStatus()
0040                 << "setting up reader for " << value.GetBranchName() << '\n';
0041       return false;
0042    }
0043    return true;
0044 }
0045
0046 void printProgress(int cur, int total){
0047    if((cur % (total/100))==0){
0048       float frac = float(cur)/float(total) + 0.01;
0049       int barWidth = 70;
0050       std::cout << "[";
0051       int pos = barWidth * frac;
0052       for (int i = 0; i < barWidth; ++i) {
0053          if (i < pos) std::cout << "=";
0054          else if (i == pos) std::cout << ">";
0055          else std::cout << " ";
0056       }
0057       std::cout << "] " << int(frac * 100.0) << " %\r";
0058       std::cout.flush();
0059       if(cur >= total*0.99) std::cout<<std::endl;
0060    }
0061 }
0062
0063 // Analyze TTree "AntiKt_r04/Data" in the file passed into the function.
0064 // Returns false in case of errors.
0065 //takes list of datafiles as entry
0066 bool readJetTree(const std::string &dataFiles = "/sphenix/tg/tg01/hf/jkvapil/JET10_r22_npile_full/condorJob/myTestJets/outputData_000*.root"){
0067    TChain *tc = new TChain("AntiKt_r04/Data"); //TTRee name
0068    tc->Add(dataFiles.data());
0069    int n_events = tc->GetEntries(); //gen total number of events
0070    tc->LoadTree(-1); //previous commands detached the Tree head, reset back to supress warning
0071    TTreeReader reader(tc);
0072    TTreeReaderValue<FullJetFinder::Container> jet_container(reader,"JetContainer");
0073
0074    TH1I *h_jet_n = new TH1I("h_jet_n","h_jet_n",2,-0.5,1.5);
0075    TH1F *h_reco_jet_pt = new TH1F("h_reco_jet_pt","h_reco_jet_pt",100,0,100);
0076    TH1F *h_reco_track_pt = new TH1F("h_reco_track_pt","h_reco_track_pt",100,0,100);
0077    TH1F *h_dca3d = new TH1F("h_dca3d","h_dca3d",100,-10,10);
0078
0079    std::cout<<"Total number of events: "<<n_events<<std::endl;
0080
0081    //main event loop
0082    while (reader.Next()){
0083       if (!CheckValue(jet_container)) return false; //NULL guards
0084       printProgress(reader.GetCurrentEntry(),n_events);
0085
0086       h_jet_n->Fill(0.0,jet_container->reco_jet_n);
0087       h_jet_n->Fill(1.0,jet_container->truth_jet_n);
0088
0089       //reco jet loop
0090       for (auto reco_jet : jet_container->recojets){
0091          h_reco_jet_pt->Fill(reco_jet.pt);
0092          //track inside jet loop
0093          for (auto chConstituent : reco_jet.chConstituents){
0094             h_reco_track_pt->Fill(chConstituent.pt);
0095        h_dca3d->Fill(chConstituent.sDCA3d);
0096          }  // for (auto chConstituent : reco_jet.chConstituents)
0097       }//f or (auto reco_jet : jet_container->recojets)
0098    } //while (reader.Next())
0099
0100    //plotting
0101    TCanvas *c = new TCanvas("c","c",2400,800);
0102    c->Divide(4,1);
0103    c->cd(1);
0104    h_jet_n->Draw();
0105    c->cd(2);
0106    h_reco_jet_pt->Draw();
0107    c->cd(3);
0108    h_reco_track_pt->Draw();
0109 c->cd(4);
0110 h_dca3d->Draw();
0111    return true;
0112 }