ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/PollackPrograms/src/ntupleProducer.cc
Revision: 1.1
Committed: Wed Feb 20 21:39:44 2013 UTC (12 years, 2 months ago) by bpollack
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:
commiting moriond analysis package

File Contents

# Content
1 #include "ntupleProducer.h"
2
3 ntupleProducer::ntupleProducer(const edm::ParameterSet& iConfig)
4 {
5 jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("JetTag");
6 metTag_ = iConfig.getUntrackedParameter<edm::InputTag>("METTag");
7 muonTag_ = iConfig.getUntrackedParameter<edm::InputTag>("MuonTag");
8 electronTag_ = iConfig.getUntrackedParameter<edm::InputTag>("ElectronTag");
9 photonTag_ = iConfig.getUntrackedParameter<edm::InputTag>("PhotonTag");
10 tauTag_ = iConfig.getUntrackedParameter<edm::InputTag>("TauTag");
11 genJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag");
12 primaryVtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag");
13
14 rhoCorrTag_ = iConfig.getUntrackedParameter<edm::InputTag>("rhoCorrTag");
15 rho25CorrTag_ = iConfig.getUntrackedParameter<edm::InputTag>("rho25CorrTag");
16 rhoMuCorrTag_ = iConfig.getUntrackedParameter<edm::InputTag>("rhoMuCorrTag");
17 hlTriggerResults_ = iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults");
18 hltProcess_ = iConfig.getUntrackedParameter<string>("hltName");
19 triggerPaths_ = iConfig.getUntrackedParameter<vector<string> >("triggers");
20
21 partFlowTag_ = iConfig.getUntrackedParameter<edm::InputTag>("partFlowTag");
22
23 saveJets_ = iConfig.getUntrackedParameter<bool>("saveJets");
24 saveElectrons_ = iConfig.getUntrackedParameter<bool>("saveElectrons");
25 saveMuons_ = iConfig.getUntrackedParameter<bool>("saveMuons");
26 saveTaus_ = iConfig.getUntrackedParameter<bool>("saveTaus");
27 savePhotons_ = iConfig.getUntrackedParameter<bool>("savePhotons");
28 saveMET_ = iConfig.getUntrackedParameter<bool>("saveMET");
29 saveGenJets_ = iConfig.getUntrackedParameter<bool>("saveGenJets");
30 saveGenParticles_ = iConfig.getUntrackedParameter<bool>("saveGenParticles");
31
32 ecalTPFilterTag_ = iConfig.getUntrackedParameter<edm::InputTag>("ecalTPFilterTag");
33 ecalBEFilterTag_ = iConfig.getUntrackedParameter<edm::InputTag>("ecalBEFilterTag");
34 hcalHBHEFilterTag_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHEFilterTag");
35 hcalLaserFilterTag_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalLaserFilterTag");
36
37 photonIsoCalcTag_ = iConfig.getParameter<edm::ParameterSet>("photonIsoCalcTag");
38 }
39
40 ntupleProducer::~ntupleProducer()
41 {
42
43 }
44
45 //
46 // member functions
47 //
48
49 // ------------ method called to for each event ------------
50 void ntupleProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
51 {
52 eventNumber = iEvent.id().event();
53 runNumber = iEvent.id().run();
54 lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
55 bunchCross = (unsigned int)iEvent.bunchCrossing();
56 isRealData = iEvent.isRealData();
57
58 edm::Handle<reco::BeamSpot> beamSpotHandle;
59 iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
60 reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
61
62 beamSpot->SetXYZ(vertexBeamSpot.x0(), vertexBeamSpot.y0(), vertexBeamSpot.z0());
63
64 int vtxCount, jetCount, jptCount, metCount, muCount, pfMuCount, eleCount, photonCount, pfPhotonCount, tauCount, genCount, genPartCount;
65 vtxCount = jetCount = jptCount = metCount = muCount = pfMuCount = eleCount = photonCount = pfPhotonCount = tauCount = genCount = genPartCount = 0;
66
67
68 /////////////////////////////////////
69 // Get PF candidates for later use //
70 /////////////////////////////////////
71
72
73 Handle<PFCandidateCollection> pfCands;
74 iEvent.getByLabel(partFlowTag_,pfCands);
75 const PFCandidateCollection thePfColl = *(pfCands.product());
76
77
78 //////////////////////////
79 //Get vertex information//
80 //////////////////////////
81
82
83 Handle<reco::VertexCollection> primaryVtcs;
84 iEvent.getByLabel(primaryVtxTag_, primaryVtcs);
85
86 for(VertexCollection::const_iterator iVtx = primaryVtcs->begin(); iVtx!= primaryVtcs->end(); ++iVtx){
87 reco::Vertex myVtx = reco::Vertex(*iVtx);
88 if(!myVtx.isValid() || myVtx.isFake()) continue;
89 TCPrimaryVtx* vtxCon = new ((*primaryVtx)[vtxCount]) TCPrimaryVtx;
90 vtxCon->SetXYZ(myVtx.x(), myVtx.y(), myVtx.z());
91 vtxCon->SetNDof(myVtx.ndof());
92 vtxCon->SetChi2(myVtx.chi2());
93 vtxCon->SetNtracks(myVtx.nTracks());
94 vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
95 vtxCon->SetIsFake(myVtx.isFake());
96 ++vtxCount;
97 }
98
99 unsigned ivtx = 0;
100 VertexRef myVtxRef(primaryVtcs, ivtx);
101
102 ///////////////////////
103 //get jet information//
104 ///////////////////////
105
106 Handle<double> rhoCorr;
107 iEvent.getByLabel(rhoCorrTag_, rhoCorr);
108 rhoFactor = (float)(*rhoCorr);
109
110 Handle<double> rho25Corr;
111 iEvent.getByLabel(rho25CorrTag_, rho25Corr);
112 rho25Factor = (float)(*rho25Corr);
113
114 Handle<double> rhoMuCorr;
115 iEvent.getByLabel(rhoMuCorrTag_, rhoMuCorr);
116 rhoMuFactor = (float)(*rhoMuCorr);
117
118 //cout<<" RHOS. In eta 4.4 = "<<rhoFactor<<" in eta25 "<<rho25Factor<<" MUs: "<<rhoMuFactor<<endl;
119
120 if(saveJets_){
121
122 edm::Handle<reco::JetTagCollection> bTagCollectionTCHE;
123 iEvent.getByLabel("trackCountingHighEffBJetTags", bTagCollectionTCHE);
124 const reco::JetTagCollection & bTagsTCHE = *(bTagCollectionTCHE.product());
125
126 Handle<vector<reco::PFJet> > jets;
127 iEvent.getByLabel(jetTag_, jets);
128
129 //Handle<vector<pat::Jet> > jets;
130 //iEvent.getByLabel(jetTag_, jets);
131
132 //for (vector<pat::Jet>::const_iterator iJet = jets->begin(); iJet != jets->end(); ++iJet) {
133
134 for (vector<reco::PFJet>::const_iterator iJet = jets->begin(); iJet != jets->end(); ++iJet) {
135
136 if (iJet->pt() < 10.) continue;
137
138 TCJet* jetCon = new ((*recoJets)[jetCount]) TCJet;
139
140 jetCon->SetPxPyPzE(iJet->px(), iJet->py(), iJet->pz(), iJet->energy());
141 jetCon->SetVtx(0., 0., 0.);
142 jetCon->SetChHadFrac(iJet->chargedHadronEnergyFraction());
143 jetCon->SetNeuHadFrac(iJet->neutralHadronEnergyFraction());
144 jetCon->SetChEmFrac(iJet->chargedEmEnergyFraction());
145 jetCon->SetNeuEmFrac(iJet->neutralEmEnergyFraction());
146 jetCon->SetNumConstit(iJet->chargedMultiplicity() + iJet->neutralMultiplicity());
147 jetCon->SetNumChPart(iJet->chargedMultiplicity());
148
149 //jetCon->SetJetFlavor(iJet->partonFlavour());
150
151 jetCon->SetUncertaintyJES(-1);
152
153 for (tag_iter iTag = bTagsTCHE.begin(); iTag != bTagsTCHE.end(); iTag++) {
154 if (sqrt(pow(iTag->first->eta() - iJet->eta(), 2) + pow(deltaPhi(iTag->first->phi(),iJet->phi()), 2)) == 0.) {
155 jetCon->SetBDiscriminatorMap("TCHE", iTag->second);
156 break;
157 }
158 }
159
160 //jetCon->SetBDiscriminatorMap("TCHE", iJet->bDiscriminator("trackCountingHighEffBJetTags"));
161 //jetCon->SetBDiscriminatorMap("TCHP", iJet->bDiscriminator("trackCountingHighPurBJetTags"));
162 //jetCon->SetBDiscriminatorMap("SSVHE", iJet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
163 //jetCon->SetBDiscriminatorMap("JPB", iJet->bDiscriminator("jetProbabilityBJetTags"));
164
165
166 /////////////////////////
167 // Associate to vertex //
168 /////////////////////////
169
170 associateJetToVertex(*iJet, primaryVtcs, jetCon);
171
172 ++jetCount;
173 }
174
175 // For VBF analysis...
176
177 /*
178 Handle<reco::JPTJetCollection> jptJets;
179 iEvent.getByLabel("ak5JPTJetsL1L2L3", jptJets);
180
181 for (reco::JPTJetCollection::const_iterator iJet = jptJets->begin(); iJet != jptJets->end(); ++iJet) {
182
183 // Perform some pre-cleaning
184 edm::RefToBase<reco::Jet> jptjetRef = iJet->getCaloJetRef();
185 reco::CaloJet const * rawCaloJet = dynamic_cast<reco::CaloJet const*>(&*jptjetRef);
186
187 if (
188 iJet->pt() < 10.
189 || (iJet->chargedEmEnergyFraction() + iJet->neutralEmEnergyFraction()) > 0.01
190 || rawCaloJet->n90() < 2
191 //|| (*jetsID)[(*iJet).getCaloJetRef()].fHPD > 0.98
192 ) continue;
193
194 TCJet* jetCon = new ((*recoJPT)[jptCount]) TCJet;
195
196 jetCon->SetPxPyPzE(iJet->px(), iJet->py(), iJet->pz(), iJet->energy());
197 jetCon->SetVtx(0., 0., 0.);
198
199 //jetCon->SetIDMap("Zch", iJet->getSpecific().Zch);
200 jetCon->SetChHadFrac(iJet->chargedHadronEnergyFraction());
201 jetCon->SetNeuHadFrac(iJet->neutralHadronEnergyFraction());
202 jetCon->SetChEmFrac(iJet->chargedEmEnergyFraction());
203 jetCon->SetNeuEmFrac(iJet->neutralEmEnergyFraction());
204 jetCon->SetNumConstit(iJet->chargedMultiplicity());// + iJet->neutralMultiplicity());
205 jetCon->SetNumChPart(iJet->chargedMultiplicity());
206
207 ++jptCount;
208 }
209 */
210 }
211
212 /////////////
213 // Get MET //
214 /////////////
215
216 if (saveMET_){
217
218 //Handle<vector<pat::MET> > MET;
219 //iEvent.getByLabel(metTag_, MET);
220 //vector<pat::MET>::const_iterator met = MET->begin();
221
222 Handle<PFMETCollection> MET;
223 iEvent.getByLabel(metTag_, MET);
224 PFMETCollection::const_iterator met = MET->begin();
225
226 if (MET->begin() != MET->end()) {
227 recoMET->SetSumEt(met->sumEt());
228 recoMET->SetMagPhi(met->et(), met->phi());
229
230 // PF specififc methods
231 recoMET->SetMuonFraction(met->MuonEtFraction());
232 recoMET->SetNeutralHadronFraction(met->NeutralHadEtFraction());
233 recoMET->SetNeutralEMFraction(met->NeutralEMFraction());
234 recoMET->SetChargedHadronFraction(met->ChargedHadEtFraction());
235 recoMET->SetChargedEMFraction(met->ChargedEMEtFraction());
236
237 }
238 }
239
240 ///////////////
241 // Get muons //
242 ///////////////
243
244
245 if (saveMuons_) {
246
247
248 Handle<vector<reco::Muon> > muons;
249 iEvent.getByLabel(muonTag_, muons);
250
251 for (vector<reco::Muon>::const_iterator iMuon = muons->begin(); iMuon != muons->end(); ++iMuon) {
252 if (!(iMuon->isGlobalMuon() && iMuon->isTrackerMuon())) continue;
253
254 TCMuon* muCon = new ((*recoMuons)[muCount]) TCMuon;
255
256 muCon->SetPxPyPzE(iMuon->px(), iMuon->py(), iMuon->pz(), iMuon->energy());
257 muCon->SetVtx(iMuon->globalTrack()->vx(),iMuon->globalTrack()->vy(),iMuon->globalTrack()->vz());
258 muCon->SetPtError(iMuon->globalTrack()->ptError());
259 muCon->SetCharge(iMuon->charge());
260
261 // Muon ID variables
262 muCon->SetIsPF(iMuon->isPFMuon());
263 muCon->SetIsGLB(iMuon->isGlobalMuon());
264 muCon->SetIsTRK(iMuon->isTrackerMuon());
265 muCon->SetNormalizedChi2(iMuon->globalTrack()->normalizedChi2());
266 muCon->SetCaloComp(iMuon->caloCompatibility());
267 muCon->SetSegComp(muon::segmentCompatibility(*iMuon));
268
269 muCon->SetNumberOfMatches(iMuon->numberOfMatches());
270 muCon->SetNumberOfValidPixelHits(iMuon->globalTrack()->hitPattern().numberOfValidPixelHits());
271 muCon->SetNumberOfValidTrackerHits(iMuon->globalTrack()->hitPattern().numberOfValidTrackerHits());
272 muCon->SetNumberOfValidMuonHits(iMuon->globalTrack()->hitPattern().numberOfValidMuonHits());
273 muCon->SetNumberOfLostPixelHits(iMuon->globalTrack()->hitPattern().numberOfLostPixelHits());
274 muCon->SetNumberOfLostTrackerHits(iMuon->globalTrack()->hitPattern().numberOfLostTrackerHits());
275
276
277 // Set isolation map values
278 // Detector-based isolation
279
280 muCon->SetIsoMap("NTracks_R03", iMuon->isolationR03().nTracks);
281 muCon->SetIsoMap("EmIso_R03", iMuon->isolationR03().emEt);
282 muCon->SetIsoMap("HadIso_R03", iMuon->isolationR03().hadEt);
283 muCon->SetIsoMap("SumPt_R03", iMuon->isolationR03().sumPt);
284
285 muCon->SetIsoMap("NTracks_R05", iMuon->isolationR05().nTracks);
286 muCon->SetIsoMap("EmIso_R05", iMuon->isolationR05().emEt);
287 muCon->SetIsoMap("HadIso_R05", iMuon->isolationR05().hadEt);
288 muCon->SetIsoMap("SumPt_R05", iMuon->isolationR05().sumPt);
289
290 // PF-based isolation
291 muCon->SetIsoMap("pfChargedPt_R03", iMuon->pfIsolationR03().sumChargedParticlePt);
292 muCon->SetIsoMap("pfChargedHadronPt_R03", iMuon->pfIsolationR03().sumChargedHadronPt);
293 muCon->SetIsoMap("pfPhotonEt_R03", iMuon->pfIsolationR03().sumPhotonEt);
294 muCon->SetIsoMap("pfNeutralHadronEt_R03", iMuon->pfIsolationR03().sumNeutralHadronEt);
295 muCon->SetIsoMap("pfPUPt_R03", iMuon->pfIsolationR03().sumPUPt);
296
297 muCon->SetIsoMap("pfChargedPt_R04", iMuon->pfIsolationR04().sumChargedParticlePt);
298 muCon->SetIsoMap("pfChargedHadronPt_R04", iMuon->pfIsolationR04().sumChargedHadronPt);
299 muCon->SetIsoMap("pfPhotonEt_R04", iMuon->pfIsolationR04().sumPhotonEt);
300 muCon->SetIsoMap("pfNeutralHadronEt_R04", iMuon->pfIsolationR04().sumNeutralHadronEt);
301 muCon->SetIsoMap("pfPUPt_R04", iMuon->pfIsolationR04().sumPUPt);
302
303
304 muCount++;
305 }
306 }
307
308
309 ///////////////////
310 // Get electrons //
311 ///////////////////
312
313
314 edm::Handle<reco::ConversionCollection> hConversions;
315 iEvent.getByLabel("allConversions", hConversions);
316
317 if (saveElectrons_) {
318
319 Handle<reco::GsfElectronCollection > electrons;
320 iEvent.getByLabel(electronTag_, electrons);
321
322 for (vector<reco::GsfElectron>::const_iterator iElectron = electrons->begin(); iElectron != electrons->end(); ++iElectron) {
323 if (iElectron->pt() < 10) continue;
324
325 TCElectron* eleCon = new ((*recoElectrons)[eleCount]) TCElectron;
326
327 // Basic physics object info
328 eleCon->SetPxPyPzE(iElectron->px(), iElectron->py(), iElectron->pz(), iElectron->energy());
329 eleCon->SetVtx(iElectron->gsfTrack()->vx(),iElectron->gsfTrack()->vy(),iElectron->gsfTrack()->vz());
330 eleCon->SetCharge(iElectron->charge());
331
332 // Fiducial variables
333 eleCon->SetIsEB(iElectron->isEB());
334 eleCon->SetIsEE(iElectron->isEE());
335 eleCon->SetIsInGap(iElectron->isGap());
336
337
338 // Electron ID variables
339 eleCon->SetHadOverEm(iElectron->hadronicOverEm());
340 eleCon->SetDphiSuperCluster(iElectron->deltaPhiSuperClusterTrackAtVtx());
341 eleCon->SetDetaSuperCluster(iElectron->deltaEtaSuperClusterTrackAtVtx());
342 eleCon->SetSigmaIetaIeta(iElectron->sigmaIetaIeta());
343 eleCon->SetFBrem(iElectron->fbrem());
344 eleCon->SetEOverP(iElectron->eSuperClusterOverP());
345 eleCon->SetSCEta(iElectron->superCluster()->eta());
346 eleCon->SetR9(iElectron->r9());
347
348 eleCon->SetPtError(iElectron->gsfTrack()->ptError());
349 eleCon->SetNormalizedChi2(iElectron->gsfTrack()->normalizedChi2());
350
351 eleCon->SetNumberOfValidPixelHits(iElectron->gsfTrack()->hitPattern().numberOfValidPixelHits());
352 eleCon->SetNumberOfValidTrackerHits(iElectron->gsfTrack()->hitPattern().numberOfValidTrackerHits());
353 eleCon->SetNumberOfLostPixelHits(iElectron->gsfTrack()->hitPattern().numberOfLostPixelHits());
354 eleCon->SetNumberOfLostTrackerHits(iElectron->gsfTrack()->hitPattern().numberOfLostTrackerHits());
355
356 eleCon->SetIdMap("fabsEPDiff",fabs((1/iElectron->ecalEnergy()) - (1/iElectron->trackMomentumAtVtx().R())));
357
358 // Electron Iso variables
359 eleCon->SetIsoMap("EmIso_R03", iElectron->dr03EcalRecHitSumEt());
360 eleCon->SetIsoMap("HadIso_R03", iElectron->dr03HcalTowerSumEt());
361 eleCon->SetIsoMap("SumPt_R03", iElectron->dr03TkSumPt());
362
363 eleCon->SetIsoMap("EmIso_R04", iElectron->dr04EcalRecHitSumEt());
364 eleCon->SetIsoMap("HadIso_R04", iElectron->dr04HcalTowerSumEt());
365 eleCon->SetIsoMap("SumPt_R04", iElectron->dr04TkSumPt());
366
367 eleCon->SetIsoMap("pfPhotonEt_R03", iElectron->pfIsolationVariables().photonIso);
368 eleCon->SetIsoMap("pfChargedHadron_R03", iElectron->pfIsolationVariables().chargedHadronIso);
369 eleCon->SetIsoMap("pfNeutralHadron_R03", iElectron->pfIsolationVariables().neutralHadronIso);
370
371 eleIsolator.fGetIsolation(&(*iElectron), &thePfColl, myVtxRef, primaryVtcs);
372 eleCon->SetIsoMap("pfChIso_R04",eleIsolator.getIsolationCharged());
373 eleCon->SetIsoMap("pfNeuIso_R04",eleIsolator.getIsolationNeutral());
374 eleCon->SetIsoMap("pfPhoIso_R04",eleIsolator.getIsolationPhoton());
375
376 // Effective area for rho PU corrections (not sure if needed)
377 float AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, iElectron->eta(), ElectronEffectiveArea::kEleEAData2012);
378 float AEff04 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso04, iElectron->eta(), ElectronEffectiveArea::kEleEAData2012);
379 eleCon->SetIsoMap("EffArea_R03", AEff03);
380 eleCon->SetIsoMap("EffArea_R04", AEff04);
381
382 // Conversion information
383 bool convVeto = !(ConversionTools::hasMatchedConversion(*iElectron,hConversions,vertexBeamSpot.position()));
384 eleCon->SetConversionVeto(convVeto);
385 eleCon->SetConversionMissHits(iElectron->gsfTrack()->trackerExpectedHitsInner().numberOfHits());
386
387 eleCon->SetIsoMap("fabsEPDiff",fabs((1/iElectron->ecalEnergy()) - (1/iElectron->trackMomentumAtVtx().R()))); //ooemoop
388
389 // EID maps for VBTF working points -- probably not needed anymore
390 //eleCon->SetCutLevel(iElectron->electronID("eidVBTF95"), 95);
391 //eleCon->SetCutLevel(iElectron->electronID("eidVBTF90"), 90);
392 //eleCon->SetCutLevel(iElectron->electronID("eidVBTF85"), 85);
393 //eleCon->SetCutLevel(iElectron->electronID("eidVBTF80"), 80);
394 //eleCon->SetCutLevel(iElectron->electronID("eidVBTF70"), 70);
395 //eleCon->SetCutLevel(iElectron->electronID("eidVBTF60"), 60);
396
397 // Add electron MVA ID and ISO when done -- needs work
398 //electronMVA(vtxCollection, iElectron);
399
400 eleCount++;
401 }
402 }
403
404
405 /////////////////
406 // Get photons //
407 /////////////////
408
409
410 if (savePhotons_) {
411
412 Handle<vector<reco::Photon> > photons;
413 iEvent.getByLabel(photonTag_, photons);
414
415 edm::Handle<reco::GsfElectronCollection> hElectrons;
416 iEvent.getByLabel("gsfElectrons", hElectrons);
417
418 for (vector<reco::Photon>::const_iterator iPhoton = photons->begin(); iPhoton != photons->end() ; ++iPhoton) {
419
420 TCPhoton* myPhoton = new ((*recoPhotons)[photonCount]) TCPhoton();
421 myPhoton->SetPxPyPzE(iPhoton->px(), iPhoton->py(), iPhoton->pz(), iPhoton->p());
422 myPhoton->SetVtx(iPhoton->vx(), iPhoton->vy(), iPhoton->vz());
423
424 // ID variables
425 myPhoton->SetHadOverEm(iPhoton->hadronicOverEm());
426 myPhoton->SetSigmaIEtaIEta(iPhoton->sigmaIetaIeta());
427 myPhoton->SetR9(iPhoton->r9());
428 myPhoton->SetTrackVeto(iPhoton->hasPixelSeed());
429
430 myPhoton->SetSCEta(iPhoton->superCluster()->eta());
431 myPhoton->SetSCPhi(iPhoton->superCluster()->phi());
432 myPhoton->SetSCEnergy(iPhoton->superCluster()->energy());
433
434 // detector-based isolation
435 myPhoton->SetIsoMap("EmIso_R03", (iPhoton->ecalRecHitSumEtConeDR03()));
436 myPhoton->SetIsoMap("HadIso_R03", (iPhoton->hcalTowerSumEtConeDR03()));
437 myPhoton->SetIsoMap("TrkIso_R03", (iPhoton->trkSumPtHollowConeDR03()));
438
439 myPhoton->SetIsoMap("EmIso_R04", (iPhoton->ecalRecHitSumEtConeDR04()));
440 myPhoton->SetIsoMap("HadIso_R04", (iPhoton->hcalTowerSumEtConeDR04()));
441 myPhoton->SetIsoMap("TrkIso_R04", (iPhoton->trkSumPtHollowConeDR04()));
442
443 // PF Iso for photons
444 phoIsolator.fGetIsolation(&(*iPhoton),&thePfColl, myVtxRef, primaryVtcs);
445 myPhoton->SetIsoMap("chIso03",phoIsolator.getIsolationCharged());
446 myPhoton->SetIsoMap("nhIso03",phoIsolator.getIsolationNeutral());
447 myPhoton->SetIsoMap("phIso03",phoIsolator.getIsolationPhoton());
448
449 // Hcal isolation for 2012
450 //myPhoton->SetIsoMap("HadIso_R03",iPhoton->hcalTowerSumEtConeDR03() +
451 // (iPhoton->hadronicOverEm() - iPhoton->hadTowOverEm())*iPhoton->superCluster()->energy()/cosh(iPhoton->superCluster()->eta()));
452 //myPhoton->SetIsoMap("HadIso_R04",iPhoton->hcalTowerSumEtConeDR04() +
453 // (iPhoton->hadronicOverEm() - iPhoton->hadTowOverEm())*iPhoton->superCluster()->energy()/cosh(iPhoton->superCluster()->eta()));
454
455
456 //Conversion info
457 bool passElectronVeto = !(ConversionTools::hasMatchedPromptElectron(iPhoton->superCluster(), hElectrons, hConversions, vertexBeamSpot.position()));
458 myPhoton->SetConversionVeto(passElectronVeto);
459
460 ++photonCount;
461 }
462 }
463
464
465 //////////////
466 // Get taus //
467 //////////////
468
469
470 if (saveTaus_) {
471
472 Handle<vector<pat::Tau> > taus;
473 iEvent.getByLabel(tauTag_, taus);
474
475 for (vector<pat::Tau>::const_iterator iTau = taus->begin(); iTau != taus->end(); ++iTau) {
476
477 if (!iTau->isPFTau()
478 or iTau->signalPFChargedHadrCands().size() < 1
479 or iTau->pt() < 10
480 ) continue;
481
482 TCTau* tauCon = new ((*recoTaus)[tauCount]) TCTau;
483
484 tauCon->SetNChHad(iTau->signalPFChargedHadrCands().size());
485 tauCon->SetNGamma (iTau->signalPFGammaCands().size());
486 tauCon->SetNNeutrHad (iTau->signalPFNeutrHadrCands().size());
487 tauCon->SetCharge(iTau->charge());
488 tauCon->SetDecayMode(iTau->decayMode());
489
490 tauCon->SetPxPyPzE(iTau->px(),iTau->py(),iTau->pz(),iTau->energy());
491 tauCon->SetCharge(iTau->charge());
492
493 if (iTau->leadPFChargedHadrCand()->trackRef().isNonnull()) {
494 tauCon->SetLeadChHadP4(iTau->leadPFChargedHadrCand()->px(),
495 iTau->leadPFChargedHadrCand()->py(),
496 iTau->leadPFChargedHadrCand()->pz(),
497 iTau->leadPFChargedHadrCand()->energy());
498
499 tauCon->SetPositionFromTrack(iTau->leadPFChargedHadrCand()->trackRef()->vx(),
500 iTau->leadPFChargedHadrCand()->trackRef()->vy(),
501 iTau->leadPFChargedHadrCand()->trackRef()->vz());
502 }
503
504
505 if (iTau->signalPFGammaCands().size()+iTau->signalPFNeutrHadrCands().size()>0)
506 tauCon->SetLeadNeutrP4(iTau->leadPFNeutralCand()->px(),
507 iTau->leadPFNeutralCand()->py(),
508 iTau->leadPFNeutralCand()->pz(),
509 iTau->leadPFNeutralCand()->energy());
510
511 tauCon->SetPositionFromTau(iTau->vx(),iTau->vy(), iTau->vz());
512
513 tauCon->SetIsoGammaEtSum(iTau->isolationPFGammaCandsEtSum());
514 tauCon->SetIsoChHadPtSum(iTau->isolationPFChargedHadrCandsPtSum());
515
516
517
518 // set the discriminators
519 // note that the strings for PAT and RECO are different. The names of the TCTau accessors are set following the RECO names
520 // the "mapping" is taken from tauTools.py
521
522 tauCon->SetHpsPFTauDiscriminationByDecayModeFinding(iTau->tauID("decayModeFinding")); // "DiscriminationByDecayModeFinding"
523
524 // isolation
525 tauCon->SetHpsPFTauDiscriminationByVLooseIsolation(iTau->tauID("byVLooseIsolation")); // "DiscriminationByVLooseIsolation"
526 tauCon->SetHpsPFTauDiscriminationByLooseIsolation (iTau->tauID("byLooseIsolation")); // "DiscriminationByLooseIsolation"
527 tauCon->SetHpsPFTauDiscriminationByMediumIsolation(iTau->tauID("byMediumIsolation")); // "DiscriminationByMediumIsolation"
528 tauCon->SetHpsPFTauDiscriminationByTightIsolation (iTau->tauID("byTightIsolation")); // "DiscriminationByTightIsolation"
529
530 // isolation with corrections
531 tauCon->SetHpsPFTauDiscriminationByVLooseIsolationDBSumPtCorr
532 (iTau->tauID("byVLooseIsolationDeltaBetaCorr")); // "DiscriminationByVLooseIsolationDBSumPtCorr"
533 tauCon->SetHpsPFTauDiscriminationByLooseIsolationDBSumPtCorr
534 (iTau->tauID("byLooseIsolationDeltaBetaCorr")); // "DiscriminationByLooseIsolationDBSumPtCorr"
535 tauCon->SetHpsPFTauDiscriminationByMediumIsolationDBSumPtCorr
536 (iTau->tauID("byMediumIsolationDeltaBetaCorr")); // "DiscriminationByMediumIsolationDBSumPtCorr"
537 tauCon->SetHpsPFTauDiscriminationByTightIsolationDBSumPtCorr
538 (iTau->tauID("byTightIsolationDeltaBetaCorr")); // "DiscriminationByTightIsolationDBSumPtCorr"
539
540 // combined isolation with corrections
541 tauCon->SetHpsPFTauDiscriminationByVLooseCombinedIsolationDBSumPtCorr
542 (iTau->tauID("byVLooseCombinedIsolationDeltaBetaCorr")); // "DiscriminationByVLooseCombinedIsolationDBSumPtCorr"
543 tauCon->SetHpsPFTauDiscriminationByLooseCombinedIsolationDBSumPtCorr
544 (iTau->tauID("byLooseCombinedIsolationDeltaBetaCorr")); // "DiscriminationByLooseCombinedIsolationDBSumPtCorr"
545 tauCon->SetHpsPFTauDiscriminationByMediumCombinedIsolationDBSumPtCorr
546 (iTau->tauID("byMediumCombinedIsolationDeltaBetaCorr")); // "DiscriminationByMediumCombinedIsolationDBSumPtCorr"
547 tauCon->SetHpsPFTauDiscriminationByTightCombinedIsolationDBSumPtCorr
548 (iTau->tauID("byTightCombinedIsolationDeltaBetaCorr")); // "DiscriminationByTightCombinedIsolationDBSumPtCorr"
549
550 // anti e/mu discriminators
551 tauCon->SetHpsPFTauDiscriminationAgainstElectronLoose (iTau->tauID("againstElectronLoose")); // "DiscriminationByLooseElectronRejection"
552 tauCon->SetHpsPFTauDiscriminationAgainstElectronMedium(iTau->tauID("againstElectronMedium")); // "DiscriminationByMediumElectronRejection"
553 tauCon->SetHpsPFTauDiscriminationAgainstElectronTight (iTau->tauID("againstElectronTight")); // "DiscriminationByTightElectronRejection"
554
555 tauCon->SetHpsPFTauDiscriminationAgainstMuonLoose (iTau->tauID("againstMuonLoose")); // "DiscriminationByLooseMuonRejection")
556 // tauCon->SetHpsPFTauDiscriminationAgainstMuonMediumt(iTau->tauID("againstMuonMedium")); // "DiscriminationByMediumMuonRejection" <- not in python
557 tauCon->SetHpsPFTauDiscriminationAgainstMuonTight (iTau->tauID("againstMuonTight")); // "DiscriminationByTightMuonRejection"
558
559 tauCount++;
560 }
561 }
562
563
564 ////////////////////////
565 // Get gen-level info //
566 ////////////////////////
567
568
569 if (!isRealData) {
570
571 Handle<GenEventInfoProduct> GenEventInfoHandle;
572 iEvent.getByLabel("generator", GenEventInfoHandle);
573
574 evtWeight = ptHat = qScale = -1;
575
576 if (GenEventInfoHandle.isValid()) {
577 //qScale = GenEventInfoHandle->qScale();
578 //evtWeight = GenEventInfoHandle->weight();
579 ptHat = (GenEventInfoHandle->hasBinningValues() ? GenEventInfoHandle->binningValues()[0] : 0.0);
580 }
581
582
583 ////////////////////
584 // PU information //
585 ////////////////////
586
587 Handle<std::vector< PileupSummaryInfo > > PUInfo;
588 iEvent.getByLabel(edm::InputTag("addPileupInfo"), PUInfo);
589 std::vector<PileupSummaryInfo>::const_iterator iPV;
590
591 for(iPV = PUInfo->begin(); iPV != PUInfo->end(); ++iPV){
592 if (iPV->getBunchCrossing() == 0){
593 nPUVertices = iPV->getPU_NumInteractions();
594 nPUVerticesTrue = iPV->getTrueNumInteractions();
595 }
596 }
597
598 //////////////////////
599 // Get genParticles //
600 //////////////////////
601
602 if (saveGenParticles_) {
603 Handle<GenParticleCollection> genParticleColl;
604 iEvent.getByLabel("genParticles", genParticleColl);
605
606 for (GenParticleCollection::const_iterator iGenPart = genParticleColl->begin(); iGenPart != genParticleColl->end(); ++iGenPart) {
607 const reco::GenParticle myParticle = reco::GenParticle(*iGenPart);
608
609 //// Leptons and photons and b's, (oh my)
610 if (
611 myParticle.pt() > 8
612 && (
613 (abs(myParticle.pdgId()) >= 11 && abs(myParticle.pdgId()) <= 16)
614 || myParticle.pdgId() == 22
615 || abs(myParticle.pdgId()) == 5
616 )
617 ) {
618
619 TCGenParticle* genCon = new ((*genParticles)[genPartCount]) TCGenParticle;
620 genCon->SetPxPyPzE(myParticle.px(), myParticle.py(), myParticle.pz(), myParticle.energy() );
621 genCon->SetVtx(myParticle.vx(), myParticle.vy(), myParticle.vz());
622 genCon->SetCharge(myParticle.charge());
623 genCon->SetPDGId(myParticle.pdgId());
624 genCon->SetMother(myParticle.mother()->pdgId());
625 genCon->SetStatus(myParticle.status());
626 if (myParticle.mother()->numberOfMothers() != 0) genCon->SetGrandmother(myParticle.mother()->mother()->pdgId());
627 ++genPartCount;
628 }
629
630 //// Z's, W's, H's, and now big juicy Gravitons
631 if (
632 abs(myParticle.pdgId()) == 23
633 || abs(myParticle.pdgId()) == 24
634 || abs(myParticle.pdgId()) == 25
635 || abs(myParticle.pdgId()) == 35
636 || abs(myParticle.pdgId()) == 36
637 || abs(myParticle.pdgId()) == 39
638 ){
639
640
641 TCGenParticle* genCon = new ((*genParticles)[genPartCount]) TCGenParticle;
642 genCon->SetPxPyPzE(myParticle.px(), myParticle.py(), myParticle.pz(), myParticle.energy() );
643 genCon->SetVtx(myParticle.vx(), myParticle.vy(), myParticle.vz() );
644 genCon->SetCharge(myParticle.charge());
645 genCon->SetPDGId(myParticle.pdgId());
646 genCon->SetMother(myParticle.mother()->pdgId());
647 genCon->SetStatus(myParticle.status());
648 ++genPartCount;
649 }
650 }
651 }
652
653
654 /////////////////
655 // Get genJets //
656 /////////////////
657
658 if (saveGenJets_) {
659
660 Handle<reco::GenJetCollection> GenJets;
661 iEvent.getByLabel(genJetTag_, GenJets);
662
663 for (GenJetCollection::const_iterator iJet = GenJets->begin(); iJet!= GenJets->end(); ++iJet) {
664 reco::GenJet myJet = reco::GenJet(*iJet);
665 if (myJet.pt() > 10) {
666 TCGenJet* jetCon = new ((*genJets)[genCount]) TCGenJet;
667 jetCon->SetPxPyPzE(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
668 jetCon->SetHadEnergy(myJet.hadEnergy());
669 jetCon->SetEmEnergy(myJet.emEnergy());
670 jetCon->SetInvEnergy(myJet.invisibleEnergy());
671 jetCon->SetAuxEnergy(myJet.auxiliaryEnergy());
672 jetCon->SetNumConstit(myJet.getGenConstituents().size());
673 jetCon->SetJetFlavor(0);
674 }
675 ++genCount;
676 }
677 }
678 }
679
680
681 ///////////////////
682 // Noise filters //
683 ///////////////////
684
685 //if (isRealData) {
686
687 myNoiseFilters.isScraping = isFilteredOutScraping(iEvent, iSetup, 10, 0.25);
688
689 Handle<bool> hcalNoiseFilterHandle;
690 iEvent.getByLabel(hcalHBHEFilterTag_, hcalNoiseFilterHandle);
691 if (hcalNoiseFilterHandle.isValid()) myNoiseFilters.isNoiseHcalHBHE = !(Bool_t)(*hcalNoiseFilterHandle);
692 else LogWarning("Filters")<<"hcal noise NOT valid ";
693
694 Handle<bool> hcalLaserFilterHandle;
695 iEvent.getByLabel(hcalLaserFilterTag_, hcalLaserFilterHandle);
696 if (hcalLaserFilterHandle.isValid()) myNoiseFilters.isNoiseHcalLaser = !(Bool_t)(*hcalLaserFilterHandle);
697 else LogWarning("Filters")<<"hcal Laser NOT valid ";
698
699 Handle<bool> ecalTPFilterHandle;
700 iEvent.getByLabel(ecalTPFilterTag_, ecalTPFilterHandle);
701 if (ecalTPFilterHandle.isValid()) myNoiseFilters.isNoiseEcalTP = !(Bool_t)(*ecalTPFilterHandle);
702 else LogWarning("Filters")<<"Ecal TP NOT valid ";
703
704 Handle<bool> ecalBEFilterHandle;
705 iEvent.getByLabel(ecalBEFilterTag_, ecalBEFilterHandle);
706 if (ecalBEFilterHandle.isValid()) myNoiseFilters.isNoiseEcalBE = !(Bool_t)(*ecalBEFilterHandle);
707 else LogWarning("Filters")<<"Ecal BE NOT valid ";
708
709 edm::Handle<BeamHaloSummary> TheBeamHaloSummary;
710 iEvent.getByLabel("BeamHaloSummary",TheBeamHaloSummary);
711 const BeamHaloSummary TheSummary = (*TheBeamHaloSummary.product() );
712 if (!TheBeamHaloSummary.isValid()) LogWarning("Filters")<<"The Summary (for CSC halo) NOT valid ";
713
714 myNoiseFilters.isCSCTightHalo = TheSummary.CSCTightHaloId();
715 myNoiseFilters.isCSCLooseHalo = TheSummary.CSCLooseHaloId();
716
717
718 //LogWarning("Filters")<<"\n csc1 "<< myNoiseFilters.isCSCTightHalo<<" csc2 "<<myNoiseFilters.isCSCLooseHalo
719 // <<" isNoiseHcal HBHE "<<myNoiseFilters.isNoiseHcalHBHE<<" laser "<<myNoiseFilters.isNoiseHcalLaser<<"\n"
720 // <<" ecal TP "<<myNoiseFilters.isNoiseEcalTP<<" ecal BE "<<myNoiseFilters.isNoiseEcalBE;
721
722 //}
723
724 ////////////////////////////
725 // get trigger information//
726 ////////////////////////////
727
728 edm::Handle<TriggerResults> hltR;
729 triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltProcess_);
730 iEvent.getByLabel(triggerResultsTag_,hltR);
731
732 const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
733 hlNames = triggerNames.triggerNames();
734
735 triggerStatus = 0x0;
736
737 for (int i=0; i < (int)hlNames.size(); ++i) {
738 if (!triggerDecision(hltR, i)) continue;
739
740 for (int j = 0; j < (int)triggerPaths_.size(); ++j){
741 if (triggerPaths_[j] == "") continue;
742
743 if (hlNames[i].compare(0, triggerPaths_[j].length(),triggerPaths_[j]) == 0) {
744 //cout << hlNames[i] << " ?= " << triggerPaths_[j] << endl;
745 triggerStatus |= ULong64_t(0x01) << j;
746
747 if (isRealData) {
748 pair<int, int> preScales;
749 preScales = hltConfig_.prescaleValues(iEvent, iSetup, hlNames[i]);
750 hltPrescale[j] = preScales.first*preScales.second;
751 } else {
752 hltPrescale[j] = 1;
753 }
754 }
755 }
756 }
757
758 ++nEvents;
759
760 if (eleCount > 0 || muCount > 0) eventTree -> Fill(); // possibly specify a cut in configuration
761
762 primaryVtx -> Clear("C");
763 recoJets -> Clear("C");
764 recoJPT -> Clear("C");
765 recoMuons -> Clear("C");
766 recoElectrons -> Clear("C");
767 recoTaus -> Clear("C");
768 recoPhotons -> Clear("C");
769 //pfPhotons -> Clear("C");
770 triggerObjects-> Clear("C");
771 genJets -> Clear("C");
772 genParticles -> Clear("C");
773 }
774
775 // ------------ method called once each job just before starting event loop ------------
776 void ntupleProducer::beginJob()
777 {
778 eventTree = fs->make<TTree>("eventTree","eventTree");
779 runTree = fs->make<TTree>("runTree","runTree");
780 jobTree = fs->make<TTree>("jobTree", "jobTree");
781
782 primaryVtx = new TClonesArray("TCPrimaryVtx");
783 recoJets = new TClonesArray("TCJet");
784 recoJPT = new TClonesArray("TCJet");
785 recoElectrons = new TClonesArray("TCElectron");
786 recoMuons = new TClonesArray("TCMuon");
787 recoTaus = new TClonesArray("TCTau");
788 recoPhotons = new TClonesArray("TCPhoton");
789 //pfPhotons = new TClonesArray("TCPhoton");
790 triggerObjects = new TClonesArray("TCTriggerObject");
791 genJets = new TClonesArray("TCGenJet");
792 genParticles = new TClonesArray("TCGenParticle");
793 beamSpot = new TVector3();
794 recoMET = 0;
795
796 eventTree->Branch("recoJets",&recoJets, 6400, 0);
797 eventTree->Branch("recoJPT",&recoJPT, 6400, 0);
798 eventTree->Branch("recoElectrons",&recoElectrons, 6400, 0);
799 eventTree->Branch("recoMuons",&recoMuons, 6400, 0);
800 eventTree->Branch("recoTaus",&recoTaus, 6400, 0);
801 eventTree->Branch("recoPhotons",&recoPhotons, 6400, 0);
802 //eventTree->Branch("pfPhotons",&pfPhotons, 6400, 0);
803 eventTree->Branch("recoMET", &recoMET, 6400, 0);
804 eventTree->Branch("triggerObjects", &triggerObjects, 6400, 0);
805 eventTree->Branch("genJets",&genJets, 6400, 0);
806 eventTree->Branch("genParticles",&genParticles, 6400, 0);
807
808 eventTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
809 eventTree->Branch("beamSpot", &beamSpot, 6400, 0);
810 eventTree->Branch("nPUVertices", &nPUVertices, "nPUVertices/I");
811 eventTree->Branch("nPUVerticesTrue", &nPUVerticesTrue, "nPUVerticesTrue/F");
812
813 eventTree->Branch("isRealData",&isRealData, "isRealData/O");
814 eventTree->Branch("runNumber",&runNumber, "runNumber/i");
815 eventTree->Branch("eventNumber",&eventNumber, "eventNumber/l");
816 eventTree->Branch("lumiSection",&lumiSection, "lumiSection/i");
817 eventTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
818
819 eventTree->Branch("ptHat",&ptHat, "ptHat/F");
820 eventTree->Branch("qScale", &qScale, "qScale/F");
821 eventTree->Branch("evtWeight", &evtWeight, "evtWeight/F");
822 eventTree->Branch("rhoFactor",&rhoFactor, "rhoFactor/F");
823 eventTree->Branch("rho25Factor",&rho25Factor, "rho25Factor/F");
824 eventTree->Branch("rhoMuFactor",&rhoMuFactor, "rhoMuFactor/F");
825 eventTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/l");
826 eventTree->Branch("hltPrescale",hltPrescale, "hltPrescale[64]/i");
827
828 eventTree->Branch("NoiseFilters", &myNoiseFilters.isScraping, "isScraping/O:isNoiseHcalHBHE:isNoiseHcalLaser:isNoiseEcalTP:isNoiseEcalBE:isCSCTightHalo:isCSCLooseHalo");
829
830 runTree->Branch("deliveredLumi",&deliveredLumi, "deliveredLumi/F");
831 runTree->Branch("recordedLumi",&recordedLumi, "recordedLumi/F");
832 runTree->Branch("runNumber",&runNumber, "runNumber/i");
833
834 jobTree->Branch("nEvents",&nEvents, "nEvents/i");
835 jobTree->Branch("triggerNames", "vector<string>", &triggerPaths_);
836
837 // Initialize HLT prescales //
838
839 for (int i = 0; i < (int)(sizeof(hltPrescale)/sizeof(int)); ++i) hltPrescale[i] = 1;
840
841 // Start counting number of events per job //
842 nEvents = 0;
843
844 // Photon Iso maker init
845 phoIsolator.initializePhotonIsolation(kTRUE);
846 phoIsolator.setConeSize(0.3);
847
848 // Initialize Electron MVA nonsense
849 eleIsolator.initializeElectronIsolation(kTRUE);
850 eleIsolator.setConeSize(0.4);
851
852 }
853
854 void ntupleProducer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
855 {
856 bool changed = true;
857 hltConfig_.init(iRun, iSetup, hltProcess_, changed);
858 deliveredLumi = 0;
859 recordedLumi = 0;
860 }
861
862 void ntupleProducer::endLuminosityBlock(const edm::LuminosityBlock& iLumi, const edm::EventSetup& iSetup)
863 {
864 //if (isRealData) {
865 if (false) {
866 edm::Handle<LumiSummary> lumiSummary;
867 iLumi.getByLabel("lumiProducer", lumiSummary);
868
869 deliveredLumi += lumiSummary->avgInsDelLumi()*93.244;
870 recordedLumi += deliveredLumi*lumiSummary->liveFrac();
871 }
872 }
873
874
875 void ntupleProducer::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup)
876 {
877 runTree->Fill();
878 }
879
880
881 void ntupleProducer::endJob()
882 {
883 cout<<nEvents<<endl;
884 jobTree->Fill();
885 }
886
887
888 bool ntupleProducer::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
889 {
890 bool triggerPassed = false;
891 if(hltR->wasrun(iTrigger) &&
892 hltR->accept(iTrigger) &&
893 !hltR->error(iTrigger) ){
894 triggerPassed = true;
895 }
896 return triggerPassed;
897 }
898
899
900 float ntupleProducer::sumPtSquared(const Vertex & v)
901 {
902 float sum = 0.;
903 float pT;
904 for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
905 pT = (**it).pt();
906 float epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
907
908 sum += pT*pT;
909 }
910 return sum;
911 }
912
913
914 bool ntupleProducer::isFilteredOutScraping( const edm::Event& iEvent, const edm::EventSetup& iSetup, int numtrack, double thresh)
915 {
916
917 bool accepted = false;
918 float fraction = 0;
919 // get GeneralTracks collection
920
921 edm::Handle<reco::TrackCollection> tkRef;
922 iEvent.getByLabel("generalTracks",tkRef);
923 const reco::TrackCollection* tkColl = tkRef.product();
924
925 int numhighpurity=0;
926 reco::TrackBase::TrackQuality _trackQuality = reco::TrackBase::qualityByName("highPurity");
927
928 if(tkColl->size()>(UInt_t)numtrack){
929 reco::TrackCollection::const_iterator itk = tkColl->begin();
930 reco::TrackCollection::const_iterator itk_e = tkColl->end();
931 for(;itk!=itk_e;++itk){
932 if(itk->quality(_trackQuality)) numhighpurity++;
933 }
934 fraction = (float)numhighpurity/(float)tkColl->size();
935 if(fraction>thresh) accepted=true;
936 } else {
937 //if less than 10 Tracks accept the event anyway
938 accepted= true;
939 }
940 return !accepted; //if filtered out it's not accepted.
941 }
942
943
944 bool ntupleProducer::associateJetToVertex(reco::PFJet inJet, Handle<reco::VertexCollection> vtxCollection, TCJet *outJet)
945 {
946 if(fabs(inJet.eta()) > 2.5){
947 outJet->SetVtxSumPtFrac(-1);
948 outJet->SetVtxSumPt(-1);
949 outJet->SetVtxTrackFrac(-1);
950 outJet->SetVtxNTracks(-1);
951 outJet->SetVtxSumPtIndex(0);
952 outJet->SetVtxCountIndex(0);
953
954 return false;
955 }
956
957 vector<float> associatedTrackSumPt;
958 vector<float> associatedTrackCount;
959 vector<const reco::Track*> jetTracks;
960 float sumTrackX, sumTrackY, sumTrackZ, sumTrackPt;
961 int nJetTracks = 0;
962 int vCount = 0;
963
964 sumTrackX = sumTrackY = sumTrackZ = sumTrackPt = 0;
965
966 //const reco::TrackRefVector &tracks = inJet.associatedTracks();
967 const reco::TrackRefVector &tracks = inJet.getTrackRefs();
968
969 for (TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack) {
970 const reco::Track &jetTrack = **iTrack;
971
972 sumTrackPt += jetTrack.pt();
973 sumTrackX += jetTrack.vx();
974 sumTrackY += jetTrack.vy();
975 sumTrackZ += jetTrack.vz();
976 jetTracks.push_back(&jetTrack);
977 ++nJetTracks;
978 }
979
980 if(jetTracks.size() == 0){
981 outJet->SetVtxSumPtFrac(-1);
982 outJet->SetVtxSumPt(0);
983 outJet->SetVtxTrackFrac(-1);
984 outJet->SetVtxNTracks(0);
985 outJet->SetVtxSumPtIndex(0);
986 outJet->SetVtxCountIndex(0);
987 outJet->SetVtx(0., 0., 0.);
988 } else {
989 outJet->SetVtx(sumTrackX/nJetTracks, sumTrackY/nJetTracks, sumTrackZ/nJetTracks);
990
991 for (VertexCollection::const_iterator iVtx = vtxCollection->begin(); iVtx!= vtxCollection->end(); ++iVtx) {
992 reco::Vertex myVtx = reco::Vertex(*iVtx);
993 if(!myVtx.isValid() || myVtx.isFake()) continue;
994 associatedTrackSumPt.push_back(0);
995 associatedTrackCount.push_back(0);
996
997 for(Vertex::trackRef_iterator iTrackRef = myVtx.tracks_begin(); iTrackRef != myVtx.tracks_end(); ++iTrackRef){
998 const edm::RefToBase<reco::Track> &myTrackRef = *iTrackRef;
999
1000 if(myTrackRef.isAvailable()){
1001 const reco::Track &myVertexTrack = *myTrackRef.get();
1002
1003 for(vector<const reco::Track*>::const_iterator iTrack = jetTracks.begin(); iTrack != jetTracks.end(); ++iTrack){
1004 if (*iTrack == &myVertexTrack) {
1005 associatedTrackSumPt.at(vCount) += myVertexTrack.pt()/sumTrackPt;
1006 associatedTrackCount.at(vCount) += 1/nJetTracks;
1007 }
1008 }
1009 }
1010 }
1011 ++vCount;
1012 }
1013
1014 float maxSumPtFraction = 0; float maxCountFraction = 0;
1015 int vtxSumPtIndex = 0; int vtxCountIndex = 0;
1016 int count = 0;
1017
1018 for (int i = 0; i < vCount; ++i) {
1019 if (associatedTrackSumPt.at(i) > maxSumPtFraction) {
1020 maxSumPtFraction = associatedTrackSumPt.at(i);
1021 vtxSumPtIndex = count + 1;
1022 }
1023 if (associatedTrackCount.at(i) > maxCountFraction) {
1024 maxCountFraction = associatedTrackCount.at(i);
1025 vtxCountIndex = count + 1;
1026 }
1027 ++count;
1028 }
1029 outJet->SetVtxSumPtFrac(maxSumPtFraction);
1030 outJet->SetVtxSumPt(sumTrackPt);
1031 outJet->SetVtxTrackFrac(maxCountFraction);
1032 outJet->SetVtxNTracks(nJetTracks);
1033 outJet->SetVtxSumPtIndex(vtxSumPtIndex);
1034 outJet->SetVtxCountIndex(vtxCountIndex);
1035 }
1036
1037 return true;
1038 }
1039
1040 //void ntupleProducer::JetSelector(Handle<reco::PFJetsCollection> inJets, TClonesArray& outJets)
1041 //{
1042 //}
1043
1044 bool ntupleProducer::electronMVA(Handle<reco::VertexCollection> vtxCollection, vector<pat::Electron>::const_iterator iElectron)
1045 {
1046 if (vtxCollection->size() != 0) return false;
1047
1048 /*
1049 edm::ESHandle<TransientTrackBuilder> builder;
1050 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
1051 TransientTrackBuilder thebuilder = *(builder.product());
1052
1053 InputTag reducedEBRecHitCollection(string("reducedEcalRecHitsEB"));
1054 InputTag reducedEERecHitCollection(string("reducedEcalRecHitsEE"));
1055
1056 EcalClusterLazyTools lazyTools(iEvent, iSetup, reducedEBRecHitCollection, reducedEERecHitCollection);
1057
1058 pv = &*vtxCollection->begin();
1059 double myMVANonTrigMethod1 = myMVANonTrig->mvaValue(*iElectron,*pv,thebuilder,lazyTools,false);
1060 double myMVATrigMethod1 = myMVATrig->mvaValue(*iElectron,*pv,thebuilder,lazyTools,false);
1061
1062 //eleCon->SetIdMap("MVATrigMethod1",myMVATrigMethod1);
1063 //eleCon->SetIdMap("MVANonTrigMethod1",myMVANonTrigMethod1);
1064
1065 //ID'd electrons to feed into MVAIso
1066
1067 InputTag gsfEleLabel(string("gsfElectrons"));
1068 Handle<reco::GsfElectronCollection> theEGammaCollection;
1069 iEvent.getByLabel(gsfEleLabel,theEGammaCollection);
1070
1071 reco::GsfElectronCollection identifiedElectrons;
1072 for (reco::GsfElectronCollection::const_iterator iE = theEGammaCollection->begin(); iE != theEGammaCollection->end(); ++iE) {
1073
1074 double electronTrackZ = 0;
1075 if (iE->gsfTrack().isNonnull()) {
1076 electronTrackZ = iE->gsfTrack()->dz(pv->position());
1077 } else if (iE->closestCtfTrackRef().isNonnull()) {
1078 electronTrackZ = iE->closestCtfTrackRef()->dz(pv->position());
1079 }
1080 if(fabs(electronTrackZ) > 0.2) continue;
1081
1082
1083 if(fabs(iE->superCluster()->eta())<1.479) {
1084 if(iE->pt() > 20) {
1085 if(
1086 iE->sigmaIetaIeta() > 0.01
1087 || fabs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.007
1088 || fabs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8
1089 || iE->hadronicOverEm() > 0.15
1090 ) continue;
1091 } else {
1092 if(iE->sigmaIetaIeta() > 0.012) continue;
1093 if(fabs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.007) continue;
1094 if(fabs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8) continue;
1095 if(iE->hadronicOverEm() > 0.15) continue;
1096 }
1097 } else {
1098 if(iE->pt() > 20) {
1099 if(iE->sigmaIetaIeta() > 0.03) continue;
1100 if(fabs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.010) continue;
1101 if(fabs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8) continue;
1102 } else {
1103 if(iE->sigmaIetaIeta() > 0.032) continue;
1104 if(fabs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.010) continue;
1105 if(fabs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8) continue;
1106 }
1107 }
1108 identifiedElectrons.push_back(*iE);
1109 }
1110
1111 Handle<reco::MuonCollection> hMuonProduct;
1112 iEvent.getByLabel("muons", hMuonProduct);
1113 const reco::MuonCollection inMuons = *(hMuonProduct.product());
1114
1115 reco::MuonCollection identifiedMuons;
1116 for (reco::MuonCollection::const_iterator iMuon = inMuons.begin(); iMuon != inMuons.end(); ++iMuon) {
1117 if (
1118 iMuon->innerTrack().isNonnull()
1119 && iMuon->isGlobalMuon()
1120 && iMuon->isTrackerMuon();
1121 && iMuon->innerTrack()->numberOfValidHits() > 11
1122 ) identifiedMuons.push_back(*iMuon);
1123 }
1124
1125 Handle<PFCandidateCollection> pfCands;
1126 iEvent.getByLabel("particleFlow", pfCands);
1127 const PFCandidateCollection pfCanIso = *(pfCands.product());
1128 double isomva = fElectronIsoMVA->mvaValue( *iElectron, *pv, pfCanIso, rhoFactor, ElectronEffectiveArea::kEleEAData2011, identifiedElectrons, identifiedMuons);
1129 */
1130
1131 return true;
1132 }
1133
1134 //define this as a plug-in
1135 DEFINE_FWK_MODULE(ntupleProducer);