ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.20
Committed: Wed Dec 10 12:47:06 2008 UTC (16 years, 4 months ago) by loizides
Content type: text/plain
Branch: MAIN
Changes since 1.19: +13 -15 lines
Log Message: 
Use GetFillHist

File Contents

# User Rev Content
1 loizides 1.20 // $Id: GeneratorMod.cc,v 1.19 2008/12/10 11:44:33 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     fMCPartName(Names::gkMCPartBrn),
17 loizides 1.8 fMCLeptonsName(ModNames::gkMCLeptonsName),
18     fMCAllLeptonsName(ModNames::gkMCAllLeptonsName),
19     fMCTausName(ModNames::gkMCTausName),
20     fMCNeutrinosName(ModNames::gkMCNeutrinosName),
21     fMCQuarksName(ModNames::gkMCQuarksName),
22     fMCqqHsName(ModNames::gkMCqqHsName),
23     fMCBosonsName(ModNames::gkMCBosonsName),
24 ceballos 1.10 fMCPhotonsName(ModNames::gkMCPhotonsName),
25 ceballos 1.14 fPtLeptonMin(0.0),
26     fEtaLeptonMax(5.0),
27     fPtPhotonMin(0.0),
28 loizides 1.15 fEtaPhotonMax(5.0),
29     fParticles(0)
30 loizides 1.1 {
31     // Constructor.
32     }
33    
34     //--------------------------------------------------------------------------------------------------
35     void GeneratorMod::Process()
36     {
37 loizides 1.5 // Process entries of the tree.
38 loizides 1.1
39 loizides 1.5 // these arrays will be filled in the loop of particles
40 loizides 1.9 MCParticleOArr *GenLeptons = new MCParticleOArr;
41 loizides 1.13 GenLeptons->SetName(fMCLeptonsName);
42 loizides 1.9 MCParticleOArr *GenAllLeptons = new MCParticleOArr;
43 loizides 1.13 GenAllLeptons->SetName(fMCAllLeptonsName);
44 loizides 1.9 MCParticleOArr *GenTaus = new MCParticleOArr;
45 loizides 1.13 GenTaus->SetName(fMCTausName);
46     GenTaus->SetOwner(kTRUE);
47 loizides 1.9 MCParticleOArr *GenNeutrinos = new MCParticleOArr;
48 loizides 1.13 GenNeutrinos->SetName(fMCNeutrinosName);
49 loizides 1.9 MCParticleOArr *GenQuarks = new MCParticleOArr;
50 loizides 1.13 GenQuarks->SetName(fMCQuarksName);
51 loizides 1.9 MCParticleOArr *GenqqHs = new MCParticleOArr;
52 loizides 1.13 GenqqHs->SetName(fMCqqHsName);
53 loizides 1.9 MCParticleOArr *GenBosons = new MCParticleOArr;
54 loizides 1.13 GenBosons->SetName(fMCBosonsName);
55 ceballos 1.10 MCParticleOArr *GenPhotons = new MCParticleOArr;
56 loizides 1.13 GenPhotons->SetName(fMCPhotonsName);
57 loizides 1.1
58 loizides 1.5 // load MCParticle branch
59     LoadBranch(fMCPartName);
60    
61 loizides 1.6 Bool_t isqqH = kFALSE;
62 loizides 1.5 for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
63 loizides 1.7 const MCParticle *p = fParticles->At(i);
64 loizides 1.5
65     if (!p->IsGenerated()) continue;
66    
67     // muons/electrons from W/Z decays
68 loizides 1.7 if ((p->Is(MCParticle::kEl) || p->Is(MCParticle::kMu)) && p->Status() == 1) {
69 ceballos 1.14 if (p->Pt() > fPtLeptonMin && p->AbsEta() < fEtaLeptonMax) {
70 loizides 1.5 GenAllLeptons->Add(p);
71     }
72 loizides 1.6 Bool_t isGoodLepton = kFALSE;
73 loizides 1.7 const MCParticle *pm = p;
74 loizides 1.5 while (pm->HasMother() && isGoodLepton == kFALSE) {
75 loizides 1.7 if (pm->Mother()->Is(MCParticle::kZ) || pm->Mother()->Is(MCParticle::kW)) {
76 loizides 1.5 GenLeptons->Add(p);
77     isGoodLepton = kTRUE;
78 loizides 1.7 break;
79 loizides 1.8 } else if (pm->Mother()->Is(MCParticle::kPi0) || pm->Mother()->Is(MCParticle::kEta)) {
80 loizides 1.7 // this is fake, but it is a trick to get rid of these cases and abort the loop
81 loizides 1.5 isGoodLepton = kTRUE;
82 loizides 1.7 break;
83     }
84     pm = pm->Mother();
85 loizides 1.1 }
86 loizides 1.5 }
87    
88 loizides 1.7 // hadronic taus
89     else if (p->Is(MCParticle::kTau) && p->Status() == 2) {
90     if (!p->HasDaughter(MCParticle::kEl) && !p->HasDaughter(MCParticle::kMu)) {
91     const MCParticle *tv = p->FindDaughter(MCParticle::kTauNu);
92     if (tv) {
93     MCParticle *pm_f = new MCParticle(*p);
94     pm_f->SetMom(p->Px()-tv->Px(), p->Py()-tv->Py(),
95     p->Pz()-tv->Pz(), p->E()-tv->E());
96     GenTaus->AddOwned(pm_f);
97     } else {
98 loizides 1.18 SendError(kWarning, "Process", "Could not find a tau neutrino!");
99 loizides 1.5 }
100 loizides 1.1 }
101 loizides 1.5 }
102 loizides 1.1
103 loizides 1.5 // neutrinos
104 loizides 1.7 else if (p->Status() == 1 && p->IsNeutrino()) {
105 loizides 1.5 GenNeutrinos->Add(p);
106     }
107 loizides 1.1
108 loizides 1.5 // quarks from W/Z decays or top particles
109 loizides 1.7 else if (p->IsQuark() && p->HasMother()) {
110     if (p->Mother()->Is(MCParticle::kZ) || p->Mother()->Is(MCParticle::kW) ||
111     p->Is(MCParticle::kTop) || p->Mother()->Is(MCParticle::kTop)) {
112 loizides 1.5 GenQuarks->Add(p);
113 loizides 1.1 }
114 loizides 1.5 }
115 loizides 1.1
116 loizides 1.5 // qqH, information about the forward jets
117 loizides 1.8 else if (isqqH == kFALSE && p->Is(MCParticle::kH)) {
118 loizides 1.5 isqqH = kTRUE;
119 loizides 1.17 const MCParticle *pq1 = fParticles->At(i-1);
120     const MCParticle *pq2 = fParticles->At(i-2);
121 loizides 1.7 if (!pq1 || !pq2) {
122     SendError(kWarning, "Process", "Could not find quark pair!");
123 loizides 1.8 } else if (pq1->IsQuark() && pq2->IsQuark() &&
124     pq1->HasMother() && pq2->HasMother() &&
125     pq1->Mother()->PdgId() == p->Mother()->PdgId() &&
126     pq2->Mother()->PdgId() == p->Mother()->PdgId()) {
127 loizides 1.5 GenqqHs->Add(pq1);
128     GenqqHs->Add(pq2);
129 loizides 1.1 }
130 loizides 1.7 if (p->Status() == 3)
131     GenBosons->Add(p); // take higgs boson in account here rather in next else if
132 loizides 1.5 }
133 loizides 1.1
134 loizides 1.5 // information about bosons: W, Z, h, Z', W', H0, A0, H+
135 loizides 1.7 else if (p->Status() == 3 &&
136     (p->Is(MCParticle::kZ) || p->Is(MCParticle::kW) || p->Is(MCParticle::kH) ||
137     p->Is(MCParticle::kZp) || p->Is(MCParticle::kZpp) ||
138     p->Is(MCParticle::kH0) || p->Is(MCParticle::kA0) || p->Is(MCParticle::kHp))) {
139 loizides 1.5 GenBosons->Add(p);
140 loizides 1.1 }
141 ceballos 1.10
142 ceballos 1.14 // photons
143 ceballos 1.10 else if (p->Status() == 1 && p->Is(MCParticle::kGamma) &&
144 ceballos 1.14 p->Pt() > fPtPhotonMin && p->AbsEta() < fEtaPhotonMax) {
145 ceballos 1.10 GenPhotons->Add(p);
146     }
147    
148 loizides 1.18 // W/Z -> lnu for Madgraph
149     if (p->IsParton() && p->NDaughters() == 2) {
150 ceballos 1.16 CompositeParticle *diBoson = new CompositeParticle();
151     diBoson->AddDaughter(p->Daughter(0));
152     diBoson->AddDaughter(p->Daughter(1));
153 loizides 1.18 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
154 loizides 1.20 if (GetFillHist())
155 loizides 1.18 hDVMass[0]->Fill(TMath::Min(diBoson->Mass(),199.999));
156 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kMu));
157     }
158 loizides 1.18 else if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
159 loizides 1.20 if (GetFillHist())
160 loizides 1.18 hDVMass[1]->Fill(TMath::Min(diBoson->Mass(),199.999));
161 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kEl));
162     }
163 loizides 1.18 else if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
164 loizides 1.20 if (GetFillHist())
165 loizides 1.18 hDVMass[2]->Fill(TMath::Min(diBoson->Mass(),199.999));
166     const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
167     if (tau->HasDaughter(MCParticle::kMu))
168     GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
169     if (tau->HasDaughter(MCParticle::kEl))
170     GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
171     }
172     else if (p->Daughter(0)->Is(MCParticle::kMu) && p->Daughter(1)->Is(MCParticle::kMu)) {
173 loizides 1.20 if (GetFillHist())
174 loizides 1.18 hDVMass[3]->Fill(TMath::Min(diBoson->Mass(),199.999));
175 ceballos 1.16 GenLeptons->Add(p->Daughter(0));
176     GenLeptons->Add(p->Daughter(1));
177     }
178 loizides 1.18 else if (p->Daughter(0)->Is(MCParticle::kEl) && p->Daughter(1)->Is(MCParticle::kEl)) {
179 loizides 1.20 if (GetFillHist())
180 loizides 1.18 hDVMass[4]->Fill(TMath::Min(diBoson->Mass(),199.999));
181 ceballos 1.16 GenLeptons->Add(p->Daughter(0));
182     GenLeptons->Add(p->Daughter(1));
183     }
184 loizides 1.18 else if (p->Daughter(0)->Is(MCParticle::kTau) && p->Daughter(1)->Is(MCParticle::kTau)) {
185 loizides 1.20 if (GetFillHist())
186 loizides 1.18 hDVMass[5]->Fill(TMath::Min(diBoson->Mass(),199.999));
187     const MCParticle *tau0 = p->Daughter(0);
188     if (tau0->HasDaughter(MCParticle::kMu))
189     GenLeptons->Add(tau0->FindDaughter(MCParticle::kMu));
190     if (tau0->HasDaughter(MCParticle::kEl))
191     GenLeptons->Add(tau0->FindDaughter(MCParticle::kEl));
192     const MCParticle *tau1 = p->Daughter(1);
193     if (tau1->HasDaughter(MCParticle::kMu))
194     GenLeptons->Add(tau1->FindDaughter(MCParticle::kMu));
195     if (tau1->HasDaughter(MCParticle::kEl))
196     GenLeptons->Add(tau1->FindDaughter(MCParticle::kEl));
197 ceballos 1.16 }
198     delete diBoson;
199     }
200    
201 loizides 1.18 // t -> lnu for Madgraph
202     if (p->Is(MCParticle::kTop)) {
203 ceballos 1.16 CompositeParticle *diBoson = new CompositeParticle();
204 loizides 1.18 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
205 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
206     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
207 loizides 1.20 if (GetFillHist())
208 loizides 1.18 hDVMass[6]->Fill(TMath::Min(diBoson->Mass(),199.999));
209 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kMu));
210     }
211 loizides 1.18 else if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
212 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
213     diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
214 loizides 1.20 if (GetFillHist())
215 loizides 1.18 hDVMass[7]->Fill(TMath::Min(diBoson->Mass(),199.999));
216 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kEl));
217     }
218 loizides 1.18 else if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
219 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
220     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
221 loizides 1.20 if (GetFillHist())
222 loizides 1.18 hDVMass[8]->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     }
229     else if (!p->HasDaughter(MCParticle::kW)) {
230     for(UInt_t nd=0; nd<p->NDaughters(); ++nd)
231     if (p->Daughter(nd)->IsNot(MCParticle::kBottom) &&
232     p->Daughter(nd)->IsNot(MCParticle::kGamma))
233     diBoson->AddDaughter(p->Daughter(nd));
234 loizides 1.20 if (GetFillHist())
235 loizides 1.18 hDVMass[9]->Fill(TMath::Min(diBoson->Mass(),199.999));
236 ceballos 1.16 }
237     delete diBoson;
238     }
239 loizides 1.5 }
240 loizides 1.1
241 loizides 1.19 // sort according to pt
242     GenLeptons->Sort();
243     GenAllLeptons->Sort();
244     GenTaus->Sort();
245     GenNeutrinos->Sort();
246     GenQuarks->Sort();
247     GenqqHs->Sort();
248     GenBosons->Sort();
249     GenPhotons->Sort();
250    
251 loizides 1.8 // add objects to this event for other modules to use
252 loizides 1.13 AddObjThisEvt(GenLeptons);
253     AddObjThisEvt(GenAllLeptons);
254     AddObjThisEvt(GenTaus);
255     AddObjThisEvt(GenNeutrinos);
256     AddObjThisEvt(GenQuarks);
257     AddObjThisEvt(GenqqHs);
258     AddObjThisEvt(GenBosons);
259     AddObjThisEvt(GenPhotons);
260 ceballos 1.16
261 loizides 1.5 // fill histograms if requested
262 loizides 1.20 if (GetFillHist()) {
263 loizides 1.5
264 loizides 1.6 // leptons
265 loizides 1.1 hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
266 loizides 1.5 Int_t idxMaxLep[2] = {-1, -1};
267     Double_t ptMaxLep[2] = {-1, -1};
268     for(UInt_t i=0; i<GenLeptons->GetEntries(); i++) {
269 loizides 1.1 hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
270 loizides 1.5 hDGenLeptons[2]->Fill(TMath::Abs(GenLeptons->At(i)->Eta()));
271     hDGenLeptons[3]->Fill(GenLeptons->At(i)->PhiDeg());
272     for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++) {
273 loizides 1.1 CompositeParticle *dilepton = new CompositeParticle();
274     dilepton->AddDaughter(GenLeptons->At(i));
275     dilepton->AddDaughter(GenLeptons->At(j));
276     hDGenLeptons[4]->Fill(dilepton->Mass());
277     delete dilepton;
278     }
279 loizides 1.6 // selecting the two highest pt leptons
280 loizides 1.5 if (GenLeptons->At(i)->Pt() > ptMaxLep[0]) {
281 loizides 1.1 ptMaxLep[1] = ptMaxLep[0];
282     idxMaxLep[1] = idxMaxLep[0];
283     ptMaxLep[0] = GenLeptons->At(i)->Pt();
284     idxMaxLep[0] = i;
285     }
286 loizides 1.5 else if (GenLeptons->At(i)->Pt() > ptMaxLep[1]) {
287 loizides 1.1 ptMaxLep[1] = GenLeptons->At(i)->Pt();
288     idxMaxLep[1] = i;
289     }
290     }
291 loizides 1.6 // looking at events with at least two leptons
292 loizides 1.5 if (ptMaxLep[0] > 0 && ptMaxLep[1] > 0) {
293     hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()),
294     TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta())),
295     4.999));
296     hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()),
297     TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta())),
298     4.999));
299     if (TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
300     TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5) {
301 loizides 1.1 hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[0])->Pt(),199.999));
302 loizides 1.5 if (GenLeptons->At(idxMaxLep[0])->Pt() > 20.0) {
303 loizides 1.1 hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[1])->Pt(),199.999));
304 loizides 1.5 if (GenLeptons->At(idxMaxLep[1])->Pt() > 10.0) {
305 loizides 1.1 CompositeParticle *dilepton = new CompositeParticle();
306     dilepton->AddDaughter(GenLeptons->At(idxMaxLep[0]));
307     dilepton->AddDaughter(GenLeptons->At(idxMaxLep[1]));
308     hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
309     hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(idxMaxLep[0])->Phi(),
310     GenLeptons->At(idxMaxLep[1])->Phi())
311     * 180./ TMath::Pi());
312     delete dilepton;
313     }
314     }
315     }
316     }
317    
318 loizides 1.6 // all leptons
319 ceballos 1.3 hDGenAllLeptons[0]->Fill(GenAllLeptons->GetEntries());
320 loizides 1.5 for(UInt_t i=0; i<GenAllLeptons->GetEntries(); i++) {
321 ceballos 1.3 hDGenAllLeptons[1]->Fill(GenAllLeptons->At(i)->Pt());
322     hDGenAllLeptons[2]->Fill(GenAllLeptons->At(i)->Eta());
323 loizides 1.5 hDGenAllLeptons[3]->Fill(GenAllLeptons->At(i)->PhiDeg());
324 ceballos 1.3 }
325    
326 loizides 1.6 // taus
327 loizides 1.1 hDGenTaus[0]->Fill(GenTaus->GetEntries());
328 loizides 1.5 for(UInt_t i=0; i<GenTaus->GetEntries(); i++) {
329 loizides 1.1 hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
330     hDGenTaus[2]->Fill(GenTaus->At(i)->Eta());
331 loizides 1.5 hDGenTaus[3]->Fill(GenTaus->At(i)->PhiDeg());
332 loizides 1.1 }
333    
334 loizides 1.6 // neutrinos
335 loizides 1.1 hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
336     CompositeParticle *neutrinoTotal = new CompositeParticle();
337 loizides 1.5 for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++) {
338     if (GenNeutrinos->At(i)->HasMother())
339 loizides 1.1 neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
340     }
341 loizides 1.5 if (GenNeutrinos->GetEntries() > 0) {
342 loizides 1.1 hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
343     hDGenNeutrinos[2]->Fill(neutrinoTotal->Eta());
344 loizides 1.5 hDGenNeutrinos[3]->Fill(neutrinoTotal->PhiDeg());
345 loizides 1.1 }
346     delete neutrinoTotal;
347    
348 loizides 1.6 // quarks
349 loizides 1.1 hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
350 loizides 1.5 for(UInt_t i=0; i<GenQuarks->GetEntries(); i++) {
351     for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++) {
352 loizides 1.1 CompositeParticle *dijet = new CompositeParticle();
353     dijet->AddDaughter(GenQuarks->At(i));
354     dijet->AddDaughter(GenQuarks->At(j));
355     hDGenQuarks[1]->Fill(dijet->Pt());
356     hDGenQuarks[2]->Fill(dijet->Mass());
357 loizides 1.5 if (TMath::Abs(GenQuarks->At(i)->Eta()) < 2.5 &&
358     TMath::Abs(GenQuarks->At(j)->Eta()) < 2.5) {
359 loizides 1.1 hDGenQuarks[3]->Fill(dijet->Pt());
360     hDGenQuarks[4]->Fill(dijet->Mass());
361     }
362     delete dijet;
363     }
364 ceballos 1.2 // b quark info
365 loizides 1.5 if (GenQuarks->At(i)->AbsPdgId() == 5) {
366 ceballos 1.2 hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
367     hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());
368     hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());
369 loizides 1.5 if (ptMaxLep[0] > 20 && ptMaxLep[1] > 15) {
370     if (TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
371     TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5) {
372 ceballos 1.2 hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());
373     hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());
374     hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());
375     }
376     }
377     }
378     // t quark info
379 loizides 1.5 else if (GenQuarks->At(i)->AbsPdgId() == 6) {
380 ceballos 1.2 hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
381     hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());
382     hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());
383     }
384     // light quark info
385     else {
386     hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
387     hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());
388     hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());
389     }
390 loizides 1.1 }
391    
392 loizides 1.6 // wbf
393 loizides 1.5 if (GenqqHs->GetEntries() == 2) {
394 loizides 1.1 hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
395     GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
396 loizides 1.5 hDGenWBF[1]->Fill(TMath::Abs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
397 loizides 1.1 hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
398     hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
399     CompositeParticle *diqq = new CompositeParticle();
400     diqq->AddDaughter(GenqqHs->At(0));
401     diqq->AddDaughter(GenqqHs->At(1));
402     hDGenWBF[4]->Fill(diqq->Mass());
403     delete diqq;
404     }
405    
406 loizides 1.6 // bosons
407 loizides 1.1 hDGenBosons[0]->Fill(GenBosons->GetEntries());
408 loizides 1.5 for(UInt_t i=0; i<GenBosons->GetEntries(); i++) {
409 loizides 1.1 hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
410     hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
411     hDGenBosons[3]->Fill(GenBosons->At(i)->Mass());
412     }
413 ceballos 1.10
414     // photons
415     hDGenPhotons[0]->Fill(GenPhotons->GetEntries());
416     for(UInt_t i=0; i<GenPhotons->GetEntries(); i++) {
417     hDGenPhotons[1]->Fill(GenPhotons->At(i)->Pt());
418     hDGenPhotons[2]->Fill(GenPhotons->At(i)->Eta());
419     }
420 loizides 1.1 }
421     }
422    
423     //--------------------------------------------------------------------------------------------------
424     void GeneratorMod::SlaveBegin()
425     {
426 loizides 1.5 // Book branch and histograms if wanted.
427    
428 loizides 1.1 ReqBranch(fMCPartName, fParticles);
429    
430 loizides 1.5 // fill histograms
431 loizides 1.20 if (GetFillHist()) {
432 loizides 1.5 char sb[1024];
433 loizides 1.6 // leptons from W
434 loizides 1.1 sprintf(sb,"hDGenLeptons_%d", 0); hDGenLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
435     sprintf(sb,"hDGenLeptons_%d", 1); hDGenLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
436     sprintf(sb,"hDGenLeptons_%d", 2); hDGenLeptons[2] = new TH1D(sb,sb,50,0.0,5.0);
437     sprintf(sb,"hDGenLeptons_%d", 3); hDGenLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
438     sprintf(sb,"hDGenLeptons_%d", 4); hDGenLeptons[4] = new TH1D(sb,sb,1000,0.0,1000.0);
439     sprintf(sb,"hDGenLeptons_%d", 5); hDGenLeptons[5] = new TH1D(sb,sb,50,0.0,5.0);
440     sprintf(sb,"hDGenLeptons_%d", 6); hDGenLeptons[6] = new TH1D(sb,sb,50,0.0,5.0);
441     sprintf(sb,"hDGenLeptons_%d", 7); hDGenLeptons[7] = new TH1D(sb,sb,100,0.0,200.0);
442     sprintf(sb,"hDGenLeptons_%d", 8); hDGenLeptons[8] = new TH1D(sb,sb,100,0.0,200.0);
443     sprintf(sb,"hDGenLeptons_%d", 9); hDGenLeptons[9] = new TH1D(sb,sb,1000,0.0,1000.0);
444     sprintf(sb,"hDGenLeptons_%d",10); hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0);
445 loizides 1.5 for(Int_t i=0; i<11; i++) AddOutput(hDGenLeptons[i]);
446 loizides 1.1
447 loizides 1.6 // all leptons
448 ceballos 1.3 sprintf(sb,"hDGenAllLeptons_%d", 0); hDGenAllLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
449     sprintf(sb,"hDGenAllLeptons_%d", 1); hDGenAllLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
450 ceballos 1.14 sprintf(sb,"hDGenAllLeptons_%d", 2); hDGenAllLeptons[2] = new TH1D(sb,sb,100,-5.0,5.0);
451 ceballos 1.3 sprintf(sb,"hDGenAllLeptons_%d", 3); hDGenAllLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
452 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenAllLeptons[i]);
453 ceballos 1.3
454 loizides 1.6 // taus
455 loizides 1.1 sprintf(sb,"hDGenTaus_%d", 0); hDGenTaus[0] = new TH1D(sb,sb,10,-0.5,9.5);
456     sprintf(sb,"hDGenTaus_%d", 1); hDGenTaus[1] = new TH1D(sb,sb,100,0.0,200.0);
457     sprintf(sb,"hDGenTaus_%d", 2); hDGenTaus[2] = new TH1D(sb,sb,100,-5.0,5.0);
458     sprintf(sb,"hDGenTaus_%d", 3); hDGenTaus[3] = new TH1D(sb,sb,90,0.0,180.0);
459 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenTaus[i]);
460 loizides 1.1
461 loizides 1.6 // neutrinos
462 loizides 1.1 sprintf(sb,"hDGenNeutrinos_%d", 0); hDGenNeutrinos[0] = new TH1D(sb,sb,10,-0.5,9.5);
463     sprintf(sb,"hDGenNeutrinos_%d", 1); hDGenNeutrinos[1] = new TH1D(sb,sb,100,0.0,200.0);
464     sprintf(sb,"hDGenNeutrinos_%d", 2); hDGenNeutrinos[2] = new TH1D(sb,sb,100,-5.0,5.0);
465     sprintf(sb,"hDGenNeutrinos_%d", 3); hDGenNeutrinos[3] = new TH1D(sb,sb,90,0.0,180.0);
466 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);
467 loizides 1.1
468 loizides 1.6 // quarks
469 loizides 1.1 sprintf(sb,"hDGenQuarks_%d", 0); hDGenQuarks[0] = new TH1D(sb,sb,10,-0.5,9.5);
470     sprintf(sb,"hDGenQuarks_%d", 1); hDGenQuarks[1] = new TH1D(sb,sb,200,0.0,400.);
471     sprintf(sb,"hDGenQuarks_%d", 2); hDGenQuarks[2] = new TH1D(sb,sb,2000,0.0,2000.);
472     sprintf(sb,"hDGenQuarks_%d", 3); hDGenQuarks[3] = new TH1D(sb,sb,200,0.0,400.);
473     sprintf(sb,"hDGenQuarks_%d", 4); hDGenQuarks[4] = new TH1D(sb,sb,2000,0.0,2000.);
474 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 5); hDGenQuarks[5] = new TH1D(sb,sb,200,0.0,400.);
475     sprintf(sb,"hDGenQuarks_%d", 6); hDGenQuarks[6] = new TH1D(sb,sb,200,-10.0,10.);
476 loizides 1.5 sprintf(sb,"hDGenQuarks_%d", 7); hDGenQuarks[7] = new TH1D(sb,sb,200,-TMath::Pi(),
477     TMath::Pi());
478 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 8); hDGenQuarks[8] = new TH1D(sb,sb,200,0.0,400.);
479     sprintf(sb,"hDGenQuarks_%d", 9); hDGenQuarks[9] = new TH1D(sb,sb,200,-10.0,10.);
480 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",10); hDGenQuarks[10] = new TH1D(sb,sb,200,-TMath::Pi(),
481     TMath::Pi());
482 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d",11); hDGenQuarks[11] = new TH1D(sb,sb,200,0.0,400.);
483     sprintf(sb,"hDGenQuarks_%d",12); hDGenQuarks[12] = new TH1D(sb,sb,200,-10.0,10.);
484 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",13); hDGenQuarks[13] = new TH1D(sb,sb,200,-TMath::Pi(),
485     TMath::Pi());
486 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d",14); hDGenQuarks[14] = new TH1D(sb,sb,200,0.0,400.);
487     sprintf(sb,"hDGenQuarks_%d",15); hDGenQuarks[15] = new TH1D(sb,sb,200,-10.0,10.);
488 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",16); hDGenQuarks[16] = new TH1D(sb,sb,200,-TMath::Pi(),
489     TMath::Pi());
490     for(Int_t i=0; i<17; i++) AddOutput(hDGenQuarks[i]);
491 loizides 1.1
492     // qqH
493     sprintf(sb,"hDGenWBF_%d", 0); hDGenWBF[0] = new TH1D(sb,sb,90,0.0,180.);
494     sprintf(sb,"hDGenWBF_%d", 1); hDGenWBF[1] = new TH1D(sb,sb,100,0.0,10.);
495     sprintf(sb,"hDGenWBF_%d", 2); hDGenWBF[2] = new TH1D(sb,sb,200,0.0,400.);
496     sprintf(sb,"hDGenWBF_%d", 3); hDGenWBF[3] = new TH1D(sb,sb,200,0.0,400.);
497     sprintf(sb,"hDGenWBF_%d", 4); hDGenWBF[4] = new TH1D(sb,sb,200,0.0,4000.);
498 loizides 1.5 for(Int_t i=0; i<5; i++) AddOutput(hDGenWBF[i]);
499 loizides 1.1
500 loizides 1.6 // bosons
501 loizides 1.1 sprintf(sb,"hDGenBosons_%d", 0); hDGenBosons[0] = new TH1D(sb,sb,10,-0.5,9.5);
502     sprintf(sb,"hDGenBosons_%d", 1); hDGenBosons[1] = new TH1D(sb,sb,200,0.0,400.0);
503     sprintf(sb,"hDGenBosons_%d", 2); hDGenBosons[2] = new TH1D(sb,sb,100,-5.0,5.0);
504     sprintf(sb,"hDGenBosons_%d", 3); hDGenBosons[3] = new TH1D(sb,sb,2000,0.0,2000.0);
505 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenBosons[i]);
506 ceballos 1.10
507     // photons
508     sprintf(sb,"hDGenPhotons_%d", 0); hDGenPhotons[0] = new TH1D(sb,sb,10,-0.5,9.5);
509     sprintf(sb,"hDGenPhotons_%d", 1); hDGenPhotons[1] = new TH1D(sb,sb,200,0.0,400.0);
510 ceballos 1.14 sprintf(sb,"hDGenPhotons_%d", 2); hDGenPhotons[2] = new TH1D(sb,sb,100,-5.0,5.0);
511 ceballos 1.10 for(Int_t i=0; i<3; i++) AddOutput(hDGenPhotons[i]);
512 ceballos 1.16
513     // auxiliar
514     sprintf(sb,"hDVMass_%d", 0); hDVMass[0] = new TH1D(sb,sb,200,0.,200.);
515     sprintf(sb,"hDVMass_%d", 1); hDVMass[1] = new TH1D(sb,sb,200,0.,200.);
516     sprintf(sb,"hDVMass_%d", 2); hDVMass[2] = new TH1D(sb,sb,200,0.,200.);
517     sprintf(sb,"hDVMass_%d", 3); hDVMass[3] = new TH1D(sb,sb,200,0.,200.);
518     sprintf(sb,"hDVMass_%d", 4); hDVMass[4] = new TH1D(sb,sb,200,0.,200.);
519     sprintf(sb,"hDVMass_%d", 5); hDVMass[5] = new TH1D(sb,sb,200,0.,200.);
520     sprintf(sb,"hDVMass_%d", 6); hDVMass[6] = new TH1D(sb,sb,200,0.,200.);
521     sprintf(sb,"hDVMass_%d", 7); hDVMass[7] = new TH1D(sb,sb,200,0.,200.);
522     sprintf(sb,"hDVMass_%d", 8); hDVMass[8] = new TH1D(sb,sb,200,0.,200.);
523     sprintf(sb,"hDVMass_%d", 9); hDVMass[9] = new TH1D(sb,sb,200,0.,200.);
524     for(Int_t i=0; i<10; i++) AddOutput(hDVMass[i]);
525 loizides 1.1 }
526     }