ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/SelMods/src/HwwExampleAnalysisMod.cc
Revision: 1.3
Committed: Tue Oct 5 06:42:48 2010 UTC (14 years, 7 months ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.2: +0 -3 lines
Log Message:
small update

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