ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.25
Committed: Tue Nov 30 14:08:14 2010 UTC (14 years, 5 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.24: +18 -8 lines
Log Message:
Add genJet container to handle extra variables

File Contents

# User Rev Content
1 naodell 1.1 // Original Author: "Nathaniel Odell"
2     // Secondary Author: Steven Won
3 andrey 1.2 // With contributions from: Andrey Pozdnyakov
4 naodell 1.1 // Created: Thurs April 22 2010
5 naodell 1.25 // $Id: MPIntuple.cc,v 1.24 2010/11/29 13:09:04 naodell Exp $
6 naodell 1.1
7     // system include files
8     #include <memory>
9     #include <string>
10    
11     // user include files
12     #include "FWCore/Framework/interface/Frameworkfwd.h"
13     #include "FWCore/Framework/interface/EDAnalyzer.h"
14     #include "FWCore/Framework/interface/ESHandle.h"
15     #include "FWCore/Framework/interface/EventSetup.h"
16     #include "FWCore/Framework/interface/Event.h"
17     #include "FWCore/Framework/interface/MakerMacros.h"
18     #include "FWCore/ParameterSet/interface/ParameterSet.h"
19     #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
20     #include "Geometry/Records/interface/CaloGeometryRecord.h"
21 naodell 1.10 #include "DataFormats/Math/interface/deltaPhi.h"
22 naodell 1.1
23 naodell 1.15 // Libraries for objects
24 naodell 1.1 #include "DataFormats/JetReco/interface/CaloJet.h"
25     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
26     #include "DataFormats/JetReco/interface/PFJet.h"
27     #include "DataFormats/JetReco/interface/PFJetCollection.h"
28     #include "DataFormats/JetReco/interface/GenJet.h"
29     #include "DataFormats/JetReco/interface/GenJetCollection.h"
30     #include "DataFormats/JetReco/interface/Jet.h"
31     #include "DataFormats/VertexReco/interface/Vertex.h"
32     #include "DataFormats/VertexReco/interface/VertexFwd.h"
33 naodell 1.7 #include "DataFormats/BTauReco/interface/JetTag.h"
34 naodell 1.15 #include "DataFormats/TrackReco/interface/Track.h"
35 naodell 1.16 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
36 naodell 1.18 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
37 naodell 1.24 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
38 naodell 1.18
39 devildog 1.11 //#include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h"
40 naodell 1.7
41 naodell 1.10 //GenParticles
42 naodell 1.24 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
43 naodell 1.1
44 naodell 1.4 // Need to get correctors
45     #include "JetMETCorrections/Objects/interface/JetCorrector.h"
46    
47     // ntuple storage classes
48     #include "TCJet.h"
49     #include "TCPrimaryVtx.h"
50 naodell 1.25 #include "TCGenJet.h"
51 naodell 1.4
52 andrey 1.2 // Need for HLT trigger info:
53     #include "FWCore/Common/interface/TriggerNames.h"
54     #include "DataFormats/Common/interface/TriggerResults.h"
55     #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
56    
57 naodell 1.4 //Root stuff
58 naodell 1.1 #include "TROOT.h"
59 naodell 1.15 #include "TH1.h"
60 naodell 1.16 #include "TH2.h"
61 naodell 1.15 #include "TProfile.h"
62 naodell 1.1 #include "TFile.h"
63     #include "TTree.h"
64     #include "TString.h"
65     #include "TObject.h"
66     #include "TObjArray.h"
67     #include "TClonesArray.h"
68     #include "TRefArray.h"
69     #include "TLorentzVector.h"
70     #include "TVector3.h"
71    
72     using namespace edm;
73     using namespace std;
74     using namespace reco;
75    
76     //
77     // class declaration
78     //
79    
80     class MPIntuple : public edm::EDAnalyzer {
81     public:
82     explicit MPIntuple(const edm::ParameterSet&);
83     ~MPIntuple();
84    
85    
86     private:
87     virtual void beginJob() ;
88     virtual void analyze(const edm::Event&, const edm::EventSetup&);
89     virtual void endJob() ;
90    
91 andrey 1.2 bool triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger);
92 devildog 1.11 double sumPtSquared(const Vertex & v);
93 andrey 1.2
94 naodell 1.1 // ----------member data ---------------------------
95    
96     edm::InputTag PFJetHandle_;
97     edm::InputTag GenJetHandle_;
98     edm::InputTag PrimaryVtxHandle_;
99 naodell 1.4 edm::InputTag triggerResultsTag_;
100 naodell 1.1 // Counting variables //
101    
102 naodell 1.8 int eventNumber, runNumber, lumiSection, bunchCross;
103 naodell 1.24 double ptHat, qScale, crossSection;
104 naodell 1.18
105 naodell 1.6 TTree* sTree;
106 naodell 1.1 TFile* ntupleFile;
107    
108 naodell 1.4 TClonesArray* jet_ak5PF;
109 naodell 1.25 TClonesArray* jet_ak5Gen;
110 naodell 1.5 TClonesArray* primaryVtx;
111 naodell 1.1
112     bool doGenJets_;
113     bool doPFJets_;
114 naodell 1.7 bool triggerHLT_;
115 naodell 1.8 bool isRealData;
116 naodell 1.1 string rootfilename;
117    
118 andrey 1.2 //Triggers
119 naodell 1.4 string hlTriggerResults_, hltName_, triggerName_;
120     TriggerNames triggerNames;
121 andrey 1.2 HLTConfigProvider hltConfig_;
122 naodell 1.4 vector<string> hlNames;
123 naodell 1.5 unsigned int triggerStatus;
124 andrey 1.2
125 naodell 1.15 //Histograms
126 naodell 1.24 TH1D * h1_fracAssociatedTracks;
127     TH1D * h1_meanJetTrackZ;
128     TH1D * h1_trackDxy;
129     TH1D * h1_jetVertexZ;
130     TH1D * h1_associatedSumPt;
131     TH1D * h1_associatedVertexIndex;
132 naodell 1.16 TH2F * h2_nAssTracksVsJetPt;
133 naodell 1.15 TProfile * p1_nVtcs;
134    
135    
136 naodell 1.1 };
137    
138 naodell 1.18 MPIntuple::MPIntuple(const edm::ParameterSet& iConfig)
139 naodell 1.1 {
140 naodell 1.18 PFJetHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag");
141     GenJetHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag");
142     PrimaryVtxHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag");
143     doGenJets_ = iConfig.getUntrackedParameter<bool>("doGenJets");
144     doPFJets_ = iConfig.getUntrackedParameter<bool>("doPFJets");
145     triggerHLT_ = iConfig.getUntrackedParameter<bool>("triggerHLT");
146     rootfilename = iConfig.getUntrackedParameter<string>("rootfilename");
147     hlTriggerResults_ = iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults");
148 naodell 1.24 hltName_ = iConfig.getUntrackedParameter<string>("hltName");
149 naodell 1.1 }
150    
151     MPIntuple::~MPIntuple()
152     {
153    
154     }
155    
156     //
157     // member functions
158     //
159    
160     // ------------ method called to for each event ------------
161     void MPIntuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
162     {
163    
164 naodell 1.6 eventNumber = iEvent.id().event();
165     runNumber = iEvent.id().run();
166 naodell 1.1 lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
167 naodell 1.8 bunchCross = (unsigned int)iEvent.bunchCrossing();
168     isRealData = iEvent.isRealData();
169    
170 naodell 1.18
171 naodell 1.16 edm::Handle<reco::BeamSpot> beamSpotHandle;
172     iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
173     reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
174    
175     int pfCount = 0;
176     int genCount = 0;
177     int vtxCount = 0;
178 naodell 1.24 double primaryVertexZ = -999;
179 naodell 1.16
180     //////////////////////////
181     //Get vertex information//
182     //////////////////////////
183    
184     Handle<reco::VertexCollection> primaryVtcs;
185     iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
186    
187     for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
188    
189     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
190     if(!myVtx.isValid() || myVtx.isFake()) continue;
191     TCPrimaryVtx* vtxCon = new ((*primaryVtx)[vtxCount]) TCPrimaryVtx;
192    
193 naodell 1.23 //cout<<myVtx.nTracks()<<endl;
194    
195 naodell 1.16 vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
196     vtxCon->SetNDof(myVtx.ndof());
197     vtxCon->SetChi2(myVtx.chi2());
198     vtxCon->SetNTrks(myVtx.tracksSize());
199     vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
200    
201     if(vtxCount == 0) primaryVertexZ = myVtx.z();
202    
203     ++vtxCount;
204    
205     }
206    
207     p1_nVtcs->Fill(runNumber, vtxCount);
208    
209     ///////////////////////
210     //get jet information//
211     ///////////////////////
212 naodell 1.1
213     if(doPFJets_){
214 naodell 1.4
215 naodell 1.7 edm::Handle<reco::JetTagCollection> bTagHandle1;
216     iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
217     const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
218     reco::JetTagCollection::const_iterator jet_it_1;
219    
220     edm::Handle<reco::JetTagCollection> bTagHandle2;
221     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
222     const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
223     reco::JetTagCollection::const_iterator jet_it_2;
224    
225 naodell 1.14
226 naodell 1.4 const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
227     const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
228    
229 naodell 1.1 Handle<reco::PFJetCollection> PFJets;
230     iEvent.getByLabel(PFJetHandle_, PFJets);
231 naodell 1.16
232 naodell 1.1 for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
233    
234     reco::PFJet myJet = reco::PFJet(*jet_iter);
235 naodell 1.22 if(myJet.pt() > 5){
236 naodell 1.7
237     for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
238     if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
239     }
240    
241     for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
242     if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
243     }
244 naodell 1.6
245 naodell 1.16 TCJet* jetCon = new ((*jet_ak5PF)[pfCount]) TCJet;
246 naodell 1.6
247     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
248 naodell 1.15 jetCon->SetVtx(-999.0, -999.0, -999.0);
249 naodell 1.4
250 naodell 1.6 jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
251     jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
252     jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
253     jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
254 naodell 1.4
255 naodell 1.6 jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
256     jetCon->SetNumChPart(myJet.chargedMultiplicity());
257 naodell 1.4
258 naodell 1.6 jetCon->SetNumChPart(myJet.chargedMultiplicity());
259 naodell 1.4
260 naodell 1.7 if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
261     if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
262    
263 naodell 1.24 double scale2 = correctorL2->correction(myJet);
264 naodell 1.7 myJet.scaleEnergy(scale2);
265 naodell 1.24 double scale3 = correctorL3->correction(myJet);
266 naodell 1.7 myJet.scaleEnergy(scale3);
267    
268     //more corrections?
269    
270 naodell 1.6 jetCon->SetJetCorr(2, scale2);
271     jetCon->SetJetCorr(3, scale3);
272 naodell 1.7
273 naodell 1.16 /////////////////////////
274     //get associated tracks//
275     /////////////////////////
276 naodell 1.15
277     const reco::TrackRefVector &tracks = myJet.getTrackRefs();
278    
279 naodell 1.24 vector<TVector3> vtxPositionCollection;
280     vector<double> associatedTrackSumPt;
281     vector<unsigned int> jetTrackAddresses;
282     double sumTrackX, sumTrackY, sumTrackZ, sumTrackIP, sumTrackPt;
283     //double meanTrackX, meanTrackY, meanTrackZ, meanTrackIP;
284     int nJetTracks, nVertexTracks, nAssociatedTracks, vertexIndex;
285 naodell 1.23 int vCount = 0;
286 naodell 1.16
287 naodell 1.24 nJetTracks = nVertexTracks = nAssociatedTracks = 0;
288     sumTrackX = sumTrackY = sumTrackZ = sumTrackIP = sumTrackPt = 0;
289     //meanTrackZ = meanTrackIP = -999;
290 naodell 1.15
291 naodell 1.23 if(myJet.pt() > 10 && fabs(myJet.eta()) < 2.4){
292 naodell 1.17
293 naodell 1.15 for(TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack){
294    
295 naodell 1.23 const reco::Track &myJetTrack = **iTrack;
296 naodell 1.20
297 naodell 1.23 sumTrackPt += myJetTrack.pt();
298 naodell 1.24 sumTrackX += myJetTrack.vx();
299     sumTrackY += myJetTrack.vy();
300 naodell 1.23 sumTrackZ += myJetTrack.vz();
301 naodell 1.21 sumTrackIP += myJetTrack.dxy(vertexBeamSpot.position());
302 naodell 1.24 jetTrackAddresses.push_back((unsigned int)&myJetTrack);
303     ++nJetTracks;
304     }
305    
306     if(nJetTracks > 0){
307    
308     jetCon->SetVtx(sumTrackX/nJetTracks, sumTrackY/nJetTracks, sumTrackZ/nJetTracks);
309    
310     h1_meanJetTrackZ->Fill(sumTrackZ/nJetTracks);
311     h1_trackDxy->Fill(sumTrackIP/nJetTracks);
312    
313 naodell 1.23 }
314 naodell 1.24 /*meanTrackX = sumTrackX/nJetTracks;
315     meanTrackY = sumTrackY/nJetTracks;
316     meanTrackZ = sumTrackZ/nJetTracks;
317     */
318    
319     if(jetTrackAddresses.size() > 0){
320 naodell 1.23
321     for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
322    
323 naodell 1.24 reco::Vertex myVtx = reco::Vertex(*vtx_iter);
324 naodell 1.23 if(!myVtx.isValid() || myVtx.isFake()) continue;
325    
326 naodell 1.24 TVector3 *iVtxPosition = new TVector3(myVtx.x(), myVtx.y(), myVtx.z());
327     vtxPositionCollection.push_back(*iVtxPosition);
328     associatedTrackSumPt.push_back(0);
329    
330 naodell 1.23 for(Vertex::trackRef_iterator iTrackRef = myVtx.tracks_begin(); iTrackRef != myVtx.tracks_end(); ++iTrackRef){
331 naodell 1.24 const edm::RefToBase<reco::Track> &myTrackRef = *iTrackRef;
332 naodell 1.23
333     if(myTrackRef.isAvailable()){
334     const reco::Track &myVertexTrack = *myTrackRef.get();
335    
336 naodell 1.24 for(vector<unsigned int>::const_iterator iTrackAddress = jetTrackAddresses.begin(); iTrackAddress != jetTrackAddresses.end(); ++iTrackAddress){
337    
338     if(*iTrackAddress == (unsigned int)&myVertexTrack){
339 naodell 1.23
340 naodell 1.24 associatedTrackSumPt.at(vCount) += myVertexTrack.pt()/sumTrackPt;
341     ++nAssociatedTracks;
342    
343     }
344     }
345 naodell 1.23 }
346 naodell 1.24 }
347     ++vCount;
348     }
349    
350     double maxSumPtFraction = 0;
351     vCount = vertexIndex = 0;
352    
353     for(vector<double>::const_iterator iTrackSumPt = associatedTrackSumPt.begin(); iTrackSumPt != associatedTrackSumPt.end(); ++iTrackSumPt) {
354    
355     if(*iTrackSumPt > maxSumPtFraction){
356    
357     maxSumPtFraction = *iTrackSumPt;
358     vertexIndex = vCount + 1;
359    
360     }
361     ++vCount;
362     }
363    
364     if(vertexIndex > 0){
365    
366     h1_jetVertexZ->Fill(vtxPositionCollection[vertexIndex-1].z());
367    
368     }
369    
370     jetCon->SetVtxSumPt(maxSumPtFraction);
371     jetCon->SetVtxIndex(vertexIndex);
372    
373     h1_fracAssociatedTracks->Fill((double)nAssociatedTracks/(double)nJetTracks);
374     h1_associatedSumPt->Fill(maxSumPtFraction);
375     h1_associatedVertexIndex->Fill(vertexIndex);
376     h2_nAssTracksVsJetPt->Fill(myJet.pt(), nAssociatedTracks);
377 naodell 1.23 }
378 naodell 1.24 }
379 naodell 1.16 ++pfCount;
380 naodell 1.4 }
381     }
382 naodell 1.1 }
383    
384 naodell 1.10 if(!isRealData){
385 naodell 1.1
386 naodell 1.24 Handle<GenEventInfoProduct> GenEventInfoHandle;
387     iEvent.getByLabel("generator", GenEventInfoHandle);
388    
389     Handle<GenRunInfoProduct> GenRunInfoHandle;
390     iEvent.getByLabel("generator", GenRunInfoHandle);
391 naodell 1.18
392 naodell 1.1 Handle<reco::GenJetCollection> GenJets;
393     iEvent.getByLabel(GenJetHandle_, GenJets);
394 naodell 1.18
395 naodell 1.24 ptHat = qScale = -1; crossSection = 0;
396 naodell 1.18
397 naodell 1.24 if(GenEventInfoHandle.isValid())
398     {
399     qScale = GenEventInfoHandle->qScale();
400     ptHat = (GenEventInfoHandle->hasBinningValues() ? GenEventInfoHandle->binningValues()[0] : 0.0);
401     }
402     if(GenRunInfoHandle.isValid())
403     {
404     crossSection = GenRunInfoHandle->crossSection();
405 naodell 1.18 }
406 naodell 1.16
407 naodell 1.1 for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
408    
409     reco::GenJet myJet = reco::GenJet(*jet_iter);
410    
411 naodell 1.4 if(myJet.pt() > 5){
412 naodell 1.1
413 naodell 1.25 TCGenJet* jetCon = new ((*jet_ak5Gen)[genCount]) TCGenJet;
414    
415     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
416    
417     jetCon->SetHadEnergy(myJet.hadEnergy());
418     jetCon->SetEmEnergy(myJet.emEnergy());
419     jetCon->SetInvEnergy(myJet.invisibleEnergy());
420     jetCon->SetAuxEnergy(myJet.auxiliaryEnergy());
421     jetCon->SetNumConstit(myJet.getGenConstituents().size());
422 naodell 1.24
423     const reco::GenParticle *myGenParticle = myJet.getGenConstituent(0);
424     const reco::Candidate *myMother = myGenParticle->mother();
425     int nGrandMas = 1;
426    
427     for(int iter = 0; nGrandMas != 0; ++iter)
428     {
429     myMother = myMother->mother();
430 naodell 1.25 nGrandMas = myMother->mother()->numberOfMothers();
431 naodell 1.24 }
432     //cout<<myJet.print()<<endl;
433 naodell 1.25 //cout<<"Jet flavor: "<< myMother->pdgId()<<endl;
434     jetCon->SetJetFlavor(myMother->pdgId());
435 naodell 1.16 ++genCount;
436 naodell 1.1
437     }
438 naodell 1.13 }
439 naodell 1.1 }
440 naodell 1.10
441 naodell 1.16 ///////////////////////////
442     //get trigger information//
443     ///////////////////////////
444 andrey 1.2
445 naodell 1.7 if(triggerHLT_){
446 andrey 1.2
447 naodell 1.7 edm::Handle<TriggerResults> hltR;
448     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
449     iEvent.getByLabel(triggerResultsTag_,hltR);
450    
451     const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
452     hlNames=triggerNames.triggerNames();
453    
454 naodell 1.24 string MPI_TriggerNames[] = {"HLT_PixelTracks_Multiplicity70",
455     "HLT_MinBiasBSC_NoBPTX",
456     "HLT_PixelTracks_Multiplicity40",
457     "HLT_L1Tech_HCAL_HF",
458     "HLT_IsoTrackHB_8E29",
459     "HLT_IsoTrackHE_8E29",
460     "HLT_L1Tech_RPC_TTU_RBst1_collisions",
461     "HLT_L1_BscMinBiasOR_BptxPlusORMinus",
462     "HLT_L1Tech_BSC_halo_forPhysicsBackground",
463     "HLT_L1Tech_BSC_HighMultiplicity",
464     "HLT_MinBiasPixel_DoubleIsoTrack5",
465     "HLT_MinBiasPixel_DoubleTrack",
466     "HLT_MinBiasPixel_SingleTrack",
467     "HLT_ZeroBiasPixel_SingleTrack",
468     "HLT_MinBiasBSC",
469     "HLT_StoppedHSCP_8E29",
470     "HLT_Jet15U_HcalNoiseFiltered",
471     "HLT_QuadJet15U",
472     "HLT_DiJetAve30U_8E29",
473     "HLT_DiJetAve15U_8E29",
474     "HLT_FwdJet20U",
475     "HLT_Jet50U",
476     "HLT_Jet30U",
477     "HLT_Jet15U",
478     "HLT_BTagMu_Jet10U",
479     "HLT_DoubleJet15U_ForwardBackward",
480     "HLT_BTagIP_Jet50U",
481     "HLT_DoubleLooseIsoTau15",
482     "HLT_SingleLooseIsoTau20",
483     "HLT_HT100U",
484     "HLT_MET100",
485     "HLT_MET45"};
486 naodell 1.7
487     bool triggerPassed = false;
488     triggerStatus = 0x0;
489 naodell 1.4
490 naodell 1.7 for (uint i=0; i<hlNames.size(); ++i) {
491    
492     triggerPassed = triggerDecision(hltR, i);
493 naodell 1.4
494 naodell 1.7 if(triggerPassed){
495 andrey 1.3
496 naodell 1.7 for (uint j = 0; j != 32; ++j){
497 naodell 1.24
498 naodell 1.7 if (hlNames[i] == MPI_TriggerNames[j])
499 naodell 1.24 {
500     //cout<<"HLTrigger name: "<<hlNames[i]<<" bit: "<<dec<<j+1<<endl;
501     triggerStatus |= 0x01 << j;
502     }
503 naodell 1.7 }
504 naodell 1.4 }
505 naodell 1.7 }
506     }
507 naodell 1.16 //cout<< "total status: "<<hex<<triggerStatus<<endl;
508 naodell 1.18
509 naodell 1.14
510    
511 naodell 1.24 if((pfCount > 1 || genCount > 1) && vtxCount > 0) sTree -> Fill();
512 naodell 1.15
513 naodell 1.4 jet_ak5PF->Clear();
514 naodell 1.5 primaryVtx->Clear();
515 naodell 1.25 jet_ak5Gen->Clear();
516 andrey 1.2
517 naodell 1.1 }
518    
519    
520     // ------------ method called once each job just before starting event loop ------------
521     void MPIntuple::beginJob()
522     {
523    
524 naodell 1.24 ntupleFile = new TFile(rootfilename.c_str(), "RECREATE");
525     sTree = new TTree("mpiTree", "Tree for Jets");
526 naodell 1.15
527 naodell 1.24 h1_fracAssociatedTracks = new TH1D("h1_fracAssociatedTracks", "Fraction of associated jet tracks", 20, 0.0, 1.0);
528     h1_trackDxy = new TH1D("h1_trackDxy", "dxy for all associated tracks", 50, -0.4, 0.4);
529     h1_meanJetTrackZ = new TH1D("h1_meanJetTrackZ", "<z> from jet tracks", 50, -25.0, 25.0);
530     h1_jetVertexZ = new TH1D("h1_jetVertexZ", "z for vertex associated to jet", 50, -25.0, 25.0);
531     h1_associatedSumPt = new TH1D("h1_associatedSumPt", "ratio of sum p_{T} between associatedTracks and jetTracks", 20, 0.0, 1.0);
532     h1_associatedVertexIndex = new TH1D("h1_associatedVertex", "number of jets associated to primary vertex", 5, -0.5, 4.5);
533     h2_nAssTracksVsJetPt = new TH2F("h2_nAssTracksVsJetPt", "nTracks vs Jet p_{T}", 50, 0, 100, 20, -0.5, 19.5);
534     p1_nVtcs = new TProfile("p1_nVtcs", "<N> per run", 5200.0, 132300.0, 137500.0, 0.0, 6.0);
535    
536     jet_ak5PF = new TClonesArray("TCJet");
537 naodell 1.25 jet_ak5Gen = new TClonesArray("TCGenJet");
538 naodell 1.24 primaryVtx = new TClonesArray("TCPrimaryVtx");
539 naodell 1.1
540 naodell 1.15
541 naodell 1.4 sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0);
542 naodell 1.25 sTree->Branch("jet_ak5Gen",&jet_ak5Gen, 6400, 0);
543 naodell 1.5 sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
544 naodell 1.1
545     sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
546     sTree->Branch("runNumber",&runNumber, "runNumber/I");
547     sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
548 naodell 1.6 sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
549 naodell 1.8 sTree->Branch("isRealData",&isRealData, "isRealData/i");
550     sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
551 naodell 1.18 sTree->Branch("ptHat",&ptHat, "ptHat/d");
552 naodell 1.24 sTree->Branch("qScale", &qScale, "qScale/d");
553     sTree->Branch("crossSection", &crossSection, "crossSection/d");
554 naodell 1.1 }
555    
556     // ------------ method called once each job just after ending the event loop ------------
557     void MPIntuple::endJob()
558     {
559 naodell 1.15
560 naodell 1.17 ntupleFile->cd();
561    
562 naodell 1.24 h1_fracAssociatedTracks->Write();
563     h1_meanJetTrackZ->Write();
564 naodell 1.16 h1_trackDxy->Write();
565 naodell 1.24 h1_jetVertexZ->Write();
566     h1_associatedSumPt->Write();
567     h1_associatedVertexIndex->Write();
568 naodell 1.16 h2_nAssTracksVsJetPt->Write();
569 naodell 1.15 p1_nVtcs->Write();
570    
571 naodell 1.1 ntupleFile->Write();
572     ntupleFile->Close();
573    
574     }
575    
576    
577 andrey 1.2 bool MPIntuple::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
578     {
579     bool triggerPassed = false;
580     if(hltR->wasrun(iTrigger) &&
581     hltR->accept(iTrigger) &&
582     !hltR->error(iTrigger) ){
583     triggerPassed = true;
584     }
585     return triggerPassed;
586 naodell 1.16 }
587    
588     double MPIntuple::sumPtSquared(const Vertex & v)
589     {
590     double sum = 0.;
591     double pT;
592     for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
593     pT = (**it).pt();
594     double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
595    
596     sum += pT*pT;
597 andrey 1.2 }
598 naodell 1.16 return sum;
599     }
600 andrey 1.2
601 naodell 1.1 //define this as a plug-in
602     DEFINE_FWK_MODULE(MPIntuple);