ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/CutFlowAnalyzer/src/CutFlowAnalyzer.cc
Revision: 1.10
Committed: Thu Jun 27 06:37:35 2013 UTC (11 years, 10 months ago) by pakhotin
Content type: text/plain
Branch: MAIN
Changes since 1.9: +114 -61 lines
Log Message:
add PF isolation

File Contents

# User Rev Content
1 pakhotin 1.1 // -*- C++ -*-
2     //
3     // Package: CutFlowAnalyzer
4     // Class: CutFlowAnalyzer
5     //
6     /**\class CutFlowAnalyzer CutFlowAnalyzer.cc AnalysisTools/CutFlowAnalyzer/src/CutFlowAnalyzer.cc
7    
8     Description: [one line class summary]
9    
10     Implementation:
11     [Notes on implementation]
12     */
13     //
14     // Original Author: Yuriy Pakhotin
15     // Created: Tue Feb 19 18:51:12 CST 2013
16 pakhotin 1.9 // $Id: CutFlowAnalyzer.cc,v 1.8 2013/06/04 05:54:05 pakhotin Exp $
17 pakhotin 1.1 //
18     //
19    
20    
21     // system include files
22     #include <memory>
23    
24 pakhotin 1.2
25 pakhotin 1.1 // user include files
26     #include "FWCore/Framework/interface/Frameworkfwd.h"
27     #include "FWCore/Framework/interface/EDAnalyzer.h"
28    
29     #include "FWCore/Framework/interface/Event.h"
30     #include "FWCore/Framework/interface/MakerMacros.h"
31    
32     #include "FWCore/ParameterSet/interface/ParameterSet.h"
33 pakhotin 1.2
34     // user include files
35     #include "TTree.h"
36     #include "TRandom3.h"
37     #include "TMath.h"
38    
39     #include "CommonTools/UtilAlgos/interface/TFileService.h"
40     #include "DataFormats/Candidate/interface/Candidate.h"
41     #include "DataFormats/DetId/interface/DetId.h"
42     #include "DataFormats/GeometrySurface/interface/Cylinder.h"
43     #include "DataFormats/GeometrySurface/interface/Plane.h"
44     #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
45     #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
46     #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
47     #include "DataFormats/MuonReco/interface/Muon.h"
48     #include "DataFormats/MuonReco/interface/MuonSegmentMatch.h"
49     #include "DataFormats/PatCandidates/interface/Muon.h"
50     #include "DataFormats/PatCandidates/interface/TriggerEvent.h"
51     #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
52     #include "DataFormats/TrackReco/interface/HitPattern.h"
53     #include "DataFormats/TrackReco/interface/TrackFwd.h"
54     #include "DataFormats/TrackReco/interface/Track.h"
55     #include "DataFormats/VertexReco/interface/VertexFwd.h"
56     #include "DataFormats/VertexReco/interface/Vertex.h"
57     #include "FWCore/Framework/interface/EDAnalyzer.h"
58     #include "FWCore/Framework/interface/ESHandle.h"
59     #include "FWCore/Framework/interface/Event.h"
60     #include "FWCore/Framework/interface/Frameworkfwd.h"
61     #include "FWCore/Framework/interface/MakerMacros.h"
62     #include "FWCore/ParameterSet/interface/ParameterSet.h"
63     #include "FWCore/ServiceRegistry/interface/Service.h"
64     #include "FWCore/Utilities/interface/InputTag.h"
65     #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
66     #include "Geometry/CommonTopologies/interface/PixelTopology.h"
67     #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
68     #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
69     #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
70     #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
71     #include "MagneticField/Engine/interface/MagneticField.h"
72     #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
73     #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
74     #include "TrackingTools/GeomPropagators/interface/Propagator.h"
75     #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
76     #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
77     #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
78     #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
79    
80     #include "AnalysisDataFormats/MuJetAnalysis/interface/MultiMuon.h"
81    
82 pakhotin 1.3 //******************************************************************************
83     // Auxiliary function: Order objects by pT
84     //******************************************************************************
85     bool PtOrder (const reco::GenParticle* p1, const reco::GenParticle* p2) { return (p1->pt() > p2->pt() ); }
86    
87     //******************************************************************************
88     // Auxiliary function: Calculate difference between two angles: -PI < phi < PI
89     //******************************************************************************
90     double My_dPhi (double phi1, double phi2) {
91     double dPhi = phi1 - phi2;
92     if (dPhi > M_PI) dPhi -= 2.*M_PI;
93     if (dPhi < -M_PI) dPhi += 2.*M_PI;
94     return dPhi;
95     }
96    
97 pakhotin 1.7 // Loose Muons
98     bool isPFMuonLoose (const reco::Muon* mu) {
99     bool isMoonLoose = false;
100     if ( fabs(mu->eta()) < 2.4
101     // && mu->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2
102     && ( mu->isTrackerMuon() || mu->isGlobalMuon() ) ) {
103     isMoonLoose = true;
104     }
105     return isMoonLoose;
106     }
107    
108     // Private ID Muons
109     bool isTrackerMuonPrivateID (const reco::Muon* mu) {
110     bool isTrackerMuonPrivateID = false;
111     if ( fabs(mu->eta()) < 2.4
112     && mu->isTrackerMuon()
113     && mu->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2
114     && mu->innerTrack()->numberOfValidHits() >= 8
115     && mu->innerTrack()->normalizedChi2() < 4. ) {
116     isTrackerMuonPrivateID = true;
117     }
118     return isTrackerMuonPrivateID;
119     }
120    
121 pakhotin 1.3 //******************************************************************************
122     // Class declaration
123     //******************************************************************************
124 pakhotin 1.1
125     class CutFlowAnalyzer : public edm::EDAnalyzer {
126 pakhotin 1.3 public:
127     explicit CutFlowAnalyzer(const edm::ParameterSet&);
128     ~CutFlowAnalyzer();
129    
130     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
131    
132     private:
133     virtual void beginJob() ;
134     virtual void analyze(const edm::Event&, const edm::EventSetup&);
135     virtual void endJob() ;
136    
137     virtual void beginRun(edm::Run const&, edm::EventSetup const&);
138     virtual void endRun(edm::Run const&, edm::EventSetup const&);
139     virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
140     virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
141 pakhotin 1.2
142 pakhotin 1.3 // ---------- Member data ----------
143 pakhotin 1.2
144 pakhotin 1.3 TTree * m_ttree; // Store some variables in this tree for later access
145    
146     Int_t m_run; // run number | these three numbers required to extract event
147     Int_t m_lumi; // lumi number | from sample (data or MC) and examine it in
148     Int_t m_event; // event number | event display
149    
150 pakhotin 1.5 Float_t b_beamSpot_x;
151     Float_t b_beamSpot_y;
152     Float_t b_beamSpot_z;
153    
154 pakhotin 1.3 Int_t m_events; // total number of events
155 pakhotin 1.4
156     Float_t m_Mu17_threshold;
157     Float_t m_Mu8_threshold;
158    
159 pakhotin 1.3 // ---------- Generator Level ----------
160    
161     edm::InputTag m_genParticles; // Label to access generator particles
162    
163     // Events counters
164     Int_t m_events4GenMu; // number of events with 4 gen muons
165     Int_t m_events1GenMu17; // n events with 1 gen muon: pT > 17 GeV, |eta| < 0.9
166     Int_t m_events2GenMu8; // n events with 2 gen muons: pT > 8 GeV, |eta| < 2.4
167     Int_t m_events3GenMu8; // n events with 3 gen muons: pT > 8 GeV, |eta| < 2.4
168     Int_t m_events4GenMu8; // n events with 4 gen muons: pT > 8 GeV, |eta| < 2.4
169    
170     // Gen branches in ROOT tree (they all start with b_)
171 pakhotin 1.5 Float_t b_genH_m;
172 pakhotin 1.3 Float_t b_genH_p;
173     Float_t b_genH_pT;
174 pakhotin 1.5 Float_t b_genH_eta;
175     Float_t b_genH_phi;
176 pakhotin 1.4 Float_t b_genH_vx;
177     Float_t b_genH_vy;
178     Float_t b_genH_vz;
179    
180 pakhotin 1.5 Float_t b_genA0_m;
181 pakhotin 1.3 Float_t b_genA0_p;
182     Float_t b_genA0_pT;
183     Float_t b_genA0_eta;
184 pakhotin 1.5 Float_t b_genA0_phi;
185     Float_t b_genA0_vx;
186     Float_t b_genA0_vy;
187     Float_t b_genA0_vz;
188    
189     Float_t b_genA1_m;
190 pakhotin 1.3 Float_t b_genA1_p;
191     Float_t b_genA1_pT;
192     Float_t b_genA1_eta;
193 pakhotin 1.5 Float_t b_genA1_phi;
194     Float_t b_genA1_vx;
195     Float_t b_genA1_vy;
196     Float_t b_genA1_vz;
197    
198     Float_t b_genA0Mu0_p;
199     Float_t b_genA0Mu1_p;
200     Float_t b_genA1Mu0_p;
201     Float_t b_genA1Mu1_p;
202 pakhotin 1.3
203     Float_t b_genA0Mu0_pT;
204     Float_t b_genA0Mu1_pT;
205     Float_t b_genA1Mu0_pT;
206     Float_t b_genA1Mu1_pT;
207 pakhotin 1.4
208 pakhotin 1.5 Float_t b_genA0Mu0_eta;
209     Float_t b_genA0Mu1_eta;
210     Float_t b_genA1Mu0_eta;
211     Float_t b_genA1Mu1_eta;
212    
213     Float_t b_genA0Mu0_phi;
214     Float_t b_genA0Mu1_phi;
215     Float_t b_genA1Mu0_phi;
216     Float_t b_genA1Mu1_phi;
217    
218 pakhotin 1.4 Float_t b_genA0Mu0_vx;
219     Float_t b_genA0Mu1_vx;
220     Float_t b_genA1Mu0_vx;
221     Float_t b_genA1Mu1_vx;
222    
223     Float_t b_genA0Mu0_vy;
224     Float_t b_genA0Mu1_vy;
225     Float_t b_genA1Mu0_vy;
226     Float_t b_genA1Mu1_vy;
227    
228     Float_t b_genA0Mu0_vz;
229     Float_t b_genA0Mu1_vz;
230     Float_t b_genA1Mu0_vz;
231     Float_t b_genA1Mu1_vz;
232    
233 pakhotin 1.5 Float_t b_genA0_vLT;
234     Float_t b_genA1_vLT;
235    
236     Float_t b_genA0_vL;
237     Float_t b_genA1_vL;
238    
239     Float_t b_genA0Mu_dEta;
240     Float_t b_genA1Mu_dEta;
241     Float_t b_genA0Mu_dPhi;
242     Float_t b_genA1Mu_dPhi;
243 pakhotin 1.3 Float_t b_genA0Mu_dR;
244     Float_t b_genA1Mu_dR;
245    
246     Float_t b_genMu0_pT;
247     Float_t b_genMu1_pT;
248     Float_t b_genMu2_pT;
249     Float_t b_genMu3_pT;
250    
251 pakhotin 1.6 Float_t b_genMu0_eta;
252     Float_t b_genMu1_eta;
253     Float_t b_genMu2_eta;
254     Float_t b_genMu3_eta;
255    
256     Float_t b_genMu0_phi;
257     Float_t b_genMu1_phi;
258     Float_t b_genMu2_phi;
259     Float_t b_genMu3_phi;
260    
261 pakhotin 1.3 // ---------- RECO Level ----------
262    
263     // Labels to access
264     edm::InputTag m_muons; // reconstructed muons
265     edm::InputTag m_muJets; // muon jets built from reconstructed muons
266    
267     // Events counters
268     Int_t m_events1SelMu17; // n events with 1 selected reco muon: pT > 17 GeV, |eta| < 0.9
269     Int_t m_events2SelMu8; // n events with 2 selected reco muon: pT > 8 GeV, |eta| < 2.4
270     Int_t m_events3SelMu8; // n events with 2 selected reco muon: pT > 8 GeV, |eta| < 2.4
271     Int_t m_events4SelMu8; // n events with 2 selected reco muon: pT > 8 GeV, |eta| < 2.4
272     Int_t m_events2MuJets; // n events with 2 muon jets
273     Int_t m_events2DiMuons; // n events with 2 dimuons (dimuon = muon jet with 2 muons)
274     Int_t m_eventsDz2DiMuonsOK; // n events with dz between dimuons less than 0.1
275     Int_t m_eventsDiMuonsHLTFired;// n events with dimuon HLT fired
276     Int_t m_eventsMassDiMuonsOK; // n events with both dimuons masses in diagonal
277     Int_t m_eventsIsoDiMuonsOK; // n events with both dimuons isolated
278     Int_t m_eventsVertexOK; // n events with one good vertex
279    
280     // Cuts thresholds from Python config file
281     Double_t m_maxIsoDiMuons;
282    
283     // Auxiliary variables
284     TRandom3 m_trandom3;
285    
286     // Reco branches in ROOT tree (they all start with b_)
287 pakhotin 1.6
288     Float_t b_selMu0_pT;
289     Float_t b_selMu1_pT;
290     Float_t b_selMu2_pT;
291     Float_t b_selMu3_pT;
292    
293     Float_t b_selMu0_eta;
294     Float_t b_selMu1_eta;
295     Float_t b_selMu2_eta;
296     Float_t b_selMu3_eta;
297    
298     Float_t b_selMu0_phi;
299     Float_t b_selMu1_phi;
300     Float_t b_selMu2_phi;
301     Float_t b_selMu3_phi;
302    
303 pakhotin 1.10 Float_t b_diMuonC_isoTk;
304     Float_t b_diMuonF_isoTk;
305     Float_t b_diMuonC_isoPF;
306     Float_t b_diMuonF_isoPF;
307    
308     Float_t b_dZDiMuonC;
309     Float_t b_dZDiMuonF;
310     Float_t b_dZDiMuons;
311    
312     // Bool branches in ROOT tree
313    
314     Bool_t b_is1SelMu17;
315     Bool_t b_is4SelMu8;
316     Bool_t b_is2MuJets;
317     Bool_t b_is2DiMuons;
318     Bool_t b_isDzDiMuonsOK;
319     Bool_t b_isDiMuonHLTFired;
320     Bool_t b_isMassDiMuonsOK;
321     Bool_t b_isIsoDiMuonsOK;
322     Bool_t b_isVLT;
323    
324 pakhotin 1.1 };
325    
326     //
327     // constants, enums and typedefs
328     //
329    
330     //
331     // static data member definitions
332     //
333    
334     //
335     // constructors and destructor
336     //
337     CutFlowAnalyzer::CutFlowAnalyzer(const edm::ParameterSet& iConfig)
338    
339     {
340 pakhotin 1.5
341 pakhotin 1.3 m_events = 0;
342     m_ttree = NULL;
343    
344     // ---------- Generator Level ----------
345     m_genParticles = iConfig.getParameter<edm::InputTag>("genParticles");
346    
347     m_events4GenMu = 0;
348     m_events1GenMu17 = 0;
349     m_events2GenMu8 = 0;
350     m_events3GenMu8 = 0;
351     m_events4GenMu8 = 0;
352    
353     // ---------- RECO Level ----------
354     m_muons = iConfig.getParameter<edm::InputTag>("muons");
355     m_muJets = iConfig.getParameter<edm::InputTag>("muJets");
356    
357     m_maxIsoDiMuons = iConfig.getParameter<double>("maxIsoDiMuons");
358 pakhotin 1.5 m_Mu17_threshold = 17.0;
359     m_Mu8_threshold = 8.0;
360 pakhotin 1.3
361     m_events1SelMu17 = 0;
362     m_events2SelMu8 = 0;
363     m_events3SelMu8 = 0;
364     m_events4SelMu8 = 0;
365     m_events2MuJets = 0;
366     m_events2DiMuons = 0;
367     m_eventsDz2DiMuonsOK = 0;
368     m_eventsDiMuonsHLTFired = 0;
369     m_eventsMassDiMuonsOK = 0;
370     m_eventsIsoDiMuonsOK = 0;
371     m_eventsVertexOK = 0;
372 pakhotin 1.1 }
373    
374    
375     CutFlowAnalyzer::~CutFlowAnalyzer()
376     {
377    
378     // do anything here that needs to be done at desctruction time
379     // (e.g. close files, deallocate resources etc.)
380    
381     }
382    
383     //
384     // member functions
385     //
386    
387     // ------------ method called for each event ------------
388     void
389     CutFlowAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
390     {
391     using namespace edm;
392 pakhotin 1.3
393     // get the run, lumi and event number
394     m_run = iEvent.id().run();
395     m_lumi = iEvent.id().luminosityBlock();
396 pakhotin 1.2 m_event = iEvent.id().event();
397    
398 pakhotin 1.5 edm::Handle<reco::BeamSpot> beamSpot;
399     iEvent.getByLabel("offlineBeamSpot",beamSpot);
400     b_beamSpot_x = beamSpot->position().x();
401     b_beamSpot_y = beamSpot->position().y();
402     b_beamSpot_z = beamSpot->position().z();
403    
404     // std::cout << "Beam spot x: " << m_beamSpot_x << " y: " << m_beamSpot_y << " z: " << m_beamSpot_z << std::endl;
405    
406 pakhotin 1.3 // count events
407     m_events++;
408     if ( !(m_events%1000) ) std::cout << "Event " << m_events << std::endl;
409    
410     //****************************************************************************
411     // GEN LEVEL
412     //****************************************************************************
413 pakhotin 1.2
414     edm::Handle<reco::GenParticleCollection> genParticles;
415 pakhotin 1.3 iEvent.getByLabel(m_genParticles, genParticles);
416 pakhotin 1.2
417 pakhotin 1.3 // Loop over all genParticles and save prompt muons from particles with codes 36 (a1) or 3000022 (gammaD) in vector genMuons
418     std::vector<const reco::GenParticle*> genH;
419     std::vector<const reco::GenParticle*> genA;
420 pakhotin 1.2 std::vector<const reco::GenParticle*> genMuons;
421 pakhotin 1.3 std::vector<const reco::Candidate*> genMuonMothers;
422     // Loop over all gen particles
423 pakhotin 1.4 int counterGenParticle = 0;
424 pakhotin 1.2 for(reco::GenParticleCollection::const_iterator iGenParticle = genParticles->begin(); iGenParticle != genParticles->end(); ++iGenParticle) {
425 pakhotin 1.4 counterGenParticle++;
426     // std::cout << counterGenParticle << " " << iGenParticle->status() << " " << iGenParticle->pdgId() << " " << iGenParticle->vx() << " " << iGenParticle->vy() << " " << iGenParticle->vz() << std::endl;
427 pakhotin 1.3 // Check if gen particle is muon (pdgId = +/-13) and stable (status = 1)
428     if ( fabs( iGenParticle->pdgId() ) == 13 && iGenParticle->status() == 1 ) {
429     // Mother of the muon can be muon. Find the last muon in this chain: genMuonCand
430     // Example: a1 -> mu+ (status = 3) mu- (status = 3)
431     // mu- (status = 3) -> mu- (status = 2) -> mu- (status = 1)
432 pakhotin 1.2 const reco::Candidate *genMuonCand = &(*iGenParticle);
433     bool isMuonMother = true;
434     while(isMuonMother) {
435     isMuonMother = false;
436     for ( size_t iMother = 0; iMother < genMuonCand->numberOfMothers(); iMother++ ) {
437     if ( fabs( genMuonCand->mother(iMother)->pdgId() ) == 13 ) {
438     isMuonMother = true;
439     genMuonCand = genMuonCand->mother(iMother);
440     }
441     }
442     }
443 pakhotin 1.3 // Loop over all real (non-muon) mothers of the muon (here we use genMuonCand)
444 pakhotin 1.2 for ( size_t iMother = 0; iMother < genMuonCand->numberOfMothers(); iMother++ ) {
445 pakhotin 1.3 // Check if mother is CP-odd Higgs (PdgId = 36) or gamma_Dark (PdgId = 3000022)
446 pakhotin 1.5 // if ( genMuonCand->mother(iMother)->pdgId() == 36 || genMuonCand->mother(iMother)->pdgId() == 3000022 || genMuonCand->mother(iMother)->pdgId() == 443 ) {
447     if ( genMuonCand->mother(iMother)->pdgId() == 36 || genMuonCand->mother(iMother)->pdgId() == 3000022 ) {
448 pakhotin 1.3 // Store the muon (stable, first in chain) into vector
449 pakhotin 1.2 genMuons.push_back(&(*iGenParticle));
450 pakhotin 1.3 // Store mother of the muon into vector. We need this to group muons into dimuons later
451     genMuonMothers.push_back(genMuonCand->mother(iMother));
452 pakhotin 1.2 }
453     }
454     }
455 pakhotin 1.5 // Check if gen particle is
456     if ( ( iGenParticle->status() == 3 && iGenParticle->pdgId() == 25 ) // decaying (status = 3) SM Higgs (pdgId = 25)
457     || ( iGenParticle->status() == 3 && iGenParticle->pdgId() == 35 ) ) { // decaying (status = 3) CP-even Higgs (pdgId = 35)
458 pakhotin 1.3 genH.push_back(&(*iGenParticle)); // Store the Higgs into vector
459     }
460 pakhotin 1.4 // Check if gen particle is
461     if ( ( iGenParticle->status() == 3 && iGenParticle->pdgId() == 36 ) // decaying (status = 3) CP-odd Higgs (pdgId = 36)
462 pakhotin 1.5 || ( iGenParticle->status() == 3 && iGenParticle->pdgId() == 3000022 ) ) { // decaying (status = 3) gamma_Dark (pdgId = 3000022)
463     // || ( iGenParticle->status() == 2 && iGenParticle->pdgId() == 443 ) ) { // decaying (status = 2) J/psi (pdgId = 443)
464 pakhotin 1.3 genA.push_back(&(*iGenParticle));
465     }
466     }
467    
468     if ( genH.size() == 1 ) {
469 pakhotin 1.5 b_genH_m = genH[0]->mass();
470 pakhotin 1.3 b_genH_p = genH[0]->p();
471     b_genH_pT = genH[0]->pt();
472     b_genH_eta = genH[0]->eta();
473 pakhotin 1.5 b_genH_phi = genH[0]->phi();
474     b_genH_vx = genH[0]->vx() - b_beamSpot_x;
475     b_genH_vy = genH[0]->vy() - b_beamSpot_y;
476     b_genH_vz = genH[0]->vz() - b_beamSpot_z;
477 pakhotin 1.3 } else {
478 pakhotin 1.4 // std::cout << "WARNING! genH.size() != 1" << std::endl;
479 pakhotin 1.3 }
480    
481     if ( genA.size() >= 2 ) {
482     // Sort genA by pT (leading pT first)
483     std::sort (genA.begin(), genA.end(), PtOrder);
484 pakhotin 1.5 b_genA0_m = genA[0]->mass();
485 pakhotin 1.3 b_genA0_p = genA[0]->p();
486     b_genA0_pT = genA[0]->pt();
487     b_genA0_eta = genA[0]->eta();
488 pakhotin 1.5 b_genA0_phi = genA[0]->phi();
489     b_genA0_vx = genA[0]->vx() - b_beamSpot_x;
490     b_genA0_vy = genA[0]->vy() - b_beamSpot_y;
491     b_genA0_vz = genA[0]->vz() - b_beamSpot_z;
492    
493     b_genA1_m = genA[1]->mass();
494 pakhotin 1.3 b_genA1_p = genA[1]->p();
495     b_genA1_pT = genA[1]->pt();
496     b_genA1_eta = genA[1]->eta();
497 pakhotin 1.5 b_genA1_phi = genA[1]->phi();
498     b_genA1_vx = genA[1]->vx() - b_beamSpot_x;
499     b_genA1_vy = genA[1]->vy() - b_beamSpot_y;
500     b_genA1_vz = genA[1]->vz() - b_beamSpot_z;
501 pakhotin 1.3 } else {
502     std::cout << "WARNING! genA.size() < 2" << std::endl;
503     }
504    
505 pakhotin 1.5 // Group muons with the same mother
506     std::vector< std::vector<const reco::GenParticle*> > genMuonGroupsUnsorted;
507     std::vector<const reco::Candidate*> genMuonGroupsUnsortedMothers;
508 pakhotin 1.3 std::vector<const reco::GenParticle*> genMuonsTMP1 = genMuons;
509     std::vector<const reco::Candidate*> genMuonMothersTMP1 = genMuonMothers;
510     unsigned int nMuonGroup = 0;
511     while ( genMuonsTMP1.size() > 0 ) {
512     std::vector<const reco::GenParticle*> genMuonsTMP2;
513     std::vector<const reco::Candidate*> genMuonMothersTMP2;
514     std::vector<const reco::GenParticle*> genMuonsSameMother;
515 pakhotin 1.5 std::vector<const reco::Candidate*> genMuonMothersSame;
516 pakhotin 1.3 for ( unsigned int j = 0; j < genMuonsTMP1.size(); j++ ) {
517     // Check if mothers are the same particle
518 pakhotin 1.5 if ( fabs( genMuonMothersTMP1[0]->pt() - genMuonMothersTMP1[j]->pt() ) < 0.000001 ) {
519 pakhotin 1.3 genMuonsSameMother.push_back( genMuonsTMP1[j] );
520     } else {
521     genMuonsTMP2.push_back( genMuonsTMP1[j] );
522     genMuonMothersTMP2.push_back( genMuonMothersTMP1[j] );
523     }
524     }
525 pakhotin 1.5 genMuonGroupsUnsorted.push_back(genMuonsSameMother);
526     genMuonGroupsUnsortedMothers.push_back(genMuonMothersTMP1[0]);
527 pakhotin 1.3 genMuonsTMP1 = genMuonsTMP2;
528     genMuonMothersTMP1 = genMuonMothersTMP2;
529     nMuonGroup++;
530 pakhotin 1.2 }
531    
532 pakhotin 1.5 // Sort muon groups to match order of genA vector
533     std::vector< std::vector<const reco::GenParticle*> > genMuonGroups;
534     std::vector<const reco::Candidate*> genMuonGroupsMothers;
535     for (unsigned int iA = 0; iA < genA.size(); iA++ ) {
536     bool isMuGroupMatchedToA = false;
537     int nMuGroup = -1;
538     for ( unsigned int iMuGroup = 0; iMuGroup < genMuonGroupsUnsortedMothers.size(); iMuGroup++ ) {
539     if ( fabs ( genA[iA]->pt() - genMuonGroupsUnsortedMothers[iMuGroup]->pt() ) < 0.000001 ) {
540     isMuGroupMatchedToA = true;
541     nMuGroup = iMuGroup;
542     break;
543     }
544     }
545     if ( isMuGroupMatchedToA && nMuGroup >= 0 ) {
546     genMuonGroups.push_back( genMuonGroupsUnsorted[nMuGroup] );
547     genMuonGroupsMothers.push_back( genMuonGroupsUnsortedMothers[nMuGroup] );
548     } else {
549     std::cout << "Error! Muon group has no matched boson A" << std::endl;
550     }
551     }
552    
553 pakhotin 1.10 bool b_isVLT = true;
554 pakhotin 1.3 if ( genMuonGroups.size() == 2 && genMuonGroups[0].size() == 2 && genMuonGroups[1].size() == 2 ) {
555     std::sort( genMuonGroups[0].begin(), genMuonGroups[0].end(), PtOrder );
556     std::sort( genMuonGroups[1].begin(), genMuonGroups[1].end(), PtOrder );
557 pakhotin 1.5
558     b_genA0Mu0_p = genMuonGroups[0][0]->p();
559     b_genA0Mu1_p = genMuonGroups[0][1]->p();
560     b_genA1Mu0_p = genMuonGroups[1][0]->p();
561     b_genA1Mu1_p = genMuonGroups[1][1]->p();
562    
563 pakhotin 1.3 b_genA0Mu0_pT = genMuonGroups[0][0]->pt();
564     b_genA0Mu1_pT = genMuonGroups[0][1]->pt();
565     b_genA1Mu0_pT = genMuonGroups[1][0]->pt();
566     b_genA1Mu1_pT = genMuonGroups[1][1]->pt();
567    
568 pakhotin 1.5 b_genA0Mu0_eta = genMuonGroups[0][0]->eta();
569     b_genA0Mu1_eta = genMuonGroups[0][1]->eta();
570     b_genA1Mu0_eta = genMuonGroups[1][0]->eta();
571     b_genA1Mu1_eta = genMuonGroups[1][1]->eta();
572    
573     b_genA0Mu0_phi = genMuonGroups[0][0]->phi();
574     b_genA0Mu1_phi = genMuonGroups[0][1]->phi();
575     b_genA1Mu0_phi = genMuonGroups[1][0]->phi();
576     b_genA1Mu1_phi = genMuonGroups[1][1]->phi();
577    
578     b_genA0Mu0_vx = genMuonGroups[0][0]->vx() - b_beamSpot_x;
579     b_genA0Mu1_vx = genMuonGroups[0][1]->vx() - b_beamSpot_x;
580     b_genA1Mu0_vx = genMuonGroups[1][0]->vx() - b_beamSpot_x;
581     b_genA1Mu1_vx = genMuonGroups[1][1]->vx() - b_beamSpot_x;
582    
583     b_genA0Mu0_vy = genMuonGroups[0][0]->vy() - b_beamSpot_y;
584     b_genA0Mu1_vy = genMuonGroups[0][1]->vy() - b_beamSpot_y;
585     b_genA1Mu0_vy = genMuonGroups[1][0]->vy() - b_beamSpot_y;
586     b_genA1Mu1_vy = genMuonGroups[1][1]->vy() - b_beamSpot_y;
587    
588     b_genA0Mu0_vz = genMuonGroups[0][0]->vz() - b_beamSpot_z;
589     b_genA0Mu1_vz = genMuonGroups[0][1]->vz() - b_beamSpot_z;
590     b_genA1Mu0_vz = genMuonGroups[1][0]->vz() - b_beamSpot_z;
591     b_genA1Mu1_vz = genMuonGroups[1][1]->vz() - b_beamSpot_z;
592    
593     double eq = 0.000001;
594     if ( fabs(b_genA0Mu0_vx - b_genA0Mu1_vx) < eq
595     && fabs(b_genA1Mu0_vx - b_genA1Mu1_vx) < eq
596     && fabs(b_genA0Mu0_vy - b_genA0Mu1_vy) < eq
597     && fabs(b_genA1Mu0_vy - b_genA1Mu1_vy) < eq
598     && fabs(b_genA0Mu0_vz - b_genA0Mu1_vz) < eq
599     && fabs(b_genA1Mu0_vz - b_genA1Mu1_vz) < eq
600     ) {
601     Float_t genA0_vLx = b_genA0Mu0_vx - b_genA0_vx;
602     Float_t genA1_vLx = b_genA1Mu0_vx - b_genA1_vx;
603    
604     Float_t genA0_vLy = b_genA0Mu0_vy - b_genA0_vy;
605     Float_t genA1_vLy = b_genA1Mu0_vy - b_genA1_vy;
606    
607     Float_t genA0_vLz = b_genA0Mu0_vz - b_genA0_vz;
608     Float_t genA1_vLz = b_genA1Mu0_vz - b_genA1_vz;
609    
610     b_genA0_vLT = sqrt( genA0_vLx * genA0_vLx + genA0_vLy * genA0_vLy );
611     b_genA1_vLT = sqrt( genA1_vLx * genA1_vLx + genA1_vLy * genA1_vLy );
612    
613 pakhotin 1.10 if ( b_genA0_vLT < 4. && b_genA0_vLT < 4. ) b_isVLT = true;
614 pakhotin 1.5
615     b_genA0_vL = sqrt( genA0_vLx * genA0_vLx + genA0_vLy * genA0_vLy + genA0_vLz * genA0_vLz );
616     b_genA1_vL = sqrt( genA1_vLx * genA1_vLx + genA1_vLy * genA1_vLy + genA1_vLz * genA1_vLz );
617     } else {
618     std::cout << "WARNING! Muon vertexes are different" << std::endl;
619     }
620    
621     b_genA0Mu_dEta = genMuonGroups[0][0]->eta() - genMuonGroups[0][1]->eta();
622     b_genA1Mu_dEta = genMuonGroups[1][0]->eta() - genMuonGroups[1][1]->eta();
623     b_genA0Mu_dPhi = My_dPhi( genMuonGroups[0][0]->phi(), genMuonGroups[0][1]->phi() );
624     b_genA1Mu_dPhi = My_dPhi( genMuonGroups[1][0]->phi(), genMuonGroups[1][1]->phi() );
625     b_genA0Mu_dR = sqrt(b_genA0Mu_dEta*b_genA0Mu_dEta + b_genA0Mu_dPhi*b_genA0Mu_dPhi);
626     b_genA1Mu_dR = sqrt(b_genA1Mu_dEta*b_genA1Mu_dEta + b_genA1Mu_dPhi*b_genA1Mu_dPhi);
627 pakhotin 1.3 }
628    
629     // Sort genMuons by pT (leading pT first)
630     if ( genMuons.size() > 1 ) std::sort( genMuons.begin(), genMuons.end(), PtOrder );
631    
632     if ( genMuons.size() == 4 ) m_events4GenMu++;
633    
634 pakhotin 1.6 if ( genMuons.size() > 0 ) {
635     b_genMu0_pT = genMuons[0]->pt();
636     b_genMu0_eta = genMuons[0]->eta();
637     b_genMu0_phi = genMuons[0]->phi();
638     } else {
639     b_genMu0_pT = -100.0;
640     b_genMu0_eta = -100.0;
641     b_genMu0_phi = -100.0;
642     }
643     if ( genMuons.size() > 1 ) {
644     b_genMu1_pT = genMuons[1]->pt();
645     b_genMu1_eta = genMuons[1]->eta();
646     b_genMu1_phi = genMuons[1]->phi();
647     } else {
648     b_genMu1_pT = -100.0;
649     b_genMu1_eta = -100.0;
650     b_genMu1_phi = -100.0;
651 pakhotin 1.2 }
652 pakhotin 1.6 if ( genMuons.size() > 2 ) {
653     b_genMu2_pT = genMuons[2]->pt();
654     b_genMu2_eta = genMuons[2]->eta();
655     b_genMu2_phi = genMuons[2]->phi();
656     } else {
657     b_genMu2_pT = -100.0;
658     b_genMu2_eta = -100.0;
659     b_genMu2_phi = -100.0;
660     }
661     if ( genMuons.size() > 3 ) {
662     b_genMu3_pT = genMuons[3]->pt();
663     b_genMu3_eta = genMuons[3]->eta();
664     b_genMu3_phi = genMuons[3]->phi();
665     } else {
666     b_genMu3_pT = -100.0;
667     b_genMu3_eta = -100.0;
668     b_genMu3_phi = -100.0;
669     }
670    
671 pakhotin 1.2
672     std::vector<const reco::GenParticle*> genMuons17;
673     std::vector<const reco::GenParticle*> genMuons8;
674    
675     for ( unsigned int i = 0; i < genMuons.size(); i++ ) {
676 pakhotin 1.4 if ( genMuons[i]->pt() > m_Mu17_threshold && fabs( genMuons[i]->eta() ) < 0.9 ) {
677 pakhotin 1.2 genMuons17.push_back(genMuons[i]);
678     }
679 pakhotin 1.4 if ( genMuons[i]->pt() > m_Mu8_threshold && fabs( genMuons[i]->eta() ) < 2.4 ) {
680 pakhotin 1.2 genMuons8.push_back(genMuons[i]);
681     }
682     }
683    
684 pakhotin 1.10 // std::cout << " b_isVLT: " << b_isVLT << std::endl;
685     if ( genMuons17.size() >=1 && b_isVLT) {
686 pakhotin 1.3 m_events1GenMu17++;
687 pakhotin 1.2 if ( genMuons8.size() >=2 ) {
688 pakhotin 1.3 m_events2GenMu8++;
689 pakhotin 1.2 }
690     if ( genMuons8.size() >=3 ) {
691 pakhotin 1.3 m_events3GenMu8++;
692 pakhotin 1.2 }
693     if ( genMuons8.size() >=4 ) {
694 pakhotin 1.3 m_events4GenMu8++;
695     }
696     }
697    
698     //****************************************************************************
699     // RECO LEVEL
700     //****************************************************************************
701    
702     edm::Handle<pat::MuonCollection> muons;
703     iEvent.getByLabel(m_muons, muons);
704    
705     std::vector<const reco::Muon*> selMuons;
706     std::vector<const reco::Muon*> selMuons8;
707     std::vector<const reco::Muon*> selMuons17;
708    
709     for (pat::MuonCollection::const_iterator iMuon = muons->begin(); iMuon != muons->end(); ++iMuon) {
710 pakhotin 1.7 if ( isPFMuonLoose( &(*iMuon) ) ) {
711     // if ( isTrackerMuonPrivateID( &(*iMuon) ) ) {
712 pakhotin 1.3 selMuons.push_back( &(*iMuon) );
713 pakhotin 1.4 if ( iMuon->pt() > m_Mu8_threshold ) {
714 pakhotin 1.3 selMuons8.push_back( &(*iMuon) );
715     }
716 pakhotin 1.4 if ( iMuon->pt() > m_Mu17_threshold && fabs(iMuon->eta()) < 0.9 ) {
717 pakhotin 1.3 selMuons17.push_back( &(*iMuon) );
718     }
719     }
720     }
721    
722 pakhotin 1.6 if ( selMuons.size() > 0 ) {
723     b_selMu0_pT = selMuons[0]->pt();
724     b_selMu0_eta = selMuons[0]->eta();
725     b_selMu0_phi = selMuons[0]->phi();
726     } else {
727     b_selMu0_pT = -100.0;
728     b_selMu0_eta = -100.0;
729     b_selMu0_phi = -100.0;
730     }
731     if ( selMuons.size() > 1 ) {
732     b_selMu1_pT = selMuons[1]->pt();
733     b_selMu1_eta = selMuons[1]->eta();
734     b_selMu1_phi = selMuons[1]->phi();
735     } else {
736     b_selMu1_pT = -100.0;
737     b_selMu1_eta = -100.0;
738     b_selMu1_phi = -100.0;
739     }
740     if ( selMuons.size() > 2 ) {
741     b_selMu2_pT = selMuons[2]->pt();
742     b_selMu2_eta = selMuons[2]->eta();
743     b_selMu2_phi = selMuons[2]->phi();
744     } else {
745     b_selMu2_pT = -100.0;
746     b_selMu2_eta = -100.0;
747     b_selMu2_phi = -100.0;
748     }
749     if ( selMuons.size() > 3 ) {
750     b_selMu3_pT = selMuons[3]->pt();
751     b_selMu3_eta = selMuons[3]->eta();
752     b_selMu3_phi = selMuons[3]->phi();
753     } else {
754     b_selMu3_pT = -100.0;
755     b_selMu3_eta = -100.0;
756     b_selMu3_phi = -100.0;
757     }
758    
759    
760 pakhotin 1.10 b_is1SelMu17 = false;
761     b_is4SelMu8 = false;
762 pakhotin 1.3 if ( selMuons17.size() >=1 ) {
763     m_events1SelMu17++;
764 pakhotin 1.10 b_is1SelMu17 = true;
765 pakhotin 1.3 if ( selMuons8.size() >=2 ) {
766     m_events2SelMu8++;
767     }
768     if ( selMuons8.size() >=3 ) {
769     m_events3SelMu8++;
770 pakhotin 1.2 }
771 pakhotin 1.3 if ( selMuons8.size() >=4 ) {
772     m_events4SelMu8++;
773 pakhotin 1.10 b_is4SelMu8 = true;
774 pakhotin 1.3 }
775     }
776    
777     edm::Handle<pat::MultiMuonCollection> muJets;
778     iEvent.getByLabel(m_muJets, muJets);
779     const pat::MultiMuon *muJetC = NULL;
780     const pat::MultiMuon *muJetF = NULL;
781     int nMuJetsContainMu17 = 0;
782     unsigned int nMuJets = muJets->size();
783 pakhotin 1.10 b_is2MuJets = false;
784 pakhotin 1.3 if ( nMuJets == 2) {
785     for ( unsigned int j = 0; j < nMuJets; j++ ) {
786     bool isMuJetContainMu17 = false;
787     for ( unsigned int m = 0; m < (*muJets)[j].numberOfDaughters(); m++ ) {
788 pakhotin 1.4 if ( (*muJets)[j].muon(m)->pt() > m_Mu17_threshold && fabs( (*muJets)[j].muon(m)->eta() ) < 0.9 ) {
789 pakhotin 1.3 isMuJetContainMu17 = true;
790     nMuJetsContainMu17++;
791     break;
792     }
793     }
794     if ( isMuJetContainMu17 ) {
795     muJetC = &((*muJets)[j]);
796     } else {
797     muJetF = &((*muJets)[j]);
798     }
799     }
800     if ( nMuJetsContainMu17 == 2 ) {
801     if (m_trandom3.Integer(2) == 0) {
802     muJetC = &((*muJets)[0]);
803     muJetF = &((*muJets)[1]);
804     } else {
805     muJetC = &((*muJets)[1]);
806     muJetF = &((*muJets)[0]);
807     }
808     }
809 pakhotin 1.10 if ( nMuJetsContainMu17 > 0 ) b_is2MuJets = true;
810 pakhotin 1.3 }
811    
812 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_isVLT) m_events2MuJets++;
813 pakhotin 1.3
814 pakhotin 1.10 b_is2DiMuons = false;
815 pakhotin 1.3 const pat::MultiMuon *diMuonC = NULL;
816     const pat::MultiMuon *diMuonF = NULL;
817     if ( muJetC != NULL && muJetF != NULL ) {
818     if ( muJetC->numberOfDaughters() == 2 && muJetF->numberOfDaughters() == 2 ) {
819     diMuonC = muJetC;
820     diMuonF = muJetF;
821 pakhotin 1.10 b_is2DiMuons = true;
822 pakhotin 1.3 }
823     }
824    
825 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_is2DiMuons && b_isVLT) m_events2DiMuons++;
826 pakhotin 1.3
827 pakhotin 1.10 b_isDzDiMuonsOK = false;
828     b_dZDiMuonC = -1000.0;
829     b_dZDiMuonF = -1000.0;
830     b_dZDiMuons = -1000.0;
831 pakhotin 1.3 if ( diMuonC != NULL && diMuonF != NULL ) {
832 pakhotin 1.10 b_dZDiMuonC = diMuonC->dz(beamSpot->position());
833     b_dZDiMuonF = diMuonF->dz(beamSpot->position());
834     b_dZDiMuons = b_dZDiMuonC - b_dZDiMuonF;
835     if ( fabs( b_dZDiMuons ) < 0.1 ) b_isDzDiMuonsOK = true;
836 pakhotin 1.3 }
837    
838 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_is2DiMuons && b_isDzDiMuonsOK && b_isVLT ) m_eventsDz2DiMuonsOK++;
839 pakhotin 1.3
840     edm::Handle<pat::TriggerEvent> triggerEvent;
841     iEvent.getByLabel("patTriggerEvent", triggerEvent);
842    
843 pakhotin 1.10 b_isDiMuonHLTFired = false;
844 pakhotin 1.3 if ( ( triggerEvent->path("HLT_Mu17_Mu8_v22") && triggerEvent->path("HLT_Mu17_Mu8_v22")->wasAccept() )
845     || ( triggerEvent->path("HLT_Mu17_Mu8_v21") && triggerEvent->path("HLT_Mu17_Mu8_v21")->wasAccept() )
846     || ( triggerEvent->path("HLT_Mu17_Mu8_v19") && triggerEvent->path("HLT_Mu17_Mu8_v19")->wasAccept() )
847     || ( triggerEvent->path("HLT_Mu17_Mu8_v18") && triggerEvent->path("HLT_Mu17_Mu8_v18")->wasAccept() )
848     || ( triggerEvent->path("HLT_Mu17_Mu8_v17") && triggerEvent->path("HLT_Mu17_Mu8_v17")->wasAccept() )
849     || ( triggerEvent->path("HLT_Mu17_Mu8_v16") && triggerEvent->path("HLT_Mu17_Mu8_v16")->wasAccept() ) ) {
850 pakhotin 1.10 b_isDiMuonHLTFired = true;
851 pakhotin 1.2 }
852    
853 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_is2DiMuons && b_isDzDiMuonsOK && b_isDiMuonHLTFired && b_isVLT ) m_eventsDiMuonsHLTFired++;
854 pakhotin 1.3
855 pakhotin 1.10 b_isMassDiMuonsOK = false;
856 pakhotin 1.3 if ( diMuonC != NULL && diMuonF != NULL ) {
857     double massC = diMuonC->mass();
858     double massF = diMuonF->mass();
859 pakhotin 1.10 if ( fabs(massC-massF) < (0.13 + 0.065*(massC+massF)/2.0) ) b_isMassDiMuonsOK = true;
860 pakhotin 1.3 }
861    
862 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_is2DiMuons && b_isDzDiMuonsOK && b_isDiMuonHLTFired && b_isMassDiMuonsOK && b_isVLT ) m_eventsMassDiMuonsOK++;
863 pakhotin 1.3
864 pakhotin 1.10 edm::Handle<reco::TrackCollection> tracks;
865     iEvent.getByLabel("generalTracks", tracks);
866 pakhotin 1.9
867     edm::Handle<reco::PFCandidateCollection> pfCandidates;
868     iEvent.getByLabel("particleFlow", pfCandidates);
869    
870 pakhotin 1.10 b_isIsoDiMuonsOK = false;
871     b_diMuonC_isoTk = -1.;
872     b_diMuonF_isoTk = -1.;
873     b_diMuonC_isoPF = -1.;
874     b_diMuonF_isoPF = -1.;
875 pakhotin 1.3 if ( diMuonC != NULL && diMuonF != NULL ) {
876 pakhotin 1.10 double diMuonC_isoTk = 0.0;
877     double diMuonF_isoTk = 0.0;
878     double diMuonC_isoPF = 0.0;
879     double diMuonF_isoPF = 0.0;
880 pakhotin 1.3 const pat::MultiMuon *diMuonTmp = NULL;
881     for ( unsigned int i = 1; i <= 2; i++ ) {
882 pakhotin 1.10 double diMuonTmp_isoTk = 0.0;
883     double diMuonTmp_isoPF = 0.0;
884 pakhotin 1.3 if ( i == 1 ) diMuonTmp = diMuonC;
885     if ( i == 2 ) diMuonTmp = diMuonF;
886     for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
887     bool trackIsMuon = false;
888     if ( diMuonTmp->sameTrack( &*track, &*(diMuonTmp->muon(0)->innerTrack()) ) || diMuonTmp->sameTrack( &*track, &*(diMuonTmp->muon(1)->innerTrack()) ) ) trackIsMuon = true;
889 pakhotin 1.10 if ( trackIsMuon == false ) {
890 pakhotin 1.3 double dPhi = My_dPhi( diMuonTmp->phi(), track->phi() );
891     double dEta = diMuonTmp->eta() - track->eta();
892 pakhotin 1.10 double dR = sqrt( dPhi*dPhi + dEta*dEta );
893 pakhotin 1.3 if ( dR < 0.4 && track->pt() > 0.5 ) {
894     double dz = fabs( track->dz(beamSpot->position()) - diMuonTmp->dz(beamSpot->position()) );
895 pakhotin 1.10 if ( dz < 0.1 ) diMuonTmp_isoTk += track->pt();
896 pakhotin 1.3 }
897     }
898     }
899 pakhotin 1.10 for (reco::PFCandidateCollection::const_iterator pfCand = pfCandidates->begin(); pfCand != pfCandidates->end(); ++pfCand) {
900     if ( pfCand->particleId() == reco::PFCandidate::h ) {
901     double dPhi = My_dPhi( diMuonTmp->phi(), pfCand->phi() );
902     double dEta = diMuonTmp->eta() - pfCand->eta();
903     double dR = sqrt( dPhi*dPhi + dEta*dEta );
904     if ( dR < 0.4 && pfCand->pt() > 0.5 ) {
905     const math::XYZPoint dummyBeamSpot = math::XYZPoint(0.0, 0.0, 0.0);
906     double dz = fabs( pfCand->vertex().z() - diMuonTmp->dz(dummyBeamSpot) );
907     if ( dz < 0.1 ) diMuonTmp_isoPF += pfCand->pt();
908     }
909     }
910     }
911     if ( i == 1 ) {
912     diMuonC_isoTk = diMuonTmp_isoTk;
913     diMuonC_isoPF = diMuonTmp_isoPF;
914     }
915     if ( i == 2 ) {
916     diMuonF_isoTk = diMuonTmp_isoTk;
917     diMuonF_isoPF = diMuonTmp_isoPF;
918     }
919 pakhotin 1.3 }
920 pakhotin 1.10 b_diMuonC_isoTk = diMuonC_isoTk;
921     b_diMuonF_isoTk = diMuonF_isoTk;
922     b_diMuonC_isoPF = diMuonC_isoPF;
923     b_diMuonF_isoPF = diMuonF_isoPF;
924     std::cout << "diMuonC_isoTk: " << diMuonC_isoTk << " diMuonC_isoPF: " << diMuonC_isoPF << " diMuonF_isoTk: " << diMuonF_isoTk << " diMuonF_isoPF: " << diMuonF_isoPF << std::endl;
925     if ( diMuonC_isoTk < m_maxIsoDiMuons && diMuonF_isoTk < m_maxIsoDiMuons ) b_isIsoDiMuonsOK = true;
926 pakhotin 1.3 }
927    
928 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_is2DiMuons && b_isDzDiMuonsOK && b_isDiMuonHLTFired && b_isMassDiMuonsOK && b_isIsoDiMuonsOK && b_isVLT ) m_eventsIsoDiMuonsOK++;
929 pakhotin 1.3
930     edm::Handle<reco::VertexCollection> primaryVertices;
931     iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
932    
933     bool isVertexOK = false;
934     for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin(); vertex != primaryVertices->end(); ++vertex) {
935 pakhotin 1.10 if (vertex->isValid() && !vertex->isFake() && vertex->tracksSize() >= 4 && fabs(vertex->z()) < 24.) isVertexOK = true;
936 pakhotin 1.3 }
937    
938 pakhotin 1.10 if ( b_is1SelMu17 && b_is4SelMu8 && b_is2MuJets && b_is2DiMuons && b_isDzDiMuonsOK && b_isDiMuonHLTFired && b_isMassDiMuonsOK && b_isIsoDiMuonsOK && isVertexOK && b_isVLT ) m_eventsVertexOK++;
939 pakhotin 1.3
940     m_ttree->Fill();
941    
942 pakhotin 1.1 }
943    
944    
945     // ------------ method called once each job just before starting event loop ------------
946     void
947     CutFlowAnalyzer::beginJob()
948     {
949 pakhotin 1.2 std::cout << "BEGIN JOB" << std::endl;
950 pakhotin 1.3
951     edm::Service<TFileService> tFileService;
952     m_ttree = tFileService->make<TTree>("Events", "Events");
953 pakhotin 1.5 // Beam spot
954     m_ttree->Branch("beamSpot_x", &b_beamSpot_x, "beamSpot_x/F");
955     m_ttree->Branch("beamSpot_y", &b_beamSpot_y, "beamSpot_y/F");
956     m_ttree->Branch("beamSpot_z", &b_beamSpot_z, "beamSpot_z/F");
957 pakhotin 1.3 // Bosons
958 pakhotin 1.5 m_ttree->Branch("genH_m", &b_genH_m, "genH_m/F");
959 pakhotin 1.3 m_ttree->Branch("genH_p", &b_genH_p, "genH_p/F");
960     m_ttree->Branch("genH_pT", &b_genH_pT, "genH_pT/F");
961 pakhotin 1.5 m_ttree->Branch("genH_eta", &b_genH_eta, "genH_eta/F");
962     m_ttree->Branch("genH_phi", &b_genH_phi, "genH_phi/F");
963 pakhotin 1.4 m_ttree->Branch("genH_vx", &b_genH_vx, "genH_vx/F");
964     m_ttree->Branch("genH_vy", &b_genH_vy, "genH_vy/F");
965     m_ttree->Branch("genH_vz", &b_genH_vz, "genH_vz/F");
966    
967 pakhotin 1.5 m_ttree->Branch("genA0_m", &b_genA0_m, "genA0_m/F");
968 pakhotin 1.3 m_ttree->Branch("genA0_p", &b_genA0_p, "genA0_p/F");
969     m_ttree->Branch("genA0_pT", &b_genA0_pT, "genA0_pT/F");
970     m_ttree->Branch("genA0_eta", &b_genA0_eta, "genA0_eta/F");
971 pakhotin 1.5 m_ttree->Branch("genA0_phi", &b_genA0_phi, "genA0_phi/F");
972     m_ttree->Branch("genA0_vx", &b_genA0_vx, "genA0_vx/F");
973     m_ttree->Branch("genA0_vy", &b_genA0_vy, "genA0_vy/F");
974     m_ttree->Branch("genA0_vz", &b_genA0_vz, "genA0_vz/F");
975    
976     m_ttree->Branch("genA1_m", &b_genA1_m, "genA1_m/F");
977 pakhotin 1.3 m_ttree->Branch("genA1_p", &b_genA1_p, "genA1_p/F");
978     m_ttree->Branch("genA1_pT", &b_genA1_pT, "genA1_pT/F");
979     m_ttree->Branch("genA1_eta", &b_genA1_eta, "genA1_eta/F");
980 pakhotin 1.5 m_ttree->Branch("genA1_phi", &b_genA1_phi, "genA1_phi/F");
981     m_ttree->Branch("genA1_vx", &b_genA1_vx, "genA1_vx/F");
982     m_ttree->Branch("genA1_vy", &b_genA1_vy, "genA1_vy/F");
983     m_ttree->Branch("genA1_vz", &b_genA1_vz, "genA1_vz/F");
984    
985 pakhotin 1.3 // Dimuons
986 pakhotin 1.5 m_ttree->Branch("genA0Mu0_p", &b_genA0Mu0_p, "genA0Mu0_p/F");
987     m_ttree->Branch("genA0Mu1_p", &b_genA0Mu1_p, "genA0Mu1_p/F");
988     m_ttree->Branch("genA1Mu0_p", &b_genA1Mu0_p, "genA1Mu0_p/F");
989     m_ttree->Branch("genA1Mu1_p", &b_genA1Mu1_p, "genA1Mu1_p/F");
990    
991 pakhotin 1.3 m_ttree->Branch("genA0Mu0_pT", &b_genA0Mu0_pT, "genA0Mu0_pT/F");
992     m_ttree->Branch("genA0Mu1_pT", &b_genA0Mu1_pT, "genA0Mu1_pT/F");
993     m_ttree->Branch("genA1Mu0_pT", &b_genA1Mu0_pT, "genA1Mu0_pT/F");
994     m_ttree->Branch("genA1Mu1_pT", &b_genA1Mu1_pT, "genA1Mu1_pT/F");
995 pakhotin 1.4
996 pakhotin 1.5 m_ttree->Branch("genA0Mu0_eta", &b_genA0Mu0_eta, "genA0Mu0_eta/F");
997     m_ttree->Branch("genA0Mu1_eta", &b_genA0Mu1_eta, "genA0Mu1_eta/F");
998     m_ttree->Branch("genA1Mu0_eta", &b_genA1Mu0_eta, "genA1Mu0_eta/F");
999     m_ttree->Branch("genA1Mu1_eta", &b_genA1Mu1_eta, "genA1Mu1_eta/F");
1000    
1001     m_ttree->Branch("genA0Mu0_phi", &b_genA0Mu0_phi, "genA0Mu0_phi/F");
1002     m_ttree->Branch("genA0Mu1_phi", &b_genA0Mu1_phi, "genA0Mu1_phi/F");
1003     m_ttree->Branch("genA1Mu0_phi", &b_genA1Mu0_phi, "genA1Mu0_phi/F");
1004     m_ttree->Branch("genA1Mu1_phi", &b_genA1Mu1_phi, "genA1Mu1_phi/F");
1005    
1006 pakhotin 1.4 m_ttree->Branch("genA0Mu0_vx", &b_genA0Mu0_vx, "genA0Mu0_vx/F");
1007     m_ttree->Branch("genA0Mu1_vx", &b_genA0Mu1_vx, "genA0Mu1_vx/F");
1008     m_ttree->Branch("genA1Mu0_vx", &b_genA1Mu0_vx, "genA1Mu0_vx/F");
1009     m_ttree->Branch("genA1Mu1_vx", &b_genA1Mu1_vx, "genA1Mu1_vx/F");
1010    
1011     m_ttree->Branch("genA0Mu0_vy", &b_genA0Mu0_vy, "genA0Mu0_vy/F");
1012     m_ttree->Branch("genA0Mu1_vy", &b_genA0Mu1_vy, "genA0Mu1_vy/F");
1013     m_ttree->Branch("genA1Mu0_vy", &b_genA1Mu0_vy, "genA1Mu0_vy/F");
1014     m_ttree->Branch("genA1Mu1_vy", &b_genA1Mu1_vy, "genA1Mu1_vy/F");
1015    
1016     m_ttree->Branch("genA0Mu0_vz", &b_genA0Mu0_vz, "genA0Mu0_vz/F");
1017     m_ttree->Branch("genA0Mu1_vz", &b_genA0Mu1_vz, "genA0Mu1_vz/F");
1018     m_ttree->Branch("genA1Mu0_vz", &b_genA1Mu0_vz, "genA1Mu0_vz/F");
1019     m_ttree->Branch("genA1Mu1_vz", &b_genA1Mu1_vz, "genA1Mu1_vz/F");
1020    
1021 pakhotin 1.5 m_ttree->Branch("genA0_vLT", &b_genA0_vLT, "genA0_vLT/F");
1022     m_ttree->Branch("genA1_vLT", &b_genA1_vLT, "genA1_vLT/F");
1023     m_ttree->Branch("genA0_vL", &b_genA0_vL, "genA0_vL/F");
1024     m_ttree->Branch("genA1_vL", &b_genA1_vL, "genA1_vL/F");
1025    
1026     m_ttree->Branch("genA0Mu_dEta", &b_genA0Mu_dEta, "genA0Mu_dEta/F");
1027     m_ttree->Branch("genA1Mu_dEta", &b_genA1Mu_dEta, "genA1Mu_dEta/F");
1028     m_ttree->Branch("genA0Mu_dPhi", &b_genA0Mu_dPhi, "genA0Mu_dPhi/F");
1029     m_ttree->Branch("genA1Mu_dPhi", &b_genA1Mu_dPhi, "genA1Mu_dPhi/F");
1030     m_ttree->Branch("genA0Mu_dR", &b_genA0Mu_dR, "genA0Mu_dR/F");
1031     m_ttree->Branch("genA1Mu_dR", &b_genA1Mu_dR, "genA1Mu_dR/F");
1032 pakhotin 1.3
1033     // Generator Muons
1034 pakhotin 1.6 m_ttree->Branch("genMu0_pT", &b_genMu0_pT, "genMu0_pT/F");
1035     m_ttree->Branch("genMu1_pT", &b_genMu1_pT, "genMu1_pT/F");
1036     m_ttree->Branch("genMu2_pT", &b_genMu2_pT, "genMu2_pT/F");
1037     m_ttree->Branch("genMu3_pT", &b_genMu3_pT, "genMu3_pT/F");
1038     m_ttree->Branch("genMu0_eta", &b_genMu0_eta, "genMu0_eta/F");
1039     m_ttree->Branch("genMu1_eta", &b_genMu1_eta, "genMu1_eta/F");
1040     m_ttree->Branch("genMu2_eta", &b_genMu2_eta, "genMu2_eta/F");
1041     m_ttree->Branch("genMu3_eta", &b_genMu3_eta, "genMu3_eta/F");
1042     m_ttree->Branch("genMu0_phi", &b_genMu0_phi, "genMu0_phi/F");
1043     m_ttree->Branch("genMu1_phi", &b_genMu1_phi, "genMu1_phi/F");
1044     m_ttree->Branch("genMu2_phi", &b_genMu2_phi, "genMu2_phi/F");
1045     m_ttree->Branch("genMu3_phi", &b_genMu3_phi, "genMu3_phi/F");
1046    
1047     // Reco Muons
1048     m_ttree->Branch("selMu0_pT", &b_selMu0_pT, "selMu0_pT/F");
1049     m_ttree->Branch("selMu1_pT", &b_selMu1_pT, "selMu1_pT/F");
1050     m_ttree->Branch("selMu2_pT", &b_selMu2_pT, "selMu2_pT/F");
1051     m_ttree->Branch("selMu3_pT", &b_selMu3_pT, "selMu3_pT/F");
1052     m_ttree->Branch("selMu0_eta", &b_selMu0_eta, "selMu0_eta/F");
1053     m_ttree->Branch("selMu1_eta", &b_selMu1_eta, "selMu1_eta/F");
1054     m_ttree->Branch("selMu2_eta", &b_selMu2_eta, "selMu2_eta/F");
1055     m_ttree->Branch("selMu3_eta", &b_selMu3_eta, "selMu3_eta/F");
1056     m_ttree->Branch("selMu0_phi", &b_selMu0_phi, "selMu0_phi/F");
1057     m_ttree->Branch("selMu1_phi", &b_selMu1_phi, "selMu1_phi/F");
1058     m_ttree->Branch("selMu2_phi", &b_selMu2_phi, "selMu2_phi/F");
1059     m_ttree->Branch("selMu3_phi", &b_selMu3_phi, "selMu3_phi/F");
1060 pakhotin 1.3
1061 pakhotin 1.10 m_ttree->Branch("dZDiMuonC", &b_dZDiMuonC, "dZDiMuonC/F");
1062     m_ttree->Branch("dZDiMuonF", &b_dZDiMuonF, "dZDiMuonF/F");
1063     m_ttree->Branch("dZDiMuons", &b_dZDiMuons, "dZDiMuons/F");
1064    
1065     m_ttree->Branch("diMuonC_isoTk", &b_diMuonC_isoTk, "diMuonC_isoTk/F");
1066     m_ttree->Branch("diMuonF_isoTk", &b_diMuonF_isoTk, "diMuonF_isoTk/F");
1067     m_ttree->Branch("diMuonC_isoPF", &b_diMuonC_isoPF, "diMuonC_isoPF/F");
1068     m_ttree->Branch("diMuonF_isoPF", &b_diMuonF_isoPF, "diMuonF_isoPF/F");
1069    
1070     m_ttree->Branch("is1SelMu17", &b_is1SelMu17, "is1SelMu17/O");
1071     m_ttree->Branch("is4SelMu8", &b_is4SelMu8, "is4SelMu8/O");
1072     m_ttree->Branch("is2MuJets", &b_is2MuJets, "is2MuJets/O");
1073     m_ttree->Branch("is2DiMuons", &b_is2DiMuons, "is2DiMuons/O");
1074     m_ttree->Branch("isDzDiMuonsOK", &b_isDzDiMuonsOK, "isDzDiMuonsOK/O");
1075     m_ttree->Branch("isDiMuonHLTFired", &b_isDiMuonHLTFired, "isDiMuonHLTFired/O");
1076     m_ttree->Branch("isMassDiMuonsOK", &b_isMassDiMuonsOK, "isMassDiMuonsOK/O");
1077     m_ttree->Branch("isIsoDiMuonsOK", &b_isIsoDiMuonsOK, "isIsoDiMuonsOK/O");
1078     m_ttree->Branch("isVLT", &b_isVLT, "isVLT/O");
1079    
1080 pakhotin 1.1 }
1081    
1082     // ------------ method called once each job just after ending the event loop ------------
1083     void
1084     CutFlowAnalyzer::endJob()
1085     {
1086 pakhotin 1.2 std::cout << "END JOB" << std::endl;
1087    
1088 pakhotin 1.3 std:: cout << "Total number of events: " << m_events << std::endl;
1089     std:: cout << "Total number of events with 4mu: " << m_events4GenMu << " fraction: " << m_events4GenMu/m_events << std::endl;
1090    
1091     std:: cout << "********** GEN **********" << std::endl;
1092     std:: cout << "Selection " << "nEv" << " \t RelEff" << " \t Eff" << std::endl;
1093     std:: cout << "pT1>17 |eta1|<0.9: " << m_events1GenMu17 << " \t" << (float)m_events1GenMu17/(float)m_events << " \t" << (float)m_events1GenMu17/(float)m_events << std::endl;
1094     std:: cout << "pT2>8 |eta2|<2.4: " << m_events2GenMu8 << " \t" << (float)m_events2GenMu8/(float)m_events1GenMu17 << " \t" << (float)m_events2GenMu8/(float)m_events << std::endl;
1095     std:: cout << "pT3>8 |eta2|<2.4: " << m_events3GenMu8 << " \t" << (float)m_events3GenMu8/(float)m_events2GenMu8 << " \t" << (float)m_events3GenMu8/(float)m_events << std::endl;
1096     std:: cout << "pT4>8 |eta2|<2.4: " << m_events4GenMu8 << " \t" << (float)m_events4GenMu8/(float)m_events3GenMu8 << " \t" << (float)m_events4GenMu8/(float)m_events << std::endl;
1097     std:: cout << "Basic MC Acceptance: " << (float)m_events4GenMu8/(float)m_events << std::endl;
1098    
1099     std:: cout << "********** RECO **********" << std::endl;
1100     std:: cout << "Selection " << "nEv" << " \t RelEff" << " \t Eff" << std::endl;
1101     std:: cout << "m_events1SelMu17: " << m_events1SelMu17 << " \t" << (float)m_events1SelMu17/(float)m_events << " \t" << (float)m_events1SelMu17/(float)m_events << std::endl;
1102     std:: cout << "m_events2SelMu8: " << m_events2SelMu8 << " \t" << (float)m_events2SelMu8/(float)m_events1SelMu17 << " \t" << (float)m_events2SelMu8/(float)m_events << std::endl;
1103     std:: cout << "m_events3SelMu8: " << m_events3SelMu8 << " \t" << (float)m_events3SelMu8/(float)m_events2SelMu8 << " \t" << (float)m_events3SelMu8/(float)m_events << std::endl;
1104     std:: cout << "m_events4SelMu8: " << m_events4SelMu8 << " \t" << (float)m_events4SelMu8/(float)m_events3SelMu8 << " \t" << (float)m_events4SelMu8/(float)m_events << std::endl;
1105    
1106     std:: cout << "Basic Acceptance: " << (float)m_events4SelMu8/(float)m_events << std::endl;
1107     std:: cout << "Basic MC Accept. a_gen: " << (float)m_events4GenMu8/(float)m_events << std::endl;
1108    
1109     std:: cout << "m_events2MuJets: " << m_events2MuJets << " \t" << (float)m_events2MuJets/(float)m_events4SelMu8 << " \t" << (float)m_events2MuJets/(float)m_events << std::endl;
1110     std:: cout << "m_events2DiMuons: " << m_events2DiMuons << " \t" << (float)m_events2DiMuons/(float)m_events2MuJets << " \t" << (float)m_events2DiMuons/(float)m_events << std::endl;
1111     std:: cout << "m_eventsDz2DiMuonsOK: " << m_eventsDz2DiMuonsOK << " \t" << (float)m_eventsDz2DiMuonsOK/(float)m_events2DiMuons << " \t" << (float)m_eventsDz2DiMuonsOK/(float)m_events << std::endl;
1112     std:: cout << "m_eventsDiMuonsHLTFired: " << m_eventsDiMuonsHLTFired << " \t" << (float)m_eventsDiMuonsHLTFired/(float)m_eventsDz2DiMuonsOK << " \t" << (float)m_eventsDiMuonsHLTFired/(float)m_events << std::endl;
1113     std:: cout << "m_eventsMassDiMuonsOK: " << m_eventsMassDiMuonsOK << " \t" << (float)m_eventsMassDiMuonsOK/(float)m_eventsDiMuonsHLTFired << " \t" << (float)m_eventsMassDiMuonsOK/(float)m_events << std::endl;
1114     std:: cout << "m_eventsIsoDiMuonsOK: " << m_eventsIsoDiMuonsOK << " \t" << (float)m_eventsIsoDiMuonsOK/(float)m_eventsMassDiMuonsOK << " \t" << (float)m_eventsIsoDiMuonsOK/(float)m_events << std::endl;
1115     std:: cout << "m_eventsVertexOK: " << m_eventsVertexOK << " \t" << (float)m_eventsVertexOK/(float)m_eventsIsoDiMuonsOK << " \t" << (float)m_eventsVertexOK/(float)m_events << std::endl;
1116    
1117     std:: cout << "Further selections: " << (float)m_eventsVertexOK/(float)m_events4SelMu8 << std::endl;
1118     std:: cout << "Full sel eff e_full: " << (float)m_eventsVertexOK/(float)m_events << std::endl;
1119     std:: cout << "e_full/a_gen: " << (float)m_eventsVertexOK/(float)m_events4GenMu8 << std::endl;
1120    
1121     std::cout << m_events << std::endl;
1122     std::cout << m_events1GenMu17 << std::endl;
1123     std::cout << m_events2GenMu8 << std::endl;
1124     std::cout << m_events3GenMu8 << std::endl;
1125     std::cout << m_events4GenMu8 << std::endl;
1126     std::cout << m_events1SelMu17 << std::endl;
1127     std::cout << m_events2SelMu8 << std::endl;
1128     std::cout << m_events3SelMu8 << std::endl;
1129     std::cout << m_events4SelMu8 << std::endl;
1130     std::cout << m_events2MuJets << std::endl;
1131     std::cout << m_events2DiMuons << std::endl;
1132     std::cout << m_eventsDz2DiMuonsOK << std::endl;
1133     std::cout << m_eventsDiMuonsHLTFired << std::endl;
1134     std::cout << m_eventsMassDiMuonsOK << std::endl;
1135     std::cout << m_eventsIsoDiMuonsOK << std::endl;
1136     std::cout << m_eventsVertexOK << std::endl;
1137    
1138    
1139 pakhotin 1.1 }
1140    
1141     // ------------ method called when starting to processes a run ------------
1142     void
1143     CutFlowAnalyzer::beginRun(edm::Run const&, edm::EventSetup const&)
1144     {
1145     }
1146    
1147     // ------------ method called when ending the processing of a run ------------
1148     void
1149     CutFlowAnalyzer::endRun(edm::Run const&, edm::EventSetup const&)
1150     {
1151     }
1152    
1153     // ------------ method called when starting to processes a luminosity block ------------
1154     void
1155     CutFlowAnalyzer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1156     {
1157     }
1158    
1159     // ------------ method called when ending the processing of a luminosity block ------------
1160     void
1161     CutFlowAnalyzer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
1162     {
1163     }
1164    
1165     // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
1166     void
1167     CutFlowAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1168     //The following says we do not know what parameters are allowed so do no validation
1169     // Please change this to state exactly what you do use, even if it is no parameters
1170     edm::ParameterSetDescription desc;
1171     desc.setUnknown();
1172     descriptions.addDefault(desc);
1173     }
1174    
1175     //define this as a plug-in
1176     DEFINE_FWK_MODULE(CutFlowAnalyzer);