ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/FitNtuple/src/FitNtuple.cc
Revision: 1.1
Committed: Sat Nov 20 03:09:31 2010 UTC (14 years, 5 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-11-19-a
Log Message:
FitNtuple for Aysen and Vadim

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