ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UfCode/UserArea/TemplateNtuple/src/TemplateNtuple.cc
Revision: 1.1
Committed: Mon Oct 25 15:29:04 2010 UTC (14 years, 6 months ago) by digiovan
Content type: text/plain
Branch: MAIN
CVS Tags: V2012-H-02, V2012-01-00, V2011-01-01, V2011-01-00, AnnaDimuon, V01-00-01, V01-00-00, HEAD
Error occurred while calculating annotation data.
Log Message:
package for Z'/boostedZ analysis

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: TemplateNtuple
4 // Class: TemplateNtuple
5 //
6 /**\class TemplateNtuple TemplateNtuple.cc UserArea/TemplateNtuple/src/TemplateNtuple.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Theodore Nicholas Kypreos,8 R-031,+41227675208,
15 // Created: Fri Apr 9 14:49:13 CEST 2010
16 // $Id$
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23 #include <vector>
24 #include <iostream>
25 #include <fstream>
26 #include <string>
27
28 #include "TLorentzVector.h"
29 #include "TTree.h"
30 #include "TFile.h"
31 #include "TH1F.h"
32 #include "TH2F.h"
33 #include "TBranch.h"
34
35 // user include files
36 #include "FWCore/Framework/interface/Frameworkfwd.h"
37 #include "FWCore/Framework/interface/EDAnalyzer.h"
38
39 #include "FWCore/Framework/interface/Event.h"
40 #include "FWCore/Framework/interface/MakerMacros.h"
41
42 #include "FWCore/ParameterSet/interface/ParameterSet.h"
43
44 #include "DataFormats/MuonReco/interface/Muon.h"
45 #include "DataFormats/MuonReco/interface/MuonFwd.h"
46 #include "DataFormats/MuonReco/interface/MuonIsolation.h"
47 #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
48 #include "DataFormats/MuonReco/interface/MuonSegmentMatch.h"
49 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
50
51 #include "DataFormats/TrackReco/interface/Track.h"
52 #include "DataFormats/TrackReco/interface/TrackFwd.h"
53 #include "FWCore/Framework/interface/ESHandle.h"
54 #include "FWCore/Framework/interface/Event.h"
55 #include "FWCore/Framework/interface/EventSetup.h"
56
57
58 #include "DataFormats/DetId/interface/DetId.h"
59 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
60 #include "DataFormats/MuonDetId/interface/DTWireId.h"
61 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
62 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
63
64 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
65 #include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h"
66
67
68 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
69
70 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
71
72 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
73 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
74 #include "MagneticField/Engine/interface/MagneticField.h"
75 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
76 //
77 // math classes
78 //
79 #include "DataFormats/Math/interface/deltaPhi.h"
80 //
81 // trigger
82 //
83 #include "DataFormats/Common/interface/TriggerResults.h"
84 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
85 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
86 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
87 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
88 #include "DataFormats/Common/interface/TriggerResults.h"
89 #include "FWCore/Common/interface/TriggerNames.h"
90
91 //
92 // vertexing
93 //
94 #include "DataFormats/VertexReco/interface/Vertex.h"
95 #include "DataFormats/VertexReco/interface/VertexFwd.h"
96 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
97 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
98 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
99 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
100 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
101 //
102 // gen particles
103 //
104 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
105
106
107 //
108 // class declaration
109 //
110 /*
111 class BetterMuon {
112
113 public:
114 reco::Muon muon;
115 std::vector<CSCSegment> segments;
116
117 };
118 */
119 class TemplateNtuple : public edm::EDAnalyzer {
120
121 public:
122 explicit TemplateNtuple(const edm::ParameterSet&);
123 ~TemplateNtuple();
124 int _numEvents;
125 edm::ESHandle<TransientTrackBuilder> transientTrackBuilder;
126
127 edm::ESHandle<GlobalTrackingGeometry> globalTrackingGeometry;
128 MuonServiceProxy* theService;
129
130
131 typedef std::pair<reco::Track,reco::Track> TrackPair;
132 typedef std::vector<TrackPair> TrackPairs;
133
134 typedef std::pair<reco::Muon,reco::Muon> MuonPair;
135 typedef std::vector<MuonPair> MuonPairs;
136
137 typedef std::vector<reco::MuonChamberMatch> MuonChamberMatches;
138 typedef std::pair<reco::MuonChamberMatch, CSCSegment> ChamberSegmentPair;
139 typedef std::vector<ChamberSegmentPair> ChamberSegmentPairs;
140
141 typedef struct {
142 float charge;
143 float pt;
144 float eta;
145 float phi;
146
147 } _TrackInfo;
148
149 _TrackInfo _true1,_true2;
150 _TrackInfo _reco1,_reco2;
151 _TrackInfo _refit1,_refit2;
152
153
154 TFile* _outFile;
155 TTree* _outTree;
156
157
158
159 private:
160 virtual void beginJob() ;
161 virtual void analyze(const edm::Event&, const edm::EventSetup&);
162 virtual void endJob() ;
163
164 edm::InputTag _muonColl;
165
166
167 std::string _getFilename;
168
169
170 MuonPairs const GetMuonPosNegPairs(reco::MuonCollection const* muons) const;
171
172 double const GetInvariantMass(MuonPair const* partpair) const;
173
174
175
176
177
178
179 void DumpMuonCollection(reco::MuonCollection const& muons);
180 void DumpTrackInfo(reco::Track const* track);
181 TLorentzVector const GetLorentzVector(TemplateNtuple::MuonPair const* partpair) const;
182 TransientVertex const GetVertexFromPair(MuonPair const& muPair) const;
183
184 // float _treeMass;
185
186
187
188
189 // ----------member data ---------------------------
190 };
191
192 //
193 // constants, enums and typedefs
194 //
195
196 //
197 // static data member definitions
198 //
199
200 //
201 // constructors and destructor
202 //
203 TemplateNtuple::TemplateNtuple(const edm::ParameterSet& iConfig):_numEvents(0)
204
205 {
206 //now do what ever initialization is needed
207 _getFilename = iConfig.getUntrackedParameter<std::string> ("getFilename" , "test.root");
208 _muonColl = iConfig.getParameter<edm::InputTag>("muonColl");
209
210 edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
211 theService = new MuonServiceProxy(serviceParameters);
212 }
213
214
215 TemplateNtuple::~TemplateNtuple()
216 {
217
218 // do anything here that needs to be done at desctruction time
219 // (e.g. close files, deallocate resources etc.)
220
221 }
222
223 #include <Geometry/Records/interface/MuonGeometryRecord.h>
224 #include <Geometry/CSCGeometry/interface/CSCGeometry.h>
225 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
226 #include <Geometry/CSCGeometry/interface/CSCLayerGeometry.h>
227
228
229
230 //
231 // member functions
232 //
233 // ------------ method called to for each event ------------
234 void
235 TemplateNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
236 {
237
238 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",transientTrackBuilder);
239
240 edm::Handle<reco::VertexCollection> offlinePrimaryVertices;
241 iEvent.getByLabel("offlinePrimaryVertices",offlinePrimaryVertices);
242
243 edm::Handle<reco::BeamSpot> beamSpot;
244 iEvent.getByLabel("offlineBeamSpot", beamSpot);
245 math::XYZPoint theVertexPoint = beamSpot->position();
246 if (!offlinePrimaryVertices->empty()) theVertexPoint = offlinePrimaryVertices->front().position();
247
248
249 edm::Handle<reco::MuonCollection> allRecoMuons;
250 iEvent.getByLabel("muons", allRecoMuons);
251
252 edm::Handle<reco::GenParticleCollection> allGenParticles;
253 iEvent.getByLabel("genParticles", allGenParticles);
254
255
256 // std::cout<<"number of gen particles: "<<allGenParticles->size()<<std::endl;
257
258 reco::GenParticleCollection theGens;
259
260 for (reco::GenParticleCollection::const_iterator gen = allGenParticles->begin(), genEnd = allGenParticles->end();
261 gen != genEnd; ++gen)
262 {
263 int id = gen->pdgId();
264 // if (abs(id) != 13 && abs(id) != 32) continue;
265 // std::cout<<gen->status()<<"\t"<<gen->pdgId()<<"\t"<<gen->mass()<<std::endl;
266 if (gen->status() == 3 && abs(id) == 13) theGens.push_back(*gen);
267 }
268
269 // std::cout<<"num gens: "<<theGens.size()<<std::endl;
270
271 if (theGens.size() != 2) return;
272 TLorentzVector vtrue1, vtrue2, vtrueMother;
273 vtrue1.SetPtEtaPhiM(theGens[0].pt(), theGens[0].eta(), theGens[0].phi(), theGens[0].mass());
274 vtrue2.SetPtEtaPhiM(theGens[1].pt(), theGens[1].eta(), theGens[1].phi(), theGens[1].mass());
275 _true1.charge = theGens[0].charge(); _true1.pt = theGens[0].pt(); _true1.eta = theGens[0].eta(); _true1.phi=theGens[0].phi();
276 _true2.charge = theGens[1].charge(); _true2.pt = theGens[1].pt(); _true2.eta = theGens[1].eta(); _true2.phi=theGens[1].phi();
277
278 // _TrackInfo _true1,_true2;
279 vtrueMother = vtrue1+vtrue2;
280 // std::cout<<"true mass = "<<vtrueMother.M()<<std::endl;
281
282
283 // reco::MuonCollection globalMuons, trackerMuons, standaloneMuons, trackerGlobalMuons, trackerGlobalMuonsBGroup;
284
285 reco::MuonCollection globalMuons;//, trackerMuons;
286
287 for (reco::MuonCollection::const_iterator muonsBegin = allRecoMuons->begin()
288 , muonsEnd = allRecoMuons->end(),muon = muonsBegin;
289 muon != muonsEnd; ++muon){
290 if (!muon->isGlobalMuon()) continue;
291 bool isGood = true;
292 isGood &= (muon->innerTrack()->pt() >= 1);
293 isGood &= (muon->innerTrack()->p() >= 2.5);
294 isGood &= (muon->innerTrack()->numberOfValidHits() >= 13);
295 isGood &= (muon->innerTrack()->normalizedChi2() <= 5);
296 isGood &= (fabs(muon->innerTrack()->dxy(theVertexPoint)) <= 0.3);
297
298 if (!isGood) continue;
299 globalMuons .push_back(*muon);
300 }
301 MuonPairs globalPosNeg= GetMuonPosNegPairs(&globalMuons);
302
303 // std::cout<<"number of pairs: "<<globalPosNeg.size()<<std::endl;
304
305 // _TrackInfo _reco1,_reco2;
306 // _TrackInfo _refit1,_refit2;
307
308 for (MuonPairs::const_iterator pair= globalPosNeg.begin(); pair != globalPosNeg.end(); ++pair){
309
310 reco::Track glb1 = *(pair->first.globalTrack());
311 reco::Track glb2 = *(pair->second.globalTrack());
312
313 _reco1.charge = glb1.charge(); _reco1.pt = glb1.pt(); _reco1.eta = glb1.eta(); _reco1.phi=glb1.phi();
314 _reco2.charge = glb2.charge(); _reco2.pt = glb2.pt(); _reco2.eta = glb2.eta(); _reco2.phi=glb2.phi();
315
316 TLorentzVector mother= GetLorentzVector(&*pair);
317 TransientVertex vtx= GetVertexFromPair(*pair);
318 // std::cout<<"reco mass = "<<mother.M()<<std::endl;
319 if (vtx.isValid() && vtx.refittedTracks().size() == 2){
320 reco::TransientTrack tt1 = vtx.refittedTracks().front();
321 reco::TransientTrack tt2 = vtx.refittedTracks().back();
322
323 GlobalVector gv1 = tt1.trajectoryStateClosestToPoint(vtx.position()).momentum();
324 GlobalVector gv2 = tt2.trajectoryStateClosestToPoint(vtx.position()).momentum();
325
326 double const MASS_MUON = 0.105658367;
327 TLorentzVector vrf1, vrf2, vrfMother;
328 vrf1.SetPtEtaPhiM(gv1.perp(), gv1.eta(), gv1.barePhi(), MASS_MUON);
329 vrf2.SetPtEtaPhiM(gv2.perp(), gv2.eta(), gv2.barePhi(), MASS_MUON);
330 _refit1.charge = tt1.charge(); _refit1.pt = gv1.perp(); _refit1.eta = gv1.eta(); _refit1.phi=gv1.barePhi();
331 _refit2.charge = tt2.charge(); _refit2.pt = gv2.perp(); _refit2.eta = gv2.eta(); _refit2.phi=gv2.barePhi();
332 vrfMother = vrf1+vrf2;
333 // std::cout<<"refit mass = "<<vrfMother.M()<<std::endl;
334
335
336
337 // vtx.isValid(),vtx.position().x(), vtx.position().y(), vtx.position().z(), vtx.totalChiSquared(), vtx.degreesOfFreedom(),
338
339 // vtrue2.SetPtEtaPhiM(.pt(), theGens[1].eta(), theGens[1].phi(), theGens[1].mass());
340
341
342 // tt1.charge(), gv1.perp(), gv1.eta(),gv1.barePhi(),
343 // tt2.charge(), gv2.perp(), gv2.eta(),gv2.barePhi()
344 }
345
346 }
347
348 _outTree->Fill();
349
350 return;
351 }
352
353 /*
354 //
355 // trigger bits
356 //
357 edm::Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
358 iEvent.getByLabel("gtDigis",L1GTRR);
359
360 ///
361 std::cout<<"fire bit 0: "<<L1GTRR->technicalTriggerWord()[0]<<std::endl;
362 std::cout<<"fire bit 36: "<<L1GTRR->technicalTriggerWord()[36]<<std::endl;
363 std::cout<<"fire bit 37: "<<L1GTRR->technicalTriggerWord()[37]<<std::endl;
364 std::cout<<"fire bit 38: "<<L1GTRR->technicalTriggerWord()[38]<<std::endl;
365 std::cout<<"fire bit 39: "<<L1GTRR->technicalTriggerWord()[39]<<std::endl;
366 std::cout<<"fire bit 40: "<<L1GTRR->technicalTriggerWord()[40]<<std::endl;
367 std::cout<<"fire bit 41: "<<L1GTRR->technicalTriggerWord()[41]<<std::endl;
368 std::cout<<"fire bit 124: "<<L1GTRR->technicalTriggerWord()[124]<<std::endl;
369
370 if (!(L1GTRR->technicalTriggerWord()[0])) return;
371 if ((L1GTRR->technicalTriggerWord()[36])) return;
372 if ((L1GTRR->technicalTriggerWord()[37])) return;
373 if ((L1GTRR->technicalTriggerWord()[38])) return;
374 if ((L1GTRR->technicalTriggerWord()[39])) return;
375 ///
376
377
378
379
380
381
382 std::cout<<"passed triggers..."<<std::endl;
383
384 bool passBit40or41 = (L1GTRR->technicalTriggerWord()[40] || L1GTRR->technicalTriggerWord()[41]);
385
386
387 edm::ESHandle<CSCGeometry> cscGeometry;
388 iSetup.get<MuonGeometryRecord>().get(cscGeometry);
389
390 // using namespace edm;
391 // edm::Handle<reco::MuonCollection> allRecoMuons;
392 iEvent.getByLabel("muons", allRecoMuons);
393 // iEvent.getByLabel("globalAndTrackerMuons2", allRecoMuons);
394
395 ++_numEvents;
396
397
398 int theRun = iEvent.id().run();
399 int theLumi = iEvent.luminosityBlock();
400 int theEvent= iEvent.id().event();
401
402 edm::Handle<CSCSegmentCollection> cscSegments;
403 iEvent.getByLabel("cscSegments", cscSegments);
404
405 iSetup.get<GlobalTrackingGeometryRecord>().get(globalTrackingGeometry);
406 theService->update(iSetup);
407
408 //
409 // parse the muon collection into global, tracker, and standalone muons
410 //
411 reco::MuonCollection globalMuons, trackerMuons, standaloneMuons, trackerGlobalMuons, trackerGlobalMuonsBGroup;
412
413 for (reco::MuonCollection::const_iterator muonsBegin = allRecoMuons->begin()
414 , muonsEnd = allRecoMuons->end(),muon = muonsBegin;
415 muon != muonsEnd; ++muon){
416 if (muon->isGlobalMuon() ) globalMuons .push_back(*muon);
417
418 if (muon->isTrackerMuon()) {
419 bool isGood = true;
420 isGood &= (muon->innerTrack()->pt() >= 1);
421 isGood &= (muon->innerTrack()->p() >= 2.5);
422 isGood &= (muon->innerTrack()->numberOfValidHits() >= 13);
423 isGood &= (muon->innerTrack()->normalizedChi2() <= 5);
424 isGood &= (fabs(muon->innerTrack()->dxy(theVertexPoint)) <= 0.3);
425 // isGood &= muon::isGoodMuon(*muon,muon::TMLastStationAngTight);
426 isGood &= muon::isGoodMuon(*muon,muon::TrackerMuonArbitrated);
427
428 if (isGood) trackerMuons .push_back(*muon);
429 // if (isGood && muon->isGlobalMuon()) trkVsGlobalPt->Fill(muon->innerTrack()->pt(), muon->globalTrack()->pt());
430 }
431
432 if (muon->isStandAloneMuon()) standaloneMuons .push_back(*muon);
433 if (muon->isGlobalMuon() || muon->isTrackerMuon()) {
434 reco::TrackRef theTrack = muon->innerTrack();
435 bool isGood = true;
436 isGood &= (theTrack->pt() >= 1);
437 isGood &= (theTrack->p() >= 2.5);
438 isGood &= (theTrack->numberOfValidHits() >= 13);
439 isGood &= (theTrack->normalizedChi2() <= 5);
440 isGood &= (fabs(theTrack->dxy(theVertexPoint)) <= 5.);
441 std::cout<<"num pixel hits: "<<theTrack->hitPattern().numberOfValidPixelHits()<<std::endl;
442 std::cout<<"d0: "<<theTrack->dxy(theVertexPoint)<<std::endl;
443 isGood &= (theTrack->hitPattern().pixelLayersWithMeasurement() >1);
444
445 // if (!isGood) continue;
446 if (muon->isTrackerMuon()){
447 isGood &= muon::isGoodMuon(*muon,muon::TrackerMuonArbitrated);
448 }
449 trackerGlobalMuons.push_back(*muon);
450 }
451 }
452
453 std::cout<<"number of selected muons: "<<trackerGlobalMuons.size()<<std::endl;
454
455
456 MuonPairs trackerGlobalMuonPairsOs = GetMuonPosNegPairs(&trackerGlobalMuons);
457
458 // MuonPairs trackerGlobalMuonPairsBGroupOs = GetMuonPosNegPairs(&trackerGlobalMuonsBGroup);
459
460
461
462 for (MuonPairs::const_iterator pair = trackerGlobalMuonPairsOs.begin();
463 pair != trackerGlobalMuonPairsOs.end(); ++pair){
464
465
466
467 if (true){
468 printf("checking some of my good times..\n");
469 MuonChamberMatches const& matches = pair->first.matches();
470 std::cout<<"number of matches = "<<pair->first.matches().size()<<std::endl;
471 for (MuonChamberMatches::const_iterator mcm = matches.begin(); mcm != matches.end(); ++mcm){
472 // if (mcm->detector() != MuonSubdetId::CSC) continue;
473 DetId const& chamberId = mcm->id;
474
475 if (chamberId.det() != DetId::Muon) continue;
476 if (chamberId.subdetId() != MuonSubdetId::CSC) continue;
477 CSCDetId id(chamberId.rawId());
478 // CSCLayerGeometry* lgeo = &*cscGeometry->layer(id)->geometry();
479 // cscGeometry->layer(id)->geometry()->nearestStrip(LocalPoint(mcm->x,mcm->y,0));
480 // long channel = chamberId.rawId()*16 + ((stripNumber-1)%16);
481 std::cout<<"CSC match:"
482 <<"\t"<<"x = "<<mcm->x
483 <<"\t"<<"y = "<<mcm->y
484 // <<"\t"<<"stripNumber = "<<stripNumber
485 // <<"\t"<<"channel = "<<channel
486 // <<"\t"<<"id.channel = "<<id.channel()
487 <<std::endl;
488 GlobalPoint gp = theService->trackingGeometry()->idToDet(id)->surface().toGlobal(LocalPoint(mcm->x,mcm->y,0));
489 std::cout<<gp<<std::endl;
490 Surface const& s = cscGeometry->idToDet(chamberId)->surface();
491 GlobalPoint gp1 = s.toGlobal(LocalPoint(mcm->x,mcm->y,0));
492 std::cout<<gp1<<std::endl;
493 }
494
495 }
496
497
498 // std::cout<<"compareDR = "<<deltaR<<"\t"<<deltaR1<<std::endl;
499
500 // double mass = GetInvariantMass(&*pair);
501 // trkOsDPhiChambArb2VsMass ->Fill(mass, deltaPhi);
502 // trkOsDRhoChambArb2VsMass ->Fill(mass,deltaRho);
503
504 bool tmLastStationAngTight = muon::isGoodMuon(pair->first,muon::TMLastStationAngTight);
505 tmLastStationAngTight &= muon::isGoodMuon(pair->second,muon::TMLastStationAngTight);
506 bool overlap = muon::overlap(pair->first, pair->second, 1, 1, true);
507 TransientVertex vtx = GetVertexFromPair(*pair);
508
509 // std::cout<<"number of refitted tracks: "<<vtx.refittedTracks().size()<<std::endl;
510 if (vtx.isValid() && vtx.refittedTracks().size() == 2){
511
512 // reco::Track const& refitTrack1 = vtx.refittedTracks().front().track();
513 // reco::Track const& refitTrack2 = vtx.refittedTracks().back().track();
514 reco::TransientTrack tt1 = vtx.refittedTracks().front();
515 reco::TransientTrack tt2 = vtx.refittedTracks().back();
516
517 GlobalVector gv1 = tt1.trajectoryStateClosestToPoint(vtx.position()).momentum();
518 GlobalVector gv2 = tt2.trajectoryStateClosestToPoint(vtx.position()).momentum();
519
520 // isGood &= (muon->innerTrack()->numberOfValidHits() >= 13);
521 // isGood &= (muon->innerTrack()->normalizedChi2() <= 5);
522 // isGood &= (fabs(muon->innerTrack()->dxy(theVertexPoint)) <= 0.3);
523 // isGood &= trk->hitPattern().pixelLayersWithMeasurement() > 1;
524 // isGood &= (trk->hitPattern().numberOfValidHits() > 11);
525 // isGood &= fabs(trk->dxy(beamSpot)) < 5.;
526 // isGood &= fabs(trk->dz(beamSpot)) < 20;
527
528
529 reco::Muon muon1 = pair->first;
530 reco::TrackRef track1 = muon1.innerTrack();
531
532 // std::cout<<"algo = "<<track1->algo()<<std::endl;
533 // std::cout<<"quality = "<<track1->qualityMask()<<std::endl;
534 // std::cout<<"high purity: = "<<track1->quality(reco::TrackBase::highPurity)<<std::endl;
535 // std::cout<<"confirmed: = "<<track1->quality(reco::TrackBase::confirmed)<<std::endl;
536 // std::cout<<"loose: = "<<track1->quality(reco::TrackBase::loose)<<std::endl;
537 printf("badger-%d:%d:%d:%d", theRun, theLumi, theEvent,passBit40or41);
538
539 printf(
540 "\t%d:%f:%f:%f"
541 ":%d:%d:%d:%d"
542 ":%d:%d"
543 ":%f:%f"
544 ":%f:%f"
545 ":%d:%d:%f"
546 ":%f"
547 ,
548 track1->charge(), track1->pt(), track1->eta(), track1->phi()
549 ,muon1.isGlobalMuon(), muon1.isTrackerMuon(),
550 muon::isGoodMuon(muon1, muon::TMLastStationAngTight),
551 muon::isGoodMuon(muon1, muon::GlobalMuonPromptTight)
552 ,track1->hitPattern().numberOfValidHits(), track1->hitPattern().pixelLayersWithMeasurement()
553 ,track1->dxy(theVertexPoint), track1->dz(theVertexPoint)
554 ,track1->chi2(), track1->ndof()
555 ,track1->algo(), track1->qualityMask(),track1->ptError(),
556
557 );
558
559 muon1 = pair->second;
560 track1 = muon1.innerTrack();
561 //
562 //
563 //
564 printf(
565 "\t%d:%f:%f:%f"
566 ":%d:%d:%d:%d"
567 ":%d:%d"
568 ":%f:%f"
569 ":%f:%f"
570 ":%d:%d:%f"
571 ":%f"
572 ,
573 track1->charge(), track1->pt(), track1->eta(), track1->phi()
574 ,muon1.isGlobalMuon(), muon1.isTrackerMuon(),
575 muon::isGoodMuon(muon1, muon::TMLastStationAngTight),
576 muon::isGoodMuon(muon1, muon::GlobalMuonPromptTight)
577 ,track1->hitPattern().numberOfValidHits(), track1->hitPattern().pixelLayersWithMeasurement()
578 ,track1->dxy(theVertexPoint), track1->dz(theVertexPoint)
579 ,track1->chi2(), track1->ndof()
580 ,track1->algo(), track1->qualityMask(),track1->ptError(),
581 );
582 //
583 // print the rest
584 //
585
586 printf(
587 "\t%d:%d:%f:%d:%f"
588 "\t%d:%f:%f:%f:%f:%f"
589 "\t%f:%f:%f"
590 "\t%f:%f:%f"
591 "\t%d:%f:%f:%f"
592 "\t%d:%f:%f:%f"
593 "\n",
594 shareChamber?1:0, tmLastStationAngTight?1:0,deltaR, overlap?1:0,deltaR1,
595 vtx.isValid(),vtx.position().x(), vtx.position().y(), vtx.position().z(), vtx.totalChiSquared(), vtx.degreesOfFreedom(),
596 theVertexPoint.x(), theVertexPoint.y(), theVertexPoint.z(),
597 beamSpot->x0(), beamSpot->y0(), beamSpot->z0(),
598 tt1.charge(), gv1.perp(), gv1.eta(),gv1.barePhi(),
599 tt2.charge(), gv2.perp(), gv2.eta(),gv2.barePhi()
600 );
601
602 }
603 }
604
605
606 for (MuonPairs::const_iterator pair = trackerGlobalMuonPairsSs.begin();
607 pair != trackerGlobalMuonPairsSs.end(); ++pair){
608
609
610 // double mass = GetInvariantMass(&*pair);
611 // trkOsDPhiChambArb2VsMass ->Fill(mass, deltaPhi);
612 // trkOsDRhoChambArb2VsMass ->Fill(mass,deltaRho);
613
614 bool tmLastStationAngTight = muon::isGoodMuon(pair->first,muon::TMLastStationAngTight);
615 tmLastStationAngTight &= muon::isGoodMuon(pair->second,muon::TMLastStationAngTight);
616 bool overlap = muon::overlap(pair->first, pair->second, 1, 1, true);
617
618 TransientVertex vtx = GetVertexFromPair(*pair);
619
620 // std::cout<<"number of refitted tracks: "<<vtx.refittedTracks().size()<<std::endl;
621 if (vtx.isValid() && vtx.refittedTracks().size() == 2){
622
623 // reco::Track const& refitTrack1 = vtx.refittedTracks().front().track();
624 // reco::Track const& refitTrack2 = vtx.refittedTracks().back().track();
625 reco::TransientTrack tt1 = vtx.refittedTracks().front();
626 reco::TransientTrack tt2 = vtx.refittedTracks().back();
627
628 GlobalVector gv1 = tt1.trajectoryStateClosestToPoint(vtx.position()).momentum();
629 GlobalVector gv2 = tt2.trajectoryStateClosestToPoint(vtx.position()).momentum();
630
631 // isGood &= (muon->innerTrack()->numberOfValidHits() >= 13);
632 // isGood &= (muon->innerTrack()->normalizedChi2() <= 5);
633 // isGood &= (fabs(muon->innerTrack()->dxy(theVertexPoint)) <= 0.3);
634 // isGood &= trk->hitPattern().pixelLayersWithMeasurement() > 1;
635 // isGood &= (trk->hitPattern().numberOfValidHits() > 11);
636 // isGood &= fabs(trk->dxy(beamSpot)) < 5.;
637 // isGood &= fabs(trk->dz(beamSpot)) < 20;
638
639 printf("badger-%d:%d:%d:%d", theRun, theLumi, theEvent,passBit40or41);
640
641 reco::Muon muon1 = pair->first;
642 reco::TrackRef track1 = muon1.innerTrack();
643 printf(
644 "\t%d:%f:%f:%f"
645 ":%d:%d:%d:%d"
646 ":%d:%d"
647 ":%f:%f"
648 ":%f:%f"
649 ":%d:%d:%f"
650 ":%f"
651 ,
652 track1->charge(), track1->pt(), track1->eta(), track1->phi()
653 ,muon1.isGlobalMuon(), muon1.isTrackerMuon(),
654 muon::isGoodMuon(muon1, muon::TMLastStationAngTight),
655 muon::isGoodMuon(muon1, muon::GlobalMuonPromptTight)
656 ,track1->hitPattern().numberOfValidHits(), track1->hitPattern().pixelLayersWithMeasurement()
657 ,track1->dxy(theVertexPoint), track1->dz(theVertexPoint)
658 ,track1->chi2(), track1->ndof()
659 ,track1->algo(), track1->qualityMask(),track1->ptError(),
660 );
661
662 muon1 = pair->second;
663 track1 = muon1.innerTrack();
664 printf(
665 "\t%d:%f:%f:%f"
666 ":%d:%d:%d:%d"
667 ":%d:%d"
668 ":%f:%f"
669 ":%f:%f"
670 ":%d:%d:%f"
671 ":%f"
672 ,
673 track1->charge(), track1->pt(), track1->eta(), track1->phi()
674 ,muon1.isGlobalMuon(), muon1.isTrackerMuon(),
675 muon::isGoodMuon(muon1, muon::TMLastStationAngTight),
676 muon::isGoodMuon(muon1, muon::GlobalMuonPromptTight)
677 ,track1->hitPattern().numberOfValidHits(), track1->hitPattern().pixelLayersWithMeasurement()
678 ,track1->dxy(theVertexPoint), track1->dz(theVertexPoint)
679 ,track1->chi2(), track1->ndof()
680 ,track1->algo(), track1->qualityMask(),track1->ptError(),
681 );
682 //
683 //
684 //
685 //
686 // print the rest
687 //
688
689 printf(
690 "\t%d:%d:%f:%d:%f"
691 "\t%d:%f:%f:%f:%f:%f"
692 "\t%f:%f:%f"
693 "\t%f:%f:%f"
694 "\t%d:%f:%f:%f"
695 "\t%d:%f:%f:%f"
696 "\n",
697 shareChamber?1:0, tmLastStationAngTight?1:0,deltaR, overlap?1:0,deltaR1,
698 vtx.isValid(),vtx.position().x(), vtx.position().y(), vtx.position().z(), vtx.totalChiSquared(), vtx.degreesOfFreedom(),
699 theVertexPoint.x(), theVertexPoint.y(), theVertexPoint.z(),
700 beamSpot->x0(), beamSpot->y0(), beamSpot->z0(),
701 tt1.charge(), gv1.perp(), gv1.eta(),gv1.barePhi(),
702 tt2.charge(), gv2.perp(), gv2.eta(),gv2.barePhi()
703 );
704
705 }
706 }
707
708 //
709 //
710 //
711 //
712 //
713 //
714
715
716 return;
717 }
718 */
719
720 // ------------ method called once each job just before starting event loop ------------
721 void
722 TemplateNtuple::beginJob()
723 {
724 // _outFile = new TFile(_getFilename.c_str(), "recreate");
725 _outFile = new TFile("badger.root","recreate");
726 _outFile->cd();
727 _outTree = new TTree("tree", "myTree");
728 // _outTree->Branch("mass", &_treeMass);
729
730 // _TrackInfo;
731 _outTree->Branch("true1" , &_true1 , "charge:pt:eta:phi");
732 _outTree->Branch("true2" , &_true2 , "charge:pt:eta:phi");
733 _outTree->Branch("reco1" , &_reco1 , "charge:pt:eta:phi");
734 _outTree->Branch("reco2" , &_reco2 , "charge:pt:eta:phi");
735 _outTree->Branch("refit1" , &_refit1 , "charge:pt:eta:phi");
736 _outTree->Branch("refit2" , &_refit2 , "charge:pt:eta:phi");
737
738
739
740 }
741
742 // ------------ method called once each job just after ending the event loop ------------
743 void
744 TemplateNtuple::endJob() {
745
746 std::cout<<"number of candidate di-muon events: "<<_numEvents<<std::endl;
747 _outFile->cd();
748 _outTree->Write();
749 // _outFile->Close();
750 }
751 //
752 //
753 //
754 TemplateNtuple::MuonPairs const TemplateNtuple::GetMuonPosNegPairs(reco::MuonCollection const* muons) const {
755
756 reco::MuonCollection posMuons, negMuons; posMuons.clear(); negMuons.clear();
757
758 for (reco::MuonCollection::const_iterator muonsBegin = muons->begin()
759 , muonsEnd = muons->end(),muon = muonsBegin;
760 muon != muonsEnd; ++muon){
761
762 if (muon->charge() > 0) posMuons.push_back(*muon);
763 else if (muon->charge() < 0) negMuons.push_back(*muon);
764 }
765
766 MuonPairs muonpairs; muonpairs.clear();
767 for (reco::MuonCollection::const_iterator pos = posMuons.begin(); pos != posMuons.end(); ++pos)
768 for (reco::MuonCollection::const_iterator neg = negMuons.begin(); neg != negMuons.end(); ++neg)
769 muonpairs.push_back(TemplateNtuple::MuonPair(*pos,*neg));
770
771 return muonpairs;
772 }
773 TLorentzVector const TemplateNtuple::GetLorentzVector(TemplateNtuple::MuonPair const* partpair) const {
774 TLorentzVector v1, v2, vsum;
775 double const MASS_MUON = 0.105658367;
776
777 v1.SetPtEtaPhiM(partpair->first.pt() , partpair->first.eta() , partpair->first.phi() , MASS_MUON);
778 v2.SetPtEtaPhiM(partpair->second.pt() , partpair->second.eta(), partpair->second.phi(), MASS_MUON);
779
780 vsum = v1+v2;
781 return vsum;
782 }
783
784 double const TemplateNtuple::GetInvariantMass(TemplateNtuple::MuonPair const* partpair) const {
785 TLorentzVector vsum = GetLorentzVector(partpair);
786 return vsum.M();
787 }
788
789 //define this as a plug-in
790 void TemplateNtuple::DumpMuonCollection(reco::MuonCollection const& muons){
791 for (reco::MuonCollection::const_iterator muonsBegin = muons.begin()
792 , muonsEnd = muons.end(),muon = muonsBegin;
793 muon != muonsEnd; ++muon){
794 }
795 }
796 //
797 //
798 //
799 void TemplateNtuple::DumpTrackInfo(reco::Track const* track) {
800
801 std::cout<<"<track>"
802 <<"\t"<<track->pt()
803 <<"\t"<<track->p()
804 <<"\t"<<track->phi()
805 // <<"\t"<<track->dxy(beamSpot->position())
806 <<"\t"<<track->numberOfValidHits()
807 <<"\t"<<track->normalizedChi2()
808 <<std::endl;
809 }
810 //
811 //
812 //
813 TransientVertex const TemplateNtuple::GetVertexFromPair(MuonPair const& muPair) const {
814 // reco::TransientTrack ttrk1 = (*transientTrackBuilder).build(*muPair.first.innerTrack());
815 // reco::TransientTrack ttrk2 = (*transientTrackBuilder).build(*muPair.second.innerTrack());
816 reco::TransientTrack ttrk1 = (*transientTrackBuilder).build(*muPair.first.globalTrack());
817 reco::TransientTrack ttrk2 = (*transientTrackBuilder).build(*muPair.second.globalTrack());
818
819 std::vector<reco::TransientTrack> t_trks; t_trks.clear();
820 t_trks.push_back(ttrk1);
821 t_trks.push_back(ttrk2);
822
823
824 KalmanVertexFitter kvf(true); // false means no smoothing which means no track re-fit
825 TransientVertex tv = kvf.vertex(t_trks);
826 return tv;
827
828 }
829
830 DEFINE_FWK_MODULE(TemplateNtuple);
831
832