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 |
|