Fork me on GitHub

Changeset d77b51d in git for examples/Example4.C


Ignore:
Timestamp:
Sep 29, 2015, 2:08:10 PM (9 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
a98c7ef
Parents:
d870fc5 (diff), 06ec139 (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 remote-tracking branch 'upstream/master'

File:
1 edited

Legend:

Unmodified
Added
Removed
  • examples/Example4.C

    rd870fc5 rd77b51d  
    11/*
    2 
    32This macro shows how to compute jet energy scale.
    43root -l examples/Example4.C'("delphes_output.root", "plots.root")'
    54
    6 The output ROOT file contains the pT(MC)/pT(Reco) distributions for various pT(Reco) and |eta| bins.
    7 The peak value of such distribution is interpreted as the jet energy correction to be applied for that given pT(Reco), |eta| bin.
    8 
    9 This can be done by modifying the "ScaleFormula" input parameter to the JetEnergyScale module in the delphes_card_XXX.tcl
    10 
    11 
     5The output ROOT file contains the pT(MC)/pT(Reco) distributions for
     6various pT(Reco) and |eta| bins. The peak value of such distribution is
     7interpreted as the jet energy correction to be applied for that
     8given pT(Reco), |eta| bin.
     9
     10This can be done by modifying the "ScaleFormula" input parameter to
     11the JetEnergyScale module in the delphes_card_XXX.tcl
    1212
    1313e.g  a smooth function:
    1414
    15 
    16   set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0 ) }
    17 
     15  set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0) }
    1816
    1917or a binned function:
    20 
    2118
    2219  set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0)  * (1.10) +
     
    2724                    (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0)               * (1.00)}
    2825
    29 
    30 Be aware that a binned jet energy scale can produce "steps" in the corrected jet pt distribution ...
    31 
    32 
    33 
     26Be aware that a binned jet energy scale can produce "steps" in the corrected
     27jet pt distribution ...
    3428*/
     29
     30#ifdef __CLING__
     31R__LOAD_LIBRARY(libDelphes)
     32#include "classes/DelphesClasses.h"
     33#include "external/ExRootAnalysis/ExRootTreeReader.h"
     34#include "external/ExRootAnalysis/ExRootResult.h"
     35#else
     36class ExRootTreeReader;
     37class ExRootResult;
     38#endif
    3539
    3640//------------------------------------------------------------------------------
     
    161165
    162166  Jet *jet, *genjet;
    163   GenParticle *part;
     167  GenParticle *particle;
    164168  TObject *object;
    165169
    166   TLorentzVector JetMom, GenJetMom, BestGenJetMom;
    167 
    168   Float_t Dr;
     170  TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
     171
     172  Float_t deltaR;
    169173  Float_t pt, eta;
    170174  Long64_t entry;
     
    177181    // Load selected branches with data from specified event
    178182    treeReader->ReadEntry(entry);
    179     //  cout<<"--  New event -- "<<endl;
    180183
    181184    if(entry%500 == 0) cout << "Event number: "<< entry <<endl;
     
    186189
    187190      jet = (Jet*) branchJet->At(i);
    188       JetMom = jet-> P4();
    189 
    190       plots->fJetPT->Fill(JetMom.Pt());
    191 
    192       Dr = 999;
     191      jetMomentum = jet->P4();
     192
     193      plots->fJetPT->Fill(jetMomentum.Pt());
     194
     195      deltaR = 999;
    193196
    194197     // Loop over all hard partons in event
    195198     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
    196199     {
    197 
    198         part = (GenParticle*) branchParticle->At(j);
    199 
    200         GenJetMom = part -> P4();
    201 
    202         //this is simply to avoid warnings from initial state particle having infite rapidity ...
    203         if(GenJetMom.Px() == 0 && GenJetMom.Py() == 0) continue;
    204 
    205         //take the closest parton candidate
    206         if( GenJetMom.DeltaR(JetMom) < Dr )
     200        particle = (GenParticle*) branchParticle->At(j);
     201
     202        genJetMomentum = particle->P4();
     203
     204        // this is simply to avoid warnings from initial state particle
     205        // having infite rapidity ...
     206        if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
     207
     208        // take the closest parton candidate
     209        if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
    207210        {
    208            Dr = GenJetMom.DeltaR(JetMom);
    209            BestGenJetMom = GenJetMom;
     211           deltaR = genJetMomentum.DeltaR(jetMomentum);
     212           bestGenJetMomentum = genJetMomentum;
    210213        }
    211 
    212214      }
    213215
    214      if(Dr < 0.3)
    215      {
    216        pt  = JetMom.Pt();
    217        eta = TMath::Abs(JetMom.Eta());
    218 
    219 
    220        if( pt > 20.0 && pt < 50.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_20_50_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    221        if( pt > 20.0 && pt < 50.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_20_50_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    222 
    223        if( pt > 50.0 && pt < 100.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_50_100_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    224        if( pt > 50.0 && pt < 100.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_50_100_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    225 
    226        if( pt > 100.0 && pt < 200.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_100_200_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    227        if( pt > 100.0 && pt < 200.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_100_200_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    228 
    229        if( pt > 200.0 && pt < 500.0 && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_200_500_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    230        if( pt > 200.0 && pt < 500.0 && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_200_500_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    231 
    232        if( pt > 500.0               && eta > 0.0 && eta < 2.5 ) plots -> fJetRes_Pt_500_inf_Eta_0_25->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    233        if( pt > 500.0               && eta > 2.5 && eta < 5.0 ) plots -> fJetRes_Pt_500_inf_Eta_25_5->Fill(BestGenJetMom.Pt()/JetMom.Pt());
    234 
    235 
    236      }
    237 
    238 
     216      if(deltaR < 0.3)
     217      {
     218        pt  = jetMomentum.Pt();
     219        eta = TMath::Abs(jetMomentum.Eta());
     220
     221
     222        if(pt > 20.0 && pt < 50.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_20_50_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     223        if(pt > 20.0 && pt < 50.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_20_50_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     224
     225        if(pt > 50.0 && pt < 100.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_50_100_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     226        if(pt > 50.0 && pt < 100.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_50_100_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     227
     228        if(pt > 100.0 && pt < 200.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_100_200_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     229        if(pt > 100.0 && pt < 200.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_100_200_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     230
     231        if(pt > 200.0 && pt < 500.0 && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_200_500_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     232        if(pt > 200.0 && pt < 500.0 && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_200_500_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     233
     234        if(pt > 500.0               && eta > 0.0 && eta < 2.5) plots -> fJetRes_Pt_500_inf_Eta_0_25->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     235        if(pt > 500.0               && eta > 2.5 && eta < 5.0) plots -> fJetRes_Pt_500_inf_Eta_25_5->Fill(bestGenJetMomentum.Pt()/jetMomentum.Pt());
     236
     237      }
    239238    }
    240239  }
    241240}
    242 
    243241
    244242//------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.