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
Log Message:
package for Z'/boostedZ analysis

File Contents

# User Rev Content
1 digiovan 1.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