ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/FakeMods/src/GenFakesMod.cc
Revision: 1.1
Committed: Tue Jun 30 10:47:17 2009 UTC (15 years, 10 months ago) by loizides
Content type: text/plain
Branch: MAIN
Log Message:
Added FakeMods.

File Contents

# User Rev Content
1 loizides 1.1 // $Id: $
2    
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     fElectronFakeableObjectsName(ModNames::gkElFakeableObjsName),
29     fMuonFakeableObjectsName(ModNames::gkMuFakeableObjsName),
30     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     const ElectronCol *ElectronFakeableObjects = 0;
55     if (!fElectronFakeableObjectsName.IsNull())
56     ElectronFakeableObjects = GetObjThisEvt<ElectronCol>(fElectronFakeableObjectsName);
57     const MuonCol *MuonFakeableObjects = 0;
58     if (!fMuonFakeableObjectsName.IsNull())
59     MuonFakeableObjects = GetObjThisEvt<MuonCol>(fMuonFakeableObjectsName);
60    
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     for (UInt_t n = 0; n < MuonFakeableObjects->GetEntries();n++) {
110    
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     Double_t deltaR = MathUtils::DeltaR(MuonFakeableObjects->At(n)->Phi(),
124     MuonFakeableObjects->At(n)->Eta(),
125     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     if (MathUtils::DeltaR(MuonFakeableObjects->At(n)->Phi(), MuonFakeableObjects->At(n)->Eta(),
139     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     MuonFakeableObjects->At(n)->Mom());
159     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     muonFakeProb = fFakeRate->MuonFakeRate(MuonFakeableObjects->At(n)->Et(),
174     MuonFakeableObjects->At(n)->Eta(),
175     MuonFakeableObjects->At(n)->Phi());
176     muonFakeProbLowError = fFakeRate->MuonFakeRateError(MuonFakeableObjects->At(n)->Et(),
177     MuonFakeableObjects->At(n)->Eta(),
178     MuonFakeableObjects->At(n)->Phi());
179     muonFakeProbHighError = fFakeRate->MuonFakeRateError(MuonFakeableObjects->At(n)->Et(),
180     MuonFakeableObjects->At(n)->Eta(),
181     MuonFakeableObjects->At(n)->Phi());
182     } 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     fakeMuonEvent->AddFakeObject(MuonFakeableObjects->At(n),kMuon,isCleanLepton,isGenMuon);
216    
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     for (UInt_t n = 0; n < ElectronFakeableObjects->GetEntries();n++) {
260    
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     Double_t deltaR = MathUtils::DeltaR(ElectronFakeableObjects->At(n)->Phi(),
274     ElectronFakeableObjects->At(n)->Eta(),
275     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     if (MathUtils::DeltaR(ElectronFakeableObjects->At(n)->Phi(),
289     ElectronFakeableObjects->At(n)->Eta(),
290     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     ElectronFakeableObjects->At(n)->Mom());
310     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     ElectronFakeableObjects->At(n)->Mom());
325     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     electronFakeProb = fFakeRate->ElectronFakeRate(ElectronFakeableObjects->At(n)->Et(),
340     ElectronFakeableObjects->At(n)->Eta(),
341     ElectronFakeableObjects->At(n)->Phi());
342     electronFakeProbLowError = fFakeRate->ElectronFakeRateError(
343     ElectronFakeableObjects->At(n)->Et(),
344     ElectronFakeableObjects->At(n)->Eta(),
345     ElectronFakeableObjects->At(n)->Phi());
346     electronFakeProbHighError = fFakeRate->ElectronFakeRateError(
347     ElectronFakeableObjects->At(n)->Et(),
348     ElectronFakeableObjects->At(n)->Eta(),
349     ElectronFakeableObjects->At(n)->Phi());
350     } 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     fakeElectronEvent->AddFakeObject(ElectronFakeableObjects->At(n),
384     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     }