Fork me on GitHub

Changes in / [2264876:7993cad] in git


Ignore:
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r2264876 r7993cad  
    137137tmp/examples/Example1.$(ObjSuf): \
    138138        examples/Example1.cpp \
    139         classes/DelphesClasses.h \
    140         external/ExRootAnalysis/ExRootTreeReader.h \
    141         external/ExRootAnalysis/ExRootTreeWriter.h \
    142         external/ExRootAnalysis/ExRootTreeBranch.h \
    143         external/ExRootAnalysis/ExRootResult.h \
    144         external/ExRootAnalysis/ExRootUtilities.h
    145 Validation$(ExeSuf): \
    146         tmp/examples/Validation.$(ObjSuf)
    147 
    148 tmp/examples/Validation.$(ObjSuf): \
    149         examples/Validation.cpp \
    150139        classes/DelphesClasses.h \
    151140        external/ExRootAnalysis/ExRootTreeReader.h \
    … …  
    161150        root2pileup$(ExeSuf) \
    162151        stdhep2pileup$(ExeSuf) \
    163         Example1$(ExeSuf) \
    164         Validation$(ExeSuf)
     152        Example1$(ExeSuf)
    165153
    166154EXECUTABLE_OBJ +=  \
    … …  
    171159        tmp/converters/root2pileup.$(ObjSuf) \
    172160        tmp/converters/stdhep2pileup.$(ObjSuf) \
    173         tmp/examples/Example1.$(ObjSuf) \
    174         tmp/examples/Validation.$(ObjSuf)
     161        tmp/examples/Example1.$(ObjSuf)
    175162
    176163DelphesHepMC$(ExeSuf): \
    … …  
    196183        classes/DelphesLHEFReader.h \
    197184        external/ExRootAnalysis/ExRootTreeWriter.h \
    198         external/ExRootAnalysis/ExRootTreeBranch.h \
    199         external/ExRootAnalysis/ExRootProgressBar.h
    200 DelphesROOT$(ExeSuf): \
    201         tmp/readers/DelphesROOT.$(ObjSuf)
    202 
    203 tmp/readers/DelphesROOT.$(ObjSuf): \
    204         readers/DelphesROOT.cpp \
    205         modules/Delphes.h \
    206         classes/DelphesStream.h \
    207         classes/DelphesClasses.h \
    208         classes/DelphesFactory.h \
    209         external/ExRootAnalysis/ExRootTreeWriter.h \
    210         external/ExRootAnalysis/ExRootTreeReader.h \
    211185        external/ExRootAnalysis/ExRootTreeBranch.h \
    212186        external/ExRootAnalysis/ExRootProgressBar.h
    … …  
    226200        DelphesHepMC$(ExeSuf) \
    227201        DelphesLHEF$(ExeSuf) \
    228         DelphesROOT$(ExeSuf) \
    229202        DelphesSTDHEP$(ExeSuf)
    230203
    … …  
    232205        tmp/readers/DelphesHepMC.$(ObjSuf) \
    233206        tmp/readers/DelphesLHEF.$(ObjSuf) \
    234         tmp/readers/DelphesROOT.$(ObjSuf) \
    235207        tmp/readers/DelphesSTDHEP.$(ObjSuf)
    236208
    … …  
    372344        modules/PdgCodeFilter.h \
    373345        modules/BeamSpotFilter.h \
    374         modules/RecoPuFilter.h \
    375346        modules/Cloner.h \
    376347        modules/Weighter.h \
    … …  
    380351        modules/VertexSorter.h \
    381352        modules/VertexFinder.h \
    382         modules/VertexFinderDA4D.h \
    383353        modules/ExampleModule.h
    384354tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    828798        external/ExRootAnalysis/ExRootFilter.h \
    829799        external/ExRootAnalysis/ExRootClassifier.h
    830 tmp/modules/RecoPuFilter.$(ObjSuf): \
    831         modules/RecoPuFilter.$(SrcSuf) \
    832         modules/RecoPuFilter.h \
    833         classes/DelphesClasses.h \
    834         classes/DelphesFactory.h \
    835         classes/DelphesFormula.h \
    836         external/ExRootAnalysis/ExRootResult.h \
    837         external/ExRootAnalysis/ExRootFilter.h \
    838         external/ExRootAnalysis/ExRootClassifier.h
    839800tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    840801        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    935896        modules/VertexFinder.$(SrcSuf) \
    936897        modules/VertexFinder.h \
    937         classes/DelphesClasses.h \
    938         classes/DelphesFactory.h \
    939         classes/DelphesFormula.h \
    940         classes/DelphesPileUpReader.h \
    941         external/ExRootAnalysis/ExRootResult.h \
    942         external/ExRootAnalysis/ExRootFilter.h \
    943         external/ExRootAnalysis/ExRootClassifier.h
    944 tmp/modules/VertexFinderDA4D.$(ObjSuf): \
    945         modules/VertexFinderDA4D.$(SrcSuf) \
    946         modules/VertexFinderDA4D.h \
    947898        classes/DelphesClasses.h \
    948899        classes/DelphesFactory.h \
    … …  
    1047998        tmp/modules/PileUpJetID.$(ObjSuf) \
    1048999        tmp/modules/PileUpMerger.$(ObjSuf) \
    1049         tmp/modules/RecoPuFilter.$(ObjSuf) \
    10501000        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    10511001        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    10601010        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
    10611011        tmp/modules/VertexFinder.$(ObjSuf) \
    1062         tmp/modules/VertexFinderDA4D.$(ObjSuf) \
    10631012        tmp/modules/VertexSorter.$(ObjSuf) \
    10641013        tmp/modules/Weighter.$(ObjSuf)
    … …  
    16661615        tmp/external/tcl/tclVar.$(ObjSuf)
    16671616
    1668 modules/VertexFinderDA4D.h: \
    1669         classes/DelphesModule.h \
    1670         classes/DelphesClasses.h
    1671         @touch $@
    1672 
    16731617modules/TrackSmearing.h: \
    16741618        classes/DelphesModule.h
    … …  
    20151959
    20161960modules/BTagging.h: \
    2017         classes/DelphesModule.h
    2018         @touch $@
    2019 
    2020 modules/RecoPuFilter.h: \
    20211961        classes/DelphesModule.h
    20221962        @touch $@

  • examples/Validation.cpp

    r2264876 r7993cad  
    2121#include <utility>
    2222#include <vector>
    23 #include <typeinfo>
    2423
    2524#include "TROOT.h"
    … …  
    2928#include "TString.h"
    3029
    31 #include "TH1.h"
    3230#include "TH2.h"
    33 #include "TMath.h"
    34 #include "TStyle.h"
    35 #include "TGraph.h"
    36 #include "TCanvas.h"
    3731#include "THStack.h"
    3832#include "TLegend.h"
    … …  
    4034#include "TClonesArray.h"
    4135#include "TLorentzVector.h"
    42 #include "TGraphErrors.h"
    43 #include "TMultiGraph.h"
    4436
    4537#include "classes/DelphesClasses.h"
    … …  
    5547//------------------------------------------------------------------------------
    5648
    57 double ptrangemin = 10;
    58 double ptrangemax = 10000;
    59 static const int Nbins = 20;
     49// Here you can put your analysis macro
    6050
    61 int objStyle = 1;
    62 int trackStyle = 7;
    63 int towerStyle = 3;
    64 
    65 Color_t objColor = kBlack;
    66 Color_t trackColor = kBlack;
    67 Color_t towerColor = kBlack;
    68 
    69 double effLegXmin = 0.22;
    70 double effLegXmax = 0.7;
    71 double effLegYmin = 0.22;
    72 double effLegYmax = 0.5;
    73 
    74 double resLegXmin = 0.62;
    75 double resLegXmax = 0.9;
    76 double resLegYmin = 0.52;
    77 double resLegYmax = 0.85;
    78 
    79 double topLeftLegXmin = 0.22;
    80 double topLeftLegXmax = 0.7;
    81 double topLeftLegYmin = 0.52;
    82 double topLeftLegYmax = 0.85;
    83 
    84 
    85 struct resolPlot
    86 {
    87   TH1 *cenResolHist;
    88   TH1 *fwdResolHist;
    89   double ptmin;
    90   double ptmax;
    91   double xmin;
    92   double xmax;
    93   TString obj;
    94 
    95   resolPlot();
    96   resolPlot(double ptdown, double ptup, TString object);
    97   void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
    98   void print(){std::cout << ptmin << std::endl;}
    99 };
    100 
    101 
    102 resolPlot::resolPlot()
    103 {
    104 }
    105 
    106 resolPlot::resolPlot(double ptdown, double ptup, TString object)
    107 {
    108   this->set(ptdown,ptup,object);
    109 }
    110 
    111 void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax)
    112 {
    113   ptmin = ptdown;
    114   ptmax = ptup;
    115   obj = object;
    116 
    117   cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 200,  xmin, xmax);
    118   fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 200,  xmin, xmax);
    119 }
    120 
    121 void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
    122 {
    123   double width;
    124   double ptdown;
    125   double ptup;
    126   resolPlot ptemp;
    127 
    128   for (int i = 0; i < Nbins; i++)
    129   {
    130     width = (ptmax - ptmin) / Nbins;
    131     ptdown = TMath::Power(10,ptmin + i * width );
    132     ptup = TMath::Power(10,ptmin + (i+1) * width );
    133     ptemp.set(ptdown, ptup, obj, xmin, xmax);
    134     histos->push_back(ptemp);
    135   }
    136 }
    137 
    138 //------------------------------------------------------------------------------
    139 
    140 class ExRootResult;
    141 class ExRootTreeReader;
    142 
    143 //------------------------------------------------------------------------------
    144 
    145 void BinLogX(TH1*h)
    146 {
    147   TAxis *axis = h->GetXaxis();
    148   int bins = axis->GetNbins();
    149 
    150   Axis_t from = axis->GetXmin();
    151   Axis_t to = axis->GetXmax();
    152   Axis_t width = (to - from) / bins;
    153   Axis_t *new_bins = new Axis_t[bins + 1];
    154 
    155   for (int i = 0; i <= bins; i++)
    156   {
    157     new_bins[i] = TMath::Power(10, from + i * width);
    158   }
    159   axis->Set(bins, new_bins);
    160   delete new_bins;
    161 }
    162 
    163 
    164 //------------------------------------------------------------------------------
    165 
    166 template<typename T>
    167 std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
    168 {
    169 
    170   cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
    171 
    172   Long64_t allEntries = treeReader->GetEntries();
    173 
    174   GenParticle *particle;
    175   T *recoObj;
    176 
    177   TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
    178 
    179   Float_t deltaR;
    180   Float_t pt, eta;
    181   Long64_t entry;
    182 
    183   Int_t i, j;
    184 
    185   TH1D *histGenPtcen = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
    186   TH1D *histRecoPtcen = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
    187   TH1D *histGenPtfwd  = new TH1D(name+" gen spectra Eta",name+" gen spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
    188   TH1D *histRecoPtfwd = new TH1D(name+" reco spectra Eta",name+" reco spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
    189 
    190 
    191   BinLogX(histGenPtcen);
    192   BinLogX(histRecoPtcen);
    193   BinLogX(histGenPtfwd);
    194   BinLogX(histRecoPtfwd);
    195 
    196   // Loop over all events
    197   for(entry = 0; entry < allEntries; ++entry)
    198   {
    199     // Load selected branches with data from specified event
    200     treeReader->ReadEntry(entry);
    201 
    202     // Loop over all generated particle in event
    203     for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
    204     {
    205 
    206       particle = (GenParticle*) branchParticle->At(i);
    207       genMomentum = particle->P4();
    208 
    209       deltaR = 999;
    210 
    211       if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
    212       {
    213 
    214         // Loop over all reco object in event
    215         for(j = 0; j < branchReco->GetEntriesFast(); ++j)
    216         {
    217           recoObj = (T*)branchReco->At(j);
    218           recoMomentum = recoObj->P4();
    219           // this is simply to avoid warnings from initial state particle
    220           // having infite rapidity ...
    221         //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
    222 
    223           // take the closest parton candidate
    224           if(TMath::Abs(pdgID) == 5)
    225           {
    226             Jet *jet = (Jet *)recoObj;
    227             if(jet->BTag != 1) continue;
    228           }
    229           if(TMath::Abs(pdgID) == 15)
    230           {
    231             Jet *jet = (Jet *)recoObj;
    232             if(jet->TauTag != 1) continue;
    233           }
    234           if(genMomentum.DeltaR(recoMomentum) < deltaR)
    235           {
    236             deltaR = genMomentum.DeltaR(recoMomentum);
    237             bestRecoMomentum = recoMomentum;
    238           }
    239         }
    240 
    241         pt  = genMomentum.Pt();
    242         eta = genMomentum.Eta();
    243 
    244         if (TMath::Abs(eta) < 1.5)
    245         {
    246           histGenPtcen->Fill(pt);
    247           if(deltaR < 0.3) { histRecoPtcen->Fill(pt); }
    248         }
    249         else if (TMath::Abs(eta) < 2.5)
    250         {
    251           histGenPtfwd->Fill(pt);
    252           if(deltaR < 0.3) { histRecoPtfwd->Fill(pt); }
    253 
    254         }
    255       }
    256     }
    257   }
    258 
    259 
    260   std::pair<TH1D*,TH1D*> histos;
    261 
    262   histRecoPtcen->Divide(histGenPtcen);
    263   histRecoPtfwd->Divide(histGenPtfwd);
    264 
    265   histos.first = histRecoPtcen;
    266   histos.second = histRecoPtfwd;
    267 
    268   return histos;
    269 }
    270 
    271 template<typename T>
    272 void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
    273 {
    274   Long64_t allEntries = treeReader->GetEntries();
    275 
    276   cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
    277 
    278   GenParticle *particle;
    279   T* recoObj;
    280 
    281   TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
    282 
    283   Float_t deltaR;
    284   Float_t pt, eta;
    285   Long64_t entry;
    286 
    287   Int_t i, j, bin;
    288 
    289   // Loop over all events
    290   for(entry = 0; entry < allEntries; ++entry)
    291   {
    292     // Load selected branches with data from specified event
    293     treeReader->ReadEntry(entry);
    294 
    295     // Loop over all reconstructed jets in event
    296     for(i = 0; i < branchReco->GetEntriesFast(); ++i)
    297     {
    298       recoObj = (T*) branchReco->At(i);
    299       recoMomentum = recoObj->P4();
    300 
    301       deltaR = 999;
    302 
    303      // Loop over all hard partons in event
    304      for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
    305      {
    306         particle = (GenParticle*) branchParticle->At(j);
    307         if (particle->PID == pdgID && particle->Status == 1)
    308         {
    309           genMomentum = particle->P4();
    310 
    311           // this is simply to avoid warnings from initial state particle
    312           // having infite rapidity ...
    313           if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
    314 
    315           // take the closest parton candidate
    316           if(genMomentum.DeltaR(recoMomentum) < deltaR)
    317           {
    318              deltaR = genMomentum.DeltaR(recoMomentum);
    319              bestGenMomentum = genMomentum;
    320           }
    321         }
    322       }
    323 
    324       if(deltaR < 0.3)
    325       {
    326         pt  = bestGenMomentum.Pt();
    327         eta = TMath::Abs(bestGenMomentum.Eta());
    328 
    329         for (bin = 0; bin < Nbins; bin++)
    330         {
    331           if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
    332           {
    333             if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
    334             else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
    335           }
    336         }
    337       }
    338     }
    339   }
    340 }
    341 void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
    342 {
    343 
    344   Long64_t allEntries = treeReader->GetEntries();
    345 
    346   cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
    347 
    348   Jet *jet, *genjet;
    349 
    350   TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
    351 
    352   Float_t deltaR;
    353   Float_t pt, eta;
    354   Long64_t entry;
    355 
    356   Int_t i, j, bin;
    357 
    358   // Loop over all events
    359   for(entry = 0; entry < allEntries; ++entry)
    360   {
    361     // Load selected branches with data from specified event
    362     treeReader->ReadEntry(entry);
    363 
    364     if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
    365 
    366     // Loop over all reconstructed jets in event
    367     for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
    368     {
    369 
    370       jet = (Jet*) branchJet->At(i);
    371       jetMomentum = jet->P4();
    372 
    373       deltaR = 999;
    374 
    375      // Loop over all hard partons in event
    376      for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
    377      {
    378         genjet = (Jet*) branchGenJet->At(j);
    379 
    380         genJetMomentum = genjet->P4();
    381 
    382         // this is simply to avoid warnings from initial state particle
    383         // having infite rapidity ...
    384         if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
    385 
    386         // take the closest parton candidate
    387         if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
    388         {
    389            deltaR = genJetMomentum.DeltaR(jetMomentum);
    390            bestGenJetMomentum = genJetMomentum;
    391         }
    392       }
    393 
    394       if(deltaR < 0.25)
    395       {
    396         pt  = genJetMomentum.Pt();
    397         eta = TMath::Abs(genJetMomentum.Eta());
    398 
    399         for (bin = 0; bin < Nbins; bin++)
    400         {
    401             if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
    402             {
    403                 histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
    404             }
    405             else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
    406             {
    407                 histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
    408             }
    409 
    410         }
    411       }
    412     }
    413   }
    414 }
    415 
    416 void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
    417 {
    418 
    419   Long64_t allEntries = treeReader->GetEntries();
    420 
    421   cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
    422 
    423   MissingET *met;
    424   ScalarHT *scalarHT;
    425 
    426   Long64_t entry;
    427 
    428   Int_t bin;
    429   Double_t ht;
    430 
    431   Jet *jet;
    432   TLorentzVector p1, p2;
    433 
    434   // Loop over all events
    435   for(entry = 0; entry < allEntries; ++entry)
    436   {
    437     // Load selected branches with data from specified event
    438     treeReader->ReadEntry(entry);
    439 
    440     if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
    441 
    442     if (branchJet->GetEntriesFast() > 1)
    443     {
    444 
    445       jet = (Jet*) branchJet->At(0);
    446       p1 = jet->P4();
    447       jet = (Jet*) branchJet->At(1);
    448       p2 = jet->P4();
    449 
    450       met = (MissingET*) branchMet->At(0);
    451       scalarHT = (ScalarHT*) branchScalarHT->At(0);
    452       ht = scalarHT->HT;
    453 
    454       if(p1.Pt() < 0.75*ht/2) continue;
    455       if(p2.Pt() < 0.75*ht/2) continue;
    456 
    457       for (bin = 0; bin < Nbins; bin++)
    458       {
    459         if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax )
    460         {
    461           histos->at(bin).cenResolHist->Fill(met->P4().Px());
    462           histos->at(bin).fwdResolHist->Fill(met->P4().Py());
    463         }
    464       }
    465     }
    466   }
    467 }
    468 
    469 
    470 std::pair<Double_t, Double_t> GausFit(TH1* hist)
    471 {
    472   TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
    473   hist->Fit("f1","RQ");
    474   Double_t sig = f1->GetParameter(2);
    475   Double_t sigErr = f1->GetParError(2);
    476   delete f1;
    477   return make_pair (sig, sigErr);
    478 }
    479 
    480 
    481 TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
    482 {
    483   Int_t bin;
    484   Int_t count = 0;
    485   TGraphErrors gr = TGraphErrors(Nbins/2);
    486   Double_t sig = 0;
    487   Double_t sigErr = 0;
    488   for (bin = 0; bin < Nbins; bin++)
    489   {
    490     if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
    491     {
    492       std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
    493       if (rms == true)
    494       {
    495         gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
    496         gr.SetPointError(count,0, sigvalues.second); // to correct
    497       }
    498       else
    499       {
    500         gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
    501         gr.SetPointError(count,0, sigvalues.second);
    502       }
    503       count++;
    504     }
    505 
    506     else if (central == false && histos->at(bin).fwdResolHist->GetEntries() > 10)
    507     {
    508       std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).fwdResolHist);
    509       if (rms == true)
    510       {
    511         gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
    512         gr.SetPointError(count,0, sigvalues.second); // to correct
    513       }
    514       else
    515       {
    516         gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
    517         gr.SetPointError(count,0, sigvalues.second);
    518       }
    519       count++;
    520     }
    521 
    522   }
    523   return gr;
    524 }
    525 
    526 
    527 //------------------------------------------------------------------------------
    528 
    529 
    530 // type 1 : object, 2 : track, 3 : tower
    531 
    532 void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
    533 {
    534 
    535   gr->SetLineWidth(2);
    536 
    537   switch ( type )
    538   {
    539     case 1:
    540       gr->SetLineColor(objColor);
    541       gr->SetLineStyle(objStyle);
    542       std::cout << "Adding " << gr->GetName() << std::endl;
    543       mg->Add(gr);
    544       leg->AddEntry(gr,"Reco","l");
    545       break;
    546 
    547     case 2:
    548       gr->SetLineColor(trackColor);
    549       gr->SetLineStyle(trackStyle);
    550       mg->Add(gr);
    551       leg->AddEntry(gr,"Track","l");
    552       break;
    553 
    554     case 3:
    555       gr->SetLineColor(towerColor);
    556       gr->SetLineStyle(towerStyle);
    557       mg->Add(gr);
    558       leg->AddEntry(gr,"Tower","l");
    559       break;
    560 
    561     case 0:
    562       gr->SetLineColor(objColor);
    563       gr->SetLineStyle(objStyle);
    564       mg->Add(gr);
    565       break;
    566 
    567     default:
    568       std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
    569       break;
    570   }
    571 }
    572 
    573 void addHist(TH1D *h, TLegend *leg, int type)
    574 {
    575   h->SetLineWidth(2);
    576 
    577   switch ( type )
    578   {
    579     case 1:
    580       h->SetLineColor(objColor);
    581       h->SetLineStyle(objStyle);
    582       leg->AddEntry(h,"Reco","l");
    583       break;
    584 
    585     case 2:
    586       h->SetLineColor(trackColor);
    587       h->SetLineStyle(trackStyle);
    588       leg->AddEntry(h,"Track","l");
    589       break;
    590 
    591     case 3:
    592       h->SetLineColor(towerColor);
    593       h->SetLineStyle(towerStyle);
    594       leg->AddEntry(h,"Tower","l");
    595       break;
    596 
    597     case 0:
    598       h->SetLineColor(objColor);
    599       h->SetLineStyle(objStyle);
    600       break;
    601 
    602     default:
    603       std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
    604       break;
    605   }
    606 }
    607 
    608 void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
    609 {
    610   mg->SetMinimum(0.);
    611   mg->SetMaximum(max);
    612   mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
    613   mg->GetYaxis()->SetTitle("resolution");
    614   mg->GetXaxis()->SetTitle("p_{T} [GeV]");
    615   mg->GetYaxis()->SetTitleSize(0.07);
    616   mg->GetXaxis()->SetTitleSize(0.07);
    617   mg->GetYaxis()->SetLabelSize(0.06);
    618   mg->GetXaxis()->SetLabelSize(0.06);
    619   mg->GetYaxis()->SetLabelOffset(0.03);
    620   mg->GetYaxis()->SetTitleOffset(1.4);
    621   mg->GetXaxis()->SetTitleOffset(1.4);
    622 
    623   mg->GetYaxis()->SetNdivisions(505);
    624 
    625   leg->SetBorderSize(0);
    626   leg->SetShadowColor(0);
    627   leg->SetFillColor(0);
    628   leg->SetFillStyle(0);
    629 
    630   gStyle->SetOptTitle(0);
    631   gPad->SetLogx();
    632   gPad->SetBottomMargin(0.2);
    633   gPad->SetLeftMargin(0.2);
    634   gPad->Modified();
    635   gPad->Update();
    636 
    637 }
    638 
    639 void DrawAxis(TH1D *h, TLegend *leg, int type)
    640 {
    641 
    642   h->GetYaxis()->SetRangeUser(0,1.0);
    643   if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
    644   else h->GetXaxis()->SetTitle("#eta");
    645   h->GetYaxis()->SetTitle("efficiency");
    646   h->GetYaxis()->SetTitleSize(0.07);
    647   h->GetXaxis()->SetTitleSize(0.07);
    648   h->GetYaxis()->SetLabelSize(0.06);
    649   h->GetXaxis()->SetLabelSize(0.06);
    650   h->GetYaxis()->SetLabelOffset(0.03);
    651   h->GetYaxis()->SetTitleOffset(1.3);
    652   h->GetXaxis()->SetTitleOffset(1.4);
    653 
    654   h->GetYaxis()->SetNdivisions(505);
    655 
    656   leg->SetBorderSize(0);
    657   leg->SetShadowColor(0);
    658   leg->SetFillColor(0);
    659   leg->SetFillStyle(0);
    660 
    661   gStyle->SetOptTitle(0);
    662   gStyle->SetOptStat(0);
    663   gPad->SetBottomMargin(0.2);
    664   gPad->SetLeftMargin(0.2);
    665 
    666   gPad->Modified();
    667   gPad->Update();
    668 
    669 }
    670 
    671 
    672 void Validation(const char *inputFileElectron, const char *inputFileMuon, const char *inputFilePhoton, const char *inputFileJet, const char *inputFileBJet, const char *inputFileTauJet, const char *outputFile)
    673 {
    674   TChain *chainElectron = new TChain("Delphes");
    675   chainElectron->Add(inputFileElectron);
    676   ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
    677 
    678   TChain *chainMuon = new TChain("Delphes");
    679   chainMuon->Add(inputFileMuon);
    680   ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
    681 
    682   TChain *chainPhoton = new TChain("Delphes");
    683   chainPhoton->Add(inputFilePhoton);
    684   ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
    685 
    686   TChain *chainJet = new TChain("Delphes");
    687   chainJet->Add(inputFileJet);
    688   ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
    689 
    690   TChain *chainBJet = new TChain("Delphes");
    691   chainBJet->Add(inputFileBJet);
    692   ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
    693 
    694   TChain *chainTauJet = new TChain("Delphes");
    695   chainTauJet->Add(inputFileTauJet);
    696   ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
    697 
    698   TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
    699   TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
    700   TClonesArray *branchTowerElectron = treeReaderElectron->UseBranch("Tower");
    701   TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
    702 
    703   TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
    704   TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
    705   TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
    706 
    707   TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
    708   TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
    709   TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
    710 
    711   TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
    712   TClonesArray *branchPFJet = treeReaderJet->UseBranch("Jet");
    713   TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
    714 
    715   TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
    716   TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
    717 
    718   TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
    719   TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
    720 
    721   TClonesArray *branchScalarHT = treeReaderJet->UseBranch("ScalarHT");
    722   TClonesArray *branchMet = treeReaderJet->UseBranch("MissingET");
    723 
    724   ///////////////
    725   // Electrons //
    726   ///////////////
    727 
    728   // Reconstruction efficiency
    729   TString elecs = "Electron";
    730   int elID = 11;
    731   std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticleElectron, "Electron", elID, treeReaderElectron);
    732 
    733   // tracking reconstruction efficiency
    734   std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrackElectron, branchParticleElectron, "electronTrack", elID, treeReaderElectron);
    735 
    736   // Tower reconstruction efficiency
    737   std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTowerElectron, branchParticleElectron, "electronTower", elID, treeReaderElectron);
    738 
    739   // Electron Energy Resolution
    740   std::vector<resolPlot> plots_el;
    741   HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
    742   GetEres<Electron>(&plots_el, branchElectron, branchParticleElectron, elID, treeReaderElectron);
    743   TGraphErrors gr_el = EresGraph(&plots_el, true);
    744   TGraphErrors gr_elFwd = EresGraph(&plots_el, false);
    745   gr_el.SetName("Electron");
    746   gr_elFwd.SetName("ElectronFwd");
    747 
    748   // Electron Track Energy Resolution
    749   std::vector<resolPlot> plots_eltrack;
    750   HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
    751   GetEres<Track>(&plots_eltrack, branchTrackElectron, branchParticleElectron, elID, treeReaderElectron);
    752   TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
    753   TGraphErrors gr_eltrackFwd = EresGraph(&plots_eltrack, false);
    754   gr_eltrack.SetName("ElectronTracks");
    755   gr_eltrackFwd.SetName("ElectronTracksFwd");
    756 
    757   // Electron Tower Energy Resolution
    758   std::vector<resolPlot> plots_eltower;
    759   HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
    760   GetEres<Tower>(&plots_eltower, branchTowerElectron, branchParticleElectron, elID, treeReaderElectron);
    761   TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
    762   TGraphErrors gr_eltowerFwd = EresGraph(&plots_eltower, false);
    763   gr_eltower.SetName("ElectronTower");
    764   gr_eltrackFwd.SetName("ElectronTracksFwd");
    765 
    766   // Canvases
    767   TString elEff = "electronEff";
    768   TCanvas *C_el1 = new TCanvas(elEff,elEff, 1600, 600);
    769   C_el1->Divide(2);
    770   C_el1->cd(1);
    771   TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
    772   leg_el1->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
    773   leg_el1->AddEntry("","","");
    774 
    775   gPad->SetLogx();
    776   histos_eltrack.first->Draw("][");
    777   addHist(histos_eltrack.first, leg_el1, 2);
    778   histos_el.first->Draw("same ][");
    779   addHist(histos_el.first, leg_el1, 1);
    780   DrawAxis(histos_eltrack.first, leg_el1, 0);
    781 
    782   leg_el1->Draw();
    783 
    784   C_el1->cd(2);
    785   TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
    786   leg_el2->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
    787   leg_el2->AddEntry("","","");
    788 
    789   gPad->SetLogx();
    790   histos_eltrack.second->Draw("][");
    791   addHist(histos_eltrack.second, leg_el2, 2);
    792   histos_el.second->Draw("same ][");
    793   addHist(histos_el.second, leg_el2, 1);
    794 
    795   DrawAxis(histos_eltrack.second, leg_el2, 0);
    796   leg_el2->Draw();
    797 
    798   C_el1->cd(0);
    799 
    800   TString elRes = "electronERes";
    801   TString elResFwd = "electronEResForward";
    802   TCanvas *C_el2 = new TCanvas(elRes,elRes, 1600, 600);
    803   C_el2->Divide(2);
    804   C_el2->cd(1);
    805   TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
    806   TLegend *leg_el = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    807   leg_el->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
    808   leg_el->AddEntry("","","");
    809 
    810   addGraph(mg_el, &gr_eltower, leg_el, 3);
    811   addGraph(mg_el, &gr_eltrack, leg_el, 2);
    812   addGraph(mg_el, &gr_el, leg_el, 1);
    813 
    814   mg_el->Draw("ACX");
    815   leg_el->Draw();
    816 
    817   DrawAxis(mg_el, leg_el, 0.1);
    818 
    819   C_el2->cd(2);
    820   TMultiGraph *mg_elFwd = new TMultiGraph(elResFwd,elResFwd);
    821   TLegend *leg_elFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    822   leg_elFwd->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
    823   leg_elFwd->AddEntry("","","");
    824 
    825   addGraph(mg_elFwd, &gr_eltowerFwd, leg_elFwd, 3);
    826   addGraph(mg_elFwd, &gr_eltrackFwd, leg_elFwd, 2);
    827   addGraph(mg_elFwd, &gr_elFwd, leg_elFwd, 1);
    828 
    829   mg_elFwd->Draw("ACX");
    830   leg_elFwd->Draw();
    831 
    832   DrawAxis(mg_elFwd, leg_elFwd, 0.2);
    833 
    834   C_el1->Print("validation.pdf(","pdf");
    835   C_el2->Print("validation.pdf","pdf");
    836 
    837   gDirectory->cd(0);
    838 
    839   ///////////
    840   // Muons //
    841   ///////////
    842 
    843   // Reconstruction efficiency
    844   int muID = 13;
    845   std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticleMuon,"Muon", muID, treeReaderMuon);
    846 
    847   // muon tracking reconstruction efficiency
    848   std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
    849 
    850   // Muon Energy Resolution
    851   std::vector<resolPlot> plots_mu;
    852   HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
    853   GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
    854   TGraphErrors gr_mu = EresGraph(&plots_mu, true);
    855   TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
    856   gr_mu.SetName("Muon");
    857   gr_muFwd.SetName("MuonFwd");
    858 
    859   // Muon Track Energy Resolution
    860   std::vector<resolPlot> plots_mutrack;
    861   HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
    862   GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
    863   TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
    864   TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
    865   gr_mutrackFwd.SetName("MuonTracksFwd");
    866 
    867   // Canvas
    868 
    869   TString muEff = "muonEff";
    870   TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600);
    871   C_mu1->Divide(2);
    872   C_mu1->cd(1);
    873   TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
    874   leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}");
    875   leg_mu1->AddEntry("","","");
    876 
    877 
    878   gPad->SetLogx();
    879   histos_mutrack.first->Draw("][");
    880   addHist(histos_mutrack.first, leg_mu1, 2);
    881   histos_mu.first->Draw("same ][");
    882   addHist(histos_mu.first, leg_mu1, 1);
    883 
    884   DrawAxis(histos_mutrack.first, leg_mu1, 0);
    885 
    886   leg_mu1->Draw();
    887 
    888   C_mu1->cd(2);
    889   TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
    890   leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
    891   leg_mu2->AddEntry("","","");
    892 
    893   gPad->SetLogx();
    894   histos_mutrack.second->Draw("][");
    895   addHist(histos_mutrack.second, leg_mu2, 2);
    896   histos_mu.second->Draw("same ][");
    897   addHist(histos_mu.second, leg_mu2, 1);
    898 
    899   DrawAxis(histos_mutrack.second, leg_mu2, 1);
    900   leg_mu2->Draw();
    901 
    902   TString muRes = "muonERes";
    903   TString muResFwd = "muonEResFwd";
    904 
    905   TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600);
    906   C_mu->Divide(2);
    907   C_mu->cd(1);
    908   TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
    909   TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
    910   leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}");
    911   leg_mu->AddEntry("","","");
    912 
    913   addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
    914   addGraph(mg_mu, &gr_mu, leg_mu, 1);
    915 
    916   mg_mu->Draw("ACX");
    917   leg_mu->Draw();
    918 
    919   DrawAxis(mg_mu, leg_mu, 0.3);
    920 
    921   C_mu->cd(2);
    922   TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd);
    923   TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
    924   leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
    925   leg_muFwd->AddEntry("","","");
    926 
    927   addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2);
    928   addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1);
    929 
    930   mg_muFwd->Draw("ACX");
    931   leg_muFwd->Draw();
    932 
    933   DrawAxis(mg_muFwd, leg_muFwd, 0.3);
    934 
    935   //C_mu1->SaveAs(muEff+".eps");
    936   //C_mu->SaveAs(muRes+".eps");
    937 
    938   C_mu1->Print("validation.pdf","pdf");
    939   C_mu->Print("validation.pdf","pdf");
    940 
    941   gDirectory->cd(0);
    942 
    943   /////////////
    944   // Photons //
    945   /////////////
    946 
    947   // Reconstruction efficiency
    948   int phID = 22;
    949   std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
    950   std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTowerPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
    951 
    952   // Photon Energy Resolution
    953   std::vector<resolPlot> plots_ph;
    954   HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
    955   GetEres<Photon>(&plots_ph, branchPhoton, branchParticlePhoton, phID, treeReaderPhoton);
    956   TGraphErrors gr_ph = EresGraph(&plots_ph, true);
    957   TGraphErrors gr_phFwd = EresGraph(&plots_ph, false);
    958   gr_ph.SetName("Photon");
    959   gr_phFwd.SetName("PhotonFwd");
    960 
    961 
    962   // Photon Tower Energy Resolution
    963   std::vector<resolPlot> plots_phtower;
    964   HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
    965   GetEres<Tower>(&plots_phtower, branchTowerPhoton, branchParticlePhoton, phID, treeReaderPhoton);
    966   TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
    967   TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false);
    968   gr_phtower.SetName("PhotonTower");
    969   gr_phtowerFwd.SetName("PhotonTowerFwd");
    970 
    971   // Canvas
    972 
    973   TString phEff = "photonEff";
    974   TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600);
    975   C_ph1->Divide(2);
    976   C_ph1->cd(1);
    977   TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
    978   leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}");
    979   leg_ph1->AddEntry("","","");
    980 
    981 
    982   gPad->SetLogx();
    983   histos_phtower.first->Draw("][");
    984   addHist(histos_phtower.first, leg_ph1, 3);
    985   histos_ph.first->Draw("same ][");
    986   addHist(histos_ph.first, leg_ph1, 1);
    987 
    988   DrawAxis(histos_phtower.first, leg_ph1, 0);
    989   leg_ph1->Draw();
    990 
    991   C_ph1->cd(2);
    992   TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
    993   leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
    994   leg_ph2->AddEntry("","","");
    995 
    996 
    997   gPad->SetLogx();
    998   histos_phtower.second->Draw("][");
    999   addHist(histos_phtower.second, leg_ph2, 3);
    1000   histos_ph.second->Draw("same ][");
    1001   addHist(histos_ph.second, leg_ph2, 1);
    1002 
    1003   DrawAxis(histos_phtower.second, leg_ph2, 1);
    1004   leg_ph2->Draw();
    1005 
    1006   C_ph1->SaveAs(phEff+".eps");
    1007 
    1008   TString phRes = "phERes";
    1009   TString phResFwd = "phEResFwd";
    1010 
    1011   TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600);
    1012   C_ph->Divide(2);
    1013   C_ph->cd(1);
    1014   TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
    1015   TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1016   leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}");
    1017   leg_ph->AddEntry("","","");
    1018 
    1019   addGraph(mg_ph, &gr_phtower, leg_ph, 3);
    1020   addGraph(mg_ph, &gr_ph, leg_ph, 1);
    1021 
    1022   mg_ph->Draw("ACX");
    1023   leg_ph->Draw();
    1024 
    1025   DrawAxis(mg_ph, leg_ph, 0.3);
    1026 
    1027   C_ph->cd(2);
    1028   TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd);
    1029   TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1030   leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
    1031   leg_phFwd->AddEntry("","","");
    1032 
    1033   addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3);
    1034   addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1);
    1035 
    1036   mg_phFwd->Draw("ACX");
    1037   leg_phFwd->Draw();
    1038 
    1039   DrawAxis(mg_phFwd, leg_phFwd, 0.3);
    1040 
    1041   C_ph->SaveAs(phRes+".eps");
    1042 
    1043   C_ph1->Print("validation.pdf","pdf");
    1044   C_ph->Print("validation.pdf","pdf");
    1045 
    1046   gDirectory->cd(0);
    1047 
    1048   //////////
    1049   // Jets //
    1050   //////////
    1051 
    1052   // BJets Reconstruction efficiency
    1053   int bID = 5;
    1054   std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFBJet, branchParticleBJet,"BTag", bID, treeReaderBJet);
    1055 
    1056   // TauJets Reconstruction efficiency
    1057   int tauID = 15;
    1058   std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFTauJet, branchParticleTauJet,"TauTag", tauID, treeReaderTauJet);
    1059 
    1060   // PFJets Energy Resolution
    1061   std::vector<resolPlot> plots_pfjets;
    1062   HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
    1063   GetJetsEres(&plots_pfjets, branchPFJet, branchGenJet, treeReaderJet);
    1064   TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
    1065   TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false);
    1066   gr_pfjets.SetName("pfJet");
    1067   gr_pfjetsFwd.SetName("pfJetFwd");
    1068 
    1069   // CaloJets Energy Resolution
    1070   std::vector<resolPlot> plots_calojets;
    1071   HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
    1072   GetJetsEres(&plots_calojets, branchCaloJet, branchGenJet, treeReaderJet);
    1073   TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
    1074   TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false);
    1075   gr_calojets.SetName("caloJet");
    1076   gr_calojetsFwd.SetName("caloJetFwd");
    1077 
    1078   // MET Resolution vs HT
    1079   std::vector<resolPlot> plots_met;
    1080   HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500);
    1081   GetMetres(&plots_met, branchScalarHT, branchMet, branchPFJet, treeReaderJet);
    1082   TGraphErrors gr_met = EresGraph(&plots_met, true);
    1083   gr_calojets.SetName("MET");
    1084 
    1085   // Canvas
    1086   TString btagEff = "btagEff";
    1087   TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600);
    1088   C_btag1->Divide(2);
    1089   C_btag1->cd(1);
    1090   TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1091   leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}");
    1092   leg_btag1->AddEntry("","","");
    1093 
    1094   gPad->SetLogx();
    1095   histos_btag.first->Draw();
    1096   addHist(histos_btag.first, leg_btag1, 0);
    1097 
    1098   DrawAxis(histos_btag.first, leg_btag1, 0);
    1099   leg_btag1->Draw();
    1100 
    1101   C_btag1->cd(2);
    1102   TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
    1103   leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}");
    1104   leg_btag2->AddEntry("","","");
    1105 
    1106   gPad->SetLogx();
    1107   histos_btag.second->Draw();
    1108   addHist(histos_btag.second, leg_btag2, 0);
    1109 
    1110   DrawAxis(histos_btag.second, leg_btag2, 0);
    1111   leg_btag2->Draw();
    1112   C_btag1->cd(0);
    1113 
    1114   TString tautagEff = "tautagEff";
    1115   TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600);
    1116   C_tautag1->Divide(2);
    1117   C_tautag1->cd(1);
    1118   TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1119   leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}");
    1120   leg_tautag1->AddEntry("","","");
    1121 
    1122   gPad->SetLogx();
    1123   histos_tautag.first->Draw();
    1124   addHist(histos_tautag.first, leg_tautag1, 0);
    1125 
    1126   DrawAxis(histos_tautag.first, leg_tautag1, 0);
    1127   leg_tautag1->Draw();
    1128 
    1129   C_tautag1->cd(2);
    1130   TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
    1131   leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}");
    1132   leg_tautag2->AddEntry("","","");
    1133 
    1134   gPad->SetLogx();
    1135   histos_tautag.second->Draw();
    1136   addHist(histos_tautag.second, leg_tautag2, 0);
    1137 
    1138   DrawAxis(histos_tautag.second, leg_tautag2, 0);
    1139   leg_tautag2->Draw();
    1140   C_tautag1->cd(0);
    1141 
    1142   TString jetRes = "jetERes";
    1143   TString jetResFwd = "jetEResFwd";
    1144   TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600);
    1145   C_jet->Divide(2);
    1146 
    1147   C_jet->cd(1);
    1148   TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
    1149   TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1150   leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}");
    1151   leg_jet->AddEntry("","","");
    1152 
    1153   addGraph(mg_jet, &gr_calojets, leg_jet, 3);
    1154   addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
    1155 
    1156   mg_jet->Draw("ACX");
    1157   leg_jet->Draw();
    1158 
    1159   DrawAxis(mg_jet, leg_jet, 0.25);
    1160 
    1161   C_jet->cd(2);
    1162   TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
    1163   TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1164   leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}");
    1165   leg_jetFwd->AddEntry("","","");
    1166 
    1167   addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3);
    1168   addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
    1169 
    1170   mg_jetFwd->Draw("ACX");
    1171   leg_jetFwd->Draw();
    1172 
    1173   DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
    1174 
    1175   C_btag1->SaveAs(btagEff+".eps");
    1176   C_jet->SaveAs(jetRes+".eps");
    1177 
    1178   TString metRes = "MetRes";
    1179   TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600);
    1180 
    1181   TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
    1182   TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax);
    1183   leg_met->SetHeader("E_{T}^{miss}");
    1184   leg_met->AddEntry("","","");
    1185 
    1186 
    1187   addGraph(mg_met, &gr_met, leg_met, 0);
    1188 
    1189   mg_met->Draw("ACX");
    1190   leg_met->Draw();
    1191 
    1192   DrawAxis(mg_met, leg_met, 300);
    1193 
    1194   mg_met->GetXaxis()->SetTitle("H_{T} [GeV]");
    1195   mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]");
    1196 
    1197   C_met->SaveAs(metRes+".eps");
    1198 
    1199   C_jet->Print("validation.pdf","pdf");
    1200   C_btag1->Print("validation.pdf","pdf");
    1201   C_tautag1->Print("validation.pdf","pdf");
    1202   C_met->Print("validation.pdf)","pdf");
    1203 
    1204   TFile *fout = new TFile(outputFile,"recreate");
    1205 
    1206   for (int bin = 0; bin < Nbins; bin++)
    1207   {
    1208     plots_el.at(bin).cenResolHist->Write();
    1209     plots_eltrack.at(bin).cenResolHist->Write();
    1210     plots_eltower.at(bin).cenResolHist->Write();
    1211     plots_el.at(bin).fwdResolHist->Write();
    1212     plots_eltrack.at(bin).fwdResolHist->Write();
    1213     plots_eltower.at(bin).fwdResolHist->Write();
    1214   }
    1215 
    1216   histos_el.first->Write();
    1217   histos_el.second->Write();
    1218   histos_eltrack.first->Write();
    1219   histos_eltrack.second->Write();
    1220   histos_eltower.first->Write();
    1221   histos_eltower.second->Write();
    1222   C_el1->Write();
    1223   C_el2->Write();
    1224 
    1225   histos_mu.first->Write();
    1226   histos_mu.second->Write();
    1227   histos_mutrack.first->Write();
    1228   histos_mutrack.second->Write();
    1229   C_mu1->Write();
    1230   C_mu->Write();
    1231 
    1232   histos_ph.first->Write();
    1233   histos_ph.second->Write();
    1234   C_ph1->Write();
    1235   C_ph->Write();
    1236 
    1237   for (int bin = 0; bin < Nbins; bin++)
    1238   {
    1239     plots_pfjets.at(bin).cenResolHist->Write();
    1240     plots_pfjets.at(bin).fwdResolHist->Write();
    1241     plots_calojets.at(bin).cenResolHist->Write();
    1242     plots_calojets.at(bin).fwdResolHist->Write();
    1243     plots_met.at(bin).cenResolHist->Write();
    1244   }
    1245   histos_btag.first->Write();
    1246   histos_btag.second->Write();
    1247   histos_tautag.first->Write();
    1248   histos_tautag.second->Write();
    1249   C_btag1->Write();
    1250   C_tautag1->Write();
    1251   C_jet->Write();
    1252   C_met->Write();
    1253 
    1254   fout->Write();
    1255 
    1256   cout << "** Exiting..." << endl;
    1257 
    1258   delete treeReaderElectron;
    1259   delete treeReaderMuon;
    1260   delete treeReaderPhoton;
    1261   delete treeReaderJet;
    1262   delete treeReaderBJet;
    1263   delete treeReaderTauJet;
    1264   delete chainElectron;
    1265   delete chainMuon;
    1266   delete chainPhoton;
    1267   delete chainJet;
    1268   delete chainBJet;
    1269   delete chainTauJet;
    1270 }
     51#include "Validation.C"
    127152
    127253//------------------------------------------------------------------------------
    … …  
    127859  if(argc != 3)
    127960  {
    1280     cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl;
    1281     cout << " input_file_electron  - input file in ROOT format ('Delphes' tree)," << endl;
    1282     cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
    1283     cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
    1284     cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
    1285     cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
    1286     cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
     61    cout << " Usage: " << appName << " input_file output_file" << endl;
     62    cout << " input_file  - input file in ROOT format ('Delphes' tree)," << endl;
    128763    cout << " output_file - output file in ROOT format" << endl;
    128864    return 1;
    … …  
    129571  TApplication app(appName, &appargc, appargv);
    129672
    1297   Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
     73  TString inputFile(argv[1]);
     74  TString outputFile(argv[2]);
     75
     76
     77//------------------------------------------------------------------------------
     78// Here you call your macro's main function
     79
     80  Validation(inputFile, outputFile);
     81
     82//------------------------------------------------------------------------------
     83
    129884}
    129985

  • readers/DelphesROOT.cpp

    r2264876 r7993cad  
    167167      modularDelphes->Clear();
    168168      treeWriter->Clear();
    169       for(Int_t entry = 0; entry < numberOfEvents && !interrupted; ++entry)
     169      for(Int_t entry = 0; entry < numberOfEvents; ++entry)
    170170      {
    171    
     171       
    172172        treeReader->ReadEntry(entry);
    173173       
Note: See TracChangeset for help on using the changeset viewer.