ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MuJetAnalysis/AnalysisTools/IsolationStudyNtuple/src/IsolationStudyNtuple.cc
Revision: 1.3
Committed: Wed Nov 10 02:42:45 2010 UTC (14 years, 5 months ago) by pivarski
Content type: text/plain
Branch: MAIN
CVS Tags: JP-2010-11-09-a
Changes since 1.2: +42 -2 lines
Log Message:
begin study of isolation

File Contents

# User Rev Content
1 pivarski 1.1 // -*- C++ -*-
2     //
3     // Package: IsolationStudyNtuple
4     // Class: IsolationStudyNtuple
5     //
6     /**\class IsolationStudyNtuple IsolationStudyNtuple.cc AnalysisTools/IsolationStudyNtuple/src/IsolationStudyNtuple.cc
7    
8     Description: [one line class summary]
9    
10     Implementation:
11     [Notes on implementation]
12     */
13     //
14     // Original Author: James Pivarski
15     // Created: Tue Nov 9 13:34:12 CST 2010
16 pivarski 1.3 // $Id: IsolationStudyNtuple.cc,v 1.2 2010/11/10 01:16:03 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 pivarski 1.2 #include "FWCore/ParameterSet/interface/ParameterSet.h"
30    
31     #include "AnalysisDataFormats/MuJetAnalysis/interface/MultiMuon.h"
32     #include "CommonTools/UtilAlgos/interface/TFileService.h"
33     #include "DataFormats/MuonReco/interface/Muon.h"
34     #include "DataFormats/PatCandidates/interface/Muon.h"
35     #include "DataFormats/TrackReco/interface/TrackFwd.h"
36     #include "DataFormats/PatCandidates/interface/TriggerEvent.h"
37     #include "DataFormats/PatCandidates/interface/TriggerPath.h"
38     #include "FWCore/ServiceRegistry/interface/Service.h"
39 pivarski 1.3 #include "DataFormats/Candidate/interface/Candidate.h"
40 pivarski 1.2 #include "TTree.h"
41 pivarski 1.1
42     //
43     // class declaration
44     //
45    
46     class IsolationStudyNtuple : public edm::EDAnalyzer {
47     public:
48     explicit IsolationStudyNtuple(const edm::ParameterSet&);
49     ~IsolationStudyNtuple();
50    
51    
52     private:
53     virtual void beginJob() ;
54     virtual void analyze(const edm::Event&, const edm::EventSetup&);
55     virtual void endJob() ;
56    
57     // ----------member data ---------------------------
58 pivarski 1.2 TTree *m_iso01th00, *m_iso03th00, *m_iso05th00, *m_iso07th00, *m_iso10th00;
59     TTree *m_iso01th01, *m_iso03th01, *m_iso05th01, *m_iso07th01, *m_iso10th01;
60     TTree *m_iso01th02, *m_iso03th02, *m_iso05th02, *m_iso07th02, *m_iso10th02;
61     TTree *m_iso01th05, *m_iso03th05, *m_iso05th05, *m_iso07th05, *m_iso10th05;
62     TTree *m_iso01th10, *m_iso03th10, *m_iso05th10, *m_iso07th10, *m_iso10th10;
63     TTree *m_iso01th20, *m_iso03th20, *m_iso05th20, *m_iso07th20, *m_iso10th20;
64 pivarski 1.3
65 pivarski 1.2 Int_t m_daughters;
66     Float_t m_mass;
67     Float_t m_dRmax;
68     Float_t m_dphimax;
69     Float_t m_detamax;
70     Float_t m_vertexProb;
71     Float_t m_pt;
72     Float_t m_phi;
73     Float_t m_eta;
74     Float_t m_pt1;
75     Float_t m_pt2;
76     Float_t m_trackiso;
77     Int_t m_numAboveThresh;
78 pivarski 1.3 Float_t m_closestOtherdR;
79     Float_t m_closestOtherdphi;
80     Float_t m_closestOtherdeta;
81     Float_t m_other_trackiso;
82     Int_t m_other_numAboveThresh;
83 pivarski 1.1 };
84    
85     //
86     // constants, enums and typedefs
87     //
88    
89     //
90     // static data member definitions
91     //
92    
93     //
94     // constructors and destructor
95     //
96     IsolationStudyNtuple::IsolationStudyNtuple(const edm::ParameterSet& iConfig)
97     {
98     //now do what ever initialization is needed
99 pivarski 1.2 edm::Service<TFileService> tFile;
100     m_iso01th00 = tFile->make<TTree>("iso01th00", "iso01th00"); m_iso03th00 = tFile->make<TTree>("iso03th00", "iso03th00"); m_iso05th00 = tFile->make<TTree>("iso05th00", "iso05th00"); m_iso07th00 = tFile->make<TTree>("iso07th00", "iso07th00"); m_iso10th00 = tFile->make<TTree>("iso10th00", "iso10th00");
101     m_iso01th01 = tFile->make<TTree>("iso01th01", "iso01th01"); m_iso03th01 = tFile->make<TTree>("iso03th01", "iso03th01"); m_iso05th01 = tFile->make<TTree>("iso05th01", "iso05th01"); m_iso07th01 = tFile->make<TTree>("iso07th01", "iso07th01"); m_iso10th01 = tFile->make<TTree>("iso10th01", "iso10th01");
102     m_iso01th02 = tFile->make<TTree>("iso01th02", "iso01th02"); m_iso03th02 = tFile->make<TTree>("iso03th02", "iso03th02"); m_iso05th02 = tFile->make<TTree>("iso05th02", "iso05th02"); m_iso07th02 = tFile->make<TTree>("iso07th02", "iso07th02"); m_iso10th02 = tFile->make<TTree>("iso10th02", "iso10th02");
103     m_iso01th05 = tFile->make<TTree>("iso01th05", "iso01th05"); m_iso03th05 = tFile->make<TTree>("iso03th05", "iso03th05"); m_iso05th05 = tFile->make<TTree>("iso05th05", "iso05th05"); m_iso07th05 = tFile->make<TTree>("iso07th05", "iso07th05"); m_iso10th05 = tFile->make<TTree>("iso10th05", "iso10th05");
104     m_iso01th10 = tFile->make<TTree>("iso01th10", "iso01th10"); m_iso03th10 = tFile->make<TTree>("iso03th10", "iso03th10"); m_iso05th10 = tFile->make<TTree>("iso05th10", "iso05th10"); m_iso07th10 = tFile->make<TTree>("iso07th10", "iso07th10"); m_iso10th10 = tFile->make<TTree>("iso10th10", "iso10th10");
105     m_iso01th20 = tFile->make<TTree>("iso01th20", "iso01th20"); m_iso03th20 = tFile->make<TTree>("iso03th20", "iso03th20"); m_iso05th20 = tFile->make<TTree>("iso05th20", "iso05th20"); m_iso07th20 = tFile->make<TTree>("iso07th20", "iso07th20"); m_iso10th20 = tFile->make<TTree>("iso10th20", "iso10th20");
106    
107     for (int cone = 1; cone <= 5; cone++) {
108     for (int thresh = 0; thresh <= 5; thresh++) {
109     TTree *ttree = NULL;
110     if (cone == 1) {
111     if (thresh == 0) ttree = m_iso01th00;
112     else if (thresh == 1) ttree = m_iso01th01;
113     else if (thresh == 2) ttree = m_iso01th02;
114     else if (thresh == 3) ttree = m_iso01th05;
115     else if (thresh == 4) ttree = m_iso01th10;
116     else if (thresh == 5) ttree = m_iso01th20;
117     else assert(false);
118     }
119     else if (cone == 2) {
120     if (thresh == 0) ttree = m_iso03th00;
121     else if (thresh == 1) ttree = m_iso03th01;
122     else if (thresh == 2) ttree = m_iso03th02;
123     else if (thresh == 3) ttree = m_iso03th05;
124     else if (thresh == 4) ttree = m_iso03th10;
125     else if (thresh == 5) ttree = m_iso03th20;
126     else assert(false);
127     }
128     else if (cone == 3) {
129     if (thresh == 0) ttree = m_iso05th00;
130     else if (thresh == 1) ttree = m_iso05th01;
131     else if (thresh == 2) ttree = m_iso05th02;
132     else if (thresh == 3) ttree = m_iso05th05;
133     else if (thresh == 4) ttree = m_iso05th10;
134     else if (thresh == 5) ttree = m_iso05th20;
135     else assert(false);
136     }
137     else if (cone == 4) {
138     if (thresh == 0) ttree = m_iso07th00;
139     else if (thresh == 1) ttree = m_iso07th01;
140     else if (thresh == 2) ttree = m_iso07th02;
141     else if (thresh == 3) ttree = m_iso07th05;
142     else if (thresh == 4) ttree = m_iso07th10;
143     else if (thresh == 5) ttree = m_iso07th20;
144     else assert(false);
145     }
146     else if (cone == 5) {
147     if (thresh == 0) ttree = m_iso10th00;
148     else if (thresh == 1) ttree = m_iso10th01;
149     else if (thresh == 2) ttree = m_iso10th02;
150     else if (thresh == 3) ttree = m_iso10th05;
151     else if (thresh == 4) ttree = m_iso10th10;
152     else if (thresh == 5) ttree = m_iso10th20;
153     else assert(false);
154     }
155     else assert(false);
156    
157     ttree->Branch("daughters", &m_daughters, "daughters/I");
158     ttree->Branch("mass", &m_mass, "mass/F");
159     ttree->Branch("dRmax", &m_dRmax, "dRmax/F");
160     ttree->Branch("dphimax", &m_dphimax, "dphimax/F");
161     ttree->Branch("detamax", &m_detamax, "detamax/F");
162     ttree->Branch("vertexProb", &m_vertexProb, "vertexProb/F");
163     ttree->Branch("pt", &m_pt, "pt/F");
164     ttree->Branch("phi", &m_phi, "phi/F");
165     ttree->Branch("eta", &m_eta, "eta/F");
166     ttree->Branch("pt1", &m_pt1, "pt1/F");
167     ttree->Branch("pt2", &m_pt2, "pt2/F");
168     ttree->Branch("trackiso", &m_trackiso, "trackiso/F");
169     ttree->Branch("numAboveThresh", &m_numAboveThresh, "numAboveThresh/I");
170 pivarski 1.3 ttree->Branch("closestOtherdR", &m_closestOtherdR, "closestOtherdR/F");
171     ttree->Branch("closestOtherdphi", &m_closestOtherdphi, "closestOtherdphi/F");
172     ttree->Branch("closestOtherdeta", &m_closestOtherdeta, "closestOtherdeta/F");
173     ttree->Branch("other_trackiso", &m_other_trackiso, "other_trackiso/F");
174     ttree->Branch("other_numAboveThresh", &m_other_numAboveThresh, "other_numAboveThresh/I");
175    
176 pivarski 1.2 }
177     }
178 pivarski 1.1 }
179    
180    
181     IsolationStudyNtuple::~IsolationStudyNtuple()
182     {
183    
184     // do anything here that needs to be done at desctruction time
185     // (e.g. close files, deallocate resources etc.)
186    
187     }
188    
189    
190     //
191     // member functions
192     //
193    
194     // ------------ method called to for each event ------------
195     void
196     IsolationStudyNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
197     {
198 pivarski 1.2 edm::Handle<pat::MultiMuonCollection> muJets;
199     iEvent.getByLabel("MuJetProducer", muJets);
200 pivarski 1.1
201 pivarski 1.2 edm::Handle<reco::TrackCollection> tracks;
202     iEvent.getByLabel("generalTracks", tracks);
203 pivarski 1.1
204 pivarski 1.2 for (int cone = 1; cone <= 5; cone++) {
205     for (int thresh = 0; thresh <= 5; thresh++) {
206     double centralCone;
207     if (cone == 1) centralCone = 0.1;
208     else if (cone == 2) centralCone = 0.3;
209     else if (cone == 3) centralCone = 0.5;
210     else if (cone == 4) centralCone = 0.7;
211     else if (cone == 5) centralCone = 1.0;
212     else assert(false);
213    
214     double centralThreshold;
215     if (thresh == 0) centralThreshold = 0.0;
216     else if (thresh == 1) centralThreshold = 0.1;
217     else if (thresh == 2) centralThreshold = 0.2;
218     else if (thresh == 3) centralThreshold = 0.5;
219     else if (thresh == 4) centralThreshold = 1.0;
220     else if (thresh == 5) centralThreshold = 2.0;
221     else assert(false);
222    
223     TTree *ttree = NULL;
224     if (cone == 1) {
225     if (thresh == 0) ttree = m_iso01th00;
226     else if (thresh == 1) ttree = m_iso01th01;
227     else if (thresh == 2) ttree = m_iso01th02;
228     else if (thresh == 3) ttree = m_iso01th05;
229     else if (thresh == 4) ttree = m_iso01th10;
230     else if (thresh == 5) ttree = m_iso01th20;
231     else assert(false);
232     }
233     else if (cone == 2) {
234     if (thresh == 0) ttree = m_iso03th00;
235     else if (thresh == 1) ttree = m_iso03th01;
236     else if (thresh == 2) ttree = m_iso03th02;
237     else if (thresh == 3) ttree = m_iso03th05;
238     else if (thresh == 4) ttree = m_iso03th10;
239     else if (thresh == 5) ttree = m_iso03th20;
240     else assert(false);
241     }
242     else if (cone == 3) {
243     if (thresh == 0) ttree = m_iso05th00;
244     else if (thresh == 1) ttree = m_iso05th01;
245     else if (thresh == 2) ttree = m_iso05th02;
246     else if (thresh == 3) ttree = m_iso05th05;
247     else if (thresh == 4) ttree = m_iso05th10;
248     else if (thresh == 5) ttree = m_iso05th20;
249     else assert(false);
250     }
251     else if (cone == 4) {
252     if (thresh == 0) ttree = m_iso07th00;
253     else if (thresh == 1) ttree = m_iso07th01;
254     else if (thresh == 2) ttree = m_iso07th02;
255     else if (thresh == 3) ttree = m_iso07th05;
256     else if (thresh == 4) ttree = m_iso07th10;
257     else if (thresh == 5) ttree = m_iso07th20;
258     else assert(false);
259     }
260     else if (cone == 5) {
261     if (thresh == 0) ttree = m_iso10th00;
262     else if (thresh == 1) ttree = m_iso10th01;
263     else if (thresh == 2) ttree = m_iso10th02;
264     else if (thresh == 3) ttree = m_iso10th05;
265     else if (thresh == 4) ttree = m_iso10th10;
266     else if (thresh == 5) ttree = m_iso10th20;
267     else assert(false);
268     }
269     else assert(false);
270    
271     for (pat::MultiMuonCollection::const_iterator muJetOriginal = muJets->begin(); muJetOriginal != muJets->end(); ++muJetOriginal) {
272     pat::MultiMuon muJet(*muJetOriginal);
273     muJet.calculateTrackIsolation(&*tracks, centralCone, 0., centralThreshold, 0.);
274     muJet.calculateNumberAboveThresholdIsolation(&*tracks, centralCone, 0., centralThreshold, 0.);
275    
276     m_daughters = muJet.numberOfDaughters();
277     m_mass = muJet.mass();
278     m_dRmax = muJet.dRmax();
279     m_dphimax = 0.;
280     m_detamax = 0.;
281     m_pt1 = 0.;
282     m_pt2 = 0.;
283     for (int i = 0; i < m_daughters; i++) {
284     for (int j = i+1; j < m_daughters; j++) {
285     if (fabs(muJet.dphi(i, j)) > fabs(m_dphimax)) m_dphimax = muJet.dphi(i, j);
286     if (fabs(muJet.deta(i, j)) > fabs(m_detamax)) m_detamax = muJet.deta(i, j);
287     }
288    
289     if (muJet.daughter(i)->pt() > m_pt1) {
290     m_pt2 = m_pt1;
291     m_pt1 = muJet.daughter(i)->pt();
292     }
293     else if (muJet.daughter(i)->pt() > m_pt2) {
294     m_pt2 = muJet.daughter(i)->pt();
295     }
296     }
297     if (muJet.vertexValid()) m_vertexProb = muJet.vertexProb();
298     else m_vertexProb = -1.;
299     m_pt = muJet.pt();
300     m_phi = muJet.phi();
301     m_eta = muJet.eta();
302     m_trackiso = muJet.centralTrackIsolation();
303     m_numAboveThresh = muJet.centralNumberAboveThreshold();
304    
305 pivarski 1.3 m_closestOtherdR = -1.;
306     m_closestOtherdphi = 0.;
307     m_closestOtherdeta = 0.;
308     m_other_trackiso = 0.;
309     m_other_numAboveThresh = 0;
310     pat::MultiMuonCollection::const_iterator closestOther = muJets->end();
311     for (pat::MultiMuonCollection::const_iterator other = muJetOriginal; other != muJets->end(); ++other) {
312     if (other != muJetOriginal) {
313     double dphi = muJetOriginal->phi() - other->phi();
314     while (dphi > M_PI) dphi -= 2.*M_PI;
315     while (dphi < -M_PI) dphi += 2.*M_PI;
316     double deta = muJetOriginal->eta() - other->eta();
317     double dR = sqrt(dphi*dphi + deta*deta);
318    
319     if (closestOther == muJets->end() || m_closestOtherdR > dR) {
320     closestOther = other;
321     m_closestOtherdR = dR;
322     m_closestOtherdphi = dphi;
323     m_closestOtherdeta = deta;
324     pat::MultiMuon mJ(*other);
325     mJ.calculateTrackIsolation(&*tracks, centralCone, 0., centralThreshold, 0.);
326     mJ.calculateNumberAboveThresholdIsolation(&*tracks, centralCone, 0., centralThreshold, 0.);
327     m_other_trackiso = mJ.centralTrackIsolation();
328     m_other_numAboveThresh = mJ.centralNumberAboveThreshold();
329     }
330     }
331     }
332    
333 pivarski 1.2 ttree->Fill();
334     }
335     }
336     }
337 pivarski 1.1 }
338    
339    
340     // ------------ method called once each job just before starting event loop ------------
341     void
342     IsolationStudyNtuple::beginJob()
343     {
344     }
345    
346     // ------------ method called once each job just after ending the event loop ------------
347     void
348     IsolationStudyNtuple::endJob() {
349     }
350    
351     //define this as a plug-in
352     DEFINE_FWK_MODULE(IsolationStudyNtuple);