1 |
– |
/* |
2 |
– |
Description: <one line class summary> |
3 |
– |
|
4 |
– |
Implementation: |
5 |
– |
n-tuple creator for jet studies |
6 |
– |
*/ |
7 |
– |
// |
1 |
|
// Original Author: "Nathaniel Odell" |
2 |
|
// Secondary Author: Steven Won |
3 |
|
// With contributions from: Andrey Pozdnyakov |
4 |
|
// Created: Thurs April 22 2010 |
5 |
|
// $Id$ |
13 |
– |
// |
14 |
– |
// |
15 |
– |
|
6 |
|
|
7 |
|
// system include files |
8 |
|
#include <memory> |
15 |
|
#include "FWCore/Framework/interface/EventSetup.h" |
16 |
|
#include "FWCore/Framework/interface/Event.h" |
17 |
|
#include "FWCore/Framework/interface/MakerMacros.h" |
28 |
– |
|
18 |
|
#include "FWCore/ParameterSet/interface/ParameterSet.h" |
30 |
– |
|
19 |
|
#include "DataFormats/GeometryVector/interface/GlobalPoint.h" |
20 |
|
#include "Geometry/Records/interface/CaloGeometryRecord.h" |
21 |
+ |
#include "DataFormats/Math/interface/deltaPhi.h" |
22 |
|
|
23 |
+ |
// Libraries for objects |
24 |
|
#include "DataFormats/JetReco/interface/CaloJet.h" |
25 |
|
#include "DataFormats/JetReco/interface/CaloJetCollection.h" |
26 |
|
#include "DataFormats/JetReco/interface/PFJet.h" |
30 |
|
#include "DataFormats/JetReco/interface/Jet.h" |
31 |
|
#include "DataFormats/VertexReco/interface/Vertex.h" |
32 |
|
#include "DataFormats/VertexReco/interface/VertexFwd.h" |
33 |
+ |
#include "DataFormats/BTauReco/interface/JetTag.h" |
34 |
+ |
#include "DataFormats/TrackReco/interface/Track.h" |
35 |
+ |
#include "DataFormats/BeamSpot/interface/BeamSpot.h" |
36 |
+ |
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" |
37 |
+ |
|
38 |
+ |
//#include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h" |
39 |
+ |
|
40 |
+ |
//GenParticles |
41 |
+ |
//#include "DataFormats/HepMCCandidate/interface/GenParticle.h" |
42 |
+ |
|
43 |
+ |
// Need to get correctors |
44 |
+ |
#include "JetMETCorrections/Objects/interface/JetCorrector.h" |
45 |
+ |
|
46 |
+ |
// ntuple storage classes |
47 |
+ |
#include "TCJet.h" |
48 |
+ |
#include "TCPrimaryVtx.h" |
49 |
|
|
50 |
|
// Need for HLT trigger info: |
51 |
|
#include "FWCore/Common/interface/TriggerNames.h" |
52 |
|
#include "DataFormats/Common/interface/TriggerResults.h" |
53 |
|
#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h" |
48 |
– |
// |
54 |
|
|
55 |
+ |
//Root stuff |
56 |
|
#include "TROOT.h" |
57 |
+ |
#include "TH1.h" |
58 |
+ |
#include "TH2.h" |
59 |
+ |
#include "TProfile.h" |
60 |
|
#include "TFile.h" |
61 |
|
#include "TTree.h" |
62 |
|
#include "TString.h" |
87 |
|
virtual void endJob() ; |
88 |
|
|
89 |
|
bool triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger); |
90 |
+ |
double sumPtSquared(const Vertex & v); |
91 |
|
|
92 |
|
// ----------member data --------------------------- |
93 |
|
|
94 |
|
edm::InputTag PFJetHandle_; |
85 |
– |
edm::InputTag CaloJetHandle_; |
86 |
– |
edm::InputTag PFJetHandle_corr; |
87 |
– |
edm::InputTag CaloJetHandle_corr; |
95 |
|
edm::InputTag GenJetHandle_; |
96 |
|
edm::InputTag PrimaryVtxHandle_; |
97 |
< |
|
91 |
< |
|
97 |
> |
edm::InputTag triggerResultsTag_; |
98 |
|
// Counting variables // |
99 |
|
|
100 |
< |
int eventNumber, runNumber, lumiSection; |
101 |
< |
|
102 |
< |
TTree *sTree; |
100 |
> |
int eventNumber, runNumber, lumiSection, bunchCross; |
101 |
> |
double ptHat, qScale; |
102 |
> |
|
103 |
> |
TTree* sTree; |
104 |
|
TFile* ntupleFile; |
105 |
|
|
106 |
< |
TClonesArray* jetP4_ak5PF; |
100 |
< |
TClonesArray* jetP4_ak5PF_corr; |
101 |
< |
TClonesArray* jetVtx3_ak5PF; |
102 |
< |
TClonesArray* jetVtx3_ak5PF_corr; |
103 |
< |
TClonesArray* jetP4_ak5Calo; |
104 |
< |
TClonesArray* jetP4_ak5Calo_corr; |
105 |
< |
TClonesArray* jetVtx3_ak5Calo; |
106 |
< |
TClonesArray* jetVtx3_ak5Calo_corr; |
106 |
> |
TClonesArray* jet_ak5PF; |
107 |
|
TClonesArray* jetP4_ak5Gen; |
108 |
< |
TClonesArray* jetVtx3_ak5Gen; |
109 |
< |
TClonesArray* PrimaryVtx3; |
108 |
> |
TClonesArray* primaryVtx; |
109 |
|
|
110 |
|
bool doGenJets_; |
112 |
– |
bool doCaloJets_; |
111 |
|
bool doPFJets_; |
112 |
< |
|
112 |
> |
bool triggerHLT_; |
113 |
> |
bool isRealData; |
114 |
|
string rootfilename; |
115 |
|
|
116 |
|
//Triggers |
117 |
< |
std::string hlTriggerResults_; |
118 |
< |
edm::TriggerNames triggerNames; |
120 |
< |
std::string hltName_; |
121 |
< |
std::string triggerName_; |
117 |
> |
string hlTriggerResults_, hltName_, triggerName_; |
118 |
> |
TriggerNames triggerNames; |
119 |
|
HLTConfigProvider hltConfig_; |
120 |
< |
edm::InputTag triggerResultsTag_; |
121 |
< |
std::vector<std::string> hlNames; |
120 |
> |
vector<string> hlNames; |
121 |
> |
unsigned int triggerStatus; |
122 |
|
|
123 |
+ |
//Histograms |
124 |
+ |
TH1F * h1_nAssociatedTracks; |
125 |
+ |
TH1F * h1_trackZ; |
126 |
+ |
TH1F * h1_trackDxy; |
127 |
+ |
TH1F * h1_allTrackDeltaZ_PV; |
128 |
+ |
TH1F * h1_allTrackDeltaZ_LJ; |
129 |
+ |
TH1F * h1_ptHat; |
130 |
+ |
TH2F * h2_nAssTracksVsJetPt; |
131 |
+ |
TProfile * p1_nVtcs; |
132 |
|
|
127 |
– |
}; |
133 |
|
|
134 |
< |
MPIntuple::MPIntuple(const edm::ParameterSet& iConfig): |
134 |
> |
}; |
135 |
|
|
136 |
< |
PFJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag")), |
132 |
< |
CaloJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("CaloJetTag")), |
133 |
< |
PFJetHandle_corr(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag_corr")), |
134 |
< |
CaloJetHandle_corr(iConfig.getUntrackedParameter<edm::InputTag>("CaloJetTag_corr")), |
135 |
< |
GenJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag")), |
136 |
< |
PrimaryVtxHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")), |
137 |
< |
doGenJets_(iConfig.getUntrackedParameter<bool>("doGenJets")), |
138 |
< |
doCaloJets_(iConfig.getUntrackedParameter<bool>("doCaloJets")), |
139 |
< |
doPFJets_(iConfig.getUntrackedParameter<bool>("doPFJets")), |
140 |
< |
rootfilename(iConfig.getUntrackedParameter<string>("rootfilename")), |
141 |
< |
hlTriggerResults_(iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")), |
142 |
< |
hltName_(iConfig.getUntrackedParameter<string>("hltName","HLT")) |
136 |
> |
MPIntuple::MPIntuple(const edm::ParameterSet& iConfig) |
137 |
|
{ |
138 |
< |
//edm::TriggerNames |
139 |
< |
// triggerNames(iConfig); |
140 |
< |
|
138 |
> |
PFJetHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag"); |
139 |
> |
GenJetHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag"); |
140 |
> |
PrimaryVtxHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag"); |
141 |
> |
doGenJets_ = iConfig.getUntrackedParameter<bool>("doGenJets"); |
142 |
> |
doPFJets_ = iConfig.getUntrackedParameter<bool>("doPFJets"); |
143 |
> |
triggerHLT_ = iConfig.getUntrackedParameter<bool>("triggerHLT"); |
144 |
> |
rootfilename = iConfig.getUntrackedParameter<string>("rootfilename"); |
145 |
> |
hlTriggerResults_ = iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults"); |
146 |
> |
hltName_ = iConfig.getUntrackedParameter<string>("hltName","RECO"); |
147 |
|
} |
148 |
|
|
149 |
– |
|
149 |
|
MPIntuple::~MPIntuple() |
150 |
|
{ |
151 |
|
|
152 |
|
} |
153 |
|
|
155 |
– |
|
154 |
|
// |
155 |
|
// member functions |
156 |
|
// |
159 |
|
void MPIntuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) |
160 |
|
{ |
161 |
|
|
162 |
< |
eventNumber = iEvent.id().event(); |
163 |
< |
runNumber = iEvent.id().run(); |
162 |
> |
eventNumber = iEvent.id().event(); |
163 |
> |
runNumber = iEvent.id().run(); |
164 |
|
lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock(); |
165 |
+ |
bunchCross = (unsigned int)iEvent.bunchCrossing(); |
166 |
+ |
isRealData = iEvent.isRealData(); |
167 |
+ |
|
168 |
+ |
|
169 |
+ |
edm::Handle<reco::BeamSpot> beamSpotHandle; |
170 |
+ |
iEvent.getByLabel("offlineBeamSpot", beamSpotHandle); |
171 |
+ |
reco::BeamSpot vertexBeamSpot = *beamSpotHandle; |
172 |
+ |
|
173 |
+ |
int pfCount = 0; |
174 |
+ |
int genCount = 0; |
175 |
+ |
int vtxCount = 0; |
176 |
+ |
float ljVertexZ = -999; |
177 |
+ |
float primaryVertexZ = -999; |
178 |
+ |
|
179 |
+ |
////////////////////////// |
180 |
+ |
//Get vertex information// |
181 |
+ |
////////////////////////// |
182 |
+ |
|
183 |
+ |
Handle<reco::VertexCollection> primaryVtcs; |
184 |
+ |
iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs); |
185 |
|
|
186 |
< |
int PFcount = 0; |
187 |
< |
int Calocount = 0; |
188 |
< |
int PFcount_corr = 0; |
189 |
< |
int Calocount_corr = 0; |
190 |
< |
int Gencount = 0; |
191 |
< |
int Vtxcount = 0; |
186 |
> |
for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){ |
187 |
> |
|
188 |
> |
reco::Vertex myVtx = reco::Vertex(*vtx_iter); |
189 |
> |
if(!myVtx.isValid() || myVtx.isFake()) continue; |
190 |
> |
TCPrimaryVtx* vtxCon = new ((*primaryVtx)[vtxCount]) TCPrimaryVtx; |
191 |
> |
|
192 |
> |
vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z()); |
193 |
> |
vtxCon->SetNDof(myVtx.ndof()); |
194 |
> |
vtxCon->SetChi2(myVtx.chi2()); |
195 |
> |
vtxCon->SetNTrks(myVtx.tracksSize()); |
196 |
> |
vtxCon->SetSumPt2Trks(sumPtSquared(myVtx)); |
197 |
|
|
198 |
< |
if(doPFJets_){ |
198 |
> |
if(vtxCount == 0) primaryVertexZ = myVtx.z(); |
199 |
> |
|
200 |
> |
++vtxCount; |
201 |
|
|
202 |
+ |
} |
203 |
+ |
|
204 |
+ |
p1_nVtcs->Fill(runNumber, vtxCount); |
205 |
+ |
|
206 |
+ |
/////////////////////// |
207 |
+ |
//get jet information// |
208 |
+ |
/////////////////////// |
209 |
+ |
|
210 |
+ |
if(doPFJets_){ |
211 |
+ |
|
212 |
+ |
edm::Handle<reco::JetTagCollection> bTagHandle1; |
213 |
+ |
iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1); |
214 |
+ |
const reco::JetTagCollection & bTags1 = *(bTagHandle1.product()); |
215 |
+ |
reco::JetTagCollection::const_iterator jet_it_1; |
216 |
+ |
|
217 |
+ |
edm::Handle<reco::JetTagCollection> bTagHandle2; |
218 |
+ |
iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2); |
219 |
+ |
const reco::JetTagCollection & bTags2 = *(bTagHandle2.product()); |
220 |
+ |
reco::JetTagCollection::const_iterator jet_it_2; |
221 |
+ |
|
222 |
+ |
|
223 |
+ |
const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup); |
224 |
+ |
const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup); |
225 |
+ |
|
226 |
|
Handle<reco::PFJetCollection> PFJets; |
227 |
|
iEvent.getByLabel(PFJetHandle_, PFJets); |
228 |
< |
|
180 |
< |
PFcount = 0; |
181 |
< |
|
228 |
> |
|
229 |
|
for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){ |
230 |
|
|
231 |
|
reco::PFJet myJet = reco::PFJet(*jet_iter); |
185 |
– |
|
186 |
– |
if(myJet.et() > 5){ |
187 |
– |
|
188 |
– |
new ((*jetP4_ak5PF)[PFcount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
189 |
– |
new ((*jetVtx3_ak5PF)[PFcount]) TVector3(myJet.vx(), myJet.vy(), myJet.vz()); |
190 |
– |
|
191 |
– |
++PFcount; |
192 |
– |
} |
193 |
– |
|
194 |
– |
} |
195 |
– |
|
196 |
– |
|
197 |
– |
PFcount_corr = 0; |
198 |
– |
|
199 |
– |
Handle<reco::PFJetCollection> PFJets_corr; |
200 |
– |
iEvent.getByLabel(PFJetHandle_corr, PFJets_corr); |
201 |
– |
|
202 |
– |
for(PFJetCollection::const_iterator jet_iter = PFJets_corr->begin(); jet_iter!= PFJets_corr->end(); ++jet_iter){ |
203 |
– |
|
204 |
– |
reco::PFJet myJet = reco::PFJet(*jet_iter); |
205 |
– |
|
206 |
– |
if(myJet.et() > 5){ |
207 |
– |
|
208 |
– |
new ((*jetP4_ak5PF_corr)[PFcount_corr]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
209 |
– |
new ((*jetVtx3_ak5PF_corr)[PFcount_corr]) TVector3(myJet.vx(), myJet.vy(), myJet.vz()); |
210 |
– |
|
211 |
– |
++PFcount_corr; |
212 |
– |
} |
213 |
– |
} |
214 |
– |
} |
215 |
– |
|
216 |
– |
if(doCaloJets_){ |
217 |
– |
|
218 |
– |
Calocount = 0; |
219 |
– |
|
220 |
– |
Handle<reco::CaloJetCollection> CaloJets; |
221 |
– |
iEvent.getByLabel(CaloJetHandle_, CaloJets); |
222 |
– |
|
223 |
– |
for(CaloJetCollection::const_iterator jet_iter = CaloJets->begin(); jet_iter!= CaloJets->end(); ++jet_iter){ |
224 |
– |
|
225 |
– |
reco::CaloJet myJet = reco::CaloJet(*jet_iter); |
226 |
– |
|
227 |
– |
if(myJet.et() > 5){ |
228 |
– |
|
229 |
– |
new ((*jetP4_ak5Calo)[Calocount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
230 |
– |
new ((*jetVtx3_ak5Calo)[Calocount]) TVector3(myJet.vx(), myJet.vy(), myJet.vz()); |
231 |
– |
|
232 |
– |
++Calocount; |
233 |
– |
|
234 |
– |
} |
235 |
– |
} |
236 |
– |
|
237 |
– |
Handle<reco::CaloJetCollection> CaloJets_corr; |
238 |
– |
iEvent.getByLabel(CaloJetHandle_corr, CaloJets_corr); |
239 |
– |
|
240 |
– |
Calocount_corr = 0; |
241 |
– |
|
242 |
– |
for(CaloJetCollection::const_iterator jet_iter = CaloJets_corr->begin(); jet_iter!= CaloJets_corr->end(); ++jet_iter){ |
243 |
– |
|
244 |
– |
reco::CaloJet myJet = reco::CaloJet(*jet_iter); |
245 |
– |
|
232 |
|
if(myJet.et() > 5){ |
233 |
< |
|
234 |
< |
new ((*jetP4_ak5Calo_corr)[Calocount_corr]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
235 |
< |
new ((*jetVtx3_ak5Calo_corr)[Calocount_corr]) TVector3(myJet.vx(), myJet.vy(), myJet.vz()); |
236 |
< |
|
237 |
< |
++Calocount_corr; |
238 |
< |
} |
239 |
< |
} |
233 |
> |
|
234 |
> |
for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) { |
235 |
> |
if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break; |
236 |
> |
} |
237 |
> |
|
238 |
> |
for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) { |
239 |
> |
if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break; |
240 |
> |
} |
241 |
> |
|
242 |
> |
TCJet* jetCon = new ((*jet_ak5PF)[pfCount]) TCJet; |
243 |
> |
|
244 |
> |
jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
245 |
> |
jetCon->SetVtx(-999.0, -999.0, -999.0); |
246 |
> |
|
247 |
> |
jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction()); |
248 |
> |
jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction()); |
249 |
> |
jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction()); |
250 |
> |
jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction()); |
251 |
> |
|
252 |
> |
jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity()); |
253 |
> |
jetCon->SetNumChPart(myJet.chargedMultiplicity()); |
254 |
> |
|
255 |
> |
jetCon->SetNumChPart(myJet.chargedMultiplicity()); |
256 |
> |
|
257 |
> |
if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second); |
258 |
> |
if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second); |
259 |
> |
|
260 |
> |
float scale2 = correctorL2->correction(myJet); |
261 |
> |
myJet.scaleEnergy(scale2); |
262 |
> |
float scale3 = correctorL3->correction(myJet); |
263 |
> |
myJet.scaleEnergy(scale3); |
264 |
> |
|
265 |
> |
//more corrections? |
266 |
> |
|
267 |
> |
jetCon->SetJetCorr(2, scale2); |
268 |
> |
jetCon->SetJetCorr(3, scale3); |
269 |
> |
|
270 |
> |
///////////////////////// |
271 |
> |
//get associated tracks// |
272 |
> |
///////////////////////// |
273 |
> |
|
274 |
> |
const reco::TrackRefVector &tracks = myJet.getTrackRefs(); |
275 |
> |
|
276 |
> |
float sumTrackZ, sumTrackIP; |
277 |
> |
float meanTrackZ, meanTrackIP; |
278 |
> |
int nAssociatedTracks = -1; |
279 |
> |
|
280 |
> |
sumTrackZ = sumTrackIP = 0; |
281 |
> |
meanTrackZ = meanTrackIP = -999; |
282 |
> |
|
283 |
> |
if(myJet.pt() > 20 && fabs(myJet.eta()) < 2.4){ |
284 |
> |
|
285 |
> |
nAssociatedTracks = 0; |
286 |
> |
|
287 |
> |
for(TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack){ |
288 |
> |
|
289 |
> |
const reco::Track& myTrack = **iTrack; |
290 |
> |
|
291 |
> |
sumTrackZ += myTrack.vz(); |
292 |
> |
sumTrackIP += myTrack.dxy(vertexBeamSpot.position()); |
293 |
> |
++nAssociatedTracks; |
294 |
> |
|
295 |
> |
} |
296 |
> |
|
297 |
> |
meanTrackZ = sumTrackZ/nAssociatedTracks; |
298 |
> |
h1_nAssociatedTracks->Fill(nAssociatedTracks); |
299 |
> |
h2_nAssTracksVsJetPt->Fill(myJet.pt(), nAssociatedTracks); |
300 |
> |
|
301 |
> |
for(TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack){ |
302 |
> |
|
303 |
> |
const reco::Track& myTrack = **iTrack; |
304 |
> |
|
305 |
> |
h1_trackZ->Fill(meanTrackZ - myTrack.vz()); |
306 |
> |
h1_trackDxy->Fill(myTrack.dxy(vertexBeamSpot.position())); |
307 |
> |
|
308 |
> |
} |
309 |
> |
|
310 |
> |
if(pfCount == 0){ |
311 |
> |
|
312 |
> |
ljVertexZ = meanTrackZ; |
313 |
> |
|
314 |
> |
}else{ |
315 |
> |
|
316 |
> |
h1_allTrackDeltaZ_LJ->Fill(ljVertexZ - meanTrackZ); |
317 |
> |
|
318 |
> |
} |
319 |
> |
|
320 |
> |
h1_allTrackDeltaZ_PV->Fill(primaryVertexZ - meanTrackZ); |
321 |
> |
|
322 |
> |
} |
323 |
> |
|
324 |
> |
jetCon->SetVtx(0, 0, meanTrackZ); |
325 |
> |
|
326 |
> |
++pfCount; |
327 |
> |
} |
328 |
> |
} |
329 |
|
} |
330 |
|
|
331 |
< |
if(doGenJets_){ |
331 |
> |
if(!isRealData){ |
332 |
|
|
333 |
+ |
Handle<GenEventInfoProduct> GenInfoHandle; |
334 |
+ |
iEvent.getByLabel("generator", GenInfoHandle); |
335 |
+ |
|
336 |
|
Handle<reco::GenJetCollection> GenJets; |
337 |
|
iEvent.getByLabel(GenJetHandle_, GenJets); |
338 |
< |
|
339 |
< |
|
338 |
> |
|
339 |
> |
ptHat = qScale = -1; |
340 |
> |
|
341 |
> |
if(GenInfoHandle.isValid()){ |
342 |
> |
|
343 |
> |
qScale = GenInfoHandle->qScale(); |
344 |
> |
ptHat = (GenInfoHandle->hasBinningValues() ? GenInfoHandle->binningValues()[0] : 0.0); |
345 |
> |
|
346 |
> |
h1_ptHat->Fill(ptHat); |
347 |
> |
|
348 |
> |
} |
349 |
> |
|
350 |
|
for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){ |
351 |
|
|
352 |
|
reco::GenJet myJet = reco::GenJet(*jet_iter); |
353 |
|
|
354 |
< |
if(myJet.et() > 5){ |
354 |
> |
if(myJet.pt() > 5){ |
355 |
|
|
356 |
< |
new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
269 |
< |
new ((*jetVtx3_ak5Gen)[Gencount]) TVector3(myJet.vx(), myJet.vy(), myJet.vz()); |
356 |
> |
new ((*jetP4_ak5Gen)[genCount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy()); |
357 |
|
|
358 |
< |
++Gencount; |
358 |
> |
++genCount; |
359 |
|
|
360 |
|
} |
361 |
< |
} |
361 |
> |
} |
362 |
|
} |
363 |
|
|
364 |
< |
|
365 |
< |
Handle<reco::VertexCollection> primaryVtx; |
366 |
< |
iEvent.getByLabel(PrimaryVtxHandle_, primaryVtx); |
364 |
> |
/////////////////////////// |
365 |
> |
//get trigger information// |
366 |
> |
/////////////////////////// |
367 |
|
|
368 |
< |
for(VertexCollection::const_iterator vtx_iter = primaryVtx->begin(); vtx_iter!= primaryVtx->end(); ++vtx_iter){ |
368 |
> |
if(triggerHLT_){ |
369 |
> |
|
370 |
> |
edm::Handle<TriggerResults> hltR; |
371 |
> |
triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_); |
372 |
> |
iEvent.getByLabel(triggerResultsTag_,hltR); |
373 |
|
|
374 |
< |
reco::Vertex myVtx = reco::Vertex(*vtx_iter); |
374 |
> |
const TriggerNames & triggerNames = iEvent.triggerNames(*hltR); |
375 |
> |
hlNames=triggerNames.triggerNames(); |
376 |
|
|
377 |
< |
new ((*PrimaryVtx3)[Vtxcount]) TVector3(myVtx.x(), myVtx.y(), myVtx.z()); |
377 |
> |
string MPI_TriggerNames[] = {"HLT_PixelTracks_Multiplicity70", "HLT_MinBiasBSC_NoBPTX", "HLT_PixelTracks_Multiplicity40","HLT_L1Tech_HCAL_HF", "HLT_IsoTrackHB_8E29", "HLT_IsoTrackHE_8E29", "HLT_L1Tech_RPC_TTU_RBst1_collisions", "HLT_L1_BscMinBiasOR_BptxPlusORMinus", "HLT_L1Tech_BSC_halo_forPhysicsBackground", "HLT_L1Tech_BSC_HighMultiplicity", "HLT_MinBiasPixel_DoubleIsoTrack5", "HLT_MinBiasPixel_DoubleTrack", "HLT_MinBiasPixel_SingleTrack", "HLT_ZeroBiasPixel_SingleTrack", "HLT_MinBiasBSC", "HLT_StoppedHSCP_8E29", "HLT_Jet15U_HcalNoiseFiltered", "HLT_QuadJet15U", "HLT_DiJetAve30U_8E29", "HLT_DiJetAve15U_8E29", "HLT_FwdJet20U", "HLT_Jet50U", "HLT_Jet30U", "HLT_Jet15U", "HLT_BTagMu_Jet10U", "HLT_DoubleJet15U_ForwardBackward", "HLT_BTagIP_Jet50U", "HLT_DoubleLooseIsoTau15", "HLT_SingleLooseIsoTau20", "HLT_HT100U", "HLT_MET100", "HLT_MET45"}; |
378 |
|
|
379 |
< |
++Vtxcount; |
379 |
> |
bool triggerPassed = false; |
380 |
> |
triggerStatus = 0x0; |
381 |
|
|
382 |
+ |
for (uint i=0; i<hlNames.size(); ++i) { |
383 |
+ |
|
384 |
+ |
triggerPassed = triggerDecision(hltR, i); |
385 |
+ |
|
386 |
+ |
if(triggerPassed){ |
387 |
+ |
|
388 |
+ |
for (uint j = 0; j != 32; ++j){ |
389 |
+ |
|
390 |
+ |
if (hlNames[i] == MPI_TriggerNames[j]) |
391 |
+ |
{ |
392 |
+ |
//cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl; |
393 |
+ |
triggerStatus |= 0x01 << j; |
394 |
+ |
|
395 |
+ |
} |
396 |
+ |
} |
397 |
+ |
} |
398 |
+ |
} |
399 |
|
} |
400 |
< |
|
291 |
< |
|
292 |
< |
if(Calocount > 0 || Calocount_corr >0 || PFcount_corr > 0 || PFcount > 0 || Gencount > 0) sTree -> Fill(); |
293 |
< |
|
294 |
< |
jetP4_ak5PF->Clear(); |
295 |
< |
jetP4_ak5PF_corr->Clear(); |
296 |
< |
jetVtx3_ak5PF->Clear(); |
297 |
< |
jetP4_ak5Calo->Clear(); |
298 |
< |
jetP4_ak5Calo_corr->Clear(); |
299 |
< |
jetVtx3_ak5Calo->Clear(); |
300 |
< |
jetP4_ak5Gen->Clear(); |
301 |
< |
jetVtx3_ak5Gen->Clear(); |
302 |
< |
|
400 |
> |
//cout<< "total status: "<<hex<<triggerStatus<<endl; |
401 |
|
|
402 |
|
|
403 |
< |
//---------- Filling HLT trigger bits! ------------ |
404 |
< |
|
307 |
< |
edm::Handle<TriggerResults> hltR; |
308 |
< |
triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_); |
309 |
< |
iEvent.getByLabel(triggerResultsTag_,hltR); |
310 |
< |
const edm::TriggerNames & triggerNames = iEvent.triggerNames(*hltR); |
311 |
< |
hlNames=triggerNames.triggerNames(); |
312 |
< |
bool triggerPassed = false; |
313 |
< |
|
314 |
< |
for (uint iT=0; iT<hlNames.size(); ++iT) |
315 |
< |
{ |
316 |
< |
triggerPassed = triggerDecision(hltR, iT); |
317 |
< |
cout<<"trigger name: "<<hlNames[iT]<<" status: "<<triggerPassed<<endl; |
318 |
< |
|
319 |
< |
// ----------->> Fill them Here! <<---------------------- |
320 |
< |
} |
321 |
< |
//---------- ------------------ --------------- |
322 |
< |
|
403 |
> |
|
404 |
> |
if((pfCount > 3 || genCount > 3) && vtxCount > 0) sTree -> Fill(); |
405 |
|
|
406 |
+ |
jet_ak5PF->Clear(); |
407 |
+ |
primaryVtx->Clear(); |
408 |
+ |
jetP4_ak5Gen->Clear(); |
409 |
|
|
410 |
|
} |
411 |
|
|
415 |
|
{ |
416 |
|
|
417 |
|
ntupleFile = new TFile(rootfilename.c_str(), "RECREATE"); |
418 |
< |
sTree = new TTree("dpsTree", "Tree for Jets"); |
419 |
< |
|
420 |
< |
jetP4_ak5PF = new TClonesArray("TLorentzVector"); |
421 |
< |
jetVtx3_ak5PF = new TClonesArray("TVector3"); |
422 |
< |
jetP4_ak5Calo = new TClonesArray("TLorentzVector"); |
423 |
< |
jetVtx3_ak5Calo = new TClonesArray("TVector3"); |
424 |
< |
jetP4_ak5PF_corr = new TClonesArray("TLorentzVector"); |
425 |
< |
jetVtx3_ak5PF_corr = new TClonesArray("TVector3"); |
426 |
< |
jetP4_ak5Calo_corr = new TClonesArray("TLorentzVector"); |
427 |
< |
jetVtx3_ak5Calo_corr = new TClonesArray("TVector3"); |
418 |
> |
sTree = new TTree("mpiTree", "Tree for Jets"); |
419 |
> |
|
420 |
> |
h1_trackZ = new TH1F("h1_trackZ", "#Delta z for associated tracks", 60, -0.3, 0.3); |
421 |
> |
h1_nAssociatedTracks = new TH1F("h1_nAssociatedTracks", "Number of tracks associated to jet", 20, -0.5, 19.5); |
422 |
> |
h1_trackDxy = new TH1F("h1_trackDxy", "dxy for all associated tracks", 50, -0.4, 0.4); |
423 |
> |
h1_allTrackDeltaZ_PV = new TH1F("h1_allTrackDeltaZ_PV", "#Delta z between jet and primary vertex", 60, -0.3, 0.3); |
424 |
> |
h1_allTrackDeltaZ_LJ = new TH1F("h1_allTrackDeltaZ_LJ", "#Delta z between leading jet and other jet vertices", 60, -0.3, 0.3); |
425 |
> |
h1_ptHat = new TH1F("h1_ptHat", "ptHat", 50, 10, 100); |
426 |
> |
h2_nAssTracksVsJetPt = new TH2F("h2_nAssTracksVsJetPt", "nTracks vs Jet p_{T}", 50, 0, 100, 20, -0.5, 19.5); |
427 |
> |
p1_nVtcs = new TProfile("p1_nVtcs", "<N> per run", 5200.0, 132300.0, 137500.0, 0.0, 6.0); |
428 |
> |
|
429 |
> |
jet_ak5PF = new TClonesArray("TCJet"); |
430 |
|
jetP4_ak5Gen = new TClonesArray("TLorentzVector"); |
431 |
< |
jetVtx3_ak5Gen = new TClonesArray("TVector3"); |
432 |
< |
PrimaryVtx3 = new TClonesArray("TVector3"); |
431 |
> |
primaryVtx = new TClonesArray("TCPrimaryVtx"); |
432 |
> |
|
433 |
|
|
434 |
< |
sTree->Branch("jetP4_ak5PF",&jetP4_ak5PF, 6400, 0); |
348 |
< |
sTree->Branch("jetP4_ak5PF_corr",&jetP4_ak5PF_corr, 6400, 0); |
349 |
< |
sTree->Branch("jetVtx3_ak5PF",&jetVtx3_ak5PF, 6400, 0); |
350 |
< |
sTree->Branch("jetVtx3_ak5PF_corr",&jetVtx3_ak5PF, 6400, 0); |
351 |
< |
sTree->Branch("jetP4_ak5Calo",&jetP4_ak5Calo, 6400, 0); |
352 |
< |
sTree->Branch("jetP4_ak5Calo_corr",&jetP4_ak5Calo_corr, 6400, 0); |
353 |
< |
sTree->Branch("jetVtx3_ak5Calo",&jetVtx3_ak5Calo, 6400, 0); |
354 |
< |
sTree->Branch("jetVtx3_ak5Calo_corr",&jetVtx3_ak5Calo, 6400, 0); |
434 |
> |
sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0); |
435 |
|
sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0); |
436 |
< |
sTree->Branch("jetVtx3_ak5Gen",&jetVtx3_ak5Gen, 6400, 0); |
357 |
< |
sTree->Branch("PrimaryVtx",&PrimaryVtx3, 6400, 0); |
436 |
> |
sTree->Branch("primaryVtx",&primaryVtx, 6400, 0); |
437 |
|
|
438 |
|
sTree->Branch("eventNumber",&eventNumber, "eventNumber/I"); |
439 |
|
sTree->Branch("runNumber",&runNumber, "runNumber/I"); |
440 |
|
sTree->Branch("lumiSection",&lumiSection, "lumiSection/I"); |
441 |
+ |
sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i"); |
442 |
+ |
sTree->Branch("isRealData",&isRealData, "isRealData/i"); |
443 |
+ |
sTree->Branch("bunchCross",&bunchCross, "bunchCross/i"); |
444 |
+ |
sTree->Branch("ptHat",&ptHat, "ptHat/d"); |
445 |
|
|
446 |
|
} |
447 |
|
|
448 |
|
// ------------ method called once each job just after ending the event loop ------------ |
449 |
|
void MPIntuple::endJob() |
450 |
|
{ |
451 |
+ |
|
452 |
+ |
ntupleFile->cd(); |
453 |
+ |
|
454 |
+ |
h1_trackZ->Write(); |
455 |
+ |
h1_nAssociatedTracks->Write(); |
456 |
+ |
h1_trackDxy->Write(); |
457 |
+ |
h1_allTrackDeltaZ_PV->Write(); |
458 |
+ |
h1_allTrackDeltaZ_LJ->Write(); |
459 |
+ |
h1_ptHat->Write(); |
460 |
+ |
h2_nAssTracksVsJetPt->Write(); |
461 |
+ |
p1_nVtcs->Write(); |
462 |
+ |
|
463 |
|
ntupleFile->Write(); |
464 |
|
ntupleFile->Close(); |
465 |
|
|
475 |
|
triggerPassed = true; |
476 |
|
} |
477 |
|
return triggerPassed; |
478 |
< |
} |
478 |
> |
} |
479 |
|
|
480 |
+ |
double MPIntuple::sumPtSquared(const Vertex & v) |
481 |
+ |
{ |
482 |
+ |
double sum = 0.; |
483 |
+ |
double pT; |
484 |
+ |
for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) { |
485 |
+ |
pT = (**it).pt(); |
486 |
+ |
double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0; |
487 |
+ |
|
488 |
+ |
sum += pT*pT; |
489 |
+ |
} |
490 |
+ |
return sum; |
491 |
+ |
} |
492 |
|
|
493 |
|
//define this as a plug-in |
494 |
|
DEFINE_FWK_MODULE(MPIntuple); |