ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/FitNtuple/src/FitNtuple.cc
Revision: 1.8
Committed: Wed Dec 1 00:32:28 2010 UTC (14 years, 5 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-11-30-a
Changes since 1.7: +228 -95 lines
Log Message:
dimuorphan to get a mass template for un-triggered muons

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.8 // $Id: FitNtuple.cc,v 1.7 2010/11/29 22:25:29 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 pivarski 1.7 #include "TRandom3.h"
51 pivarski 1.1
52     //
53     // class declaration
54     //
55    
56     class FitNtuple : public edm::EDAnalyzer {
57     public:
58     explicit FitNtuple(const edm::ParameterSet&);
59     ~FitNtuple();
60    
61    
62     private:
63     virtual void beginJob() ;
64     virtual void analyze(const edm::Event&, const edm::EventSetup&);
65     virtual void endJob() ;
66    
67     // ----------member data ---------------------------
68    
69     // input parameters
70     edm::InputTag m_src;
71 pivarski 1.8 edm::InputTag m_srcOrphans;
72 pivarski 1.1 bool m_getQscale;
73     bool m_getAlternating;
74    
75 pivarski 1.7 TRandom3 m_trandom3;
76    
77 pivarski 1.1 // all ntuples need these variables
78     Float_t m_qscale; // also known as "pt-hat" ($\hat{p}_T$)
79     Int_t m_run; // run and event number so that we can look at
80     Int_t m_event; // interesting cases in the event display
81     Int_t m_trigger; // satisfied trigger? 0 = failed all
82     // 1 = passed HLT_Mu5
83     // 2 = passed HLT_Mu9
84     // 4 = passed HLT_Mu11
85 pivarski 1.3 Float_t m_muon1pt;
86     Float_t m_muon1eta;
87     Float_t m_muon2pt;
88     Float_t m_muon2eta;
89     Float_t m_muon3pt;
90     Float_t m_muon3eta;
91     Float_t m_muon4pt;
92     Float_t m_muon4eta;
93     Float_t m_muontrigpt;
94     Float_t m_muontrigeta;
95 pivarski 1.1
96 pivarski 1.7 // for scaling to the Z peak
97     TTree *m_alldimuons;
98     Float_t m_alldimuons_mass;
99     Float_t m_alldimuons_iso;
100    
101 pivarski 1.1 // control sample "lowdimuon" (single mu-jet with two muons)
102     TTree *m_lowdimuon;
103     Float_t m_lowdimuon_genmass;
104     Float_t m_lowdimuon_mass;
105     Float_t m_lowdimuon_pt;
106     Float_t m_lowdimuon_eta;
107     Float_t m_lowdimuon_phi;
108     Float_t m_lowdimuon_dr;
109     Float_t m_lowdimuon_pluspx;
110     Float_t m_lowdimuon_pluspy;
111     Float_t m_lowdimuon_pluspz;
112     Float_t m_lowdimuon_minuspx;
113     Float_t m_lowdimuon_minuspy;
114     Float_t m_lowdimuon_minuspz;
115     Float_t m_lowdimuon_vprob;
116     Float_t m_lowdimuon_vx;
117     Float_t m_lowdimuon_vy;
118     Float_t m_lowdimuon_vz;
119     Float_t m_lowdimuon_iso;
120 pivarski 1.2 Float_t m_lowdimuon_compluspx;
121     Float_t m_lowdimuon_compluspy;
122     Float_t m_lowdimuon_compluspz;
123 pivarski 1.6 Float_t m_lowdimuon_comrotpluspx;
124     Float_t m_lowdimuon_comrotpluspy;
125     Float_t m_lowdimuon_comrotpluspz;
126     Int_t m_lowdimuon_plusmatches;
127     Float_t m_lowdimuon_plusst1x;
128     Float_t m_lowdimuon_plusst1xsig;
129     Float_t m_lowdimuon_plusst1dxdz;
130     Float_t m_lowdimuon_plusst1dxdzsig;
131     Int_t m_lowdimuon_plushits;
132     Float_t m_lowdimuon_plusnormchi2;
133     Int_t m_lowdimuon_minusmatches;
134     Float_t m_lowdimuon_minusst1x;
135     Float_t m_lowdimuon_minusst1xsig;
136     Float_t m_lowdimuon_minusst1dxdz;
137     Float_t m_lowdimuon_minusst1dxdzsig;
138     Int_t m_lowdimuon_minushits;
139     Float_t m_lowdimuon_minusnormchi2;
140    
141 pivarski 1.1 // signal region "highdimuon" (single mu-jet with two muons)
142     TTree *m_highdimuon;
143    
144     // control sample "muplustrack" (single mu-jet with three muons and the closest track)
145     TTree *m_muplustrack;
146    
147     // signal region "quadmuon" (single mu-jet with four muons)
148     TTree *m_quadmuon;
149    
150 pivarski 1.8 // control sample "dimuorphan" (single, two-muon mu-jet and a single muon)
151     TTree *m_dimuorphan;
152     Float_t m_dimuorphan_deltaphi;
153     Float_t m_dimuorphan_orphanpt;
154     Float_t m_dimuorphan_orphaneta;
155     Float_t m_dimuorphan_orphanphi;
156     Int_t m_dimuorphan_containstrig;
157     Float_t m_dimuorphan_mass;
158     Float_t m_dimuorphan_pt;
159     Float_t m_dimuorphan_eta;
160     Float_t m_dimuorphan_phi;
161     Float_t m_dimuorphan_dr;
162     Float_t m_dimuorphan_vprob;
163     Float_t m_dimuorphan_lxy;
164     Float_t m_dimuorphan_lxyz;
165     Float_t m_dimuorphan_iso;
166    
167 pivarski 1.1 // signal region "dimudimu" (two mu-jets with two muons each)
168     TTree *m_dimudimu;
169 pivarski 1.7 Float_t m_dimudimu_wholemass;
170     Float_t m_dimudimu_wholept;
171     Float_t m_dimudimu_deltaphi;
172 pivarski 1.8 Int_t m_dimudimu_containstrig;
173     Float_t m_dimudimu_massC;
174     Float_t m_dimudimu_ptC;
175     Float_t m_dimudimu_etaC;
176     Float_t m_dimudimu_phiC;
177     Float_t m_dimudimu_drC;
178     Float_t m_dimudimu_vprobC;
179     Float_t m_dimudimu_lxyC;
180     Float_t m_dimudimu_lxyzC;
181     Float_t m_dimudimu_isoC;
182     Float_t m_dimudimu_massF;
183     Float_t m_dimudimu_ptF;
184     Float_t m_dimudimu_etaF;
185     Float_t m_dimudimu_phiF;
186     Float_t m_dimudimu_drF;
187     Float_t m_dimudimu_vprobF;
188     Float_t m_dimudimu_lxyF;
189     Float_t m_dimudimu_lxyzF;
190     Float_t m_dimudimu_isoF;
191 pivarski 1.1
192     // signal region "dimuquadmu" (one two-muon jet and one four-muon jet)
193     TTree *m_dimuquadmu;
194    
195     // signal region "quadmuquadmu" (two mu-jets with four muons each)
196     TTree *m_quadmuquadmu;
197    
198     };
199    
200     //
201     // constants, enums and typedefs
202     //
203    
204     //
205     // static data member definitions
206     //
207    
208     //
209     // constructors and destructor
210     //
211     FitNtuple::FitNtuple(const edm::ParameterSet& iConfig)
212     : m_src(iConfig.getParameter<edm::InputTag>("src"))
213 pivarski 1.8 , m_srcOrphans(iConfig.getParameter<edm::InputTag>("srcOrphans"))
214 pivarski 1.1 , m_getQscale(iConfig.getParameter<bool>("getQscale"))
215     , m_getAlternating(iConfig.getParameter<bool>("getAlternating"))
216     {
217     //now do what ever initialization is needed
218     edm::Service<TFileService> tFile;
219    
220 pivarski 1.7 m_alldimuons = tFile->make<TTree>("alldimuons", "alldimuons");
221     m_alldimuons->Branch("trigger", &m_trigger, "trigger/I");
222     m_alldimuons->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
223     m_alldimuons->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
224     m_alldimuons->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
225     m_alldimuons->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
226     m_alldimuons->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
227     m_alldimuons->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
228     m_alldimuons->Branch("mass", &m_alldimuons_mass, "mass/F");
229     m_alldimuons->Branch("iso", &m_alldimuons_iso, "iso/F");
230    
231     m_lowdimuon = tFile->make<TTree>("lowdimuon", "lowdimuon");
232 pivarski 1.1 m_lowdimuon->Branch("qscale", &m_qscale, "qscale/F");
233     m_lowdimuon->Branch("run", &m_run, "run/I");
234     m_lowdimuon->Branch("event", &m_event, "event/I");
235     m_lowdimuon->Branch("trigger", &m_trigger, "trigger/I");
236 pivarski 1.3 m_lowdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
237     m_lowdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
238     m_lowdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
239     m_lowdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
240     m_lowdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
241     m_lowdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
242 pivarski 1.1 m_lowdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
243     m_lowdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
244     m_lowdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
245     m_lowdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
246     m_lowdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
247     m_lowdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
248     m_lowdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
249     m_lowdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
250     m_lowdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
251     m_lowdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
252     m_lowdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
253     m_lowdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
254     m_lowdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
255     m_lowdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
256     m_lowdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
257     m_lowdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
258     m_lowdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
259 pivarski 1.2 m_lowdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
260     m_lowdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
261     m_lowdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
262 pivarski 1.6 m_lowdimuon->Branch("comrotpluspx", &m_lowdimuon_comrotpluspx, "comrotpluspx/F");
263     m_lowdimuon->Branch("comrotpluspy", &m_lowdimuon_comrotpluspy, "comrotpluspy/F");
264     m_lowdimuon->Branch("comrotpluspz", &m_lowdimuon_comrotpluspz, "comrotpluspz/F");
265     m_lowdimuon->Branch("plusmatches", &m_lowdimuon_plusmatches, "plusmatches/I");
266     m_lowdimuon->Branch("plusst1x", &m_lowdimuon_plusst1x, "plusst1x/F");
267     m_lowdimuon->Branch("plusst1xsig", &m_lowdimuon_plusst1xsig, "plusst1xsig/F");
268     m_lowdimuon->Branch("plusst1dxdz", &m_lowdimuon_plusst1dxdz, "plusst1dxdz/F");
269     m_lowdimuon->Branch("plusst1dxdzsig", &m_lowdimuon_plusst1dxdzsig, "plusst1dxdzsig/F");
270     m_lowdimuon->Branch("plushits", &m_lowdimuon_plushits, "plushits/I");
271     m_lowdimuon->Branch("plusnormchi2", &m_lowdimuon_plusnormchi2, "plusnormchi2/F");
272     m_lowdimuon->Branch("minusmatches", &m_lowdimuon_minusmatches, "minusmatches/I");
273     m_lowdimuon->Branch("minusst1x", &m_lowdimuon_minusst1x, "minusst1x/F");
274     m_lowdimuon->Branch("minusst1xsig", &m_lowdimuon_minusst1xsig, "minusst1xsig/F");
275     m_lowdimuon->Branch("minusst1dxdz", &m_lowdimuon_minusst1dxdz, "minusst1dxdz/F");
276     m_lowdimuon->Branch("minusst1dxdzsig", &m_lowdimuon_minusst1dxdzsig, "minusst1dxdzsig/F");
277     m_lowdimuon->Branch("minushits", &m_lowdimuon_minushits, "minushits/I");
278     m_lowdimuon->Branch("minusnormchi2", &m_lowdimuon_minusnormchi2, "minusnormchi2/F");
279 pivarski 1.1
280     m_highdimuon = tFile->make<TTree>("highdimuon", "highdimuon");;
281     m_highdimuon->Branch("qscale", &m_qscale, "qscale/F");
282     m_highdimuon->Branch("run", &m_run, "run/I");
283     m_highdimuon->Branch("event", &m_event, "event/I");
284     m_highdimuon->Branch("trigger", &m_trigger, "trigger/I");
285 pivarski 1.3 m_highdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
286     m_highdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
287     m_highdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
288     m_highdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
289     m_highdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
290     m_highdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
291 pivarski 1.1
292     m_muplustrack = tFile->make<TTree>("muplustrack", "muplustrack");;
293     m_muplustrack->Branch("qscale", &m_qscale, "qscale/F");
294     m_muplustrack->Branch("run", &m_run, "run/I");
295     m_muplustrack->Branch("event", &m_event, "event/I");
296     m_muplustrack->Branch("trigger", &m_trigger, "trigger/I");
297 pivarski 1.3 m_muplustrack->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
298     m_muplustrack->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
299     m_muplustrack->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
300     m_muplustrack->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
301     m_muplustrack->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
302     m_muplustrack->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
303     m_muplustrack->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
304     m_muplustrack->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
305     m_muplustrack->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
306     m_muplustrack->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
307 pivarski 1.1
308     m_quadmuon = tFile->make<TTree>("quadmuon", "quadmuon");;
309     m_quadmuon->Branch("qscale", &m_qscale, "qscale/F");
310     m_quadmuon->Branch("run", &m_run, "run/I");
311     m_quadmuon->Branch("event", &m_event, "event/I");
312     m_quadmuon->Branch("trigger", &m_trigger, "trigger/I");
313 pivarski 1.3 m_quadmuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
314     m_quadmuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
315     m_quadmuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
316     m_quadmuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
317     m_quadmuon->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
318     m_quadmuon->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
319     m_quadmuon->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
320     m_quadmuon->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
321     m_quadmuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
322     m_quadmuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
323 pivarski 1.1
324 pivarski 1.8 m_dimuorphan = tFile->make<TTree>("dimuorphan", "dimuorphan");;
325     m_dimuorphan->Branch("qscale", &m_qscale, "qscale/F");
326     m_dimuorphan->Branch("run", &m_run, "run/I");
327     m_dimuorphan->Branch("event", &m_event, "event/I");
328     m_dimuorphan->Branch("trigger", &m_trigger, "trigger/I");
329     m_dimuorphan->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
330     m_dimuorphan->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
331     m_dimuorphan->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
332     m_dimuorphan->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
333     m_dimuorphan->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
334     m_dimuorphan->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
335     m_dimuorphan->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
336     m_dimuorphan->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
337     m_dimuorphan->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
338     m_dimuorphan->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
339     m_dimuorphan->Branch("deltaphi", &m_dimuorphan_deltaphi, "deltaphi/F");
340     m_dimuorphan->Branch("orphanpt", &m_dimuorphan_orphanpt, "orphanpt/F");
341     m_dimuorphan->Branch("orphaneta", &m_dimuorphan_orphaneta, "orphaneta/F");
342     m_dimuorphan->Branch("orphanphi", &m_dimuorphan_orphanphi, "orphanphi/F");
343     m_dimuorphan->Branch("containstrig", &m_dimuorphan_containstrig, "containstrig/I");
344     m_dimuorphan->Branch("mass", &m_dimuorphan_mass, "mass/F");
345     m_dimuorphan->Branch("pt", &m_dimuorphan_pt, "pt/F");
346     m_dimuorphan->Branch("eta", &m_dimuorphan_eta, "eta/F");
347     m_dimuorphan->Branch("phi", &m_dimuorphan_phi, "phi/F");
348     m_dimuorphan->Branch("dr", &m_dimuorphan_dr, "dr/F");
349     m_dimuorphan->Branch("vprob", &m_dimuorphan_vprob, "vprob/F");
350     m_dimuorphan->Branch("lxy", &m_dimuorphan_lxy, "lxy/F");
351     m_dimuorphan->Branch("lxyz", &m_dimuorphan_lxyz, "lxyz/F");
352     m_dimuorphan->Branch("iso", &m_dimuorphan_iso, "iso/F");
353    
354 pivarski 1.1 m_dimudimu = tFile->make<TTree>("dimudimu", "dimudimu");;
355     m_dimudimu->Branch("qscale", &m_qscale, "qscale/F");
356     m_dimudimu->Branch("run", &m_run, "run/I");
357     m_dimudimu->Branch("event", &m_event, "event/I");
358     m_dimudimu->Branch("trigger", &m_trigger, "trigger/I");
359 pivarski 1.3 m_dimudimu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
360     m_dimudimu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
361     m_dimudimu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
362     m_dimudimu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
363     m_dimudimu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
364     m_dimudimu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
365     m_dimudimu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
366     m_dimudimu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
367     m_dimudimu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
368     m_dimudimu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
369 pivarski 1.7 m_dimudimu->Branch("wholemass", &m_dimudimu_wholemass, "wholemass/F");
370     m_dimudimu->Branch("wholept", &m_dimudimu_wholept, "wholept/F");
371     m_dimudimu->Branch("deltaphi", &m_dimudimu_deltaphi, "deltaphi/F");
372 pivarski 1.8 m_dimudimu->Branch("containstrig", &m_dimudimu_containstrig, "containstrig/I");
373     m_dimudimu->Branch("massC", &m_dimudimu_massC, "massC/F");
374     m_dimudimu->Branch("ptC", &m_dimudimu_ptC, "ptC/F");
375     m_dimudimu->Branch("etaC", &m_dimudimu_etaC, "etaC/F");
376     m_dimudimu->Branch("phiC", &m_dimudimu_phiC, "phiC/F");
377     m_dimudimu->Branch("drC", &m_dimudimu_drC, "drC/F");
378     m_dimudimu->Branch("vprobC", &m_dimudimu_vprobC, "vprobC/F");
379     m_dimudimu->Branch("lxyC", &m_dimudimu_lxyC, "lxyC/F");
380     m_dimudimu->Branch("lxyzC", &m_dimudimu_lxyzC, "lxyzC/F");
381     m_dimudimu->Branch("isoC", &m_dimudimu_isoC, "isoC/F");
382     m_dimudimu->Branch("massF", &m_dimudimu_massF, "massF/F");
383     m_dimudimu->Branch("ptF", &m_dimudimu_ptF, "ptF/F");
384     m_dimudimu->Branch("etaF", &m_dimudimu_etaF, "etaF/F");
385     m_dimudimu->Branch("phiF", &m_dimudimu_phiF, "phiF/F");
386     m_dimudimu->Branch("drF", &m_dimudimu_drF, "drF/F");
387     m_dimudimu->Branch("vprobF", &m_dimudimu_vprobF, "vprobF/F");
388     m_dimudimu->Branch("lxyF", &m_dimudimu_lxyF, "lxyF/F");
389     m_dimudimu->Branch("lxyzF", &m_dimudimu_lxyzF, "lxyzF/F");
390     m_dimudimu->Branch("isoF", &m_dimudimu_isoF, "isoF/F");
391 pivarski 1.1
392     m_dimuquadmu = tFile->make<TTree>("dimuquadmu", "dimuquadmu");;
393     m_dimuquadmu->Branch("qscale", &m_qscale, "qscale/F");
394     m_dimuquadmu->Branch("run", &m_run, "run/I");
395     m_dimuquadmu->Branch("event", &m_event, "event/I");
396     m_dimuquadmu->Branch("trigger", &m_trigger, "trigger/I");
397 pivarski 1.3 m_dimuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
398     m_dimuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
399     m_dimuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
400     m_dimuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
401     m_dimuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
402     m_dimuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
403     m_dimuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
404     m_dimuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
405     m_dimuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
406     m_dimuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
407 pivarski 1.1
408     m_quadmuquadmu = tFile->make<TTree>("quadmuquadmu", "quadmuquadmu");;
409     m_quadmuquadmu->Branch("qscale", &m_qscale, "qscale/F");
410     m_quadmuquadmu->Branch("run", &m_run, "run/I");
411     m_quadmuquadmu->Branch("event", &m_event, "event/I");
412     m_quadmuquadmu->Branch("trigger", &m_trigger, "trigger/I");
413 pivarski 1.3 m_quadmuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
414     m_quadmuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
415     m_quadmuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
416     m_quadmuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
417     m_quadmuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
418     m_quadmuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
419     m_quadmuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
420     m_quadmuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
421     m_quadmuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
422     m_quadmuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
423 pivarski 1.1 }
424    
425    
426     FitNtuple::~FitNtuple()
427     {
428    
429     // do anything here that needs to be done at desctruction time
430     // (e.g. close files, deallocate resources etc.)
431    
432     }
433    
434    
435     //
436     // member functions
437     //
438    
439     // ------------ method called to for each event ------------
440     void
441     FitNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
442     {
443     // the QCD scale only exists for certain kinds of Monte Carlo, but if it's available, get it
444     m_qscale = 0.;
445     if (m_getQscale) {
446     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
447     iEvent.getByLabel("generator", genEventInfoProduct);
448     m_qscale = genEventInfoProduct->qScale();
449     }
450    
451     // the pair-gun Monte Carlos have a (never used) feature called
452     // alteration; the real events are the ones in which particleNumber == 0
453     bool alternating = true;
454     if (m_getAlternating) {
455     edm::Handle<unsigned int> particleNumber;
456     iEvent.getByLabel("generator", "particleNumber", particleNumber);
457     alternating = (*particleNumber == 0);
458     }
459     if (!alternating) return;
460    
461     // get the run and event number
462     m_run = iEvent.id().run();
463     m_event = iEvent.id().event();
464    
465     // mu-jets (muons grouped by mass and vertex compatibility)
466     edm::Handle<pat::MultiMuonCollection> muJets;
467     iEvent.getByLabel(m_src, muJets);
468    
469 pivarski 1.8 // orphans (muons not found in any group)
470     edm::Handle<pat::MuonCollection> orphans;
471     iEvent.getByLabel(m_srcOrphans, orphans);
472 pivarski 1.1
473 pivarski 1.3 // find the top four muons in the event (-1000. if not found)
474     edm::Handle<pat::MuonCollection> muons;
475     iEvent.getByLabel("cleanPatMuons", muons);
476     m_muon1pt = -1000.; m_muon1eta = -1000.;
477     m_muon2pt = -1000.; m_muon2eta = -1000.;
478     m_muon3pt = -1000.; m_muon3eta = -1000.;
479     m_muon4pt = -1000.; m_muon4eta = -1000.;
480     m_muontrigpt = -1000.; m_muontrigeta = -1000.;
481 pivarski 1.8 pat::MuonCollection::const_iterator muontrig = muons->end();
482 pivarski 1.3 for (pat::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
483     if (fabs(muon->eta()) < 2.4 && muon->isTrackerMuon() && muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2 && muon->innerTrack()->numberOfValidHits() >= 8 && muon->innerTrack()->normalizedChi2() < 4.) {
484     if (muon->pt() > m_muon1pt) {
485     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
486     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
487     m_muon2pt = m_muon1pt; m_muon2eta = m_muon1eta;
488     m_muon1pt = muon->pt(); m_muon1eta = muon->eta();
489     }
490     else if (muon->pt() > m_muon2pt) {
491     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
492     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
493     m_muon2pt = muon->pt(); m_muon2eta = muon->eta();
494     }
495     else if (muon->pt() > m_muon3pt) {
496     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
497     m_muon3pt = muon->pt(); m_muon3eta = muon->eta();
498     }
499     else if (muon->pt() > m_muon4pt) {
500     m_muon4pt = muon->pt(); m_muon4eta = muon->eta();
501     }
502    
503     // special muon within a more limited |eta| range, to guarantee the trigger
504 pivarski 1.8 if (fabs(muon->eta()) < 0.9) {
505 pivarski 1.3 if (muon->pt() > m_muontrigpt) {
506     m_muontrigpt = muon->pt(); m_muontrigeta = muon->eta();
507     }
508     }
509 pivarski 1.8
510     // most central muon with pT > 12 GeV/c
511     if (muon->pt() > 12.) {
512     if (muontrig == muons->end() || fabs(muon->eta()) < fabs(muontrig->eta())) {
513     muontrig = muon;
514     }
515     }
516 pivarski 1.3 }
517     }
518 pivarski 1.1
519     // // all tracker-tracks
520     // edm::Handle<reco::TrackCollection> tracks;
521     // iEvent.getByLabel("generalTracks", tracks);
522    
523     // // generator-level 4-vectors
524     // edm::Handle<reco::GenParticleCollection> genParticles;
525     // iEvent.getByLabel("genParticles", genParticles);
526    
527     // find the closest primary vertex (in Z) to the first muJet with a valid vertex
528     edm::Handle<reco::VertexCollection> primaryVertices;
529     iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
530     reco::VertexCollection::const_iterator closestPrimaryVertex = primaryVertices->end();
531     if (muJets->size() > 0) {
532     pat::MultiMuonCollection::const_iterator muJet0 = muJets->end();
533     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
534     if (muJet->vertexValid()) {
535     muJet0 = muJet;
536     break;
537     }
538     }
539    
540     if (muJet0 != muJets->end()) {
541     for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin(); vertex != primaryVertices->end(); ++vertex) {
542     if (vertex->isValid() && !vertex->isFake() && vertex->tracksSize() > 3 && fabs(vertex->z()) < 24.) {
543     if (closestPrimaryVertex == primaryVertices->end() || fabs(vertex->z() - muJet0->vertexPoint().z()) < fabs(closestPrimaryVertex->z() - muJet0->vertexPoint().z())) {
544     closestPrimaryVertex = vertex;
545     }
546     } // end vertex quality cuts
547     } // end loop over primary vertices
548     } // end if muJet0 exists
549     } // end if muJets->size > 0
550    
551     // find out which trigger bits were fired
552     edm::Handle<pat::TriggerEvent> triggerEvent;
553     iEvent.getByLabel("patTriggerEvent", triggerEvent);
554     m_trigger = 0;
555     if (triggerEvent->path("HLT_Mu5") && triggerEvent->path("HLT_Mu5")->wasAccept()) m_trigger += 1;
556     if (triggerEvent->path("HLT_Mu9") && triggerEvent->path("HLT_Mu9")->wasAccept()) m_trigger += 2;
557     if (triggerEvent->path("HLT_Mu11") && triggerEvent->path("HLT_Mu11")->wasAccept()) m_trigger += 4;
558    
559 pivarski 1.7 ////////////////////////////////////////////////////////// alldimuons (special case for scaling to the Z)
560     if (muons->size() == 2) {
561     double total_energy = (*muons)[0].energy() + (*muons)[1].energy();
562     double total_px = (*muons)[0].px() + (*muons)[1].px();
563     double total_py = (*muons)[0].py() + (*muons)[1].py();
564     double total_pz = (*muons)[0].pz() + (*muons)[1].pz();
565     m_alldimuons_mass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
566     if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2) {
567     m_alldimuons_iso = (*muJets)[0].centralTrackIsolation();
568     }
569     else {
570     m_alldimuons_iso = -1.;
571     }
572     m_alldimuons->Fill();
573     }
574    
575 pivarski 1.1 ////////////////////////////////////////////////////////// lowdimuon:
576 pivarski 1.8 if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() < 80. && orphans->size() == 0) {
577 pivarski 1.1 pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
578    
579     // generator-level mass using matched genParticles (for resolution of fit peak)
580     m_lowdimuon_genmass = -1000.;
581 pivarski 1.5 if (muJet->muon(0)->genParticlesSize() == 1 && muJet->muon(1)->genParticlesSize() == 1) {
582 pivarski 1.4 const reco::GenParticle *mu0 = muJet->muon(0)->genParticle();
583     const reco::GenParticle *mu1 = muJet->muon(1)->genParticle();
584 pivarski 1.1
585     double total_energy = mu0->energy() + mu1->energy();
586     double total_px = mu0->px() + mu1->px();
587     double total_py = mu0->py() + mu1->py();
588     double total_pz = mu0->pz() + mu1->pz();
589     m_lowdimuon_genmass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
590     }
591    
592     m_lowdimuon_mass = muJet->mass();
593     m_lowdimuon_pt = muJet->pt();
594     m_lowdimuon_eta = muJet->eta();
595     m_lowdimuon_phi = muJet->phi();
596     m_lowdimuon_dr = muJet->dRmax();
597 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
598     m_lowdimuon_pluspx = muJet->daughter(0)->px();
599     m_lowdimuon_pluspy = muJet->daughter(0)->py();
600     m_lowdimuon_pluspz = muJet->daughter(0)->pz();
601     m_lowdimuon_minuspx = muJet->daughter(1)->px();
602     m_lowdimuon_minuspy = muJet->daughter(1)->py();
603     m_lowdimuon_minuspz = muJet->daughter(1)->pz();
604     }
605     else {
606     m_lowdimuon_pluspx = muJet->daughter(1)->px();
607     m_lowdimuon_pluspy = muJet->daughter(1)->py();
608     m_lowdimuon_pluspz = muJet->daughter(1)->pz();
609     m_lowdimuon_minuspx = muJet->daughter(0)->px();
610     m_lowdimuon_minuspy = muJet->daughter(0)->py();
611     m_lowdimuon_minuspz = muJet->daughter(0)->pz();
612     }
613 pivarski 1.1 m_lowdimuon_vprob = -1000.;
614     m_lowdimuon_vx = -1000.;
615     m_lowdimuon_vy = -1000.;
616     m_lowdimuon_vz = -1000.;
617    
618 pivarski 1.2 GlobalVector complus;
619 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0);
620     else complus = muJet->daughterCOM(1);
621 pivarski 1.2 m_lowdimuon_compluspx = complus.x();
622     m_lowdimuon_compluspy = complus.y();
623     m_lowdimuon_compluspz = complus.z();
624    
625 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0);
626     else complus = muJet->daughterCOMrot(1);
627     m_lowdimuon_comrotpluspx = complus.x();
628     m_lowdimuon_comrotpluspy = complus.y();
629     m_lowdimuon_comrotpluspz = complus.z();
630    
631 pivarski 1.2 // replace all values with vertex-updated values if vertex fitting succeeded
632 pivarski 1.1 if (muJet->vertexValid()) {
633     m_lowdimuon_mass = muJet->vertexMass();
634     m_lowdimuon_pt = muJet->vertexMomentum().perp();
635     m_lowdimuon_eta = muJet->vertexMomentum().eta();
636     m_lowdimuon_phi = muJet->vertexMomentum().phi();
637     m_lowdimuon_dr = muJet->dRmax(true);
638 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
639     m_lowdimuon_pluspx = muJet->vertexMomentum(0).x();
640     m_lowdimuon_pluspy = muJet->vertexMomentum(0).y();
641     m_lowdimuon_pluspz = muJet->vertexMomentum(0).z();
642     m_lowdimuon_minuspx = muJet->vertexMomentum(1).x();
643     m_lowdimuon_minuspy = muJet->vertexMomentum(1).y();
644     m_lowdimuon_minuspz = muJet->vertexMomentum(1).z();
645     }
646     else {
647     m_lowdimuon_pluspx = muJet->vertexMomentum(1).x();
648     m_lowdimuon_pluspy = muJet->vertexMomentum(1).y();
649     m_lowdimuon_pluspz = muJet->vertexMomentum(1).z();
650     m_lowdimuon_minuspx = muJet->vertexMomentum(0).x();
651     m_lowdimuon_minuspy = muJet->vertexMomentum(0).y();
652     m_lowdimuon_minuspz = muJet->vertexMomentum(0).z();
653     }
654 pivarski 1.1 m_lowdimuon_vprob = muJet->vertexProb();
655    
656     if (closestPrimaryVertex != primaryVertices->end()) {
657     m_lowdimuon_vx = muJet->vertexPoint().x() - closestPrimaryVertex->x();
658     m_lowdimuon_vy = muJet->vertexPoint().y() - closestPrimaryVertex->y();
659     m_lowdimuon_vz = muJet->vertexPoint().z() - closestPrimaryVertex->z();
660     }
661    
662 pivarski 1.2 GlobalVector complus;
663 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0, true);
664     else complus = muJet->daughterCOM(1, true);
665 pivarski 1.2 m_lowdimuon_compluspx = complus.x();
666     m_lowdimuon_compluspy = complus.y();
667     m_lowdimuon_compluspz = complus.z();
668 pivarski 1.6
669     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0, true);
670     else complus = muJet->daughterCOMrot(1, true);
671     m_lowdimuon_comrotpluspx = complus.x();
672     m_lowdimuon_comrotpluspy = complus.y();
673     m_lowdimuon_comrotpluspz = complus.z();
674    
675 pivarski 1.1 } // end of replacements with quantities measured at the vertex
676    
677     m_lowdimuon_iso = muJet->centralTrackIsolation();
678    
679 pivarski 1.6 std::vector<reco::MuonChamberMatch> plusmatches, minusmatches;
680     if (muJet->daughter(0)->charge() > 0) {
681     m_lowdimuon_plusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
682     m_lowdimuon_minusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
683     plusmatches = muJet->muon(0)->matches();
684     minusmatches = muJet->muon(1)->matches();
685     m_lowdimuon_plushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
686     m_lowdimuon_minushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
687     m_lowdimuon_plusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
688     m_lowdimuon_minusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
689     }
690     else {
691     m_lowdimuon_plusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
692     m_lowdimuon_minusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
693     plusmatches = muJet->muon(1)->matches();
694     minusmatches = muJet->muon(0)->matches();
695     m_lowdimuon_plushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
696     m_lowdimuon_minushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
697     m_lowdimuon_plusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
698     m_lowdimuon_minusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
699     }
700    
701     m_lowdimuon_plusst1x = -1000.;
702     m_lowdimuon_plusst1xsig = -1000.;
703     m_lowdimuon_plusst1dxdz = -1000.;
704     m_lowdimuon_plusst1dxdzsig = -1000.;
705     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = plusmatches.begin(); chamberMatch != plusmatches.end(); ++chamberMatch) {
706     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
707     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
708     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
709     if (chamberMatch->station() == 1) {
710     m_lowdimuon_plusst1x = (segmentMatch->x - chamberMatch->x);
711     m_lowdimuon_plusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
712     m_lowdimuon_plusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
713     m_lowdimuon_plusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
714     }
715     }
716     }
717     }
718    
719     m_lowdimuon_minusst1x = -1000.;
720     m_lowdimuon_minusst1xsig = -1000.;
721     m_lowdimuon_minusst1dxdz = -1000.;
722     m_lowdimuon_minusst1dxdzsig = -1000.;
723     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = minusmatches.begin(); chamberMatch != minusmatches.end(); ++chamberMatch) {
724     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
725     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
726     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
727     if (chamberMatch->station() == 1) {
728     m_lowdimuon_minusst1x = (segmentMatch->x - chamberMatch->x);
729     m_lowdimuon_minusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
730     m_lowdimuon_minusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
731     m_lowdimuon_minusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
732     }
733     }
734     }
735     }
736    
737 pivarski 1.1 m_lowdimuon->Fill();
738     }
739    
740     ////////////////////////////////////////////////////////// highdimuon
741     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() > 80.) {
742     // ...
743     // m_highdimuon->Fill();
744     }
745    
746     ////////////////////////////////////////////////////////// muplustrack
747     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 3) {
748     // ...
749     // m_muplustrack->Fill();
750     }
751    
752     ////////////////////////////////////////////////////////// quadmuon
753     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 4) {
754     // ...
755     // m_quadmuon->Fill();
756     }
757    
758 pivarski 1.8 ////////////////////////////////////////////////////////// dimuorphan
759     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && orphans->size() == 1) {
760     pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
761     pat::MuonCollection::const_iterator orphan = orphans->begin();
762    
763     m_dimuorphan_deltaphi = muJet->phi() - orphan->phi();
764     while (m_dimuorphan_deltaphi > M_PI) m_dimuorphan_deltaphi -= 2.*M_PI;
765     while (m_dimuorphan_deltaphi < -M_PI) m_dimuorphan_deltaphi += 2.*M_PI;
766    
767     m_dimuorphan_orphanpt = orphan->pt();
768     m_dimuorphan_orphaneta = orphan->eta();
769     m_dimuorphan_orphanphi = orphan->phi();
770    
771     const pat::Muon *muon0 = muJet->muon(0);
772     const pat::Muon *muon1 = muJet->muon(1);
773     if (muontrig != muons->end() && muontrig->innerTrack().isAvailable() && ((muon0->innerTrack().isAvailable() && muJet->sameTrack(&*(muontrig->innerTrack()), &*(muon0->innerTrack()))) || (muon1->innerTrack().isAvailable() && muJet->sameTrack(&*(muontrig->innerTrack()), &*(muon1->innerTrack()))))) {
774     m_dimuorphan_containstrig = 1;
775     }
776     else {
777     m_dimuorphan_containstrig = 0;
778     }
779    
780     m_dimuorphan_mass = muJet->mass();
781     m_dimuorphan_pt = muJet->pt();
782     m_dimuorphan_eta = muJet->eta();
783     m_dimuorphan_phi = muJet->phi();
784     m_dimuorphan_dr = muJet->dRmax();
785     m_dimuorphan_vprob = -1.;
786     m_dimuorphan_lxy = -1000.;
787     m_dimuorphan_lxyz = -1000.;
788     m_dimuorphan_iso = muJet->centralTrackIsolation();
789    
790     if (muJet->vertexValid()) {
791     m_dimuorphan_deltaphi = muJet->vertexMomentum().phi() - orphan->phi();
792     while (m_dimuorphan_deltaphi > M_PI) m_dimuorphan_deltaphi -= 2.*M_PI;
793     while (m_dimuorphan_deltaphi < -M_PI) m_dimuorphan_deltaphi += 2.*M_PI;
794    
795     m_dimuorphan_mass = muJet->vertexMass();
796     m_dimuorphan_pt = muJet->vertexMomentum().perp();
797     m_dimuorphan_eta = muJet->vertexMomentum().eta();
798     m_dimuorphan_phi = muJet->vertexMomentum().phi();
799     m_dimuorphan_dr = muJet->dRmax(true);
800     m_dimuorphan_vprob = muJet->vertexProb();
801     m_dimuorphan_lxy = muJet->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
802     m_dimuorphan_lxyz = muJet->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
803     }
804    
805     m_dimuorphan->Fill();
806     }
807    
808 pivarski 1.1 ////////////////////////////////////////////////////////// dimudimu
809     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 2) {
810 pivarski 1.8 pat::MultiMuonCollection::const_iterator muJetC = muJets->end();
811     pat::MultiMuonCollection::const_iterator muJetF = muJets->end();
812     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
813     const pat::Muon *muon0 = muJet->muon(0);
814     const pat::Muon *muon1 = muJet->muon(1);
815     if (muontrig != muons->end() && muontrig->innerTrack().isAvailable() && ((muon0->innerTrack().isAvailable() && muJet->sameTrack(&*(muontrig->innerTrack()), &*(muon0->innerTrack()))) || (muon1->innerTrack().isAvailable() && muJet->sameTrack(&*(muontrig->innerTrack()), &*(muon1->innerTrack()))))) {
816     muJetC = muJet;
817     }
818     else {
819     muJetF = muJet;
820     }
821     }
822     if (muJetC != muJets->end() && muJetF != muJets->end()) {
823     m_dimudimu_containstrig = 1;
824     }
825     else {
826     m_dimudimu_containstrig = 0;
827     pat::MultiMuonCollection::const_iterator muJet0 = muJets->begin();
828     pat::MultiMuonCollection::const_iterator muJet1 = muJets->begin(); ++muJet1;
829     if (fabs(muJet0->eta()) < fabs(muJet1->eta())) {
830     muJetC = muJet0;
831     muJetF = muJet1;
832     }
833     else {
834     muJetC = muJet1;
835     muJetF = muJet0;
836     }
837     }
838    
839     double total_energy = muJetC->energy() + muJetF->energy();
840     double total_px = muJetC->px() + muJetF->px();
841     double total_py = muJetC->py() + muJetF->py();
842     double total_pz = muJetC->pz() + muJetF->pz();
843 pivarski 1.7
844     m_dimudimu_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
845     m_dimudimu_wholept = sqrt(total_px*total_px + total_py*total_py);
846    
847 pivarski 1.8 m_dimudimu_deltaphi = muJetC->phi() - muJetF->phi();
848 pivarski 1.7 while (m_dimudimu_deltaphi > M_PI) m_dimudimu_deltaphi -= 2.*M_PI;
849     while (m_dimudimu_deltaphi < -M_PI) m_dimudimu_deltaphi += 2.*M_PI;
850    
851 pivarski 1.8 m_dimudimu_massC = muJetC->mass();
852     m_dimudimu_ptC = muJetC->pt();
853     m_dimudimu_etaC = muJetC->eta();
854     m_dimudimu_phiC = muJetC->phi();
855     m_dimudimu_drC = muJetC->dRmax();
856     m_dimudimu_vprobC = -1.;
857     m_dimudimu_lxyC = -1000.;
858     m_dimudimu_lxyzC = -1000.;
859     m_dimudimu_isoC = muJetC->centralTrackIsolation();
860    
861     m_dimudimu_massF = muJetF->mass();
862     m_dimudimu_ptF = muJetF->pt();
863     m_dimudimu_etaF = muJetF->eta();
864     m_dimudimu_phiF = muJetF->phi();
865     m_dimudimu_drF = muJetF->dRmax();
866     m_dimudimu_vprobF = -1.;
867     m_dimudimu_lxyF = -1000.;
868     m_dimudimu_lxyzF = -1000.;
869     m_dimudimu_isoF = muJetF->centralTrackIsolation();
870    
871     if (muJetC->vertexValid() && muJetF->vertexValid()) {
872     total_energy = muJetC->vertexP4().energy() + muJetF->vertexP4().energy();
873     total_px = muJetC->vertexMomentum().x() + muJetF->vertexMomentum().x();
874     total_py = muJetC->vertexMomentum().y() + muJetF->vertexMomentum().y();
875     total_pz = muJetC->vertexMomentum().z() + muJetF->vertexMomentum().z();
876 pivarski 1.7
877     m_dimudimu_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
878     m_dimudimu_wholept = sqrt(total_px*total_px + total_py*total_py);
879    
880 pivarski 1.8 m_dimudimu_deltaphi = muJetC->vertexMomentum().phi() - muJetF->vertexMomentum().phi();
881 pivarski 1.7 while (m_dimudimu_deltaphi > M_PI) m_dimudimu_deltaphi -= 2.*M_PI;
882     while (m_dimudimu_deltaphi < -M_PI) m_dimudimu_deltaphi += 2.*M_PI;
883    
884 pivarski 1.8 m_dimudimu_massC = muJetC->vertexMass();
885     m_dimudimu_ptC = muJetC->vertexMomentum().perp();
886     m_dimudimu_etaC = muJetC->vertexMomentum().eta();
887     m_dimudimu_phiC = muJetC->vertexMomentum().phi();
888     m_dimudimu_drC = muJetC->dRmax(true);
889     m_dimudimu_vprobC = muJetC->vertexProb();
890     m_dimudimu_lxyC = muJetC->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
891     m_dimudimu_lxyzC = muJetC->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
892    
893     m_dimudimu_massF = muJetF->vertexMass();
894     m_dimudimu_ptF = muJetF->vertexMomentum().perp();
895     m_dimudimu_etaF = muJetF->vertexMomentum().eta();
896     m_dimudimu_phiF = muJetF->vertexMomentum().phi();
897     m_dimudimu_drF = muJetF->dRmax(true);
898     m_dimudimu_vprobF = muJetF->vertexProb();
899     m_dimudimu_lxyF = muJetF->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
900     m_dimudimu_lxyzF = muJetF->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
901 pivarski 1.7 }
902    
903     m_dimudimu->Fill();
904 pivarski 1.1 }
905    
906     ////////////////////////////////////////////////////////// dimuquadmu
907     else if (muJets->size() == 2 && (
908     ((*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 4) ||
909     ((*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 2)
910     )) {
911     // ...
912     // m_dimuquadmu->Fill();
913     }
914    
915     ////////////////////////////////////////////////////////// quadmuquadmu
916     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 4) {
917     // ...
918     // m_quadmuquadmu->Fill();
919     }
920     }
921    
922    
923     // ------------ method called once each job just before starting event loop ------------
924     void
925     FitNtuple::beginJob()
926     {
927     }
928    
929     // ------------ method called once each job just after ending the event loop ------------
930     void
931     FitNtuple::endJob() {
932     }
933    
934     //define this as a plug-in
935     DEFINE_FWK_MODULE(FitNtuple);