ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitAna/TreeMod/src/ObjectCleaningMod.cc
Revision: 1.7
Committed: Tue Oct 14 05:12:49 2008 UTC (16 years, 7 months ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.6: +1 -1 lines
State: FILE REMOVED
Log Message:
Moved to MitPhysics/Mods

File Contents

# Content
1 // $Id: ObjectCleaningMod.cc,v 1.6 2008/09/30 16:38:20 sixie Exp $
2
3 #include "MitAna/TreeMod/interface/ObjectCleaningMod.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 "MitAna/Utils/interface/IsolationTools.h"
9 #include "MitCommon/MathTools/interface/MathUtils.h"
10 #include "MitAna/Utils/interface/IsolationTools.h"
11
12 using namespace mithep;
13
14 ClassImp(mithep::ObjectCleaningMod)
15
16 //--------------------------------------------------------------------------------------------------
17 ObjectCleaningMod::ObjectCleaningMod(const char *name, const char *title) :
18 BaseMod(name,title),
19 fPrintDebug(false),
20 fMCPartName(Names::gkMCPartBrn),
21 fMuonName(Names::gkMuonBrn),
22 fElectronName(Names::gkElectronBrn),
23 fJetName(Names::gkCaloJetBrn),
24 fSCJetName(Names::gkSC5JetBrn),
25 fMetName(Names::gkCaloMetBrn),
26 fGoodElectronsName("GoodElectrons"),
27 fGoodMuonsName("GoodMuons"),
28 fGoodCentralJetsName("GoodCentralJets" ),
29 fGoodForwardJetsName("GoodForwardJets" ),
30 fGoodCentralSCJetsName("GoodCentralSCJets" ),
31 fGoodForwardSCJetsName("GoodForwardSCJets" ),
32 fGenLeptonsName("GenLeptons"),
33 fParticles(0),
34 fMuons(0),
35 fElectrons(0)
36 {
37 // Constructor.
38 }
39
40 //--------------------------------------------------------------------------------------------------
41 void ObjectCleaningMod::Begin()
42 {
43 // Run startup code on the client machine. For this module, we dont do
44 // anything here.
45 }
46
47 //--------------------------------------------------------------------------------------------------
48 void ObjectCleaningMod::Process()
49 {
50 // Process entries of the tree. For this module, we just load the branches and
51 //output debug info or not
52
53 fNEventsProcessed++;
54
55 if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
56 time_t systime;
57 systime = time(NULL);
58
59 cerr << endl << "ObjectCleaningMod : Process Event " << fNEventsProcessed << " Time: " << ctime(&systime) << endl;
60 }
61
62 //Muons
63 ObjArray<Muon> *GoodMuons = new ObjArray<Muon>;
64 LoadBranch(fMuonName);
65
66 for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
67 Muon *mu = fMuons->At(i);
68 fAllMuonPtHist->Fill(mu->Pt());
69 fAllMuonEtaHist->Fill(mu->Eta());
70
71 Double_t MuonClass = -1;
72 //Find Class of Muons: Global(0), Standalone(1), Tracker(2)
73 //Double_t MuonClass = -1;
74
75 if (mu->GlobalTrk())
76 MuonClass = 0;
77 else if (mu->StandaloneTrk())
78
79 MuonClass = 1;
80 else if (mu->TrackerTrk())
81 MuonClass = 2;
82
83 //Define the ID Cuts
84 const int nCuts = 4;
85 double cutValue[nCuts] = {0.1, 3.0, 3.0, 1.5 };
86 bool passCut[nCuts] = {false, false, false, false};
87
88 double muonD0 = -0.05;
89
90 muonD0 = mu->BestTrk()->D0();
91 if(muonD0 < cutValue[0] && MuonClass == 0 )
92 passCut[0] = true;
93 if(mu->IsoR03SumPt() < cutValue[1]) passCut[1] = true;
94 if(mu->IsoR03EmEt() +
95 mu->IsoR03HadEt() < cutValue[2]) passCut[2] = true;
96 if(mu->Pt() > 5)
97 passCut[3] = true;
98
99 // if(mu->Trk()->NHits() > 6)
100 // passCut[4] = true;
101
102 // if(mu->Trk()->Chi2() / mu->Trk()->Ndof() < 5 )
103 // passCut[5] = true;
104
105 // Final decision
106 bool allCuts = true;
107 for(int c=0; c<nCuts; c++) {
108 allCuts = allCuts & passCut[c];
109 }
110
111 //make muon ID selection histograms
112 fMuonSelection->Fill(-1);
113 //Fill the rest of the muon selection histograms
114 for (int k=0;k<3;k++) {
115 bool pass = true;
116 for (int p=0;p<=k;p++)
117 pass = (pass && passCut[p]);
118
119 if (pass) {
120 fMuonSelection->Fill(k);
121 }
122 }
123 //Fill histogram for the Good Muons
124 if ( allCuts
125 && abs(mu->Eta()) < 2.5
126 ) {
127 fGoodMuonPtHist->Fill(mu->Pt());
128 fGoodMuonEtaHist->Fill(mu->Eta());
129 GoodMuons->Add(mu);
130 }
131 }
132
133
134 //Get Electrons
135 LoadBranch(fElectronName);
136
137 //we have to use a vector first and then fill the ObjArray with the vector
138 //contents because ObJArray does not allow removal of duplicates.
139 vector<Electron*> GoodElectronsVector;
140 vector<Electron*> GoodLooseElectronsVector;
141 vector<Electron*> GoodTightElectronsVector;
142 vector<Electron*> GoodLikelihoodIDElectronsVector;
143
144 for (UInt_t i=0; i<fElectrons->GetEntries(); ++i) {
145 Electron *e = fElectrons->At(i);
146 fAllElectronPtHist->Fill(e->Pt());
147 fAllElectronEtaHist->Fill(e->Eta());
148
149 //from RecoEgamma/ElectronIdentification/src/CurBasedElectronID.cc : CMSSW 2_1_0
150 Int_t InEndcapOrBarrel = (fabs(e->Eta()) < 1.479)? 0 : 1;
151 double sigmaEtaEta = e->CovEtaEta();
152 //have to correct sigma_etaeta dependance on eta in the endcap
153 if (InEndcapOrBarrel == 1) {
154 sigmaEtaEta = sigmaEtaEta - 0.02*(fabs(e->Eta()) - 2.3);
155 e->SetCovEtaEta(e->CovEtaEta() - 0.02*(fabs(e->Eta()) - 2.3));
156 }
157
158 //Calculate the electron category for determining which cuts to make
159 Int_t ElectronCategory;
160 Double_t fBrem = (e->PIn() - e->POut())/e->PIn();
161 if((fabs(e->Eta())<1.479 && fBrem<0.06) || (fabs(e->Eta())>1.479 && fBrem<0.1))
162 ElectronCategory=1;
163 else if (e->ESuperClusterOverP() < 1.2 && e->ESuperClusterOverP() > 0.8)
164 ElectronCategory=0;
165 else
166 ElectronCategory=2;
167
168 // Final decision
169 bool allCuts = true;
170 allCuts = e->PassTightID();
171
172 //Check whether it overlaps with a good muon.
173 bool isMuonOverlap = false;
174 for (UInt_t j=0; j<GoodMuons->GetEntries();j++) {
175 double deltaR = MathUtils::DeltaR(GoodMuons->At(j)->Mom(), e->Mom());
176 if (deltaR < 0.1) {
177 isMuonOverlap = true;
178 break;
179 }
180 }
181
182 //Check whether it overlaps with another electron candidate
183 bool isElectronOverlap = false;
184 for (UInt_t j=0; j<GoodElectronsVector.size(); j++) {
185 double deltaR = MathUtils::DeltaR(GoodElectronsVector[j]->Mom(), e->Mom());
186
187 if (deltaR < 0.1) {
188 isElectronOverlap = true;
189
190 //if there's an overlap and the new one being considered has E/P closer to 1.0
191 //then replace the overlapping one with this new one because it's better.
192 //Once we get super cluster info, we should make sure the superclusters match
193 //before declaring it is a duplicate
194 //can have one SC matched to many tracks. or one track matched to many SC's.
195 //we should cover all the cases. see SUSYAnalyzer...
196 //Si: This will be unnecessary very soon. It is to be removed in reconstruction.
197 if (abs(GoodElectronsVector[j]->ESuperClusterOverP() - 1) >
198 abs(e->ESuperClusterOverP() - 1)) {
199 GoodElectronsVector[j] = e;
200 }
201 break;
202 }
203 }
204
205 //These are Good Electrons
206 if ( allCuts && !isMuonOverlap && !isElectronOverlap
207 ) {
208 fGoodElectronPtHist->Fill(e->Pt());
209 fGoodElectronEtaHist->Fill(e->Eta());
210 fGoodElectronClassification->Fill(e->Classification());
211 GoodElectronsVector.push_back(fElectrons->At(i));
212 }
213
214 }
215 //fill the electron ObjArray with the contents of the vector
216 //this is necessary because I want to swap out the duplicates, can't be done with ObjArray...
217 ObjArray<Electron> *GoodElectrons = new ObjArray<Electron>;
218 for (UInt_t j=0; j<GoodElectronsVector.size(); j++)
219 GoodElectrons->Add(GoodElectronsVector[j]);
220
221 //Get Jet info
222 ObjArray<Jet> *GoodCentralJets = new ObjArray<Jet>;
223 ObjArray<Jet> *GoodForwardJets = new ObjArray<Jet>;
224 LoadBranch(fJetName);
225 for (UInt_t i=0; i<fJets->GetEntries(); ++i) {
226 Jet *jet = fJets->At(i);
227
228 bool isElectronOverlap = false;
229
230 //Check for overlap with an electron
231 for (UInt_t j=0; j<GoodElectrons->GetEntries(); j++) {
232 double deltaR = MathUtils::DeltaR(GoodElectrons->At(j)->Mom(),jet->Mom());
233 if (deltaR < 0.1) {
234 isElectronOverlap = true;
235 break;
236 }
237 }
238 if (isElectronOverlap) continue; //skip this jet if it was an overlap
239
240
241 fAllJetPtHist->Fill(jet->Pt());
242 fAllJetEtaHist->Fill(jet->Eta());
243
244 const int nCuts = 4;
245 double cutValue[nCuts] = {1.0, 15., 2.5, 0.2};
246 bool passCut[nCuts] = {false, false, false, false};
247
248 passCut[0] = (!isElectronOverlap); //This is supposed to check e/ph overlap
249 if(jet->Et() > cutValue[1]) passCut[1] = true;
250 if(fabs(jet->Eta()) < cutValue[2]) passCut[2] = true;
251 if(jet->Alpha() > cutValue[3] ||
252 jet->Et() > 20.)
253 passCut[3] = true;
254
255 // Final decision
256 bool passAllCentralJetCuts = true;
257 bool passAllForwardJetCuts = true;
258 for(int i=0; i<nCuts; i++) {
259 passAllCentralJetCuts = passAllCentralJetCuts && passCut[i];
260 passAllForwardJetCuts = passAllForwardJetCuts && ((i==2)?(!passCut[i]):passCut[i]);
261 }
262
263 //make electron ID selection histogram
264 fCentralJetSelection->Fill(-1);
265
266 for (int k=0;k<nCuts;k++) {
267 bool pass = true;
268 for (int p=0;p<=k;p++)
269 pass = (pass && passCut[p]);
270
271 if (pass) {
272 fCentralJetSelection->Fill(k);
273 }
274 }
275
276 //Save Good Jets
277 if (passAllCentralJetCuts) {
278 GoodCentralJets->Add(jet);
279 }
280
281 if(passAllForwardJetCuts)
282 GoodForwardJets->Add(jet);
283
284 } //for all jets
285
286
287 //Get Siscone 5 Jet info
288 ObjArray<Jet> *GoodCentralSCJets = new ObjArray<Jet>;
289 ObjArray<Jet> *GoodForwardSCJets = new ObjArray<Jet>;
290 LoadBranch(fSCJetName);
291 for (UInt_t i=0; i<fSCJets->GetEntries(); ++i) {
292 Jet *jet = fSCJets->At(i);
293
294 bool isElectronOverlap = false;
295
296 //Check for overlap with an electron
297 for (UInt_t j=0; j<GoodElectrons->GetEntries(); j++) {
298 double deltaR = MathUtils::DeltaR(GoodElectrons->At(j)->Mom(),jet->Mom());
299 if (deltaR < 0.1) {
300 isElectronOverlap = true;
301 break;
302 }
303 }
304 if (isElectronOverlap) continue; //skip this jet if it was an overlap
305
306 const int nCuts = 4;
307 double cutValue[nCuts] = {1.0, 15., 2.5, 0.2};
308 bool passCut[nCuts] = {false, false, false, false};
309
310 passCut[0] = (!isElectronOverlap); //This is supposed to check e/ph overlap
311 if(jet->Et() > cutValue[1]) passCut[1] = true;
312 if(fabs(jet->Eta()) < cutValue[2]) passCut[2] = true;
313 if(jet->Alpha() > cutValue[3] ||
314 jet->Et() > 20.)
315 passCut[3] = true;
316
317 // Final decision
318 bool passAllCentralJetCuts = true;
319 bool passAllForwardJetCuts = true;
320 for(int i=0; i<nCuts; i++) {
321 passAllCentralJetCuts = passAllCentralJetCuts && passCut[i];
322 passAllForwardJetCuts = passAllForwardJetCuts && ((i==2)?(!passCut[i]):passCut[i]);
323 }
324
325 for (int k=0;k<nCuts;k++) {
326 bool pass = true;
327 for (int p=0;p<=k;p++)
328 pass = (pass && passCut[p]);
329
330 }
331
332 //Save Good SCJets
333 if (passAllCentralJetCuts) {
334 GoodCentralSCJets->Add(jet);
335 }
336
337 if(passAllForwardJetCuts)
338 GoodForwardSCJets->Add(jet);
339
340 } //for all jets
341
342 //Get MET
343 LoadBranch(fMetName);
344
345 Met *met = fMet->At(0); //define the met
346 fMetPtHist->Fill(met->Pt());
347 fMetPhiHist->Fill(met->Phi());
348
349 //Final Summary Debug Output
350 if ( fPrintDebug ) {
351 cerr << "Event Dump: " << fNEventsProcessed << endl;
352
353 //print out event content to text
354 cerr << "Electrons" << endl;
355 for (UInt_t i = 0; i < GoodElectrons->GetEntries(); i++) {
356 cerr << i << " " << GoodElectrons->At(i)->Pt() << " " << GoodElectrons->At(i)->Eta()
357 << " " << GoodElectrons->At(i)->Phi() << " "
358 << GoodElectrons->At(i)->ESuperClusterOverP() << endl;
359 }
360
361 cerr << "Muons" << endl;
362 for (UInt_t i = 0; i < GoodMuons->GetEntries(); i++) {
363 cerr << i << " " << GoodMuons->At(i)->Pt() << " " << GoodMuons->At(i)->Eta()
364 << " " << GoodMuons->At(i)->Phi() << endl;
365 }
366
367 cerr << "Central Jets" << endl;
368 for (UInt_t i = 0; i < GoodCentralJets->GetEntries(); i++) {
369 cerr << i << " " << GoodCentralJets->At(i)->Pt() << " "
370 << GoodCentralJets->At(i)->Eta() << " " << GoodCentralJets->At(i)->Phi() << endl;
371 }
372
373 cerr << "MET" << endl;
374 cerr << met->Pt() << " " << met->Eta() << " "
375 << met->Phi() << endl;
376 }
377
378 //Save Objects for Other Modules to use
379 AddObjThisEvt(GoodElectrons, fGoodElectronsName.Data());
380 AddObjThisEvt(GoodMuons, fGoodMuonsName.Data());
381 AddObjThisEvt(GoodCentralJets, fGoodCentralJetsName.Data());
382 AddObjThisEvt(GoodForwardJets, fGoodForwardJetsName.Data());
383 AddObjThisEvt(GoodCentralSCJets, fGoodCentralSCJetsName.Data());
384 AddObjThisEvt(GoodForwardSCJets, fGoodForwardSCJetsName.Data());
385
386 }
387
388
389 //--------------------------------------------------------------------------------------------------
390 void ObjectCleaningMod::SlaveBegin()
391 {
392 // Run startup code on the computer (slave) doing the actual analysis. Here,
393 // we typically initialize histograms and other analysis objects and request
394 // branches. For this module, we request a branch of the MitTree.
395
396 ReqBranch(fMCPartName, fParticles);
397 ReqBranch(fMuonName, fMuons);
398 ReqBranch(fElectronName, fElectrons);
399 ReqBranch(fJetName, fJets);
400 ReqBranch(fSCJetName, fSCJets);
401 ReqBranch(fMetName, fMet);
402
403
404 fAllMuonPtHist = new TH1D("hAllMuonPtHist",";p_{t};#",200,0.,200.);
405 fAllMuonEtaHist = new TH1D("hAllMuonEtaHist",";#eta;#",100,-5.,5.);
406 fGoodMuonPtHist = new TH1D("hGoodMuonPtHist",";p_{t};#",200,0.,200.);
407 fGoodMuonEtaHist = new TH1D("hGoodMuonEtaHist",";#eta;#",21,-5.,5.);
408 fMuonSelection = new TH1D("hMuonSelection", ";MuonSelection;#",4,-1.5,2.5 ) ;
409 AddOutput(fAllMuonPtHist);
410 AddOutput(fAllMuonEtaHist);
411 AddOutput(fGoodMuonPtHist);
412 AddOutput(fGoodMuonEtaHist);
413 AddOutput(fMuonSelection);
414
415 fAllElectronPtHist = new TH1D("hAllElectronPtHist",";p_{t};#",100,0.,200.);
416 fAllElectronEtaHist = new TH1D("hAllElectronEtaHist",";#eta;#",100,-5.,5.);
417 fGoodElectronPtHist = new TH1D("hGoodElectronPtHist",";p_{t};#",25,0.,200.);
418 fGoodElectronEtaHist = new TH1D("hGoodElectronEtaHist",";#eta;#",21,-5.,5.);
419 fGoodElectronClassification = new TH1D("hGoodElectronClassification",
420 ";Good Electron Classification;#",51,0,50 ) ;
421 AddOutput(fAllElectronPtHist);
422 AddOutput(fAllElectronEtaHist);
423 AddOutput(fGoodElectronPtHist);
424 AddOutput(fGoodElectronEtaHist);
425 AddOutput(fGoodElectronClassification);
426
427 //Jet Plots
428 fAllJetPtHist = new TH1D("hAllJetPtHist",";All Jet p_{t};#",100,0.,200.);
429 fAllJetEtaHist = new TH1D("hAllJetEtaHist",";All Jet #eta;#",160,-8.,8.);
430 fCentralJetSelection = new TH1D("hCentralJetSelection",
431 ";CentralJetSelection;#",5,-1.5,3.5 ) ;
432 AddOutput(fAllJetPtHist);
433 AddOutput(fAllJetEtaHist);
434 AddOutput(fCentralJetSelection);
435
436 //MET Plots
437 fMetPtHist = new TH1D("hMetPtHist",";p_{t};#",30,0.,300.);
438 fMetPhiHist = new TH1D("hMetPhiHist",";#phi;#",28,-3.5,3.5);
439 AddOutput(fMetPtHist);
440 AddOutput(fMetPhiHist);
441
442
443 }
444
445 //--------------------------------------------------------------------------------------------------
446 void ObjectCleaningMod::SlaveTerminate()
447 {
448 // Run finishing code on the computer (slave) that did the analysis. For this
449 // module, we dont do anything here.
450
451 }
452
453 //--------------------------------------------------------------------------------------------------
454 void ObjectCleaningMod::Terminate()
455 {
456 // Run finishing code on the client computer. For this module, we dont do
457 // anything here.
458 }