ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/FitNtuple/src/FitNtuple.cc
Revision: 1.12
Committed: Sat Dec 11 14:47:16 2010 UTC (14 years, 4 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-12-11-a
Changes since 1.11: +4 -4 lines
Log Message:
fix HLT_Mu15_v1 name

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.12 // $Id: FitNtuple.cc,v 1.11 2010/12/09 11:38:17 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     bool m_getQscale;
71     bool m_getAlternating;
72 pivarski 1.10 bool m_fillGenLevel;
73 pivarski 1.1
74 pivarski 1.7 TRandom3 m_trandom3;
75    
76 pivarski 1.1 // all ntuples need these variables
77     Float_t m_qscale; // also known as "pt-hat" ($\hat{p}_T$)
78 pivarski 1.10 Int_t m_genLevel; // gen-level information about all mu-jets
79     // 1 = b-jets
80     // 2 = resonances
81     // 4 = any decays-in-flight
82     // 8 = any non-muons
83 pivarski 1.1 Int_t m_run; // run and event number so that we can look at
84     Int_t m_event; // interesting cases in the event display
85     Int_t m_trigger; // satisfied trigger? 0 = failed all
86     // 1 = passed HLT_Mu5
87     // 2 = passed HLT_Mu9
88     // 4 = passed HLT_Mu11
89 pivarski 1.12 // 8 = passed HLT_Mu15_v1
90 pivarski 1.3 Float_t m_muon1pt;
91     Float_t m_muon1eta;
92     Float_t m_muon2pt;
93     Float_t m_muon2eta;
94     Float_t m_muon3pt;
95     Float_t m_muon3eta;
96     Float_t m_muon4pt;
97     Float_t m_muon4eta;
98     Float_t m_muontrigpt;
99     Float_t m_muontrigeta;
100 pivarski 1.1
101 pivarski 1.7 // for scaling to the Z peak
102     TTree *m_alldimuons;
103     Float_t m_alldimuons_mass;
104     Float_t m_alldimuons_iso;
105    
106 pivarski 1.1 // control sample "lowdimuon" (single mu-jet with two muons)
107     TTree *m_lowdimuon;
108 pivarski 1.10 Int_t m_lowdimuon_containstrig;
109 pivarski 1.1 Float_t m_lowdimuon_genmass;
110     Float_t m_lowdimuon_mass;
111     Float_t m_lowdimuon_pt;
112     Float_t m_lowdimuon_eta;
113     Float_t m_lowdimuon_phi;
114     Float_t m_lowdimuon_dr;
115     Float_t m_lowdimuon_pluspx;
116     Float_t m_lowdimuon_pluspy;
117     Float_t m_lowdimuon_pluspz;
118     Float_t m_lowdimuon_minuspx;
119     Float_t m_lowdimuon_minuspy;
120     Float_t m_lowdimuon_minuspz;
121     Float_t m_lowdimuon_vprob;
122     Float_t m_lowdimuon_vx;
123     Float_t m_lowdimuon_vy;
124     Float_t m_lowdimuon_vz;
125     Float_t m_lowdimuon_iso;
126 pivarski 1.2 Float_t m_lowdimuon_compluspx;
127     Float_t m_lowdimuon_compluspy;
128     Float_t m_lowdimuon_compluspz;
129 pivarski 1.6 Float_t m_lowdimuon_comrotpluspx;
130     Float_t m_lowdimuon_comrotpluspy;
131     Float_t m_lowdimuon_comrotpluspz;
132     Int_t m_lowdimuon_plusmatches;
133     Float_t m_lowdimuon_plusst1x;
134     Float_t m_lowdimuon_plusst1xsig;
135     Float_t m_lowdimuon_plusst1dxdz;
136     Float_t m_lowdimuon_plusst1dxdzsig;
137     Int_t m_lowdimuon_plushits;
138     Float_t m_lowdimuon_plusnormchi2;
139     Int_t m_lowdimuon_minusmatches;
140     Float_t m_lowdimuon_minusst1x;
141     Float_t m_lowdimuon_minusst1xsig;
142     Float_t m_lowdimuon_minusst1dxdz;
143     Float_t m_lowdimuon_minusst1dxdzsig;
144     Int_t m_lowdimuon_minushits;
145     Float_t m_lowdimuon_minusnormchi2;
146    
147 pivarski 1.1 // signal region "highdimuon" (single mu-jet with two muons)
148     TTree *m_highdimuon;
149    
150 pivarski 1.10 // control sample "mujetplustracks" (single mu-jet with 2-3 muons and 1-2 random, nearby tracks)
151     TTree *m_mujetplustracks;
152     Int_t m_mujetplustracks_containstrig;
153     Int_t m_mujetplustracks_muons;
154     Int_t m_mujetplustracks_fakes;
155 pivarski 1.11 Int_t m_mujetplustracks_charge;
156 pivarski 1.10 Float_t m_mujetplustracks_mass;
157 pivarski 1.11 Float_t m_mujetplustracks_massa;
158     Float_t m_mujetplustracks_massb;
159 pivarski 1.10 Float_t m_mujetplustracks_pt;
160     Float_t m_mujetplustracks_eta;
161     Float_t m_mujetplustracks_phi;
162     Float_t m_mujetplustracks_dr;
163     Float_t m_mujetplustracks_vprob;
164     Float_t m_mujetplustracks_lxy;
165     Float_t m_mujetplustracks_lxyz;
166     Float_t m_mujetplustracks_iso;
167     Float_t m_mujetplustracks_mupt1;
168     Float_t m_mujetplustracks_mupt2;
169     Float_t m_mujetplustracks_mupt3;
170     Float_t m_mujetplustracks_fakept1;
171     Float_t m_mujetplustracks_fakept2;
172 pivarski 1.1
173     // signal region "quadmuon" (single mu-jet with four muons)
174     TTree *m_quadmuon;
175    
176 pivarski 1.8 // control sample "dimuorphan" (single, two-muon mu-jet and a single muon)
177     TTree *m_dimuorphan;
178     Float_t m_dimuorphan_deltaphi;
179     Float_t m_dimuorphan_orphanpt;
180     Float_t m_dimuorphan_orphaneta;
181     Float_t m_dimuorphan_orphanphi;
182 pivarski 1.9 Int_t m_dimuorphan_orphanisglobal;
183     Int_t m_dimuorphan_orphanmatches;
184     Int_t m_dimuorphan_orphanhits;
185     Float_t m_dimuorphan_orphanchi2;
186     Float_t m_dimuorphan_orphanst1x;
187     Float_t m_dimuorphan_orphanst1xsig;
188     Float_t m_dimuorphan_orphanst1dxdz;
189     Float_t m_dimuorphan_orphanst1dxdzsig;
190 pivarski 1.8 Int_t m_dimuorphan_containstrig;
191     Float_t m_dimuorphan_mass;
192     Float_t m_dimuorphan_pt;
193     Float_t m_dimuorphan_eta;
194     Float_t m_dimuorphan_phi;
195     Float_t m_dimuorphan_dr;
196     Float_t m_dimuorphan_vprob;
197     Float_t m_dimuorphan_lxy;
198     Float_t m_dimuorphan_lxyz;
199     Float_t m_dimuorphan_iso;
200    
201 pivarski 1.1 // signal region "dimudimu" (two mu-jets with two muons each)
202     TTree *m_dimudimu;
203 pivarski 1.7 Float_t m_dimudimu_wholemass;
204     Float_t m_dimudimu_wholept;
205     Float_t m_dimudimu_deltaphi;
206 pivarski 1.8 Int_t m_dimudimu_containstrig;
207     Float_t m_dimudimu_massC;
208     Float_t m_dimudimu_ptC;
209     Float_t m_dimudimu_etaC;
210     Float_t m_dimudimu_phiC;
211     Float_t m_dimudimu_drC;
212     Float_t m_dimudimu_vprobC;
213     Float_t m_dimudimu_lxyC;
214     Float_t m_dimudimu_lxyzC;
215     Float_t m_dimudimu_isoC;
216     Float_t m_dimudimu_massF;
217     Float_t m_dimudimu_ptF;
218     Float_t m_dimudimu_etaF;
219     Float_t m_dimudimu_phiF;
220     Float_t m_dimudimu_drF;
221     Float_t m_dimudimu_vprobF;
222     Float_t m_dimudimu_lxyF;
223     Float_t m_dimudimu_lxyzF;
224     Float_t m_dimudimu_isoF;
225 pivarski 1.1
226     // signal region "dimuquadmu" (one two-muon jet and one four-muon jet)
227     TTree *m_dimuquadmu;
228    
229     // signal region "quadmuquadmu" (two mu-jets with four muons each)
230     TTree *m_quadmuquadmu;
231    
232     };
233    
234     //
235     // constants, enums and typedefs
236     //
237    
238     //
239     // static data member definitions
240     //
241    
242     //
243     // constructors and destructor
244     //
245     FitNtuple::FitNtuple(const edm::ParameterSet& iConfig)
246 pivarski 1.10 : m_getQscale(iConfig.getParameter<bool>("getQscale"))
247 pivarski 1.1 , m_getAlternating(iConfig.getParameter<bool>("getAlternating"))
248 pivarski 1.10 , m_fillGenLevel(iConfig.getParameter<bool>("fillGenLevel"))
249 pivarski 1.1 {
250     //now do what ever initialization is needed
251     edm::Service<TFileService> tFile;
252    
253 pivarski 1.7 m_alldimuons = tFile->make<TTree>("alldimuons", "alldimuons");
254     m_alldimuons->Branch("trigger", &m_trigger, "trigger/I");
255     m_alldimuons->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
256     m_alldimuons->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
257     m_alldimuons->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
258     m_alldimuons->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
259     m_alldimuons->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
260     m_alldimuons->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
261     m_alldimuons->Branch("mass", &m_alldimuons_mass, "mass/F");
262     m_alldimuons->Branch("iso", &m_alldimuons_iso, "iso/F");
263    
264     m_lowdimuon = tFile->make<TTree>("lowdimuon", "lowdimuon");
265 pivarski 1.1 m_lowdimuon->Branch("qscale", &m_qscale, "qscale/F");
266 pivarski 1.10 m_lowdimuon->Branch("genLevel", &m_genLevel, "genLevel/I");
267 pivarski 1.1 m_lowdimuon->Branch("run", &m_run, "run/I");
268     m_lowdimuon->Branch("event", &m_event, "event/I");
269     m_lowdimuon->Branch("trigger", &m_trigger, "trigger/I");
270 pivarski 1.3 m_lowdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
271     m_lowdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
272     m_lowdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
273     m_lowdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
274     m_lowdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
275     m_lowdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
276 pivarski 1.10 m_lowdimuon->Branch("containstrig", &m_lowdimuon_containstrig, "containstrig/I");
277 pivarski 1.1 m_lowdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
278     m_lowdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
279     m_lowdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
280     m_lowdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
281     m_lowdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
282     m_lowdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
283     m_lowdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
284     m_lowdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
285     m_lowdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
286     m_lowdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
287     m_lowdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
288     m_lowdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
289     m_lowdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
290     m_lowdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
291     m_lowdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
292     m_lowdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
293     m_lowdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
294 pivarski 1.2 m_lowdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
295     m_lowdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
296     m_lowdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
297 pivarski 1.6 m_lowdimuon->Branch("comrotpluspx", &m_lowdimuon_comrotpluspx, "comrotpluspx/F");
298     m_lowdimuon->Branch("comrotpluspy", &m_lowdimuon_comrotpluspy, "comrotpluspy/F");
299     m_lowdimuon->Branch("comrotpluspz", &m_lowdimuon_comrotpluspz, "comrotpluspz/F");
300     m_lowdimuon->Branch("plusmatches", &m_lowdimuon_plusmatches, "plusmatches/I");
301     m_lowdimuon->Branch("plusst1x", &m_lowdimuon_plusst1x, "plusst1x/F");
302     m_lowdimuon->Branch("plusst1xsig", &m_lowdimuon_plusst1xsig, "plusst1xsig/F");
303     m_lowdimuon->Branch("plusst1dxdz", &m_lowdimuon_plusst1dxdz, "plusst1dxdz/F");
304     m_lowdimuon->Branch("plusst1dxdzsig", &m_lowdimuon_plusst1dxdzsig, "plusst1dxdzsig/F");
305     m_lowdimuon->Branch("plushits", &m_lowdimuon_plushits, "plushits/I");
306     m_lowdimuon->Branch("plusnormchi2", &m_lowdimuon_plusnormchi2, "plusnormchi2/F");
307     m_lowdimuon->Branch("minusmatches", &m_lowdimuon_minusmatches, "minusmatches/I");
308     m_lowdimuon->Branch("minusst1x", &m_lowdimuon_minusst1x, "minusst1x/F");
309     m_lowdimuon->Branch("minusst1xsig", &m_lowdimuon_minusst1xsig, "minusst1xsig/F");
310     m_lowdimuon->Branch("minusst1dxdz", &m_lowdimuon_minusst1dxdz, "minusst1dxdz/F");
311     m_lowdimuon->Branch("minusst1dxdzsig", &m_lowdimuon_minusst1dxdzsig, "minusst1dxdzsig/F");
312     m_lowdimuon->Branch("minushits", &m_lowdimuon_minushits, "minushits/I");
313     m_lowdimuon->Branch("minusnormchi2", &m_lowdimuon_minusnormchi2, "minusnormchi2/F");
314 pivarski 1.1
315 pivarski 1.10 m_highdimuon = tFile->make<TTree>("highdimuon", "highdimuon");
316 pivarski 1.1 m_highdimuon->Branch("qscale", &m_qscale, "qscale/F");
317 pivarski 1.10 m_highdimuon->Branch("genLevel", &m_genLevel, "genLevel/I");
318 pivarski 1.1 m_highdimuon->Branch("run", &m_run, "run/I");
319     m_highdimuon->Branch("event", &m_event, "event/I");
320     m_highdimuon->Branch("trigger", &m_trigger, "trigger/I");
321 pivarski 1.3 m_highdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
322     m_highdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
323     m_highdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
324     m_highdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
325     m_highdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
326     m_highdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
327 pivarski 1.1
328 pivarski 1.10 m_mujetplustracks = tFile->make<TTree>("mujetplustracks", "mujetplustracks");
329     m_mujetplustracks->Branch("qscale", &m_qscale, "qscale/F");
330     m_mujetplustracks->Branch("run", &m_run, "run/I");
331     m_mujetplustracks->Branch("event", &m_event, "event/I");
332     m_mujetplustracks->Branch("trigger", &m_trigger, "trigger/I");
333     m_mujetplustracks->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
334     m_mujetplustracks->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
335     m_mujetplustracks->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
336     m_mujetplustracks->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
337     m_mujetplustracks->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
338     m_mujetplustracks->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
339     m_mujetplustracks->Branch("containstrig", &m_mujetplustracks_containstrig, "containstrig/I");
340     m_mujetplustracks->Branch("muons", &m_mujetplustracks_muons, "muons/I");
341     m_mujetplustracks->Branch("fakes", &m_mujetplustracks_fakes, "fakes/I");
342 pivarski 1.11 m_mujetplustracks->Branch("charge", &m_mujetplustracks_charge, "charge/I");
343 pivarski 1.10 m_mujetplustracks->Branch("mass", &m_mujetplustracks_mass, "mass/F");
344 pivarski 1.11 m_mujetplustracks->Branch("massa", &m_mujetplustracks_massa, "massa/F");
345     m_mujetplustracks->Branch("massb", &m_mujetplustracks_massb, "massb/F");
346 pivarski 1.10 m_mujetplustracks->Branch("pt", &m_mujetplustracks_pt, "pt/F");
347     m_mujetplustracks->Branch("eta", &m_mujetplustracks_eta, "eta/F");
348     m_mujetplustracks->Branch("phi", &m_mujetplustracks_phi, "phi/F");
349     m_mujetplustracks->Branch("dr", &m_mujetplustracks_dr, "dr/F");
350     m_mujetplustracks->Branch("vprob", &m_mujetplustracks_vprob, "vprob/F");
351     m_mujetplustracks->Branch("lxy", &m_mujetplustracks_lxy, "lxy/F");
352     m_mujetplustracks->Branch("lxyz", &m_mujetplustracks_lxyz, "lxyz/F");
353     m_mujetplustracks->Branch("iso", &m_mujetplustracks_iso, "iso/F");
354     m_mujetplustracks->Branch("mupt1", &m_mujetplustracks_mupt1, "mupt1/F");
355     m_mujetplustracks->Branch("mupt2", &m_mujetplustracks_mupt2, "mupt2/F");
356     m_mujetplustracks->Branch("mupt3", &m_mujetplustracks_mupt3, "mupt3/F");
357     m_mujetplustracks->Branch("fakept1", &m_mujetplustracks_fakept1, "fakept1/F");
358     m_mujetplustracks->Branch("fakept2", &m_mujetplustracks_fakept2, "fakept2/F");
359 pivarski 1.1
360 pivarski 1.10 m_quadmuon = tFile->make<TTree>("quadmuon", "quadmuon");
361 pivarski 1.1 m_quadmuon->Branch("qscale", &m_qscale, "qscale/F");
362 pivarski 1.10 m_quadmuon->Branch("genLevel", &m_genLevel, "genLevel/I");
363 pivarski 1.1 m_quadmuon->Branch("run", &m_run, "run/I");
364     m_quadmuon->Branch("event", &m_event, "event/I");
365     m_quadmuon->Branch("trigger", &m_trigger, "trigger/I");
366 pivarski 1.3 m_quadmuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
367     m_quadmuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
368     m_quadmuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
369     m_quadmuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
370     m_quadmuon->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
371     m_quadmuon->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
372     m_quadmuon->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
373     m_quadmuon->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
374     m_quadmuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
375     m_quadmuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
376 pivarski 1.1
377 pivarski 1.10 m_dimuorphan = tFile->make<TTree>("dimuorphan", "dimuorphan");
378 pivarski 1.8 m_dimuorphan->Branch("qscale", &m_qscale, "qscale/F");
379 pivarski 1.10 m_dimuorphan->Branch("genLevel", &m_genLevel, "genLevel/I");
380 pivarski 1.8 m_dimuorphan->Branch("run", &m_run, "run/I");
381     m_dimuorphan->Branch("event", &m_event, "event/I");
382     m_dimuorphan->Branch("trigger", &m_trigger, "trigger/I");
383     m_dimuorphan->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
384     m_dimuorphan->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
385     m_dimuorphan->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
386     m_dimuorphan->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
387     m_dimuorphan->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
388     m_dimuorphan->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
389     m_dimuorphan->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
390     m_dimuorphan->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
391     m_dimuorphan->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
392     m_dimuorphan->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
393     m_dimuorphan->Branch("deltaphi", &m_dimuorphan_deltaphi, "deltaphi/F");
394     m_dimuorphan->Branch("orphanpt", &m_dimuorphan_orphanpt, "orphanpt/F");
395     m_dimuorphan->Branch("orphaneta", &m_dimuorphan_orphaneta, "orphaneta/F");
396     m_dimuorphan->Branch("orphanphi", &m_dimuorphan_orphanphi, "orphanphi/F");
397 pivarski 1.9 m_dimuorphan->Branch("orphanisglobal", &m_dimuorphan_orphanisglobal, "orphanisglobal/I");
398     m_dimuorphan->Branch("orphanmatches", &m_dimuorphan_orphanmatches, "orphan/I");
399     m_dimuorphan->Branch("orphanhits", &m_dimuorphan_orphanhits, "orphanhits/I");
400     m_dimuorphan->Branch("orphanchi2", &m_dimuorphan_orphanchi2, "orphanchi2/F");
401     m_dimuorphan->Branch("orphanst1x", &m_dimuorphan_orphanst1x, "orphanst1x/F");
402     m_dimuorphan->Branch("orphanst1xsig", &m_dimuorphan_orphanst1xsig, "orphanst1xsig/F");
403     m_dimuorphan->Branch("orphanst1dxdz", &m_dimuorphan_orphanst1dxdz, "orphanst1dxdz/F");
404     m_dimuorphan->Branch("orphanst1dxdzsig", &m_dimuorphan_orphanst1dxdzsig, "orphanst1dxdzsig/F");
405 pivarski 1.8 m_dimuorphan->Branch("containstrig", &m_dimuorphan_containstrig, "containstrig/I");
406     m_dimuorphan->Branch("mass", &m_dimuorphan_mass, "mass/F");
407     m_dimuorphan->Branch("pt", &m_dimuorphan_pt, "pt/F");
408     m_dimuorphan->Branch("eta", &m_dimuorphan_eta, "eta/F");
409     m_dimuorphan->Branch("phi", &m_dimuorphan_phi, "phi/F");
410     m_dimuorphan->Branch("dr", &m_dimuorphan_dr, "dr/F");
411     m_dimuorphan->Branch("vprob", &m_dimuorphan_vprob, "vprob/F");
412     m_dimuorphan->Branch("lxy", &m_dimuorphan_lxy, "lxy/F");
413     m_dimuorphan->Branch("lxyz", &m_dimuorphan_lxyz, "lxyz/F");
414     m_dimuorphan->Branch("iso", &m_dimuorphan_iso, "iso/F");
415    
416 pivarski 1.10 m_dimudimu = tFile->make<TTree>("dimudimu", "dimudimu");
417 pivarski 1.1 m_dimudimu->Branch("qscale", &m_qscale, "qscale/F");
418 pivarski 1.10 m_dimudimu->Branch("genLevel", &m_genLevel, "genLevel/I");
419 pivarski 1.1 m_dimudimu->Branch("run", &m_run, "run/I");
420     m_dimudimu->Branch("event", &m_event, "event/I");
421     m_dimudimu->Branch("trigger", &m_trigger, "trigger/I");
422 pivarski 1.3 m_dimudimu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
423     m_dimudimu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
424     m_dimudimu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
425     m_dimudimu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
426     m_dimudimu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
427     m_dimudimu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
428     m_dimudimu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
429     m_dimudimu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
430     m_dimudimu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
431     m_dimudimu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
432 pivarski 1.7 m_dimudimu->Branch("wholemass", &m_dimudimu_wholemass, "wholemass/F");
433     m_dimudimu->Branch("wholept", &m_dimudimu_wholept, "wholept/F");
434     m_dimudimu->Branch("deltaphi", &m_dimudimu_deltaphi, "deltaphi/F");
435 pivarski 1.8 m_dimudimu->Branch("containstrig", &m_dimudimu_containstrig, "containstrig/I");
436     m_dimudimu->Branch("massC", &m_dimudimu_massC, "massC/F");
437     m_dimudimu->Branch("ptC", &m_dimudimu_ptC, "ptC/F");
438     m_dimudimu->Branch("etaC", &m_dimudimu_etaC, "etaC/F");
439     m_dimudimu->Branch("phiC", &m_dimudimu_phiC, "phiC/F");
440     m_dimudimu->Branch("drC", &m_dimudimu_drC, "drC/F");
441     m_dimudimu->Branch("vprobC", &m_dimudimu_vprobC, "vprobC/F");
442     m_dimudimu->Branch("lxyC", &m_dimudimu_lxyC, "lxyC/F");
443     m_dimudimu->Branch("lxyzC", &m_dimudimu_lxyzC, "lxyzC/F");
444     m_dimudimu->Branch("isoC", &m_dimudimu_isoC, "isoC/F");
445     m_dimudimu->Branch("massF", &m_dimudimu_massF, "massF/F");
446     m_dimudimu->Branch("ptF", &m_dimudimu_ptF, "ptF/F");
447     m_dimudimu->Branch("etaF", &m_dimudimu_etaF, "etaF/F");
448     m_dimudimu->Branch("phiF", &m_dimudimu_phiF, "phiF/F");
449     m_dimudimu->Branch("drF", &m_dimudimu_drF, "drF/F");
450     m_dimudimu->Branch("vprobF", &m_dimudimu_vprobF, "vprobF/F");
451     m_dimudimu->Branch("lxyF", &m_dimudimu_lxyF, "lxyF/F");
452     m_dimudimu->Branch("lxyzF", &m_dimudimu_lxyzF, "lxyzF/F");
453     m_dimudimu->Branch("isoF", &m_dimudimu_isoF, "isoF/F");
454 pivarski 1.1
455 pivarski 1.10 m_dimuquadmu = tFile->make<TTree>("dimuquadmu", "dimuquadmu");
456 pivarski 1.1 m_dimuquadmu->Branch("qscale", &m_qscale, "qscale/F");
457 pivarski 1.10 m_dimuquadmu->Branch("genLevel", &m_genLevel, "genLevel/I");
458 pivarski 1.1 m_dimuquadmu->Branch("run", &m_run, "run/I");
459     m_dimuquadmu->Branch("event", &m_event, "event/I");
460     m_dimuquadmu->Branch("trigger", &m_trigger, "trigger/I");
461 pivarski 1.3 m_dimuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
462     m_dimuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
463     m_dimuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
464     m_dimuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
465     m_dimuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
466     m_dimuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
467     m_dimuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
468     m_dimuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
469     m_dimuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
470     m_dimuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
471 pivarski 1.1
472 pivarski 1.10 m_quadmuquadmu = tFile->make<TTree>("quadmuquadmu", "quadmuquadmu");
473 pivarski 1.1 m_quadmuquadmu->Branch("qscale", &m_qscale, "qscale/F");
474 pivarski 1.10 m_quadmuquadmu->Branch("genLevel", &m_genLevel, "genLevel/I");
475 pivarski 1.1 m_quadmuquadmu->Branch("run", &m_run, "run/I");
476     m_quadmuquadmu->Branch("event", &m_event, "event/I");
477     m_quadmuquadmu->Branch("trigger", &m_trigger, "trigger/I");
478 pivarski 1.3 m_quadmuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
479     m_quadmuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
480     m_quadmuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
481     m_quadmuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
482     m_quadmuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
483     m_quadmuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
484     m_quadmuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
485     m_quadmuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
486     m_quadmuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
487     m_quadmuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
488 pivarski 1.1 }
489    
490    
491     FitNtuple::~FitNtuple()
492     {
493    
494     // do anything here that needs to be done at desctruction time
495     // (e.g. close files, deallocate resources etc.)
496    
497     }
498    
499    
500     //
501     // member functions
502     //
503    
504     // ------------ method called to for each event ------------
505     void
506     FitNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
507     {
508     // the QCD scale only exists for certain kinds of Monte Carlo, but if it's available, get it
509     m_qscale = 0.;
510     if (m_getQscale) {
511     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
512     iEvent.getByLabel("generator", genEventInfoProduct);
513     m_qscale = genEventInfoProduct->qScale();
514     }
515    
516     // the pair-gun Monte Carlos have a (never used) feature called
517     // alteration; the real events are the ones in which particleNumber == 0
518     bool alternating = true;
519     if (m_getAlternating) {
520     edm::Handle<unsigned int> particleNumber;
521     iEvent.getByLabel("generator", "particleNumber", particleNumber);
522     alternating = (*particleNumber == 0);
523     }
524     if (!alternating) return;
525    
526     // get the run and event number
527     m_run = iEvent.id().run();
528     m_event = iEvent.id().event();
529    
530     // mu-jets (muons grouped by mass and vertex compatibility)
531     edm::Handle<pat::MultiMuonCollection> muJets;
532 pivarski 1.10 iEvent.getByLabel("MuJetProducer", muJets);
533 pivarski 1.1
534 pivarski 1.8 // orphans (muons not found in any group)
535     edm::Handle<pat::MuonCollection> orphans;
536 pivarski 1.10 iEvent.getByLabel("MuJetProducer", "Orphans", orphans);
537    
538     // mu-jets plus at most two tracks as fake muons
539     edm::Handle<pat::MultiMuonCollection> muJetPlusTracks;
540     iEvent.getByLabel("MuJetPlusTracks", muJetPlusTracks);
541 pivarski 1.1
542 pivarski 1.3 // find the top four muons in the event (-1000. if not found)
543     edm::Handle<pat::MuonCollection> muons;
544 pivarski 1.10 iEvent.getByLabel("cleanPatMuonsTriggerMatch", muons);
545 pivarski 1.3 m_muon1pt = -1000.; m_muon1eta = -1000.;
546     m_muon2pt = -1000.; m_muon2eta = -1000.;
547     m_muon3pt = -1000.; m_muon3eta = -1000.;
548     m_muon4pt = -1000.; m_muon4eta = -1000.;
549     m_muontrigpt = -1000.; m_muontrigeta = -1000.;
550 pivarski 1.10 std::vector<pat::MuonCollection::const_iterator> muontrig;
551    
552 pivarski 1.3 for (pat::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
553     if (fabs(muon->eta()) < 2.4 && muon->isTrackerMuon() && muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2 && muon->innerTrack()->numberOfValidHits() >= 8 && muon->innerTrack()->normalizedChi2() < 4.) {
554     if (muon->pt() > m_muon1pt) {
555     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
556     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
557     m_muon2pt = m_muon1pt; m_muon2eta = m_muon1eta;
558     m_muon1pt = muon->pt(); m_muon1eta = muon->eta();
559     }
560     else if (muon->pt() > m_muon2pt) {
561     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
562     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
563     m_muon2pt = muon->pt(); m_muon2eta = muon->eta();
564     }
565     else if (muon->pt() > m_muon3pt) {
566     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
567     m_muon3pt = muon->pt(); m_muon3eta = muon->eta();
568     }
569     else if (muon->pt() > m_muon4pt) {
570     m_muon4pt = muon->pt(); m_muon4eta = muon->eta();
571     }
572    
573     // special muon within a more limited |eta| range, to guarantee the trigger
574 pivarski 1.8 if (fabs(muon->eta()) < 0.9) {
575 pivarski 1.3 if (muon->pt() > m_muontrigpt) {
576     m_muontrigpt = muon->pt(); m_muontrigeta = muon->eta();
577     }
578     }
579 pivarski 1.8
580 pivarski 1.10 // all muons with pT > 15 GeV/c and |eta| < 0.9, matched to a trigger object
581     if (muon->pt() > 15. && fabs(muon->eta()) < 0.9) {
582     const pat::TriggerObjectStandAlone *mu3 = muon->triggerObjectMatchByPath("HLT_Mu3");
583     const pat::TriggerObjectStandAlone *mu5 = muon->triggerObjectMatchByPath("HLT_Mu5");
584     const pat::TriggerObjectStandAlone *mu9 = muon->triggerObjectMatchByPath("HLT_Mu9");
585     const pat::TriggerObjectStandAlone *mu11 = muon->triggerObjectMatchByPath("HLT_Mu11");
586 pivarski 1.12 const pat::TriggerObjectStandAlone *mu15 = muon->triggerObjectMatchByPath("HLT_Mu15_v1");
587 pivarski 1.10
588     if ((mu3 != NULL && mu3->collection() == std::string("hltL3MuonCandidates::HLT") && mu3->pt() > 15.) ||
589     (mu5 != NULL && mu5->collection() == std::string("hltL3MuonCandidates::HLT") && mu5->pt() > 15.) ||
590     (mu9 != NULL && mu9->collection() == std::string("hltL3MuonCandidates::HLT") && mu9->pt() > 15.) ||
591     (mu11 != NULL && mu11->collection() == std::string("hltL3MuonCandidates::HLT") && mu11->pt() > 15.) ||
592     (mu15 != NULL && mu15->collection() == std::string("hltL3MuonCandidates::HLT") && mu15->pt() > 15.)) {
593     muontrig.push_back(muon);
594 pivarski 1.8 }
595     }
596 pivarski 1.3 }
597     }
598 pivarski 1.1
599     // // all tracker-tracks
600     // edm::Handle<reco::TrackCollection> tracks;
601     // iEvent.getByLabel("generalTracks", tracks);
602    
603     // // generator-level 4-vectors
604     // edm::Handle<reco::GenParticleCollection> genParticles;
605     // iEvent.getByLabel("genParticles", genParticles);
606    
607     // find the closest primary vertex (in Z) to the first muJet with a valid vertex
608     edm::Handle<reco::VertexCollection> primaryVertices;
609     iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
610     reco::VertexCollection::const_iterator closestPrimaryVertex = primaryVertices->end();
611     if (muJets->size() > 0) {
612     pat::MultiMuonCollection::const_iterator muJet0 = muJets->end();
613     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
614     if (muJet->vertexValid()) {
615     muJet0 = muJet;
616     break;
617     }
618     }
619    
620     if (muJet0 != muJets->end()) {
621     for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin(); vertex != primaryVertices->end(); ++vertex) {
622     if (vertex->isValid() && !vertex->isFake() && vertex->tracksSize() > 3 && fabs(vertex->z()) < 24.) {
623     if (closestPrimaryVertex == primaryVertices->end() || fabs(vertex->z() - muJet0->vertexPoint().z()) < fabs(closestPrimaryVertex->z() - muJet0->vertexPoint().z())) {
624     closestPrimaryVertex = vertex;
625     }
626     } // end vertex quality cuts
627     } // end loop over primary vertices
628     } // end if muJet0 exists
629     } // end if muJets->size > 0
630    
631     // find out which trigger bits were fired
632     edm::Handle<pat::TriggerEvent> triggerEvent;
633     iEvent.getByLabel("patTriggerEvent", triggerEvent);
634     m_trigger = 0;
635     if (triggerEvent->path("HLT_Mu5") && triggerEvent->path("HLT_Mu5")->wasAccept()) m_trigger += 1;
636     if (triggerEvent->path("HLT_Mu9") && triggerEvent->path("HLT_Mu9")->wasAccept()) m_trigger += 2;
637     if (triggerEvent->path("HLT_Mu11") && triggerEvent->path("HLT_Mu11")->wasAccept()) m_trigger += 4;
638 pivarski 1.12 if (triggerEvent->path("HLT_Mu15_v1") && triggerEvent->path("HLT_Mu15_v1")->wasAccept()) m_trigger += 8;
639 pivarski 1.10
640     // get some gen-level information
641     m_genLevel = -1;
642     if (m_fillGenLevel) {
643     bool any_missing_mcmatch = false;
644     bool any_decays_in_flight = false;
645     unsigned int bjets = 0;
646     unsigned int resonances = 0;
647     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
648     for (unsigned int i = 0; i < muJet->numberOfDaughters(); i++) {
649     const pat::Muon *muon = dynamic_cast<const pat::Muon*>(muJet->daughter(i));
650    
651     if (muon->genParticlesSize() == 0) any_missing_mcmatch = true;
652    
653     else if (muon->genParticle()->mother()) {
654     int mother = muon->genParticle()->mother()->pdgId();
655     if (abs(mother) == 13 && muon->genParticle()->mother()->mother()) muon->genParticle()->mother()->mother()->pdgId(); // don't let muons be their own mothers
656    
657     if (abs(mother) == 211 /* charged pion */ ||
658     abs(mother) == 321 /* charged kaon */ ||
659     (3000 <= abs(mother) && abs(mother) < 4000) /* strange baryon */ ) any_decays_in_flight = true;
660     }
661     }
662    
663     if (muJet->genParticlesSize() == 1) {
664     for (const reco::Candidate *particle = muJet->genParticle(); particle != NULL; particle = particle->mother()) {
665     int pdgId = particle->pdgId();
666     if (abs(pdgId) == 511 ||
667     abs(pdgId) == 521 ||
668     abs(pdgId) == 513 ||
669     abs(pdgId) == 523 ||
670     abs(pdgId) == 515 ||
671     abs(pdgId) == 525 ||
672     abs(pdgId) == 531 ||
673     abs(pdgId) == 533 ||
674     abs(pdgId) == 535 ||
675     abs(pdgId) == 541 ||
676     abs(pdgId) == 543 ||
677     abs(pdgId) == 545 ||
678     (5000 <= abs(pdgId) && abs(pdgId) < 6000)) {
679     bjets++;
680     break;
681     }
682     }
683    
684     if (muJet->genParticle()->pdgId() == 113 /* rho */ ||
685     muJet->genParticle()->pdgId() == 221 /* eta */ ||
686     muJet->genParticle()->pdgId() == 223 /* omega(782) */ ||
687     muJet->genParticle()->pdgId() == 333 /* phi(1020) */ ||
688     muJet->genParticle()->pdgId() == 443 /* J/psi */ ||
689     muJet->genParticle()->pdgId() == 100443 /* psi(2S) */ ) {
690     resonances++;
691     }
692     }
693     }
694    
695     m_genLevel = 0;
696     if (bjets > 0 && bjets == muJets->size()) m_genLevel += 1;
697     if (resonances > 0 && resonances == muJets->size()) m_genLevel += 2;
698     if (any_decays_in_flight) m_genLevel += 4;
699     if (any_missing_mcmatch) m_genLevel += 8;
700     }
701 pivarski 1.1
702 pivarski 1.7 ////////////////////////////////////////////////////////// alldimuons (special case for scaling to the Z)
703     if (muons->size() == 2) {
704     double total_energy = (*muons)[0].energy() + (*muons)[1].energy();
705     double total_px = (*muons)[0].px() + (*muons)[1].px();
706     double total_py = (*muons)[0].py() + (*muons)[1].py();
707     double total_pz = (*muons)[0].pz() + (*muons)[1].pz();
708     m_alldimuons_mass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
709     if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2) {
710     m_alldimuons_iso = (*muJets)[0].centralTrackIsolation();
711     }
712     else {
713     m_alldimuons_iso = -1.;
714     }
715     m_alldimuons->Fill();
716     }
717    
718 pivarski 1.1 ////////////////////////////////////////////////////////// lowdimuon:
719 pivarski 1.8 if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() < 80. && orphans->size() == 0) {
720 pivarski 1.1 pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
721    
722 pivarski 1.10 m_lowdimuon_containstrig = 0;
723     for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin(); iter != muontrig.end(); ++iter) {
724     if (muJet->muon(0)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
725     muJet->sameTrack(&*(muJet->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
726     m_lowdimuon_containstrig = 1;
727     }
728     if (muJet->muon(1)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
729     muJet->sameTrack(&*(muJet->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
730     m_lowdimuon_containstrig = 1;
731     }
732     }
733    
734 pivarski 1.1 // generator-level mass using matched genParticles (for resolution of fit peak)
735     m_lowdimuon_genmass = -1000.;
736 pivarski 1.5 if (muJet->muon(0)->genParticlesSize() == 1 && muJet->muon(1)->genParticlesSize() == 1) {
737 pivarski 1.4 const reco::GenParticle *mu0 = muJet->muon(0)->genParticle();
738     const reco::GenParticle *mu1 = muJet->muon(1)->genParticle();
739 pivarski 1.1
740     double total_energy = mu0->energy() + mu1->energy();
741     double total_px = mu0->px() + mu1->px();
742     double total_py = mu0->py() + mu1->py();
743     double total_pz = mu0->pz() + mu1->pz();
744     m_lowdimuon_genmass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
745     }
746    
747     m_lowdimuon_mass = muJet->mass();
748     m_lowdimuon_pt = muJet->pt();
749     m_lowdimuon_eta = muJet->eta();
750     m_lowdimuon_phi = muJet->phi();
751     m_lowdimuon_dr = muJet->dRmax();
752 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
753     m_lowdimuon_pluspx = muJet->daughter(0)->px();
754     m_lowdimuon_pluspy = muJet->daughter(0)->py();
755     m_lowdimuon_pluspz = muJet->daughter(0)->pz();
756     m_lowdimuon_minuspx = muJet->daughter(1)->px();
757     m_lowdimuon_minuspy = muJet->daughter(1)->py();
758     m_lowdimuon_minuspz = muJet->daughter(1)->pz();
759     }
760     else {
761     m_lowdimuon_pluspx = muJet->daughter(1)->px();
762     m_lowdimuon_pluspy = muJet->daughter(1)->py();
763     m_lowdimuon_pluspz = muJet->daughter(1)->pz();
764     m_lowdimuon_minuspx = muJet->daughter(0)->px();
765     m_lowdimuon_minuspy = muJet->daughter(0)->py();
766     m_lowdimuon_minuspz = muJet->daughter(0)->pz();
767     }
768 pivarski 1.1 m_lowdimuon_vprob = -1000.;
769     m_lowdimuon_vx = -1000.;
770     m_lowdimuon_vy = -1000.;
771     m_lowdimuon_vz = -1000.;
772    
773 pivarski 1.2 GlobalVector complus;
774 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0);
775     else complus = muJet->daughterCOM(1);
776 pivarski 1.2 m_lowdimuon_compluspx = complus.x();
777     m_lowdimuon_compluspy = complus.y();
778     m_lowdimuon_compluspz = complus.z();
779    
780 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0);
781     else complus = muJet->daughterCOMrot(1);
782     m_lowdimuon_comrotpluspx = complus.x();
783     m_lowdimuon_comrotpluspy = complus.y();
784     m_lowdimuon_comrotpluspz = complus.z();
785    
786 pivarski 1.2 // replace all values with vertex-updated values if vertex fitting succeeded
787 pivarski 1.1 if (muJet->vertexValid()) {
788     m_lowdimuon_mass = muJet->vertexMass();
789     m_lowdimuon_pt = muJet->vertexMomentum().perp();
790     m_lowdimuon_eta = muJet->vertexMomentum().eta();
791     m_lowdimuon_phi = muJet->vertexMomentum().phi();
792     m_lowdimuon_dr = muJet->dRmax(true);
793 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
794     m_lowdimuon_pluspx = muJet->vertexMomentum(0).x();
795     m_lowdimuon_pluspy = muJet->vertexMomentum(0).y();
796     m_lowdimuon_pluspz = muJet->vertexMomentum(0).z();
797     m_lowdimuon_minuspx = muJet->vertexMomentum(1).x();
798     m_lowdimuon_minuspy = muJet->vertexMomentum(1).y();
799     m_lowdimuon_minuspz = muJet->vertexMomentum(1).z();
800     }
801     else {
802     m_lowdimuon_pluspx = muJet->vertexMomentum(1).x();
803     m_lowdimuon_pluspy = muJet->vertexMomentum(1).y();
804     m_lowdimuon_pluspz = muJet->vertexMomentum(1).z();
805     m_lowdimuon_minuspx = muJet->vertexMomentum(0).x();
806     m_lowdimuon_minuspy = muJet->vertexMomentum(0).y();
807     m_lowdimuon_minuspz = muJet->vertexMomentum(0).z();
808     }
809 pivarski 1.1 m_lowdimuon_vprob = muJet->vertexProb();
810    
811     if (closestPrimaryVertex != primaryVertices->end()) {
812     m_lowdimuon_vx = muJet->vertexPoint().x() - closestPrimaryVertex->x();
813     m_lowdimuon_vy = muJet->vertexPoint().y() - closestPrimaryVertex->y();
814     m_lowdimuon_vz = muJet->vertexPoint().z() - closestPrimaryVertex->z();
815     }
816    
817 pivarski 1.2 GlobalVector complus;
818 pivarski 1.6 if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOM(0, true);
819     else complus = muJet->daughterCOM(1, true);
820 pivarski 1.2 m_lowdimuon_compluspx = complus.x();
821     m_lowdimuon_compluspy = complus.y();
822     m_lowdimuon_compluspz = complus.z();
823 pivarski 1.6
824     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0, true);
825     else complus = muJet->daughterCOMrot(1, true);
826     m_lowdimuon_comrotpluspx = complus.x();
827     m_lowdimuon_comrotpluspy = complus.y();
828     m_lowdimuon_comrotpluspz = complus.z();
829    
830 pivarski 1.1 } // end of replacements with quantities measured at the vertex
831    
832     m_lowdimuon_iso = muJet->centralTrackIsolation();
833    
834 pivarski 1.6 std::vector<reco::MuonChamberMatch> plusmatches, minusmatches;
835     if (muJet->daughter(0)->charge() > 0) {
836     m_lowdimuon_plusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
837     m_lowdimuon_minusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
838     plusmatches = muJet->muon(0)->matches();
839     minusmatches = muJet->muon(1)->matches();
840     m_lowdimuon_plushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
841     m_lowdimuon_minushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
842     m_lowdimuon_plusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
843     m_lowdimuon_minusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
844     }
845     else {
846     m_lowdimuon_plusmatches = muJet->muon(1)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
847     m_lowdimuon_minusmatches = muJet->muon(0)->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
848     plusmatches = muJet->muon(1)->matches();
849     minusmatches = muJet->muon(0)->matches();
850     m_lowdimuon_plushits = muJet->muon(1)->innerTrack()->numberOfValidHits();
851     m_lowdimuon_minushits = muJet->muon(0)->innerTrack()->numberOfValidHits();
852     m_lowdimuon_plusnormchi2 = muJet->muon(1)->innerTrack()->normalizedChi2();
853     m_lowdimuon_minusnormchi2 = muJet->muon(0)->innerTrack()->normalizedChi2();
854     }
855    
856     m_lowdimuon_plusst1x = -1000.;
857     m_lowdimuon_plusst1xsig = -1000.;
858     m_lowdimuon_plusst1dxdz = -1000.;
859     m_lowdimuon_plusst1dxdzsig = -1000.;
860     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = plusmatches.begin(); chamberMatch != plusmatches.end(); ++chamberMatch) {
861     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
862     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
863     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
864     if (chamberMatch->station() == 1) {
865     m_lowdimuon_plusst1x = (segmentMatch->x - chamberMatch->x);
866     m_lowdimuon_plusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
867     m_lowdimuon_plusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
868     m_lowdimuon_plusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
869     }
870     }
871     }
872     }
873    
874     m_lowdimuon_minusst1x = -1000.;
875     m_lowdimuon_minusst1xsig = -1000.;
876     m_lowdimuon_minusst1dxdz = -1000.;
877     m_lowdimuon_minusst1dxdzsig = -1000.;
878     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = minusmatches.begin(); chamberMatch != minusmatches.end(); ++chamberMatch) {
879     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
880     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
881     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
882     if (chamberMatch->station() == 1) {
883     m_lowdimuon_minusst1x = (segmentMatch->x - chamberMatch->x);
884     m_lowdimuon_minusst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
885     m_lowdimuon_minusst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
886     m_lowdimuon_minusst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
887     }
888     }
889     }
890     }
891    
892 pivarski 1.1 m_lowdimuon->Fill();
893     }
894    
895     ////////////////////////////////////////////////////////// highdimuon
896     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() > 80.) {
897     // ...
898     // m_highdimuon->Fill();
899     }
900    
901     ////////////////////////////////////////////////////////// quadmuon
902     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 4) {
903     // ...
904     // m_quadmuon->Fill();
905     }
906    
907 pivarski 1.8 ////////////////////////////////////////////////////////// dimuorphan
908     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && orphans->size() == 1) {
909     pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
910     pat::MuonCollection::const_iterator orphan = orphans->begin();
911    
912     m_dimuorphan_deltaphi = muJet->phi() - orphan->phi();
913     while (m_dimuorphan_deltaphi > M_PI) m_dimuorphan_deltaphi -= 2.*M_PI;
914     while (m_dimuorphan_deltaphi < -M_PI) m_dimuorphan_deltaphi += 2.*M_PI;
915    
916     m_dimuorphan_orphanpt = orphan->pt();
917     m_dimuorphan_orphaneta = orphan->eta();
918     m_dimuorphan_orphanphi = orphan->phi();
919 pivarski 1.9 m_dimuorphan_orphanisglobal = (orphan->isGlobalMuon() ? 1 : 0);
920     m_dimuorphan_orphanmatches = orphan->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
921     m_dimuorphan_orphanhits = (orphan->innerTrack().isAvailable() ? orphan->innerTrack()->numberOfValidHits(): -1);
922     m_dimuorphan_orphanchi2 = (orphan->innerTrack().isAvailable() ? orphan->innerTrack()->normalizedChi2(): -1.);
923     m_dimuorphan_orphanst1x = -1000.;
924     m_dimuorphan_orphanst1xsig = -1000.;
925     m_dimuorphan_orphanst1dxdz = -1000.;
926     m_dimuorphan_orphanst1dxdzsig = -1000.;
927     for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = orphan->matches().begin(); chamberMatch != orphan->matches().end(); ++chamberMatch) {
928     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin(); segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch) {
929     if (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
930     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
931     if (chamberMatch->station() == 1) {
932     m_dimuorphan_orphanst1x = (segmentMatch->x - chamberMatch->x);
933     m_dimuorphan_orphanst1xsig = (segmentMatch->x - chamberMatch->x)/sqrt(pow(segmentMatch->xErr, 2) + pow(chamberMatch->xErr, 2));
934     m_dimuorphan_orphanst1dxdz = (segmentMatch->dXdZ - chamberMatch->dXdZ);
935     m_dimuorphan_orphanst1dxdzsig = (segmentMatch->dXdZ - chamberMatch->dXdZ)/sqrt(pow(segmentMatch->dXdZErr, 2) + pow(chamberMatch->dXdZErr, 2));
936     }
937     }
938     }
939     }
940 pivarski 1.8
941 pivarski 1.10 m_dimuorphan_containstrig = 0;
942     for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin(); iter != muontrig.end(); ++iter) {
943     if (orphan->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
944     muJet->sameTrack(&*(orphan->innerTrack()), &*((*iter)->innerTrack()))) {
945     m_dimuorphan_containstrig = 1;
946     }
947 pivarski 1.8 }
948 pivarski 1.10
949 pivarski 1.8 m_dimuorphan_mass = muJet->mass();
950     m_dimuorphan_pt = muJet->pt();
951     m_dimuorphan_eta = muJet->eta();
952     m_dimuorphan_phi = muJet->phi();
953     m_dimuorphan_dr = muJet->dRmax();
954     m_dimuorphan_vprob = -1.;
955     m_dimuorphan_lxy = -1000.;
956     m_dimuorphan_lxyz = -1000.;
957     m_dimuorphan_iso = muJet->centralTrackIsolation();
958    
959     if (muJet->vertexValid()) {
960     m_dimuorphan_deltaphi = muJet->vertexMomentum().phi() - orphan->phi();
961     while (m_dimuorphan_deltaphi > M_PI) m_dimuorphan_deltaphi -= 2.*M_PI;
962     while (m_dimuorphan_deltaphi < -M_PI) m_dimuorphan_deltaphi += 2.*M_PI;
963    
964     m_dimuorphan_mass = muJet->vertexMass();
965     m_dimuorphan_pt = muJet->vertexMomentum().perp();
966     m_dimuorphan_eta = muJet->vertexMomentum().eta();
967     m_dimuorphan_phi = muJet->vertexMomentum().phi();
968     m_dimuorphan_dr = muJet->dRmax(true);
969     m_dimuorphan_vprob = muJet->vertexProb();
970 pivarski 1.10 if (closestPrimaryVertex != primaryVertices->end()) {
971     m_dimuorphan_lxy = muJet->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
972     m_dimuorphan_lxyz = muJet->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
973     }
974 pivarski 1.8 }
975    
976     m_dimuorphan->Fill();
977     }
978    
979 pivarski 1.1 ////////////////////////////////////////////////////////// dimudimu
980     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 2) {
981 pivarski 1.10 const pat::MultiMuon *muJet0 = &((*muJets)[0]);
982     const pat::MultiMuon *muJet1 = &((*muJets)[1]);
983    
984     bool muJet0_canBeC = false;
985     bool muJet1_canBeC = false;
986     for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin(); iter != muontrig.end(); ++iter) {
987     if (muJet0->muon(0)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
988     muJet0->sameTrack(&*(muJet0->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
989     muJet0_canBeC = true;
990     }
991     if (muJet0->muon(1)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
992     muJet0->sameTrack(&*(muJet0->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
993     muJet0_canBeC = true;
994     }
995    
996     if (muJet1->muon(0)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
997     muJet1->sameTrack(&*(muJet1->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
998     muJet1_canBeC = true;
999 pivarski 1.8 }
1000 pivarski 1.10 if (muJet1->muon(1)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
1001     muJet1->sameTrack(&*(muJet1->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
1002     muJet1_canBeC = true;
1003 pivarski 1.8 }
1004     }
1005 pivarski 1.10
1006     const pat::MultiMuon *muJetC = NULL;
1007     const pat::MultiMuon *muJetF = NULL;
1008     if (muJet0_canBeC && muJet1_canBeC) {
1009     if (m_trandom3.Integer(2) == 0) {
1010 pivarski 1.8 muJetC = muJet0;
1011     muJetF = muJet1;
1012     }
1013     else {
1014     muJetC = muJet1;
1015     muJetF = muJet0;
1016     }
1017     }
1018 pivarski 1.10 else if (muJet0_canBeC) {
1019     muJetC = muJet0;
1020     muJetF = muJet1;
1021     }
1022     else if (muJet1_canBeC) {
1023     muJetC = muJet1;
1024     muJetF = muJet0;
1025     }
1026 pivarski 1.8
1027 pivarski 1.10 if (muJetC != NULL && muJetF != NULL) {
1028     double total_energy = muJetC->energy() + muJetF->energy();
1029     double total_px = muJetC->px() + muJetF->px();
1030     double total_py = muJetC->py() + muJetF->py();
1031     double total_pz = muJetC->pz() + muJetF->pz();
1032 pivarski 1.7
1033     m_dimudimu_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
1034     m_dimudimu_wholept = sqrt(total_px*total_px + total_py*total_py);
1035    
1036 pivarski 1.10 m_dimudimu_deltaphi = muJetC->phi() - muJetF->phi();
1037 pivarski 1.7 while (m_dimudimu_deltaphi > M_PI) m_dimudimu_deltaphi -= 2.*M_PI;
1038     while (m_dimudimu_deltaphi < -M_PI) m_dimudimu_deltaphi += 2.*M_PI;
1039    
1040 pivarski 1.10 m_dimudimu_massC = muJetC->mass();
1041     m_dimudimu_ptC = muJetC->pt();
1042     m_dimudimu_etaC = muJetC->eta();
1043     m_dimudimu_phiC = muJetC->phi();
1044     m_dimudimu_drC = muJetC->dRmax();
1045     m_dimudimu_vprobC = -1.;
1046     m_dimudimu_lxyC = -1000.;
1047     m_dimudimu_lxyzC = -1000.;
1048     m_dimudimu_isoC = muJetC->centralTrackIsolation();
1049    
1050     m_dimudimu_massF = muJetF->mass();
1051     m_dimudimu_ptF = muJetF->pt();
1052     m_dimudimu_etaF = muJetF->eta();
1053     m_dimudimu_phiF = muJetF->phi();
1054     m_dimudimu_drF = muJetF->dRmax();
1055     m_dimudimu_vprobF = -1.;
1056     m_dimudimu_lxyF = -1000.;
1057     m_dimudimu_lxyzF = -1000.;
1058     m_dimudimu_isoF = muJetF->centralTrackIsolation();
1059    
1060     if (muJetC->vertexValid() && muJetF->vertexValid()) {
1061     total_energy = muJetC->vertexP4().energy() + muJetF->vertexP4().energy();
1062     total_px = muJetC->vertexMomentum().x() + muJetF->vertexMomentum().x();
1063     total_py = muJetC->vertexMomentum().y() + muJetF->vertexMomentum().y();
1064     total_pz = muJetC->vertexMomentum().z() + muJetF->vertexMomentum().z();
1065    
1066     m_dimudimu_wholemass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
1067     m_dimudimu_wholept = sqrt(total_px*total_px + total_py*total_py);
1068    
1069     m_dimudimu_deltaphi = muJetC->vertexMomentum().phi() - muJetF->vertexMomentum().phi();
1070     while (m_dimudimu_deltaphi > M_PI) m_dimudimu_deltaphi -= 2.*M_PI;
1071     while (m_dimudimu_deltaphi < -M_PI) m_dimudimu_deltaphi += 2.*M_PI;
1072    
1073     m_dimudimu_massC = muJetC->vertexMass();
1074     m_dimudimu_ptC = muJetC->vertexMomentum().perp();
1075     m_dimudimu_etaC = muJetC->vertexMomentum().eta();
1076     m_dimudimu_phiC = muJetC->vertexMomentum().phi();
1077     m_dimudimu_drC = muJetC->dRmax(true);
1078     m_dimudimu_vprobC = muJetC->vertexProb();
1079     if (closestPrimaryVertex != primaryVertices->end()) {
1080     m_dimudimu_lxyC = muJetC->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
1081     m_dimudimu_lxyzC = muJetC->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
1082     }
1083    
1084     m_dimudimu_massF = muJetF->vertexMass();
1085     m_dimudimu_ptF = muJetF->vertexMomentum().perp();
1086     m_dimudimu_etaF = muJetF->vertexMomentum().eta();
1087     m_dimudimu_phiF = muJetF->vertexMomentum().phi();
1088     m_dimudimu_drF = muJetF->dRmax(true);
1089     m_dimudimu_vprobF = muJetF->vertexProb();
1090     if (closestPrimaryVertex != primaryVertices->end()) {
1091     m_dimudimu_lxyF = muJetF->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
1092     m_dimudimu_lxyzF = muJetF->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
1093     }
1094     }
1095    
1096     m_dimudimu->Fill();
1097 pivarski 1.7 }
1098 pivarski 1.1 }
1099    
1100     ////////////////////////////////////////////////////////// dimuquadmu
1101     else if (muJets->size() == 2 && (
1102     ((*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 4) ||
1103     ((*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 2)
1104     )) {
1105     // ...
1106     // m_dimuquadmu->Fill();
1107     }
1108    
1109     ////////////////////////////////////////////////////////// quadmuquadmu
1110     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 4) {
1111     // ...
1112     // m_quadmuquadmu->Fill();
1113     }
1114 pivarski 1.10
1115     ////////////////////////////////////////////////////////// mujetplustracks
1116     if (muJetPlusTracks->size() == 1 && ((*muJetPlusTracks)[0].userInt("numReal") == 2 || (*muJetPlusTracks)[0].userInt("numReal") == 3) && (*muJetPlusTracks)[0].userInt("numFakes") > 0) {
1117     const pat::MultiMuon *muJet = &((*muJetPlusTracks)[0]);
1118    
1119     m_mujetplustracks_containstrig = 0;
1120     for (std::vector<pat::MuonCollection::const_iterator>::const_iterator iter = muontrig.begin(); iter != muontrig.end(); ++iter) {
1121     if (muJet->muon(0)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
1122     muJet->sameTrack(&*(muJet->muon(0)->innerTrack()), &*((*iter)->innerTrack()))) {
1123     m_mujetplustracks_containstrig = 1;
1124     }
1125     if (muJet->muon(1)->innerTrack().isAvailable() && (*iter)->innerTrack().isAvailable() &&
1126     muJet->sameTrack(&*(muJet->muon(1)->innerTrack()), &*((*iter)->innerTrack()))) {
1127     m_mujetplustracks_containstrig = 1;
1128     }
1129     }
1130    
1131     m_mujetplustracks_muons = muJet->userInt("numReal");
1132     m_mujetplustracks_fakes = muJet->userInt("numFakes");
1133     m_mujetplustracks_mass = muJet->mass();
1134 pivarski 1.11 m_mujetplustracks_charge = muJet->charge();
1135     m_mujetplustracks_massa = -1.;
1136     m_mujetplustracks_massb = -1.;
1137     std::vector<double> masses = muJet->consistentPairMasses();
1138     if (masses.size() >= 1) m_mujetplustracks_massa = masses[0];
1139     if (masses.size() >= 2) m_mujetplustracks_massb = masses[1];
1140    
1141 pivarski 1.10 m_mujetplustracks_pt = muJet->pt();
1142     m_mujetplustracks_eta = muJet->eta();
1143     m_mujetplustracks_phi = muJet->phi();
1144     m_mujetplustracks_dr = muJet->dRmax();
1145     m_mujetplustracks_vprob = -1.;
1146     m_mujetplustracks_lxy = -1000.;
1147     m_mujetplustracks_lxyz = -1000.;
1148     m_mujetplustracks_iso = muJet->centralTrackIsolation();
1149    
1150     m_mujetplustracks_mupt1 = -1.;
1151     m_mujetplustracks_mupt2 = -1.;
1152     m_mujetplustracks_mupt3 = -1.;
1153     m_mujetplustracks_fakept1 = -1.;
1154     m_mujetplustracks_fakept2 = -1.;
1155     for (unsigned int i = 0; i < muJet->numberOfDaughters(); i++) {
1156     if (muJet->muon(i)->matches().size() > 0) {
1157     if (muJet->muon(i)->pt() > m_mujetplustracks_mupt1) {
1158     m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
1159     m_mujetplustracks_mupt2 = m_mujetplustracks_mupt1;
1160     m_mujetplustracks_mupt1 = muJet->muon(i)->pt();
1161     }
1162     else if (muJet->muon(i)->pt() > m_mujetplustracks_mupt2) {
1163     m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
1164     m_mujetplustracks_mupt2 = muJet->muon(i)->pt();
1165     }
1166     else if (muJet->muon(i)->pt() > m_mujetplustracks_mupt3) {
1167     m_mujetplustracks_mupt3 = muJet->muon(i)->pt();
1168     }
1169     }
1170     else {
1171     if (muJet->muon(i)->pt() > m_mujetplustracks_fakept1) {
1172     m_mujetplustracks_fakept2 = m_mujetplustracks_fakept1;
1173     m_mujetplustracks_fakept1 = muJet->muon(i)->pt();
1174     }
1175     else if (muJet->muon(i)->pt() > m_mujetplustracks_fakept1) {
1176     m_mujetplustracks_fakept2 = muJet->muon(i)->pt();
1177     }
1178     }
1179     }
1180    
1181     if (muJet->vertexValid()) {
1182     m_mujetplustracks_mass = muJet->vertexMass();
1183 pivarski 1.11 m_mujetplustracks_massa = -1.;
1184     m_mujetplustracks_massb = -1.;
1185     std::vector<double> masses = muJet->consistentPairMasses(true);
1186     if (masses.size() >= 1) m_mujetplustracks_massa = masses[0];
1187     if (masses.size() >= 2) m_mujetplustracks_massb = masses[1];
1188    
1189 pivarski 1.10 m_mujetplustracks_pt = muJet->vertexMomentum().perp();
1190     m_mujetplustracks_eta = muJet->vertexMomentum().eta();
1191     m_mujetplustracks_phi = muJet->vertexMomentum().phi();
1192     m_mujetplustracks_dr = muJet->dRmax(true);
1193     m_mujetplustracks_vprob = muJet->vertexProb();
1194     if (closestPrimaryVertex != primaryVertices->end()) {
1195     m_mujetplustracks_lxy = muJet->lxy(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
1196     m_mujetplustracks_lxyz = muJet->lxyz(GlobalPoint(closestPrimaryVertex->x(), closestPrimaryVertex->y(), closestPrimaryVertex->z()));
1197     }
1198    
1199     m_mujetplustracks_mupt1 = -1.;
1200     m_mujetplustracks_mupt2 = -1.;
1201     m_mujetplustracks_mupt3 = -1.;
1202     m_mujetplustracks_fakept1 = -1.;
1203     m_mujetplustracks_fakept2 = -1.;
1204     for (unsigned int i = 0; i < muJet->numberOfDaughters(); i++) {
1205     if (muJet->muon(i)->matches().size() > 0) {
1206     if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_mupt1) {
1207     m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
1208     m_mujetplustracks_mupt2 = m_mujetplustracks_mupt1;
1209     m_mujetplustracks_mupt1 = muJet->vertexMomentum(i).perp();
1210     }
1211     else if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_mupt2) {
1212     m_mujetplustracks_mupt3 = m_mujetplustracks_mupt2;
1213     m_mujetplustracks_mupt2 = muJet->vertexMomentum(i).perp();
1214     }
1215     else if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_mupt3) {
1216     m_mujetplustracks_mupt3 = muJet->vertexMomentum(i).perp();
1217     }
1218     }
1219     else {
1220     if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_fakept1) {
1221     m_mujetplustracks_fakept2 = m_mujetplustracks_fakept1;
1222     m_mujetplustracks_fakept1 = muJet->vertexMomentum(i).perp();
1223     }
1224     else if (muJet->vertexMomentum(i).perp() > m_mujetplustracks_fakept1) {
1225     m_mujetplustracks_fakept2 = muJet->vertexMomentum(i).perp();
1226     }
1227     }
1228     }
1229     }
1230    
1231     m_mujetplustracks->Fill();
1232     }
1233 pivarski 1.1 }
1234    
1235    
1236     // ------------ method called once each job just before starting event loop ------------
1237     void
1238     FitNtuple::beginJob()
1239     {
1240     }
1241    
1242     // ------------ method called once each job just after ending the event loop ------------
1243     void
1244     FitNtuple::endJob() {
1245     }
1246    
1247     //define this as a plug-in
1248     DEFINE_FWK_MODULE(FitNtuple);