ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/CutFlowAnalyzer/src/CutFlowAnalyzer.cc
Revision: 1.6
Committed: Fri Apr 26 08:17:35 2013 UTC (12 years ago) by pakhotin
Content type: text/plain
Branch: MAIN
Changes since 1.5: +127 -10 lines
Log Message:
*** empty log message ***

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