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.8 by naodell, Fri Jun 11 23:01:47 2010 UTC vs.
Revision 1.21 by naodell, Wed Nov 3 09:53:54 2010 UTC

# Line 18 | Line 18
18   #include "FWCore/ParameterSet/interface/ParameterSet.h"
19   #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
20   #include "Geometry/Records/interface/CaloGeometryRecord.h"
21 + #include "DataFormats/Math/interface/deltaPhi.h"
22  
23 < // Jet and vertex functions
23 > // Libraries for objects
24   #include "DataFormats/JetReco/interface/CaloJet.h"
25   #include "DataFormats/JetReco/interface/CaloJetCollection.h"
26   #include "DataFormats/JetReco/interface/PFJet.h"
# Line 30 | Line 31
31   #include "DataFormats/VertexReco/interface/Vertex.h"
32   #include "DataFormats/VertexReco/interface/VertexFwd.h"
33   #include "DataFormats/BTauReco/interface/JetTag.h"
34 + #include "DataFormats/TrackReco/interface/Track.h"
35 + #include "DataFormats/BeamSpot/interface/BeamSpot.h"
36 + #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
37  
38 < #include "DataFormats/Math/interface/deltaPhi.h"
38 > //#include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h"
39 >
40 > //GenParticles
41 > //#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
42  
43   // Need to get correctors
44   #include "JetMETCorrections/Objects/interface/JetCorrector.h"
# Line 47 | Line 54
54  
55   //Root  stuff
56   #include "TROOT.h"
57 + #include "TH1.h"
58 + #include "TH2.h"
59 + #include "TProfile.h"
60   #include "TFile.h"
61   #include "TTree.h"
62   #include "TString.h"
# Line 77 | Line 87 | class MPIntuple : public edm::EDAnalyzer
87        virtual void endJob() ;
88  
89    bool triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger);
90 +  double sumPtSquared(const Vertex & v);
91  
92        // ----------member data ---------------------------
93  
# Line 84 | Line 95 | class MPIntuple : public edm::EDAnalyzer
95    edm::InputTag GenJetHandle_;
96    edm::InputTag PrimaryVtxHandle_;
97    edm::InputTag triggerResultsTag_;
87
88
98    //  Counting variables   //
99  
100    int eventNumber, runNumber, lumiSection, bunchCross;
101 <
101 >  double ptHat, qScale;
102 >  
103    TTree* sTree;
104    TFile* ntupleFile;
105  
# Line 110 | Line 120 | class MPIntuple : public edm::EDAnalyzer
120    vector<string>  hlNames;
121    unsigned int triggerStatus;
122  
123 +  //Histograms
124 +  TH1F * h1_nAssociatedTracks;
125 +  TH1F * h1_trackZ;
126 +  TH1F * h1_trackDxy;
127 +  TH1F * h1_allTrackDeltaZ_PV;
128 +  TH1F * h1_allTrackDeltaZ_LJ;
129 +  TH2F * h2_nAssTracksVsJetPt;
130 +  TProfile * p1_nVtcs;
131  
114 };
132  
133 < MPIntuple::MPIntuple(const edm::ParameterSet& iConfig):
133 > };
134  
135 <  PFJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag")),
119 <  GenJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag")),
120 <  PrimaryVtxHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")),
121 <  doGenJets_(iConfig.getUntrackedParameter<bool>("doGenJets")),
122 <  doPFJets_(iConfig.getUntrackedParameter<bool>("doPFJets")),
123 <  triggerHLT_(iConfig.getUntrackedParameter<bool>("triggerHLT")),
124 <  rootfilename(iConfig.getUntrackedParameter<string>("rootfilename")),
125 <  hlTriggerResults_(iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")),
126 <  hltName_(iConfig.getUntrackedParameter<string>("hltName","HLT"))
135 > MPIntuple::MPIntuple(const edm::ParameterSet& iConfig)
136   {
137 <  //edm::TriggerNames
138 <  // triggerNames(iConfig);
139 <
137 >  PFJetHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag");
138 >  GenJetHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag");
139 >  PrimaryVtxHandle_ = iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag");
140 >  doGenJets_ = iConfig.getUntrackedParameter<bool>("doGenJets");
141 >  doPFJets_ = iConfig.getUntrackedParameter<bool>("doPFJets");
142 >  triggerHLT_ = iConfig.getUntrackedParameter<bool>("triggerHLT");
143 >  rootfilename = iConfig.getUntrackedParameter<string>("rootfilename");
144 >  hlTriggerResults_ = iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults");
145 >  hltName_ = iConfig.getUntrackedParameter<string>("hltName","HLT");
146   }
147  
133
148   MPIntuple::~MPIntuple()
149   {
150  
# Line 150 | Line 164 | void MPIntuple::analyze(const edm::Event
164    bunchCross   = (unsigned int)iEvent.bunchCrossing();
165    isRealData   = iEvent.isRealData();
166  
167 <  int PFcount  = 0;
168 <  int Gencount = 0;
169 <  int Vtxcount = 0;
167 >
168 >  edm::Handle<reco::BeamSpot> beamSpotHandle;
169 >  iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
170 >  reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
171 >
172 >  int pfCount   = 0;
173 >  int genCount  = 0;
174 >  int vtxCount  = 0;
175 >  float ljVertexZ = -999;
176 >  float primaryVertexZ = -999;
177 >
178 >  //////////////////////////
179 >  //Get vertex information//
180 >  //////////////////////////
181 >
182 >  Handle<reco::VertexCollection> primaryVtcs;
183 >  iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
184 >  
185 >  for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
186 >    
187 >    reco::Vertex myVtx = reco::Vertex(*vtx_iter);
188 >    if(!myVtx.isValid() || myVtx.isFake()) continue;
189 >    TCPrimaryVtx* vtxCon = new ((*primaryVtx)[vtxCount]) TCPrimaryVtx;
190 >      
191 >    cout<<myVtx.nTracks()<<endl;
192 >
193 >    vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
194 >    vtxCon->SetNDof(myVtx.ndof());
195 >    vtxCon->SetChi2(myVtx.chi2());
196 >    vtxCon->SetNTrks(myVtx.tracksSize());
197 >    vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
198 >
199 >    if(vtxCount == 0) primaryVertexZ = myVtx.z();
200 >
201 >    ++vtxCount;
202 >    
203 >  }
204 >
205 >  p1_nVtcs->Fill(runNumber, vtxCount);
206 >
207 >  ///////////////////////
208 >  //get jet information//
209 >  ///////////////////////
210  
211    if(doPFJets_){
212  
# Line 166 | Line 220 | void MPIntuple::analyze(const edm::Event
220      const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
221      reco::JetTagCollection::const_iterator jet_it_2;
222  
223 +
224      const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
225      const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
226  
227      Handle<reco::PFJetCollection> PFJets;
228      iEvent.getByLabel(PFJetHandle_, PFJets);
229 <    
175 <    PFcount = 0;
176 <    
229 >        
230      for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
231        
232        reco::PFJet myJet = reco::PFJet(*jet_iter);
180
181 //      float scale2 = correctorL2->correction(myJet);
182 //      float scale3 = correctorL3->correction(myJet);
183
233        if(myJet.et() > 5){
234  
235          for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
# Line 191 | Line 240 | void MPIntuple::analyze(const edm::Event
240             if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
241          }
242                  
243 <        TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
243 >        TCJet* jetCon = new ((*jet_ak5PF)[pfCount]) TCJet;
244        
245          jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
246 <        jetCon->SetVtx(myJet.vx(), myJet.vy(), myJet.vz());
246 >        jetCon->SetVtx(-999.0, -999.0, -999.0);
247  
248          jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
249          jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
# Line 219 | Line 268 | void MPIntuple::analyze(const edm::Event
268          jetCon->SetJetCorr(2, scale2);
269          jetCon->SetJetCorr(3, scale3);
270  
271 < //if(( (lumiSection==666690) && (eventNumber==2439)) ||
272 < //   ( (lumiSection==666702) && (eventNumber==3629)) ||
273 < //   ( (lumiSection==666702) && (eventNumber==3681)) ||
274 < //   ( (lumiSection==666709) && (eventNumber==4320)) ||
275 < //   ( (lumiSection==666709) && (eventNumber==4328)) ||
276 < //   ( (lumiSection==666724) && (eventNumber==5855)) ||
277 < //   ( (lumiSection==666742) && (eventNumber==7642)) ||
278 < //   ( (lumiSection==666754) && (eventNumber==8869)) ||
279 < //   ( (lumiSection==666670) && (eventNumber==466)) ||
280 < //   ( (lumiSection==666675) && (eventNumber==919)) ||
281 < //   ( (lumiSection==666718) && (eventNumber==5235)))
282 < //{
234 < //      std::cout << runNumber << " " << lumiSection << "  " << eventNumber << "\n";
235 < //      std::cout << "(" << myJet.pt() << "," << myJet.eta() << "," << myJet.phi() << ", [" << scale2 << ", " << scale3 << "]) " << std::endl;
236 < //      std::cout << "(" << jetCon->Pt(3) << "," << jetCon->P4(3).Eta() << "," << jetCon->P4(3).Phi() << "," << ", [" << jetCon->JetCorr(2) << ", " <<  jetCon->JetCorr(3)  << "]) " << std::endl;
237 < //}
271 >        /////////////////////////
272 >        //get associated tracks//
273 >        /////////////////////////
274 >
275 >        const reco::TrackRefVector &tracks = myJet.getTrackRefs();
276 >
277 >        float sumTrackZ, sumTrackIP;
278 >        float meanTrackZ, meanTrackIP;
279 >        int   nAssociatedTracks = -1;
280 >
281 >        sumTrackZ  = sumTrackIP  = 0;
282 >        meanTrackZ = meanTrackIP = -999;
283  
284 +        if(myJet.pt() > 20 && fabs(myJet.eta()) < 2.4){
285  
286 <        ++PFcount;
286 >          nAssociatedTracks = 0;
287 >
288 >          for(TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack){
289 >
290 >            const reco::Track    &myJetTrack = **iTrack;
291 >            //const reco::TrackRef &myTrackRef = *iTrack;
292 >
293 >            sumTrackZ += myJetTrack.vz();
294 >            sumTrackIP += myJetTrack.dxy(vertexBeamSpot.position());
295 >            ++nAssociatedTracks;
296 >
297 >            for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
298 >              
299 >              reco::Vertex myVtx = reco::Vertex(*vtx_iter);
300 >              if(!myVtx.isValid() || myVtx.isFake()) continue;
301 >              
302 >              for(Vertex::trackRef_iterator iTrackRef = myVtx.tracks_begin(); iTrackRef != myVtx.tracks_end(); ++iTrackRef){
303 >                const edm::RefToBase<reco::Track> &myTrackRef = *iTrackRef;
304 >
305 >                if(myTrackRef.isAvailable()){
306 >
307 >                  const reco::Track &myVertexTrack = *myTrackRef.get();
308 >                                
309 >                  cout<<"Vertex track z = "<<myVertexTrack.vz()<<endl;
310 >                  cout<<"Jet track z = "<<myJetTrack.vz()<<endl;
311 >                
312 >                }
313 >              }      
314 >            }              
315 >          }
316 >          
317 >          meanTrackZ = sumTrackZ/nAssociatedTracks;
318 >          h1_nAssociatedTracks->Fill(nAssociatedTracks);
319 >          h2_nAssTracksVsJetPt->Fill(myJet.pt(), nAssociatedTracks);
320 >          
321 >          for(TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); ++iTrack){
322 >            
323 >            const reco::Track &myJetTrack = **iTrack;
324 >            
325 >            if(nAssociatedTracks == 3) h1_trackZ->Fill(meanTrackZ - myJetTrack.vz());
326 >            h1_trackDxy->Fill(myJetTrack.dxy(vertexBeamSpot.position()));
327 >            
328 >          }    
329 >          
330 >          if(pfCount == 0){
331 >            
332 >            ljVertexZ = meanTrackZ;
333 >            
334 >          }else{
335 >            
336 >            h1_allTrackDeltaZ_LJ->Fill(ljVertexZ - meanTrackZ);
337 >            
338 >          }
339 >          
340 >          h1_allTrackDeltaZ_PV->Fill(primaryVertexZ - meanTrackZ);
341 >          
342 >        }
343 >        
344 >        jetCon->SetVtx(0, 0, meanTrackZ);      
345 >        
346 >        ++pfCount;
347        }      
348      }  
349    }
350    
351 <  if(doGenJets_){
351 >  if(!isRealData){
352      
353 +    Handle<GenEventInfoProduct> GenInfoHandle;
354 +    iEvent.getByLabel("generator", GenInfoHandle);
355 +
356      Handle<reco::GenJetCollection> GenJets;
357      iEvent.getByLabel(GenJetHandle_, GenJets);
358 <    
359 <    
358 >
359 >    ptHat = qScale = -1;
360 >
361 >    if(GenInfoHandle.isValid()){
362 >
363 >      qScale = GenInfoHandle->qScale();
364 >      ptHat  = (GenInfoHandle->hasBinningValues() ? GenInfoHandle->binningValues()[0] : 0.0);
365 >
366 >    }
367 >        
368      for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
369        
370        reco::GenJet myJet = reco::GenJet(*jet_iter);
371        
372        if(myJet.pt() > 5){
373          
374 <        new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
374 >        new ((*jetP4_ak5Gen)[genCount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
375          
376 <        ++Gencount;
376 >        ++genCount;
377          
378        }
379 <    }
263 <  }
264 <
265 <
266 <  Handle<reco::VertexCollection> primaryVtcs;
267 <  iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
268 <
269 <  for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
270 <    
271 <    reco::Vertex myVtx = reco::Vertex(*vtx_iter);
272 <    
273 <    TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
274 <      
275 <    vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
276 <    vtxCon->SetNDof(myVtx.ndof());
277 <    vtxCon->SetChi2(myVtx.chi2());
278 <            
279 <    ++Vtxcount;
280 <    
379 >    }    
380    }
381 <
381 >  
382 >  ///////////////////////////  
383 >  //get trigger information//
384 >  ///////////////////////////
385  
386    if(triggerHLT_){
387  
286    //----------  Filling HLT trigger bits!  ------------
287
388      edm::Handle<TriggerResults> hltR;
389      triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
390      iEvent.getByLabel(triggerResultsTag_,hltR);
# Line 292 | Line 392 | void MPIntuple::analyze(const edm::Event
392      const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
393      hlNames=triggerNames.triggerNames();
394      
395 <    string MPI_TriggerNames[32] = {"HLT_PixelTracks_Multiplicity70", "HLT_MinBiasBSC_NoBPTX", "HLT_PixelTracks_Multiplicity40","HLT_L1Tech_HCAL_HF", "HLT_IsoTrackHB_8E29", "HLT_IsoTrackHE_8E29", "HLT_L1Tech_RPC_TTU_RBst1_collisions", "HLT_L1_BscMinBiasOR_BptxPlusORMinus", "HLT_L1Tech_BSC_halo_forPhysicsBackground", "HLT_L1Tech_BSC_HighMultiplicity", "HLT_MinBiasPixel_DoubleIsoTrack5", "HLT_MinBiasPixel_DoubleTrack", "HLT_MinBiasPixel_SingleTrack", "HLT_ZeroBiasPixel_SingleTrack", "HLT_MinBiasBSC", "HLT_StoppedHSCP_8E29", "HLT_Jet15U_HcalNoiseFiltered", "HLT_QuadJet15U", "HLT_DiJetAve30U_8E29", "HLT_DiJetAve15U_8E29", "HLT_FwdJet20U", "HLT_Jet50U", "HLT_Jet30U", "HLT_Jet15U", "HLT_BTagMu_Jet10U", "HLT_DoubleJet15U_ForwardBackward", "HLT_BTagIP_Jet50U", "HLT_DoubleLooseIsoTau15", "HLT_SingleLooseIsoTau20", "HLT_HT100U", "HLT_MET100", "HLT_MET45"};
395 >    string MPI_TriggerNames[] = {"HLT_PixelTracks_Multiplicity70", "HLT_MinBiasBSC_NoBPTX", "HLT_PixelTracks_Multiplicity40","HLT_L1Tech_HCAL_HF", "HLT_IsoTrackHB_8E29", "HLT_IsoTrackHE_8E29", "HLT_L1Tech_RPC_TTU_RBst1_collisions", "HLT_L1_BscMinBiasOR_BptxPlusORMinus", "HLT_L1Tech_BSC_halo_forPhysicsBackground", "HLT_L1Tech_BSC_HighMultiplicity", "HLT_MinBiasPixel_DoubleIsoTrack5", "HLT_MinBiasPixel_DoubleTrack", "HLT_MinBiasPixel_SingleTrack", "HLT_ZeroBiasPixel_SingleTrack", "HLT_MinBiasBSC", "HLT_StoppedHSCP_8E29", "HLT_Jet15U_HcalNoiseFiltered", "HLT_QuadJet15U", "HLT_DiJetAve30U_8E29", "HLT_DiJetAve15U_8E29", "HLT_FwdJet20U", "HLT_Jet50U", "HLT_Jet30U", "HLT_Jet15U", "HLT_BTagMu_Jet10U", "HLT_DoubleJet15U_ForwardBackward", "HLT_BTagIP_Jet50U", "HLT_DoubleLooseIsoTau15", "HLT_SingleLooseIsoTau20", "HLT_HT100U", "HLT_MET100", "HLT_MET45"};
396      
397      bool triggerPassed = false;
398      triggerStatus = 0x0;
# Line 307 | Line 407 | void MPIntuple::analyze(const edm::Event
407            
408            if (hlNames[i] == MPI_TriggerNames[j])
409              {
410 <              //            cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
410 >              //cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
411                triggerStatus |= 0x01 << j;
412                
413              }
# Line 315 | Line 415 | void MPIntuple::analyze(const edm::Event
415        }
416      }
417    }
418 <  //  cout<< "total status: "<<hex<<triggerStatus<<endl;
419 <  
420 <  if(PFcount > 0 || Gencount > 0 || Vtxcount > 0);  sTree -> Fill();
421 <  
418 >  //cout<< "total status: "<<hex<<triggerStatus<<endl;
419 >
420 >
421 >
422 >  if((pfCount > 3 || genCount > 3) && vtxCount > 0) sTree -> Fill();
423 >
424    jet_ak5PF->Clear();
425    primaryVtx->Clear();
426    jetP4_ak5Gen->Clear();
# Line 332 | Line 434 | void  MPIntuple::beginJob()
434    
435    ntupleFile           = new TFile(rootfilename.c_str(), "RECREATE");
436    sTree                = new TTree("mpiTree", "Tree for Jets");
437 <  
437 >
438 >  h1_trackZ            = new TH1F("h1_trackZ", "#Delta z for associated tracks", 60, -0.3, 0.3);
439 >  h1_nAssociatedTracks = new TH1F("h1_nAssociatedTracks", "Number of tracks associated to jet", 20, -0.5, 19.5);
440 >  h1_trackDxy          = new TH1F("h1_trackDxy", "dxy for all associated tracks", 50, -0.4, 0.4);
441 >  h1_allTrackDeltaZ_PV = new TH1F("h1_allTrackDeltaZ_PV", "#Delta z between jet and primary vertex", 60, -0.3, 0.3);
442 >  h1_allTrackDeltaZ_LJ = new TH1F("h1_allTrackDeltaZ_LJ", "#Delta z between leading jet and other jet vertices", 60, -0.3, 0.3);
443 >  h2_nAssTracksVsJetPt = new TH2F("h2_nAssTracksVsJetPt", "nTracks vs Jet p_{T}", 50, 0, 100, 20, -0.5, 19.5);
444 >  p1_nVtcs             = new TProfile("p1_nVtcs", "<N> per run", 5200.0, 132300.0, 137500.0, 0.0, 6.0);
445 >
446    jet_ak5PF            = new TClonesArray("TCJet");
447    jetP4_ak5Gen         = new TClonesArray("TLorentzVector");
448    primaryVtx           = new TClonesArray("TCPrimaryVtx");
449  
450 +
451    sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0);
452    sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0);
453    sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
# Line 347 | Line 458 | void  MPIntuple::beginJob()
458    sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
459    sTree->Branch("isRealData",&isRealData, "isRealData/i");
460    sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
461 +  sTree->Branch("ptHat",&ptHat, "ptHat/d");
462  
463   }
464  
465   // ------------ method called once each job just after ending the event loop  ------------
466   void MPIntuple::endJob()
467   {
468 +
469 +  ntupleFile->cd();
470 +
471 +  h1_trackZ->Write();
472 +  h1_nAssociatedTracks->Write();
473 +  h1_trackDxy->Write();
474 +  h1_allTrackDeltaZ_PV->Write();
475 +  h1_allTrackDeltaZ_LJ->Write();
476 +  h2_nAssTracksVsJetPt->Write();
477 +  p1_nVtcs->Write();
478 +
479    ntupleFile->Write();
480    ntupleFile->Close();
481  
# Line 368 | Line 491 | bool MPIntuple::triggerDecision(edm::Han
491        triggerPassed = true;
492      }
493      return triggerPassed;
494 <  }
494 > }
495 >
496 > double MPIntuple::sumPtSquared(const Vertex & v)
497 > {
498 >  double sum = 0.;
499 >  double pT;
500 >  for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
501 >    pT = (**it).pt();
502 >    double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
503  
504 +    sum += pT*pT;
505 +  }
506 +  return sum;
507 + }
508  
509   //define this as a plug-in
510   DEFINE_FWK_MODULE(MPIntuple);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines