ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/FakeMods/src/GenFakesMod.cc
Revision: 1.3
Committed: Mon Jul 13 11:27:13 2009 UTC (15 years, 9 months ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_010
Changes since 1.2: +12 -7 lines
Log Message:
Fix compiler warnings... Still this code needs further cleanup.

File Contents

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