ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/SelMods/src/HwwExampleAnalysisMod.cc
Revision: 1.5
Committed: Wed Oct 20 02:44:02 2010 UTC (14 years, 6 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_015b, Mit_015a
Changes since 1.4: +2 -3 lines
Log Message:
fixing vertex issues

File Contents

# Content
1 // $Id $
2
3 #include "MitPhysics/SelMods/interface/HwwExampleAnalysisMod.h"
4 #include <TH1D.h>
5 #include <TH2D.h>
6 #include <TParameter.h>
7 #include "MitAna/DataUtil/interface/Debug.h"
8 #include "MitAna/DataTree/interface/Names.h"
9 #include "MitPhysics/Init/interface/ModNames.h"
10 #include "MitAna/DataCont/interface/ObjArray.h"
11 #include "MitCommon/MathTools/interface/MathUtils.h"
12 #include "MitAna/DataTree/interface/ParticleCol.h"
13 #include "TFile.h"
14 #include "TTree.h"
15
16 using namespace mithep;
17 ClassImp(mithep::HwwExampleAnalysisMod)
18
19 //--------------------------------------------------------------------------------------------------
20 HwwExampleAnalysisMod::HwwExampleAnalysisMod(const char *name, const char *title) :
21 BaseMod(name,title),
22 fMuonBranchName(Names::gkMuonBrn),
23 fMetName("NoDefaultNameSet"),
24 fCleanJetsName("NoDefaultNameSet"),
25 fVertexName(ModNames::gkGoodVertexesName),
26 fMuons(0),
27 fMet(0),
28 fVertices(0)
29 {
30 // Constructor.
31 }
32
33 //--------------------------------------------------------------------------------------------------
34 void HwwExampleAnalysisMod::Begin()
35 {
36 // Run startup code on the client machine. For this module, we dont do
37 // anything here.
38 }
39
40 //--------------------------------------------------------------------------------------------------
41 void HwwExampleAnalysisMod::SlaveBegin()
42 {
43 // Run startup code on the computer (slave) doing the actual analysis. Here,
44 // we typically initialize histograms and other analysis objects and request
45 // branches. For this module, we request a branch of the MitTree.
46
47 // Load Branches
48 ReqBranch(fMuonBranchName, fMuons);
49
50 //Create your histograms here
51
52 //*************************************************************************************************
53 // Selection Histograms
54 //*************************************************************************************************
55 AddTH1(fHWWSelection,"hHWWSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
56 AddTH1(fHWWToEESelection,"hHWWToEESelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
57 AddTH1(fHWWToMuMuSelection,"hHWWToMuMuSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
58 AddTH1(fHWWToEMuSelection,"hHWWToEMuSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
59
60 //***********************************************************************************************
61 // Histograms after preselection
62 //***********************************************************************************************
63 AddTH1(fLeptonEta ,"hLeptonEta",";LeptonEta;Number of Events",100,-5.,5.0);
64 AddTH1(fLeptonPtMax ,"hLeptonPtMax",";Lepton P_t Max;Number of Events",150,0.,150.);
65 AddTH1(fLeptonPtMin ,"hLeptonPtMin",";Lepton P_t Min;Number of Events",150,0.,150.);
66 AddTH1(fMetPtHist ,"hMetPtHist",";Met;Number of Events",150,0.,300.);
67 AddTH1(fMetPhiHist ,"hMetPhiHist",";#phi;Number of Events",28,-3.5,3.5);
68 AddTH1(fUncorrMetPtHist ,"hUncorrMetPtHist",";Met;Number of Events",150,0.,300.);
69 AddTH1(fUncorrMetPhiHist ,"hUncorrMetPhiHist",";#phi;Number of Events",28,-3.5,3.5);
70 AddTH1(fDeltaPhiLeptons ,"hDeltaPhiLeptons",";#Delta#phi_{ll};Number of Events",90,0,180);
71 AddTH1(fDeltaEtaLeptons ,"hDeltaEtaLeptons",";#Delta#eta_{ll};Number of Events",100,-50.,5.0);
72 AddTH1(fDileptonMass ,"hDileptonMass",";Mass_{ll};Number of Events",150,0.,300.);
73
74 //***********************************************************************************************
75 // N-1 Histograms
76 //***********************************************************************************************
77 //All events
78 AddTH1(fLeptonPtMax_NMinusOne ,"hLeptonPtMax_NMinusOne",
79 ";Lepton P_t Max;Number of Events",150,0.,150.);
80 AddTH1(fLeptonPtMin_NMinusOne ,"hLeptonPtMin_NMinusOne",
81 ";Lepton P_t Min;Number of Events",150,0.,150.);
82 AddTH1(fMetPtHist_NMinusOne ,"hMetPtHist_NMinusOne",
83 ";Met;Number of Events",150,0.,300.);
84 AddTH1(fMetPhiHist_NMinusOne ,"hMetPhiHist_NMinusOne",
85 ";#phi;Number of Events",28,-3.5,3.5);
86 AddTH1(fMETdeltaPhilEtHist_NMinusOne ,"hMETdeltaPhilEtHist_NMinusOne",
87 ";METdeltaPhilEtHist;Number of Events",150,0.,300.);
88 AddTH1(fNCentralJets_NMinusOne ,"hNCentralJets_NMinusOne",
89 ";Number of Central Jets;Number of Events",6,-0.5,5.5);
90 AddTH1(fNSoftMuonsHist_NMinusOne ,"hNSoftMuonsHist_NMinusOne",
91 ";Number of Dirty Muons; Number of Events",6,-0.5,5.5);
92 AddTH1(fDeltaPhiLeptons_NMinusOne ,"hDeltaPhiLeptons_NMinusOne",
93 ";#Delta#phi_{ll};Number of Events",90,0,180);
94 AddTH1(fDeltaEtaLeptons_NMinusOne ,"hDeltaEtaLeptons_NMinusOne",
95 ";#Delta#eta_{ll};Number of Events",100,-5.0,5.0);
96 AddTH1(fDileptonMass_NMinusOne ,"hDileptonMass_NMinusOne",
97 ";Mass_{ll};Number of Events",150,0.,300.);
98 AddTH1(fMinDeltaPhiLeptonMet_NMinusOne ,"hMinDeltaPhiLeptonMet_NMinusOne",
99 ";Min #Delta#phi_{l,Met};Number of Events",90,0.,180);
100
101 //***********************************************************************************************
102 // After all cuts Histograms
103 //***********************************************************************************************
104 AddTH1(fMinDeltaPhiLeptonMet_afterCuts ,"hMinDeltaPhiLeptonMet_afterCuts",
105 ";Min #Delta#phi_{l,Met};Number of Events",90,0.,180);
106 AddTH1(fMtLepton1_afterCuts ,"hMtLepton1_afterCuts",
107 ";M_t (Lepton1,Met);Number of Events",100,0.,200.);
108 AddTH1(fMtLepton2_afterCuts ,"hMtLepton2_afterCuts",
109 ";M_t (Lepton2,Met);Number of Events",100,0.,200.);
110 AddTH1(fMtHiggs_afterCuts ,"hMtHiggs_afterCuts",
111 ";M_t (l1+l2+Met);Number of Events",150,0.,300.);
112 AddTH1(fLeptonPtPlusMet_afterCuts ,"hLeptonPtPlusMet_afterCuts",
113 ";LeptonPtPlusMet;Number of Events",150,0., 300.);
114
115 }
116
117 //--------------------------------------------------------------------------------------------------
118 void HwwExampleAnalysisMod::Process()
119 {
120 // Process entries of the tree. For this module, we just load the branches and
121 LoadBranch(fMuonBranchName);
122
123 //Obtain all the good objects from the event cleaning module
124 fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
125 ObjArray<Muon> *CleanMuons = dynamic_cast<ObjArray<Muon>* >(FindObjThisEvt(ModNames::gkCleanMuonsName));
126 ParticleOArr *CleanLeptons = dynamic_cast<mithep::ParticleOArr*>
127 (FindObjThisEvt(ModNames::gkMergedLeptonsName));
128 ObjArray<Jet> *CleanJets = dynamic_cast<ObjArray<Jet>* >
129 (FindObjThisEvt(fCleanJetsName.Data()));
130 TParameter<Double_t> *NNLOWeight = GetObjThisEvt<TParameter<Double_t> >("NNLOWeight");
131
132 MetCol *met = dynamic_cast<ObjArray<Met>* >(FindObjThisEvt(fMetName));
133 const Met *caloMet = 0;
134 if (met) {
135 caloMet = met->At(0);
136 } else {
137 cout << "Error: Met Collection " << fMetName << " could not be loaded.\n";
138 return;
139 }
140
141 //***********************************************************************************************
142 //Kinematic PreSelection
143 //***********************************************************************************************
144 // At least two leptons in the event
145 if (CleanLeptons->GetEntries() < 2) return;
146 // Pt1 > 20 && Pt2 > 10
147 if(CleanLeptons->At(0)->Pt() <= 20 || CleanLeptons->At(1)->Pt() <= 10) return;
148 // opposite charge leptons
149 if(CleanLeptons->At(0)->Charge() * CleanLeptons->At(1)->Charge() > 0) return;
150
151 CompositeParticle *dilepton = new CompositeParticle();
152 dilepton->AddDaughter(CleanLeptons->At(0));
153 dilepton->AddDaughter(CleanLeptons->At(1));
154
155 //***********************************************************************************************
156 //Get Dirty Muons: Non-isolated Muons (exclude the clean muons)
157 //***********************************************************************************************
158 ObjArray<Muon> *SoftMuons = new ObjArray<Muon>;
159 for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
160 const Muon *mu = fMuons->At(i);
161 if(!MuonTools::PassSoftMuonCut(mu, fVertices)) continue;
162
163 bool isCleanMuon = kFALSE;
164 for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
165 if(fMuons->At(i) == CleanMuons->At(j) &&
166 CleanMuons->At(j)->Pt() > 10) isCleanMuon = kTRUE;
167 }
168 if(isCleanMuon == kFALSE) SoftMuons->Add(mu);
169 }
170
171 //***********************************************************************************************
172 //Define Event Variables
173 //***********************************************************************************************
174 //delta phi between the 2 leptons in degrees
175 double deltaPhiLeptons = MathUtils::DeltaPhi(CleanLeptons->At(0)->Phi(),
176 CleanLeptons->At(1)->Phi())* 180.0 / TMath::Pi();
177
178 double deltaEtaLeptons = abs(CleanLeptons->At(0)->Eta() - CleanLeptons->At(1)->Eta()) * 180.0 / TMath::Pi();
179
180 double deltaPhiDileptonMet = MathUtils::DeltaPhi(caloMet->Phi(),
181 dilepton->Phi())*180.0 / TMath::Pi();
182
183 double mtHiggs = TMath::Sqrt(2.0*dilepton->Pt() * caloMet->Pt()*
184 (1.0 - cos(deltaPhiDileptonMet * TMath::Pi() / 180.0)));
185
186 //angle between MET and closest lepton
187 double deltaPhiMetLepton[2] = {MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(0)->Phi()),
188 MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(1)->Phi())};
189
190 double mTW[2] = {TMath::Sqrt(2.0*CleanLeptons->At(0)->Pt()*caloMet->Pt()*
191 (1.0 - cos(deltaPhiMetLepton[0]))),
192 TMath::Sqrt(2.0*CleanLeptons->At(1)->Pt()*caloMet->Pt()*
193 (1.0 - cos(deltaPhiMetLepton[1])))};
194
195 double minDeltaPhiMetLepton = (deltaPhiMetLepton[0] < deltaPhiMetLepton[1])?
196 deltaPhiMetLepton[0]:deltaPhiMetLepton[1];
197
198 double METdeltaPhilEt = caloMet->Pt();
199 if(minDeltaPhiMetLepton < TMath::Pi()/2.)
200 METdeltaPhilEt = METdeltaPhilEt * sin(minDeltaPhiMetLepton);
201
202 //count the number of central Jets for vetoing
203 int nCentralJets = 0;
204 for(UInt_t i=0; i<CleanJets->GetEntries(); i++){
205 if(TMath::Abs(CleanJets->At(i)->Eta()) < 5.0 &&
206 CleanJets->At(i)->Pt() > 25.0){
207 nCentralJets++;
208 }
209 }
210
211 //Lepton Type
212 int finalstateType = -1;
213 if (CleanLeptons->At(0)->ObjType() == kMuon && CleanLeptons->At(1)->ObjType() == kMuon ){ // mumu
214 finalstateType = 10;
215 } else if(CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kElectron ){ // ee
216 finalstateType = 11;
217 } else if((CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kMuon) ||
218 (CleanLeptons->At(1)->ObjType() == kElectron && CleanLeptons->At(0)->ObjType() == kMuon)) {
219 finalstateType = 12;
220 } else {
221 cerr << "Error: finalstate lepton type not supported\n";
222 }
223
224 //*********************************************************************************************
225 //Define Cuts
226 //*********************************************************************************************
227 const int nCuts = 7;
228 bool passCut[nCuts] = {false, false, false, false, false, false, false};
229
230 if(CleanLeptons->At(0)->Pt() > 20.0 &&
231 CleanLeptons->At(1)->Pt() >= 20.0) passCut[0] = true;
232
233 if(caloMet->Pt() > 20.0) passCut[1] = true;
234
235 if(dilepton->Mass() > 12.0) passCut[2] = true;
236
237 if(nCentralJets < 1) passCut[5] = true;
238
239 if(CleanLeptons->GetEntries() == 2 &&
240 SoftMuons->GetEntries() == 0) passCut[6] = true;
241
242 if (finalstateType == 10 || finalstateType == 11){ // mumu/ee
243 if(fabs(dilepton->Mass()-91.1876) > 15.0) passCut[3] = true;
244 if(METdeltaPhilEt > 35) passCut[4] = true;
245 }
246 else if(finalstateType == 12) { // emu
247 passCut[3] = true;
248 if(METdeltaPhilEt > 20) passCut[4] = true;
249 }
250
251 //*********************************************************************************************
252 //Make Selection Histograms. Number of events passing each level of cut
253 //*********************************************************************************************
254 bool passAllCuts = true;
255 for(int c=0; c<nCuts; c++) passAllCuts = passAllCuts & passCut[c];
256
257 //Cut Selection Histograms
258 fHWWSelection->Fill(-1,NNLOWeight->GetVal());
259 if (finalstateType == 10 )
260 fHWWToMuMuSelection->Fill(-1,NNLOWeight->GetVal());
261 else if(finalstateType == 11 )
262 fHWWToEESelection->Fill(-1,NNLOWeight->GetVal());
263 else if(finalstateType == 12 )
264 fHWWToEMuSelection->Fill(-1,NNLOWeight->GetVal());
265
266 for (int k=0;k<nCuts;k++) {
267 bool pass = true;
268 bool passPreviousCut = true;
269 for (int p=0;p<=k;p++) {
270 pass = (pass && passCut[p]);
271 if (p<k)
272 passPreviousCut = (passPreviousCut&& passCut[p]);
273 }
274
275 if (pass) {
276 fHWWSelection->Fill(k,NNLOWeight->GetVal());
277 if (finalstateType == 10 )
278 fHWWToMuMuSelection->Fill(k,NNLOWeight->GetVal());
279 else if(finalstateType == 11)
280 fHWWToEESelection->Fill(k,NNLOWeight->GetVal());
281 else if(finalstateType == 12)
282 fHWWToEMuSelection->Fill(k,NNLOWeight->GetVal());
283 }
284 }
285
286 //*****************************************************************************************
287 //Make Preselection Histograms
288 //*****************************************************************************************
289 fLeptonEta->Fill(CleanLeptons->At(0)->Eta(),NNLOWeight->GetVal());
290 fLeptonEta->Fill(CleanLeptons->At(1)->Eta(),NNLOWeight->GetVal());
291 fLeptonPtMax->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
292 fLeptonPtMin->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
293 fMetPtHist->Fill(caloMet->Pt(),NNLOWeight->GetVal());
294 fMetPhiHist->Fill(caloMet->Phi(),NNLOWeight->GetVal());
295 fDeltaPhiLeptons->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
296 fDeltaEtaLeptons->Fill(deltaEtaLeptons,NNLOWeight->GetVal());
297 fDileptonMass->Fill(dilepton->Mass(),NNLOWeight->GetVal());
298
299 //*********************************************************************************************
300 // N-1 Histograms
301 //*********************************************************************************************
302 bool pass;;
303
304 //N Jet Veto
305 pass = true;
306 for (int k=0;k<nCuts;k++) {
307 if (k != 5) {
308 pass = (pass && passCut[k]);
309 }
310 }
311 if (pass) {
312 fNCentralJets_NMinusOne->Fill(nCentralJets,NNLOWeight->GetVal());
313 }
314
315 // Final Met Cut
316 pass = true;
317 for (int k=0;k<nCuts;k++) {
318 if (k != 4) {
319 pass = (pass && passCut[k]);
320 }
321 }
322 if (pass) {
323 fMetPtHist_NMinusOne->Fill(caloMet->Pt(),NNLOWeight->GetVal());
324 }
325
326 // dilepton mass
327 pass = true;
328 for (int k=0;k<nCuts;k++) {
329 if (k != 2 && k != 3)
330 pass = (pass && passCut[k]);
331 }
332 if (pass) {
333 fDileptonMass_NMinusOne->Fill(dilepton->Mass(),NNLOWeight->GetVal());
334 }
335
336 // Lepton Pt Max, Lepton Pt Min, DeltaPhiLeptons
337 pass = true;
338 for (int k=0;k<nCuts;k++) {
339 pass = (pass && passCut[k]);
340 }
341 if (pass) {
342 fLeptonPtMax_NMinusOne->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
343 fLeptonPtMin_NMinusOne->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
344 fDeltaPhiLeptons_NMinusOne->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
345 }
346
347 // NSoftMuons
348 pass = true;
349 for (int k=0;k<nCuts;k++) {
350 if (k != 6)
351 pass = (pass && passCut[k]);
352 }
353 if (pass) {
354 fNSoftMuonsHist_NMinusOne->Fill(SoftMuons->GetEntries(),NNLOWeight->GetVal());
355 }
356
357 //*********************************************************************************************
358 //Plots after all Cuts
359 //*********************************************************************************************
360 if (passAllCuts) {
361 fMinDeltaPhiLeptonMet_afterCuts->Fill(minDeltaPhiMetLepton,NNLOWeight->GetVal());
362 fMtLepton1_afterCuts->Fill(mTW[0],NNLOWeight->GetVal());
363 fMtLepton2_afterCuts->Fill(mTW[1],NNLOWeight->GetVal());
364 fMtHiggs_afterCuts->Fill(mtHiggs,NNLOWeight->GetVal());
365 fLeptonPtPlusMet_afterCuts->Fill(CleanLeptons->At(0)->Pt()+CleanLeptons->At(1)->Pt()+caloMet->Pt(),NNLOWeight->GetVal());
366 }
367
368 delete dilepton;
369 delete SoftMuons;
370 return;
371 }
372
373 //--------------------------------------------------------------------------------------------------
374 void HwwExampleAnalysisMod::SlaveTerminate()
375 {
376
377 // Run finishing code on the computer (slave) that did the analysis. For this
378 // module, we dont do anything here.
379
380 }
381 //--------------------------------------------------------------------------------------------------
382 void HwwExampleAnalysisMod::Terminate()
383 {
384 // Run finishing code on the client computer. For this module, we dont do
385 // anything here.
386
387 }