ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.39
Committed: Thu Apr 23 10:50:59 2009 UTC (16 years ago) by phedex
Content type: text/plain
Branch: MAIN
Changes since 1.38: +18 -8 lines
Log Message: 
correct madgraph Z boson section. we were using daughter0 and daughter1 to find leptons. sometimes the leptons are not the first 2 entries.

File Contents

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