ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
(Generate patch)

Comparing UserCode/MPIAnalyzer/src/MPIntuple.cc (file contents):
Revision 1.3 by andrey, Sat May 22 00:51:59 2010 UTC vs.
Revision 1.26 by naodell, Mon Dec 6 19:36:19 2010 UTC

# Line 1 | Line 1
1 /*
2 Description: <one line class summary>
3
4 Implementation:
5     n-tuple creator for jet studies
6 */
7 //
1   // Original Author:  "Nathaniel Odell"
2   // Secondary Author: Steven Won
3   // With contributions from: Andrey Pozdnyakov
4   //         Created:  Thurs April 22  2010
5   // $Id$
13 //
14 //
15
6  
7   // system include files
8   #include <memory>
# Line 25 | Line 15
15   #include "FWCore/Framework/interface/EventSetup.h"
16   #include "FWCore/Framework/interface/Event.h"
17   #include "FWCore/Framework/interface/MakerMacros.h"
18 <
18 > #include "FWCore/Framework/interface/LuminosityBlock.h"
19   #include "FWCore/ParameterSet/interface/ParameterSet.h"
20  
31 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
21   #include "Geometry/Records/interface/CaloGeometryRecord.h"
22 <
22 > #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
23 > #include "DataFormats/GeometryVector/interface/GlobalVector.h"
24 > #include "DataFormats/GeometryVector/interface/LocalPoint.h"
25 > #include "DataFormats/GeometryVector/interface/LocalVector.h"
26 > #include "DataFormats/Math/interface/Point3D.h"
27 > #include "DataFormats/Math/interface/Vector3D.h"
28 > #include "DataFormats/Math/interface/LorentzVector.h"
29 > #include "DataFormats/Math/interface/deltaPhi.h"
30 >
31 > // Libraries for objects
32 > #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
33 > #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
34 > #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
35   #include "DataFormats/JetReco/interface/CaloJet.h"
36   #include "DataFormats/JetReco/interface/CaloJetCollection.h"
36 #include "DataFormats/JetReco/interface/PFJet.h"
37 #include "DataFormats/JetReco/interface/PFJetCollection.h"
37   #include "DataFormats/JetReco/interface/GenJet.h"
38   #include "DataFormats/JetReco/interface/GenJetCollection.h"
39 < #include "DataFormats/JetReco/interface/Jet.h"
39 > #include "DataFormats/JetReco/interface/PFJet.h"
40 > #include "DataFormats/JetReco/interface/PFJetCollection.h"
41 > #include "DataFormats/MuonReco/interface/Muon.h"
42 > #include "DataFormats/MuonReco/interface/MuonFwd.h"
43 > #include "DataFormats/MuonReco/interface/MuonSelectors.h"
44 > #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
45 > #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
46 > #include "DataFormats/EgammaCandidates/interface/Photon.h"
47 > #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
48 > #include "DataFormats/TrackReco/interface/Track.h"
49 > #include "DataFormats/TrackReco/interface/TrackFwd.h"
50 > #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
51   #include "DataFormats/VertexReco/interface/Vertex.h"
52   #include "DataFormats/VertexReco/interface/VertexFwd.h"
53 + #include "DataFormats/BTauReco/interface/JetTag.h"
54 + #include "DataFormats/BeamSpot/interface/BeamSpot.h"
55 + #include "DataFormats/Luminosity/interface/LumiSummary.h"
56 + #include "DataFormats/Common/interface/ValueMap.h"
57 + #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
58 + #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
59 + #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
60 +
61 + //#include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h"
62 +
63 + // Need to get corrections
64 + #include "JetMETCorrections/Objects/interface/JetCorrector.h"
65 +
66 + // ntuple storage classes
67 + #include "TCJet.h"
68 + #include "TCPrimaryVtx.h"
69 + #include "TCGenJet.h"
70 + #include "TCMuon.h"
71 + #include "TCElectron.h"
72  
73   // Need for HLT trigger info:
74   #include "FWCore/Common/interface/TriggerNames.h"
75   #include "DataFormats/Common/interface/TriggerResults.h"
76   #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
48 //
77  
78 + //Root  stuff
79   #include "TROOT.h"
80 + #include "TH1.h"
81 + #include "TH2.h"
82 + #include "TProfile.h"
83   #include "TFile.h"
84   #include "TTree.h"
85   #include "TString.h"
# Line 74 | Line 106 | class MPIntuple : public edm::EDAnalyzer
106  
107     private:
108        virtual void beginJob() ;
109 +      virtual void beginRun(const edm::Run&, const edm::EventSetup&) ;
110        virtual void analyze(const edm::Event&, const edm::EventSetup&);
111 +      virtual void endLuminosityBlock(const edm::LuminosityBlock&,const edm::EventSetup&);
112 +      virtual void endRun(const edm::Run&, const edm::EventSetup&);
113        virtual void endJob() ;
114  
115 <  bool triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger);
115 >      bool triggerDecision(edm::Handle<edm::TriggerResults>& hltR, int iTrigger);
116 >      double sumPtSquared(const Vertex& v);
117  
118        // ----------member data ---------------------------
119  
84  edm::InputTag PFJetHandle_;
85  edm::InputTag CaloJetHandle_;
86  edm::InputTag PFJetHandle_corr;
87  edm::InputTag CaloJetHandle_corr;
88  edm::InputTag GenJetHandle_;
89  edm::InputTag PrimaryVtxHandle_;
90
91
92  //  Counting variables   //
93
94  int eventNumber, runNumber, lumiSection;
95
96  TTree *sTree;
97  TFile* ntupleFile;
98
99  TClonesArray* jetP4_ak5PF;
100  TClonesArray* jetP4_ak5PF_corr;
101  TClonesArray* jetVtx3_ak5PF;
102  TClonesArray* jetVtx3_ak5PF_corr;
103  TClonesArray* jetP4_ak5Calo;
104  TClonesArray* jetP4_ak5Calo_corr;
105  TClonesArray* jetVtx3_ak5Calo;
106  TClonesArray* jetVtx3_ak5Calo_corr;
107  TClonesArray* jetP4_ak5Gen;
108  TClonesArray* jetVtx3_ak5Gen;
109  TClonesArray* PrimaryVtx3;
110
111  bool doGenJets_;
112  bool doCaloJets_;
113  bool doPFJets_;
114
115  string rootfilename;
116
117  //Triggers
118  std::string hlTriggerResults_;
119  edm::TriggerNames triggerNames;
120  std::string hltName_;
121  std::string  triggerName_;
122  HLTConfigProvider hltConfig_;
123  edm::InputTag triggerResultsTag_;
124  std::vector<std::string>  hlNames;
125
120  
121 +      int eventNumber, runNumber, lumiSection, bunchCross;
122 +      double ptHat, qScale, crossSection;
123 +  
124 +      TTree* sTree;
125 +      TFile* ntupleFile;
126 +      edm::InputTag recoJetTag_;
127 +      edm::InputTag genJetTag_;
128 +      edm::InputTag muonTag_;
129 +      edm::InputTag electronTag_;
130 +      edm::InputTag primaryVtxTag_;
131 +      edm::InputTag triggerResultsTag_;
132 +      edm::InputTag electronIDMap_;
133 +      bool doGenJets_;
134 +      bool doPFJets_;
135 +      bool triggerHLT_;
136 +      bool isRealData;
137 +      string rootfilename;
138 +
139 +      TClonesArray* recoJets;
140 +      TClonesArray* genJets;
141 +      TClonesArray* primaryVtx;
142 +      TClonesArray* recoMuons;
143 +      TClonesArray* recoElectrons;
144 +
145 +      //Triggers
146 +      string hlTriggerResults_, hltProcess_, triggerName_;
147 +      TriggerNames triggerNames;
148 +      HLTConfigProvider hltConfig_;
149 +      vector<string>  hlNames;
150 +      unsigned int triggerStatus;
151 +
152 +      //Histograms
153 +      TH1D * h1_fracAssociatedTracks;
154 +      TH1D * h1_meanJetTrackZ;
155 +      TH1D * h1_trackDxy;
156 +      TH1D * h1_jetVertexZ;
157 +      TH1D * h1_associatedSumPt;
158 +      TH1D * h1_associatedVertexIndex;
159 +      TH2F * h2_nAssTracksVsJetPt;
160 +      TProfile * p1_nVtcs;
161   };
162  
163 < MPIntuple::MPIntuple(const edm::ParameterSet& iConfig):
130 <
131 <  PFJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag")),
132 <  CaloJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("CaloJetTag")),
133 <  PFJetHandle_corr(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag_corr")),
134 <  CaloJetHandle_corr(iConfig.getUntrackedParameter<edm::InputTag>("CaloJetTag_corr")),
135 <  GenJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag")),
136 <  PrimaryVtxHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")),
137 <  doGenJets_(iConfig.getUntrackedParameter<bool>("doGenJets")),
138 <  doCaloJets_(iConfig.getUntrackedParameter<bool>("doCaloJets")),
139 <  doPFJets_(iConfig.getUntrackedParameter<bool>("doPFJets")),
140 <  rootfilename(iConfig.getUntrackedParameter<string>("rootfilename")),
141 <  hlTriggerResults_(iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")),
142 <  hltName_(iConfig.getUntrackedParameter<string>("hltName","HLT"))
163 > MPIntuple::MPIntuple(const edm::ParameterSet& iConfig)
164   {
165 <  //edm::TriggerNames
166 <  // triggerNames(iConfig);
167 <
165 >  recoJetTag_       = iConfig.getUntrackedParameter<edm::InputTag>("RecoJetTag");
166 >  muonTag_          = iConfig.getUntrackedParameter<edm::InputTag>("MuonTag");
167 >  electronTag_      = iConfig.getUntrackedParameter<edm::InputTag>("ElectronTag");
168 >  genJetTag_        = iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag");
169 >  primaryVtxTag_    = iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag");
170 >  electronIDMap_    = iConfig.getParameter<edm::InputTag>("electronIDMap");
171 >  doGenJets_        = iConfig.getUntrackedParameter<bool>("doGenJets");
172 >  doPFJets_         = iConfig.getUntrackedParameter<bool>("doPFJets");
173 >  triggerHLT_       = iConfig.getUntrackedParameter<bool>("triggerHLT");
174 >  hlTriggerResults_ = iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults");
175 >  hltProcess_          = iConfig.getUntrackedParameter<string>("hltName");
176 >  rootfilename      = iConfig.getUntrackedParameter<string>("rootfilename");
177   }
178  
149
179   MPIntuple::~MPIntuple()
180   {
181  
182   }
183  
155
184   //
185   // member functions
186   //
# Line 161 | Line 189 | MPIntuple::~MPIntuple()
189   void MPIntuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
190   {
191  
192 <  eventNumber = iEvent.id().event();
193 <  runNumber   = iEvent.id().run();
192 >  eventNumber  = iEvent.id().event();
193 >  runNumber    = iEvent.id().run();
194    lumiSection  = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
195 +  bunchCross   = (unsigned int)iEvent.bunchCrossing();
196 +  isRealData   = iEvent.isRealData();
197 +
198 +
199 +  edm::Handle<reco::BeamSpot> beamSpotHandle;
200 +  iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
201 +  reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
202 +
203 +  int pfCount   = 0;
204 +  int genCount  = 0;
205 +  int vtxCount  = 0;
206 +  double primaryVertexZ = -999;
207 +
208 +  //////////////////////////
209 +  //Get vertex information//
210 +  //////////////////////////
211 +
212 +  Handle<reco::VertexCollection> primaryVtcs;
213 +  iEvent.getByLabel(primaryVtxTag_, primaryVtcs);
214    
215 <  int PFcount = 0;
216 <  int Calocount = 0;
217 <  int PFcount_corr = 0;
218 <  int Calocount_corr = 0;
219 <  int Gencount = 0;
220 <  int Vtxcount = 0;
215 >  for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
216 >    
217 >    reco::Vertex myVtx = reco::Vertex(*vtx_iter);
218 >    if(!myVtx.isValid() || myVtx.isFake()) continue;
219 >    TCPrimaryVtx* vtxCon = new ((*primaryVtx)[vtxCount]) TCPrimaryVtx;
220 >      
221 >    //cout<<myVtx.nTracks()<<endl;
222  
223 <  if(doPFJets_){
223 >    vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
224 >    vtxCon->SetNDof(myVtx.ndof());
225 >    vtxCon->SetChi2(myVtx.chi2());
226 >    vtxCon->SetNTrks(myVtx.tracksSize());
227 >    vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
228 >
229 >    if(vtxCount == 0) primaryVertexZ = myVtx.z();
230 >
231 >    ++vtxCount;
232      
233 +  }
234 +
235 +  p1_nVtcs->Fill(runNumber, vtxCount);
236 +
237 +  ///////////////////////
238 +  //get jet information//
239 +  ///////////////////////
240 +
241 +  if(doPFJets_){
242 +
243 +    edm::Handle<reco::JetTagCollection> bTagHandle1;
244 +    iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
245 +    const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
246 +    reco::JetTagCollection::const_iterator jet_it_1;
247 +
248 +    edm::Handle<reco::JetTagCollection> bTagHandle2;
249 +    iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
250 +    const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
251 +    reco::JetTagCollection::const_iterator jet_it_2;
252 +
253 +    const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
254 +    const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
255 +
256      Handle<reco::PFJetCollection> PFJets;
257 <    iEvent.getByLabel(PFJetHandle_, PFJets);
258 <    
180 <    PFcount = 0;
181 <    
257 >    iEvent.getByLabel(recoJetTag_, PFJets);
258 >        
259      for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
260        
261        reco::PFJet myJet = reco::PFJet(*jet_iter);
262 <      
263 <      if(myJet.et() > 5){
264 <        
265 <        new ((*jetP4_ak5PF)[PFcount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
266 <        new ((*jetVtx3_ak5PF)[PFcount]) TVector3(myJet.vx(), myJet.vy(), myJet.vz());
267 <        
268 <        ++PFcount;
269 <      }
270 <      
271 <    }
272 <    
273 <    
274 <    PFcount_corr = 0;
275 <    
276 <    Handle<reco::PFJetCollection> PFJets_corr;
277 <    iEvent.getByLabel(PFJetHandle_corr, PFJets_corr);
278 <    
279 <    for(PFJetCollection::const_iterator jet_iter = PFJets_corr->begin(); jet_iter!= PFJets_corr->end(); ++jet_iter){
280 <      
281 <      reco::PFJet myJet = reco::PFJet(*jet_iter);
282 <      
283 <      if(myJet.et() > 5){
284 <        
285 <        new ((*jetP4_ak5PF_corr)[PFcount_corr]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
286 <        new ((*jetVtx3_ak5PF_corr)[PFcount_corr]) TVector3(myJet.vx(), myJet.vy(), myJet.vz());
287 <        
288 <        ++PFcount_corr;
289 <      }
290 <    }
262 >      if(myJet.pt() > 5){
263 >
264 >        for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
265 >           if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
266 >        }
267 >
268 >        for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
269 >           if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
270 >        }
271 >                
272 >        TCJet* jetCon = new ((*recoJets)[pfCount]) TCJet;
273 >      
274 >        jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
275 >        jetCon->SetVtx(-999.0, -999.0, -999.0);
276 >        jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
277 >        jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
278 >        jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
279 >        jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
280 >        jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
281 >        jetCon->SetNumChPart(myJet.chargedMultiplicity());
282 >
283 >        if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
284 >        if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
285 >
286 >        double scale2 = correctorL2->correction(myJet);
287 >        myJet.scaleEnergy(scale2);
288 >        double scale3 = correctorL3->correction(myJet);
289 >        myJet.scaleEnergy(scale3);
290 >
291 >        //more corrections?
292 >
293 >        jetCon->SetJetCorr(2, scale2);
294 >        jetCon->SetJetCorr(3, scale3);
295 >
296 >        /////////////////////////
297 >        //get associated tracks//
298 >        /////////////////////////
299 >
300 >        const reco::TrackRefVector &tracks = myJet.getTrackRefs();
301 >
302 >        vector<TVector3> vtxPositionCollection;
303 >        vector<double>  associatedTrackSumPt;
304 >        vector<unsigned int> jetTrackAddresses;
305 >        double sumTrackX, sumTrackY, sumTrackZ, sumTrackIP, sumTrackPt;
306 >        int   nJetTracks, nVertexTracks, nAssociatedTracks, vertexIndex;
307 >        int   vCount = 0;
308 >
309 >        nJetTracks = nVertexTracks = nAssociatedTracks = 0;
310 >        sumTrackX = sumTrackY = sumTrackZ  = sumTrackIP  = sumTrackPt = 0;
311 >
312 >        if(myJet.pt() > 10 && fabs(myJet.eta()) < 2.4){
313 >
314 >          for(TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack){
315 >            const reco::Track &myJetTrack = **iTrack;
316 >
317 >            sumTrackPt += myJetTrack.pt();
318 >            sumTrackX  += myJetTrack.vx();
319 >            sumTrackY  += myJetTrack.vy();            
320 >            sumTrackZ  += myJetTrack.vz();
321 >            sumTrackIP += myJetTrack.dxy(vertexBeamSpot.position());
322 >            jetTrackAddresses.push_back((unsigned int)&myJetTrack);
323 >            ++nJetTracks;
324 >          }
325 >
326 >          if (nJetTracks > 0) {
327 >            jetCon->SetVtx(sumTrackX/nJetTracks, sumTrackY/nJetTracks, sumTrackZ/nJetTracks);          
328 >            h1_meanJetTrackZ->Fill(sumTrackZ/nJetTracks);
329 >            h1_trackDxy->Fill(sumTrackIP/nJetTracks);
330 >          }
331 >          if(jetTrackAddresses.size() > 0){
332 >
333 >            for (VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter) {      
334 >              reco::Vertex myVtx = reco::Vertex(*vtx_iter);
335 >              if(!myVtx.isValid() || myVtx.isFake()) continue;
336 >              TVector3 *iVtxPosition = new TVector3(myVtx.x(), myVtx.y(), myVtx.z());
337 >              vtxPositionCollection.push_back(*iVtxPosition);
338 >              associatedTrackSumPt.push_back(0);            
339 >            
340 >              for(Vertex::trackRef_iterator iTrackRef = myVtx.tracks_begin(); iTrackRef != myVtx.tracks_end(); ++iTrackRef){
341 >                const edm::RefToBase<reco::Track> &myTrackRef = *iTrackRef;
342 >                if(myTrackRef.isAvailable()){
343 >                  const reco::Track &myVertexTrack = *myTrackRef.get();        
344 >
345 >                  for(vector<unsigned int>::const_iterator iTrackAddress = jetTrackAddresses.begin(); iTrackAddress != jetTrackAddresses.end(); ++iTrackAddress){
346 >                    if (*iTrackAddress == (unsigned int)&myVertexTrack) {
347 >                      associatedTrackSumPt.at(vCount) += myVertexTrack.pt()/sumTrackPt;
348 >                      ++nAssociatedTracks;
349 >                    }
350 >                  }
351 >                }
352 >              }
353 >              ++vCount;  
354 >            }
355 >                    
356 >            double maxSumPtFraction = 0;
357 >            vCount = vertexIndex = 0;
358 >
359 >            for (vector<double>::const_iterator iTrackSumPt = associatedTrackSumPt.begin(); iTrackSumPt != associatedTrackSumPt.end(); ++iTrackSumPt) {
360 >              if (*iTrackSumPt > maxSumPtFraction) {
361 >                maxSumPtFraction = *iTrackSumPt;  
362 >                vertexIndex      = vCount + 1;
363 >              }
364 >              ++vCount;
365 >            }
366 >            if (vertexIndex > 0) {
367 >              h1_jetVertexZ->Fill(vtxPositionCollection[vertexIndex-1].z());
368 >            }
369 >            jetCon->SetVtxSumPt(maxSumPtFraction);
370 >            jetCon->SetVtxIndex(vertexIndex);
371 >
372 >            h1_fracAssociatedTracks->Fill((double)nAssociatedTracks/(double)nJetTracks);
373 >            h1_associatedSumPt->Fill(maxSumPtFraction);
374 >            h1_associatedVertexIndex->Fill(vertexIndex);
375 >            h2_nAssTracksVsJetPt->Fill(myJet.pt(), nAssociatedTracks);
376 >          }
377 >        }        
378 >        ++pfCount;
379 >      }      
380 >    }  
381    }
382 <  
383 <  if(doCaloJets_){
384 <    
385 <    Calocount = 0;
386 <    
387 <    Handle<reco::CaloJetCollection> CaloJets;
388 <    iEvent.getByLabel(CaloJetHandle_, CaloJets);
389 <    
390 <    for(CaloJetCollection::const_iterator jet_iter = CaloJets->begin(); jet_iter!= CaloJets->end(); ++jet_iter){
391 <      
392 <      reco::CaloJet myJet = reco::CaloJet(*jet_iter);
393 <      
394 <      if(myJet.et() > 5){
395 <        
396 <        new ((*jetP4_ak5Calo)[Calocount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
397 <        new ((*jetVtx3_ak5Calo)[Calocount]) TVector3(myJet.vx(), myJet.vy(), myJet.vz());
398 <        
399 <        ++Calocount;
400 <        
401 <      }
402 <    }
403 <    
404 <    Handle<reco::CaloJetCollection> CaloJets_corr;
405 <    iEvent.getByLabel(CaloJetHandle_corr, CaloJets_corr);
406 <    
407 <    Calocount_corr = 0;
408 <    
409 <    for(CaloJetCollection::const_iterator jet_iter = CaloJets_corr->begin(); jet_iter!= CaloJets_corr->end(); ++jet_iter){
410 <      
411 <      reco::CaloJet myJet = reco::CaloJet(*jet_iter);
412 <      
413 <      if(myJet.et() > 5){
414 <        
248 <        new ((*jetP4_ak5Calo_corr)[Calocount_corr]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
249 <        new ((*jetVtx3_ak5Calo_corr)[Calocount_corr]) TVector3(myJet.vx(), myJet.vy(), myJet.vz());
250 <        
251 <        ++Calocount_corr;
252 <      }
382 >
383 >
384 >  ///////////////////
385 >  // Get muon info //
386 >  ///////////////////
387 >
388 >
389 >  Handle<MuonCollection> muons;
390 >  iEvent.getByLabel(muonTag_,muons);
391 >
392 >  int muCount = 0;
393 >  for (MuonCollection::const_iterator mu = muons->begin(); mu != muons->end(); ++mu) {
394 >    if (mu->isGlobalMuon() && mu->pt() > 8){
395 >      TCMuon* lepCon = new ((*recoMuons)[muCount]) TCMuon;
396 >      lepCon->Setp4(mu->px(), mu->py(), mu->pz(), mu->p());
397 >      lepCon->SetVtx(mu->globalTrack()->vx(),mu->globalTrack()->vy(),mu->globalTrack()->vz());
398 >      lepCon->SetEta(mu->eta());
399 >      lepCon->SetPhi(mu->phi());
400 >      lepCon->SetCharge(mu->charge());
401 >      lepCon->SetisGLB(mu->isGlobalMuon());
402 >      lepCon->SetisTRK(mu->isTrackerMuon());
403 >      lepCon->Setdxy(mu->globalTrack()->dxy(vertexBeamSpot.position()));
404 >      lepCon->SetnPXLHits(mu->globalTrack()->hitPattern().numberOfValidPixelHits());
405 >      lepCon->SetnTRKHits(mu->globalTrack()->hitPattern().numberOfValidTrackerHits());
406 >      lepCon->SetnValidMuHits(mu->globalTrack()->hitPattern().numberOfValidMuonHits());
407 >      lepCon->SetnMatchSeg(mu->numberOfMatches());
408 >      lepCon->SetNormChi2(mu->globalTrack()->normalizedChi2());
409 >      lepCon->SetCaloComp(mu->caloCompatibility());
410 >      lepCon->SetSegComp(muon::segmentCompatibility(*mu));
411 >      lepCon->SetEMIso(mu->isolationR03().emEt);
412 >      lepCon->SetHADIso(mu->isolationR03().hadEt);
413 >      lepCon->SetTRKIso(mu->isolationR03().sumPt);
414 >      muCount++;
415      }
416    }
417 +
418    
419 <  if(doGenJets_){
420 <    
421 <    Handle<reco::GenJetCollection> GenJets;
422 <    iEvent.getByLabel(GenJetHandle_, GenJets);
423 <    
424 <    
425 <    for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
426 <      
427 <      reco::GenJet myJet = reco::GenJet(*jet_iter);
428 <      
429 <      if(myJet.et() > 5){
430 <        
431 <        new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
432 <        new ((*jetVtx3_ak5Gen)[Gencount]) TVector3(myJet.vx(), myJet.vy(), myJet.vz());
433 <        
434 <        ++Gencount;
435 <        
436 <      }
419 >  ///////////////////////
420 >  // Get electron info //
421 >  ///////////////////////
422 >
423 >
424 >  Handle<edm::ValueMap<float> > eIDValueMap;
425 >  iEvent.getByLabel( electronIDMap_ , eIDValueMap );
426 >  const edm::ValueMap<float> & eIDmap = * eIDValueMap ;
427 >
428 >  Handle<GsfElectronCollection> electrons;
429 >  iEvent.getByLabel(electronTag_,electrons);
430 >
431 >  int elCount = 0;
432 >  for (unsigned int i = 0; i < electrons->size(); i++){
433 >    edm::Ref<reco::GsfElectronCollection> electronRef(electrons,i);
434 >    if ( eIDmap[electronRef] == 7 && electronRef->pt() > 15){
435 >      TCElectron* lepCon = new ((*recoElectrons)[elCount]) TCElectron;
436 >      lepCon->Setp4(electronRef->px(),electronRef->py(),electronRef->pz(),electronRef->p());
437 >      lepCon->SetVtx(electronRef->gsfTrack()->vx(),electronRef->gsfTrack()->vy(),electronRef->gsfTrack()->vz());
438 >      lepCon->SetCharge(electronRef->charge());
439 >      lepCon->SetEta(electronRef->eta());
440 >      lepCon->SetPhi(electronRef->phi());
441 >      lepCon->Setdxy(electronRef->gsfTrack()->dxy(vertexBeamSpot.position()));
442 >      lepCon->SetNormChi2(electronRef->gsfTrack()->normalizedChi2());
443 >      lepCon->SetEMIso(electronRef->dr03EcalRecHitSumEt());
444 >      lepCon->SetHADIso(electronRef->dr03HcalTowerSumEt());
445 >      lepCon->SetTRKIso(electronRef->dr03TkSumPt());
446 >      elCount++;
447      }
448    }
449 <  
277 <  
278 <  Handle<reco::VertexCollection> primaryVtx;
279 <  iEvent.getByLabel(PrimaryVtxHandle_, primaryVtx);
449 >      
450  
451 <  for(VertexCollection::const_iterator vtx_iter = primaryVtx->begin(); vtx_iter!= primaryVtx->end(); ++vtx_iter){
452 <    
453 <    reco::Vertex myVtx = reco::Vertex(*vtx_iter);
284 <    
285 <    new ((*PrimaryVtx3)[Vtxcount]) TVector3(myVtx.x(), myVtx.y(), myVtx.z());
286 <    
287 <    ++Vtxcount;
288 <    
289 <  }
290 <  
291 <  
292 <  if(Calocount > 0 || Calocount_corr >0 || PFcount_corr > 0 || PFcount > 0 || Gencount > 0)  sTree -> Fill();
451 >  ////////////////////////
452 >  // Get gen-level info //
453 >  ////////////////////////
454  
294  jetP4_ak5PF->Clear();
295  jetP4_ak5PF_corr->Clear();
296  jetVtx3_ak5PF->Clear();
297  jetP4_ak5Calo->Clear();
298  jetP4_ak5Calo_corr->Clear();
299  jetVtx3_ak5Calo->Clear();
300  jetP4_ak5Gen->Clear();
301  jetVtx3_ak5Gen->Clear();
302  
455  
456 +  if (!isRealData) {
457 +    
458 +    Handle<GenEventInfoProduct> GenEventInfoHandle;
459 +    iEvent.getByLabel("generator", GenEventInfoHandle);
460 +    Handle<GenRunInfoProduct> GenRunInfoHandle;
461 +    iEvent.getByLabel("generator", GenRunInfoHandle);
462 +    Handle<reco::GenJetCollection> GenJets;
463 +    iEvent.getByLabel(genJetTag_, GenJets);
464  
465 <  //----------  Filling HLT trigger bits!  ------------
465 >    ptHat = qScale = -1; crossSection = 0;
466  
467 <  edm::Handle<TriggerResults> hltR;
468 <  triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
469 <  iEvent.getByLabel(triggerResultsTag_,hltR);
470 <  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
471 <  hlNames=triggerNames.triggerNames();
472 <  bool triggerPassed = false;
467 >    if (GenEventInfoHandle.isValid()) {
468 >      qScale       = GenEventInfoHandle->qScale();
469 >      ptHat        = (GenEventInfoHandle->hasBinningValues() ? GenEventInfoHandle->binningValues()[0] : 0.0);
470 >    }
471 >   if (GenRunInfoHandle.isValid()) {
472 >      crossSection = GenRunInfoHandle->crossSection();
473 >    }
474 >        
475 >    for (GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter) {
476 >      reco::GenJet myJet = reco::GenJet(*jet_iter);      
477 >      if (myJet.pt() > 5) {    
478 >
479 >        TCGenJet* jetCon = new ((*genJets)[genCount]) TCGenJet;
480 >        jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
481 >        jetCon->SetHadEnergy(myJet.hadEnergy());
482 >        jetCon->SetEmEnergy(myJet.emEnergy());
483 >        jetCon->SetInvEnergy(myJet.invisibleEnergy());
484 >        jetCon->SetAuxEnergy(myJet.auxiliaryEnergy());
485 >        jetCon->SetNumConstit(myJet.getGenConstituents().size());
486 >
487 >        /*const reco::GenParticle *myGenParticle = myJet.getGenConstituent(0);
488 >
489 >        if(myGenParticle->numberOfMothers() == 1)
490 >        {
491 >          const reco::Candidate *myMother = myGenParticle->mother();
492 >          int   nGrandMas = myMother->numberOfMothers();
493 >        
494 >          for(int iter = 0; nGrandMas != 0 ; ++iter)
495 >          {
496 >            myMother  = myMother->mother();
497 >            nGrandMas = myMother->mother()->numberOfMothers();            
498 >            //cout<<myMother->pdgId()<<endl;
499 >          }
500 >          //if(nGrandMas == 0) jetCon->SetJetFlavor((int)myMother->pdgId()); else jetCon->SetJetFlavor(0);
501 >        }*/
502 >        ++genCount;    
503 >      }
504 >    }    
505 >  }
506    
507 <  for (uint iT=0; iT<hlNames.size(); ++iT)
508 <    {
509 <      triggerPassed = triggerDecision(hltR, iT);
507 >  ///////////////////////////  
508 >  //get trigger information//
509 >  ///////////////////////////
510 >
511 >  if (triggerHLT_) {
512 >
513 >    edm::Handle<TriggerResults> hltR;
514 >    triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltProcess_);
515 >    iEvent.getByLabel(triggerResultsTag_,hltR);
516  
517 < if (hlNames[iT]=="HLT_PixelTracks_Multiplicity70")
518 < {
320 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
321 < //And fill the bit here
322 < }
517 >    const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
518 >    hlNames=triggerNames.triggerNames();  
519  
520 < if (hlNames[iT]=="HLT_MinBiasBSC_NoBPTX")
521 < {
522 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
523 < //And fill the bit here
524 < }
520 >    bool triggerPassed = false;
521 >    triggerStatus = 0x0;    
522 >    string MPI_TriggerNames[] = {"HLT_QuadJet15U",
523 >                                 "HLT_DiJetAve50U",
524 >                                 "HLT_DiJetAve30U",
525 >                                 "HLT_DiJetAve15U",
526 >                                 "HLT_FwdJet20U",
527 >                                 "HLT_Jet100U",
528 >                                 "HLT_Jet70U",
529 >                                 "HLT_Jet50U",
530 >                                 "HLT_Jet30U",
531 >                                 "HLT_Jet15U",
532 >                                 "HLT_BTagMu_Jet10U",
533 >                                 "HLT_DoubleJet15U_ForwardBackward",
534 >                                 "HLT_BTagIP_Jet50U",
535 >                                 "HLT_DoubleLooseIsoTau15",
536 >                                 "HLT_SingleLooseIsoTau20",
537 >                                 "HLT_HT100U",
538 >                                 "HLT_MET100",
539 >                                 "HLT_MET45"};
540 >
541 >    for (int i=0; i < (int)hlNames.size(); ++i) {      
542 >      triggerPassed = triggerDecision(hltR, i);      
543 >      if(triggerPassed){        
544 >        for (int j = 0; j < 17; ++j){
545 >          if (hlNames[i] == MPI_TriggerNames[j]) {
546 >            triggerStatus |= 0x01 << j;      
547 >          }
548 >        }
549 >      }
550 >    }
551 >  }
552 >
553 >  if((pfCount > 1 || genCount > 1) && vtxCount > 0) sTree -> Fill();
554  
555 < if (hlNames[iT]=="HLT_PixelTracks_Multiplicity40")
556 < {
557 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
558 < }
559 < if (hlNames[iT]=="HLT_L1Tech_HCAL_HF")
335 < {
336 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
337 < }
338 < if (hlNames[iT]=="HLT_IsoTrackHB_8E29")
339 < {
340 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
341 < }
342 < if (hlNames[iT]=="HLT_IsoTrackHE_8E29")
343 < {
344 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
345 < }
346 < if (hlNames[iT]=="HLT_L1Tech_RPC_TTU_RBst1_collisions")
347 < {
348 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
349 < }
350 < if (hlNames[iT]=="HLT_L1_BscMinBiasOR_BptxPlusORMinus")
351 < {
352 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
353 < }
354 < if (hlNames[iT]=="HLT_L1Tech_BSC_halo_forPhysicsBackground")
355 < {
356 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
357 < }
358 < if (hlNames[iT]=="HLT_L1Tech_BSC_HighMultiplicity")
359 < {
360 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
361 < }
362 < if (hlNames[iT]=="HLT_MinBiasPixel_DoubleIsoTrack5")
363 < {
364 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
365 < }
366 < if (hlNames[iT]=="HLT_MinBiasPixel_DoubleTrack")
367 < {
368 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
369 < }
370 < if (hlNames[iT]=="HLT_MinBiasPixel_SingleTrack")
371 < {
372 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
373 < }
374 < if (hlNames[iT]=="HLT_ZeroBiasPixel_SingleTrack")
375 < {
376 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
377 < }
378 < if (hlNames[iT]=="HLT_MinBiasBSC")
379 < {
380 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
381 < }
382 < if (hlNames[iT]=="HLT_StoppedHSCP_8E29")
383 < {
384 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
385 < }
386 < if (hlNames[iT]=="HLT_Jet15U_HcalNoiseFiltered")
387 < {
388 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
389 < }
390 < if (hlNames[iT]=="HLT_QuadJet15U")
391 < {
392 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
393 < }
394 < if (hlNames[iT]=="HLT_DiJetAve30U_8E29")
395 < {
396 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
397 < }
398 < if (hlNames[iT]=="HLT_DiJetAve15U_8E29")
399 < {
400 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
401 < }
402 < if (hlNames[iT]=="HLT_FwdJet20U")
403 < {
404 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
555 >  recoJets->Clear();
556 >  recoElectrons->Clear();
557 >  recoMuons->Clear();
558 >  primaryVtx->Clear();
559 >  genJets->Clear();
560   }
561 < if (hlNames[iT]=="HLT_Jet50U")
562 < {
563 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
564 < }
565 < if (hlNames[iT]=="HLT_Jet30U")
566 < {
567 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
568 < }
569 < if (hlNames[iT]=="HLT_Jet15U")
570 < {
571 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
572 < }
573 < if (hlNames[iT]=="HLT_BTagMu_Jet10U")
574 < {
575 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
576 < }
577 < if (hlNames[iT]=="HLT_DoubleJet15U_ForwardBackward")
578 < {
579 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
580 < }
581 < if (hlNames[iT]=="HLT_BTagIP_Jet50U")
582 < {
583 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
584 < }
585 < if (hlNames[iT]=="HLT_DoubleLooseIsoTau15")
586 < {
587 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
588 < }
589 < if (hlNames[iT]=="HLT_SingleLooseIsoTau20")
590 < {
591 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
592 < }
593 < if (hlNames[iT]=="HLT_HT100U")
594 < {
595 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
596 < }
442 < if (hlNames[iT]=="HLT_MET100")
443 < {
444 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
561 >
562 > // ------------ method called once each job just before starting event loop  ------------
563 > void  MPIntuple::beginJob()
564 > {  
565 >  ntupleFile               = new TFile(rootfilename.c_str(), "RECREATE");
566 >  sTree                    = new TTree("mpiTree", "Tree for Jets");
567 >
568 >  h1_fracAssociatedTracks  = new TH1D("h1_fracAssociatedTracks", "Fraction of associated jet tracks", 20, 0.0, 1.0);
569 >  h1_trackDxy              = new TH1D("h1_trackDxy", "dxy for all associated tracks", 50, -0.4, 0.4);
570 >  h1_meanJetTrackZ         = new TH1D("h1_meanJetTrackZ", "<z> from jet tracks", 50, -25.0, 25.0);
571 >  h1_jetVertexZ            = new TH1D("h1_jetVertexZ", "z for vertex associated to jet", 50, -25.0, 25.0);
572 >  h1_associatedSumPt       = new TH1D("h1_associatedSumPt", "ratio of sum p_{T} between associatedTracks and jetTracks", 20, 0.0, 1.0);
573 >  h1_associatedVertexIndex = new TH1D("h1_associatedVertex", "number of jets associated to primary vertex", 5, -0.5, 4.5);
574 >  h2_nAssTracksVsJetPt     = new TH2F("h2_nAssTracksVsJetPt", "nTracks vs Jet p_{T}", 50, 0, 100, 20, -0.5, 19.5);
575 >  p1_nVtcs                 = new TProfile("p1_nVtcs", "<N> per run", 5200.0, 132300.0, 137500.0, 0.0, 6.0);
576 >
577 >  recoJets                 = new TClonesArray("TCJet");
578 >  recoElectrons            = new TClonesArray("TCElectron");
579 >  recoMuons                = new TClonesArray("TCMuon");
580 >  genJets                  = new TClonesArray("TCGenJet");
581 >  primaryVtx               = new TClonesArray("TCPrimaryVtx");
582 >
583 >  sTree->Branch("recoJets",&recoJets, 6400, 0);
584 >  sTree->Branch("recoElectrons",&recoElectrons, 6400, 0);
585 >  sTree->Branch("recoMuons",&recoMuons, 6400, 0);
586 >  sTree->Branch("genJets",&genJets, 6400, 0);
587 >  sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
588 >  sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
589 >  sTree->Branch("runNumber",&runNumber, "runNumber/I");
590 >  sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
591 >  sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
592 >  sTree->Branch("isRealData",&isRealData, "isRealData/i");
593 >  sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
594 >  sTree->Branch("ptHat",&ptHat, "ptHat/d");
595 >  sTree->Branch("qScale", &qScale, "qScale/d");
596 >  sTree->Branch("crossSection", &crossSection, "crossSection/d");
597   }
598 < if (hlNames[iT]=="HLT_MET45")
598 >
599 > void MPIntuple::beginRun(const edm::Run& iRun, const edm::EventSetup& iEvent)
600   {
601 <      cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
601 >  bool changed = true;
602 >  hltConfig_.init(iRun, iEvent, hltProcess_, changed);
603 >
604   }
450        
451      
452      // ----------->> Fill them Here! <<----------------------
453    }
454  //----------         ------------------   ---------------  
455
605  
606 + void MPIntuple::endLuminosityBlock(const edm::LuminosityBlock& iLumi, const edm::EventSetup& iEvent)
607 + {
608 +  edm::Handle<LumiSummary> lumiSummary;
609 +  iLumi.getByLabel("lumiProducer", lumiSummary);
610  
611 +  
612   }
613  
614 <
461 < // ------------ method called once each job just before starting event loop  ------------
462 < void  MPIntuple::beginJob()
614 > void MPIntuple::endRun(const edm::Run& iRun, const edm::EventSetup& iEvent)
615   {
464  
465  ntupleFile           = new TFile(rootfilename.c_str(), "RECREATE");
466  sTree                = new TTree("dpsTree", "Tree for Jets");
467  
468  jetP4_ak5PF          = new TClonesArray("TLorentzVector");
469  jetVtx3_ak5PF        = new TClonesArray("TVector3");
470  jetP4_ak5Calo        = new TClonesArray("TLorentzVector");
471  jetVtx3_ak5Calo      = new TClonesArray("TVector3");
472  jetP4_ak5PF_corr     = new TClonesArray("TLorentzVector");
473  jetVtx3_ak5PF_corr   = new TClonesArray("TVector3");
474  jetP4_ak5Calo_corr   = new TClonesArray("TLorentzVector");
475  jetVtx3_ak5Calo_corr = new TClonesArray("TVector3");
476  jetP4_ak5Gen         = new TClonesArray("TLorentzVector");
477  jetVtx3_ak5Gen       = new TClonesArray("TVector3");
478  PrimaryVtx3          = new TClonesArray("TVector3");
479
480  sTree->Branch("jetP4_ak5PF",&jetP4_ak5PF, 6400, 0);
481  sTree->Branch("jetP4_ak5PF_corr",&jetP4_ak5PF_corr, 6400, 0);
482  sTree->Branch("jetVtx3_ak5PF",&jetVtx3_ak5PF, 6400, 0);
483  sTree->Branch("jetVtx3_ak5PF_corr",&jetVtx3_ak5PF, 6400, 0);
484  sTree->Branch("jetP4_ak5Calo",&jetP4_ak5Calo, 6400, 0);
485  sTree->Branch("jetP4_ak5Calo_corr",&jetP4_ak5Calo_corr, 6400, 0);
486  sTree->Branch("jetVtx3_ak5Calo",&jetVtx3_ak5Calo, 6400, 0);
487  sTree->Branch("jetVtx3_ak5Calo_corr",&jetVtx3_ak5Calo, 6400, 0);
488  sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0);
489  sTree->Branch("jetVtx3_ak5Gen",&jetVtx3_ak5Gen, 6400, 0);
490  sTree->Branch("PrimaryVtx",&PrimaryVtx3, 6400, 0);
491
492  sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
493  sTree->Branch("runNumber",&runNumber, "runNumber/I");
494  sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
616  
617   }
497
618   // ------------ method called once each job just after ending the event loop  ------------
619   void MPIntuple::endJob()
620   {
621 +
622 +  ntupleFile->cd();
623 +
624 +  h1_fracAssociatedTracks->Write();
625 +  h1_meanJetTrackZ->Write();
626 +  h1_trackDxy->Write();
627 +  h1_jetVertexZ->Write();
628 +  h1_associatedSumPt->Write();
629 +  h1_associatedVertexIndex->Write();
630 +  h2_nAssTracksVsJetPt->Write();
631 +  p1_nVtcs->Write();
632 +
633    ntupleFile->Write();
634    ntupleFile->Close();
635  
# Line 513 | Line 645 | bool MPIntuple::triggerDecision(edm::Han
645        triggerPassed = true;
646      }
647      return triggerPassed;
648 <  }
648 > }
649 >
650 > double MPIntuple::sumPtSquared(const Vertex & v)
651 > {
652 >  double sum = 0.;
653 >  double pT;
654 >  for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
655 >    pT = (**it).pt();
656 >    double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
657  
658 +    sum += pT*pT;
659 +  }
660 +  return sum;
661 + }
662  
663   //define this as a plug-in
664   DEFINE_FWK_MODULE(MPIntuple);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines