Changes in examples/Example4.C [c9803f7:fbb617b] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/Example4.C
rc9803f7 rfbb617b 1 1 /* 2 2 3 This macro shows how to compute jet energy scale. 3 4 root -l examples/Example4.C'("delphes_output.root", "plots.root")' 4 5 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 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 12 12 13 13 e.g a smooth function: 14 14 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 16 18 17 19 or a binned function: 20 18 21 19 22 set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0) * (1.10) + … … 24 27 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0) * (1.00)} 25 28 26 Be aware that a binned jet energy scale can produce "steps" in the corrected 27 jet pt distribution ... 29 30 Be aware that a binned jet energy scale can produce "steps" in the corrected jet pt distribution ... 31 32 33 28 34 */ 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 #else36 class ExRootTreeReader;37 class ExRootResult;38 #endif39 35 40 36 //------------------------------------------------------------------------------ … … 165 161 166 162 Jet *jet, *genjet; 167 GenParticle *part icle;163 GenParticle *part; 168 164 TObject *object; 169 165 170 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;171 172 Float_t deltaR;166 TLorentzVector JetMom, GenJetMom, BestGenJetMom; 167 168 Float_t Dr; 173 169 Float_t pt, eta; 174 170 Long64_t entry; … … 181 177 // Load selected branches with data from specified event 182 178 treeReader->ReadEntry(entry); 179 // cout<<"-- New event -- "<<endl; 183 180 184 181 if(entry%500 == 0) cout << "Event number: "<< entry <<endl; … … 189 186 190 187 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; 196 193 197 194 // Loop over all hard partons in event 198 195 for(j = 0; j < branchParticle->GetEntriesFast(); ++j) 199 196 { 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 // 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 ) 210 207 { 211 deltaR = genJetMomentum.DeltaR(jetMomentum);212 bestGenJetMomentum = genJetMomentum;208 Dr = GenJetMom.DeltaR(JetMom); 209 BestGenJetMom = GenJetMom; 213 210 } 211 214 212 } 215 213 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 238 239 } 239 240 } 240 241 } 242 241 243 242 244 //------------------------------------------------------------------------------
Note:
See TracChangeset
for help on using the changeset viewer.