Changeset d77b51d in git for examples/Example4.C
- Timestamp:
- Sep 29, 2015, 2:08:10 PM (9 years ago)
- 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. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/Example4.C
rd870fc5 rd77b51d 1 1 /* 2 3 2 This macro shows how to compute jet energy scale. 4 3 root -l examples/Example4.C'("delphes_output.root", "plots.root")' 5 4 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 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 12 12 13 13 e.g a smooth function: 14 14 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) } 18 16 19 17 or a binned function: 20 21 18 22 19 set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0) * (1.10) + … … 27 24 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0) * (1.00)} 28 25 29 30 Be aware that a binned jet energy scale can produce "steps" in the corrected jet pt distribution ... 31 32 33 26 Be aware that a binned jet energy scale can produce "steps" in the corrected 27 jet pt distribution ... 34 28 */ 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 35 39 36 40 //------------------------------------------------------------------------------ … … 161 165 162 166 Jet *jet, *genjet; 163 GenParticle *part ;167 GenParticle *particle; 164 168 TObject *object; 165 169 166 TLorentzVector JetMom, GenJetMom, BestGenJetMom;167 168 Float_t Dr;170 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum; 171 172 Float_t deltaR; 169 173 Float_t pt, eta; 170 174 Long64_t entry; … … 177 181 // Load selected branches with data from specified event 178 182 treeReader->ReadEntry(entry); 179 // cout<<"-- New event -- "<<endl;180 183 181 184 if(entry%500 == 0) cout << "Event number: "<< entry <<endl; … … 186 189 187 190 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; 193 196 194 197 // Loop over all hard partons in event 195 198 for(j = 0; j < branchParticle->GetEntriesFast(); ++j) 196 199 { 197 198 part = (GenParticle*) branchParticle->At(j); 199 200 GenJetMom = part -> P4(); 201 202 //this is simply to avoid warnings from initial state particlehaving infite rapidity ...203 if( GenJetMom.Px() == 0 && GenJetMom.Py() == 0) continue;204 205 // take the closest parton candidate206 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) 207 210 { 208 Dr = GenJetMom.DeltaR(JetMom);209 BestGenJetMom = GenJetMom;211 deltaR = genJetMomentum.DeltaR(jetMomentum); 212 bestGenJetMomentum = genJetMomentum; 210 213 } 211 212 214 } 213 215 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 } 239 238 } 240 239 } 241 240 } 242 243 241 244 242 //------------------------------------------------------------------------------
Note:
See TracChangeset
for help on using the changeset viewer.