ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.24
Committed: Fri Jan 23 08:50:04 2009 UTC (16 years, 3 months ago) by loizides
Content type: text/plain
Branch: MAIN
Changes since 1.23: +5 -5 lines
Log Message:
Remove printout.

File Contents

# User Rev Content
1 loizides 1.24 // $Id: GeneratorMod.cc,v 1.23 2009/01/22 22:25:27 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     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.21 if (pm->PdgId() == 92) // string reached, terminate loop
76     break;
77 loizides 1.7 if (pm->Mother()->Is(MCParticle::kZ) || pm->Mother()->Is(MCParticle::kW)) {
78 loizides 1.5 GenLeptons->Add(p);
79     isGoodLepton = kTRUE;
80 loizides 1.7 break;
81 loizides 1.8 } else if (pm->Mother()->Is(MCParticle::kPi0) || pm->Mother()->Is(MCParticle::kEta)) {
82 loizides 1.7 // this is fake, but it is a trick to get rid of these cases and abort the loop
83     break;
84     }
85     pm = pm->Mother();
86 loizides 1.1 }
87 loizides 1.5 }
88    
89 loizides 1.7 // hadronic taus
90     else if (p->Is(MCParticle::kTau) && p->Status() == 2) {
91     if (!p->HasDaughter(MCParticle::kEl) && !p->HasDaughter(MCParticle::kMu)) {
92     const MCParticle *tv = p->FindDaughter(MCParticle::kTauNu);
93     if (tv) {
94     MCParticle *pm_f = new MCParticle(*p);
95     pm_f->SetMom(p->Px()-tv->Px(), p->Py()-tv->Py(),
96     p->Pz()-tv->Pz(), p->E()-tv->E());
97     GenTaus->AddOwned(pm_f);
98     } else {
99 loizides 1.18 SendError(kWarning, "Process", "Could not find a tau neutrino!");
100 loizides 1.5 }
101 loizides 1.1 }
102 loizides 1.5 }
103 loizides 1.1
104 loizides 1.5 // neutrinos
105 loizides 1.7 else if (p->Status() == 1 && p->IsNeutrino()) {
106 loizides 1.5 GenNeutrinos->Add(p);
107     }
108 loizides 1.1
109 loizides 1.5 // quarks from W/Z decays or top particles
110 loizides 1.7 else if (p->IsQuark() && p->HasMother()) {
111     if (p->Mother()->Is(MCParticle::kZ) || p->Mother()->Is(MCParticle::kW) ||
112     p->Is(MCParticle::kTop) || p->Mother()->Is(MCParticle::kTop)) {
113 loizides 1.5 GenQuarks->Add(p);
114 loizides 1.1 }
115 loizides 1.5 }
116 loizides 1.1
117 loizides 1.5 // qqH, information about the forward jets
118 loizides 1.8 else if (isqqH == kFALSE && p->Is(MCParticle::kH)) {
119 loizides 1.5 isqqH = kTRUE;
120 loizides 1.17 const MCParticle *pq1 = fParticles->At(i-1);
121     const MCParticle *pq2 = fParticles->At(i-2);
122 loizides 1.7 if (!pq1 || !pq2) {
123     SendError(kWarning, "Process", "Could not find quark pair!");
124 loizides 1.8 } else if (pq1->IsQuark() && pq2->IsQuark() &&
125     pq1->HasMother() && pq2->HasMother() &&
126 ceballos 1.23 pq1->Mother()->PdgId() == pq2->Mother()->PdgId()) {
127 loizides 1.5 GenqqHs->Add(pq1);
128     GenqqHs->Add(pq2);
129 ceballos 1.23 } else if (pq1->IsQuark() && pq2->IsQuark() &&
130     pq1->HasMother() && pq2->HasMother()) {
131 loizides 1.24 SendError(kWarning, "Process", "Mother q1(%d - %d) != mother q2 (%d - %d)!",
132     pq1->PdgId(),pq1->Mother()->PdgId(),
133     pq2->PdgId(),pq2->Mother()->PdgId());
134 loizides 1.1 }
135 ceballos 1.23
136 loizides 1.7 if (p->Status() == 3)
137 loizides 1.24 GenBosons->Add(p); // take higgs boson in account here rather in next else if
138 loizides 1.5 }
139 loizides 1.1
140 loizides 1.5 // information about bosons: W, Z, h, Z', W', H0, A0, H+
141 loizides 1.7 else if (p->Status() == 3 &&
142     (p->Is(MCParticle::kZ) || p->Is(MCParticle::kW) || p->Is(MCParticle::kH) ||
143     p->Is(MCParticle::kZp) || p->Is(MCParticle::kZpp) ||
144     p->Is(MCParticle::kH0) || p->Is(MCParticle::kA0) || p->Is(MCParticle::kHp))) {
145 loizides 1.5 GenBosons->Add(p);
146 loizides 1.1 }
147 ceballos 1.10
148 ceballos 1.14 // photons
149 ceballos 1.10 else if (p->Status() == 1 && p->Is(MCParticle::kGamma) &&
150 ceballos 1.14 p->Pt() > fPtPhotonMin && p->AbsEta() < fEtaPhotonMax) {
151 ceballos 1.10 GenPhotons->Add(p);
152     }
153    
154 loizides 1.18 // W/Z -> lnu for Madgraph
155     if (p->IsParton() && p->NDaughters() == 2) {
156 ceballos 1.16 CompositeParticle *diBoson = new CompositeParticle();
157     diBoson->AddDaughter(p->Daughter(0));
158     diBoson->AddDaughter(p->Daughter(1));
159 loizides 1.18 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
160 loizides 1.20 if (GetFillHist())
161 loizides 1.18 hDVMass[0]->Fill(TMath::Min(diBoson->Mass(),199.999));
162 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kMu));
163     }
164 loizides 1.18 else if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
165 loizides 1.20 if (GetFillHist())
166 loizides 1.18 hDVMass[1]->Fill(TMath::Min(diBoson->Mass(),199.999));
167 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kEl));
168     }
169 loizides 1.18 else if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
170 loizides 1.20 if (GetFillHist())
171 loizides 1.18 hDVMass[2]->Fill(TMath::Min(diBoson->Mass(),199.999));
172     const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
173     if (tau->HasDaughter(MCParticle::kMu))
174     GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
175     if (tau->HasDaughter(MCParticle::kEl))
176     GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
177     }
178     else if (p->Daughter(0)->Is(MCParticle::kMu) && p->Daughter(1)->Is(MCParticle::kMu)) {
179 loizides 1.20 if (GetFillHist())
180 loizides 1.18 hDVMass[3]->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::kEl) && p->Daughter(1)->Is(MCParticle::kEl)) {
185 loizides 1.20 if (GetFillHist())
186 loizides 1.18 hDVMass[4]->Fill(TMath::Min(diBoson->Mass(),199.999));
187 ceballos 1.16 GenLeptons->Add(p->Daughter(0));
188     GenLeptons->Add(p->Daughter(1));
189     }
190 loizides 1.18 else if (p->Daughter(0)->Is(MCParticle::kTau) && p->Daughter(1)->Is(MCParticle::kTau)) {
191 loizides 1.20 if (GetFillHist())
192 loizides 1.18 hDVMass[5]->Fill(TMath::Min(diBoson->Mass(),199.999));
193     const MCParticle *tau0 = p->Daughter(0);
194     if (tau0->HasDaughter(MCParticle::kMu))
195     GenLeptons->Add(tau0->FindDaughter(MCParticle::kMu));
196     if (tau0->HasDaughter(MCParticle::kEl))
197     GenLeptons->Add(tau0->FindDaughter(MCParticle::kEl));
198     const MCParticle *tau1 = p->Daughter(1);
199     if (tau1->HasDaughter(MCParticle::kMu))
200     GenLeptons->Add(tau1->FindDaughter(MCParticle::kMu));
201     if (tau1->HasDaughter(MCParticle::kEl))
202     GenLeptons->Add(tau1->FindDaughter(MCParticle::kEl));
203 ceballos 1.16 }
204     delete diBoson;
205     }
206    
207 loizides 1.18 // t -> lnu for Madgraph
208     if (p->Is(MCParticle::kTop)) {
209 ceballos 1.16 CompositeParticle *diBoson = new CompositeParticle();
210 loizides 1.18 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
211 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
212     diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
213 loizides 1.20 if (GetFillHist())
214 loizides 1.18 hDVMass[6]->Fill(TMath::Min(diBoson->Mass(),199.999));
215 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kMu));
216     }
217 loizides 1.18 else if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
218 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
219     diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
220 loizides 1.20 if (GetFillHist())
221 loizides 1.18 hDVMass[7]->Fill(TMath::Min(diBoson->Mass(),199.999));
222 ceballos 1.16 GenLeptons->Add(p->FindDaughter(MCParticle::kEl));
223     }
224 loizides 1.18 else if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
225 ceballos 1.16 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
226     diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
227 loizides 1.20 if (GetFillHist())
228 loizides 1.18 hDVMass[8]->Fill(TMath::Min(diBoson->Mass(),199.999));
229     const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
230     if (tau->HasDaughter(MCParticle::kMu))
231     GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
232     if (tau->HasDaughter(MCParticle::kEl))
233     GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
234     }
235     else if (!p->HasDaughter(MCParticle::kW)) {
236     for(UInt_t nd=0; nd<p->NDaughters(); ++nd)
237     if (p->Daughter(nd)->IsNot(MCParticle::kBottom) &&
238     p->Daughter(nd)->IsNot(MCParticle::kGamma))
239     diBoson->AddDaughter(p->Daughter(nd));
240 loizides 1.20 if (GetFillHist())
241 loizides 1.18 hDVMass[9]->Fill(TMath::Min(diBoson->Mass(),199.999));
242 ceballos 1.16 }
243     delete diBoson;
244     }
245 loizides 1.5 }
246 loizides 1.1
247 loizides 1.19 // sort according to pt
248     GenLeptons->Sort();
249     GenAllLeptons->Sort();
250     GenTaus->Sort();
251     GenNeutrinos->Sort();
252     GenQuarks->Sort();
253     GenqqHs->Sort();
254     GenBosons->Sort();
255     GenPhotons->Sort();
256    
257 loizides 1.8 // add objects to this event for other modules to use
258 loizides 1.13 AddObjThisEvt(GenLeptons);
259     AddObjThisEvt(GenAllLeptons);
260     AddObjThisEvt(GenTaus);
261     AddObjThisEvt(GenNeutrinos);
262     AddObjThisEvt(GenQuarks);
263     AddObjThisEvt(GenqqHs);
264     AddObjThisEvt(GenBosons);
265     AddObjThisEvt(GenPhotons);
266 ceballos 1.16
267 loizides 1.5 // fill histograms if requested
268 loizides 1.20 if (GetFillHist()) {
269 loizides 1.5
270 loizides 1.6 // leptons
271 loizides 1.1 hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
272 loizides 1.5 for(UInt_t i=0; i<GenLeptons->GetEntries(); i++) {
273 loizides 1.1 hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
274 loizides 1.5 hDGenLeptons[2]->Fill(TMath::Abs(GenLeptons->At(i)->Eta()));
275     hDGenLeptons[3]->Fill(GenLeptons->At(i)->PhiDeg());
276     for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++) {
277 loizides 1.1 CompositeParticle *dilepton = new CompositeParticle();
278     dilepton->AddDaughter(GenLeptons->At(i));
279     dilepton->AddDaughter(GenLeptons->At(j));
280     hDGenLeptons[4]->Fill(dilepton->Mass());
281     delete dilepton;
282     }
283 ceballos 1.22 }
284     // looking at events with two leptons
285     if (GenLeptons->GetEntries() == 2) {
286     hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(TMath::Abs(GenLeptons->At(0)->Eta()),
287     TMath::Abs(GenLeptons->At(1)->Eta())),
288 loizides 1.5 4.999));
289 ceballos 1.22 hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(TMath::Abs(GenLeptons->At(0)->Eta()),
290     TMath::Abs(GenLeptons->At(1)->Eta())),
291 loizides 1.5 4.999));
292 ceballos 1.22 if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
293     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5) {
294     hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
295     if (GenLeptons->At(0)->Pt() > 20.0) {
296     hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
297     if (GenLeptons->At(1)->Pt() > 10.0) {
298 loizides 1.1 CompositeParticle *dilepton = new CompositeParticle();
299 ceballos 1.22 dilepton->AddDaughter(GenLeptons->At(0));
300     dilepton->AddDaughter(GenLeptons->At(1));
301 loizides 1.1 hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
302 ceballos 1.22 if(dilepton->Mass() > 12.0){
303     hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(0)->Phi(),
304     GenLeptons->At(1)->Phi())
305     * 180./ TMath::Pi());
306     hDGenLeptons[11]->Fill(MathUtils::DeltaR(GenLeptons->At(0)->Eta(), GenLeptons->At(0)->Phi(),
307     GenLeptons->At(1)->Eta(), GenLeptons->At(1)->Phi()));
308     }
309 loizides 1.1 delete dilepton;
310     }
311     }
312     }
313     }
314 ceballos 1.22 // looking at events with three leptons
315     if (GenLeptons->GetEntries() == 3) {
316     if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
317     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5 &&
318     TMath::Abs(GenLeptons->At(2)->Eta()) < 2.5) {
319     hDGenLeptons[12]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
320     if (GenLeptons->At(0)->Pt() > 20.0) {
321     hDGenLeptons[13]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
322     hDGenLeptons[14]->Fill(TMath::Min(GenLeptons->At(2)->Pt(),199.999));
323     if (GenLeptons->At(1)->Pt() > 10.0 && GenLeptons->At(2)->Pt() > 10.0) {
324     CompositeParticle *dilepton01 = new CompositeParticle();
325     dilepton01->AddDaughter(GenLeptons->At(0));
326     dilepton01->AddDaughter(GenLeptons->At(1));
327     CompositeParticle *dilepton02 = new CompositeParticle();
328     dilepton02->AddDaughter(GenLeptons->At(0));
329     dilepton02->AddDaughter(GenLeptons->At(2));
330     CompositeParticle *dilepton12 = new CompositeParticle();
331     dilepton12->AddDaughter(GenLeptons->At(1));
332     dilepton12->AddDaughter(GenLeptons->At(2));
333     hDGenLeptons[15]->Fill(TMath::Min(dilepton01->Mass(),999.999));
334     hDGenLeptons[15]->Fill(TMath::Min(dilepton02->Mass(),999.999));
335     hDGenLeptons[15]->Fill(TMath::Min(dilepton12->Mass(),999.999));
336     CompositeParticle *trilepton = new CompositeParticle();
337     trilepton->AddDaughter(GenLeptons->At(0));
338     trilepton->AddDaughter(GenLeptons->At(1));
339     trilepton->AddDaughter(GenLeptons->At(2));
340     hDGenLeptons[16]->Fill(TMath::Min(trilepton->Mass(),999.999));
341     double deltaR[3] = {MathUtils::DeltaR(GenLeptons->At(0)->Eta(), GenLeptons->At(0)->Phi(),
342     GenLeptons->At(1)->Eta(), GenLeptons->At(1)->Phi()),
343     MathUtils::DeltaR(GenLeptons->At(0)->Eta(), GenLeptons->At(0)->Phi(),
344     GenLeptons->At(2)->Eta(), GenLeptons->At(2)->Phi()),
345     MathUtils::DeltaR(GenLeptons->At(1)->Eta(), GenLeptons->At(1)->Phi(),
346     GenLeptons->At(2)->Eta(), GenLeptons->At(2)->Phi())};
347     double deltaRMin = deltaR[0];
348     for(int i=1; i<3; i++) if(deltaRMin > deltaR[i]) deltaRMin = deltaR[i];
349     hDGenLeptons[17]->Fill(deltaRMin);
350    
351     delete dilepton01;
352     delete dilepton02;
353     delete dilepton12;
354     delete trilepton;
355     }
356     }
357     }
358     }
359     // looking at events with four leptons
360     if (GenLeptons->GetEntries() == 4) {
361     if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
362     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5 &&
363     TMath::Abs(GenLeptons->At(2)->Eta()) < 2.5 &&
364     TMath::Abs(GenLeptons->At(3)->Eta()) < 2.5) {
365     hDGenLeptons[18]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
366     if (GenLeptons->At(0)->Pt() > 20.0) {
367     hDGenLeptons[19]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
368     hDGenLeptons[20]->Fill(TMath::Min(GenLeptons->At(2)->Pt(),199.999));
369     hDGenLeptons[21]->Fill(TMath::Min(GenLeptons->At(3)->Pt(),199.999));
370     if (GenLeptons->At(1)->Pt() > 10.0 && GenLeptons->At(2)->Pt() > 10.0 &&
371     GenLeptons->At(3)->Pt() > 10.0) {
372     CompositeParticle *dilepton01 = new CompositeParticle();
373     dilepton01->AddDaughter(GenLeptons->At(0));
374     dilepton01->AddDaughter(GenLeptons->At(1));
375     CompositeParticle *dilepton02 = new CompositeParticle();
376     dilepton02->AddDaughter(GenLeptons->At(0));
377     dilepton02->AddDaughter(GenLeptons->At(2));
378     CompositeParticle *dilepton03 = new CompositeParticle();
379     dilepton03->AddDaughter(GenLeptons->At(0));
380     dilepton03->AddDaughter(GenLeptons->At(3));
381     CompositeParticle *dilepton12 = new CompositeParticle();
382     dilepton12->AddDaughter(GenLeptons->At(1));
383     dilepton12->AddDaughter(GenLeptons->At(2));
384     CompositeParticle *dilepton13 = new CompositeParticle();
385     dilepton13->AddDaughter(GenLeptons->At(1));
386     dilepton13->AddDaughter(GenLeptons->At(3));
387     CompositeParticle *dilepton23 = new CompositeParticle();
388     dilepton23->AddDaughter(GenLeptons->At(2));
389     dilepton23->AddDaughter(GenLeptons->At(3));
390     hDGenLeptons[22]->Fill(TMath::Min(dilepton01->Mass(),999.999));
391     hDGenLeptons[22]->Fill(TMath::Min(dilepton02->Mass(),999.999));
392     hDGenLeptons[22]->Fill(TMath::Min(dilepton03->Mass(),999.999));
393     hDGenLeptons[22]->Fill(TMath::Min(dilepton12->Mass(),999.999));
394     hDGenLeptons[22]->Fill(TMath::Min(dilepton13->Mass(),999.999));
395     hDGenLeptons[22]->Fill(TMath::Min(dilepton23->Mass(),999.999));
396     CompositeParticle *fourlepton = new CompositeParticle();
397     fourlepton->AddDaughter(GenLeptons->At(0));
398     fourlepton->AddDaughter(GenLeptons->At(1));
399     fourlepton->AddDaughter(GenLeptons->At(2));
400     fourlepton->AddDaughter(GenLeptons->At(3));
401     hDGenLeptons[23]->Fill(TMath::Min(fourlepton->Mass(),999.999));
402     double deltaR[6] = {MathUtils::DeltaR(GenLeptons->At(0)->Eta(), GenLeptons->At(0)->Phi(),
403     GenLeptons->At(1)->Eta(), GenLeptons->At(1)->Phi()),
404     MathUtils::DeltaR(GenLeptons->At(0)->Eta(), GenLeptons->At(0)->Phi(),
405     GenLeptons->At(2)->Eta(), GenLeptons->At(2)->Phi()),
406     MathUtils::DeltaR(GenLeptons->At(0)->Eta(), GenLeptons->At(0)->Phi(),
407     GenLeptons->At(3)->Eta(), GenLeptons->At(3)->Phi()),
408     MathUtils::DeltaR(GenLeptons->At(1)->Eta(), GenLeptons->At(1)->Phi(),
409     GenLeptons->At(2)->Eta(), GenLeptons->At(2)->Phi()),
410     MathUtils::DeltaR(GenLeptons->At(1)->Eta(), GenLeptons->At(1)->Phi(),
411     GenLeptons->At(3)->Eta(), GenLeptons->At(3)->Phi()),
412     MathUtils::DeltaR(GenLeptons->At(2)->Eta(), GenLeptons->At(2)->Phi(),
413     GenLeptons->At(3)->Eta(), GenLeptons->At(3)->Phi())};
414     double deltaRMin = deltaR[0];
415     for(int i=1; i<6; i++) if(deltaRMin > deltaR[i]) deltaRMin = deltaR[i];
416     hDGenLeptons[24]->Fill(deltaRMin);
417    
418     delete dilepton01;
419     delete dilepton02;
420     delete dilepton03;
421     delete dilepton12;
422     delete dilepton13;
423     delete dilepton23;
424     delete fourlepton;
425     }
426     }
427     }
428     }
429 loizides 1.1
430 loizides 1.6 // all leptons
431 ceballos 1.3 hDGenAllLeptons[0]->Fill(GenAllLeptons->GetEntries());
432 loizides 1.5 for(UInt_t i=0; i<GenAllLeptons->GetEntries(); i++) {
433 ceballos 1.3 hDGenAllLeptons[1]->Fill(GenAllLeptons->At(i)->Pt());
434     hDGenAllLeptons[2]->Fill(GenAllLeptons->At(i)->Eta());
435 loizides 1.5 hDGenAllLeptons[3]->Fill(GenAllLeptons->At(i)->PhiDeg());
436 ceballos 1.3 }
437    
438 loizides 1.6 // taus
439 loizides 1.1 hDGenTaus[0]->Fill(GenTaus->GetEntries());
440 loizides 1.5 for(UInt_t i=0; i<GenTaus->GetEntries(); i++) {
441 loizides 1.1 hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
442     hDGenTaus[2]->Fill(GenTaus->At(i)->Eta());
443 loizides 1.5 hDGenTaus[3]->Fill(GenTaus->At(i)->PhiDeg());
444 loizides 1.1 }
445    
446 loizides 1.6 // neutrinos
447 loizides 1.1 hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
448     CompositeParticle *neutrinoTotal = new CompositeParticle();
449 loizides 1.5 for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++) {
450     if (GenNeutrinos->At(i)->HasMother())
451 loizides 1.1 neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
452     }
453 loizides 1.5 if (GenNeutrinos->GetEntries() > 0) {
454 loizides 1.1 hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
455     hDGenNeutrinos[2]->Fill(neutrinoTotal->Eta());
456 loizides 1.5 hDGenNeutrinos[3]->Fill(neutrinoTotal->PhiDeg());
457 loizides 1.1 }
458     delete neutrinoTotal;
459    
460 loizides 1.6 // quarks
461 loizides 1.1 hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
462 loizides 1.5 for(UInt_t i=0; i<GenQuarks->GetEntries(); i++) {
463     for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++) {
464 loizides 1.1 CompositeParticle *dijet = new CompositeParticle();
465     dijet->AddDaughter(GenQuarks->At(i));
466     dijet->AddDaughter(GenQuarks->At(j));
467     hDGenQuarks[1]->Fill(dijet->Pt());
468     hDGenQuarks[2]->Fill(dijet->Mass());
469 loizides 1.5 if (TMath::Abs(GenQuarks->At(i)->Eta()) < 2.5 &&
470     TMath::Abs(GenQuarks->At(j)->Eta()) < 2.5) {
471 loizides 1.1 hDGenQuarks[3]->Fill(dijet->Pt());
472     hDGenQuarks[4]->Fill(dijet->Mass());
473     }
474     delete dijet;
475     }
476 ceballos 1.2 // b quark info
477 loizides 1.5 if (GenQuarks->At(i)->AbsPdgId() == 5) {
478 ceballos 1.2 hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
479     hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());
480     hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());
481 ceballos 1.22 if (GenLeptons->GetEntries() >= 2 &&
482     GenLeptons->At(0)->Pt() > 20 &&
483     GenLeptons->At(1)->Pt() > 15) {
484     if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
485     TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5) {
486 ceballos 1.2 hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());
487     hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());
488     hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());
489     }
490     }
491     }
492     // t quark info
493 loizides 1.5 else if (GenQuarks->At(i)->AbsPdgId() == 6) {
494 ceballos 1.2 hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
495     hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());
496     hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());
497     }
498     // light quark info
499     else {
500     hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
501     hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());
502     hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());
503     }
504 loizides 1.1 }
505    
506 loizides 1.6 // wbf
507 loizides 1.5 if (GenqqHs->GetEntries() == 2) {
508 loizides 1.1 hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
509     GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
510 loizides 1.5 hDGenWBF[1]->Fill(TMath::Abs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
511 loizides 1.1 hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
512     hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
513     CompositeParticle *diqq = new CompositeParticle();
514     diqq->AddDaughter(GenqqHs->At(0));
515     diqq->AddDaughter(GenqqHs->At(1));
516     hDGenWBF[4]->Fill(diqq->Mass());
517     delete diqq;
518     }
519    
520 loizides 1.6 // bosons
521 loizides 1.1 hDGenBosons[0]->Fill(GenBosons->GetEntries());
522 loizides 1.5 for(UInt_t i=0; i<GenBosons->GetEntries(); i++) {
523 loizides 1.1 hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
524     hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
525     hDGenBosons[3]->Fill(GenBosons->At(i)->Mass());
526     }
527 ceballos 1.10
528     // photons
529     hDGenPhotons[0]->Fill(GenPhotons->GetEntries());
530     for(UInt_t i=0; i<GenPhotons->GetEntries(); i++) {
531     hDGenPhotons[1]->Fill(GenPhotons->At(i)->Pt());
532     hDGenPhotons[2]->Fill(GenPhotons->At(i)->Eta());
533     }
534 loizides 1.1 }
535     }
536    
537     //--------------------------------------------------------------------------------------------------
538     void GeneratorMod::SlaveBegin()
539     {
540 loizides 1.5 // Book branch and histograms if wanted.
541    
542 loizides 1.1 ReqBranch(fMCPartName, fParticles);
543    
544 loizides 1.5 // fill histograms
545 loizides 1.20 if (GetFillHist()) {
546 loizides 1.5 char sb[1024];
547 loizides 1.6 // leptons from W
548 loizides 1.1 sprintf(sb,"hDGenLeptons_%d", 0); hDGenLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
549     sprintf(sb,"hDGenLeptons_%d", 1); hDGenLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
550     sprintf(sb,"hDGenLeptons_%d", 2); hDGenLeptons[2] = new TH1D(sb,sb,50,0.0,5.0);
551     sprintf(sb,"hDGenLeptons_%d", 3); hDGenLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
552     sprintf(sb,"hDGenLeptons_%d", 4); hDGenLeptons[4] = new TH1D(sb,sb,1000,0.0,1000.0);
553     sprintf(sb,"hDGenLeptons_%d", 5); hDGenLeptons[5] = new TH1D(sb,sb,50,0.0,5.0);
554     sprintf(sb,"hDGenLeptons_%d", 6); hDGenLeptons[6] = new TH1D(sb,sb,50,0.0,5.0);
555     sprintf(sb,"hDGenLeptons_%d", 7); hDGenLeptons[7] = new TH1D(sb,sb,100,0.0,200.0);
556     sprintf(sb,"hDGenLeptons_%d", 8); hDGenLeptons[8] = new TH1D(sb,sb,100,0.0,200.0);
557     sprintf(sb,"hDGenLeptons_%d", 9); hDGenLeptons[9] = new TH1D(sb,sb,1000,0.0,1000.0);
558     sprintf(sb,"hDGenLeptons_%d",10); hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0);
559 ceballos 1.22 sprintf(sb,"hDGenLeptons_%d",11); hDGenLeptons[11] = new TH1D(sb,sb,100,0.0,5.0);
560     sprintf(sb,"hDGenLeptons_%d",12); hDGenLeptons[12] = new TH1D(sb,sb,100,0.0,200.0);
561     sprintf(sb,"hDGenLeptons_%d",13); hDGenLeptons[13] = new TH1D(sb,sb,100,0.0,200.0);
562     sprintf(sb,"hDGenLeptons_%d",14); hDGenLeptons[14] = new TH1D(sb,sb,100,0.0,200.0);
563     sprintf(sb,"hDGenLeptons_%d",15); hDGenLeptons[15] = new TH1D(sb,sb,1000,0.0,1000.0);
564     sprintf(sb,"hDGenLeptons_%d",16); hDGenLeptons[16] = new TH1D(sb,sb,1000,0.0,1000.0);
565     sprintf(sb,"hDGenLeptons_%d",17); hDGenLeptons[17] = new TH1D(sb,sb,100,0.0,5.0);
566     sprintf(sb,"hDGenLeptons_%d",18); hDGenLeptons[18] = new TH1D(sb,sb,100,0.0,200.0);
567     sprintf(sb,"hDGenLeptons_%d",19); hDGenLeptons[19] = new TH1D(sb,sb,100,0.0,200.0);
568     sprintf(sb,"hDGenLeptons_%d",20); hDGenLeptons[20] = new TH1D(sb,sb,100,0.0,200.0);
569     sprintf(sb,"hDGenLeptons_%d",21); hDGenLeptons[21] = new TH1D(sb,sb,100,0.0,200.0);
570     sprintf(sb,"hDGenLeptons_%d",22); hDGenLeptons[22] = new TH1D(sb,sb,1000,0.0,1000.0);
571     sprintf(sb,"hDGenLeptons_%d",23); hDGenLeptons[23] = new TH1D(sb,sb,1000,0.0,1000.0);
572     sprintf(sb,"hDGenLeptons_%d",24); hDGenLeptons[24] = new TH1D(sb,sb,100,0.0,5.0);
573     for(Int_t i=0; i<25; i++) AddOutput(hDGenLeptons[i]);
574 loizides 1.1
575 loizides 1.6 // all leptons
576 ceballos 1.3 sprintf(sb,"hDGenAllLeptons_%d", 0); hDGenAllLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
577     sprintf(sb,"hDGenAllLeptons_%d", 1); hDGenAllLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
578 ceballos 1.14 sprintf(sb,"hDGenAllLeptons_%d", 2); hDGenAllLeptons[2] = new TH1D(sb,sb,100,-5.0,5.0);
579 ceballos 1.3 sprintf(sb,"hDGenAllLeptons_%d", 3); hDGenAllLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
580 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenAllLeptons[i]);
581 ceballos 1.3
582 loizides 1.6 // taus
583 loizides 1.1 sprintf(sb,"hDGenTaus_%d", 0); hDGenTaus[0] = new TH1D(sb,sb,10,-0.5,9.5);
584     sprintf(sb,"hDGenTaus_%d", 1); hDGenTaus[1] = new TH1D(sb,sb,100,0.0,200.0);
585     sprintf(sb,"hDGenTaus_%d", 2); hDGenTaus[2] = new TH1D(sb,sb,100,-5.0,5.0);
586     sprintf(sb,"hDGenTaus_%d", 3); hDGenTaus[3] = new TH1D(sb,sb,90,0.0,180.0);
587 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenTaus[i]);
588 loizides 1.1
589 loizides 1.6 // neutrinos
590 loizides 1.1 sprintf(sb,"hDGenNeutrinos_%d", 0); hDGenNeutrinos[0] = new TH1D(sb,sb,10,-0.5,9.5);
591     sprintf(sb,"hDGenNeutrinos_%d", 1); hDGenNeutrinos[1] = new TH1D(sb,sb,100,0.0,200.0);
592     sprintf(sb,"hDGenNeutrinos_%d", 2); hDGenNeutrinos[2] = new TH1D(sb,sb,100,-5.0,5.0);
593     sprintf(sb,"hDGenNeutrinos_%d", 3); hDGenNeutrinos[3] = new TH1D(sb,sb,90,0.0,180.0);
594 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);
595 loizides 1.1
596 loizides 1.6 // quarks
597 loizides 1.1 sprintf(sb,"hDGenQuarks_%d", 0); hDGenQuarks[0] = new TH1D(sb,sb,10,-0.5,9.5);
598     sprintf(sb,"hDGenQuarks_%d", 1); hDGenQuarks[1] = new TH1D(sb,sb,200,0.0,400.);
599     sprintf(sb,"hDGenQuarks_%d", 2); hDGenQuarks[2] = new TH1D(sb,sb,2000,0.0,2000.);
600     sprintf(sb,"hDGenQuarks_%d", 3); hDGenQuarks[3] = new TH1D(sb,sb,200,0.0,400.);
601     sprintf(sb,"hDGenQuarks_%d", 4); hDGenQuarks[4] = new TH1D(sb,sb,2000,0.0,2000.);
602 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 5); hDGenQuarks[5] = new TH1D(sb,sb,200,0.0,400.);
603     sprintf(sb,"hDGenQuarks_%d", 6); hDGenQuarks[6] = new TH1D(sb,sb,200,-10.0,10.);
604 loizides 1.5 sprintf(sb,"hDGenQuarks_%d", 7); hDGenQuarks[7] = new TH1D(sb,sb,200,-TMath::Pi(),
605     TMath::Pi());
606 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 8); hDGenQuarks[8] = new TH1D(sb,sb,200,0.0,400.);
607     sprintf(sb,"hDGenQuarks_%d", 9); hDGenQuarks[9] = new TH1D(sb,sb,200,-10.0,10.);
608 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",10); hDGenQuarks[10] = new TH1D(sb,sb,200,-TMath::Pi(),
609     TMath::Pi());
610 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d",11); hDGenQuarks[11] = new TH1D(sb,sb,200,0.0,400.);
611     sprintf(sb,"hDGenQuarks_%d",12); hDGenQuarks[12] = new TH1D(sb,sb,200,-10.0,10.);
612 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",13); hDGenQuarks[13] = new TH1D(sb,sb,200,-TMath::Pi(),
613     TMath::Pi());
614 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d",14); hDGenQuarks[14] = new TH1D(sb,sb,200,0.0,400.);
615     sprintf(sb,"hDGenQuarks_%d",15); hDGenQuarks[15] = new TH1D(sb,sb,200,-10.0,10.);
616 loizides 1.5 sprintf(sb,"hDGenQuarks_%d",16); hDGenQuarks[16] = new TH1D(sb,sb,200,-TMath::Pi(),
617     TMath::Pi());
618     for(Int_t i=0; i<17; i++) AddOutput(hDGenQuarks[i]);
619 loizides 1.1
620     // qqH
621     sprintf(sb,"hDGenWBF_%d", 0); hDGenWBF[0] = new TH1D(sb,sb,90,0.0,180.);
622     sprintf(sb,"hDGenWBF_%d", 1); hDGenWBF[1] = new TH1D(sb,sb,100,0.0,10.);
623     sprintf(sb,"hDGenWBF_%d", 2); hDGenWBF[2] = new TH1D(sb,sb,200,0.0,400.);
624     sprintf(sb,"hDGenWBF_%d", 3); hDGenWBF[3] = new TH1D(sb,sb,200,0.0,400.);
625     sprintf(sb,"hDGenWBF_%d", 4); hDGenWBF[4] = new TH1D(sb,sb,200,0.0,4000.);
626 loizides 1.5 for(Int_t i=0; i<5; i++) AddOutput(hDGenWBF[i]);
627 loizides 1.1
628 loizides 1.6 // bosons
629 loizides 1.1 sprintf(sb,"hDGenBosons_%d", 0); hDGenBosons[0] = new TH1D(sb,sb,10,-0.5,9.5);
630     sprintf(sb,"hDGenBosons_%d", 1); hDGenBosons[1] = new TH1D(sb,sb,200,0.0,400.0);
631     sprintf(sb,"hDGenBosons_%d", 2); hDGenBosons[2] = new TH1D(sb,sb,100,-5.0,5.0);
632     sprintf(sb,"hDGenBosons_%d", 3); hDGenBosons[3] = new TH1D(sb,sb,2000,0.0,2000.0);
633 loizides 1.5 for(Int_t i=0; i<4; i++) AddOutput(hDGenBosons[i]);
634 ceballos 1.10
635     // photons
636     sprintf(sb,"hDGenPhotons_%d", 0); hDGenPhotons[0] = new TH1D(sb,sb,10,-0.5,9.5);
637     sprintf(sb,"hDGenPhotons_%d", 1); hDGenPhotons[1] = new TH1D(sb,sb,200,0.0,400.0);
638 ceballos 1.14 sprintf(sb,"hDGenPhotons_%d", 2); hDGenPhotons[2] = new TH1D(sb,sb,100,-5.0,5.0);
639 ceballos 1.10 for(Int_t i=0; i<3; i++) AddOutput(hDGenPhotons[i]);
640 ceballos 1.16
641     // auxiliar
642     sprintf(sb,"hDVMass_%d", 0); hDVMass[0] = new TH1D(sb,sb,200,0.,200.);
643     sprintf(sb,"hDVMass_%d", 1); hDVMass[1] = new TH1D(sb,sb,200,0.,200.);
644     sprintf(sb,"hDVMass_%d", 2); hDVMass[2] = new TH1D(sb,sb,200,0.,200.);
645     sprintf(sb,"hDVMass_%d", 3); hDVMass[3] = new TH1D(sb,sb,200,0.,200.);
646     sprintf(sb,"hDVMass_%d", 4); hDVMass[4] = new TH1D(sb,sb,200,0.,200.);
647     sprintf(sb,"hDVMass_%d", 5); hDVMass[5] = new TH1D(sb,sb,200,0.,200.);
648     sprintf(sb,"hDVMass_%d", 6); hDVMass[6] = new TH1D(sb,sb,200,0.,200.);
649     sprintf(sb,"hDVMass_%d", 7); hDVMass[7] = new TH1D(sb,sb,200,0.,200.);
650     sprintf(sb,"hDVMass_%d", 8); hDVMass[8] = new TH1D(sb,sb,200,0.,200.);
651     sprintf(sb,"hDVMass_%d", 9); hDVMass[9] = new TH1D(sb,sb,200,0.,200.);
652     for(Int_t i=0; i<10; i++) AddOutput(hDVMass[i]);
653 loizides 1.1 }
654     }