ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/CutFlowAnalyzer/src/CutFlowAnalyzer.cc
Revision: 1.7
Committed: Fri Apr 26 20:54:57 2013 UTC (12 years ago) by pakhotin
Content type: text/plain
Branch: MAIN
Changes since 1.6: +27 -12 lines
Log Message:
version used for note

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