ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SelMods/src/ttEvtSelMod.cc
Revision: 1.1
Committed: Tue Oct 14 06:13:54 2008 UTC (16 years, 6 months ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
Log Message:
Start of MitPhysics

File Contents

# Content
1 // $Id: ttEvtSelMod.cc,v 1.3 2008/10/10 10:54:13 ceballos Exp $
2
3 #include "MitPhysics/SelMods/interface/ttEvtSelMod.h"
4 #include <TH1D.h>
5 #include <TH2D.h>
6 #include "MitAna/DataTree/interface/Names.h"
7 #include "MitAna/DataCont/interface/ObjArray.h"
8 #include "MitCommon/MathTools/interface/MathUtils.h"
9
10 using namespace mithep;
11 ClassImp(mithep::ttEvtSelMod)
12
13 //--------------------------------------------------------------------------------------------------
14 ttEvtSelMod::ttEvtSelMod(const char *name, const char *title) :
15 BaseMod(name,title),
16 fPrintDebug(false),
17 fMetName(Names::gkCaloMetBrn),
18 fMuonName(Names::gkMuonBrn),
19 fCleanJetsName(Names::gkCleanJetsName),
20 fMCLeptonsName(Names::gkMCLeptonsName),
21 fMCQuarksName(Names::gkMCQuarksName),
22 fMet(0),
23 fMuons(0),
24 fNEventsProcessed(0)
25 {
26 // Constructor.
27 }
28
29 //--------------------------------------------------------------------------------------------------
30 void ttEvtSelMod::Begin()
31 {
32 // Run startup code on the client machine. For this module, we dont do
33 // anything here.
34 }
35
36 //--------------------------------------------------------------------------------------------------
37 void ttEvtSelMod::Process()
38 {
39 // Process entries of the tree. For this module, we just load the branches and
40 fNEventsProcessed++;
41
42 if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
43 time_t systime;
44 systime = time(NULL);
45 cerr << endl << "ttEvtSelMod : Process Event " << fNEventsProcessed << " Time: " << ctime(&systime) << endl;
46 }
47
48 //Get Generator Level information for matching
49 //ObjArray<MCParticle> *GenLeptons = dynamic_cast<ObjArray<MCParticle>* > (FindObjThisEvt(fMCLeptonsName.Data()));
50 ObjArray<MCParticle> *GenQuarks = dynamic_cast<ObjArray<MCParticle>* > (FindObjThisEvt(fMCQuarksName.Data()));
51
52 //Obtain all the good objects from the event cleaning module
53 ObjArray<Electron> *CleanElectrons = dynamic_cast<ObjArray<Electron>* >
54 (FindObjThisEvt(Names::gkCleanElectronsName));
55 ObjArray<Muon> *CleanMuons = dynamic_cast<ObjArray<Muon>* >
56 (FindObjThisEvt(Names::gkCleanMuonsName));
57 ObjArray<Jet> *CleanJets = dynamic_cast<ObjArray<Jet>* >
58 (FindObjThisEvt(fCleanJetsName.Data()));
59
60 LoadBranch(fMetName);
61 Met *caloMet = fMet->At(0);
62
63 vector<ChargedParticle*> leptons;
64 vector<string> leptonType;
65
66 LoadBranch(fMuonName);
67 ObjArray<Muon> *DirtyMuons = new ObjArray<Muon>;
68 for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
69 Muon *mu = fMuons->At(i);
70 if(!mu->GlobalTrk()) continue;
71 if(mu->Pt() < 5.0) continue;
72
73 bool isCleanMuon = false;
74 for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
75 if(fMuons->At(i) == CleanMuons->At(j)) isCleanMuon = true;
76 }
77 if(isCleanMuon == false) DirtyMuons->Add(mu);
78 }
79
80 // Make lepton vector from muons and electrons
81 for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
82 leptons.push_back(CleanMuons->At(j));
83 leptonType.push_back("mu");
84 }
85 for (UInt_t j=0; j<CleanElectrons->GetEntries(); j++) {
86 leptons.push_back(CleanElectrons->At(j));
87 leptonType.push_back("e");
88 }
89
90 // Sort the Leptons by Pt
91 for(UInt_t i=0; i<leptons.size(); i++){
92 for(UInt_t j=i+1; j<leptons.size(); j++){
93 if(leptons[i]->Pt() < leptons[j]->Pt()) {
94 //swap i and j
95 ChargedParticle* templepton = leptons[i];
96 leptons[i] = leptons[j];
97 leptons[j] = templepton;
98 string templeptonType = leptonType[i];
99 leptonType[i] = leptonType[j];
100 leptonType[j] = templeptonType;
101 }
102 }
103 }
104 if (fPrintDebug) {
105 cout << "Check Lepton Sort\n";
106 for(UInt_t i=0; i<leptons.size(); i++)
107 cout << leptons[i]->Pt() << endl;
108 }
109
110 // Pt requirements && MET requirements
111 if (leptons.size() >= 2 && caloMet->Pt() >= 0 &&
112 leptons[0]->Pt() > 20 && leptons[1]->Pt() > 15){
113
114 CompositeParticle *dilepton = new CompositeParticle();
115 dilepton->AddDaughter(leptons[0]);
116 dilepton->AddDaughter(leptons[1]);
117 // Charge of the leptons should be opposite && minimum mass requirement
118 if (dilepton->Charge() == 0 && dilepton->Mass() > 12){
119
120 // Sort and count the number of central Jets for vetoing
121 vector<Jet*> sortedJets;
122 int nCentralJets = 0;
123 for(UInt_t i=0; i<CleanJets->GetEntries(); i++){
124 if(fabs(CleanJets->At(i)->Eta()) < 2.5){
125 nCentralJets++;
126 Jet* jet_f = new Jet(CleanJets->At(i)->Px()*CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale(),
127 CleanJets->At(i)->Py()*CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale(),
128 CleanJets->At(i)->Pz()*CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale(),
129 CleanJets->At(i)->E() *CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale());
130 sortedJets.push_back(jet_f);
131 }
132 }
133 for(UInt_t i=0; i<sortedJets.size(); i++){
134 for(UInt_t j=i+1; j<sortedJets.size(); j++){
135 if(sortedJets[i]->Pt() < sortedJets[j]->Pt()) {
136 //swap i and j
137 Jet* tempjet = sortedJets[i];
138 sortedJets[i] = sortedJets[j];
139 sortedJets[j] = tempjet;
140 }
141 }
142 }
143
144 // Delta phi between the 2 leptons in degrees
145 double deltaPhiLeptons = MathUtils::DeltaPhi(leptons[0]->Phi(),
146 leptons[1]->Phi())* 180./TMath::Pi();
147
148 double deltaEtaLeptons = abs(leptons[0]->Eta() - leptons[1]->Eta());
149
150 double deltaPhiDileptonMet = MathUtils::DeltaPhi(caloMet->Phi(),
151 dilepton->Phi())* 180./TMath::Pi();
152
153 double mtHiggs = TMath::Sqrt(2.0*dilepton->Pt() * caloMet->Pt()*
154 (1.0 - cos(deltaPhiDileptonMet * TMath::Pi()/180.0)));
155
156 // Angle between MET and closest lepton
157 double deltaPhiMetLepton[2] = {MathUtils::DeltaPhi(caloMet->Phi(), leptons[0]->Phi()),
158 MathUtils::DeltaPhi(caloMet->Phi(), leptons[1]->Phi())};
159
160 double minDeltaPhiMetLepton = (deltaPhiMetLepton[0] < deltaPhiMetLepton[1])?
161 deltaPhiMetLepton[0]:deltaPhiMetLepton[1];
162 minDeltaPhiMetLepton = minDeltaPhiMetLepton * 180./TMath::Pi();
163
164 double METdeltaPhilEt = caloMet->Pt();
165 if(minDeltaPhiMetLepton < 90.0)
166 METdeltaPhilEt = METdeltaPhilEt * sin(minDeltaPhiMetLepton * TMath::Pi()/180.);
167
168 double mTW[2] = {TMath::Sqrt(2.0*leptons[0]->Pt()*caloMet->Pt()*
169 (1.0 - cos(deltaPhiMetLepton[0]))),
170 TMath::Sqrt(2.0*leptons[1]->Pt()*caloMet->Pt()*
171 (1.0 - cos(deltaPhiMetLepton[1])))};
172
173 int pairType = -1;
174 if (leptonType[0] == "mu" && leptonType[1] == "mu" )
175 pairType = 0;
176 else if(leptonType[0] == "e" && leptonType[1] == "e")
177 pairType = 1;
178 else if((leptonType[0] == "e" && leptonType[1] == "mu") ||
179 (leptonType[0] == "mu" && leptonType[1] == "e"))
180 pairType = 2;
181 else {
182 cout << "Hey, this is not possible, leptonTypes: "
183 << leptonType[0] << " - "
184 << leptonType[1] << endl;
185 }
186
187 hDttPresel[ 0+100*pairType]->Fill((double)leptons.size());
188 // No more than 2 isolated good leptons
189 if (leptons.size() > 2) return;
190
191 hDttPresel[ 1+100*pairType]->Fill(caloMet->Pt());
192 hDttPresel[ 2+100*pairType]->Fill(leptons[0]->Pt());
193 hDttPresel[ 3+100*pairType]->Fill(leptons[1]->Pt());
194 hDttPresel[ 4+100*pairType]->Fill(dilepton->Mass());
195 hDttPresel[ 5+100*pairType]->Fill((double)nCentralJets);
196 hDttPresel[ 6+100*pairType]->Fill(deltaEtaLeptons);
197 hDttPresel[ 7+100*pairType]->Fill(deltaPhiLeptons);
198 hDttPresel[ 8+100*pairType]->Fill(mtHiggs);
199 hDttPresel[ 9+100*pairType]->Fill(minDeltaPhiMetLepton);
200 hDttPresel[10+100*pairType]->Fill(METdeltaPhilEt);
201 hDttPresel[11+100*pairType]->Fill(caloMet->MetSig());
202 hDttPresel[12+100*pairType]->Fill(caloMet->SumEt());
203 hDttPresel[13+100*pairType]->Fill(TMath::Max(mTW[0],mTW[1]));
204 hDttPresel[14+100*pairType]->Fill(TMath::Min(mTW[0],mTW[1]));
205 hDttPresel[15+100*pairType]->Fill(leptons[0]->Pt()+leptons[1]->Pt()+caloMet->Pt());
206 for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
207 hDttPresel[16+100*pairType]->Fill(TMath::Min(fabs(CleanMuons->At(j)->GlobalTrk()->D0()),0.3999));
208 }
209 for (UInt_t j=0; j<CleanElectrons->GetEntries(); j++) {
210 hDttPresel[17+100*pairType]->Fill(TMath::Min(fabs(CleanElectrons->At(j)->GsfTrk()->D0()),0.3999));
211 }
212
213 // Study for genuine tt->WbWb->2l2nubb
214 if(nCentralJets >= 2){
215 hDttPresel[18+100*pairType]->Fill(TMath::Max(fabs(sortedJets[0]->Eta()),
216 fabs(sortedJets[1]->Eta())));
217 hDttPresel[19+100*pairType]->Fill(TMath::Min(fabs(sortedJets[0]->Eta()),
218 fabs(sortedJets[1]->Eta())));
219 hDttPresel[20+100*pairType]->Fill(MathUtils::DeltaPhi(sortedJets[0]->Phi(),sortedJets[1]->Phi())
220 * 180.0 / TMath::Pi());
221 hDttPresel[21+100*pairType]->Fill(fabs(sortedJets[0]->Eta()-sortedJets[1]->Eta()));
222 hDttPresel[22+100*pairType]->Fill(MathUtils::DeltaR(sortedJets[0]->Phi(),sortedJets[0]->Eta(),
223 sortedJets[1]->Phi(),sortedJets[1]->Eta()));
224 hDttPresel[23+100*pairType]->Fill(TMath::Max(fabs(sortedJets[0]->Pt()),
225 fabs(sortedJets[1]->Pt())));
226 hDttPresel[24+100*pairType]->Fill(TMath::Min(fabs(sortedJets[0]->Pt()),
227 fabs(sortedJets[1]->Pt())));
228 CompositeParticle *dijet = new CompositeParticle();
229 dijet->AddDaughter(sortedJets[0]);
230 dijet->AddDaughter(sortedJets[1]);
231 hDttPresel[25+100*pairType]->Fill(dijet->Mass());
232 hDttPresel[26+100*pairType]->Fill(dijet->Pt());
233
234 // Study energy scale corrections
235 if(GenQuarks->GetEntries() >= 2){
236 double deltarAux[2];
237 double deltaPt;
238 deltarAux[0] = MathUtils::DeltaR(sortedJets[0]->Phi(),sortedJets[0]->Eta(),
239 GenQuarks->At(0)->Phi(),GenQuarks->At(0)->Eta());
240 deltarAux[1] = MathUtils::DeltaR(sortedJets[0]->Phi(),sortedJets[0]->Eta(),
241 GenQuarks->At(1)->Phi(),GenQuarks->At(1)->Eta());
242 hDttPresel[27+100*pairType]->Fill(TMath::Min(deltarAux[0],deltarAux[1]));
243 deltaPt = -999;
244 if (deltarAux[0] < deltarAux[1] && deltarAux[0] < 1.0){
245 deltaPt = (sortedJets[0]->Pt()-GenQuarks->At(0)->Pt())/sortedJets[0]->Pt();
246 }
247 else if(deltarAux[1] < deltarAux[0] && deltarAux[1] < 1.0){
248 deltaPt = (sortedJets[0]->Pt()-GenQuarks->At(1)->Pt())/sortedJets[0]->Pt();
249 }
250 if(deltaPt != -999) hDttPresel[28+100*pairType]->Fill(deltaPt);
251 deltarAux[0] = MathUtils::DeltaR(sortedJets[1]->Phi(),sortedJets[1]->Eta(),
252 GenQuarks->At(0)->Phi(),GenQuarks->At(0)->Eta());
253 deltarAux[1] = MathUtils::DeltaR(sortedJets[1]->Phi(),sortedJets[1]->Eta(),
254 GenQuarks->At(1)->Phi(),GenQuarks->At(1)->Eta());
255 hDttPresel[29+100*pairType]->Fill(TMath::Min(deltarAux[0],deltarAux[1]));
256 deltaPt = -999;
257 if (deltarAux[0] < deltarAux[1] && deltarAux[0] < 1.0){
258 deltaPt = (sortedJets[1]->Pt()-GenQuarks->At(0)->Pt())/sortedJets[1]->Pt();
259 }
260 else if(deltarAux[1] < deltarAux[0] && deltarAux[1] < 1.0){
261 deltaPt = (sortedJets[1]->Pt()-GenQuarks->At(1)->Pt())/sortedJets[1]->Pt();
262 }
263 if(deltaPt != -999) hDttPresel[30+100*pairType]->Fill(deltaPt);
264 }
265 delete dijet;
266 } // Njets >=2
267
268 // Study for additional muons
269 hDttPresel[31+100*pairType]->Fill((double)DirtyMuons->GetEntries());
270 for (UInt_t i=0; i<DirtyMuons->GetEntries(); ++i) {
271 Muon *mu = DirtyMuons->At(i);
272 hDttPresel[32+100*pairType]->Fill(mu->Pt());
273 hDttPresel[33+100*pairType]->Fill(TMath::Min(fabs(mu->GlobalTrk()->D0()),0.3999));
274 }
275 if(DirtyMuons->GetEntries() > 0){
276 hDttPresel[34+100*pairType]->Fill((double)nCentralJets);
277 }
278
279 // Study dijet mass
280 if(nCentralJets >= 2){
281 double dijetMass = 999.;
282 int indDiJet[2] = {-1, -1};
283 for(UInt_t i=0; i<sortedJets.size(); i++){
284 for(UInt_t j=i+1; j<sortedJets.size(); j++){
285 CompositeParticle *dijetTemp = new CompositeParticle();
286 dijetTemp->AddDaughter(sortedJets[i]);
287 dijetTemp->AddDaughter(sortedJets[j]);
288 if(fabs(dijetTemp->Mass()-91.18) < fabs(dijetMass-91.18)){
289 dijetMass = dijetTemp->Mass();
290 indDiJet[0] = i;
291 indDiJet[1] = j;
292 }
293 delete dijetTemp;
294 }
295 }
296 hDttPresel[35+100*pairType]->Fill(dijetMass);
297 if(indDiJet[0] == 0 && indDiJet[1] == 1)
298 hDttPresel[36+100*pairType]->Fill(dijetMass);
299 else
300 hDttPresel[37+100*pairType]->Fill(dijetMass);
301 hDttPresel[38+100*pairType]->Fill((double)(indDiJet[0]+10*indDiJet[1]));
302 }
303 for(UInt_t i=0; i<sortedJets.size(); i++) delete sortedJets[i];
304 } // q1+q2==0 && mass_ll>12
305 delete dilepton;
306 } // Minimun Pt, Nleptons and MET requirements
307 delete DirtyMuons;
308 leptons.clear();
309 }
310 //--------------------------------------------------------------------------------------------------
311 void ttEvtSelMod::SlaveBegin()
312 {
313 // Run startup code on the computer (slave) doing the actual analysis. Here,
314 // we typically initialize histograms and other analysis objects and request
315 // branches. For this module, we request a branch of the MitTree.
316
317 ReqBranch(fMetName, fMet);
318 ReqBranch(fMuonName, fMuons);
319
320 char sb[200];
321 for(int j=0; j<3; j++){
322 int ind = 100 * j;
323 sprintf(sb,"hDttPresel_%d",ind+0); hDttPresel[ind+0] = new TH1D(sb,sb,10,-0.5,9.5);
324 sprintf(sb,"hDttPresel_%d",ind+1); hDttPresel[ind+1] = new TH1D(sb,sb,150,0.0,300.);
325 sprintf(sb,"hDttPresel_%d",ind+2); hDttPresel[ind+2] = new TH1D(sb,sb,100,0.0,200.);
326 sprintf(sb,"hDttPresel_%d",ind+3); hDttPresel[ind+3] = new TH1D(sb,sb,100,0.0,200.);
327 sprintf(sb,"hDttPresel_%d",ind+4); hDttPresel[ind+4] = new TH1D(sb,sb,150,0.0,300.);
328 sprintf(sb,"hDttPresel_%d",ind+5); hDttPresel[ind+5] = new TH1D(sb,sb,10,-0.5,9.5);
329 sprintf(sb,"hDttPresel_%d",ind+6); hDttPresel[ind+6] = new TH1D(sb,sb,50,0.0,5.);
330 sprintf(sb,"hDttPresel_%d",ind+7); hDttPresel[ind+7] = new TH1D(sb,sb,90,0.0,180.);
331 sprintf(sb,"hDttPresel_%d",ind+8); hDttPresel[ind+8] = new TH1D(sb,sb,150,0.0,300.);
332 sprintf(sb,"hDttPresel_%d",ind+9); hDttPresel[ind+9] = new TH1D(sb,sb,90,0.0,180.);
333 sprintf(sb,"hDttPresel_%d",ind+10); hDttPresel[ind+10] = new TH1D(sb,sb,150,0.0,300.);
334 sprintf(sb,"hDttPresel_%d",ind+11); hDttPresel[ind+11] = new TH1D(sb,sb,100,0.0,20.);
335 sprintf(sb,"hDttPresel_%d",ind+12); hDttPresel[ind+12] = new TH1D(sb,sb,200,0.0,800.);
336 sprintf(sb,"hDttPresel_%d",ind+13); hDttPresel[ind+13] = new TH1D(sb,sb,200,0.0,300.);
337 sprintf(sb,"hDttPresel_%d",ind+14); hDttPresel[ind+14] = new TH1D(sb,sb,200,0.0,300.);
338 sprintf(sb,"hDttPresel_%d",ind+15); hDttPresel[ind+15] = new TH1D(sb,sb,250,0.0,500.);
339 sprintf(sb,"hDttPresel_%d",ind+16); hDttPresel[ind+16] = new TH1D(sb,sb,200,0.0,0.4);
340 sprintf(sb,"hDttPresel_%d",ind+17); hDttPresel[ind+17] = new TH1D(sb,sb,200,0.0,0.4);
341 sprintf(sb,"hDttPresel_%d",ind+18); hDttPresel[ind+18] = new TH1D(sb,sb,25,0.0,2.5);
342 sprintf(sb,"hDttPresel_%d",ind+19); hDttPresel[ind+19] = new TH1D(sb,sb,25,0.0,2.5);
343 sprintf(sb,"hDttPresel_%d",ind+20); hDttPresel[ind+20] = new TH1D(sb,sb,90,0.0,180.);
344 sprintf(sb,"hDttPresel_%d",ind+21); hDttPresel[ind+21] = new TH1D(sb,sb,50,0.0, 5.0);
345 sprintf(sb,"hDttPresel_%d",ind+22); hDttPresel[ind+22] = new TH1D(sb,sb,50,0.0,10.0);
346 sprintf(sb,"hDttPresel_%d",ind+23); hDttPresel[ind+23] = new TH1D(sb,sb,100,0.0,200.0);
347 sprintf(sb,"hDttPresel_%d",ind+24); hDttPresel[ind+24] = new TH1D(sb,sb,100,0.0,200.0);
348 sprintf(sb,"hDttPresel_%d",ind+25); hDttPresel[ind+25] = new TH1D(sb,sb,100,0.0,400.0);
349 sprintf(sb,"hDttPresel_%d",ind+26); hDttPresel[ind+26] = new TH1D(sb,sb,100,0.0,400.0);
350 sprintf(sb,"hDttPresel_%d",ind+27); hDttPresel[ind+27] = new TH1D(sb,sb,200,0.0,10.0);
351 sprintf(sb,"hDttPresel_%d",ind+28); hDttPresel[ind+28] = new TH1D(sb,sb,80,-2.0,2.0);
352 sprintf(sb,"hDttPresel_%d",ind+29); hDttPresel[ind+29] = new TH1D(sb,sb,200,0.0,10.0);
353 sprintf(sb,"hDttPresel_%d",ind+30); hDttPresel[ind+30] = new TH1D(sb,sb,80,-2.0,2.0);
354 sprintf(sb,"hDttPresel_%d",ind+31); hDttPresel[ind+31] = new TH1D(sb,sb,10,-0.5,9.5);
355 sprintf(sb,"hDttPresel_%d",ind+32); hDttPresel[ind+32] = new TH1D(sb,sb,200,0.0,50.0);
356 sprintf(sb,"hDttPresel_%d",ind+33); hDttPresel[ind+33] = new TH1D(sb,sb,200,0.0,0.4);
357 sprintf(sb,"hDttPresel_%d",ind+34); hDttPresel[ind+34] = new TH1D(sb,sb,10,-0.5,9.5);
358 sprintf(sb,"hDttPresel_%d",ind+35); hDttPresel[ind+35] = new TH1D(sb,sb,100,0.0,400.0);
359 sprintf(sb,"hDttPresel_%d",ind+36); hDttPresel[ind+36] = new TH1D(sb,sb,100,0.0,400.0);
360 sprintf(sb,"hDttPresel_%d",ind+37); hDttPresel[ind+37] = new TH1D(sb,sb,100,0.0,400.0);
361 sprintf(sb,"hDttPresel_%d",ind+38); hDttPresel[ind+38] = new TH1D(sb,sb,100,-0.5,99.5);
362 }
363
364 for(int i=0; i<39; i++){
365 for(int j=0; j<3; j++){
366 AddOutput(hDttPresel[i+j*100]);
367 }
368 }
369 }
370
371 //--------------------------------------------------------------------------------------------------
372 void ttEvtSelMod::SlaveTerminate()
373 {
374 // Run finishing code on the computer (slave) that did the analysis
375 }
376
377 //--------------------------------------------------------------------------------------------------
378 void ttEvtSelMod::Terminate()
379 {
380 // Run finishing code on the client computer
381 }