ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/FitNtuple/src/FitNtuple.cc
Revision: 1.2
Committed: Sat Nov 20 03:40:56 2010 UTC (14 years, 5 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-11-19-b
Changes since 1.1: +22 -5 lines
Log Message:
for safety, put all COM components into the ntuple

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.2 // $Id: FitNtuple.cc,v 1.1 2010/11/20 03:09:31 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    
79     // control sample "lowdimuon" (single mu-jet with two muons)
80     TTree *m_lowdimuon;
81     Float_t m_lowdimuon_genmass;
82     Float_t m_lowdimuon_mass;
83     Float_t m_lowdimuon_pt;
84     Float_t m_lowdimuon_eta;
85     Float_t m_lowdimuon_phi;
86     Float_t m_lowdimuon_dr;
87     Float_t m_lowdimuon_pluspx;
88     Float_t m_lowdimuon_pluspy;
89     Float_t m_lowdimuon_pluspz;
90     Float_t m_lowdimuon_minuspx;
91     Float_t m_lowdimuon_minuspy;
92     Float_t m_lowdimuon_minuspz;
93     Float_t m_lowdimuon_vprob;
94     Float_t m_lowdimuon_vx;
95     Float_t m_lowdimuon_vy;
96     Float_t m_lowdimuon_vz;
97     Float_t m_lowdimuon_iso;
98 pivarski 1.2 Float_t m_lowdimuon_compluspx;
99     Float_t m_lowdimuon_compluspy;
100     Float_t m_lowdimuon_compluspz;
101 pivarski 1.1
102     // signal region "highdimuon" (single mu-jet with two muons)
103     TTree *m_highdimuon;
104    
105     // control sample "muplustrack" (single mu-jet with three muons and the closest track)
106     TTree *m_muplustrack;
107    
108     // signal region "quadmuon" (single mu-jet with four muons)
109     TTree *m_quadmuon;
110    
111     // signal region "dimudimu" (two mu-jets with two muons each)
112     TTree *m_dimudimu;
113    
114     // signal region "dimuquadmu" (one two-muon jet and one four-muon jet)
115     TTree *m_dimuquadmu;
116    
117     // signal region "quadmuquadmu" (two mu-jets with four muons each)
118     TTree *m_quadmuquadmu;
119    
120     };
121    
122     //
123     // constants, enums and typedefs
124     //
125    
126     //
127     // static data member definitions
128     //
129    
130     //
131     // constructors and destructor
132     //
133     FitNtuple::FitNtuple(const edm::ParameterSet& iConfig)
134     : m_src(iConfig.getParameter<edm::InputTag>("src"))
135     , m_getQscale(iConfig.getParameter<bool>("getQscale"))
136     , m_getAlternating(iConfig.getParameter<bool>("getAlternating"))
137     {
138     //now do what ever initialization is needed
139     edm::Service<TFileService> tFile;
140    
141     m_lowdimuon = tFile->make<TTree>("lowdimuon", "lowdimuon");;
142     m_lowdimuon->Branch("qscale", &m_qscale, "qscale/F");
143     m_lowdimuon->Branch("run", &m_run, "run/I");
144     m_lowdimuon->Branch("event", &m_event, "event/I");
145     m_lowdimuon->Branch("trigger", &m_trigger, "trigger/I");
146     m_lowdimuon->Branch("genmass", &m_lowdimuon_genmass, "genmass/F");
147     m_lowdimuon->Branch("mass", &m_lowdimuon_mass, "mass/F");
148     m_lowdimuon->Branch("pt", &m_lowdimuon_pt, "pt/F");
149     m_lowdimuon->Branch("eta", &m_lowdimuon_eta, "eta/F");
150     m_lowdimuon->Branch("phi", &m_lowdimuon_phi, "phi/F");
151     m_lowdimuon->Branch("dr", &m_lowdimuon_dr, "dr/F");
152     m_lowdimuon->Branch("pluspx", &m_lowdimuon_pluspx, "pluspx/F");
153     m_lowdimuon->Branch("pluspy", &m_lowdimuon_pluspy, "pluspy/F");
154     m_lowdimuon->Branch("pluspz", &m_lowdimuon_pluspz, "pluspz/F");
155     m_lowdimuon->Branch("minuspx", &m_lowdimuon_minuspx, "minuspx/F");
156     m_lowdimuon->Branch("minuspy", &m_lowdimuon_minuspy, "minuspy/F");
157     m_lowdimuon->Branch("minuspz", &m_lowdimuon_minuspz, "minuspz/F");
158     m_lowdimuon->Branch("vprob", &m_lowdimuon_vprob, "vprob/F");
159     m_lowdimuon->Branch("vx", &m_lowdimuon_vx, "vx/F");
160     m_lowdimuon->Branch("vy", &m_lowdimuon_vy, "vy/F");
161     m_lowdimuon->Branch("vz", &m_lowdimuon_vz, "vz/F");
162     m_lowdimuon->Branch("iso", &m_lowdimuon_iso, "iso/F");
163 pivarski 1.2 m_lowdimuon->Branch("compluspx", &m_lowdimuon_compluspx, "compluspx/F");
164     m_lowdimuon->Branch("compluspy", &m_lowdimuon_compluspy, "compluspy/F");
165     m_lowdimuon->Branch("compluspz", &m_lowdimuon_compluspz, "compluspz/F");
166 pivarski 1.1
167     m_highdimuon = tFile->make<TTree>("highdimuon", "highdimuon");;
168     m_highdimuon->Branch("qscale", &m_qscale, "qscale/F");
169     m_highdimuon->Branch("run", &m_run, "run/I");
170     m_highdimuon->Branch("event", &m_event, "event/I");
171     m_highdimuon->Branch("trigger", &m_trigger, "trigger/I");
172    
173     m_muplustrack = tFile->make<TTree>("muplustrack", "muplustrack");;
174     m_muplustrack->Branch("qscale", &m_qscale, "qscale/F");
175     m_muplustrack->Branch("run", &m_run, "run/I");
176     m_muplustrack->Branch("event", &m_event, "event/I");
177     m_muplustrack->Branch("trigger", &m_trigger, "trigger/I");
178    
179     m_quadmuon = tFile->make<TTree>("quadmuon", "quadmuon");;
180     m_quadmuon->Branch("qscale", &m_qscale, "qscale/F");
181     m_quadmuon->Branch("run", &m_run, "run/I");
182     m_quadmuon->Branch("event", &m_event, "event/I");
183     m_quadmuon->Branch("trigger", &m_trigger, "trigger/I");
184    
185     m_dimudimu = tFile->make<TTree>("dimudimu", "dimudimu");;
186     m_dimudimu->Branch("qscale", &m_qscale, "qscale/F");
187     m_dimudimu->Branch("run", &m_run, "run/I");
188     m_dimudimu->Branch("event", &m_event, "event/I");
189     m_dimudimu->Branch("trigger", &m_trigger, "trigger/I");
190    
191     m_dimuquadmu = tFile->make<TTree>("dimuquadmu", "dimuquadmu");;
192     m_dimuquadmu->Branch("qscale", &m_qscale, "qscale/F");
193     m_dimuquadmu->Branch("run", &m_run, "run/I");
194     m_dimuquadmu->Branch("event", &m_event, "event/I");
195     m_dimuquadmu->Branch("trigger", &m_trigger, "trigger/I");
196    
197     m_quadmuquadmu = tFile->make<TTree>("quadmuquadmu", "quadmuquadmu");;
198     m_quadmuquadmu->Branch("qscale", &m_qscale, "qscale/F");
199     m_quadmuquadmu->Branch("run", &m_run, "run/I");
200     m_quadmuquadmu->Branch("event", &m_event, "event/I");
201     m_quadmuquadmu->Branch("trigger", &m_trigger, "trigger/I");
202     }
203    
204    
205     FitNtuple::~FitNtuple()
206     {
207    
208     // do anything here that needs to be done at desctruction time
209     // (e.g. close files, deallocate resources etc.)
210    
211     }
212    
213    
214     //
215     // member functions
216     //
217    
218     // ------------ method called to for each event ------------
219     void
220     FitNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
221     {
222     // the QCD scale only exists for certain kinds of Monte Carlo, but if it's available, get it
223     m_qscale = 0.;
224     if (m_getQscale) {
225     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
226     iEvent.getByLabel("generator", genEventInfoProduct);
227     m_qscale = genEventInfoProduct->qScale();
228     }
229    
230     // the pair-gun Monte Carlos have a (never used) feature called
231     // alteration; the real events are the ones in which particleNumber == 0
232     bool alternating = true;
233     if (m_getAlternating) {
234     edm::Handle<unsigned int> particleNumber;
235     iEvent.getByLabel("generator", "particleNumber", particleNumber);
236     alternating = (*particleNumber == 0);
237     }
238     if (!alternating) return;
239    
240     // get the run and event number
241     m_run = iEvent.id().run();
242     m_event = iEvent.id().event();
243    
244     // mu-jets (muons grouped by mass and vertex compatibility)
245     edm::Handle<pat::MultiMuonCollection> muJets;
246     iEvent.getByLabel(m_src, muJets);
247    
248     // // orphans (muons not found in any group)
249     // edm::Handle<pat::MuonCollection> orphans;
250     // iEvent.getByLabel(m_src, "Orphans", orphans);
251    
252     // // all muons
253     // edm::Handle<pat::MuonCollection> muons;
254     // iEvent.getByLabel("cleanPatMuons", muons);
255    
256     // // all tracker-tracks
257     // edm::Handle<reco::TrackCollection> tracks;
258     // iEvent.getByLabel("generalTracks", tracks);
259    
260     // // generator-level 4-vectors
261     // edm::Handle<reco::GenParticleCollection> genParticles;
262     // iEvent.getByLabel("genParticles", genParticles);
263    
264     // find the closest primary vertex (in Z) to the first muJet with a valid vertex
265     edm::Handle<reco::VertexCollection> primaryVertices;
266     iEvent.getByLabel("offlinePrimaryVertices", primaryVertices);
267     reco::VertexCollection::const_iterator closestPrimaryVertex = primaryVertices->end();
268     if (muJets->size() > 0) {
269     pat::MultiMuonCollection::const_iterator muJet0 = muJets->end();
270     for (pat::MultiMuonCollection::const_iterator muJet = muJets->begin(); muJet != muJets->end(); ++muJet) {
271     if (muJet->vertexValid()) {
272     muJet0 = muJet;
273     break;
274     }
275     }
276    
277     if (muJet0 != muJets->end()) {
278     for (reco::VertexCollection::const_iterator vertex = primaryVertices->begin(); vertex != primaryVertices->end(); ++vertex) {
279     if (vertex->isValid() && !vertex->isFake() && vertex->tracksSize() > 3 && fabs(vertex->z()) < 24.) {
280     if (closestPrimaryVertex == primaryVertices->end() || fabs(vertex->z() - muJet0->vertexPoint().z()) < fabs(closestPrimaryVertex->z() - muJet0->vertexPoint().z())) {
281     closestPrimaryVertex = vertex;
282     }
283     } // end vertex quality cuts
284     } // end loop over primary vertices
285     } // end if muJet0 exists
286     } // end if muJets->size > 0
287    
288     // find out which trigger bits were fired
289     edm::Handle<pat::TriggerEvent> triggerEvent;
290     iEvent.getByLabel("patTriggerEvent", triggerEvent);
291     m_trigger = 0;
292     if (triggerEvent->path("HLT_Mu5") && triggerEvent->path("HLT_Mu5")->wasAccept()) m_trigger += 1;
293     if (triggerEvent->path("HLT_Mu9") && triggerEvent->path("HLT_Mu9")->wasAccept()) m_trigger += 2;
294     if (triggerEvent->path("HLT_Mu11") && triggerEvent->path("HLT_Mu11")->wasAccept()) m_trigger += 4;
295    
296     ////////////////////////////////////////////////////////// lowdimuon:
297     if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() < 80.) {
298     pat::MultiMuonCollection::const_iterator muJet = muJets->begin();
299    
300     // generator-level mass using matched genParticles (for resolution of fit peak)
301     m_lowdimuon_genmass = -1000.;
302     if (dynamic_cast<const pat::Muon*>(muJet->daughter(0))->genParticlesSize() == 1 || dynamic_cast<const pat::Muon*>(muJet->daughter(1))->genParticlesSize() == 1) {
303     const reco::GenParticle *mu0 = dynamic_cast<const pat::Muon*>(muJet->daughter(0))->genParticle();
304     const reco::GenParticle *mu1 = dynamic_cast<const pat::Muon*>(muJet->daughter(1))->genParticle();
305    
306     double total_energy = mu0->energy() + mu1->energy();
307     double total_px = mu0->px() + mu1->px();
308     double total_py = mu0->py() + mu1->py();
309     double total_pz = mu0->pz() + mu1->pz();
310     m_lowdimuon_genmass = sqrt(total_energy*total_energy - total_px*total_px - total_py*total_py - total_pz*total_pz);
311     }
312    
313     m_lowdimuon_mass = muJet->mass();
314     m_lowdimuon_pt = muJet->pt();
315     m_lowdimuon_eta = muJet->eta();
316     m_lowdimuon_phi = muJet->phi();
317     m_lowdimuon_dr = muJet->dRmax();
318     m_lowdimuon_pluspx = muJet->daughter(0)->px();
319     m_lowdimuon_pluspy = muJet->daughter(0)->py();
320     m_lowdimuon_pluspz = muJet->daughter(0)->pz();
321     m_lowdimuon_minuspx = muJet->daughter(1)->px();
322     m_lowdimuon_minuspy = muJet->daughter(1)->py();
323     m_lowdimuon_minuspz = muJet->daughter(1)->pz();
324     m_lowdimuon_vprob = -1000.;
325     m_lowdimuon_vx = -1000.;
326     m_lowdimuon_vy = -1000.;
327     m_lowdimuon_vz = -1000.;
328    
329 pivarski 1.2 GlobalVector complus;
330     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0);
331     else complus = muJet->daughterCOMrot(1);
332     m_lowdimuon_compluspx = complus.x();
333     m_lowdimuon_compluspy = complus.y();
334     m_lowdimuon_compluspz = complus.z();
335    
336     // replace all values with vertex-updated values if vertex fitting succeeded
337 pivarski 1.1 if (muJet->vertexValid()) {
338     m_lowdimuon_mass = muJet->vertexMass();
339     m_lowdimuon_pt = muJet->vertexMomentum().perp();
340     m_lowdimuon_eta = muJet->vertexMomentum().eta();
341     m_lowdimuon_phi = muJet->vertexMomentum().phi();
342     m_lowdimuon_dr = muJet->dRmax(true);
343     m_lowdimuon_pluspx = muJet->vertexMomentum(0).x();
344     m_lowdimuon_pluspy = muJet->vertexMomentum(0).y();
345     m_lowdimuon_pluspz = muJet->vertexMomentum(0).z();
346     m_lowdimuon_minuspx = muJet->vertexMomentum(1).x();
347     m_lowdimuon_minuspy = muJet->vertexMomentum(1).y();
348     m_lowdimuon_minuspz = muJet->vertexMomentum(1).z();
349     m_lowdimuon_vprob = muJet->vertexProb();
350    
351     if (closestPrimaryVertex != primaryVertices->end()) {
352     m_lowdimuon_vx = muJet->vertexPoint().x() - closestPrimaryVertex->x();
353     m_lowdimuon_vy = muJet->vertexPoint().y() - closestPrimaryVertex->y();
354     m_lowdimuon_vz = muJet->vertexPoint().z() - closestPrimaryVertex->z();
355     }
356    
357 pivarski 1.2 GlobalVector complus;
358     if (muJet->daughter(0)->charge() > 0) complus = muJet->daughterCOMrot(0, true);
359     else complus = muJet->daughterCOMrot(1, true);
360     m_lowdimuon_compluspx = complus.x();
361     m_lowdimuon_compluspy = complus.y();
362     m_lowdimuon_compluspz = complus.z();
363    
364 pivarski 1.1 } // end of replacements with quantities measured at the vertex
365    
366     m_lowdimuon_iso = muJet->centralTrackIsolation();
367    
368     m_lowdimuon->Fill();
369     }
370    
371     ////////////////////////////////////////////////////////// highdimuon
372     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[0].pt() > 80.) {
373     // ...
374     // m_highdimuon->Fill();
375     }
376    
377     ////////////////////////////////////////////////////////// muplustrack
378     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 3) {
379     // ...
380     // m_muplustrack->Fill();
381     }
382    
383     ////////////////////////////////////////////////////////// quadmuon
384     else if (muJets->size() == 1 && (*muJets)[0].numberOfDaughters() == 4) {
385     // ...
386     // m_quadmuon->Fill();
387     }
388    
389     ////////////////////////////////////////////////////////// dimudimu
390     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 2) {
391     // ...
392     // m_dimudimu->Fill();
393     }
394    
395     ////////////////////////////////////////////////////////// dimuquadmu
396     else if (muJets->size() == 2 && (
397     ((*muJets)[0].numberOfDaughters() == 2 && (*muJets)[1].numberOfDaughters() == 4) ||
398     ((*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 2)
399     )) {
400     // ...
401     // m_dimuquadmu->Fill();
402     }
403    
404     ////////////////////////////////////////////////////////// quadmuquadmu
405     else if (muJets->size() == 2 && (*muJets)[0].numberOfDaughters() == 4 && (*muJets)[1].numberOfDaughters() == 4) {
406     // ...
407     // m_quadmuquadmu->Fill();
408     }
409     }
410    
411    
412     // ------------ method called once each job just before starting event loop ------------
413     void
414     FitNtuple::beginJob()
415     {
416     }
417    
418     // ------------ method called once each job just after ending the event loop ------------
419     void
420     FitNtuple::endJob() {
421     }
422    
423     //define this as a plug-in
424     DEFINE_FWK_MODULE(FitNtuple);