Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • examples/Example4.C

    rc9803f7 rfbb617b  
    11/*
     2
    23This macro shows how to compute jet energy scale.
    34root -l examples/Example4.C'("delphes_output.root", "plots.root")'
    45
    5 The output ROOT file contains the pT(MC)/pT(Reco) distributions for
    6 various pT(Reco) and |eta| bins. The peak value of such distribution is
    7 interpreted as the jet energy correction to be applied for that
    8 given pT(Reco), |eta| bin.
    9 
    10 This can be done by modifying the "ScaleFormula" input parameter to
    11 the JetEnergyScale module in the delphes_card_XXX.tcl
     6The output ROOT file contains the pT(MC)/pT(Reco) distributions for various pT(Reco) and |eta| bins.
     7The peak value of such distribution is interpreted as the jet energy correction to be applied for that given pT(Reco), |eta| bin.
     8
     9This can be done by modifying the "ScaleFormula" input parameter to the JetEnergyScale module in the delphes_card_XXX.tcl
     10
     11
    1212
    1313e.g  a smooth function:
    1414
    15   set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0) }
     15
     16  set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0 ) }
     17
    1618
    1719or a binned function:
     20
    1821
    1922  set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0)  * (1.10) +
     
    2427                    (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0)               * (1.00)}
    2528
    26 Be aware that a binned jet energy scale can produce "steps" in the corrected
    27 jet pt distribution ...
     29
     30Be aware that a binned jet energy scale can produce "steps" in the corrected jet pt distribution ...
     31
     32
     33
    2834*/
    29 
    30 #ifdef __CLING__
    31 R__LOAD_LIBRARY(libDelphes)
    32 #include "classes/DelphesClasses.h"
    33 #include "external/ExRootAnalysis/ExRootTreeReader.h"
    34 #include "external/ExRootAnalysis/ExRootResult.h"
    35 #else
    36 class ExRootTreeReader;
    37 class ExRootResult;
    38 #endif
    3935
    4036//------------------------------------------------------------------------------
     
    165161
    166162  Jet *jet, *genjet;
    167   GenParticle *particle;
     163  GenParticle *part;
    168164  TObject *object;
    169165
    170   TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
    171 
    172   Float_t deltaR;
     166  TLorentzVector JetMom, GenJetMom, BestGenJetMom;
     167
     168  Float_t Dr;
    173169  Float_t pt, eta;
    174170  Long64_t entry;
     
    181177    // Load selected branches with data from specified event
    182178    treeReader->ReadEntry(entry);
     179    //  cout<<"--  New event -- "<<endl;
    183180
    184181    if(entry%500 == 0) cout << "Event number: "<< entry <<endl;
     
    189186
    190187      jet = (Jet*) branchJet->At(i);
    191       jetMomentum = jet->P4();
    192 
    193       plots->fJetPT->Fill(jetMomentum.Pt());
    194 
    195       deltaR = 999;
     188      JetMom = jet-> P4();
     189
     190      plots->fJetPT->Fill(JetMom.Pt());
     191
     192      Dr = 999;
    196193
    197194     // Loop over all hard partons in event
    198195     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
    199196     {
    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)
     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 )
    210207        {
    211            deltaR = genJetMomentum.DeltaR(jetMomentum);
    212            bestGenJetMomentum = genJetMomentum;
     208           Dr = GenJetMom.DeltaR(JetMom);
     209           BestGenJetMom = GenJetMom;
    213210        }
     211
    214212      }
    215213
    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       }
     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
    238239    }
    239240  }
    240241}
     242
    241243
    242244//------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.