ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/FakeMods/src/GenFakesMod.cc
Revision: 1.2
Committed: Thu Jul 2 12:17:32 2009 UTC (15 years, 10 months ago) by phedex
Content type: text/plain
Branch: MAIN
Changes since 1.1: +41 -41 lines
Log Message:
Shorten ElectronFakeableObjects to ElFakeableObjects, and MuonFakeableObjects to MuFakeableObj, consistently everywhere.

File Contents

# User Rev Content
1 phedex 1.2 // $Id: GenFakesMod.cc,v 1.1 2009/06/30 10:47:17 loizides Exp $
2 loizides 1.1
3     #include "MitPhysics/FakeMods/interface/GenFakesMod.h"
4     #include "MitCommon/MathTools/interface/MathUtils.h"
5     #include "MitAna/DataUtil/interface/Debug.h"
6     #include "MitPhysics/Init/interface/ModNames.h"
7     #include "MitPhysics/FakeMods/interface/FakeEventHeader.h"
8    
9     using namespace mithep;
10    
11     ClassImp(mithep::GenFakesMod)
12    
13     //--------------------------------------------------------------------------------------------------
14     GenFakesMod::GenFakesMod(const char *name, const char *title) :
15     BaseMod(name,title),
16     fElectronFRFilename("InputRequired"),
17     fMuonFRFilename("InputRequired"),
18     fUse2DFakeRate(false),
19     fUseFitFunction(false),
20     fElectronFRFunctionName("InputRequired"),
21     fMuonFRFunctionName("InputRequired"),
22     fElectronFRHistName("InputRequired"),
23     fMuonFRHistName("InputRequired"),
24     fCleanElectronsName(ModNames::gkCleanElectronsName),
25     fCleanMuonsName(ModNames::gkCleanMuonsName),
26     fCleanPhotonsName(ModNames::gkCleanPhotonsName),
27     fCleanJetsName(ModNames::gkCleanJetsName),
28 phedex 1.2 fElFakeableObjsName(ModNames::gkElFakeableObjsName),
29     fMuFakeableObjsName(ModNames::gkMuFakeableObjsName),
30 loizides 1.1 fMCLeptonsName(ModNames::gkMCLeptonsName),
31     fMCTausName(ModNames::gkMCTausName),
32     fFakeEventHeadersName(ModNames::gkFakeEventHeadersName)
33     {
34     // Constructor.
35     }
36    
37     //--------------------------------------------------------------------------------------------------
38     void GenFakesMod::LoadFakeRate()
39     {
40     //Load FakeRate Probabilities.
41     fFakeRate = new FakeRate(fElectronFRFilename,fMuonFRFilename,fElectronFRFunctionName,
42     fMuonFRFunctionName,fElectronFRHistName,
43     fMuonFRHistName,fUse2DFakeRate, fUseFitFunction );
44    
45     }
46    
47    
48     //--------------------------------------------------------------------------------------------------
49     void GenFakesMod::Process()
50     {
51     // Process entries of the tree.
52    
53     // get input Fakeable object collections
54 phedex 1.2 const ElectronCol *ElFakeableObjs = 0;
55     if (!fElFakeableObjsName.IsNull())
56     ElFakeableObjs = GetObjThisEvt<ElectronCol>(fElFakeableObjsName);
57     const MuonCol *MuFakeableObjs = 0;
58     if (!fMuFakeableObjsName.IsNull())
59     MuFakeableObjs = GetObjThisEvt<MuonCol>(fMuFakeableObjsName);
60 loizides 1.1
61     // get input clean object collections
62     const ElectronCol *CleanElectrons = 0;
63     if (!fCleanElectronsName.IsNull())
64     CleanElectrons = GetObjThisEvt<ElectronCol>(fCleanElectronsName);
65     const MuonCol *CleanMuons = 0;
66     if (!fCleanMuonsName.IsNull())
67     CleanMuons = GetObjThisEvt<MuonCol>(fCleanMuonsName);
68     const PhotonCol *CleanPhotons = 0;
69     if (!fCleanPhotonsName.IsNull())
70     CleanPhotons = GetObjThisEvt<PhotonCol>(fCleanPhotonsName);
71     const JetCol *CleanJets = 0;
72     if (!fCleanJetsName.IsNull())
73     CleanJets = GetObjThisEvt<JetCol>(fCleanJetsName);
74     mithep::ParticleOArr *CleanLeptons = dynamic_cast<mithep::ParticleOArr*>
75     (FindObjThisEvt(ModNames::gkMergedLeptonsName));
76    
77     //get monte carlo collections
78     const MCParticleCol *GenLeptons = 0;
79     if (!fMCLeptonsName.IsNull())
80     GenLeptons = GetObjThisEvt<MCParticleCol>(fMCLeptonsName);
81     const MCParticleCol *GenTaus = 0;
82     if (!fMCTausName.IsNull())
83     GenTaus = GetObjThisEvt<MCParticleCol>(fMCTausName);
84     ObjArray<MCParticle> *GenLeptonsAndTaus = new ObjArray<MCParticle>;
85     if (GenLeptons) {
86     for (UInt_t i=0; i<GenLeptons->GetEntries(); i++)
87     GenLeptonsAndTaus->Add(GenLeptons->At(i));
88     }
89     if (GenTaus) {
90     for (UInt_t i=0; i<GenTaus->GetEntries(); i++)
91     GenLeptonsAndTaus->Add(GenTaus->At(i));
92     }
93    
94     // create final output collection
95     ObjArray <FakeEventHeader> *FakeEventHeaders = new ObjArray <FakeEventHeader> ;
96     FakeEventHeaders->SetOwner(kTRUE);
97    
98     //initialize with one fake event containing no fake objects and all jets.
99     FakeEventHeader *initialFakeEvent = new FakeEventHeader();
100     for (UInt_t j=0;j<CleanJets->GetEntries();j++)
101     initialFakeEvent->AddJets(CleanJets->At(j));
102    
103     FakeEventHeaders->AddOwned(initialFakeEvent);
104    
105     // *****************************************************************************************
106     // Fake into Muons
107     // Loop through all Muon Fakeable objects and consider the fake possibility.
108     // *****************************************************************************************
109 phedex 1.2 for (UInt_t n = 0; n < MuFakeableObjs->GetEntries();n++) {
110 loizides 1.1
111     //make temporary fake event headers array
112     ObjArray <FakeEventHeader> *tempFakeEventHeaders = new ObjArray <FakeEventHeader> ;
113     tempFakeEventHeaders->SetOwner(kTRUE);
114    
115     //loop over all fake events generated so far - and perform an additional fake if necessary
116     for (UInt_t i=0; i<FakeEventHeaders->GetEntries();i++) {
117    
118     // *****************************************************************************************
119     // Determine if the fakeable object was a clean lepton
120     // *****************************************************************************************
121     Bool_t isCleanLepton = false;
122     for (UInt_t j = 0; j < CleanLeptons->GetEntries() ; j++) {
123 phedex 1.2 Double_t deltaR = MathUtils::DeltaR(MuFakeableObjs->At(n)->Phi(),
124     MuFakeableObjs->At(n)->Eta(),
125 loizides 1.1 CleanLeptons->At(j)->Phi(), CleanLeptons->At(j)->Eta());
126    
127     if (deltaR < 0.3) {
128     isCleanLepton = true;
129     break;
130     }
131     }
132    
133     // *****************************************************************************************
134     // Determine if the fakeable object was a real muon from Monte Carlo
135     // *****************************************************************************************
136     Bool_t isGenMuon = false;
137     for (UInt_t l=0; l<GenLeptonsAndTaus->GetEntries(); l++) {
138 phedex 1.2 if (MathUtils::DeltaR(MuFakeableObjs->At(n)->Phi(), MuFakeableObjs->At(n)->Eta(),
139 loizides 1.1 GenLeptonsAndTaus->At(l)->Phi(),
140     GenLeptonsAndTaus->At(l)->Eta()) < 0.3) {
141     isGenMuon = true;
142     break;
143     }
144     }
145    
146     //this is used to determine the weight of the unfaked event.
147     double totalCumulativeFakeProbability = 0.0;
148    
149     // *****************************************************************************************
150     // Perform Muon Fake
151     // *****************************************************************************************
152    
153     //match fake to one of the jets
154     int fakeToJetMatch = -1; //index of the jet that matches to the fake
155     double minDR = 5000;
156     for (UInt_t jj=0;jj<FakeEventHeaders->At(i)->NJets();jj++) {
157     Double_t deltaR = MathUtils::DeltaR(FakeEventHeaders->At(i)->UnfakedJet(jj)->Mom(),
158 phedex 1.2 MuFakeableObjs->At(n)->Mom());
159 loizides 1.1 if (deltaR < minDR) {
160     minDR = deltaR;
161     fakeToJetMatch = jj;
162     }
163     }
164     if (!(minDR < 0.5)) {
165     fakeToJetMatch = -1;
166     }
167    
168     //Obtain the muon FakeRate
169     Double_t muonFakeProb = 0.0;
170     Double_t muonFakeProbLowError = 0.0;
171     Double_t muonFakeProbHighError = 0.0;
172     if(fFakeRate) {
173 phedex 1.2 muonFakeProb = fFakeRate->MuonFakeRate(MuFakeableObjs->At(n)->Et(),
174     MuFakeableObjs->At(n)->Eta(),
175     MuFakeableObjs->At(n)->Phi());
176     muonFakeProbLowError = fFakeRate->MuonFakeRateError(MuFakeableObjs->At(n)->Et(),
177     MuFakeableObjs->At(n)->Eta(),
178     MuFakeableObjs->At(n)->Phi());
179     muonFakeProbHighError = fFakeRate->MuonFakeRateError(MuFakeableObjs->At(n)->Et(),
180     MuFakeableObjs->At(n)->Eta(),
181     MuFakeableObjs->At(n)->Phi());
182 loizides 1.1 } else {
183     cerr << "Error: fFakeRate is a NULL pointer.\n";
184     assert(false);
185     }
186    
187     //only fake into a muon if the fakeable object did not match to a clean lepton
188     if (!isCleanLepton) {
189    
190     //create new fake event header
191     FakeEventHeader *fakeMuonEvent = new FakeEventHeader();
192     fakeMuonEvent->SetWeight(FakeEventHeaders->At(i)->Weight() * muonFakeProb);
193     fakeMuonEvent->SetWeightLowError(FakeEventHeaders->At(i)->Weight()*muonFakeProb*
194     TMath::Sqrt((FakeEventHeaders->At(i)->WeightLowError()/
195     FakeEventHeaders->At(i)->Weight())*
196     (FakeEventHeaders->At(i)->WeightLowError()/
197     FakeEventHeaders->At(i)->Weight()) +
198     (muonFakeProbLowError/muonFakeProb)*
199     (muonFakeProbLowError/muonFakeProb)
200     ));
201     fakeMuonEvent->SetWeightHighError(FakeEventHeaders->At(i)->Weight()*muonFakeProb*
202     TMath::Sqrt((FakeEventHeaders->At(i)->WeightHighError()/
203     FakeEventHeaders->At(i)->Weight())*
204     (FakeEventHeaders->At(i)->WeightHighError()/
205     FakeEventHeaders->At(i)->Weight()) +
206     (muonFakeProbHighError/muonFakeProb)*
207     (muonFakeProbHighError/muonFakeProb)
208     ));
209    
210     //add all previous fakes
211     for (UInt_t f=0;f<FakeEventHeaders->At(i)->FakeObjsSize();f++) {
212     fakeMuonEvent->AddFakeObject(FakeEventHeaders->At(i)->FakeObj(f));
213     }
214     //add new fake
215 phedex 1.2 fakeMuonEvent->AddFakeObject(MuFakeableObjs->At(n),kMuon,isCleanLepton,isGenMuon);
216 loizides 1.1
217     //add all previous jets except the one matching to the new fake
218     for (UInt_t jj=0;jj<FakeEventHeaders->At(i)->NJets();jj++) {
219     if (jj != fakeToJetMatch)
220     fakeMuonEvent->AddJets(FakeEventHeaders->At(i)->UnfakedJet(jj));
221     }
222    
223     //add fake event to the temporary fake event header array
224     tempFakeEventHeaders->AddOwned(fakeMuonEvent);
225    
226     //increase cumulative fake probability
227     totalCumulativeFakeProbability += muonFakeProb;
228     }
229    
230     // *****************************************************************************************
231     // Do not fake into Muon
232     // *****************************************************************************************
233     //create new fake event header
234     FakeEventHeader *notFakeEvent = new FakeEventHeader();
235     notFakeEvent->SetWeight(FakeEventHeaders->At(i)->Weight() *
236     (1-totalCumulativeFakeProbability));
237     //add previous fakes
238     for (UInt_t f=0;f<FakeEventHeaders->At(i)->FakeObjsSize();f++) {
239     notFakeEvent->AddFakeObject(FakeEventHeaders->At(i)->FakeObj(f));
240     }
241     //add previous jets
242     for (UInt_t jj=0;jj<FakeEventHeaders->At(i)->NJets();jj++) {
243     notFakeEvent->AddJets(FakeEventHeaders->At(i)->UnfakedJet(jj));
244     }
245     tempFakeEventHeaders->AddOwned(notFakeEvent);
246    
247     } //loop over all current fake event headers
248    
249     //replace current fake event headers with the new temporary ones.
250     delete FakeEventHeaders;
251     FakeEventHeaders = tempFakeEventHeaders;
252     } //loop over all muon fakeable objects
253    
254    
255     // *****************************************************************************************
256     // Fake into Electrons
257     // Loop through all Electron Fakeable objects and consider the fake possibility.
258     // *****************************************************************************************
259 phedex 1.2 for (UInt_t n = 0; n < ElFakeableObjs->GetEntries();n++) {
260 loizides 1.1
261     //make temporary fake event headers array
262     ObjArray <FakeEventHeader> *tempFakeEventHeaders = new ObjArray <FakeEventHeader> ;
263     tempFakeEventHeaders->SetOwner(kTRUE);
264    
265     //loop over all fake events generated so far - and perform an additional fake if necessary
266     for (UInt_t i=0; i<FakeEventHeaders->GetEntries();i++) {
267    
268     // *****************************************************************************************
269     // Determine if the fakeable object was a clean lepton
270     // *****************************************************************************************
271     Bool_t isCleanLepton = false;
272     for (UInt_t j = 0; j < CleanLeptons->GetEntries() ; j++) {
273 phedex 1.2 Double_t deltaR = MathUtils::DeltaR(ElFakeableObjs->At(n)->Phi(),
274     ElFakeableObjs->At(n)->Eta(),
275 loizides 1.1 CleanLeptons->At(j)->Phi(), CleanLeptons->At(j)->Eta());
276    
277     if (deltaR < 0.3) {
278     isCleanLepton = true;
279     break;
280     }
281     }
282    
283     // *****************************************************************************************
284     // Determine if the fakeable object was a real electron from Monte Carlo
285     // *****************************************************************************************
286     Bool_t isGenElectron = false;
287     for (UInt_t l=0; l<GenLeptonsAndTaus->GetEntries(); l++) {
288 phedex 1.2 if (MathUtils::DeltaR(ElFakeableObjs->At(n)->Phi(),
289     ElFakeableObjs->At(n)->Eta(),
290 loizides 1.1 GenLeptonsAndTaus->At(l)->Phi(),
291     GenLeptonsAndTaus->At(l)->Eta()) < 0.3) {
292     isGenElectron = true;
293     }
294     }
295    
296    
297     //this is used to determine the weight of the unfaked event.
298     double totalCumulativeFakeProbability = 0.0;
299    
300     // *****************************************************************************************
301     // Determine if the current electron fakeable object already corresponds to one of the
302     // fake muons already in the FakeEventHeader, determined based on deltaR proximity.
303     // If the current electron fakeable object corresponds to one of the fake muon, then
304     // we do not allow it to fake an electron, since one denominator cannot fake two leptons.
305     // *****************************************************************************************
306     Bool_t alreadyFaked = false;
307     for (UInt_t f = 0; f < FakeEventHeaders->At(i)->FakeObjsSize() ; f++) {
308     double deltaR = MathUtils::DeltaR(FakeEventHeaders->At(i)->FakeObj(f)->Mom(),
309 phedex 1.2 ElFakeableObjs->At(n)->Mom());
310 loizides 1.1 if (deltaR < 0.3)
311     alreadyFaked = true;
312     }
313     if (!alreadyFaked) {
314    
315     // *****************************************************************************************
316     // Perform Electron Fake
317     // *****************************************************************************************
318    
319     //match fake to one of the jets
320     int fakeToJetMatch = -1; //index of the jet that matches to the fake
321     double minDR = 5000;
322     for (UInt_t jj=0;jj<FakeEventHeaders->At(i)->NJets();jj++) {
323     Double_t deltaR = MathUtils::DeltaR(FakeEventHeaders->At(i)->UnfakedJet(jj)->Mom(),
324 phedex 1.2 ElFakeableObjs->At(n)->Mom());
325 loizides 1.1 if (deltaR < minDR) {
326     minDR = deltaR;
327     fakeToJetMatch = jj;
328     }
329     }
330     if (!(minDR < 0.5)) {
331     fakeToJetMatch = -1;
332     }
333    
334     //Obtain the electron FakeRate
335     Double_t electronFakeProb = 0.0;
336     Double_t electronFakeProbLowError = 0.0;
337     Double_t electronFakeProbHighError = 0.0;
338     if(fFakeRate) {
339 phedex 1.2 electronFakeProb = fFakeRate->ElectronFakeRate(ElFakeableObjs->At(n)->Et(),
340     ElFakeableObjs->At(n)->Eta(),
341     ElFakeableObjs->At(n)->Phi());
342 loizides 1.1 electronFakeProbLowError = fFakeRate->ElectronFakeRateError(
343 phedex 1.2 ElFakeableObjs->At(n)->Et(),
344     ElFakeableObjs->At(n)->Eta(),
345     ElFakeableObjs->At(n)->Phi());
346 loizides 1.1 electronFakeProbHighError = fFakeRate->ElectronFakeRateError(
347 phedex 1.2 ElFakeableObjs->At(n)->Et(),
348     ElFakeableObjs->At(n)->Eta(),
349     ElFakeableObjs->At(n)->Phi());
350 loizides 1.1 } else {
351     cerr << "Error: fFakeRate is a NULL pointer.\n";
352     assert(false);
353     }
354    
355     //only fake into a muon if the fakeable object did not match to a clean lepton
356     if (!isCleanLepton) {
357     //create new fake event header
358     FakeEventHeader *fakeElectronEvent = new FakeEventHeader();
359     fakeElectronEvent->SetWeight(FakeEventHeaders->At(i)->Weight() * electronFakeProb);
360     fakeElectronEvent->SetWeightLowError(FakeEventHeaders->At(i)->Weight()*electronFakeProb*
361     TMath::Sqrt((FakeEventHeaders->At(i)->WeightLowError()/
362     FakeEventHeaders->At(i)->Weight())*
363     (FakeEventHeaders->At(i)->WeightLowError()/
364     FakeEventHeaders->At(i)->Weight()) +
365     (electronFakeProbLowError/electronFakeProb)*
366     (electronFakeProbLowError/electronFakeProb)
367     ));
368     fakeElectronEvent->SetWeightHighError(FakeEventHeaders->At(i)->Weight()*electronFakeProb*
369     TMath::Sqrt((FakeEventHeaders->At(i)->WeightHighError()/
370     FakeEventHeaders->At(i)->Weight())*
371     (FakeEventHeaders->At(i)->WeightHighError()/
372     FakeEventHeaders->At(i)->Weight()) +
373     (electronFakeProbHighError/electronFakeProb)*
374     (electronFakeProbHighError/electronFakeProb)
375     ));
376    
377    
378     //add previous fakes
379     for (UInt_t f=0;f<FakeEventHeaders->At(i)->FakeObjsSize();f++) {
380     fakeElectronEvent->AddFakeObject(FakeEventHeaders->At(i)->FakeObj(f));
381     }
382     //add the new fake
383 phedex 1.2 fakeElectronEvent->AddFakeObject(ElFakeableObjs->At(n),
384 loizides 1.1 kElectron,isCleanLepton,isGenElectron);
385     //add previous jets that do not match to the new fake
386     for (UInt_t jj=0;jj<FakeEventHeaders->At(i)->NJets();jj++) {
387     if (jj != fakeToJetMatch)
388     fakeElectronEvent->AddJets(FakeEventHeaders->At(i)->UnfakedJet(jj));
389     }
390    
391     //add fake event to the temporary fake event header array
392     tempFakeEventHeaders->AddOwned(fakeElectronEvent);
393     //increase cumulative fake probability
394     totalCumulativeFakeProbability += electronFakeProb;
395     }
396     }
397    
398     // *****************************************************************************************
399     // Do not fake into anything
400     // *****************************************************************************************
401     //create new fake event header
402     FakeEventHeader *notFakeEvent = new FakeEventHeader();
403     notFakeEvent->SetWeight(FakeEventHeaders->At(i)->Weight() *
404     (1-totalCumulativeFakeProbability));
405     //add previous fakes
406     for (UInt_t f=0;f<FakeEventHeaders->At(i)->FakeObjsSize();f++) {
407     notFakeEvent->AddFakeObject(FakeEventHeaders->At(i)->FakeObj(f));
408     }
409     //add previous jets
410     for (UInt_t jj=0;jj<FakeEventHeaders->At(i)->NJets();jj++) {
411     notFakeEvent->AddJets(FakeEventHeaders->At(i)->UnfakedJet(jj));
412     }
413     tempFakeEventHeaders->AddOwned(notFakeEvent);
414    
415     } //for all current fake event headers
416    
417     //replace current fake event headers with the new temporary ones.
418     delete FakeEventHeaders;
419     FakeEventHeaders = tempFakeEventHeaders;
420     } //loop over all fakeable objects
421    
422     // export FakeEventHeaders for other modules to use
423     FakeEventHeaders->SetName(fFakeEventHeadersName);
424     AddObjThisEvt(FakeEventHeaders);
425    
426     //delete temporary collections
427     delete GenLeptonsAndTaus;
428     }