ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.43
Committed: Mon Jun 15 15:00:21 2009 UTC (15 years, 10 months ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_009c, Mit_009b
Changes since 1.42: +3 -1 lines
Log Message: 
Added proper fwd defs plus split up complilation of MitAna/DataTree LinkDefs.

File Contents

# User Rev Content
1 loizides 1.43 // $Id: GeneratorMod.cc,v 1.42 2009/06/11 13:16:53 ceballos Exp $
2 loizides 1.1
3     #include "MitPhysics/Mods/interface/GeneratorMod.h"
4     #include "MitCommon/MathTools/interface/MathUtils.h"
5 loizides 1.43 #include "MitAna/DataTree/interface/MetCol.h"
6     #include "MitAna/DataTree/interface/MCParticleCol.h"
7 loizides 1.8 #include "MitPhysics/Init/interface/ModNames.h"
8 loizides 1.1 #include <TH1D.h>
9     #include <TH2D.h>
10    
11     using namespace mithep;
12    
13     ClassImp(mithep::GeneratorMod)
14    
15     //--------------------------------------------------------------------------------------------------
16 loizides 1.8 GeneratorMod::GeneratorMod(const char *name, const char *title) :
17 loizides 1.1 BaseMod(name,title),
18 ceballos 1.32 fPrintDebug(kFALSE),
19 loizides 1.1 fMCPartName(Names::gkMCPartBrn),
20 ceballos 1.30 fMCMETName(ModNames::gkMCMETName),
21 loizides 1.8 fMCLeptonsName(ModNames::gkMCLeptonsName),
22     fMCAllLeptonsName(ModNames::gkMCAllLeptonsName),
23     fMCTausName(ModNames::gkMCTausName),
24     fMCNeutrinosName(ModNames::gkMCNeutrinosName),
25     fMCQuarksName(ModNames::gkMCQuarksName),
26     fMCqqHsName(ModNames::gkMCqqHsName),
27     fMCBosonsName(ModNames::gkMCBosonsName),
28 ceballos 1.10 fMCPhotonsName(ModNames::gkMCPhotonsName),
29 ceballos 1.32 fMCRadPhotonsName(ModNames::gkMCRadPhotonsName),
30     fMCISRPhotonsName(ModNames::gkMCISRPhotonsName),
31 ceballos 1.14 fPtLeptonMin(0.0),
32     fEtaLeptonMax(5.0),
33     fPtPhotonMin(0.0),
34 loizides 1.15 fEtaPhotonMax(5.0),
35 ceballos 1.32 fPtRadPhotonMin(0.0),
36     fEtaRadPhotonMax(5.0),
37 loizides 1.27 fPdgIdCut(0),
38     fMassMinCut(-FLT_MAX),
39     fMassMaxCut(FLT_MAX),
40 ceballos 1.34 fApplyISRFilter(kFALSE),
41 loizides 1.15 fParticles(0)
42 loizides 1.1 {
43 ceballos 1.42 // Constructor
44 loizides 1.1 }
45    
46     //--------------------------------------------------------------------------------------------------
47     void GeneratorMod::Process()
48     {
49 loizides 1.5 // Process entries of the tree.
50 loizides 1.1
51 loizides 1.5 // these arrays will be filled in the loop of particles
52 ceballos 1.30 MetOArr *GenMet = new MetOArr;
53     GenMet->SetName(fMCMETName);
54     GenMet->SetOwner(kTRUE);
55 loizides 1.9 MCParticleOArr *GenLeptons = new MCParticleOArr;
56 loizides 1.13 GenLeptons->SetName(fMCLeptonsName);
57 loizides 1.9 MCParticleOArr *GenAllLeptons = new MCParticleOArr;
58 loizides 1.13 GenAllLeptons->SetName(fMCAllLeptonsName);
59 loizides 1.9 MCParticleOArr *GenTaus = new MCParticleOArr;
60 loizides 1.13 GenTaus->SetName(fMCTausName);
61     GenTaus->SetOwner(kTRUE);
62 loizides 1.9 MCParticleOArr *GenNeutrinos = new MCParticleOArr;
63 loizides 1.13 GenNeutrinos->SetName(fMCNeutrinosName);
64 loizides 1.9 MCParticleOArr *GenQuarks = new MCParticleOArr;
65 loizides 1.13 GenQuarks->SetName(fMCQuarksName);
66 loizides 1.9 MCParticleOArr *GenqqHs = new MCParticleOArr;
67 loizides 1.13 GenqqHs->SetName(fMCqqHsName);
68 loizides 1.9 MCParticleOArr *GenBosons = new MCParticleOArr;
69 loizides 1.13 GenBosons->SetName(fMCBosonsName);
70 ceballos 1.10 MCParticleOArr *GenPhotons = new MCParticleOArr;
71 loizides 1.13 GenPhotons->SetName(fMCPhotonsName);
72 ceballos 1.32 MCParticleOArr *GenRadPhotons = new MCParticleOArr;
73     GenRadPhotons->SetName(fMCRadPhotonsName);
74     MCParticleOArr *GenISRPhotons = new MCParticleOArr;
75     GenISRPhotons->SetName(fMCISRPhotonsName);
76    
77 ceballos 1.42 MCParticleOArr *GenTempMG0 = new MCParticleOArr;
78    
79 loizides 1.40 if(fPrintDebug)
80     printf("\n************ Next Event ************\n\n");
81 loizides 1.1
82 loizides 1.5 // load MCParticle branch
83 loizides 1.41 LoadEventObject(fMCPartName, fParticles);
84 loizides 1.5
85 ceballos 1.42 Bool_t isOld = kFALSE;
86     Int_t sumV[2] = {0, 0}; // W, Z
87     Int_t sumVVFlavor[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
88 ceballos 1.30 Double_t totalMET[3] = {0.0, 0.0, 0.0};
89 loizides 1.6 Bool_t isqqH = kFALSE;
90 loizides 1.5 for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
91 loizides 1.7 const MCParticle *p = fParticles->At(i);
92 ceballos 1.32
93 loizides 1.40 if(fPrintDebug)
94     p->Print("l");
95 ceballos 1.32
96 loizides 1.40 // rad photons
97 ceballos 1.34 if(p->Is(MCParticle::kGamma) && p->HasMother() && p->Mother()->Status() == 3 &&
98 ceballos 1.33 (p->Mother()->Is(MCParticle::kEl) || p->Mother()->Is(MCParticle::kMu) ||
99 ceballos 1.34 p->Mother()->Is(MCParticle::kTau)) &&
100 ceballos 1.32 p->Pt() > fPtRadPhotonMin && p->AbsEta() < fEtaRadPhotonMax) {
101     GenRadPhotons->Add(p);
102     }
103    
104     // ISR photons
105     if(p->Is(MCParticle::kGamma) && p->HasMother() && p->Mother()->IsQuark()) {
106     GenISRPhotons->Add(p);
107     }
108    
109 ceballos 1.30 // MET computation at generation level
110     if (p->Status() == 1 && !p->IsNeutrino()) {
111     totalMET[0] = totalMET[0] + p->Px();
112     totalMET[1] = totalMET[1] + p->Py();
113     totalMET[2] = totalMET[2] + p->Pz();
114     }
115 loizides 1.5
116     if (!p->IsGenerated()) continue;
117    
118     // muons/electrons from W/Z decays
119 loizides 1.7 if ((p->Is(MCParticle::kEl) || p->Is(MCParticle::kMu)) && p->Status() == 1) {
120 ceballos 1.14 if (p->Pt() > fPtLeptonMin && p->AbsEta() < fEtaLeptonMax) {
121 loizides 1.5 GenAllLeptons->Add(p);
122     }
123 loizides 1.6 Bool_t isGoodLepton = kFALSE;
124 loizides 1.7 const MCParticle *pm = p;
125 loizides 1.5 while (pm->HasMother() && isGoodLepton == kFALSE) {
126 loizides 1.21 if (pm->PdgId() == 92) // string reached, terminate loop
127     break;
128 loizides 1.29 if (pm->Mother()->Is(MCParticle::kZ) || pm->Mother()->Is(MCParticle::kW) ||
129     pm->Mother()->Is(MCParticle::kZp) || pm->Mother()->Is(MCParticle::kWp) ||
130     pm->Mother()->Is(MCParticle::kH)) {
131 loizides 1.5 GenLeptons->Add(p);
132     isGoodLepton = kTRUE;
133 loizides 1.7 break;
134 loizides 1.8 } else if (pm->Mother()->Is(MCParticle::kPi0) || pm->Mother()->Is(MCParticle::kEta)) {
135 loizides 1.7 // this is fake, but it is a trick to get rid of these cases and abort the loop
136     break;
137     }
138     pm = pm->Mother();
139 loizides 1.1 }
140 loizides 1.5 }
141    
142 loizides 1.7 // hadronic taus
143     else if (p->Is(MCParticle::kTau) && p->Status() == 2) {
144     if (!p->HasDaughter(MCParticle::kEl) && !p->HasDaughter(MCParticle::kMu)) {
145     const MCParticle *tv = p->FindDaughter(MCParticle::kTauNu);
146     if (tv) {
147     MCParticle *pm_f = new MCParticle(*p);
148     pm_f->SetMom(p->Px()-tv->Px(), p->Py()-tv->Py(),
149     p->Pz()-tv->Pz(), p->E()-tv->E());
150     GenTaus->AddOwned(pm_f);
151     } else {
152 loizides 1.18 SendError(kWarning, "Process", "Could not find a tau neutrino!");
153 loizides 1.5 }
154 loizides 1.1 }
155 loizides 1.5 }
156 loizides 1.1
157 loizides 1.5 // neutrinos
158 loizides 1.7 else if (p->Status() == 1 && p->IsNeutrino()) {
159 loizides 1.5 GenNeutrinos->Add(p);
160     }
161 loizides 1.1
162 loizides 1.5 // quarks from W/Z decays or top particles
163 loizides 1.7 else if (p->IsQuark() && p->HasMother()) {
164     if (p->Mother()->Is(MCParticle::kZ) || p->Mother()->Is(MCParticle::kW) ||
165     p->Is(MCParticle::kTop) || p->Mother()->Is(MCParticle::kTop)) {
166 loizides 1.5 GenQuarks->Add(p);
167 loizides 1.1 }
168 loizides 1.5 }
169 loizides 1.1
170 loizides 1.5 // qqH, information about the forward jets
171 loizides 1.8 else if (isqqH == kFALSE && p->Is(MCParticle::kH)) {
172 loizides 1.5 isqqH = kTRUE;
173 loizides 1.17 const MCParticle *pq1 = fParticles->At(i-1);
174     const MCParticle *pq2 = fParticles->At(i-2);
175 loizides 1.7 if (!pq1 || !pq2) {
176     SendError(kWarning, "Process", "Could not find quark pair!");
177 loizides 1.8 } else if (pq1->IsQuark() && pq2->IsQuark() &&
178     pq1->HasMother() && pq2->HasMother() &&
179 ceballos 1.36 pq1->Mother() == pq2->Mother()) {
180 loizides 1.5 GenqqHs->Add(pq1);
181     GenqqHs->Add(pq2);
182 loizides 1.1 }
183 ceballos 1.23
184 loizides 1.7 if (p->Status() == 3)
185 loizides 1.24 GenBosons->Add(p); // take higgs boson in account here rather in next else if
186 loizides 1.5 }
187 loizides 1.1
188 loizides 1.5 // information about bosons: W, Z, h, Z', W', H0, A0, H+
189 loizides 1.7 else if (p->Status() == 3 &&
190     (p->Is(MCParticle::kZ) || p->Is(MCParticle::kW) || p->Is(MCParticle::kH) ||
191     p->Is(MCParticle::kZp) || p->Is(MCParticle::kZpp) ||
192     p->Is(MCParticle::kH0) || p->Is(MCParticle::kA0) || p->Is(MCParticle::kHp))) {
193 loizides 1.5 GenBosons->Add(p);
194 ceballos 1.42 if (p->Is(MCParticle::kW)) sumV[0]++;
195     else if(p->Is(MCParticle::kZ)) sumV[1]++;
196     if (p->Is(MCParticle::kW) && p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu))
197     sumVVFlavor[0]++;
198     else if(p->Is(MCParticle::kW) && p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu))
199     sumVVFlavor[1]++;
200     else if(p->Is(MCParticle::kW) && p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu))
201     sumVVFlavor[2]++;
202    
203     else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kMu,kTRUE) && p->HasDaughter(-1*MCParticle::kMu,kTRUE))
204     sumVVFlavor[3]++;
205     else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kEl,kTRUE) && p->HasDaughter(-1*MCParticle::kEl,kTRUE))
206     sumVVFlavor[4]++;
207     else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kTau,kTRUE) && p->HasDaughter(-1*MCParticle::kTau,kTRUE))
208     sumVVFlavor[5]++;
209    
210     else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kMuNu,kTRUE) && p->HasDaughter(-1*MCParticle::kMuNu,kTRUE))
211     sumVVFlavor[6]++;
212     else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kElNu,kTRUE) && p->HasDaughter(-1*MCParticle::kElNu,kTRUE))
213     sumVVFlavor[7]++;
214     else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kTauNu,kTRUE) && p->HasDaughter(-1*MCParticle::kTauNu,kTRUE))
215     sumVVFlavor[8]++;
216 loizides 1.1 }
217 ceballos 1.10
218 ceballos 1.14 // photons
219 ceballos 1.10 else if (p->Status() == 1 && p->Is(MCParticle::kGamma) &&
220 ceballos 1.14 p->Pt() > fPtPhotonMin && p->AbsEta() < fEtaPhotonMax) {
221 ceballos 1.10 GenPhotons->Add(p);
222     }
223    
224 loizides 1.18 // W/Z -> lnu for Madgraph
225 loizides 1.28 if (p->IsParton() && p->NDaughters() >= 2) {
226 ceballos 1.16 CompositeParticle *diBoson = new CompositeParticle();
227 loizides 1.18 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
228 ceballos 1.42 isOld = kFALSE;
229     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
230     if(p->FindDaughter(MCParticle::kMu) == GenTempMG0->At(nl)) {
231     isOld = kTRUE;
232     break;
233     }
234     }
235     if(isOld == kFALSE){
236     GenTempMG0->Add(p->FindDaughter(MCParticle::kMu));
237     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
238     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
239     sumV[0]++;
240     sumVVFlavor[0]++;
241     if (GetFillHist())
242     hDVMass[0]->Fill(TMath::Min(diBoson->Mass(),199.999));
243     const MCParticle *tmp_mu = p->FindDaughter(MCParticle::kMu);
244     while (tmp_mu->HasDaughter(MCParticle::kMu) &&
245     tmp_mu->FindDaughter(MCParticle::kMu)->IsGenerated())
246     tmp_mu = tmp_mu->FindDaughter(MCParticle::kMu);
247    
248     GenLeptons->Add(tmp_mu);
249     }
250     }
251     if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
252     isOld = kFALSE;
253     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
254     if(p->FindDaughter(MCParticle::kEl) == GenTempMG0->At(nl)) {
255     isOld = kTRUE;
256     break;
257     }
258     }
259     if(isOld == kFALSE){
260     GenTempMG0->Add(p->FindDaughter(MCParticle::kEl));
261     diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
262     diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
263     sumV[0]++;
264     sumVVFlavor[1]++;
265     if (GetFillHist())
266     hDVMass[1]->Fill(TMath::Min(diBoson->Mass(),199.999));
267     const MCParticle *tmp_e = p->FindDaughter(MCParticle::kEl);
268     while (tmp_e->HasDaughter(MCParticle::kEl) &&
269     tmp_e->FindDaughter(MCParticle::kEl)->IsGenerated())
270     tmp_e = tmp_e->FindDaughter(MCParticle::kEl);
271     GenLeptons->Add(tmp_e);
272     }
273     }
274     if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
275     isOld = kFALSE;
276     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
277     if(p->FindDaughter(MCParticle::kTau) == GenTempMG0->At(nl)) {
278     isOld = kTRUE;
279     break;
280     }
281     }
282     if(isOld == kFALSE){
283     GenTempMG0->Add(p->FindDaughter(MCParticle::kTau));
284     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
285     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
286     sumV[0]++;
287     sumVVFlavor[2]++;
288     if (GetFillHist())
289     hDVMass[2]->Fill(TMath::Min(diBoson->Mass(),199.999));
290     const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
291     if (tau->HasDaughter(MCParticle::kMu))
292     GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
293     if (tau->HasDaughter(MCParticle::kEl))
294     GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
295     if (tau->HasDaughter(MCParticle::kTau)) {
296     const MCParticle *tau_second = tau->FindDaughter(MCParticle::kTau);
297     if (tau_second->HasDaughter(MCParticle::kMu))
298     GenLeptons->Add(tau_second->FindDaughter(MCParticle::kMu));
299     if (tau_second->HasDaughter(MCParticle::kEl))
300     GenLeptons->Add(tau_second->FindDaughter(MCParticle::kEl));
301     }
302     }
303     }
304     if (p->HasDaughter(MCParticle::kMu,kTRUE) && p->HasDaughter(-1*MCParticle::kMu,kTRUE)) {
305     isOld = kFALSE;
306     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
307     if(p->FindDaughter(MCParticle::kMu,kTRUE) == GenTempMG0->At(nl)) {
308     isOld = kTRUE;
309     break;
310     }
311     }
312     if(isOld == kFALSE){
313     GenTempMG0->Add(p->FindDaughter(MCParticle::kMu,kTRUE));
314     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu,kTRUE));
315     diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kMu,kTRUE));
316     sumV[1]++;
317     sumVVFlavor[3]++;
318     if (GetFillHist())
319     hDVMass[3]->Fill(TMath::Min(diBoson->Mass(),199.999));
320     const MCParticle *tmp_mu0 = p->FindDaughter(MCParticle::kMu,kTRUE);
321     while (tmp_mu0->HasDaughter(MCParticle::kMu) &&
322     tmp_mu0->FindDaughter(MCParticle::kMu)->IsGenerated())
323     tmp_mu0 = tmp_mu0->FindDaughter(MCParticle::kMu);
324     const MCParticle *tmp_mu1 = p->FindDaughter(-1*MCParticle::kMu,kTRUE);
325     while (tmp_mu1->HasDaughter(MCParticle::kMu) &&
326     tmp_mu1->FindDaughter(MCParticle::kMu)->IsGenerated())
327     tmp_mu1 = tmp_mu1->FindDaughter(MCParticle::kMu);
328     GenLeptons->Add(tmp_mu0);
329     GenLeptons->Add(tmp_mu1);
330     }
331     }
332     if (p->HasDaughter(MCParticle::kEl,kTRUE) && p->HasDaughter(-1*MCParticle::kEl,kTRUE)) {
333     isOld = kFALSE;
334     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
335     if(p->FindDaughter(MCParticle::kEl,kTRUE) == GenTempMG0->At(nl)) {
336     isOld = kTRUE;
337     break;
338     }
339     }
340     if(isOld == kFALSE){
341     GenTempMG0->Add(p->FindDaughter(MCParticle::kEl,kTRUE));
342     diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl,kTRUE));
343     diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kEl,kTRUE));
344     sumV[1]++;
345     sumVVFlavor[4]++;
346     if (GetFillHist())
347     hDVMass[4]->Fill(TMath::Min(diBoson->Mass(),199.999));
348     const MCParticle *tmp_e0 = p->Daughter(0);
349     while (tmp_e0->HasDaughter(MCParticle::kEl) &&
350     tmp_e0->FindDaughter(MCParticle::kEl)->IsGenerated())
351     tmp_e0 = tmp_e0->FindDaughter(MCParticle::kEl);
352     const MCParticle *tmp_e1 = p->Daughter(1);
353     while (tmp_e1->HasDaughter(MCParticle::kEl) &&
354     tmp_e1->FindDaughter(MCParticle::kEl)->IsGenerated())
355     tmp_e1 = tmp_e1->FindDaughter(MCParticle::kEl);
356     GenLeptons->Add(tmp_e0);
357     GenLeptons->Add(tmp_e1);
358     }
359     }
360     if (p->HasDaughter(MCParticle::kTau,kTRUE) && p->HasDaughter(-1*MCParticle::kTau,kTRUE)) {
361     isOld = kFALSE;
362     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
363     if(p->FindDaughter(MCParticle::kTau,kTRUE) == GenTempMG0->At(nl)) {
364     isOld = kTRUE;
365     break;
366     }
367     }
368     if(isOld == kFALSE){
369     GenTempMG0->Add(p->FindDaughter(MCParticle::kTau,kTRUE));
370     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau,kTRUE));
371     diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kTau,kTRUE));
372     sumV[1]++;
373     sumVVFlavor[5]++;
374     if (GetFillHist())
375     hDVMass[5]->Fill(TMath::Min(diBoson->Mass(),199.999));
376     const MCParticle *tau0 = p->Daughter(0);
377     if (tau0->HasDaughter(MCParticle::kMu))
378     GenLeptons->Add(tau0->FindDaughter(MCParticle::kMu));
379     if (tau0->HasDaughter(MCParticle::kEl))
380     GenLeptons->Add(tau0->FindDaughter(MCParticle::kEl));
381     const MCParticle *tau1 = p->Daughter(1);
382     if (tau1->HasDaughter(MCParticle::kMu))
383     GenLeptons->Add(tau1->FindDaughter(MCParticle::kMu));
384     if (tau1->HasDaughter(MCParticle::kEl))
385     GenLeptons->Add(tau1->FindDaughter(MCParticle::kEl));
386     if (tau0->HasDaughter(MCParticle::kTau)) {
387     const MCParticle *tau0_second = tau0->FindDaughter(MCParticle::kTau);
388     if (tau0_second->HasDaughter(MCParticle::kMu))
389     GenLeptons->Add(tau0_second->FindDaughter(MCParticle::kMu));
390     if (tau0_second->HasDaughter(MCParticle::kEl))
391     GenLeptons->Add(tau0_second->FindDaughter(MCParticle::kEl));
392     }
393     if (tau1->HasDaughter(MCParticle::kTau)) {
394     const MCParticle *tau1_second = tau1->FindDaughter(MCParticle::kTau);
395     if (tau1_second->HasDaughter(MCParticle::kMu))
396     GenLeptons->Add(tau1_second->FindDaughter(MCParticle::kMu));
397     if (tau1_second->HasDaughter(MCParticle::kEl))
398     GenLeptons->Add(tau1_second->FindDaughter(MCParticle::kEl));
399     }
400     }
401 ceballos 1.16 }
402 ceballos 1.42 if (p->HasDaughter(MCParticle::kMuNu,kTRUE) && p->HasDaughter(-1*MCParticle::kMuNu,kTRUE)) {
403     isOld = kFALSE;
404     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
405     if(p->FindDaughter(MCParticle::kMuNu,kTRUE) == GenTempMG0->At(nl)) {
406     isOld = kTRUE;
407     break;
408     }
409     }
410     if(isOld == kFALSE){
411     GenTempMG0->Add(p->FindDaughter(MCParticle::kMuNu,kTRUE));
412     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu,kTRUE));
413     diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kMuNu,kTRUE));
414     sumV[1]++;
415     sumVVFlavor[6]++;
416     if (GetFillHist())
417     hDVMass[6]->Fill(TMath::Min(diBoson->Mass(),199.999));
418     }
419 ceballos 1.16 }
420 ceballos 1.42 if (p->HasDaughter(MCParticle::kElNu,kTRUE) && p->HasDaughter(-1*MCParticle::kElNu,kTRUE)) {
421     isOld = kFALSE;
422     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
423     if(p->FindDaughter(MCParticle::kElNu,kTRUE) == GenTempMG0->At(nl)) {
424     isOld = kTRUE;
425     break;
426     }
427     }
428     if(isOld == kFALSE){
429     GenTempMG0->Add(p->FindDaughter(MCParticle::kElNu,kTRUE));
430     diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu,kTRUE));
431     diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kElNu,kTRUE));
432     sumV[1]++;
433     sumVVFlavor[7]++;
434     if (GetFillHist())
435     hDVMass[7]->Fill(TMath::Min(diBoson->Mass(),199.999));
436     }
437 loizides 1.18 }
438 ceballos 1.42 if (p->HasDaughter(MCParticle::kTauNu,kTRUE) && p->HasDaughter(-1*MCParticle::kTauNu,kTRUE)) {
439     isOld = kFALSE;
440     for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
441     if(p->FindDaughter(MCParticle::kTauNu,kTRUE) == GenTempMG0->At(nl)) {
442     isOld = kTRUE;
443     break;
444     }
445     }
446     if(isOld == kFALSE){
447     GenTempMG0->Add(p->FindDaughter(MCParticle::kTauNu,kTRUE));
448     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu,kTRUE));
449     diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kTauNu,kTRUE));
450     sumV[1]++;
451     sumVVFlavor[8]++;
452     if (GetFillHist())
453     hDVMass[8]->Fill(TMath::Min(diBoson->Mass(),199.999));
454     }
455 ceballos 1.16 }
456     delete diBoson;
457     }
458    
459 loizides 1.18 // t -> lnu for Madgraph
460     if (p->Is(MCParticle::kTop)) {
461 ceballos 1.16 CompositeParticle *diBoson = new CompositeParticle();
462 loizides 1.18 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
463 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
464     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
465 loizides 1.20 if (GetFillHist())
466 ceballos 1.42 hDVMass[9]->Fill(TMath::Min(diBoson->Mass(),199.999));
467 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kMu));
468     }
469 loizides 1.18 else if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
470 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
471     diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
472 loizides 1.20 if (GetFillHist())
473 ceballos 1.42 hDVMass[10]->Fill(TMath::Min(diBoson->Mass(),199.999));
474 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kEl));
475     }
476 loizides 1.18 else if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
477 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
478     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
479 loizides 1.20 if (GetFillHist())
480 ceballos 1.42 hDVMass[11]->Fill(TMath::Min(diBoson->Mass(),199.999));
481 loizides 1.18 const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
482     if (tau->HasDaughter(MCParticle::kMu))
483     GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
484     if (tau->HasDaughter(MCParticle::kEl))
485     GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
486     }
487     else if (!p->HasDaughter(MCParticle::kW)) {
488     for(UInt_t nd=0; nd<p->NDaughters(); ++nd)
489     if (p->Daughter(nd)->IsNot(MCParticle::kBottom) &&
490     p->Daughter(nd)->IsNot(MCParticle::kGamma))
491     diBoson->AddDaughter(p->Daughter(nd));
492 loizides 1.20 if (GetFillHist())
493 ceballos 1.42 hDVMass[12]->Fill(TMath::Min(diBoson->Mass(),199.999));
494 ceballos 1.16 }
495     delete diBoson;
496     }
497 ceballos 1.26
498 loizides 1.27 // mass cut for given pid
499     if(fPdgIdCut && p->Is(fPdgIdCut) &&
500 ceballos 1.26 (p->Mass() < fMassMinCut || p->Mass() > fMassMaxCut)) {
501     SkipEvent();
502     return;
503 sixie 1.35 }
504 loizides 1.27 } // end loop of particles
505 loizides 1.1
506 ceballos 1.42 delete GenTempMG0;
507    
508 ceballos 1.30 Met *theMET = new Met(totalMET[0], totalMET[1]);
509     theMET->SetElongitudinal(totalMET[2]);
510     GenMet->AddOwned(theMET);
511    
512 loizides 1.19 // sort according to pt
513     GenLeptons->Sort();
514     GenAllLeptons->Sort();
515     GenTaus->Sort();
516     GenNeutrinos->Sort();
517     GenQuarks->Sort();
518     GenqqHs->Sort();
519     GenBosons->Sort();
520     GenPhotons->Sort();
521 ceballos 1.32 GenRadPhotons->Sort();
522     GenISRPhotons->Sort();
523 loizides 1.19
524 loizides 1.8 // add objects to this event for other modules to use
525 ceballos 1.30 AddObjThisEvt(GenMet);
526 loizides 1.13 AddObjThisEvt(GenLeptons);
527     AddObjThisEvt(GenAllLeptons);
528     AddObjThisEvt(GenTaus);
529     AddObjThisEvt(GenNeutrinos);
530     AddObjThisEvt(GenQuarks);
531     AddObjThisEvt(GenqqHs);
532     AddObjThisEvt(GenBosons);
533     AddObjThisEvt(GenPhotons);
534 ceballos 1.32 AddObjThisEvt(GenRadPhotons);
535     AddObjThisEvt(GenISRPhotons);
536    
537 ceballos 1.42 //if(GenqqHs->GetEntries() == 0) {
538     // SkipEvent();
539     // return;
540     //}
541    
542 loizides 1.5 // fill histograms if requested
543 loizides 1.20 if (GetFillHist()) {
544 loizides 1.5
545 ceballos 1.30 // MET
546     hDGenMet[0]->Fill(GenMet->At(0)->Pt());
547     hDGenMet[1]->Fill(GenMet->At(0)->Px());
548     hDGenMet[2]->Fill(GenMet->At(0)->Py());
549     hDGenMet[3]->Fill(GenMet->At(0)->Elongitudinal());
550    
551 loizides 1.6 // leptons
552 loizides 1.1 hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
553 loizides 1.5 for(UInt_t i=0; i<GenLeptons->GetEntries(); i++) {
554 loizides 1.1 hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
555 loizides 1.5 hDGenLeptons[2]->Fill(TMath::Abs(GenLeptons->At(i)->Eta()));
556     hDGenLeptons[3]->Fill(GenLeptons->At(i)->PhiDeg());
557     for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++) {
558 loizides 1.1 CompositeParticle *dilepton = new CompositeParticle();
559     dilepton->AddDaughter(GenLeptons->At(i));
560     dilepton->AddDaughter(GenLeptons->At(j));
561     hDGenLeptons[4]->Fill(dilepton->Mass());
562     delete dilepton;
563     }
564 ceballos 1.22 }
565     // looking at events with two leptons
566     if (GenLeptons->GetEntries() == 2) {
567     hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(TMath::Abs(GenLeptons->At(0)->Eta()),
568     TMath::Abs(GenLeptons->At(1)->Eta())),
569 loizides 1.5 4.999));
570 ceballos 1.22 hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(TMath::Abs(GenLeptons->At(0)->Eta()),
571     TMath::Abs(GenLeptons->At(1)->Eta())),
572 loizides 1.5 4.999));
573 ceballos 1.22 if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
574     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5) {
575     hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
576     if (GenLeptons->At(0)->Pt() > 20.0) {
577     hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
578     if (GenLeptons->At(1)->Pt() > 10.0) {
579 loizides 1.1 CompositeParticle *dilepton = new CompositeParticle();
580 ceballos 1.22 dilepton->AddDaughter(GenLeptons->At(0));
581     dilepton->AddDaughter(GenLeptons->At(1));
582 loizides 1.1 hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
583 ceballos 1.22 if(dilepton->Mass() > 12.0){
584     hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(0)->Phi(),
585     GenLeptons->At(1)->Phi())
586     * 180./ TMath::Pi());
587 loizides 1.40 hDGenLeptons[11]->Fill(MathUtils::DeltaR(*GenLeptons->At(0),
588     *GenLeptons->At(1)));
589 ceballos 1.22 }
590 loizides 1.1 delete dilepton;
591     }
592     }
593     }
594     }
595 ceballos 1.22 // looking at events with three leptons
596     if (GenLeptons->GetEntries() == 3) {
597     if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
598     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5 &&
599     TMath::Abs(GenLeptons->At(2)->Eta()) < 2.5) {
600     hDGenLeptons[12]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
601     if (GenLeptons->At(0)->Pt() > 20.0) {
602     hDGenLeptons[13]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
603     hDGenLeptons[14]->Fill(TMath::Min(GenLeptons->At(2)->Pt(),199.999));
604     if (GenLeptons->At(1)->Pt() > 10.0 && GenLeptons->At(2)->Pt() > 10.0) {
605     CompositeParticle *dilepton01 = new CompositeParticle();
606     dilepton01->AddDaughter(GenLeptons->At(0));
607     dilepton01->AddDaughter(GenLeptons->At(1));
608     CompositeParticle *dilepton02 = new CompositeParticle();
609     dilepton02->AddDaughter(GenLeptons->At(0));
610     dilepton02->AddDaughter(GenLeptons->At(2));
611     CompositeParticle *dilepton12 = new CompositeParticle();
612     dilepton12->AddDaughter(GenLeptons->At(1));
613     dilepton12->AddDaughter(GenLeptons->At(2));
614     hDGenLeptons[15]->Fill(TMath::Min(dilepton01->Mass(),999.999));
615     hDGenLeptons[15]->Fill(TMath::Min(dilepton02->Mass(),999.999));
616     hDGenLeptons[15]->Fill(TMath::Min(dilepton12->Mass(),999.999));
617     CompositeParticle *trilepton = new CompositeParticle();
618     trilepton->AddDaughter(GenLeptons->At(0));
619     trilepton->AddDaughter(GenLeptons->At(1));
620     trilepton->AddDaughter(GenLeptons->At(2));
621     hDGenLeptons[16]->Fill(TMath::Min(trilepton->Mass(),999.999));
622 loizides 1.40 Double_t deltaR[3] = {MathUtils::DeltaR(*GenLeptons->At(0),
623     *GenLeptons->At(1)),
624     MathUtils::DeltaR(*GenLeptons->At(0),
625     *GenLeptons->At(2)),
626     MathUtils::DeltaR(*GenLeptons->At(1),
627     *GenLeptons->At(2))};
628 loizides 1.27 Double_t deltaRMin = deltaR[0];
629 loizides 1.40 for(Int_t i=1; i<3; i++)
630     if(deltaRMin > deltaR[i])
631     deltaRMin = deltaR[i];
632 ceballos 1.22 hDGenLeptons[17]->Fill(deltaRMin);
633    
634     delete dilepton01;
635     delete dilepton02;
636     delete dilepton12;
637     delete trilepton;
638     }
639     }
640     }
641     }
642     // looking at events with four leptons
643     if (GenLeptons->GetEntries() == 4) {
644     if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
645     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5 &&
646     TMath::Abs(GenLeptons->At(2)->Eta()) < 2.5 &&
647     TMath::Abs(GenLeptons->At(3)->Eta()) < 2.5) {
648     hDGenLeptons[18]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
649     if (GenLeptons->At(0)->Pt() > 20.0) {
650     hDGenLeptons[19]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
651     hDGenLeptons[20]->Fill(TMath::Min(GenLeptons->At(2)->Pt(),199.999));
652     hDGenLeptons[21]->Fill(TMath::Min(GenLeptons->At(3)->Pt(),199.999));
653     if (GenLeptons->At(1)->Pt() > 10.0 && GenLeptons->At(2)->Pt() > 10.0 &&
654     GenLeptons->At(3)->Pt() > 10.0) {
655     CompositeParticle *dilepton01 = new CompositeParticle();
656     dilepton01->AddDaughter(GenLeptons->At(0));
657     dilepton01->AddDaughter(GenLeptons->At(1));
658     CompositeParticle *dilepton02 = new CompositeParticle();
659     dilepton02->AddDaughter(GenLeptons->At(0));
660     dilepton02->AddDaughter(GenLeptons->At(2));
661     CompositeParticle *dilepton03 = new CompositeParticle();
662     dilepton03->AddDaughter(GenLeptons->At(0));
663     dilepton03->AddDaughter(GenLeptons->At(3));
664     CompositeParticle *dilepton12 = new CompositeParticle();
665     dilepton12->AddDaughter(GenLeptons->At(1));
666     dilepton12->AddDaughter(GenLeptons->At(2));
667     CompositeParticle *dilepton13 = new CompositeParticle();
668     dilepton13->AddDaughter(GenLeptons->At(1));
669     dilepton13->AddDaughter(GenLeptons->At(3));
670     CompositeParticle *dilepton23 = new CompositeParticle();
671     dilepton23->AddDaughter(GenLeptons->At(2));
672     dilepton23->AddDaughter(GenLeptons->At(3));
673     hDGenLeptons[22]->Fill(TMath::Min(dilepton01->Mass(),999.999));
674     hDGenLeptons[22]->Fill(TMath::Min(dilepton02->Mass(),999.999));
675     hDGenLeptons[22]->Fill(TMath::Min(dilepton03->Mass(),999.999));
676     hDGenLeptons[22]->Fill(TMath::Min(dilepton12->Mass(),999.999));
677     hDGenLeptons[22]->Fill(TMath::Min(dilepton13->Mass(),999.999));
678     hDGenLeptons[22]->Fill(TMath::Min(dilepton23->Mass(),999.999));
679     CompositeParticle *fourlepton = new CompositeParticle();
680     fourlepton->AddDaughter(GenLeptons->At(0));
681     fourlepton->AddDaughter(GenLeptons->At(1));
682     fourlepton->AddDaughter(GenLeptons->At(2));
683     fourlepton->AddDaughter(GenLeptons->At(3));
684     hDGenLeptons[23]->Fill(TMath::Min(fourlepton->Mass(),999.999));
685 loizides 1.40 Double_t deltaR[6] = {MathUtils::DeltaR(*GenLeptons->At(0),
686     *GenLeptons->At(1)),
687     MathUtils::DeltaR(*GenLeptons->At(0),
688     *GenLeptons->At(2)),
689     MathUtils::DeltaR(*GenLeptons->At(0),
690     *GenLeptons->At(3)),
691     MathUtils::DeltaR(*GenLeptons->At(1),
692     *GenLeptons->At(2)),
693     MathUtils::DeltaR(*GenLeptons->At(1),
694     *GenLeptons->At(3)),
695     MathUtils::DeltaR(*GenLeptons->At(2),
696     *GenLeptons->At(3))};
697 loizides 1.27 Double_t deltaRMin = deltaR[0];
698 loizides 1.40 for(Int_t i=1; i<6; i++)
699     if(deltaRMin > deltaR[i])
700     deltaRMin = deltaR[i];
701 ceballos 1.22 hDGenLeptons[24]->Fill(deltaRMin);
702    
703     delete dilepton01;
704     delete dilepton02;
705     delete dilepton03;
706     delete dilepton12;
707     delete dilepton13;
708     delete dilepton23;
709     delete fourlepton;
710     }
711     }
712     }
713     }
714 loizides 1.1
715 loizides 1.6 // all leptons
716 ceballos 1.3 hDGenAllLeptons[0]->Fill(GenAllLeptons->GetEntries());
717 loizides 1.5 for(UInt_t i=0; i<GenAllLeptons->GetEntries(); i++) {
718 ceballos 1.3 hDGenAllLeptons[1]->Fill(GenAllLeptons->At(i)->Pt());
719     hDGenAllLeptons[2]->Fill(GenAllLeptons->At(i)->Eta());
720 loizides 1.5 hDGenAllLeptons[3]->Fill(GenAllLeptons->At(i)->PhiDeg());
721 ceballos 1.3 }
722    
723 loizides 1.6 // taus
724 loizides 1.1 hDGenTaus[0]->Fill(GenTaus->GetEntries());
725 loizides 1.5 for(UInt_t i=0; i<GenTaus->GetEntries(); i++) {
726 loizides 1.1 hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
727     hDGenTaus[2]->Fill(GenTaus->At(i)->Eta());
728 loizides 1.5 hDGenTaus[3]->Fill(GenTaus->At(i)->PhiDeg());
729 loizides 1.1 }
730    
731 loizides 1.6 // neutrinos
732 loizides 1.1 hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
733     CompositeParticle *neutrinoTotal = new CompositeParticle();
734 loizides 1.5 for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++) {
735     if (GenNeutrinos->At(i)->HasMother())
736 loizides 1.1 neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
737     }
738 loizides 1.5 if (GenNeutrinos->GetEntries() > 0) {
739 loizides 1.1 hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
740     hDGenNeutrinos[2]->Fill(neutrinoTotal->Eta());
741 loizides 1.5 hDGenNeutrinos[3]->Fill(neutrinoTotal->PhiDeg());
742 loizides 1.1 }
743     delete neutrinoTotal;
744    
745 loizides 1.6 // quarks
746 loizides 1.1 hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
747 loizides 1.5 for(UInt_t i=0; i<GenQuarks->GetEntries(); i++) {
748     for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++) {
749 loizides 1.1 CompositeParticle *dijet = new CompositeParticle();
750     dijet->AddDaughter(GenQuarks->At(i));
751     dijet->AddDaughter(GenQuarks->At(j));
752     hDGenQuarks[1]->Fill(dijet->Pt());
753     hDGenQuarks[2]->Fill(dijet->Mass());
754 loizides 1.5 if (TMath::Abs(GenQuarks->At(i)->Eta()) < 2.5 &&
755     TMath::Abs(GenQuarks->At(j)->Eta()) < 2.5) {
756 loizides 1.1 hDGenQuarks[3]->Fill(dijet->Pt());
757     hDGenQuarks[4]->Fill(dijet->Mass());
758     }
759     delete dijet;
760     }
761 ceballos 1.2 // b quark info
762 loizides 1.5 if (GenQuarks->At(i)->AbsPdgId() == 5) {
763 ceballos 1.2 hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
764     hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());
765     hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());
766 ceballos 1.22 if (GenLeptons->GetEntries() >= 2 &&
767     GenLeptons->At(0)->Pt() > 20 &&
768     GenLeptons->At(1)->Pt() > 15) {
769     if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
770     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5) {
771 ceballos 1.2 hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());
772     hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());
773     hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());
774     }
775     }
776     }
777     // t quark info
778 loizides 1.5 else if (GenQuarks->At(i)->AbsPdgId() == 6) {
779 ceballos 1.2 hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
780     hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());
781     hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());
782     }
783     // light quark info
784     else {
785     hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
786     hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());
787     hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());
788     }
789 loizides 1.1 }
790    
791 loizides 1.6 // wbf
792 loizides 1.5 if (GenqqHs->GetEntries() == 2) {
793 loizides 1.1 hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
794     GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
795 loizides 1.5 hDGenWBF[1]->Fill(TMath::Abs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
796 loizides 1.1 hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
797     hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
798     CompositeParticle *diqq = new CompositeParticle();
799     diqq->AddDaughter(GenqqHs->At(0));
800     diqq->AddDaughter(GenqqHs->At(1));
801     hDGenWBF[4]->Fill(diqq->Mass());
802     delete diqq;
803     }
804    
805 loizides 1.6 // bosons
806 loizides 1.1 hDGenBosons[0]->Fill(GenBosons->GetEntries());
807 loizides 1.5 for(UInt_t i=0; i<GenBosons->GetEntries(); i++) {
808 loizides 1.1 hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
809     hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
810 ceballos 1.42 hDGenBosons[3]->Fill(TMath::Min(GenBosons->At(i)->Mass(),1999.999));
811     hDGenBosons[4]->Fill(TMath::Min(GenBosons->At(i)->Mass(),199.999));
812     if(GenBosons->At(i)->Is(MCParticle::kW))
813     hDGenBosons[5]->Fill(TMath::Min(GenBosons->At(i)->Mass(),199.999));
814     if(GenBosons->At(i)->Is(MCParticle::kZ))
815     hDGenBosons[6]->Fill(TMath::Min(GenBosons->At(i)->Mass(),199.999));
816     }
817     if(sumV[0] >= 4) printf("More than 3 W bosons (%d)\n",sumV[0]);
818     if(sumV[1] >= 4) printf("More than 3 Z bosons (%d)\n",sumV[1]);
819     hDGenBosons[7]->Fill(TMath::Min((double)(sumV[0] + 4*sumV[1]),12.4999));
820 ceballos 1.10
821     // photons
822     hDGenPhotons[0]->Fill(GenPhotons->GetEntries());
823     for(UInt_t i=0; i<GenPhotons->GetEntries(); i++) {
824     hDGenPhotons[1]->Fill(GenPhotons->At(i)->Pt());
825     hDGenPhotons[2]->Fill(GenPhotons->At(i)->Eta());
826     }
827 ceballos 1.32
828     // Rad photons
829     hDGenRadPhotons[0]->Fill(GenRadPhotons->GetEntries());
830     for(UInt_t i=0; i<GenRadPhotons->GetEntries(); i++) {
831     hDGenRadPhotons[1]->Fill(TMath::Min(GenRadPhotons->At(i)->Pt(),199.999));
832     hDGenRadPhotons[2]->Fill(TMath::Min(GenRadPhotons->At(i)->AbsEta(),4.999));
833     hDGenRadPhotons[3]->Fill(TMath::Min((double)GenRadPhotons->At(i)->Mother()->Status(),19.499));
834 loizides 1.40 hDGenRadPhotons[4]->Fill(GenRadPhotons->At(i)->IsGenerated() +
835     2*GenRadPhotons->At(i)->IsSimulated());
836 ceballos 1.32 hDGenRadPhotons[5]->Fill(TMath::Min(
837 loizides 1.40 MathUtils::DeltaR(*GenRadPhotons->At(i),
838     *GenRadPhotons->At(i)->Mother()),
839     4.999));
840 ceballos 1.32 Int_t Mother = 0;
841     if(GenRadPhotons->At(i)->Mother()->Is(MCParticle::kMu)) Mother = 1;
842     hDGenRadPhotons[6]->Fill(Mother);
843     }
844    
845     // ISR photons
846     hDGenISRPhotons[0]->Fill(GenISRPhotons->GetEntries());
847     for(UInt_t i=0; i<GenISRPhotons->GetEntries(); i++) {
848     hDGenISRPhotons[1]->Fill(TMath::Min(GenISRPhotons->At(i)->Pt(),199.999));
849     hDGenISRPhotons[2]->Fill(TMath::Min(GenISRPhotons->At(i)->AbsEta(),4.999));
850 loizides 1.40 hDGenISRPhotons[3]->Fill(TMath::Min((Double_t)GenISRPhotons->At(i)->Mother()->Status(),
851     19.499));
852     hDGenISRPhotons[4]->Fill(GenISRPhotons->At(i)->IsGenerated() +
853     2*GenISRPhotons->At(i)->IsSimulated());
854 ceballos 1.32 hDGenISRPhotons[5]->Fill(TMath::Min(
855 loizides 1.40 MathUtils::DeltaR(*GenISRPhotons->At(i),
856     *GenISRPhotons->At(i)->Mother()),4.999));
857 ceballos 1.32 }
858 loizides 1.1 }
859 ceballos 1.34
860     // Apply ISR filter (but filling all histograms)
861 ceballos 1.37 if(fApplyISRFilter == kTRUE && GenISRPhotons->GetEntries() > 0 &&
862 ceballos 1.38 GenISRPhotons->At(0)->Pt() > 15.0){
863 ceballos 1.34 SkipEvent();
864     }
865 loizides 1.1 }
866    
867     //--------------------------------------------------------------------------------------------------
868     void GeneratorMod::SlaveBegin()
869     {
870 loizides 1.5 // Book branch and histograms if wanted.
871    
872 loizides 1.41 ReqEventObject(fMCPartName, fParticles, kTRUE);
873 loizides 1.1
874 loizides 1.5 // fill histograms
875 loizides 1.20 if (GetFillHist()) {
876 loizides 1.5 char sb[1024];
877 ceballos 1.30 // MET
878     sprintf(sb,"hDGenMet_%d", 0); hDGenMet[0] = new TH1D(sb,sb,200,0,200);
879     sprintf(sb,"hDGenMet_%d", 1); hDGenMet[1] = new TH1D(sb,sb,400,-200,200);
880     sprintf(sb,"hDGenMet_%d", 2); hDGenMet[2] = new TH1D(sb,sb,400,-200,200);
881     sprintf(sb,"hDGenMet_%d", 3); hDGenMet[3] = new TH1D(sb,sb,400,-1000,1000);
882     for(Int_t i=0; i<4; i++) AddOutput(hDGenMet[i]);
883    
884 loizides 1.6 // leptons from W
885 loizides 1.1 sprintf(sb,"hDGenLeptons_%d", 0); hDGenLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
886     sprintf(sb,"hDGenLeptons_%d", 1); hDGenLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
887     sprintf(sb,"hDGenLeptons_%d", 2); hDGenLeptons[2] = new TH1D(sb,sb,50,0.0,5.0);
888     sprintf(sb,"hDGenLeptons_%d", 3); hDGenLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
889     sprintf(sb,"hDGenLeptons_%d", 4); hDGenLeptons[4] = new TH1D(sb,sb,1000,0.0,1000.0);
890     sprintf(sb,"hDGenLeptons_%d", 5); hDGenLeptons[5] = new TH1D(sb,sb,50,0.0,5.0);
891     sprintf(sb,"hDGenLeptons_%d", 6); hDGenLeptons[6] = new TH1D(sb,sb,50,0.0,5.0);
892     sprintf(sb,"hDGenLeptons_%d", 7); hDGenLeptons[7] = new TH1D(sb,sb,100,0.0,200.0);
893     sprintf(sb,"hDGenLeptons_%d", 8); hDGenLeptons[8] = new TH1D(sb,sb,100,0.0,200.0);
894     sprintf(sb,"hDGenLeptons_%d", 9); hDGenLeptons[9] = new TH1D(sb,sb,1000,0.0,1000.0);
895     sprintf(sb,"hDGenLeptons_%d",10); hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0);
896 ceballos 1.22 sprintf(sb,"hDGenLeptons_%d",11); hDGenLeptons[11] = new TH1D(sb,sb,100,0.0,5.0);
897     sprintf(sb,"hDGenLeptons_%d",12); hDGenLeptons[12] = new TH1D(sb,sb,100,0.0,200.0);
898     sprintf(sb,"hDGenLeptons_%d",13); hDGenLeptons[13] = new TH1D(sb,sb,100,0.0,200.0);
899     sprintf(sb,"hDGenLeptons_%d",14); hDGenLeptons[14] = new TH1D(sb,sb,100,0.0,200.0);
900     sprintf(sb,"hDGenLeptons_%d",15); hDGenLeptons[15] = new TH1D(sb,sb,1000,0.0,1000.0);
901     sprintf(sb,"hDGenLeptons_%d",16); hDGenLeptons[16] = new TH1D(sb,sb,1000,0.0,1000.0);
902     sprintf(sb,"hDGenLeptons_%d",17); hDGenLeptons[17] = new TH1D(sb,sb,100,0.0,5.0);
903     sprintf(sb,"hDGenLeptons_%d",18); hDGenLeptons[18] = new TH1D(sb,sb,100,0.0,200.0);
904     sprintf(sb,"hDGenLeptons_%d",19); hDGenLeptons[19] = new TH1D(sb,sb,100,0.0,200.0);
905     sprintf(sb,"hDGenLeptons_%d",20); hDGenLeptons[20] = new TH1D(sb,sb,100,0.0,200.0);
906     sprintf(sb,"hDGenLeptons_%d",21); hDGenLeptons[21] = new TH1D(sb,sb,100,0.0,200.0);
907     sprintf(sb,"hDGenLeptons_%d",22); hDGenLeptons[22] = new TH1D(sb,sb,1000,0.0,1000.0);
908     sprintf(sb,"hDGenLeptons_%d",23); hDGenLeptons[23] = new TH1D(sb,sb,1000,0.0,1000.0);
909     sprintf(sb,"hDGenLeptons_%d",24); hDGenLeptons[24] = new TH1D(sb,sb,100,0.0,5.0);
910     for(Int_t i=0; i<25; i++) AddOutput(hDGenLeptons[i]);
911 loizides 1.1
912 loizides 1.6 // all leptons
913 ceballos 1.3 sprintf(sb,"hDGenAllLeptons_%d", 0); hDGenAllLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
914     sprintf(sb,"hDGenAllLeptons_%d", 1); hDGenAllLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
915 ceballos 1.14 sprintf(sb,"hDGenAllLeptons_%d", 2); hDGenAllLeptons[2] = new TH1D(sb,sb,100,-5.0,5.0);
916 ceballos 1.3 sprintf(sb,"hDGenAllLeptons_%d", 3); hDGenAllLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
917 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenAllLeptons[i]);
918 ceballos 1.3
919 loizides 1.6 // taus
920 loizides 1.1 sprintf(sb,"hDGenTaus_%d", 0); hDGenTaus[0] = new TH1D(sb,sb,10,-0.5,9.5);
921     sprintf(sb,"hDGenTaus_%d", 1); hDGenTaus[1] = new TH1D(sb,sb,100,0.0,200.0);
922     sprintf(sb,"hDGenTaus_%d", 2); hDGenTaus[2] = new TH1D(sb,sb,100,-5.0,5.0);
923     sprintf(sb,"hDGenTaus_%d", 3); hDGenTaus[3] = new TH1D(sb,sb,90,0.0,180.0);
924 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenTaus[i]);
925 loizides 1.1
926 loizides 1.6 // neutrinos
927 loizides 1.1 sprintf(sb,"hDGenNeutrinos_%d", 0); hDGenNeutrinos[0] = new TH1D(sb,sb,10,-0.5,9.5);
928     sprintf(sb,"hDGenNeutrinos_%d", 1); hDGenNeutrinos[1] = new TH1D(sb,sb,100,0.0,200.0);
929     sprintf(sb,"hDGenNeutrinos_%d", 2); hDGenNeutrinos[2] = new TH1D(sb,sb,100,-5.0,5.0);
930     sprintf(sb,"hDGenNeutrinos_%d", 3); hDGenNeutrinos[3] = new TH1D(sb,sb,90,0.0,180.0);
931 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);
932 loizides 1.1
933 loizides 1.6 // quarks
934 loizides 1.1 sprintf(sb,"hDGenQuarks_%d", 0); hDGenQuarks[0] = new TH1D(sb,sb,10,-0.5,9.5);
935     sprintf(sb,"hDGenQuarks_%d", 1); hDGenQuarks[1] = new TH1D(sb,sb,200,0.0,400.);
936     sprintf(sb,"hDGenQuarks_%d", 2); hDGenQuarks[2] = new TH1D(sb,sb,2000,0.0,2000.);
937     sprintf(sb,"hDGenQuarks_%d", 3); hDGenQuarks[3] = new TH1D(sb,sb,200,0.0,400.);
938     sprintf(sb,"hDGenQuarks_%d", 4); hDGenQuarks[4] = new TH1D(sb,sb,2000,0.0,2000.);
939 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 5); hDGenQuarks[5] = new TH1D(sb,sb,200,0.0,400.);
940     sprintf(sb,"hDGenQuarks_%d", 6); hDGenQuarks[6] = new TH1D(sb,sb,200,-10.0,10.);
941 loizides 1.5 sprintf(sb,"hDGenQuarks_%d", 7); hDGenQuarks[7] = new TH1D(sb,sb,200,-TMath::Pi(),
942     TMath::Pi());
943 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 8); hDGenQuarks[8] = new TH1D(sb,sb,200,0.0,400.);
944     sprintf(sb,"hDGenQuarks_%d", 9); hDGenQuarks[9] = new TH1D(sb,sb,200,-10.0,10.);
945 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",10); hDGenQuarks[10] = new TH1D(sb,sb,200,-TMath::Pi(),
946     TMath::Pi());
947 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d",11); hDGenQuarks[11] = new TH1D(sb,sb,200,0.0,400.);
948     sprintf(sb,"hDGenQuarks_%d",12); hDGenQuarks[12] = new TH1D(sb,sb,200,-10.0,10.);
949 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",13); hDGenQuarks[13] = new TH1D(sb,sb,200,-TMath::Pi(),
950     TMath::Pi());
951 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d",14); hDGenQuarks[14] = new TH1D(sb,sb,200,0.0,400.);
952     sprintf(sb,"hDGenQuarks_%d",15); hDGenQuarks[15] = new TH1D(sb,sb,200,-10.0,10.);
953 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",16); hDGenQuarks[16] = new TH1D(sb,sb,200,-TMath::Pi(),
954     TMath::Pi());
955     for(Int_t i=0; i<17; i++) AddOutput(hDGenQuarks[i]);
956 loizides 1.1
957     // qqH
958     sprintf(sb,"hDGenWBF_%d", 0); hDGenWBF[0] = new TH1D(sb,sb,90,0.0,180.);
959     sprintf(sb,"hDGenWBF_%d", 1); hDGenWBF[1] = new TH1D(sb,sb,100,0.0,10.);
960     sprintf(sb,"hDGenWBF_%d", 2); hDGenWBF[2] = new TH1D(sb,sb,200,0.0,400.);
961     sprintf(sb,"hDGenWBF_%d", 3); hDGenWBF[3] = new TH1D(sb,sb,200,0.0,400.);
962     sprintf(sb,"hDGenWBF_%d", 4); hDGenWBF[4] = new TH1D(sb,sb,200,0.0,4000.);
963 loizides 1.5 for(Int_t i=0; i<5; i++) AddOutput(hDGenWBF[i]);
964 loizides 1.1
965 loizides 1.6 // bosons
966 loizides 1.1 sprintf(sb,"hDGenBosons_%d", 0); hDGenBosons[0] = new TH1D(sb,sb,10,-0.5,9.5);
967     sprintf(sb,"hDGenBosons_%d", 1); hDGenBosons[1] = new TH1D(sb,sb,200,0.0,400.0);
968     sprintf(sb,"hDGenBosons_%d", 2); hDGenBosons[2] = new TH1D(sb,sb,100,-5.0,5.0);
969 ceballos 1.42 sprintf(sb,"hDGenBosons_%d", 3); hDGenBosons[3] = new TH1D(sb,sb,2000,0.0,2000.0);
970     sprintf(sb,"hDGenBosons_%d", 4); hDGenBosons[4] = new TH1D(sb,sb,200,0.0,200.0);
971     sprintf(sb,"hDGenBosons_%d", 5); hDGenBosons[5] = new TH1D(sb,sb,200,0.0,200.0);
972     sprintf(sb,"hDGenBosons_%d", 6); hDGenBosons[6] = new TH1D(sb,sb,200,0.0,200.0);
973     sprintf(sb,"hDGenBosons_%d", 7); hDGenBosons[7] = new TH1D(sb,sb,13,-0.5,12.5);
974     for(Int_t i=0; i<8; i++) AddOutput(hDGenBosons[i]);
975 ceballos 1.10
976     // photons
977     sprintf(sb,"hDGenPhotons_%d", 0); hDGenPhotons[0] = new TH1D(sb,sb,10,-0.5,9.5);
978     sprintf(sb,"hDGenPhotons_%d", 1); hDGenPhotons[1] = new TH1D(sb,sb,200,0.0,400.0);
979 ceballos 1.14 sprintf(sb,"hDGenPhotons_%d", 2); hDGenPhotons[2] = new TH1D(sb,sb,100,-5.0,5.0);
980 ceballos 1.10 for(Int_t i=0; i<3; i++) AddOutput(hDGenPhotons[i]);
981 ceballos 1.16
982 loizides 1.40 // rad photons
983 ceballos 1.32 sprintf(sb,"hDGenRadPhotons_%d", 0); hDGenRadPhotons[0] = new TH1D(sb,sb,10,-0.5,9.5);
984     sprintf(sb,"hDGenRadPhotons_%d", 1); hDGenRadPhotons[1] = new TH1D(sb,sb,400,0.0,200.0);
985     sprintf(sb,"hDGenRadPhotons_%d", 2); hDGenRadPhotons[2] = new TH1D(sb,sb,100,0.0,5.0);
986     sprintf(sb,"hDGenRadPhotons_%d", 3); hDGenRadPhotons[3] = new TH1D(sb,sb,20,-0.5,19.5);
987     sprintf(sb,"hDGenRadPhotons_%d", 4); hDGenRadPhotons[4] = new TH1D(sb,sb,4,-0.5,3.5);
988     sprintf(sb,"hDGenRadPhotons_%d", 5); hDGenRadPhotons[5] = new TH1D(sb,sb,500,0.0,5.0);
989     sprintf(sb,"hDGenRadPhotons_%d", 6); hDGenRadPhotons[6] = new TH1D(sb,sb,2,-0.5,1.5);
990     for(Int_t i=0; i<7; i++) AddOutput(hDGenRadPhotons[i]);
991    
992     // ISR photons
993     sprintf(sb,"hDGenISRPhotons_%d", 0); hDGenISRPhotons[0] = new TH1D(sb,sb,10,-0.5,9.5);
994     sprintf(sb,"hDGenISRPhotons_%d", 1); hDGenISRPhotons[1] = new TH1D(sb,sb,400,0.0,200.0);
995     sprintf(sb,"hDGenISRPhotons_%d", 2); hDGenISRPhotons[2] = new TH1D(sb,sb,100,0.0,5.0);
996     sprintf(sb,"hDGenISRPhotons_%d", 3); hDGenISRPhotons[3] = new TH1D(sb,sb,20,-0.5,19.5);
997     sprintf(sb,"hDGenISRPhotons_%d", 4); hDGenISRPhotons[4] = new TH1D(sb,sb,4,-0.5,3.5);
998     sprintf(sb,"hDGenISRPhotons_%d", 5); hDGenISRPhotons[5] = new TH1D(sb,sb,500,0.0,5.0);
999     for(Int_t i=0; i<6; i++) AddOutput(hDGenISRPhotons[i]);
1000    
1001 ceballos 1.42 // auxiliar for MG studies
1002 ceballos 1.16 sprintf(sb,"hDVMass_%d", 0); hDVMass[0] = new TH1D(sb,sb,200,0.,200.);
1003     sprintf(sb,"hDVMass_%d", 1); hDVMass[1] = new TH1D(sb,sb,200,0.,200.);
1004     sprintf(sb,"hDVMass_%d", 2); hDVMass[2] = new TH1D(sb,sb,200,0.,200.);
1005     sprintf(sb,"hDVMass_%d", 3); hDVMass[3] = new TH1D(sb,sb,200,0.,200.);
1006     sprintf(sb,"hDVMass_%d", 4); hDVMass[4] = new TH1D(sb,sb,200,0.,200.);
1007     sprintf(sb,"hDVMass_%d", 5); hDVMass[5] = new TH1D(sb,sb,200,0.,200.);
1008     sprintf(sb,"hDVMass_%d", 6); hDVMass[6] = new TH1D(sb,sb,200,0.,200.);
1009     sprintf(sb,"hDVMass_%d", 7); hDVMass[7] = new TH1D(sb,sb,200,0.,200.);
1010     sprintf(sb,"hDVMass_%d", 8); hDVMass[8] = new TH1D(sb,sb,200,0.,200.);
1011     sprintf(sb,"hDVMass_%d", 9); hDVMass[9] = new TH1D(sb,sb,200,0.,200.);
1012 ceballos 1.42 sprintf(sb,"hDVMass_%d",10); hDVMass[10] = new TH1D(sb,sb,200,0.,200.);
1013     sprintf(sb,"hDVMass_%d",11); hDVMass[11] = new TH1D(sb,sb,200,0.,200.);
1014     sprintf(sb,"hDVMass_%d",12); hDVMass[12] = new TH1D(sb,sb,200,0.,200.);
1015     for(Int_t i=0; i<13; i++) AddOutput(hDVMass[i]);
1016    
1017 loizides 1.1 }
1018     }