ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/devildog/SWonAnalysis3/Thesis/src/Thesis.cc
Revision: 1.4
Committed: Thu Apr 28 21:55:41 2011 UTC (14 years ago) by devildog
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +201 -129 lines
Log Message:
Changes for 41X, fix b-tag bug, fix jet vertex association

File Contents

# User Rev Content
1 devildog 1.1 // system include files
2     #include <memory>
3    
4     // user include files
5     #include "FWCore/Framework/interface/Frameworkfwd.h"
6     #include "FWCore/Framework/interface/EDAnalyzer.h"
7     #include "FWCore/Framework/interface/Event.h"
8     #include "FWCore/Framework/interface/MakerMacros.h"
9     #include "FWCore/ParameterSet/interface/ParameterSet.h"
10     #include "FWCore/Framework/interface/LuminosityBlock.h"
11     #include "DataFormats/BeamSpot/interface/BeamSpot.h"
12     #include "FWCore/Common/interface/TriggerNames.h"
13     #include "DataFormats/Common/interface/TriggerResults.h"
14     #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
15    
16     #include "DataFormats/Common/interface/View.h"
17     #include "DataFormats/TrackReco/interface/TrackFwd.h"
18     #include "DataFormats/TrackReco/interface/Track.h"
19     #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
20     #include "DataFormats/MuonReco/interface/Muon.h"
21     #include "DataFormats/MuonReco/interface/MuonFwd.h"
22     #include "DataFormats/JetReco/interface/Jet.h"
23     #include "DataFormats/EgammaCandidates/interface/Electron.h"
24     #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
25     #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
26     #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
27     #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
28    
29     #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
30     #include "Geometry/Records/interface/CaloGeometryRecord.h"
31     #include "DataFormats/Math/interface/deltaPhi.h"
32     #include "DataFormats/GeometryVector/interface/VectorUtil.h"
33    
34     #include "DataFormats/JetReco/interface/CaloJet.h"
35     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
36     #include "DataFormats/JetReco/interface/PFJet.h"
37     #include "DataFormats/JetReco/interface/PFJetCollection.h"
38     #include "DataFormats/METReco/interface/PFMETFwd.h"
39     #include "DataFormats/METReco/interface/PFMET.h"
40     #include "DataFormats/METReco/interface/PFMETCollection.h"
41     #include "DataFormats/TauReco/interface/PFTau.h"
42     #include "DataFormats/TauReco/interface/PFTauFwd.h"
43     #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
44     #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
45     #include "DataFormats/BeamSpot/interface/BeamSpot.h"
46     #include "DataFormats/JetReco/interface/GenJet.h"
47     #include "DataFormats/JetReco/interface/GenJetCollection.h"
48     #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
49     #include "DataFormats/JetReco/interface/Jet.h"
50     #include "DataFormats/VertexReco/interface/Vertex.h"
51     #include "DataFormats/VertexReco/interface/VertexFwd.h"
52     #include "DataFormats/BTauReco/interface/JetTag.h"
53     #include "JetMETCorrections/Objects/interface/JetCorrector.h"
54     #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
55     #include "DataFormats/Common/interface/ValueMap.h"
56     #include "DataFormats/Math/interface/deltaPhi.h"
57    
58 devildog 1.4 #include "DataFormats/PatCandidates/interface/Electron.h"
59     #include "DataFormats/PatCandidates/interface/Muon.h"
60     #include "DataFormats/PatCandidates/interface/MET.h"
61     #include "DataFormats/PatCandidates/interface/Jet.h"
62    
63 devildog 1.1 #include "TH1.h"
64     #include "TH2.h"
65     #include "TFile.h"
66     #include "TTree.h"
67     #include "TCJet.h"
68     #include "TCPrimaryVtx.h"
69     #include "TCMET.h"
70     #include "TCTrig.h"
71     #include "TCElectron.h"
72     #include "TCMuon.h"
73     #include "TCTau.h"
74     #include "TClonesArray.h"
75 devildog 1.2 #include "TGenParticle.h"
76     #include "TGenMET.h"
77 devildog 1.1 #include "TLorentzVector.h"
78 devildog 1.2 #include "TVector3.h"
79 devildog 1.1 #include <string>
80    
81     using namespace edm;
82     using namespace std;
83     using namespace reco;
84    
85     class Thesis : public edm::EDAnalyzer {
86     public:
87     explicit Thesis(const edm::ParameterSet&);
88     ~Thesis();
89    
90     private:
91     virtual void beginJob() ;
92     virtual void analyze(const edm::Event&, const edm::EventSetup&);
93     virtual void endJob() ;
94     virtual void beginRun(const edm::Run&, const edm::EventSetup&) ;
95     double sumPtSquared(const Vertex& v);
96     bool triggerDecision(edm::Handle<edm::TriggerResults>& hltR, int iTrigger);
97    
98     int eventNumber, runNumber, lumiSection, bunchCross;
99     TCMET MET;
100     TCTrig triggers;
101     TTree * sTree;
102     TFile * ntupleFile;
103     TClonesArray * jet_ak5PF;
104     TClonesArray * primaryVtx;
105     TClonesArray * electrons;
106     TClonesArray * muons;
107     TClonesArray * taus;
108 devildog 1.2 TClonesArray * genParticles;
109     TVector3 * beamSpot;
110 devildog 1.4
111 devildog 1.1 edm::InputTag PFJetHandle_;
112     edm::InputTag PrimaryVtxHandle_;
113     edm::InputTag muonTag_;
114     edm::InputTag electronTag_;
115     edm::InputTag electronIDMap95_;
116     edm::InputTag electronIDMap90_;
117     edm::InputTag electronIDMap85_;
118     edm::InputTag electronIDMap80_;
119     edm::InputTag electronIDMap70_;
120     edm::InputTag electronIDMap60_;
121    
122     edm::InputTag triggerResultsTag_;
123     edm::InputTag tauJetSource_;
124    
125     edm::InputTag tauDiscrByLeadTrackFinding_;
126     edm::InputTag tauDiscrByLeadTrackPtCut_;
127     edm::InputTag tauDiscrByTrackIso_;
128     edm::InputTag tauDiscrByEcalIso_;
129     edm::InputTag tauDiscrAgainstElectrons_;
130     edm::InputTag tauDiscrAgainstMuons_;
131     edm::InputTag tauDiscrTaNC_;
132    
133     std::string rootfilename;
134     bool doMCMatching;
135     bool keepAllEvents;
136     bool triggerHLT_;
137 devildog 1.2 int isMC;
138    
139 devildog 1.1 string hlTriggerResults_, hltProcess_, triggerName_;
140     TriggerNames triggerNames;
141     HLTConfigProvider hltConfig_;
142     std::vector<string> hlNames;
143     };
144    
145     Thesis::Thesis(const edm::ParameterSet& iConfig):
146     PFJetHandle_ (iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag")),
147     PrimaryVtxHandle_ (iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")),
148     muonTag_ (iConfig.getParameter<edm::InputTag>("muonTag")),
149     electronTag_ (iConfig.getParameter<edm::InputTag>("electronTag")),
150     electronIDMap95_ (iConfig.getParameter<edm::InputTag>("electronIDMap95")),
151     electronIDMap90_ (iConfig.getParameter<edm::InputTag>("electronIDMap90")),
152     electronIDMap85_ (iConfig.getParameter<edm::InputTag>("electronIDMap85")),
153     electronIDMap80_ (iConfig.getParameter<edm::InputTag>("electronIDMap80")),
154     electronIDMap70_ (iConfig.getParameter<edm::InputTag>("electronIDMap70")),
155     electronIDMap60_ (iConfig.getParameter<edm::InputTag>("electronIDMap60")),
156     rootfilename (iConfig.getUntrackedParameter<std::string>("rootfilename")),
157     doMCMatching (iConfig.getParameter<bool>("doMCMatch")),
158     keepAllEvents (iConfig.getParameter<bool>("keepAllEvents")),
159     triggerHLT_ (iConfig.getUntrackedParameter<bool>("triggerHLT")),
160     hlTriggerResults_ (iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")),
161     hltProcess_ (iConfig.getUntrackedParameter<string>("hltName"))
162     { }
163    
164     Thesis::~Thesis()
165     { }
166    
167     // ------------ method called to for each event ------------
168     void
169     Thesis::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
170     {
171     using namespace edm;
172     using namespace reco;
173    
174     eventNumber = iEvent.id().event();
175     runNumber = iEvent.id().run();
176     lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
177     bunchCross = (unsigned int)iEvent.bunchCrossing();
178 devildog 1.2 isMC = 0;
179    
180     if(doMCMatching){
181     int MCcount = 0;
182     isMC = 1;
183     edm::Handle<reco::GenParticleCollection> pGenPart;
184     iEvent.getByLabel("genParticles", pGenPart);
185     reco::GenParticleCollection::const_iterator pGen_;
186     for(pGen_ = pGenPart->begin(); pGen_ != pGenPart->end(); pGen_++) {
187     if( (abs(pGen_->pdgId()) == 11) || //electron
188     (abs(pGen_->pdgId()) == 12) || //e v
189     (abs(pGen_->pdgId()) == 13) || //mu
190     (abs(pGen_->pdgId()) == 14) || //mu v
191     (abs(pGen_->pdgId()) == 15) || //tau
192     (abs(pGen_->pdgId()) == 16) || //tau v
193     (abs(pGen_->pdgId()) == 5) || //b
194     (abs(pGen_->pdgId()) == 6) || //t
195     (abs(pGen_->pdgId()) == 37) || //H+/-
196     (abs(pGen_->pdgId()) == 24) //W+/-
197     ) {
198     TGenParticle * mcCon = new ((*genParticles)[MCcount]) TGenParticle;
199     mcCon->SetPDGId(pGen_->pdgId());
200     mcCon->SetP4(pGen_->px(), pGen_->py(), pGen_->pz(),pGen_->energy());
201     mcCon->SetCharge(pGen_->charge());
202     if(pGen_->numberOfMothers() != 0) mcCon->SetMother(pGen_->mother(0)->pdgId());
203     MCcount++;
204     }
205     }
206     }
207 devildog 1.1
208     edm::Handle<reco::VertexCollection> primaryVtcs;
209     iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
210    
211     int vtxCount = 0;
212     for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
213     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
214     if(!myVtx.isValid() || myVtx.isFake()) continue;
215     TCPrimaryVtx* vtxCon = new ((*primaryVtx)[vtxCount]) TCPrimaryVtx;
216     vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
217     vtxCon->SetNDof(myVtx.ndof());
218     vtxCon->SetChi2(myVtx.chi2());
219     vtxCon->SetNTrks(myVtx.tracksSize());
220     vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
221     // if(vtxCount == 0) primaryVertexZ = myVtx.z(); //Shouldn't be needed, using VertexFilter upstream ofhe analyzer
222     ++vtxCount;
223     }
224    
225     edm::Handle<reco::BeamSpot> beamSpotHandle;
226     iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
227     reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
228    
229     edm::Handle<reco::PFCandidateCollection> PFCandCollection;
230     iEvent.getByLabel("particleFlow", PFCandCollection);
231    
232     int muoncount = 0;
233     int passmuons = 0;
234    
235     //Muons
236 devildog 1.2 beamSpot->SetX(0);
237     beamSpot->SetY(0);
238     beamSpot->SetZ(0);
239    
240 devildog 1.1 math::XYZPoint BSPoint(vertexBeamSpot.x0(),vertexBeamSpot.y0(),vertexBeamSpot.z0());
241     edm::Handle<MuonCollection> muColl;
242     iEvent.getByLabel(muonTag_,muColl);
243 devildog 1.2
244     beamSpot->SetX(vertexBeamSpot.x0());
245     beamSpot->SetY(vertexBeamSpot.y0());
246     beamSpot->SetZ(vertexBeamSpot.z0());
247 devildog 1.4
248     // cout << "Muons... ";
249    
250     for(reco::MuonCollection::const_iterator mu_ = muColl->begin(); mu_ != muColl->end(); mu_++){
251 devildog 1.1 reco::Muon myMuon = reco::Muon(*mu_);
252     if(!(myMuon.isGlobalMuon() && myMuon.isTrackerMuon())) continue;
253     TCMuon * muCon = new ((*muons)[muoncount]) TCMuon;
254     muCon->SetP4(myMuon.px(), myMuon.py(), myMuon.pz(), myMuon.energy());
255     // muCon->SetP4(myMuon.track().px(), myMuon.track().py(), myMuon.track().pz(), myMuon.track().energy());
256     muCon->SetPosition(mu_->globalTrack()->vx(),mu_->globalTrack()->vy(),mu_->globalTrack()->vz());
257     muCon->SetCharge(myMuon.charge());
258     muCon->SetsumPt03(myMuon.isolationR03().sumPt);
259     muCon->SetemEt03(myMuon.isolationR03().emEt);
260     muCon->SethadEt03(myMuon.isolationR03().hadEt);
261     muCon->SetnTracks03(myMuon.isolationR03().nTracks);
262     muCon->SetnumberOfValidPixelHits(myMuon.globalTrack()->hitPattern().numberOfValidPixelHits());
263     muCon->SetnumberOfValidTrackerHits(myMuon.globalTrack()->hitPattern().numberOfValidTrackerHits());
264     muCon->SetnumberOfValidMuonHits(myMuon.globalTrack()->hitPattern().numberOfValidMuonHits());
265     muCon->SetnumberOfMissingPixelHits(myMuon.globalTrack()->hitPattern().numberOfLostPixelHits());
266     muCon->SetnumberOfMissingTrackerHits(myMuon.globalTrack()->hitPattern().numberOfLostTrackerHits());
267     muCon->SetnormalizedChi2(myMuon.globalTrack()->normalizedChi2());
268     float sumPt3 = 0;
269     float gamma3 = 0;
270     float neutral3 = 0;
271     float sumPt4 = 0;
272     float gamma4 = 0;
273     float neutral4 = 0;
274     float sumPt5 = 0;
275     float gamma5 = 0;
276     float neutral5 = 0;
277     for(reco::PFCandidateCollection::const_iterator PF_ = PFCandCollection->begin(); PF_ != PFCandCollection->end(); ++PF_) {
278     reco::PFCandidate pf = reco::PFCandidate(*PF_);
279     double deltaR = Geom::deltaR(myMuon.momentum(),pf.momentum());
280     if (deltaR > 0.5) continue;
281     if (pf.particleId()==reco::PFCandidate::h){// || pf.particleId()==reco::PFCandidate::e || pf.particleId()==reco::PFCandidate::mu ) {
282 devildog 1.3 if(deltaR > 0.01) {
283     if(!pf.trackRef().isNull()) {
284     if((pf.pt() > .5) && (fabs(pf.trackRef()->dz(BSPoint) - myMuon.innerTrack()->dz(BSPoint)) < 0.2 ) ){
285     sumPt5 += pf.pt();
286     if (deltaR < 0.4) sumPt4 += pf.pt();
287     if (deltaR < 0.3) sumPt3 += pf.pt();
288     }
289 devildog 1.1 }
290     }
291     }else if(pf.particleId()==reco::PFCandidate::e){
292     if(deltaR > 0.01) {
293 devildog 1.3 if(!pf.gsfTrackRef().isNull()) {
294     if((pf.pt() > .5) && (fabs(pf.gsfTrackRef()->dz(BSPoint) - myMuon.innerTrack()->dz(BSPoint)) < 0.2 ) ){
295     sumPt5 += pf.pt();
296     if (deltaR < 0.4) sumPt4 += pf.pt();
297     if (deltaR < 0.3) sumPt3 += pf.pt();
298     }
299 devildog 1.1 }
300     }
301     }else if(pf.particleId()==reco::PFCandidate::mu){
302     if(deltaR > 0.01) {
303 devildog 1.3 if(!pf.muonRef().isNull()) {
304     if((pf.pt() > .5) && (fabs(pf.muonRef()->innerTrack()->dz(BSPoint) - myMuon.innerTrack()->dz(BSPoint)) < 0.2 ) ){
305     sumPt5 += pf.pt();
306     if (deltaR < 0.4) sumPt4 += pf.pt();
307     if (deltaR < 0.3) sumPt3 += pf.pt();
308     }
309 devildog 1.1 }
310     }
311     }else if(pf.particleId()==reco::PFCandidate::h0 || pf.particleId()==reco::PFCandidate::h_HF) {
312     if(deltaR > 0.08) {
313     if(pf.pt() > 1) {
314     neutral5 += pf.pt();
315     if (deltaR < 0.4) neutral4 += pf.pt();
316     if (deltaR < 0.3) neutral3 += pf.pt();
317     }
318     }
319     }else if(pf.particleId()==reco::PFCandidate::gamma || pf.particleId()==reco::PFCandidate::egamma_HF) {
320     if(deltaR > 0.05) {
321     if(pf.pt() > 1) {
322     gamma5 += pf.pt();
323     if (deltaR < 0.4) gamma4 += pf.pt();
324     if (deltaR < 0.3) gamma3 += pf.pt();
325     }
326     }
327     }
328     }
329     muCon->SetPFSumPt(0.5, sumPt5);
330     muCon->SetPFSumPt(0.4, sumPt4);
331     muCon->SetPFSumPt(0.3, sumPt3);
332     muCon->SetPFEGamma(0.5, gamma5);
333     muCon->SetPFEGamma(0.4, gamma4);
334     muCon->SetPFEGamma(0.3, gamma3);
335     muCon->SetPFENeutral(0.5, neutral5);
336     muCon->SetPFENeutral(0.4, neutral4);
337     muCon->SetPFENeutral(0.3, neutral3);
338    
339     muCon->SetDxy(myMuon.innerTrack()->dxy(BSPoint));
340     muCon->SetDz(myMuon.innerTrack()->dz(BSPoint));
341    
342     if(
343     (myMuon.innerTrack()->hitPattern().pixelLayersWithMeasurement() > 0 ) &&
344     (myMuon.globalTrack()->hitPattern().numberOfValidTrackerHits() > 10) &&
345     (myMuon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0) &&
346     (myMuon.globalTrack()->normalizedChi2() < 10) &&
347     (myMuon.numberOfMatches() > 1) &&
348     // (fabs(myMuon.innerTrack()->dxy(BSPoint)) < 0.1) && //WAS 0.02
349     // (fabs(myMuon.innerTrack()->dz(BSPoint)) < 5) && //WAS 1
350 devildog 1.2 (myMuon.pt() > 12) &&
351 devildog 1.1 (fabs(myMuon.eta()) < 2.5) && //WAS 2.1
352     ((myMuon.isolationR03().emEt + myMuon.isolationR03().hadEt + myMuon.isolationR03().sumPt) / myMuon.pt() < 0.5 )
353     ) passmuons++;
354     muoncount++;
355     }
356    
357     int passelectrons = 0;
358     //Electrons
359 devildog 1.4 // cout << "Electrons... ";
360 devildog 1.1 unsigned int electroncount = 0;
361     edm::Handle<GsfElectronCollection> eColl;
362     iEvent.getByLabel(electronTag_,eColl);
363    
364     edm::Handle<edm::ValueMap<float> > eIDValueMap95;
365     iEvent.getByLabel( electronIDMap95_ , eIDValueMap95 );
366     const edm::ValueMap<float> & eIDmap95 = * eIDValueMap95 ;
367    
368     edm::Handle<edm::ValueMap<float> > eIDValueMap90;
369     iEvent.getByLabel( electronIDMap90_ , eIDValueMap90 );
370     const edm::ValueMap<float> & eIDmap90 = * eIDValueMap90 ;
371    
372     edm::Handle<edm::ValueMap<float> > eIDValueMap85;
373     iEvent.getByLabel( electronIDMap85_ , eIDValueMap85 );
374     const edm::ValueMap<float> & eIDmap85 = * eIDValueMap85 ;
375    
376     edm::Handle<edm::ValueMap<float> > eIDValueMap80;
377     iEvent.getByLabel( electronIDMap80_ , eIDValueMap80 );
378     const edm::ValueMap<float> & eIDmap80 = * eIDValueMap80 ;
379    
380     edm::Handle<edm::ValueMap<float> > eIDValueMap70;
381     iEvent.getByLabel( electronIDMap70_ , eIDValueMap70 );
382     const edm::ValueMap<float> & eIDmap70 = * eIDValueMap70 ;
383    
384     edm::Handle<edm::ValueMap<float> > eIDValueMap60;
385     iEvent.getByLabel( electronIDMap60_ , eIDValueMap60 );
386     const edm::ValueMap<float> & eIDmap60 = * eIDValueMap60 ;
387    
388    
389     for(reco::GsfElectronCollection::const_iterator e_ = eColl->begin(); e_ != eColl->end(); e_++)
390     {
391     edm::Ref<reco::GsfElectronCollection> electronRef(eColl,electroncount);
392     int cuts95 = eIDmap95[electronRef];
393     int cuts90 = eIDmap90[electronRef];
394     int cuts85 = eIDmap85[electronRef];
395     int cuts80 = eIDmap80[electronRef];
396     int cuts70 = eIDmap70[electronRef];
397     int cuts60 = eIDmap60[electronRef];
398     reco::GsfElectron myElectron = reco::GsfElectron(*e_);
399     TCElectron * eCon = new ((*electrons)[electroncount]) TCElectron;
400     eCon->SetP4(myElectron.px(), myElectron.py(), myElectron.pz(), myElectron.energy());
401     eCon->SetPosition(electronRef->gsfTrack()->vx(),electronRef->gsfTrack()->vy(),electronRef->gsfTrack()->vz());
402    
403     eCon->SetnumberOfValidPixelHits(myElectron.gsfTrack()->hitPattern().numberOfValidPixelHits());
404     eCon->SetnumberOfValidTrackerHits(myElectron.gsfTrack()->hitPattern().numberOfValidTrackerHits());
405     eCon->SetnumberOfMissingPixelHits(myElectron.gsfTrack()->hitPattern().numberOfLostPixelHits());
406     eCon->SetnumberOfMissingTrackerHits(myElectron.gsfTrack()->hitPattern().numberOfLostTrackerHits());
407    
408     eCon->SetisEB(myElectron.isEB());
409     eCon->SetisEE(myElectron.isEE());
410     eCon->SetisInGap(myElectron.isGap());
411    
412     eCon->SetCharge(myElectron.charge());
413     eCon->SetsumPt03(myElectron.dr03TkSumPt());
414     eCon->SetemEt03(myElectron.dr03EcalRecHitSumEt());
415     eCon->SethadEt03(myElectron.dr03HcalTowerSumEt());
416     eCon->SetCutLevel(cuts95, 95);
417     eCon->SetCutLevel(cuts90, 90);
418     eCon->SetCutLevel(cuts85, 85);
419     eCon->SetCutLevel(cuts80, 80);
420     eCon->SetCutLevel(cuts70, 70);
421     eCon->SetCutLevel(cuts60, 60);
422     float sumPt3 = 0;
423     float gamma3 = 0;
424     float neutral3 = 0;
425     float sumPt4 = 0;
426     float gamma4 = 0;
427     float neutral4 = 0;
428     float sumPt5 = 0;
429     float gamma5 = 0;
430     float neutral5 = 0;
431     for(reco::PFCandidateCollection::const_iterator PF_ = PFCandCollection->begin(); PF_ != PFCandCollection->end(); ++PF_) {
432     reco::PFCandidate pf = reco::PFCandidate(*PF_);
433     double deltaR = Geom::deltaR(myElectron.momentum(),pf.momentum());
434     if (deltaR > 0.5) continue;
435    
436     if (pf.particleId()==reco::PFCandidate::h){// || pf.particleId()==reco::PFCandidate::e || pf.particleId()==reco::PFCandidate::mu ) {
437     if(deltaR > 0.01) {
438 devildog 1.3 if(!pf.trackRef().isNull()) {
439     if((pf.pt() > .5) && (fabs(pf.trackRef()->dz(BSPoint) - myElectron.gsfTrack()->dz(BSPoint)) < 0.2 ) ){
440     sumPt5 += pf.pt();
441     if (deltaR < 0.4) sumPt4 += pf.pt();
442     if (deltaR < 0.3) sumPt3 += pf.pt();
443     }
444 devildog 1.1 }
445     }
446     }else if(pf.particleId()==reco::PFCandidate::e){
447     if(deltaR > 0.01) {
448 devildog 1.3 if(!pf.gsfTrackRef().isNull()) {
449     if((pf.pt() > .5) && (fabs(pf.gsfTrackRef()->dz(BSPoint) - myElectron.gsfTrack()->dz(BSPoint)) < 0.2 ) ){
450     sumPt5 += pf.pt();
451     if (deltaR < 0.4) sumPt4 += pf.pt();
452     if (deltaR < 0.3) sumPt3 += pf.pt();
453     }
454 devildog 1.1 }
455     }
456     }else if(pf.particleId()==reco::PFCandidate::mu){
457     if(deltaR > 0.01) {
458 devildog 1.3 if(!pf.muonRef().isNull()) {
459     if((pf.pt() > .5) && (fabs(pf.muonRef()->innerTrack()->dz(BSPoint) - myElectron.gsfTrack()->dz(BSPoint)) < 0.2 ) ){
460     sumPt5 += pf.pt();
461     if (deltaR < 0.4) sumPt4 += pf.pt();
462     if (deltaR < 0.3) sumPt3 += pf.pt();
463     }
464 devildog 1.1 }
465     }
466     }else if(pf.particleId()==reco::PFCandidate::h0 || pf.particleId()==reco::PFCandidate::h_HF) {
467     if(deltaR > 0.08) {
468     if(pf.pt() > 1) {
469     neutral5 += pf.pt();
470     if (deltaR < 0.4) neutral4 += pf.pt();
471     if (deltaR < 0.3) neutral3 += pf.pt();
472     }
473     }
474     }else if(pf.particleId()==reco::PFCandidate::gamma || pf.particleId()==reco::PFCandidate::egamma_HF) {
475     if(deltaR > 0.05) {
476     if(pf.pt() > 1) {
477     gamma5 += pf.pt();
478     if (deltaR < 0.4) gamma4 += pf.pt();
479     if (deltaR < 0.3) gamma3 += pf.pt();
480    
481     }
482     }
483     }
484     }
485     eCon->SetPFSumPt(0.5, sumPt5);
486     eCon->SetPFSumPt(0.4, sumPt4);
487     eCon->SetPFSumPt(0.3, sumPt3);
488     eCon->SetPFEGamma(0.5, gamma5);
489     eCon->SetPFEGamma(0.4, gamma4);
490     eCon->SetPFEGamma(0.3, gamma3);
491     eCon->SetPFENeutral(0.5, neutral5);
492     eCon->SetPFENeutral(0.4, neutral4);
493     eCon->SetPFENeutral(0.3, neutral3);
494 devildog 1.2
495     eCon->SetDxy(myElectron.gsfTrack()->dxy(BSPoint));
496     eCon->SetDz(myElectron.gsfTrack()->dz(BSPoint));
497    
498 devildog 1.1 if((myElectron.et() > 15) && (fabs(myElectron.eta() < 2.5)) && (cuts95 == 7)) passelectrons++;
499     //use this wider cut for making ntuples
500    
501     electroncount++;
502     }
503    
504     //Jets
505 devildog 1.4 // cout << "Jets... ";
506    
507 devildog 1.1 edm::Handle<reco::JetTagCollection> bTagHandle1;
508     iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
509     const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
510 devildog 1.4
511     edm::Handle<reco::JetTagCollection> bTagHandle2;
512     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
513     const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
514    
515     edm::Handle<reco::JetTagCollection> bTagHandle3;
516     iEvent.getByLabel("combinedSecondaryVertexMVABJetTags", bTagHandle3);
517     const reco::JetTagCollection & bTags3 = *(bTagHandle3.product());
518    
519     reco::JetTagCollection::const_iterator jet_it_1, jet_it_2, jet_it_3;
520 devildog 1.1
521 devildog 1.4
522    
523     //const JetCorrector* correctorL1 = JetCorrector::getJetCorrector("ak5PFL1Offset",iSetup);
524     const JetCorrector* correctorL2 = JetCorrector::getJetCorrector("ak5PFL2Relative",iSetup);
525     const JetCorrector* correctorL3 = JetCorrector::getJetCorrector("ak5PFL3Absolute",iSetup);
526     // const JetCorrector* correctorRes = JetCorrector::getJetCorrector("ak5PFResidual", iSetup);
527 devildog 1.1
528     Handle<reco::PFJetCollection> PFJets;
529     iEvent.getByLabel(PFJetHandle_, PFJets);
530    
531     int PFcount = 0;
532    
533     for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
534    
535     reco::PFJet myJet = reco::PFJet(*jet_iter);
536 devildog 1.4 reco::PFJet corJet = reco::PFJet(*jet_iter);
537    
538     float scale2 = correctorL2->correction(corJet);
539     corJet.scaleEnergy(scale2);
540     float scale3 = correctorL3->correction(corJet);
541     corJet.scaleEnergy(scale3);
542    
543     if (corJet.pt() > 5) {
544     TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
545    
546     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
547     jetCon->SetVtx(-999.0, -999.0, -999.0);
548     jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
549     jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
550     jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
551     jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
552     jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
553     jetCon->SetNumChPart(myJet.chargedMultiplicity());
554    
555    
556     for (jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
557     if (sqrt(pow(jet_it_1->first->eta() - corJet.eta(), 2) + pow(deltaPhi(jet_it_1->first->phi(),corJet.phi()), 2)) == 0.) {
558     jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
559     break; }
560     }
561    
562     for (jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
563     if (sqrt(pow(jet_it_2->first->eta() - corJet.eta(), 2) + pow(deltaPhi(jet_it_2->first->phi(),corJet.phi()), 2)) == 0.) {
564     jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
565     break; }
566     }
567 devildog 1.1
568 devildog 1.4 for (jet_it_3 = bTags3.begin(); jet_it_3 != bTags3.end(); jet_it_3++) {
569     if (sqrt(pow(jet_it_3->first->eta() - corJet.eta(), 2) + pow(deltaPhi(jet_it_3->first->phi(),corJet.phi()), 2)) == 0.) {
570     jetCon->SetBDiscrSecVtxMVA(jet_it_3->second);
571     break; }
572 devildog 1.1 }
573    
574    
575 devildog 1.4 jetCon->SetJetCorr(2, scale2);
576     jetCon->SetJetCorr(3, scale3);
577    
578     /////////////////////////
579     //get associated tracks//
580     /////////////////////////
581    
582     const reco::TrackRefVector &tracks = myJet.getTrackRefs();
583    
584     vector<TVector3> vtxPositionCollection;
585     vector<float> associatedTrackSumPt;
586     vector<const reco::Track*> jetTrackAddresses;
587     // vector<unsigned int> jetTrackAddresses;
588     float sumTrackX, sumTrackY, sumTrackZ, sumTrackIP, sumTrackPt;
589     int nJetTracks, nVertexTracks, nAssociatedTracks, vertexIndex;
590     int vCount = 0;
591    
592     nJetTracks = nVertexTracks = nAssociatedTracks = 0;
593     sumTrackX = sumTrackY = sumTrackZ = sumTrackIP = sumTrackPt = 0;
594    
595     if(fabs(myJet.eta()) < 2.5){
596    
597     for (TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack) {
598     const reco::Track &myJetTrack = **iTrack;
599    
600     sumTrackPt += myJetTrack.pt();
601     sumTrackX += myJetTrack.vx();
602     sumTrackY += myJetTrack.vy();
603     sumTrackZ += myJetTrack.vz();
604     sumTrackIP += myJetTrack.dxy(vertexBeamSpot.position());
605     jetTrackAddresses.push_back(&myJetTrack);
606     // jetTrackAddresses.push_back((unsigned int)&myJetTrack);
607     ++nJetTracks;
608     }
609    
610     if(jetTrackAddresses.size() > 0){
611    
612     for (VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter) {
613     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
614     if(!myVtx.isValid() || myVtx.isFake()) continue;
615     TVector3 *iVtxPosition = new TVector3(myVtx.x(), myVtx.y(), myVtx.z());
616     vtxPositionCollection.push_back(*iVtxPosition);
617     associatedTrackSumPt.push_back(0);
618    
619     for(Vertex::trackRef_iterator iTrackRef = myVtx.tracks_begin(); iTrackRef != myVtx.tracks_end(); ++iTrackRef){
620     const edm::RefToBase<reco::Track> &myTrackRef = *iTrackRef;
621     if(myTrackRef.isAvailable()){
622     const reco::Track &myVertexTrack = *myTrackRef.get();
623    
624     for(vector<const reco::Track*>::const_iterator iTrackAddress = jetTrackAddresses.begin();
625     iTrackAddress != jetTrackAddresses.end(); ++iTrackAddress){
626     if (*iTrackAddress == &myVertexTrack) {
627     associatedTrackSumPt.at(vCount) += myVertexTrack.pt()/sumTrackPt;
628     ++nAssociatedTracks;
629     }
630     }
631     }
632     }
633     ++vCount;
634     }
635    
636     float maxSumPtFraction = 0;
637     vCount = vertexIndex = 0;
638    
639     for (vector<float>::const_iterator iTrackSumPt = associatedTrackSumPt.begin(); iTrackSumPt != associatedTrackSumPt.end(); ++iTrackSumPt) {
640     if (*iTrackSumPt > maxSumPtFraction) {
641     maxSumPtFraction = *iTrackSumPt;
642     vertexIndex = vCount + 1;
643     }
644     ++vCount;
645 devildog 1.1 }
646 devildog 1.4
647     jetCon->SetVtxSumPtFrac(maxSumPtFraction);
648     jetCon->SetVtxSumPt(sumTrackPt);
649     jetCon->SetVtxTrackFrac((float)nAssociatedTracks/(float)nJetTracks);
650     jetCon->SetVtxNTracks(nJetTracks);
651     jetCon->SetVtxIndex(vertexIndex);
652     }
653 devildog 1.1 }
654 devildog 1.4 ++PFcount;
655 devildog 1.1 }
656     }
657    
658     //MET
659     edm::Handle<reco::PFMETCollection> pfMET;
660     iEvent.getByLabel("pfMet", pfMET);
661    
662     if (pfMET.isValid())
663     {
664     MET.SetEt(pfMET->begin()->et());
665     MET.SetPhi(pfMET->begin()->phi());
666     }
667    
668     //Taus
669 devildog 1.4 // cout << "Taus... ";
670 devildog 1.1
671     edm::Handle<reco::PFTauCollection> PFtauHPS;
672     iEvent.getByLabel("hpsPFTauProducer", PFtauHPS);
673 devildog 1.4 edm::Handle<reco::PFTauDiscriminator> hps11;
674     iEvent.getByLabel("hpsPFTauDiscriminationByLooseElectronRejection", hps11);
675     edm::Handle<reco::PFTauDiscriminator> hps21;
676     iEvent.getByLabel("hpsPFTauDiscriminationByLooseMuonRejection", hps21);
677     edm::Handle<reco::PFTauDiscriminator> hps12;
678     iEvent.getByLabel("hpsPFTauDiscriminationByMediumElectronRejection", hps12);
679     // edm::Handle<reco::PFTauDiscriminator> hps22;
680     // iEvent.getByLabel("hpsPFTauDiscriminationByMediumMuonRejection", hps22);
681     edm::Handle<reco::PFTauDiscriminator> hps13;
682     iEvent.getByLabel("hpsPFTauDiscriminationByTightElectronRejection", hps13);
683     edm::Handle<reco::PFTauDiscriminator> hps23;
684     iEvent.getByLabel("hpsPFTauDiscriminationByTightMuonRejection", hps23);
685 devildog 1.1 edm::Handle<reco::PFTauDiscriminator> hps3;
686     iEvent.getByLabel("hpsPFTauDiscriminationByDecayModeFinding", hps3);
687     edm::Handle<reco::PFTauDiscriminator> hps4;
688     iEvent.getByLabel("hpsPFTauDiscriminationByLooseIsolation", hps4);
689     edm::Handle<reco::PFTauDiscriminator> hps5;
690     iEvent.getByLabel("hpsPFTauDiscriminationByMediumIsolation", hps5);
691     edm::Handle<reco::PFTauDiscriminator> hps6;
692     iEvent.getByLabel("hpsPFTauDiscriminationByTightIsolation", hps6);
693    
694     int iTau = 0;
695     for(PFTauCollection::const_iterator tau_iter = PFtauHPS->begin(); tau_iter!= PFtauHPS->end(); ++tau_iter){
696     reco::PFTau myTau = reco::PFTau(*tau_iter);
697     reco::PFTauRef theTauJetRef(PFtauHPS, iTau);
698     TCTau* mTau = new((*taus)[iTau]) TCTau;
699     mTau->SetAlgorithm(2);
700     mTau->SetP4(myTau.px(), myTau.py(), myTau.pz(), myTau.energy());
701 devildog 1.4 mTau->SethpsPFTauDiscriminationAgainstElectronLoose((*hps11)[theTauJetRef]);
702     mTau->SethpsPFTauDiscriminationAgainstMuonLoose((*hps21)[theTauJetRef]);
703     mTau->SethpsPFTauDiscriminationAgainstElectronMedium((*hps12)[theTauJetRef]);
704     // mTau->SethpsPFTauDiscriminationAgainstMuonMedium((*hps22)[theTauJetRef]);
705     mTau->SethpsPFTauDiscriminationAgainstElectronTight((*hps13)[theTauJetRef]);
706     mTau->SethpsPFTauDiscriminationAgainstMuonTight((*hps23)[theTauJetRef]);
707 devildog 1.1 mTau->SethpsPFTauDiscriminationByDecayModeFinding((*hps3)[theTauJetRef]);
708     mTau->SethpsPFTauDiscriminationByLooseIsolation((*hps4)[theTauJetRef]);
709     mTau->SethpsPFTauDiscriminationByMediumIsolation((*hps5)[theTauJetRef]);
710     mTau->SethpsPFTauDiscriminationByTightIsolation((*hps6)[theTauJetRef]);
711     mTau->SetCharge(myTau.charge());
712     mTau->SetNTracks(myTau.signalPFChargedHadrCands().size());
713     mTau->SetNConst(myTau.signalPFCands().size());
714 devildog 1.4 mTau->SetDecayMode((int)myTau.decayMode());
715    
716     if(!myTau.signalPFChargedHadrCands().isNull()) {
717     for(PFCandidateRefVector::const_iterator chit = myTau.signalPFChargedHadrCands().begin();
718     chit != myTau.signalPFChargedHadrCands().end(); chit++){
719     TLorentzVector temp((*chit)->px(),(*chit)->py(),(*chit)->pz(),(*chit)->energy());
720     mTau->AddChargedTrack(temp);
721     }
722     }
723     if(!myTau.signalPFNeutrHadrCands().isNull()) {
724     for(PFCandidateRefVector::const_iterator h0it = myTau.signalPFNeutrHadrCands().begin();
725     h0it != myTau.signalPFNeutrHadrCands().end(); h0it++){
726     TLorentzVector temp((*h0it)->px(),(*h0it)->py(),(*h0it)->pz(),(*h0it)->energy());
727     mTau->AddH0Track(temp);
728     }
729     }
730     if(!myTau.signalPFGammaCands().isNull()) {
731     for(PFCandidateRefVector::const_iterator egit = myTau.signalPFGammaCands().begin();
732     egit != myTau.signalPFGammaCands().end(); egit++){
733     TLorentzVector temp((*egit)->px(),(*egit)->py(),(*egit)->pz(),(*egit)->energy());
734     mTau->AddEMTrack(temp);
735     }
736     }
737    
738 devildog 1.1 iTau++;
739     }
740    
741    
742     //Need to store trigger information here.
743     //
744 devildog 1.4 // cout << "Triggers...\n ";
745 devildog 1.1
746 devildog 1.4 edm::Handle<TriggerResults> hltR;
747     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltProcess_);
748     iEvent.getByLabel(triggerResultsTag_,hltR);
749 devildog 1.1
750 devildog 1.4 const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
751     hlNames=triggerNames.triggerNames();
752 devildog 1.1
753 devildog 1.4 bool triggerPassed = false;
754     for(int i=0; i < (int)hlNames.size(); ++i) {
755 devildog 1.1 triggerPassed = triggerDecision(hltR, i);
756     if(triggerPassed){
757     int thisPrescale = hltConfig_.prescaleValue(iEvent, iSetup,hlNames[i]);
758     triggers.addTrigger(hlNames[i],thisPrescale);
759     }
760 devildog 1.4 }
761 devildog 1.1
762    
763    
764     if(
765     ( passmuons >= 1 ) ||
766     ( passelectrons >= 1 ) ||
767     ( keepAllEvents )
768     )
769     sTree->Fill();
770    
771     jet_ak5PF->Clear();
772     primaryVtx->Clear();
773     electrons->Clear();
774     muons->Clear();
775     taus->Clear();
776 devildog 1.4 genParticles->Clear();
777 devildog 1.1 triggers.clearTriggers();
778     }
779    
780     void Thesis::beginRun(const edm::Run& iRun, const edm::EventSetup& iEvent)
781     {
782     bool changed = true;
783     hltConfig_.init(iRun, iEvent, hltProcess_, changed);
784     }
785    
786    
787     // ------------ method called once each job just before starting event loop ------------
788     void
789     Thesis::beginJob()
790     {
791     ntupleFile = new TFile(rootfilename.c_str(), "RECREATE");
792     sTree = new TTree("myLepTree", "Tree for Leptons and Jets");
793 devildog 1.3 beamSpot = new TVector3();
794 devildog 1.1
795     jet_ak5PF = new TClonesArray("TCJet");
796     primaryVtx = new TClonesArray("TCPrimaryVtx");
797     electrons = new TClonesArray("TCElectron");
798     muons = new TClonesArray("TCMuon");
799     taus = new TClonesArray("TCTau");
800 devildog 1.2 genParticles = new TClonesArray("TGenParticle");
801 devildog 1.1
802 devildog 1.3 sTree->Branch("jet_ak5PF", &jet_ak5PF, 3200, 0);
803     sTree->Branch("primaryVtx", &primaryVtx, 3200, 0);
804     sTree->Branch("electrons", &electrons, 3200, 0);
805     sTree->Branch("muons", &muons, 3200, 0);
806     sTree->Branch("MET", &MET, 3200, 0);
807     sTree->Branch("taus", &taus, 3200, 0);
808     sTree->Branch("genParticles", &genParticles, 3200, 0);
809     sTree->Branch("triggers", &triggers, 3200,0);
810     sTree->Branch("beamSpot", &beamSpot, 3200,0);
811    
812 devildog 1.4 sTree->Branch("isMC", &isMC, "isMC/i");
813     sTree->Branch("eventNumber",&eventNumber, "eventNumber/i");
814     sTree->Branch("runNumber",&runNumber, "runNumber/i");
815     sTree->Branch("lumiSection",&lumiSection, "lumiSection/i");
816     sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
817    
818 devildog 1.1 }
819    
820     // ------------ method called once each job just after ending the event loop ------------
821     void
822     Thesis::endJob()
823     {
824     ntupleFile->Write();
825     ntupleFile->Close();
826     }
827    
828     double Thesis::sumPtSquared(const Vertex & v)
829     {
830     double sum = 0.;
831     double pT;
832     for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
833     pT = (**it).pt();
834     double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
835    
836     sum += pT*pT;
837     }
838     return sum;
839     }
840    
841     bool Thesis::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
842     {
843     bool triggerPassed = false;
844     if(hltR->wasrun(iTrigger) &&
845     hltR->accept(iTrigger) &&
846     !hltR->error(iTrigger) ){
847     triggerPassed = true;
848     }
849     return triggerPassed;
850     }
851    
852     //define this as a plug-in
853     DEFINE_FWK_MODULE(Thesis);