ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/FitNtuple/src/FitNtuple.cc
Revision: 1.6
Committed: Thu Nov 25 16:42:05 2010 UTC (14 years, 5 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-11-25-a
Changes since 1.5: +114 -6 lines
Log Message:
new Drell-Yan sample (not fully propagated through yet)

File Contents

# User Rev Content
1 pivarski 1.1 // -*- C++ -*-
2     //
3     // Package: FitNtuple
4     // Class: FitNtuple
5     //
6     /**\class FitNtuple FitNtuple.cc AnalysisTools/FitNtuple/src/FitNtuple.cc
7    
8     Description: [one line class summary]
9    
10     Implementation:
11     [Notes on implementation]
12     */
13     //
14     // Original Author: James Pivarski
15     // Created: Fri Nov 19 17:04:33 CST 2010
16 pivarski 1.6 // $Id: FitNtuple.cc,v 1.5 2010/11/22 21:17:43 pivarski Exp $
17 pivarski 1.1 //
18     //
19    
20    
21     // system include files
22     #include <memory>
23    
24     // user include files
25     #include "FWCore/Framework/interface/Frameworkfwd.h"
26     #include "FWCore/Framework/interface/EDAnalyzer.h"
27     #include "FWCore/Framework/interface/Event.h"
28     #include "FWCore/Framework/interface/MakerMacros.h"
29     #include "FWCore/ParameterSet/interface/ParameterSet.h"
30     #include "FWCore/Utilities/interface/InputTag.h"
31    
32     #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
33     #include "DataFormats/Candidate/interface/Candidate.h"
34     #include "AnalysisDataFormats/MuJetAnalysis/interface/MultiMuon.h"
35     #include "DataFormats/PatCandidates/interface/Muon.h"
36     #include "DataFormats/PatCandidates/interface/TriggerEvent.h"
37     #include "DataFormats/TrackReco/interface/Track.h"
38     #include "DataFormats/TrackReco/interface/TrackFwd.h"
39     #include "DataFormats/VertexReco/interface/Vertex.h"
40     #include "DataFormats/VertexReco/interface/VertexFwd.h"
41     #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
42     #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
43 pivarski 1.6 #include "DataFormats/MuonReco/interface/Muon.h"
44     #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
45     #include "DataFormats/MuonReco/interface/MuonSegmentMatch.h"
46 pivarski 1.1
47     #include "CommonTools/UtilAlgos/interface/TFileService.h"
48     #include "FWCore/ServiceRegistry/interface/Service.h"
49     #include "TTree.h"
50    
51     //
52     // class declaration
53     //
54    
55     class FitNtuple : public edm::EDAnalyzer {
56     public:
57     explicit FitNtuple(const edm::ParameterSet&);
58     ~FitNtuple();
59    
60    
61     private:
62     virtual void beginJob() ;
63     virtual void analyze(const edm::Event&, const edm::EventSetup&);
64     virtual void endJob() ;
65    
66     // ----------member data ---------------------------
67    
68     // input parameters
69     edm::InputTag m_src;
70     bool m_getQscale;
71     bool m_getAlternating;
72    
73     // all ntuples need these variables
74     Float_t m_qscale; // also known as "pt-hat" ($\hat{p}_T$)
75     Int_t m_run; // run and event number so that we can look at
76     Int_t m_event; // interesting cases in the event display
77     Int_t m_trigger; // satisfied trigger? 0 = failed all
78     // 1 = passed HLT_Mu5
79     // 2 = passed HLT_Mu9
80     // 4 = passed HLT_Mu11
81 pivarski 1.3 Float_t m_muon1pt;
82     Float_t m_muon1eta;
83     Float_t m_muon2pt;
84     Float_t m_muon2eta;
85     Float_t m_muon3pt;
86     Float_t m_muon3eta;
87     Float_t m_muon4pt;
88     Float_t m_muon4eta;
89     Float_t m_muontrigpt;
90     Float_t m_muontrigeta;
91 pivarski 1.1
92     // control sample "lowdimuon" (single mu-jet with two muons)
93     TTree *m_lowdimuon;
94     Float_t m_lowdimuon_genmass;
95     Float_t m_lowdimuon_mass;
96     Float_t m_lowdimuon_pt;
97     Float_t m_lowdimuon_eta;
98     Float_t m_lowdimuon_phi;
99     Float_t m_lowdimuon_dr;
100     Float_t m_lowdimuon_pluspx;
101     Float_t m_lowdimuon_pluspy;
102     Float_t m_lowdimuon_pluspz;
103     Float_t m_lowdimuon_minuspx;
104     Float_t m_lowdimuon_minuspy;
105     Float_t m_lowdimuon_minuspz;
106     Float_t m_lowdimuon_vprob;
107     Float_t m_lowdimuon_vx;
108     Float_t m_lowdimuon_vy;
109     Float_t m_lowdimuon_vz;
110     Float_t m_lowdimuon_iso;
111 pivarski 1.2 Float_t m_lowdimuon_compluspx;
112     Float_t m_lowdimuon_compluspy;
113     Float_t m_lowdimuon_compluspz;
114 pivarski 1.6 Float_t m_lowdimuon_comrotpluspx;
115     Float_t m_lowdimuon_comrotpluspy;
116     Float_t m_lowdimuon_comrotpluspz;
117     Int_t m_lowdimuon_plusmatches;
118     Float_t m_lowdimuon_plusst1x;
119     Float_t m_lowdimuon_plusst1xsig;
120     Float_t m_lowdimuon_plusst1dxdz;
121     Float_t m_lowdimuon_plusst1dxdzsig;
122     Int_t m_lowdimuon_plushits;
123     Float_t m_lowdimuon_plusnormchi2;
124     Int_t m_lowdimuon_minusmatches;
125     Float_t m_lowdimuon_minusst1x;
126     Float_t m_lowdimuon_minusst1xsig;
127     Float_t m_lowdimuon_minusst1dxdz;
128     Float_t m_lowdimuon_minusst1dxdzsig;
129     Int_t m_lowdimuon_minushits;
130     Float_t m_lowdimuon_minusnormchi2;
131    
132 pivarski 1.1 // signal region "highdimuon" (single mu-jet with two muons)
133     TTree *m_highdimuon;
134    
135     // control sample "muplustrack" (single mu-jet with three muons and the closest track)
136     TTree *m_muplustrack;
137    
138     // signal region "quadmuon" (single mu-jet with four muons)
139     TTree *m_quadmuon;
140    
141     // signal region "dimudimu" (two mu-jets with two muons each)
142     TTree *m_dimudimu;
143    
144     // signal region "dimuquadmu" (one two-muon jet and one four-muon jet)
145     TTree *m_dimuquadmu;
146    
147     // signal region "quadmuquadmu" (two mu-jets with four muons each)
148     TTree *m_quadmuquadmu;
149    
150     };
151    
152     //
153     // constants, enums and typedefs
154     //
155    
156     //
157     // static data member definitions
158     //
159    
160     //
161     // constructors and destructor
162     //
163     FitNtuple::FitNtuple(const edm::ParameterSet& iConfig)
164     : m_src(iConfig.getParameter<edm::InputTag>("src"))
165     , m_getQscale(iConfig.getParameter<bool>("getQscale"))
166     , m_getAlternating(iConfig.getParameter<bool>("getAlternating"))
167     {
168     //now do what ever initialization is needed
169     edm::Service<TFileService> tFile;
170    
171     m_lowdimuon = tFile->make<TTree>("lowdimuon", "lowdimuon");;
172     m_lowdimuon->Branch("qscale", &m_qscale, "qscale/F");
173     m_lowdimuon->Branch("run", &m_run, "run/I");
174     m_lowdimuon->Branch("event", &m_event, "event/I");
175     m_lowdimuon->Branch("trigger", &m_trigger, "trigger/I");
176 pivarski 1.3 m_lowdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
177     m_lowdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
178     m_lowdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
179     m_lowdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
180     m_lowdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
181     m_lowdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
182 pivarski 1.1 m_lowdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
183     m_lowdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
184     m_lowdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
185     m_lowdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
186     m_lowdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
187     m_lowdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
188     m_lowdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
189     m_lowdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
190     m_lowdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
191     m_lowdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
192     m_lowdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
193     m_lowdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
194     m_lowdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
195     m_lowdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
196     m_lowdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
197     m_lowdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
198     m_lowdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
199 pivarski 1.2 m_lowdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
200     m_lowdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
201     m_lowdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
202 pivarski 1.6 m_lowdimuon->Branch("comrotpluspx", &m_lowdimuon_comrotpluspx, "comrotpluspx/F");
203     m_lowdimuon->Branch("comrotpluspy", &m_lowdimuon_comrotpluspy, "comrotpluspy/F");
204     m_lowdimuon->Branch("comrotpluspz", &m_lowdimuon_comrotpluspz, "comrotpluspz/F");
205     m_lowdimuon->Branch("plusmatches", &m_lowdimuon_plusmatches, "plusmatches/I");
206     m_lowdimuon->Branch("plusst1x", &m_lowdimuon_plusst1x, "plusst1x/F");
207     m_lowdimuon->Branch("plusst1xsig", &m_lowdimuon_plusst1xsig, "plusst1xsig/F");
208     m_lowdimuon->Branch("plusst1dxdz", &m_lowdimuon_plusst1dxdz, "plusst1dxdz/F");
209     m_lowdimuon->Branch("plusst1dxdzsig", &m_lowdimuon_plusst1dxdzsig, "plusst1dxdzsig/F");
210     m_lowdimuon->Branch("plushits", &m_lowdimuon_plushits, "plushits/I");
211     m_lowdimuon->Branch("plusnormchi2", &m_lowdimuon_plusnormchi2, "plusnormchi2/F");
212     m_lowdimuon->Branch("minusmatches", &m_lowdimuon_minusmatches, "minusmatches/I");
213     m_lowdimuon->Branch("minusst1x", &m_lowdimuon_minusst1x, "minusst1x/F");
214     m_lowdimuon->Branch("minusst1xsig", &m_lowdimuon_minusst1xsig, "minusst1xsig/F");
215     m_lowdimuon->Branch("minusst1dxdz", &m_lowdimuon_minusst1dxdz, "minusst1dxdz/F");
216     m_lowdimuon->Branch("minusst1dxdzsig", &m_lowdimuon_minusst1dxdzsig, "minusst1dxdzsig/F");
217     m_lowdimuon->Branch("minushits", &m_lowdimuon_minushits, "minushits/I");
218     m_lowdimuon->Branch("minusnormchi2", &m_lowdimuon_minusnormchi2, "minusnormchi2/F");
219 pivarski 1.1
220     m_highdimuon = tFile->make<TTree>("highdimuon", "highdimuon");;
221     m_highdimuon->Branch("qscale", &m_qscale, "qscale/F");
222     m_highdimuon->Branch("run", &m_run, "run/I");
223     m_highdimuon->Branch("event", &m_event, "event/I");
224     m_highdimuon->Branch("trigger", &m_trigger, "trigger/I");
225 pivarski 1.3 m_highdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
226     m_highdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
227     m_highdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
228     m_highdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
229     m_highdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
230     m_highdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
231 pivarski 1.1
232     m_muplustrack = tFile->make<TTree>("muplustrack", "muplustrack");;
233     m_muplustrack->Branch("qscale", &m_qscale, "qscale/F");
234     m_muplustrack->Branch("run", &m_run, "run/I");
235     m_muplustrack->Branch("event", &m_event, "event/I");
236     m_muplustrack->Branch("trigger", &m_trigger, "trigger/I");
237 pivarski 1.3 m_muplustrack->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
238     m_muplustrack->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
239     m_muplustrack->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
240     m_muplustrack->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
241     m_muplustrack->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
242     m_muplustrack->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
243     m_muplustrack->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
244     m_muplustrack->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
245     m_muplustrack->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
246     m_muplustrack->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
247 pivarski 1.1
248     m_quadmuon = tFile->make<TTree>("quadmuon", "quadmuon");;
249     m_quadmuon->Branch("qscale", &m_qscale, "qscale/F");
250     m_quadmuon->Branch("run", &m_run, "run/I");
251     m_quadmuon->Branch("event", &m_event, "event/I");
252     m_quadmuon->Branch("trigger", &m_trigger, "trigger/I");
253 pivarski 1.3 m_quadmuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
254     m_quadmuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
255     m_quadmuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
256     m_quadmuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
257     m_quadmuon->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
258     m_quadmuon->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
259     m_quadmuon->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
260     m_quadmuon->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
261     m_quadmuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
262     m_quadmuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
263 pivarski 1.1
264     m_dimudimu = tFile->make<TTree>("dimudimu", "dimudimu");;
265     m_dimudimu->Branch("qscale", &m_qscale, "qscale/F");
266     m_dimudimu->Branch("run", &m_run, "run/I");
267     m_dimudimu->Branch("event", &m_event, "event/I");
268     m_dimudimu->Branch("trigger", &m_trigger, "trigger/I");
269 pivarski 1.3 m_dimudimu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
270     m_dimudimu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
271     m_dimudimu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
272     m_dimudimu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
273     m_dimudimu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
274     m_dimudimu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
275     m_dimudimu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
276     m_dimudimu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
277     m_dimudimu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
278     m_dimudimu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
279 pivarski 1.1
280     m_dimuquadmu = tFile->make<TTree>("dimuquadmu", "dimuquadmu");;
281     m_dimuquadmu->Branch("qscale", &m_qscale, "qscale/F");
282     m_dimuquadmu->Branch("run", &m_run, "run/I");
283     m_dimuquadmu->Branch("event", &m_event, "event/I");
284     m_dimuquadmu->Branch("trigger", &m_trigger, "trigger/I");
285 pivarski 1.3 m_dimuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
286     m_dimuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
287     m_dimuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
288     m_dimuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
289     m_dimuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
290     m_dimuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
291     m_dimuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
292     m_dimuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
293     m_dimuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
294     m_dimuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
295 pivarski 1.1
296     m_quadmuquadmu = tFile->make<TTree>("quadmuquadmu", "quadmuquadmu");;
297     m_quadmuquadmu->Branch("qscale", &m_qscale, "qscale/F");
298     m_quadmuquadmu->Branch("run", &m_run, "run/I");
299     m_quadmuquadmu->Branch("event", &m_event, "event/I");
300     m_quadmuquadmu->Branch("trigger", &m_trigger, "trigger/I");
301 pivarski 1.3 m_quadmuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
302     m_quadmuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
303     m_quadmuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
304     m_quadmuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
305     m_quadmuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
306     m_quadmuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
307     m_quadmuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
308     m_quadmuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
309     m_quadmuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
310     m_quadmuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
311 pivarski 1.1 }
312    
313    
314     FitNtuple::~FitNtuple()
315     {
316    
317     // do anything here that needs to be done at desctruction time
318     // (e.g. close files, deallocate resources etc.)
319    
320     }
321    
322    
323     //
324     // member functions
325     //
326    
327     // ------------ method called to for each event ------------
328     void
329     FitNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
330     {
331     // the QCD scale only exists for certain kinds of Monte Carlo, but if it's available, get it
332     m_qscale = 0.;
333     if (m_getQscale) {
334     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
335     iEvent.getByLabel("generator", genEventInfoProduct);
336     m_qscale = genEventInfoProduct->qScale();
337     }
338    
339     // the pair-gun Monte Carlos have a (never used) feature called
340     // alteration; the real events are the ones in which particleNumber == 0
341     bool alternating = true;
342     if (m_getAlternating) {
343     edm::Handle<unsigned int> particleNumber;
344     iEvent.getByLabel("generator", "particleNumber", particleNumber);
345     alternating = (*particleNumber == 0);
346     }
347     if (!alternating) return;
348    
349     // get the run and event number
350     m_run = iEvent.id().run();
351     m_event = iEvent.id().event();
352    
353     // mu-jets (muons grouped by mass and vertex compatibility)
354     edm::Handle<pat::MultiMuonCollection> muJets;
355     iEvent.getByLabel(m_src, muJets);
356    
357     // // orphans (muons not found in any group)
358     // edm::Handle<pat::MuonCollection> orphans;
359     // iEvent.getByLabel(m_src, "Orphans", orphans);
360    
361 pivarski 1.3 // find the top four muons in the event (-1000. if not found)
362     edm::Handle<pat::MuonCollection> muons;
363     iEvent.getByLabel("cleanPatMuons", muons);
364     m_muon1pt = -1000.; m_muon1eta = -1000.;
365     m_muon2pt = -1000.; m_muon2eta = -1000.;
366     m_muon3pt = -1000.; m_muon3eta = -1000.;
367     m_muon4pt = -1000.; m_muon4eta = -1000.;
368     m_muontrigpt = -1000.; m_muontrigeta = -1000.;
369     for (pat::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
370     if (fabs(muon->eta()) < 2.4 && muon->isTrackerMuon() && muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2 && muon->innerTrack()->numberOfValidHits() >= 8 && muon->innerTrack()->normalizedChi2() < 4.) {
371     if (muon->pt() > m_muon1pt) {
372     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
373     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
374     m_muon2pt = m_muon1pt; m_muon2eta = m_muon1eta;
375     m_muon1pt = muon->pt(); m_muon1eta = muon->eta();
376     }
377     else if (muon->pt() > m_muon2pt) {
378     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
379     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
380     m_muon2pt = muon->pt(); m_muon2eta = muon->eta();
381     }
382     else if (muon->pt() > m_muon3pt) {
383     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
384     m_muon3pt = muon->pt(); m_muon3eta = muon->eta();
385     }
386     else if (muon->pt() > m_muon4pt) {
387     m_muon4pt = muon->pt(); m_muon4eta = muon->eta();
388     }
389    
390     // special muon within a more limited |eta| range, to guarantee the trigger
391     if (fabs(muon->eta()) < 1.0) {
392     if (muon->pt() > m_muontrigpt) {
393     m_muontrigpt = muon->pt(); m_muontrigeta = muon->eta();
394     }
395     }
396     }
397     }
398 pivarski 1.1
399     // // all tracker-tracks
400     // edm::Handle<reco::TrackCollection> tracks;
401     // iEvent.getByLabel("generalTracks", tracks);
402    
403     // // generator-level 4-vectors
404     // edm::Handle<reco::GenParticleCollection> genParticles;
405     // iEvent.getByLabel("genParticles", genParticles);
406    
407     // find the closest primary vertex (in Z) to the first muJet with a valid vertex
408     edm::Handle<reco::VertexCollection> primaryVertices;
409     iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
410     reco::VertexCollection::const_iterator closestPrimaryVertex = primaryVertices->end();
411     if (muJets->size() > 0) {
412     pat::MultiMuonCollection::const_iterator muJet0 = muJets->end();
413     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
414     if (muJet->vertexValid()) {
415     muJet0 = muJet;
416     break;
417     }
418     }
419    
420     if (muJet0 != muJets->end()) {
421     for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin(); vertex != primaryVertices->end(); ++vertex) {
422     if (vertex->isValid() && !vertex->isFake() && vertex->tracksSize() > 3 && fabs(vertex->z()) < 24.) {
423     if (closestPrimaryVertex == primaryVertices->end() || fabs(vertex->z() - muJet0->vertexPoint().z()) < fabs(closestPrimaryVertex->z() - muJet0->vertexPoint().z())) {
424     closestPrimaryVertex = vertex;
425     }
426     } // end vertex quality cuts
427     } // end loop over primary vertices
428     } // end if muJet0 exists
429     } // end if muJets->size > 0
430    
431     // find out which trigger bits were fired
432     edm::Handle<pat::TriggerEvent> triggerEvent;
433     iEvent.getByLabel("patTriggerEvent", triggerEvent);
434     m_trigger = 0;
435     if (triggerEvent->path("HLT_Mu5") && triggerEvent->path("HLT_Mu5")->wasAccept()) m_trigger += 1;
436     if (triggerEvent->path("HLT_Mu9") && triggerEvent->path("HLT_Mu9")->wasAccept()) m_trigger += 2;
437     if (triggerEvent->path("HLT_Mu11") && triggerEvent->path("HLT_Mu11")->wasAccept()) m_trigger += 4;
438    
439     ////////////////////////////////////////////////////////// lowdimuon:
440     if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() < 80.) {
441     pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
442    
443     // generator-level mass using matched genParticles (for resolution of fit peak)
444     m_lowdimuon_genmass = -1000.;
445 pivarski 1.5 if (muJet->muon(0)->genParticlesSize() == 1 && muJet->muon(1)->genParticlesSize() == 1) {
446 pivarski 1.4 const reco::GenParticle *mu0 = muJet->muon(0)->genParticle();
447     const reco::GenParticle *mu1 = muJet->muon(1)->genParticle();
448 pivarski 1.1
449     double total_energy = mu0->energy() + mu1->energy();
450     double total_px = mu0->px() + mu1->px();
451     double total_py = mu0->py() + mu1->py();
452     double total_pz = mu0->pz() + mu1->pz();
453     m_lowdimuon_genmass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
454     }
455    
456     m_lowdimuon_mass = muJet->mass();
457     m_lowdimuon_pt = muJet->pt();
458     m_lowdimuon_eta = muJet->eta();
459     m_lowdimuon_phi = muJet->phi();
460     m_lowdimuon_dr = muJet->dRmax();
461 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
462     m_lowdimuon_pluspx = muJet->daughter(0)->px();
463     m_lowdimuon_pluspy = muJet->daughter(0)->py();
464     m_lowdimuon_pluspz = muJet->daughter(0)->pz();
465     m_lowdimuon_minuspx = muJet->daughter(1)->px();
466     m_lowdimuon_minuspy = muJet->daughter(1)->py();
467     m_lowdimuon_minuspz = muJet->daughter(1)->pz();
468     }
469     else {
470     m_lowdimuon_pluspx = muJet->daughter(1)->px();
471     m_lowdimuon_pluspy = muJet->daughter(1)->py();
472     m_lowdimuon_pluspz = muJet->daughter(1)->pz();
473     m_lowdimuon_minuspx = muJet->daughter(0)->px();
474     m_lowdimuon_minuspy = muJet->daughter(0)->py();
475     m_lowdimuon_minuspz = muJet->daughter(0)->pz();
476     }
477 pivarski 1.1 m_lowdimuon_vprob = -1000.;
478     m_lowdimuon_vx = -1000.;
479     m_lowdimuon_vy = -1000.;
480     m_lowdimuon_vz = -1000.;
481    
482 pivarski 1.2 GlobalVector complus;
483 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0);
484     else complus = muJet->daughterCOM(1);
485 pivarski 1.2 m_lowdimuon_compluspx = complus.x();
486     m_lowdimuon_compluspy = complus.y();
487     m_lowdimuon_compluspz = complus.z();
488    
489 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0);
490     else complus = muJet->daughterCOMrot(1);
491     m_lowdimuon_comrotpluspx = complus.x();
492     m_lowdimuon_comrotpluspy = complus.y();
493     m_lowdimuon_comrotpluspz = complus.z();
494    
495 pivarski 1.2 // replace all values with vertex-updated values if vertex fitting succeeded
496 pivarski 1.1 if (muJet->vertexValid()) {
497     m_lowdimuon_mass = muJet->vertexMass();
498     m_lowdimuon_pt = muJet->vertexMomentum().perp();
499     m_lowdimuon_eta = muJet->vertexMomentum().eta();
500     m_lowdimuon_phi = muJet->vertexMomentum().phi();
501     m_lowdimuon_dr = muJet->dRmax(true);
502 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
503     m_lowdimuon_pluspx = muJet->vertexMomentum(0).x();
504     m_lowdimuon_pluspy = muJet->vertexMomentum(0).y();
505     m_lowdimuon_pluspz = muJet->vertexMomentum(0).z();
506     m_lowdimuon_minuspx = muJet->vertexMomentum(1).x();
507     m_lowdimuon_minuspy = muJet->vertexMomentum(1).y();
508     m_lowdimuon_minuspz = muJet->vertexMomentum(1).z();
509     }
510     else {
511     m_lowdimuon_pluspx = muJet->vertexMomentum(1).x();
512     m_lowdimuon_pluspy = muJet->vertexMomentum(1).y();
513     m_lowdimuon_pluspz = muJet->vertexMomentum(1).z();
514     m_lowdimuon_minuspx = muJet->vertexMomentum(0).x();
515     m_lowdimuon_minuspy = muJet->vertexMomentum(0).y();
516     m_lowdimuon_minuspz = muJet->vertexMomentum(0).z();
517     }
518 pivarski 1.1 m_lowdimuon_vprob = muJet->vertexProb();
519    
520     if (closestPrimaryVertex != primaryVertices->end()) {
521     m_lowdimuon_vx = muJet->vertexPoint().x() - closestPrimaryVertex->x();
522     m_lowdimuon_vy = muJet->vertexPoint().y() - closestPrimaryVertex->y();
523     m_lowdimuon_vz = muJet->vertexPoint().z() - closestPrimaryVertex->z();
524     }
525    
526 pivarski 1.2 GlobalVector complus;
527 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0, true);
528     else complus = muJet->daughterCOM(1, true);
529 pivarski 1.2 m_lowdimuon_compluspx = complus.x();
530     m_lowdimuon_compluspy = complus.y();
531     m_lowdimuon_compluspz = complus.z();
532 pivarski 1.6
533     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0, true);
534     else complus = muJet->daughterCOMrot(1, true);
535     m_lowdimuon_comrotpluspx = complus.x();
536     m_lowdimuon_comrotpluspy = complus.y();
537     m_lowdimuon_comrotpluspz = complus.z();
538    
539 pivarski 1.1 } // end of replacements with quantities measured at the vertex
540    
541     m_lowdimuon_iso = muJet->centralTrackIsolation();
542    
543 pivarski 1.6 std::vector<reco::MuonChamberMatch> plusmatches, minusmatches;
544     if (muJet->daughter(0)->charge() > 0) {
545     m_lowdimuon_plusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
546     m_lowdimuon_minusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
547     plusmatches = muJet->muon(0)->matches();
548     minusmatches = muJet->muon(1)->matches();
549     m_lowdimuon_plushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
550     m_lowdimuon_minushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
551     m_lowdimuon_plusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
552     m_lowdimuon_minusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
553     }
554     else {
555     m_lowdimuon_plusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
556     m_lowdimuon_minusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
557     plusmatches = muJet->muon(1)->matches();
558     minusmatches = muJet->muon(0)->matches();
559     m_lowdimuon_plushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
560     m_lowdimuon_minushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
561     m_lowdimuon_plusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
562     m_lowdimuon_minusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
563     }
564    
565     m_lowdimuon_plusst1x = -1000.;
566     m_lowdimuon_plusst1xsig = -1000.;
567     m_lowdimuon_plusst1dxdz = -1000.;
568     m_lowdimuon_plusst1dxdzsig = -1000.;
569     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = plusmatches.begin(); chamberMatch != plusmatches.end(); ++chamberMatch) {
570     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
571     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
572     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
573     if (chamberMatch->station() == 1) {
574     m_lowdimuon_plusst1x = (segmentMatch->x - chamberMatch->x);
575     m_lowdimuon_plusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
576     m_lowdimuon_plusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
577     m_lowdimuon_plusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
578     }
579     }
580     }
581     }
582    
583     m_lowdimuon_minusst1x = -1000.;
584     m_lowdimuon_minusst1xsig = -1000.;
585     m_lowdimuon_minusst1dxdz = -1000.;
586     m_lowdimuon_minusst1dxdzsig = -1000.;
587     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = minusmatches.begin(); chamberMatch != minusmatches.end(); ++chamberMatch) {
588     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
589     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
590     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
591     if (chamberMatch->station() == 1) {
592     m_lowdimuon_minusst1x = (segmentMatch->x - chamberMatch->x);
593     m_lowdimuon_minusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
594     m_lowdimuon_minusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
595     m_lowdimuon_minusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
596     }
597     }
598     }
599     }
600    
601 pivarski 1.1 m_lowdimuon->Fill();
602     }
603    
604     ////////////////////////////////////////////////////////// highdimuon
605     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() > 80.) {
606     // ...
607     // m_highdimuon->Fill();
608     }
609    
610     ////////////////////////////////////////////////////////// muplustrack
611     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 3) {
612     // ...
613     // m_muplustrack->Fill();
614     }
615    
616     ////////////////////////////////////////////////////////// quadmuon
617     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 4) {
618     // ...
619     // m_quadmuon->Fill();
620     }
621    
622     ////////////////////////////////////////////////////////// dimudimu
623     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 2) {
624     // ...
625     // m_dimudimu->Fill();
626     }
627    
628     ////////////////////////////////////////////////////////// dimuquadmu
629     else if (muJets->size() == 2 && (
630     ((*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 4) ||
631     ((*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 2)
632     )) {
633     // ...
634     // m_dimuquadmu->Fill();
635     }
636    
637     ////////////////////////////////////////////////////////// quadmuquadmu
638     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 4) {
639     // ...
640     // m_quadmuquadmu->Fill();
641     }
642     }
643    
644    
645     // ------------ method called once each job just before starting event loop ------------
646     void
647     FitNtuple::beginJob()
648     {
649     }
650    
651     // ------------ method called once each job just after ending the event loop ------------
652     void
653     FitNtuple::endJob() {
654     }
655    
656     //define this as a plug-in
657     DEFINE_FWK_MODULE(FitNtuple);