Fork me on GitHub

Changeset 2264876 in git


Ignore:
Timestamp:
Aug 25, 2016, 4:27:14 PM (8 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
bc58cf5
Parents:
7993cad (diff), 1408174 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'dev_01' of github.com:delphes/delphes into dev_01

Files:
1 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r7993cad r2264876  
    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
     145Validation$(ExeSuf): \
     146        tmp/examples/Validation.$(ObjSuf)
     147
     148tmp/examples/Validation.$(ObjSuf): \
     149        examples/Validation.cpp \
    139150        classes/DelphesClasses.h \
    140151        external/ExRootAnalysis/ExRootTreeReader.h \
    … …  
    150161        root2pileup$(ExeSuf) \
    151162        stdhep2pileup$(ExeSuf) \
    152         Example1$(ExeSuf)
     163        Example1$(ExeSuf) \
     164        Validation$(ExeSuf)
    153165
    154166EXECUTABLE_OBJ +=  \
    … …  
    159171        tmp/converters/root2pileup.$(ObjSuf) \
    160172        tmp/converters/stdhep2pileup.$(ObjSuf) \
    161         tmp/examples/Example1.$(ObjSuf)
     173        tmp/examples/Example1.$(ObjSuf) \
     174        tmp/examples/Validation.$(ObjSuf)
    162175
    163176DelphesHepMC$(ExeSuf): \
    … …  
    183196        classes/DelphesLHEFReader.h \
    184197        external/ExRootAnalysis/ExRootTreeWriter.h \
     198        external/ExRootAnalysis/ExRootTreeBranch.h \
     199        external/ExRootAnalysis/ExRootProgressBar.h
     200DelphesROOT$(ExeSuf): \
     201        tmp/readers/DelphesROOT.$(ObjSuf)
     202
     203tmp/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 \
    185211        external/ExRootAnalysis/ExRootTreeBranch.h \
    186212        external/ExRootAnalysis/ExRootProgressBar.h
    … …  
    200226        DelphesHepMC$(ExeSuf) \
    201227        DelphesLHEF$(ExeSuf) \
     228        DelphesROOT$(ExeSuf) \
    202229        DelphesSTDHEP$(ExeSuf)
    203230
    … …  
    205232        tmp/readers/DelphesHepMC.$(ObjSuf) \
    206233        tmp/readers/DelphesLHEF.$(ObjSuf) \
     234        tmp/readers/DelphesROOT.$(ObjSuf) \
    207235        tmp/readers/DelphesSTDHEP.$(ObjSuf)
    208236
    … …  
    344372        modules/PdgCodeFilter.h \
    345373        modules/BeamSpotFilter.h \
     374        modules/RecoPuFilter.h \
    346375        modules/Cloner.h \
    347376        modules/Weighter.h \
    … …  
    351380        modules/VertexSorter.h \
    352381        modules/VertexFinder.h \
     382        modules/VertexFinderDA4D.h \
    353383        modules/ExampleModule.h
    354384tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    798828        external/ExRootAnalysis/ExRootFilter.h \
    799829        external/ExRootAnalysis/ExRootClassifier.h
     830tmp/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
    800839tmp/modules/SimpleCalorimeter.$(ObjSuf): \
    801840        modules/SimpleCalorimeter.$(SrcSuf) \
    … …  
    896935        modules/VertexFinder.$(SrcSuf) \
    897936        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
     944tmp/modules/VertexFinderDA4D.$(ObjSuf): \
     945        modules/VertexFinderDA4D.$(SrcSuf) \
     946        modules/VertexFinderDA4D.h \
    898947        classes/DelphesClasses.h \
    899948        classes/DelphesFactory.h \
    … …  
    9981047        tmp/modules/PileUpJetID.$(ObjSuf) \
    9991048        tmp/modules/PileUpMerger.$(ObjSuf) \
     1049        tmp/modules/RecoPuFilter.$(ObjSuf) \
    10001050        tmp/modules/SimpleCalorimeter.$(ObjSuf) \
    10011051        tmp/modules/StatusPidFilter.$(ObjSuf) \
    … …  
    10101060        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
    10111061        tmp/modules/VertexFinder.$(ObjSuf) \
     1062        tmp/modules/VertexFinderDA4D.$(ObjSuf) \
    10121063        tmp/modules/VertexSorter.$(ObjSuf) \
    10131064        tmp/modules/Weighter.$(ObjSuf)
    … …  
    16151666        tmp/external/tcl/tclVar.$(ObjSuf)
    16161667
     1668modules/VertexFinderDA4D.h: \
     1669        classes/DelphesModule.h \
     1670        classes/DelphesClasses.h
     1671        @touch $@
     1672
    16171673modules/TrackSmearing.h: \
    16181674        classes/DelphesModule.h
    … …  
    19592015
    19602016modules/BTagging.h: \
     2017        classes/DelphesModule.h
     2018        @touch $@
     2019
     2020modules/RecoPuFilter.h: \
    19612021        classes/DelphesModule.h
    19622022        @touch $@

  • examples/Validation.cpp

    r7993cad r2264876  
    2121#include <utility>
    2222#include <vector>
     23#include <typeinfo>
    2324
    2425#include "TROOT.h"
    … …  
    2829#include "TString.h"
    2930
     31#include "TH1.h"
    3032#include "TH2.h"
     33#include "TMath.h"
     34#include "TStyle.h"
     35#include "TGraph.h"
     36#include "TCanvas.h"
    3137#include "THStack.h"
    3238#include "TLegend.h"
    … …  
    3440#include "TClonesArray.h"
    3541#include "TLorentzVector.h"
     42#include "TGraphErrors.h"
     43#include "TMultiGraph.h"
    3644
    3745#include "classes/DelphesClasses.h"
    … …  
    4755//------------------------------------------------------------------------------
    4856
    49 // Here you can put your analysis macro
    50 
    51 #include "Validation.C"
     57double ptrangemin = 10;
     58double ptrangemax = 10000;
     59static const int Nbins = 20;
     60
     61int objStyle = 1;
     62int trackStyle = 7;
     63int towerStyle = 3;
     64
     65Color_t objColor = kBlack;
     66Color_t trackColor = kBlack;
     67Color_t towerColor = kBlack;
     68
     69double effLegXmin = 0.22;
     70double effLegXmax = 0.7;
     71double effLegYmin = 0.22;
     72double effLegYmax = 0.5;
     73
     74double resLegXmin = 0.62;
     75double resLegXmax = 0.9;
     76double resLegYmin = 0.52;
     77double resLegYmax = 0.85;
     78
     79double topLeftLegXmin = 0.22;
     80double topLeftLegXmax = 0.7;
     81double topLeftLegYmin = 0.52;
     82double topLeftLegYmax = 0.85;
     83
     84
     85struct 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
     102resolPlot::resolPlot()
     103{
     104}
     105
     106resolPlot::resolPlot(double ptdown, double ptup, TString object)
     107{
     108  this->set(ptdown,ptup,object);
     109}
     110
     111void 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
     121void 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
     140class ExRootResult;
     141class ExRootTreeReader;
     142
     143//------------------------------------------------------------------------------
     144
     145void 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
     166template<typename T>
     167std::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
     271template<typename T>
     272void 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}
     341void 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
     416void 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
     470std::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
     481TGraphErrors 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
     532void 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
     573void 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
     608void 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
     639void 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
     672void 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}
    521271
    531272//------------------------------------------------------------------------------
    … …  
    591278  if(argc != 3)
    601279  {
    61     cout << " Usage: " << appName << " input_file output_file" << endl;
    62     cout << " input_file  - input file in ROOT format ('Delphes' tree)," << endl;
     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;
    631287    cout << " output_file - output file in ROOT format" << endl;
    641288    return 1;
    … …  
    711295  TApplication app(appName, &appargc, appargv);
    721296
    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 
    84 }
    85 
    86 
     1297  Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
     1298}
     1299
     1300

  • readers/DelphesROOT.cpp

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