ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.42
Committed: Thu Jun 11 13:16:53 2009 UTC (15 years, 10 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_009a
Changes since 1.41: +286 -109 lines
Log Message: 
new stuff to fix MG

File Contents

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