ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/FitNtuple/src/FitNtuple.cc
Revision: 1.4
Committed: Sat Nov 20 17:16:11 2010 UTC (14 years, 5 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-11-20-c
Changes since 1.3: +4 -4 lines
Log Message:
use MuJet's muon() method to do the Candidate -> pat::Muon cast (I had forgotten about that feature)

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.4 // $Id: FitNtuple.cc,v 1.3 2010/11/20 15:43:06 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    
44     #include "CommonTools/UtilAlgos/interface/TFileService.h"
45     #include "FWCore/ServiceRegistry/interface/Service.h"
46     #include "TTree.h"
47    
48     //
49     // class declaration
50     //
51    
52     class FitNtuple : public edm::EDAnalyzer {
53     public:
54     explicit FitNtuple(const edm::ParameterSet&);
55     ~FitNtuple();
56    
57    
58     private:
59     virtual void beginJob() ;
60     virtual void analyze(const edm::Event&, const edm::EventSetup&);
61     virtual void endJob() ;
62    
63     // ----------member data ---------------------------
64    
65     // input parameters
66     edm::InputTag m_src;
67     bool m_getQscale;
68     bool m_getAlternating;
69    
70     // all ntuples need these variables
71     Float_t m_qscale; // also known as "pt-hat" ($\hat{p}_T$)
72     Int_t m_run; // run and event number so that we can look at
73     Int_t m_event; // interesting cases in the event display
74     Int_t m_trigger; // satisfied trigger? 0 = failed all
75     // 1 = passed HLT_Mu5
76     // 2 = passed HLT_Mu9
77     // 4 = passed HLT_Mu11
78 pivarski 1.3 Float_t m_muon1pt;
79     Float_t m_muon1eta;
80     Float_t m_muon2pt;
81     Float_t m_muon2eta;
82     Float_t m_muon3pt;
83     Float_t m_muon3eta;
84     Float_t m_muon4pt;
85     Float_t m_muon4eta;
86     Float_t m_muontrigpt;
87     Float_t m_muontrigeta;
88 pivarski 1.1
89     // control sample "lowdimuon" (single mu-jet with two muons)
90     TTree *m_lowdimuon;
91     Float_t m_lowdimuon_genmass;
92     Float_t m_lowdimuon_mass;
93     Float_t m_lowdimuon_pt;
94     Float_t m_lowdimuon_eta;
95     Float_t m_lowdimuon_phi;
96     Float_t m_lowdimuon_dr;
97     Float_t m_lowdimuon_pluspx;
98     Float_t m_lowdimuon_pluspy;
99     Float_t m_lowdimuon_pluspz;
100     Float_t m_lowdimuon_minuspx;
101     Float_t m_lowdimuon_minuspy;
102     Float_t m_lowdimuon_minuspz;
103     Float_t m_lowdimuon_vprob;
104     Float_t m_lowdimuon_vx;
105     Float_t m_lowdimuon_vy;
106     Float_t m_lowdimuon_vz;
107     Float_t m_lowdimuon_iso;
108 pivarski 1.2 Float_t m_lowdimuon_compluspx;
109     Float_t m_lowdimuon_compluspy;
110     Float_t m_lowdimuon_compluspz;
111 pivarski 1.1
112     // signal region "highdimuon" (single mu-jet with two muons)
113     TTree *m_highdimuon;
114    
115     // control sample "muplustrack" (single mu-jet with three muons and the closest track)
116     TTree *m_muplustrack;
117    
118     // signal region "quadmuon" (single mu-jet with four muons)
119     TTree *m_quadmuon;
120    
121     // signal region "dimudimu" (two mu-jets with two muons each)
122     TTree *m_dimudimu;
123    
124     // signal region "dimuquadmu" (one two-muon jet and one four-muon jet)
125     TTree *m_dimuquadmu;
126    
127     // signal region "quadmuquadmu" (two mu-jets with four muons each)
128     TTree *m_quadmuquadmu;
129    
130     };
131    
132     //
133     // constants, enums and typedefs
134     //
135    
136     //
137     // static data member definitions
138     //
139    
140     //
141     // constructors and destructor
142     //
143     FitNtuple::FitNtuple(const edm::ParameterSet& iConfig)
144     : m_src(iConfig.getParameter<edm::InputTag>("src"))
145     , m_getQscale(iConfig.getParameter<bool>("getQscale"))
146     , m_getAlternating(iConfig.getParameter<bool>("getAlternating"))
147     {
148     //now do what ever initialization is needed
149     edm::Service<TFileService> tFile;
150    
151     m_lowdimuon = tFile->make<TTree>("lowdimuon", "lowdimuon");;
152     m_lowdimuon->Branch("qscale", &m_qscale, "qscale/F");
153     m_lowdimuon->Branch("run", &m_run, "run/I");
154     m_lowdimuon->Branch("event", &m_event, "event/I");
155     m_lowdimuon->Branch("trigger", &m_trigger, "trigger/I");
156 pivarski 1.3 m_lowdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
157     m_lowdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
158     m_lowdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
159     m_lowdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
160     m_lowdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
161     m_lowdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
162 pivarski 1.1 m_lowdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
163     m_lowdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
164     m_lowdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
165     m_lowdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
166     m_lowdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
167     m_lowdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
168     m_lowdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
169     m_lowdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
170     m_lowdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
171     m_lowdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
172     m_lowdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
173     m_lowdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
174     m_lowdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
175     m_lowdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
176     m_lowdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
177     m_lowdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
178     m_lowdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
179 pivarski 1.2 m_lowdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
180     m_lowdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
181     m_lowdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
182 pivarski 1.1
183     m_highdimuon = tFile->make<TTree>("highdimuon", "highdimuon");;
184     m_highdimuon->Branch("qscale", &m_qscale, "qscale/F");
185     m_highdimuon->Branch("run", &m_run, "run/I");
186     m_highdimuon->Branch("event", &m_event, "event/I");
187     m_highdimuon->Branch("trigger", &m_trigger, "trigger/I");
188 pivarski 1.3 m_highdimuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
189     m_highdimuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
190     m_highdimuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
191     m_highdimuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
192     m_highdimuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
193     m_highdimuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
194 pivarski 1.1
195     m_muplustrack = tFile->make<TTree>("muplustrack", "muplustrack");;
196     m_muplustrack->Branch("qscale", &m_qscale, "qscale/F");
197     m_muplustrack->Branch("run", &m_run, "run/I");
198     m_muplustrack->Branch("event", &m_event, "event/I");
199     m_muplustrack->Branch("trigger", &m_trigger, "trigger/I");
200 pivarski 1.3 m_muplustrack->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
201     m_muplustrack->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
202     m_muplustrack->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
203     m_muplustrack->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
204     m_muplustrack->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
205     m_muplustrack->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
206     m_muplustrack->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
207     m_muplustrack->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
208     m_muplustrack->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
209     m_muplustrack->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
210 pivarski 1.1
211     m_quadmuon = tFile->make<TTree>("quadmuon", "quadmuon");;
212     m_quadmuon->Branch("qscale", &m_qscale, "qscale/F");
213     m_quadmuon->Branch("run", &m_run, "run/I");
214     m_quadmuon->Branch("event", &m_event, "event/I");
215     m_quadmuon->Branch("trigger", &m_trigger, "trigger/I");
216 pivarski 1.3 m_quadmuon->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
217     m_quadmuon->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
218     m_quadmuon->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
219     m_quadmuon->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
220     m_quadmuon->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
221     m_quadmuon->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
222     m_quadmuon->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
223     m_quadmuon->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
224     m_quadmuon->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
225     m_quadmuon->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
226 pivarski 1.1
227     m_dimudimu = tFile->make<TTree>("dimudimu", "dimudimu");;
228     m_dimudimu->Branch("qscale", &m_qscale, "qscale/F");
229     m_dimudimu->Branch("run", &m_run, "run/I");
230     m_dimudimu->Branch("event", &m_event, "event/I");
231     m_dimudimu->Branch("trigger", &m_trigger, "trigger/I");
232 pivarski 1.3 m_dimudimu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
233     m_dimudimu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
234     m_dimudimu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
235     m_dimudimu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
236     m_dimudimu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
237     m_dimudimu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
238     m_dimudimu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
239     m_dimudimu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
240     m_dimudimu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
241     m_dimudimu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
242 pivarski 1.1
243     m_dimuquadmu = tFile->make<TTree>("dimuquadmu", "dimuquadmu");;
244     m_dimuquadmu->Branch("qscale", &m_qscale, "qscale/F");
245     m_dimuquadmu->Branch("run", &m_run, "run/I");
246     m_dimuquadmu->Branch("event", &m_event, "event/I");
247     m_dimuquadmu->Branch("trigger", &m_trigger, "trigger/I");
248 pivarski 1.3 m_dimuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
249     m_dimuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
250     m_dimuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
251     m_dimuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
252     m_dimuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
253     m_dimuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
254     m_dimuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
255     m_dimuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
256     m_dimuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
257     m_dimuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
258 pivarski 1.1
259     m_quadmuquadmu = tFile->make<TTree>("quadmuquadmu", "quadmuquadmu");;
260     m_quadmuquadmu->Branch("qscale", &m_qscale, "qscale/F");
261     m_quadmuquadmu->Branch("run", &m_run, "run/I");
262     m_quadmuquadmu->Branch("event", &m_event, "event/I");
263     m_quadmuquadmu->Branch("trigger", &m_trigger, "trigger/I");
264 pivarski 1.3 m_quadmuquadmu->Branch("muon1pt", &m_muon1pt, "muon1pt/F");
265     m_quadmuquadmu->Branch("muon1eta", &m_muon1eta, "muon1eta/F");
266     m_quadmuquadmu->Branch("muon2pt", &m_muon2pt, "muon2pt/F");
267     m_quadmuquadmu->Branch("muon2eta", &m_muon2eta, "muon2eta/F");
268     m_quadmuquadmu->Branch("muon3pt", &m_muon3pt, "muon3pt/F");
269     m_quadmuquadmu->Branch("muon3eta", &m_muon3eta, "muon3eta/F");
270     m_quadmuquadmu->Branch("muon4pt", &m_muon4pt, "muon4pt/F");
271     m_quadmuquadmu->Branch("muon4eta", &m_muon4eta, "muon4eta/F");
272     m_quadmuquadmu->Branch("muontrigpt", &m_muontrigpt, "muontrigpt/F");
273     m_quadmuquadmu->Branch("muontrigeta", &m_muontrigeta, "muontrigeta/F");
274 pivarski 1.1 }
275    
276    
277     FitNtuple::~FitNtuple()
278     {
279    
280     // do anything here that needs to be done at desctruction time
281     // (e.g. close files, deallocate resources etc.)
282    
283     }
284    
285    
286     //
287     // member functions
288     //
289    
290     // ------------ method called to for each event ------------
291     void
292     FitNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
293     {
294     // the QCD scale only exists for certain kinds of Monte Carlo, but if it's available, get it
295     m_qscale = 0.;
296     if (m_getQscale) {
297     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
298     iEvent.getByLabel("generator", genEventInfoProduct);
299     m_qscale = genEventInfoProduct->qScale();
300     }
301    
302     // the pair-gun Monte Carlos have a (never used) feature called
303     // alteration; the real events are the ones in which particleNumber == 0
304     bool alternating = true;
305     if (m_getAlternating) {
306     edm::Handle<unsigned int> particleNumber;
307     iEvent.getByLabel("generator", "particleNumber", particleNumber);
308     alternating = (*particleNumber == 0);
309     }
310     if (!alternating) return;
311    
312     // get the run and event number
313     m_run = iEvent.id().run();
314     m_event = iEvent.id().event();
315    
316     // mu-jets (muons grouped by mass and vertex compatibility)
317     edm::Handle<pat::MultiMuonCollection> muJets;
318     iEvent.getByLabel(m_src, muJets);
319    
320     // // orphans (muons not found in any group)
321     // edm::Handle<pat::MuonCollection> orphans;
322     // iEvent.getByLabel(m_src, "Orphans", orphans);
323    
324 pivarski 1.3 // find the top four muons in the event (-1000. if not found)
325     edm::Handle<pat::MuonCollection> muons;
326     iEvent.getByLabel("cleanPatMuons", muons);
327     m_muon1pt = -1000.; m_muon1eta = -1000.;
328     m_muon2pt = -1000.; m_muon2eta = -1000.;
329     m_muon3pt = -1000.; m_muon3eta = -1000.;
330     m_muon4pt = -1000.; m_muon4eta = -1000.;
331     m_muontrigpt = -1000.; m_muontrigeta = -1000.;
332     for (pat::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
333     if (fabs(muon->eta()) < 2.4 && muon->isTrackerMuon() && muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) >= 2 && muon->innerTrack()->numberOfValidHits() >= 8 && muon->innerTrack()->normalizedChi2() < 4.) {
334     if (muon->pt() > m_muon1pt) {
335     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
336     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
337     m_muon2pt = m_muon1pt; m_muon2eta = m_muon1eta;
338     m_muon1pt = muon->pt(); m_muon1eta = muon->eta();
339     }
340     else if (muon->pt() > m_muon2pt) {
341     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
342     m_muon3pt = m_muon2pt; m_muon3eta = m_muon2eta;
343     m_muon2pt = muon->pt(); m_muon2eta = muon->eta();
344     }
345     else if (muon->pt() > m_muon3pt) {
346     m_muon4pt = m_muon3pt; m_muon4eta = m_muon3eta;
347     m_muon3pt = muon->pt(); m_muon3eta = muon->eta();
348     }
349     else if (muon->pt() > m_muon4pt) {
350     m_muon4pt = muon->pt(); m_muon4eta = muon->eta();
351     }
352    
353     // special muon within a more limited |eta| range, to guarantee the trigger
354     if (fabs(muon->eta()) < 1.0) {
355     if (muon->pt() > m_muontrigpt) {
356     m_muontrigpt = muon->pt(); m_muontrigeta = muon->eta();
357     }
358     }
359     }
360     }
361 pivarski 1.1
362     // // all tracker-tracks
363     // edm::Handle<reco::TrackCollection> tracks;
364     // iEvent.getByLabel("generalTracks", tracks);
365    
366     // // generator-level 4-vectors
367     // edm::Handle<reco::GenParticleCollection> genParticles;
368     // iEvent.getByLabel("genParticles", genParticles);
369    
370     // find the closest primary vertex (in Z) to the first muJet with a valid vertex
371     edm::Handle<reco::VertexCollection> primaryVertices;
372     iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
373     reco::VertexCollection::const_iterator closestPrimaryVertex = primaryVertices->end();
374     if (muJets->size() > 0) {
375     pat::MultiMuonCollection::const_iterator muJet0 = muJets->end();
376     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
377     if (muJet->vertexValid()) {
378     muJet0 = muJet;
379     break;
380     }
381     }
382    
383     if (muJet0 != muJets->end()) {
384     for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin(); vertex != primaryVertices->end(); ++vertex) {
385     if (vertex->isValid() && !vertex->isFake() && vertex->tracksSize() > 3 && fabs(vertex->z()) < 24.) {
386     if (closestPrimaryVertex == primaryVertices->end() || fabs(vertex->z() - muJet0->vertexPoint().z()) < fabs(closestPrimaryVertex->z() - muJet0->vertexPoint().z())) {
387     closestPrimaryVertex = vertex;
388     }
389     } // end vertex quality cuts
390     } // end loop over primary vertices
391     } // end if muJet0 exists
392     } // end if muJets->size > 0
393    
394     // find out which trigger bits were fired
395     edm::Handle<pat::TriggerEvent> triggerEvent;
396     iEvent.getByLabel("patTriggerEvent", triggerEvent);
397     m_trigger = 0;
398     if (triggerEvent->path("HLT_Mu5") && triggerEvent->path("HLT_Mu5")->wasAccept()) m_trigger += 1;
399     if (triggerEvent->path("HLT_Mu9") && triggerEvent->path("HLT_Mu9")->wasAccept()) m_trigger += 2;
400     if (triggerEvent->path("HLT_Mu11") && triggerEvent->path("HLT_Mu11")->wasAccept()) m_trigger += 4;
401    
402     ////////////////////////////////////////////////////////// lowdimuon:
403     if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() < 80.) {
404     pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
405    
406     // generator-level mass using matched genParticles (for resolution of fit peak)
407     m_lowdimuon_genmass = -1000.;
408 pivarski 1.4 if (muJet->muon(0)->genParticlesSize() == 1 || muJet->muon(1)->genParticlesSize() == 1) {
409     const reco::GenParticle *mu0 = muJet->muon(0)->genParticle();
410     const reco::GenParticle *mu1 = muJet->muon(1)->genParticle();
411 pivarski 1.1
412     double total_energy = mu0->energy() + mu1->energy();
413     double total_px = mu0->px() + mu1->px();
414     double total_py = mu0->py() + mu1->py();
415     double total_pz = mu0->pz() + mu1->pz();
416     m_lowdimuon_genmass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
417     }
418    
419     m_lowdimuon_mass = muJet->mass();
420     m_lowdimuon_pt = muJet->pt();
421     m_lowdimuon_eta = muJet->eta();
422     m_lowdimuon_phi = muJet->phi();
423     m_lowdimuon_dr = muJet->dRmax();
424 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
425     m_lowdimuon_pluspx = muJet->daughter(0)->px();
426     m_lowdimuon_pluspy = muJet->daughter(0)->py();
427     m_lowdimuon_pluspz = muJet->daughter(0)->pz();
428     m_lowdimuon_minuspx = muJet->daughter(1)->px();
429     m_lowdimuon_minuspy = muJet->daughter(1)->py();
430     m_lowdimuon_minuspz = muJet->daughter(1)->pz();
431     }
432     else {
433     m_lowdimuon_pluspx = muJet->daughter(1)->px();
434     m_lowdimuon_pluspy = muJet->daughter(1)->py();
435     m_lowdimuon_pluspz = muJet->daughter(1)->pz();
436     m_lowdimuon_minuspx = muJet->daughter(0)->px();
437     m_lowdimuon_minuspy = muJet->daughter(0)->py();
438     m_lowdimuon_minuspz = muJet->daughter(0)->pz();
439     }
440 pivarski 1.1 m_lowdimuon_vprob = -1000.;
441     m_lowdimuon_vx = -1000.;
442     m_lowdimuon_vy = -1000.;
443     m_lowdimuon_vz = -1000.;
444    
445 pivarski 1.2 GlobalVector complus;
446     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0);
447     else complus = muJet->daughterCOMrot(1);
448     m_lowdimuon_compluspx = complus.x();
449     m_lowdimuon_compluspy = complus.y();
450     m_lowdimuon_compluspz = complus.z();
451    
452     // replace all values with vertex-updated values if vertex fitting succeeded
453 pivarski 1.1 if (muJet->vertexValid()) {
454     m_lowdimuon_mass = muJet->vertexMass();
455     m_lowdimuon_pt = muJet->vertexMomentum().perp();
456     m_lowdimuon_eta = muJet->vertexMomentum().eta();
457     m_lowdimuon_phi = muJet->vertexMomentum().phi();
458     m_lowdimuon_dr = muJet->dRmax(true);
459 pivarski 1.3 if (muJet->daughter(0)->charge() > 0) {
460     m_lowdimuon_pluspx = muJet->vertexMomentum(0).x();
461     m_lowdimuon_pluspy = muJet->vertexMomentum(0).y();
462     m_lowdimuon_pluspz = muJet->vertexMomentum(0).z();
463     m_lowdimuon_minuspx = muJet->vertexMomentum(1).x();
464     m_lowdimuon_minuspy = muJet->vertexMomentum(1).y();
465     m_lowdimuon_minuspz = muJet->vertexMomentum(1).z();
466     }
467     else {
468     m_lowdimuon_pluspx = muJet->vertexMomentum(1).x();
469     m_lowdimuon_pluspy = muJet->vertexMomentum(1).y();
470     m_lowdimuon_pluspz = muJet->vertexMomentum(1).z();
471     m_lowdimuon_minuspx = muJet->vertexMomentum(0).x();
472     m_lowdimuon_minuspy = muJet->vertexMomentum(0).y();
473     m_lowdimuon_minuspz = muJet->vertexMomentum(0).z();
474     }
475 pivarski 1.1 m_lowdimuon_vprob = muJet->vertexProb();
476    
477     if (closestPrimaryVertex != primaryVertices->end()) {
478     m_lowdimuon_vx = muJet->vertexPoint().x() - closestPrimaryVertex->x();
479     m_lowdimuon_vy = muJet->vertexPoint().y() - closestPrimaryVertex->y();
480     m_lowdimuon_vz = muJet->vertexPoint().z() - closestPrimaryVertex->z();
481     }
482    
483 pivarski 1.2 GlobalVector complus;
484     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0, true);
485     else complus = muJet->daughterCOMrot(1, true);
486     m_lowdimuon_compluspx = complus.x();
487     m_lowdimuon_compluspy = complus.y();
488     m_lowdimuon_compluspz = complus.z();
489    
490 pivarski 1.1 } // end of replacements with quantities measured at the vertex
491    
492     m_lowdimuon_iso = muJet->centralTrackIsolation();
493    
494     m_lowdimuon->Fill();
495     }
496    
497     ////////////////////////////////////////////////////////// highdimuon
498     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() > 80.) {
499     // ...
500     // m_highdimuon->Fill();
501     }
502    
503     ////////////////////////////////////////////////////////// muplustrack
504     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 3) {
505     // ...
506     // m_muplustrack->Fill();
507     }
508    
509     ////////////////////////////////////////////////////////// quadmuon
510     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 4) {
511     // ...
512     // m_quadmuon->Fill();
513     }
514    
515     ////////////////////////////////////////////////////////// dimudimu
516     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 2) {
517     // ...
518     // m_dimudimu->Fill();
519     }
520    
521     ////////////////////////////////////////////////////////// dimuquadmu
522     else if (muJets->size() == 2 && (
523     ((*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 4) ||
524     ((*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 2)
525     )) {
526     // ...
527     // m_dimuquadmu->Fill();
528     }
529    
530     ////////////////////////////////////////////////////////// quadmuquadmu
531     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 4) {
532     // ...
533     // m_quadmuquadmu->Fill();
534     }
535     }
536    
537    
538     // ------------ method called once each job just before starting event loop ------------
539     void
540     FitNtuple::beginJob()
541     {
542     }
543    
544     // ------------ method called once each job just after ending the event loop ------------
545     void
546     FitNtuple::endJob() {
547     }
548    
549     //define this as a plug-in
550     DEFINE_FWK_MODULE(FitNtuple);