ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.62
Committed: Mon Sep 27 15:49:24 2010 UTC (14 years, 7 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_015b, Mit_015a, Mit_015
Changes since 1.61: +14 -6 lines
Log Message:
small fix

File Contents

# Content
1 // $Id: GeneratorMod.cc,v 1.61 2010/09/27 12:42:20 ceballos Exp $
2
3 #include "MitPhysics/Mods/interface/GeneratorMod.h"
4 #include "MitCommon/MathTools/interface/MathUtils.h"
5 #include "MitAna/DataTree/interface/MetCol.h"
6 #include "MitAna/DataTree/interface/MCParticleCol.h"
7 #include "MitPhysics/Init/interface/ModNames.h"
8 #include <TH1D.h>
9 #include <TH2D.h>
10 #include <TParameter.h>
11
12 using namespace mithep;
13
14 ClassImp(mithep::GeneratorMod)
15
16 //--------------------------------------------------------------------------------------------------
17 GeneratorMod::GeneratorMod(const char *name, const char *title) :
18 BaseMod(name,title),
19 fIsData(kFALSE),
20 fPrintDebug(kFALSE),
21 fCopyArrays(kFALSE),
22 fMCPartName(Names::gkMCPartBrn),
23 fMCMETName(ModNames::gkMCMETName),
24 fMCLeptonsName(ModNames::gkMCLeptonsName),
25 fMCAllLeptonsName(ModNames::gkMCAllLeptonsName),
26 fMCTausName(ModNames::gkMCTausName),
27 fMCNeutrinosName(ModNames::gkMCNeutrinosName),
28 fMCQuarksName(ModNames::gkMCQuarksName),
29 fMCqqHsName(ModNames::gkMCqqHsName),
30 fMCBosonsName(ModNames::gkMCBosonsName),
31 fMCPhotonsName(ModNames::gkMCPhotonsName),
32 fMCRadPhotonsName(ModNames::gkMCRadPhotonsName),
33 fMCISRPhotonsName(ModNames::gkMCISRPhotonsName),
34 fPtLeptonMin(0.0),
35 fEtaLeptonMax(5.0),
36 fPtPhotonMin(0.0),
37 fEtaPhotonMax(5.0),
38 fPtRadPhotonMin(0.0),
39 fEtaRadPhotonMax(5.0),
40 fPdgIdCut(0),
41 fMassMinCut(-FLT_MAX),
42 fMassMaxCut(FLT_MAX),
43 fApplyISRFilter(kFALSE),
44 fParticles(0)
45 {
46 // Constructor
47 fGenLeptons = new MCParticleArr();
48 fGenLeptons->SetName(TString("Pub") + ModNames::gkMCLeptonsName);
49 fGenAllLeptons = new MCParticleArr();
50 fGenAllLeptons->SetName(TString("Pub") + ModNames::gkMCAllLeptonsName);
51 fGenTaus = new MCParticleArr();
52 fGenTaus->SetName(TString("Pub") + ModNames::gkMCTausName);
53 fGenNeutrinos = new MCParticleArr();
54 fGenNeutrinos->SetName(TString("Pub") + ModNames::gkMCNeutrinosName);
55 fGenQuarks = new MCParticleArr();
56 fGenQuarks->SetName(TString("Pub") + ModNames::gkMCQuarksName);
57 fGenqqHs = new MCParticleArr();
58 fGenqqHs->SetName(TString("Pub") + ModNames::gkMCqqHsName);
59 fGenBosons = new MCParticleArr();
60 fGenBosons->SetName(TString("Pub") + ModNames::gkMCBosonsName);
61 fGenPhotons = new MCParticleArr();
62 fGenPhotons->SetName(TString("Pub") + ModNames::gkMCPhotonsName);
63 fGenRadPhotons = new MCParticleArr();
64 fGenRadPhotons->SetName(TString("Pub") + ModNames::gkMCRadPhotonsName);
65 fGenISRPhotons = new MCParticleArr();
66 fGenISRPhotons->SetName(TString("Pub") + ModNames::gkMCISRPhotonsName);
67 }
68
69
70 //--------------------------------------------------------------------------------------------------
71 GeneratorMod::~GeneratorMod()
72 {
73 // Destructor.
74
75 delete fGenLeptons;
76 delete fGenAllLeptons;
77 delete fGenTaus;
78 delete fGenNeutrinos;
79 delete fGenQuarks;
80 delete fGenqqHs;
81 delete fGenBosons;
82 delete fGenPhotons;
83 delete fGenRadPhotons;
84 delete fGenISRPhotons;
85 }
86
87 //--------------------------------------------------------------------------------------------------
88 void GeneratorMod::Process()
89 {
90 // Process entries of the tree.
91
92 // these arrays will be filled in the loop of particles
93 MetOArr *GenMet = new MetOArr;
94 GenMet->SetName(fMCMETName);
95 GenMet->SetOwner(kTRUE);
96 MCParticleOArr *GenLeptons = new MCParticleOArr;
97 GenLeptons->SetName(fMCLeptonsName);
98 MCParticleOArr *GenAllLeptons = new MCParticleOArr;
99 GenAllLeptons->SetName(fMCAllLeptonsName);
100 MCParticleOArr *GenTaus = new MCParticleOArr;
101 GenTaus->SetName(fMCTausName);
102 GenTaus->SetOwner(kTRUE);
103 MCParticleOArr *GenNeutrinos = new MCParticleOArr;
104 GenNeutrinos->SetName(fMCNeutrinosName);
105 MCParticleOArr *GenQuarks = new MCParticleOArr;
106 GenQuarks->SetName(fMCQuarksName);
107 MCParticleOArr *GenqqHs = new MCParticleOArr;
108 GenqqHs->SetName(fMCqqHsName);
109 MCParticleOArr *GenBosons = new MCParticleOArr;
110 GenBosons->SetName(fMCBosonsName);
111 MCParticleOArr *GenPhotons = new MCParticleOArr;
112 GenPhotons->SetName(fMCPhotonsName);
113 MCParticleOArr *GenRadPhotons = new MCParticleOArr;
114 GenRadPhotons->SetName(fMCRadPhotonsName);
115 MCParticleOArr *GenISRPhotons = new MCParticleOArr;
116 GenISRPhotons->SetName(fMCISRPhotonsName);
117
118 MCParticleOArr *GenTempMG0 = new MCParticleOArr;
119
120 Bool_t isOld = kFALSE;
121 Int_t sumV[2] = {0, 0}; // W, Z
122 Int_t sumVVFlavor[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
123 Double_t totalMET[3] = {0.0, 0.0, 0.0};
124 Bool_t isqqH = kFALSE;
125 Double_t ptMin = 999999.;
126 if(fIsData == kFALSE){
127 if (fPrintDebug)
128 printf("\n************ Next Event ************\n\n");
129
130 // load MCParticle branch
131 LoadEventObject(fMCPartName, fParticles);
132
133 for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
134 const MCParticle *p = fParticles->At(i);
135
136 if(fPrintDebug)
137 p->Print("l");
138
139 // rad photons, includes gamma from WWGamma vertex.
140 if( p->Is(MCParticle::kGamma) && p->Status() == 1 &&
141 p->Pt() > fPtRadPhotonMin && p->AbsEta() < fEtaRadPhotonMax &&
142 p->DistinctMother() && p->DistinctMother()->Status() == 3 &&
143 (p->DistinctMother()->Is(MCParticle::kEl) || p->DistinctMother()->Is(MCParticle::kMu) ||
144 p->DistinctMother()->Is(MCParticle::kTau) || p->DistinctMother()->Is(MCParticle::kW))
145 ) {
146 CompositeParticle *object = new CompositeParticle();
147 object->AddDaughter(p);
148 object->AddDaughter(p->DistinctMother());
149 if(object->Mass() > 1.0 || p->DistinctMother()->Is(MCParticle::kW)) GenRadPhotons->Add(p);
150 delete object;
151 }
152
153 // ISR photons
154 if( p->Is(MCParticle::kGamma) && p->Status() == 1 &&
155 p->Pt() > fPtRadPhotonMin && p->AbsEta() < fEtaRadPhotonMax &&
156 p->DistinctMother() && p->DistinctMother()->IsParton()
157 ) {
158 GenISRPhotons->Add(p);
159 }
160
161 // MET computation at generation level
162 if (p->Status() == 1 && !p->IsNeutrino()) {
163 totalMET[0] = totalMET[0] + p->Px();
164 totalMET[1] = totalMET[1] + p->Py();
165 totalMET[2] = totalMET[2] + p->Pz();
166 }
167
168 if (!p->IsGenerated()) continue;
169
170 if ((((p->Is(MCParticle::kEl) || p->Is(MCParticle::kMu)) &&
171 !p->HasMother(MCParticle::kTau, kFALSE)) || p->Is(MCParticle::kTau)) &&
172 p->Status() == 3 &&
173 (p->HasMother(MCParticle::kW, kFALSE) || p->HasMother(MCParticle::kZ, kFALSE))) {
174 if(p->Pt() < ptMin) ptMin = p->Pt();
175 }
176
177 // all muons/electrons
178 if ((p->Is(MCParticle::kEl) || p->Is(MCParticle::kMu)) && p->Status() == 1) {
179 if (p->Pt() > fPtLeptonMin && p->AbsEta() < fEtaLeptonMax) {
180 GenAllLeptons->Add(p);
181 }
182 Bool_t isGoodLepton = kFALSE;
183 const MCParticle *pm = p;
184 while (pm->HasMother() && isGoodLepton == kFALSE) {
185 if (pm->PdgId() == 92) // string reached, terminate loop
186 break;
187 if (pm->Mother()->Is(MCParticle::kZ) || pm->Mother()->Is(MCParticle::kW) ||
188 pm->Mother()->Is(MCParticle::kZp) || pm->Mother()->Is(MCParticle::kWp) ||
189 pm->Mother()->Is(MCParticle::kH) || pm->Mother()->Is(MCParticle::kH0) ||
190 pm->Mother()->Is(MCParticle::kA0) || pm->Mother()->Is(MCParticle::kHp)) {
191 GenLeptons->Add(p);
192 isGoodLepton = kTRUE;
193 break;
194 } else if (pm->Mother()->Is(MCParticle::kPi0) || pm->Mother()->Is(MCParticle::kEta)) {
195 // this is fake, but it is a trick to get rid of these cases and abort the loop
196 break;
197 }
198 pm = pm->Mother();
199 }
200 }
201
202 // hadronic taus
203 else if (p->Is(MCParticle::kTau) && p->Status() == 2) {
204 if (!p->HasDaughter(MCParticle::kEl) && !p->HasDaughter(MCParticle::kMu)) {
205 const MCParticle *tv = p->FindDaughter(MCParticle::kTauNu);
206 if (tv) {
207 MCParticle *pm_f = new MCParticle(*p);
208 pm_f->SetMom(p->Px()-tv->Px(), p->Py()-tv->Py(),
209 p->Pz()-tv->Pz(), p->E()-tv->E());
210 GenTaus->AddOwned(pm_f);
211 } else {
212 //SendError(kWarning, "Process", "Could not find a tau neutrino!");
213 }
214 }
215 }
216
217 // neutrinos
218 else if (p->Status() == 1 && p->IsNeutrino()) {
219 GenNeutrinos->Add(p);
220 }
221
222 // quarks from W/Z decays or top particles
223 else if (p->IsQuark() && p->HasMother()) {
224 if (p->Mother()->Is(MCParticle::kZ) || p->Mother()->Is(MCParticle::kW) ||
225 p->Is(MCParticle::kTop) || p->Mother()->Is(MCParticle::kTop)) {
226 GenQuarks->Add(p);
227 }
228 }
229
230 // qqH, information about the forward jets
231 else if (isqqH == kFALSE && p->Is(MCParticle::kH)) {
232 isqqH = kTRUE;
233 const MCParticle *pq1 = fParticles->At(i-1);
234 const MCParticle *pq2 = fParticles->At(i-2);
235 if (!pq1 || !pq2) {
236 SendError(kWarning, "Process", "Could not find quark pair!");
237 } else if (pq1->IsQuark() && pq2->IsQuark() &&
238 pq1->HasMother() && pq2->HasMother() &&
239 pq1->Mother() == pq2->Mother()) {
240 GenqqHs->Add(pq1);
241 GenqqHs->Add(pq2);
242 }
243
244 if (p->Status() == 3)
245 GenBosons->Add(p); // take higgs boson in account here rather in next else if
246 }
247
248 // information about bosons: W, Z, h, Z', W', H0, A0, H+
249 else if ((p->Status() == 3 &&
250 (p->Is(MCParticle::kZ) || p->Is(MCParticle::kW) || p->Is(MCParticle::kH) ||
251 p->Is(MCParticle::kZp) || p->Is(MCParticle::kZpp) ||
252 p->Is(MCParticle::kH0) || p->Is(MCParticle::kA0) || p->Is(MCParticle::kHp))) ||
253 (p->Status() == 2 &&
254 (p->Is(MCParticle::kJPsi) || p->Is(MCParticle::kUpsilon)))) {
255 GenBosons->Add(p);
256 if (p->Is(MCParticle::kW)) sumV[0]++;
257 else if(p->Is(MCParticle::kZ)) sumV[1]++;
258 if (p->Is(MCParticle::kW) && p->HasDaughter(MCParticle::kMu) &&
259 p->HasDaughter(MCParticle::kMuNu))
260 sumVVFlavor[0]++;
261 else if(p->Is(MCParticle::kW) && p->HasDaughter(MCParticle::kEl) &&
262 p->HasDaughter(MCParticle::kElNu))
263 sumVVFlavor[1]++;
264 else if(p->Is(MCParticle::kW) && p->HasDaughter(MCParticle::kTau) &&
265 p->HasDaughter(MCParticle::kTauNu))
266 sumVVFlavor[2]++;
267 else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kMu,kTRUE) &&
268 p->HasDaughter(-1*MCParticle::kMu,kTRUE))
269 sumVVFlavor[3]++;
270 else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kEl,kTRUE) &&
271 p->HasDaughter(-1*MCParticle::kEl,kTRUE))
272 sumVVFlavor[4]++;
273 else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kTau,kTRUE) &&
274 p->HasDaughter(-1*MCParticle::kTau,kTRUE))
275 sumVVFlavor[5]++;
276 else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kMuNu,kTRUE) &&
277 p->HasDaughter(-1*MCParticle::kMuNu,kTRUE))
278 sumVVFlavor[6]++;
279 else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kElNu,kTRUE) &&
280 p->HasDaughter(-1*MCParticle::kElNu,kTRUE))
281 sumVVFlavor[7]++;
282 else if(p->Is(MCParticle::kZ) && p->HasDaughter(MCParticle::kTauNu,kTRUE) &&
283 p->HasDaughter(-1*MCParticle::kTauNu,kTRUE))
284 sumVVFlavor[8]++;
285 }
286
287 // photons
288 else if (p->Status() == 1 && p->Is(MCParticle::kGamma) &&
289 p->Pt() > fPtPhotonMin && p->AbsEta() < fEtaPhotonMax) {
290 GenPhotons->Add(p);
291 }
292
293 // W/Z -> lnu for Madgraph
294 if (p->IsParton() && p->NDaughters() >= 2) {
295 CompositeParticle *diBoson = new CompositeParticle();
296 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
297 isOld = kFALSE;
298 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
299 if(p->FindDaughter(MCParticle::kMu) == GenTempMG0->At(nl)) {
300 isOld = kTRUE;
301 break;
302 }
303 }
304 if(isOld == kFALSE){
305 GenTempMG0->Add(p->FindDaughter(MCParticle::kMu));
306 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
307 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
308 sumV[0]++;
309 sumVVFlavor[0]++;
310 if (GetFillHist())
311 hDVMass[0]->Fill(TMath::Min(diBoson->Mass(),199.999));
312 const MCParticle *tmp_mu = p->FindDaughter(MCParticle::kMu);
313 while (tmp_mu->HasDaughter(MCParticle::kMu) &&
314 tmp_mu->FindDaughter(MCParticle::kMu)->IsGenerated())
315 tmp_mu = tmp_mu->FindDaughter(MCParticle::kMu);
316
317 GenLeptons->Add(tmp_mu);
318 }
319 }
320 if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
321 isOld = kFALSE;
322 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
323 if(p->FindDaughter(MCParticle::kEl) == GenTempMG0->At(nl)) {
324 isOld = kTRUE;
325 break;
326 }
327 }
328 if(isOld == kFALSE){
329 GenTempMG0->Add(p->FindDaughter(MCParticle::kEl));
330 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
331 diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
332 sumV[0]++;
333 sumVVFlavor[1]++;
334 if (GetFillHist())
335 hDVMass[1]->Fill(TMath::Min(diBoson->Mass(),199.999));
336 const MCParticle *tmp_e = p->FindDaughter(MCParticle::kEl);
337 while (tmp_e->HasDaughter(MCParticle::kEl) &&
338 tmp_e->FindDaughter(MCParticle::kEl)->IsGenerated())
339 tmp_e = tmp_e->FindDaughter(MCParticle::kEl);
340 GenLeptons->Add(tmp_e);
341 }
342 }
343 if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
344 isOld = kFALSE;
345 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
346 if(p->FindDaughter(MCParticle::kTau) == GenTempMG0->At(nl)) {
347 isOld = kTRUE;
348 break;
349 }
350 }
351 if(isOld == kFALSE){
352 GenTempMG0->Add(p->FindDaughter(MCParticle::kTau));
353 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
354 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
355 sumV[0]++;
356 sumVVFlavor[2]++;
357 if (GetFillHist())
358 hDVMass[2]->Fill(TMath::Min(diBoson->Mass(),199.999));
359 const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
360 if (tau->HasDaughter(MCParticle::kMu))
361 GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
362 if (tau->HasDaughter(MCParticle::kEl))
363 GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
364 if (tau->HasDaughter(MCParticle::kTau)) {
365 const MCParticle *tau_second = tau->FindDaughter(MCParticle::kTau);
366 if (tau_second->HasDaughter(MCParticle::kMu))
367 GenLeptons->Add(tau_second->FindDaughter(MCParticle::kMu));
368 if (tau_second->HasDaughter(MCParticle::kEl))
369 GenLeptons->Add(tau_second->FindDaughter(MCParticle::kEl));
370 }
371 }
372 }
373 if (p->HasDaughter(MCParticle::kMu,kTRUE) && p->HasDaughter(-1*MCParticle::kMu,kTRUE)) {
374 isOld = kFALSE;
375 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
376 if(p->FindDaughter(MCParticle::kMu,kTRUE) == GenTempMG0->At(nl)) {
377 isOld = kTRUE;
378 break;
379 }
380 }
381 if(isOld == kFALSE){
382 GenTempMG0->Add(p->FindDaughter(MCParticle::kMu,kTRUE));
383 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu,kTRUE));
384 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kMu,kTRUE));
385 sumV[1]++;
386 sumVVFlavor[3]++;
387 if (GetFillHist())
388 hDVMass[3]->Fill(TMath::Min(diBoson->Mass(),199.999));
389 const MCParticle *tmp_mu0 = p->FindDaughter(MCParticle::kMu,kTRUE);
390 while (tmp_mu0->HasDaughter(MCParticle::kMu) &&
391 tmp_mu0->FindDaughter(MCParticle::kMu)->IsGenerated())
392 tmp_mu0 = tmp_mu0->FindDaughter(MCParticle::kMu);
393 const MCParticle *tmp_mu1 = p->FindDaughter(-1*MCParticle::kMu,kTRUE);
394 while (tmp_mu1->HasDaughter(MCParticle::kMu) &&
395 tmp_mu1->FindDaughter(MCParticle::kMu)->IsGenerated())
396 tmp_mu1 = tmp_mu1->FindDaughter(MCParticle::kMu);
397 GenLeptons->Add(tmp_mu0);
398 GenLeptons->Add(tmp_mu1);
399 }
400 }
401 if (p->HasDaughter(MCParticle::kEl,kTRUE) && p->HasDaughter(-1*MCParticle::kEl,kTRUE)) {
402 isOld = kFALSE;
403 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
404 if(p->FindDaughter(MCParticle::kEl,kTRUE) == GenTempMG0->At(nl)) {
405 isOld = kTRUE;
406 break;
407 }
408 }
409 if(isOld == kFALSE){
410 GenTempMG0->Add(p->FindDaughter(MCParticle::kEl,kTRUE));
411 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl,kTRUE));
412 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kEl,kTRUE));
413 sumV[1]++;
414 sumVVFlavor[4]++;
415 if (GetFillHist())
416 hDVMass[4]->Fill(TMath::Min(diBoson->Mass(),199.999));
417 const MCParticle *tmp_e0 = p->Daughter(0);
418 while (tmp_e0->HasDaughter(MCParticle::kEl) &&
419 tmp_e0->FindDaughter(MCParticle::kEl)->IsGenerated())
420 tmp_e0 = tmp_e0->FindDaughter(MCParticle::kEl);
421 const MCParticle *tmp_e1 = p->Daughter(1);
422 while (tmp_e1->HasDaughter(MCParticle::kEl) &&
423 tmp_e1->FindDaughter(MCParticle::kEl)->IsGenerated())
424 tmp_e1 = tmp_e1->FindDaughter(MCParticle::kEl);
425 GenLeptons->Add(tmp_e0);
426 GenLeptons->Add(tmp_e1);
427 }
428 }
429 if (p->HasDaughter(MCParticle::kTau,kTRUE) && p->HasDaughter(-1*MCParticle::kTau,kTRUE)) {
430 isOld = kFALSE;
431 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
432 if(p->FindDaughter(MCParticle::kTau,kTRUE) == GenTempMG0->At(nl)) {
433 isOld = kTRUE;
434 break;
435 }
436 }
437 if(isOld == kFALSE){
438 GenTempMG0->Add(p->FindDaughter(MCParticle::kTau,kTRUE));
439 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau,kTRUE));
440 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kTau,kTRUE));
441 sumV[1]++;
442 sumVVFlavor[5]++;
443 if (GetFillHist())
444 hDVMass[5]->Fill(TMath::Min(diBoson->Mass(),199.999));
445 const MCParticle *tau0 = p->Daughter(0);
446 if (tau0->HasDaughter(MCParticle::kMu))
447 GenLeptons->Add(tau0->FindDaughter(MCParticle::kMu));
448 if (tau0->HasDaughter(MCParticle::kEl))
449 GenLeptons->Add(tau0->FindDaughter(MCParticle::kEl));
450 const MCParticle *tau1 = p->Daughter(1);
451 if (tau1->HasDaughter(MCParticle::kMu))
452 GenLeptons->Add(tau1->FindDaughter(MCParticle::kMu));
453 if (tau1->HasDaughter(MCParticle::kEl))
454 GenLeptons->Add(tau1->FindDaughter(MCParticle::kEl));
455 if (tau0->HasDaughter(MCParticle::kTau)) {
456 const MCParticle *tau0_second = tau0->FindDaughter(MCParticle::kTau);
457 if (tau0_second->HasDaughter(MCParticle::kMu))
458 GenLeptons->Add(tau0_second->FindDaughter(MCParticle::kMu));
459 if (tau0_second->HasDaughter(MCParticle::kEl))
460 GenLeptons->Add(tau0_second->FindDaughter(MCParticle::kEl));
461 }
462 if (tau1->HasDaughter(MCParticle::kTau)) {
463 const MCParticle *tau1_second = tau1->FindDaughter(MCParticle::kTau);
464 if (tau1_second->HasDaughter(MCParticle::kMu))
465 GenLeptons->Add(tau1_second->FindDaughter(MCParticle::kMu));
466 if (tau1_second->HasDaughter(MCParticle::kEl))
467 GenLeptons->Add(tau1_second->FindDaughter(MCParticle::kEl));
468 }
469 }
470 }
471 if (p->HasDaughter(MCParticle::kMuNu,kTRUE) && p->HasDaughter(-1*MCParticle::kMuNu,kTRUE)) {
472 isOld = kFALSE;
473 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
474 if(p->FindDaughter(MCParticle::kMuNu,kTRUE) == GenTempMG0->At(nl)) {
475 isOld = kTRUE;
476 break;
477 }
478 }
479 if(isOld == kFALSE){
480 GenTempMG0->Add(p->FindDaughter(MCParticle::kMuNu,kTRUE));
481 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu,kTRUE));
482 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kMuNu,kTRUE));
483 sumV[1]++;
484 sumVVFlavor[6]++;
485 if (GetFillHist())
486 hDVMass[6]->Fill(TMath::Min(diBoson->Mass(),199.999));
487 }
488 }
489 if (p->HasDaughter(MCParticle::kElNu,kTRUE) && p->HasDaughter(-1*MCParticle::kElNu,kTRUE)) {
490 isOld = kFALSE;
491 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
492 if(p->FindDaughter(MCParticle::kElNu,kTRUE) == GenTempMG0->At(nl)) {
493 isOld = kTRUE;
494 break;
495 }
496 }
497 if(isOld == kFALSE){
498 GenTempMG0->Add(p->FindDaughter(MCParticle::kElNu,kTRUE));
499 diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu,kTRUE));
500 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kElNu,kTRUE));
501 sumV[1]++;
502 sumVVFlavor[7]++;
503 if (GetFillHist())
504 hDVMass[7]->Fill(TMath::Min(diBoson->Mass(),199.999));
505 }
506 }
507 if (p->HasDaughter(MCParticle::kTauNu,kTRUE) && p->HasDaughter(-1*MCParticle::kTauNu,kTRUE)) {
508 isOld = kFALSE;
509 for(UInt_t nl = 0; nl < GenTempMG0->GetEntries(); nl++){
510 if(p->FindDaughter(MCParticle::kTauNu,kTRUE) == GenTempMG0->At(nl)) {
511 isOld = kTRUE;
512 break;
513 }
514 }
515 if(isOld == kFALSE){
516 GenTempMG0->Add(p->FindDaughter(MCParticle::kTauNu,kTRUE));
517 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu,kTRUE));
518 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kTauNu,kTRUE));
519 sumV[1]++;
520 sumVVFlavor[8]++;
521 if (GetFillHist())
522 hDVMass[8]->Fill(TMath::Min(diBoson->Mass(),199.999));
523 }
524 }
525 delete diBoson;
526 }
527
528 // t -> lnu for Madgraph
529 if (p->Is(MCParticle::kTop)) {
530 CompositeParticle *diBoson = new CompositeParticle();
531 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
532 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
533 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
534 if (GetFillHist())
535 hDVMass[9]->Fill(TMath::Min(diBoson->Mass(),199.999));
536 GenLeptons->Add(p->FindDaughter(MCParticle::kMu));
537 }
538 else if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
539 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
540 diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
541 if (GetFillHist())
542 hDVMass[10]->Fill(TMath::Min(diBoson->Mass(),199.999));
543 GenLeptons->Add(p->FindDaughter(MCParticle::kEl));
544 }
545 else if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
546 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
547 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
548 if (GetFillHist())
549 hDVMass[11]->Fill(TMath::Min(diBoson->Mass(),199.999));
550 const MCParticle *tau = p->FindDaughter(MCParticle::kTau);
551 if (tau->HasDaughter(MCParticle::kMu))
552 GenLeptons->Add(tau->FindDaughter(MCParticle::kMu));
553 if (tau->HasDaughter(MCParticle::kEl))
554 GenLeptons->Add(tau->FindDaughter(MCParticle::kEl));
555 }
556 else if (!p->HasDaughter(MCParticle::kW)) {
557 for(UInt_t nd=0; nd<p->NDaughters(); ++nd)
558 if (p->Daughter(nd)->IsNot(MCParticle::kBottom) &&
559 p->Daughter(nd)->IsNot(MCParticle::kGamma))
560 diBoson->AddDaughter(p->Daughter(nd));
561 if (GetFillHist())
562 hDVMass[12]->Fill(TMath::Min(diBoson->Mass(),199.999));
563 }
564 delete diBoson;
565 }
566
567 // mass cut for given pid
568 if(fPdgIdCut && p->Is(fPdgIdCut) &&
569 (p->Mass() < fMassMinCut || p->Mass() > fMassMaxCut)) {
570 SkipEvent();
571 return;
572 }
573 } // end loop of particles
574
575 delete GenTempMG0;
576
577 // --------------------------------
578 // Begin special study about VVjets
579 // --------------------------------
580 if(sumV[0] + 4*sumV[1] == 2 || sumV[0] + 4*sumV[1] == 5 || sumV[0] + 4*sumV[1] == 8){
581 MCParticleOArr *GenTempMG1 = new MCParticleOArr;
582 Double_t diBosonMass[2] = {0., 0.};
583 for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
584 const MCParticle *p = fParticles->At(i);
585
586 if (p->IsParton() && p->NDaughters() >= 2) {
587 CompositeParticle *diBoson = new CompositeParticle();
588 if (p->HasDaughter(MCParticle::kMu) && p->HasDaughter(MCParticle::kMuNu)) {
589 isOld = kFALSE;
590 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
591 if(p->FindDaughter(MCParticle::kMu) == GenTempMG1->At(nl)) {
592 isOld = kTRUE;
593 break;
594 }
595 }
596 if(isOld == kFALSE){
597 GenTempMG1->Add(p->FindDaughter(MCParticle::kMu));
598 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu));
599 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu));
600 if (GetFillHist() && sumV[0] + 4*sumV[1] == 2)
601 hDVVMass[0]->Fill(TMath::Min(diBoson->Mass(),199.999));
602 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
603 hDVVMass[1]->Fill(TMath::Min(diBoson->Mass(),199.999));
604 }
605 }
606 if (p->HasDaughter(MCParticle::kEl) && p->HasDaughter(MCParticle::kElNu)) {
607 isOld = kFALSE;
608 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
609 if(p->FindDaughter(MCParticle::kEl) == GenTempMG1->At(nl)) {
610 isOld = kTRUE;
611 break;
612 }
613 }
614 if(isOld == kFALSE){
615 GenTempMG1->Add(p->FindDaughter(MCParticle::kEl));
616 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl));
617 diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu));
618 if (GetFillHist() && sumV[0] + 4*sumV[1] == 2)
619 hDVVMass[2]->Fill(TMath::Min(diBoson->Mass(),199.999));
620 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
621 hDVVMass[3]->Fill(TMath::Min(diBoson->Mass(),199.999));
622 }
623 }
624 if (p->HasDaughter(MCParticle::kTau) && p->HasDaughter(MCParticle::kTauNu)) {
625 isOld = kFALSE;
626 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
627 if(p->FindDaughter(MCParticle::kTau) == GenTempMG1->At(nl)) {
628 isOld = kTRUE;
629 break;
630 }
631 }
632 if(isOld == kFALSE){
633 GenTempMG1->Add(p->FindDaughter(MCParticle::kTau));
634 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau));
635 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu));
636 if (GetFillHist() && sumV[0] + 4*sumV[1] == 2)
637 hDVVMass[4]->Fill(TMath::Min(diBoson->Mass(),199.999));
638 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
639 hDVVMass[5]->Fill(TMath::Min(diBoson->Mass(),199.999));
640 }
641 }
642 if (p->HasDaughter(MCParticle::kMu,kTRUE) && p->HasDaughter(-1*MCParticle::kMu,kTRUE)) {
643 isOld = kFALSE;
644 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
645 if(p->FindDaughter(MCParticle::kMu,kTRUE) == GenTempMG1->At(nl)) {
646 isOld = kTRUE;
647 break;
648 }
649 }
650 if(isOld == kFALSE){
651 GenTempMG1->Add(p->FindDaughter(MCParticle::kMu,kTRUE));
652 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMu,kTRUE));
653 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kMu,kTRUE));
654 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
655 hDVVMass[6]->Fill(TMath::Min(diBoson->Mass(),199.999));
656 if (GetFillHist() && sumV[0] + 4*sumV[1] == 8)
657 hDVVMass[7]->Fill(TMath::Min(diBoson->Mass(),199.999));
658 }
659 }
660 if (p->HasDaughter(MCParticle::kEl,kTRUE) && p->HasDaughter(-1*MCParticle::kEl,kTRUE)) {
661 isOld = kFALSE;
662 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
663 if(p->FindDaughter(MCParticle::kEl,kTRUE) == GenTempMG1->At(nl)) {
664 isOld = kTRUE;
665 break;
666 }
667 }
668 if(isOld == kFALSE){
669 GenTempMG1->Add(p->FindDaughter(MCParticle::kEl,kTRUE));
670 diBoson->AddDaughter(p->FindDaughter(MCParticle::kEl,kTRUE));
671 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kEl,kTRUE));
672 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
673 hDVVMass[8]->Fill(TMath::Min(diBoson->Mass(),199.999));
674 if (GetFillHist() && sumV[0] + 4*sumV[1] == 8)
675 hDVVMass[9]->Fill(TMath::Min(diBoson->Mass(),199.999));
676 }
677 }
678 if (p->HasDaughter(MCParticle::kTau,kTRUE) && p->HasDaughter(-1*MCParticle::kTau,kTRUE)) {
679 isOld = kFALSE;
680 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
681 if(p->FindDaughter(MCParticle::kTau,kTRUE) == GenTempMG1->At(nl)) {
682 isOld = kTRUE;
683 break;
684 }
685 }
686 if(isOld == kFALSE){
687 GenTempMG1->Add(p->FindDaughter(MCParticle::kTau,kTRUE));
688 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTau,kTRUE));
689 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kTau,kTRUE));
690 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
691 hDVVMass[10]->Fill(TMath::Min(diBoson->Mass(),199.999));
692 if (GetFillHist() && sumV[0] + 4*sumV[1] == 8)
693 hDVVMass[11]->Fill(TMath::Min(diBoson->Mass(),199.999));
694 }
695 }
696 if (p->HasDaughter(MCParticle::kMuNu,kTRUE) && p->HasDaughter(-1*MCParticle::kMuNu,kTRUE)) {
697 isOld = kFALSE;
698 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
699 if(p->FindDaughter(MCParticle::kMuNu,kTRUE) == GenTempMG1->At(nl)) {
700 isOld = kTRUE;
701 break;
702 }
703 }
704 if(isOld == kFALSE){
705 GenTempMG1->Add(p->FindDaughter(MCParticle::kMuNu,kTRUE));
706 diBoson->AddDaughter(p->FindDaughter(MCParticle::kMuNu,kTRUE));
707 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kMuNu,kTRUE));
708 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
709 hDVVMass[12]->Fill(TMath::Min(diBoson->Mass(),199.999));
710 if (GetFillHist() && sumV[0] + 4*sumV[1] == 8)
711 hDVVMass[13]->Fill(TMath::Min(diBoson->Mass(),199.999));
712 }
713 }
714 if (p->HasDaughter(MCParticle::kElNu,kTRUE) && p->HasDaughter(-1*MCParticle::kElNu,kTRUE)) {
715 isOld = kFALSE;
716 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
717 if(p->FindDaughter(MCParticle::kElNu,kTRUE) == GenTempMG1->At(nl)) {
718 isOld = kTRUE;
719 break;
720 }
721 }
722 if(isOld == kFALSE){
723 GenTempMG1->Add(p->FindDaughter(MCParticle::kElNu,kTRUE));
724 diBoson->AddDaughter(p->FindDaughter(MCParticle::kElNu,kTRUE));
725 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kElNu,kTRUE));
726 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
727 hDVVMass[14]->Fill(TMath::Min(diBoson->Mass(),199.999));
728 if (GetFillHist() && sumV[0] + 4*sumV[1] == 8)
729 hDVVMass[15]->Fill(TMath::Min(diBoson->Mass(),199.999));
730 }
731 }
732 if (p->HasDaughter(MCParticle::kTauNu,kTRUE) &&
733 p->HasDaughter(-1*MCParticle::kTauNu,kTRUE)) {
734 isOld = kFALSE;
735 for(UInt_t nl = 0; nl < GenTempMG1->GetEntries(); nl++){
736 if(p->FindDaughter(MCParticle::kTauNu,kTRUE) == GenTempMG1->At(nl)) {
737 isOld = kTRUE;
738 break;
739 }
740 }
741 if(isOld == kFALSE){
742 GenTempMG1->Add(p->FindDaughter(MCParticle::kTauNu,kTRUE));
743 diBoson->AddDaughter(p->FindDaughter(MCParticle::kTauNu,kTRUE));
744 diBoson->AddDaughter(p->FindDaughter(-1*MCParticle::kTauNu,kTRUE));
745 if (GetFillHist() && sumV[0] + 4*sumV[1] == 5)
746 hDVVMass[16]->Fill(TMath::Min(diBoson->Mass(),199.999));
747 if (GetFillHist() && sumV[0] + 4*sumV[1] == 8)
748 hDVVMass[17]->Fill(TMath::Min(diBoson->Mass(),199.999));
749 }
750 }
751 if (diBoson && diBosonMass[0] <= 0) diBosonMass[0] = diBoson->Mass();
752 else if(diBoson && diBosonMass[1] <= 0) diBosonMass[1] = diBoson->Mass();
753 delete diBoson;
754 }
755 else if (p->Status() == 3 && (p->Is(MCParticle::kZ) || p->Is(MCParticle::kW))) {
756 if (diBosonMass[0] <= 0) diBosonMass[0] = p->Mass();
757 else if(diBosonMass[1] <= 0) diBosonMass[1] = p->Mass();
758 if (GetFillHist()) {
759 if (sumV[0] + 4*sumV[1] == 2 && p->Is(MCParticle::kW) &&
760 p->HasDaughter(MCParticle::kMu) &&
761 p->HasDaughter(MCParticle::kMuNu))
762 hDVVMass[0]->Fill(TMath::Min(p->Mass(),199.999));
763 else if(sumV[0] + 4*sumV[1] == 2 && p->Is(MCParticle::kW) &&
764 p->HasDaughter(MCParticle::kEl) &&
765 p->HasDaughter(MCParticle::kElNu))
766 hDVVMass[2]->Fill(TMath::Min(p->Mass(),199.999));
767 else if(sumV[0] + 4*sumV[1] == 2 && p->Is(MCParticle::kW) &&
768 p->HasDaughter(MCParticle::kTau) &&
769 p->HasDaughter(MCParticle::kTauNu))
770 hDVVMass[4]->Fill(TMath::Min(p->Mass(),199.999));
771 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kW) &&
772 p->HasDaughter(MCParticle::kMu) &&
773 p->HasDaughter(MCParticle::kMuNu))
774 hDVVMass[1]->Fill(TMath::Min(p->Mass(),199.999));
775 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kW) &&
776 p->HasDaughter(MCParticle::kEl) &&
777 p->HasDaughter(MCParticle::kElNu))
778 hDVVMass[3]->Fill(TMath::Min(p->Mass(),199.999));
779 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kW) &&
780 p->HasDaughter(MCParticle::kTau) &&
781 p->HasDaughter(MCParticle::kTauNu))
782 hDVVMass[5]->Fill(TMath::Min(p->Mass(),199.999));
783 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kZ) &&
784 p->HasDaughter(MCParticle::kMu,kTRUE) &&
785 p->HasDaughter(-1*MCParticle::kMu,kTRUE))
786 hDVVMass[6]->Fill(TMath::Min(p->Mass(),199.999));
787 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kZ) &&
788 p->HasDaughter(MCParticle::kEl,kTRUE) &&
789 p->HasDaughter(-1*MCParticle::kEl,kTRUE))
790 hDVVMass[8]->Fill(TMath::Min(p->Mass(),199.999));
791 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kZ) &&
792 p->HasDaughter(MCParticle::kTau,kTRUE) &&
793 p->HasDaughter(-1*MCParticle::kTau,kTRUE))
794 hDVVMass[10]->Fill(TMath::Min(p->Mass(),199.999));
795 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kZ) &&
796 p->HasDaughter(MCParticle::kMuNu,kTRUE) &&
797 p->HasDaughter(-1*MCParticle::kMuNu,kTRUE))
798 hDVVMass[12]->Fill(TMath::Min(p->Mass(),199.999));
799 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kZ) &&
800 p->HasDaughter(MCParticle::kElNu,kTRUE) &&
801 p->HasDaughter(-1*MCParticle::kElNu,kTRUE))
802 hDVVMass[14]->Fill(TMath::Min(p->Mass(),199.999));
803 else if(sumV[0] + 4*sumV[1] == 5 && p->Is(MCParticle::kZ) &&
804 p->HasDaughter(MCParticle::kTauNu,kTRUE) &&
805 p->HasDaughter(-1*MCParticle::kTauNu,kTRUE))
806 hDVVMass[16]->Fill(TMath::Min(p->Mass(),199.999));
807 else if(sumV[0] + 4*sumV[1] == 8 && p->Is(MCParticle::kZ) &&
808 p->HasDaughter(MCParticle::kMu,kTRUE) &&
809 p->HasDaughter(-1*MCParticle::kMu,kTRUE))
810 hDVVMass[7]->Fill(TMath::Min(p->Mass(),199.999));
811 else if(sumV[0] + 4*sumV[1] == 8 && p->Is(MCParticle::kZ) &&
812 p->HasDaughter(MCParticle::kEl,kTRUE) &&
813 p->HasDaughter(-1*MCParticle::kEl,kTRUE))
814 hDVVMass[9]->Fill(TMath::Min(p->Mass(),199.999));
815 else if(sumV[0] + 4*sumV[1] == 8 && p->Is(MCParticle::kZ) &&
816 p->HasDaughter(MCParticle::kTau,kTRUE) &&
817 p->HasDaughter(-1*MCParticle::kTau,kTRUE))
818 hDVVMass[11]->Fill(TMath::Min(p->Mass(),199.999));
819 else if(sumV[0] + 4*sumV[1] == 8 && p->Is(MCParticle::kZ) &&
820 p->HasDaughter(MCParticle::kMuNu,kTRUE) &&
821 p->HasDaughter(-1*MCParticle::kMuNu,kTRUE))
822 hDVVMass[13]->Fill(TMath::Min(p->Mass(),199.999));
823 else if(sumV[0] + 4*sumV[1] == 8 && p->Is(MCParticle::kZ) &&
824 p->HasDaughter(MCParticle::kElNu,kTRUE) &&
825 p->HasDaughter(-1*MCParticle::kElNu,kTRUE))
826 hDVVMass[15]->Fill(TMath::Min(p->Mass(),199.999));
827 else if(sumV[0] + 4*sumV[1] == 8 && p->Is(MCParticle::kZ) &&
828 p->HasDaughter(MCParticle::kTauNu,kTRUE) &&
829 p->HasDaughter(-1*MCParticle::kTauNu,kTRUE))
830 hDVVMass[17]->Fill(TMath::Min(p->Mass(),199.999));
831 }
832 }
833 } // end loop of particles
834 if (GetFillHist()) {
835 if(diBosonMass[0] > 70 && diBosonMass[0] < 110 &&
836 diBosonMass[1] > 70 && diBosonMass[1] < 110){
837 if(sumV[0] + 4*sumV[1] == 2){
838 if (sumVVFlavor[0] == 2) hDVVMass[18]->Fill(0.);
839 else if(sumVVFlavor[1] == 2) hDVVMass[18]->Fill(1.);
840 else if(sumVVFlavor[2] == 2) hDVVMass[18]->Fill(2.);
841 else if(sumVVFlavor[0] == 1 && sumVVFlavor[1] == 1) hDVVMass[18]->Fill(3.);
842 else if(sumVVFlavor[0] == 1 && sumVVFlavor[2] == 1) hDVVMass[18]->Fill(4.);
843 else if(sumVVFlavor[1] == 1 && sumVVFlavor[2] == 1) hDVVMass[18]->Fill(5.);
844 else hDVVMass[18]->Fill(6.);
845 }
846 if(sumV[0] + 4*sumV[1] == 5){
847 if (sumVVFlavor[3] == 1 && sumVVFlavor[0] == 1) hDVVMass[19]->Fill(0.);
848 else if(sumVVFlavor[3] == 1 && sumVVFlavor[1] == 1) hDVVMass[19]->Fill(1.);
849 else if(sumVVFlavor[3] == 1 && sumVVFlavor[2] == 1) hDVVMass[19]->Fill(2.);
850 else if(sumVVFlavor[4] == 1 && sumVVFlavor[0] == 1) hDVVMass[19]->Fill(3.);
851 else if(sumVVFlavor[4] == 1 && sumVVFlavor[1] == 1) hDVVMass[19]->Fill(4.);
852 else if(sumVVFlavor[4] == 1 && sumVVFlavor[2] == 1) hDVVMass[19]->Fill(5.);
853 else if(sumVVFlavor[5] == 1 && sumVVFlavor[0] == 1) hDVVMass[19]->Fill(6.);
854 else if(sumVVFlavor[5] == 1 && sumVVFlavor[1] == 1) hDVVMass[19]->Fill(7.);
855 else if(sumVVFlavor[5] == 1 && sumVVFlavor[2] == 1) hDVVMass[19]->Fill(8.);
856 else { hDVVMass[19]->Fill(9.);
857 /*for(int i=0; i<9; i++) cout << sumVVFlavor[i] << " " ;cout << endl;*/}
858 }
859 if(sumV[0] + 4*sumV[1] == 8 &&
860 sumVVFlavor[3] + sumVVFlavor[4] + sumVVFlavor[5] == 2){
861 if (sumVVFlavor[3] == 2) hDVVMass[20]->Fill(0.);
862 else if(sumVVFlavor[4] == 2) hDVVMass[20]->Fill(1.);
863 else if(sumVVFlavor[5] == 2) hDVVMass[20]->Fill(2.);
864 else if(sumVVFlavor[3] == 1 && sumVVFlavor[4] == 1) hDVVMass[20]->Fill(3.);
865 else if(sumVVFlavor[3] == 1 && sumVVFlavor[5] == 1) hDVVMass[20]->Fill(4.);
866 else if(sumVVFlavor[4] == 1 && sumVVFlavor[5] == 1) hDVVMass[20]->Fill(5.);
867 else hDVVMass[20]->Fill(6.);
868 }
869 else if(sumV[0] + 4*sumV[1] == 8){
870 if (sumVVFlavor[6] == 2) hDVVMass[21]->Fill(0.);
871 else if(sumVVFlavor[7] == 2) hDVVMass[21]->Fill(1.);
872 else if(sumVVFlavor[8] == 2) hDVVMass[21]->Fill(2.);
873 else if(sumVVFlavor[3] == 1 && sumVVFlavor[6] == 1) hDVVMass[21]->Fill(3.);
874 else if(sumVVFlavor[3] == 1 && sumVVFlavor[7] == 1) hDVVMass[21]->Fill(4.);
875 else if(sumVVFlavor[3] == 1 && sumVVFlavor[8] == 1) hDVVMass[21]->Fill(5.);
876 else if(sumVVFlavor[4] == 1 && sumVVFlavor[6] == 1) hDVVMass[21]->Fill(6.);
877 else if(sumVVFlavor[4] == 1 && sumVVFlavor[7] == 1) hDVVMass[21]->Fill(7.);
878 else if(sumVVFlavor[4] == 1 && sumVVFlavor[8] == 1) hDVVMass[21]->Fill(8.);
879 else if(sumVVFlavor[5] == 1 && sumVVFlavor[6] == 1) hDVVMass[21]->Fill(9.);
880 else if(sumVVFlavor[5] == 1 && sumVVFlavor[7] == 1) hDVVMass[21]->Fill(10.);
881 else if(sumVVFlavor[5] == 1 && sumVVFlavor[8] == 1) hDVVMass[21]->Fill(11.);
882 else if(sumVVFlavor[6] == 1 && sumVVFlavor[7] == 1) hDVVMass[21]->Fill(12.);
883 else if(sumVVFlavor[6] == 1 && sumVVFlavor[8] == 1) hDVVMass[21]->Fill(13.);
884 else if(sumVVFlavor[7] == 1 && sumVVFlavor[8] == 1) hDVVMass[21]->Fill(14.);
885 else hDVVMass[21]->Fill(15.);
886 }
887 } // 60<mV1/2<120
888 if(sumV[0] + 4*sumV[1] == 2)
889 hDVVMass[22]->Fill(TMath::Min(TMath::Min(diBosonMass[0],diBosonMass[1]),199.999));
890 if(sumV[0] + 4*sumV[1] == 2)
891 hDVVMass[23]->Fill(TMath::Min(TMath::Max(diBosonMass[0],diBosonMass[1]),199.999));
892 if(sumV[0] + 4*sumV[1] == 5)
893 hDVVMass[24]->Fill(TMath::Min(TMath::Min(diBosonMass[0],diBosonMass[1]),199.999));
894 if(sumV[0] + 4*sumV[1] == 5)
895 hDVVMass[25]->Fill(TMath::Min(TMath::Max(diBosonMass[0],diBosonMass[1]),199.999));
896 if(sumV[0] + 4*sumV[1] == 8)
897 hDVVMass[26]->Fill(TMath::Min(TMath::Min(diBosonMass[0],diBosonMass[1]),199.999));
898 if(sumV[0] + 4*sumV[1] == 8)
899 hDVVMass[27]->Fill(TMath::Min(TMath::Max(diBosonMass[0],diBosonMass[1]),199.999));
900 }
901 delete GenTempMG1;
902 } // WW, WZ or ZZ
903 // --------------------------------
904 // End special study about VVjets
905 // --------------------------------
906
907 Met *theMET = new Met(-totalMET[0], -totalMET[1]);
908 theMET->SetElongitudinal(-totalMET[2]);
909 GenMet->AddOwned(theMET);
910
911 // sort according to pt
912 GenLeptons->Sort();
913 GenAllLeptons->Sort();
914 GenTaus->Sort();
915 GenNeutrinos->Sort();
916 GenQuarks->Sort();
917 GenqqHs->Sort();
918 GenBosons->Sort();
919 GenPhotons->Sort();
920 GenRadPhotons->Sort();
921 GenISRPhotons->Sort();
922 } // Only for Monte Carlo
923
924 // add objects to this event for other modules to use
925 AddObjThisEvt(GenMet);
926 AddObjThisEvt(GenLeptons);
927 AddObjThisEvt(GenAllLeptons);
928 AddObjThisEvt(GenTaus);
929 AddObjThisEvt(GenNeutrinos);
930 AddObjThisEvt(GenQuarks);
931 AddObjThisEvt(GenqqHs);
932 AddObjThisEvt(GenBosons);
933 AddObjThisEvt(GenPhotons);
934 AddObjThisEvt(GenRadPhotons);
935 AddObjThisEvt(GenISRPhotons);
936
937 // --------------------------------
938 // Copy these Collections into the Arrays for Publication for Output Module
939 // --------------------------------
940 fGenLeptons->Delete();
941 fGenAllLeptons->Delete();
942 fGenTaus->Delete();
943 fGenNeutrinos->Delete();
944 fGenQuarks->Delete();
945 fGenqqHs->Delete();
946 fGenBosons->Delete();
947 fGenPhotons->Delete();
948 fGenRadPhotons->Delete();
949 fGenISRPhotons->Delete();
950
951 for (UInt_t i=0; i < GenLeptons->GetEntries(); ++i) {
952 mithep::MCParticle *genParticle = fGenLeptons->Allocate();
953 new (genParticle) mithep::MCParticle(GenLeptons->At(i)->Px(),GenLeptons->At(i)->Py(),
954 GenLeptons->At(i)->Pz(),GenLeptons->At(i)->E(),
955 GenLeptons->At(i)->PdgId(),GenLeptons->At(i)->Status());
956 }
957 for (UInt_t i=0; i < GenAllLeptons->GetEntries(); ++i) {
958 mithep::MCParticle *genParticle = fGenAllLeptons->Allocate();
959 new (genParticle) mithep::MCParticle(GenAllLeptons->At(i)->Px(),GenAllLeptons->At(i)->Py(),
960 GenAllLeptons->At(i)->Pz(),GenAllLeptons->At(i)->E(),
961 GenAllLeptons->At(i)->PdgId(),GenAllLeptons->At(i)->Status());
962 }
963 for (UInt_t i=0; i < GenTaus->GetEntries(); ++i) {
964 mithep::MCParticle *genParticle = fGenTaus->Allocate();
965 new (genParticle) mithep::MCParticle(GenTaus->At(i)->Px(),GenTaus->At(i)->Py(),
966 GenTaus->At(i)->Pz(),GenTaus->At(i)->E(),
967 GenTaus->At(i)->PdgId(),GenTaus->At(i)->Status());
968 }
969 for (UInt_t i=0; i < GenNeutrinos->GetEntries(); ++i) {
970 mithep::MCParticle *genParticle = fGenNeutrinos->Allocate();
971 new (genParticle) mithep::MCParticle(GenNeutrinos->At(i)->Px(),GenNeutrinos->At(i)->Py(),
972 GenNeutrinos->At(i)->Pz(),GenNeutrinos->At(i)->E(),
973 GenNeutrinos->At(i)->PdgId(),GenNeutrinos->At(i)->Status());
974 }
975 for (UInt_t i=0; i < GenQuarks->GetEntries(); ++i) {
976 mithep::MCParticle *genParticle = fGenQuarks->Allocate();
977 new (genParticle) mithep::MCParticle(GenQuarks->At(i)->Px(),GenQuarks->At(i)->Py(),
978 GenQuarks->At(i)->Pz(),GenQuarks->At(i)->E(),
979 GenQuarks->At(i)->PdgId(),GenQuarks->At(i)->Status());
980 }
981 for (UInt_t i=0; i < GenqqHs->GetEntries(); ++i) {
982 mithep::MCParticle *genParticle = fGenqqHs->Allocate();
983 new (genParticle) mithep::MCParticle(GenqqHs->At(i)->Px(),GenqqHs->At(i)->Py(),
984 GenqqHs->At(i)->Pz(),GenqqHs->At(i)->E(),
985 GenqqHs->At(i)->PdgId(),GenqqHs->At(i)->Status());
986 }
987
988 if (fCopyArrays) {
989 // --------------------------------
990 // Copy these Collections into the Arrays for Publication for Output Module
991 // --------------------------------
992 fGenLeptons->Delete();
993 fGenAllLeptons->Delete();
994 fGenTaus->Delete();
995 fGenNeutrinos->Delete();
996 fGenQuarks->Delete();
997 fGenqqHs->Delete();
998 fGenBosons->Delete();
999 fGenPhotons->Delete();
1000 fGenRadPhotons->Delete();
1001 fGenISRPhotons->Delete();
1002
1003 for (UInt_t i=0; i < GenLeptons->GetEntries(); ++i) {
1004 mithep::MCParticle *genParticle = fGenLeptons->Allocate();
1005 new (genParticle) mithep::MCParticle(GenLeptons->At(i)->Px(),GenLeptons->At(i)->Py(),
1006 GenLeptons->At(i)->Pz(),GenLeptons->At(i)->E(),
1007 GenLeptons->At(i)->PdgId(),GenLeptons->At(i)->Status());
1008 }
1009 for (UInt_t i=0; i < GenAllLeptons->GetEntries(); ++i) {
1010 mithep::MCParticle *genParticle = fGenAllLeptons->Allocate();
1011 new (genParticle) mithep::MCParticle(GenAllLeptons->At(i)->Px(),GenAllLeptons->At(i)->Py(),
1012 GenAllLeptons->At(i)->Pz(),GenAllLeptons->At(i)->E(),
1013 GenAllLeptons->At(i)->PdgId(),
1014 GenAllLeptons->At(i)->Status());
1015 }
1016 for (UInt_t i=0; i < GenTaus->GetEntries(); ++i) {
1017 mithep::MCParticle *genParticle = fGenTaus->Allocate();
1018 new (genParticle) mithep::MCParticle(GenTaus->At(i)->Px(),GenTaus->At(i)->Py(),
1019 GenTaus->At(i)->Pz(),GenTaus->At(i)->E(),
1020 GenTaus->At(i)->PdgId(),
1021 GenTaus->At(i)->Status());
1022 }
1023 for (UInt_t i=0; i < GenNeutrinos->GetEntries(); ++i) {
1024 mithep::MCParticle *genParticle = fGenNeutrinos->Allocate();
1025 new (genParticle) mithep::MCParticle(GenNeutrinos->At(i)->Px(),GenNeutrinos->At(i)->Py(),
1026 GenNeutrinos->At(i)->Pz(),GenNeutrinos->At(i)->E(),
1027 GenNeutrinos->At(i)->PdgId(),
1028 GenNeutrinos->At(i)->Status());
1029 }
1030 for (UInt_t i=0; i < GenQuarks->GetEntries(); ++i) {
1031 mithep::MCParticle *genParticle = fGenQuarks->Allocate();
1032 new (genParticle) mithep::MCParticle(GenQuarks->At(i)->Px(),GenQuarks->At(i)->Py(),
1033 GenQuarks->At(i)->Pz(),GenQuarks->At(i)->E(),
1034 GenQuarks->At(i)->PdgId(),
1035 GenQuarks->At(i)->Status());
1036 }
1037 for (UInt_t i=0; i < GenqqHs->GetEntries(); ++i) {
1038 mithep::MCParticle *genParticle = fGenqqHs->Allocate();
1039 new (genParticle) mithep::MCParticle(GenqqHs->At(i)->Px(),GenqqHs->At(i)->Py(),
1040 GenqqHs->At(i)->Pz(),GenqqHs->At(i)->E(),
1041 GenqqHs->At(i)->PdgId(),
1042 GenqqHs->At(i)->Status());
1043 }
1044 for (UInt_t i=0; i < GenBosons->GetEntries(); ++i) {
1045 mithep::MCParticle *genParticle = fGenBosons->Allocate();
1046 new (genParticle) mithep::MCParticle(GenBosons->At(i)->Px(),GenBosons->At(i)->Py(),
1047 GenBosons->At(i)->Pz(),GenBosons->At(i)->E(),
1048 GenBosons->At(i)->PdgId(),
1049 GenBosons->At(i)->Status());
1050 }
1051 for (UInt_t i=0; i < GenPhotons->GetEntries(); ++i) {
1052 mithep::MCParticle *genParticle = fGenPhotons->Allocate();
1053 new (genParticle) mithep::MCParticle(GenPhotons->At(i)->Px(),GenPhotons->At(i)->Py(),
1054 GenPhotons->At(i)->Pz(),GenPhotons->At(i)->E(),
1055 GenPhotons->At(i)->PdgId(),
1056 GenPhotons->At(i)->Status());
1057 }
1058 for (UInt_t i=0; i < GenRadPhotons->GetEntries(); ++i) {
1059 mithep::MCParticle *genParticle = fGenRadPhotons->Allocate();
1060 new (genParticle) mithep::MCParticle(GenRadPhotons->At(i)->Px(),GenRadPhotons->At(i)->Py(),
1061 GenRadPhotons->At(i)->Pz(),GenRadPhotons->At(i)->E(),
1062 GenRadPhotons->At(i)->PdgId(),
1063 GenRadPhotons->At(i)->Status());
1064 }
1065 for (UInt_t i=0; i < GenISRPhotons->GetEntries(); ++i) {
1066 mithep::MCParticle *genParticle = fGenISRPhotons->Allocate();
1067 new (genParticle) mithep::MCParticle(GenISRPhotons->At(i)->Px(),GenISRPhotons->At(i)->Py(),
1068 GenISRPhotons->At(i)->Pz(),GenISRPhotons->At(i)->E(),
1069 GenISRPhotons->At(i)->PdgId(),
1070 GenISRPhotons->At(i)->Status());
1071 }
1072 }
1073
1074 for (UInt_t i=0; i < GenBosons->GetEntries(); ++i) {
1075 mithep::MCParticle *genParticle = fGenBosons->Allocate();
1076 new (genParticle) mithep::MCParticle(GenBosons->At(i)->Px(),GenBosons->At(i)->Py(),
1077 GenBosons->At(i)->Pz(),GenBosons->At(i)->E(),
1078 GenBosons->At(i)->PdgId(),GenBosons->At(i)->Status());
1079 }
1080 for (UInt_t i=0; i < GenPhotons->GetEntries(); ++i) {
1081 mithep::MCParticle *genParticle = fGenPhotons->Allocate();
1082 new (genParticle) mithep::MCParticle(GenPhotons->At(i)->Px(),GenPhotons->At(i)->Py(),
1083 GenPhotons->At(i)->Pz(),GenPhotons->At(i)->E(),
1084 GenPhotons->At(i)->PdgId(),GenPhotons->At(i)->Status());
1085 }
1086 for (UInt_t i=0; i < GenRadPhotons->GetEntries(); ++i) {
1087 mithep::MCParticle *genParticle = fGenRadPhotons->Allocate();
1088 new (genParticle) mithep::MCParticle(GenRadPhotons->At(i)->Px(),GenRadPhotons->At(i)->Py(),
1089 GenRadPhotons->At(i)->Pz(),GenRadPhotons->At(i)->E(),
1090 GenRadPhotons->At(i)->PdgId(),GenRadPhotons->At(i)->Status());
1091 }
1092 for (UInt_t i=0; i < GenISRPhotons->GetEntries(); ++i) {
1093 mithep::MCParticle *genParticle = fGenISRPhotons->Allocate();
1094 new (genParticle) mithep::MCParticle(GenISRPhotons->At(i)->Px(),GenISRPhotons->At(i)->Py(),
1095 GenISRPhotons->At(i)->Pz(),GenISRPhotons->At(i)->E(),
1096 GenISRPhotons->At(i)->PdgId(),GenISRPhotons->At(i)->Status());
1097 }
1098
1099 //bool theZCut = (GenLeptons->GetEntries() > 0 || GenTaus->GetEntries() > 0) && GenQuarks->GetEntries() == 0;
1100 //if(!theZCut){
1101 // SkipEvent();
1102 // return;
1103 //}
1104
1105 // fill histograms if requested
1106 if (GetFillHist()) {
1107 // MET
1108 hDGenMet[0]->Fill(GenMet->At(0)->Pt());
1109 hDGenMet[1]->Fill(GenMet->At(0)->Px());
1110 hDGenMet[2]->Fill(GenMet->At(0)->Py());
1111 hDGenMet[3]->Fill(GenMet->At(0)->Elongitudinal());
1112
1113 // pt min for leptons from W/Z
1114 hDGenPtMin->Fill(TMath::Min(ptMin, 199.999));
1115 //if(ptMin < 10){
1116 // SkipEvent();
1117 // return;
1118 //}
1119
1120 // leptons
1121 hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
1122 for(UInt_t i=0; i<GenLeptons->GetEntries(); i++) {
1123 hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
1124 hDGenLeptons[2]->Fill(TMath::Abs(GenLeptons->At(i)->Eta()));
1125 hDGenLeptons[3]->Fill(GenLeptons->At(i)->PhiDeg());
1126 for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++) {
1127 CompositeParticle *dilepton = new CompositeParticle();
1128 dilepton->AddDaughter(GenLeptons->At(i));
1129 dilepton->AddDaughter(GenLeptons->At(j));
1130 hDGenLeptons[4]->Fill(dilepton->Mass());
1131 delete dilepton;
1132 }
1133 }
1134 // looking at events with two leptons
1135 if (GenLeptons->GetEntries() == 2) {
1136 hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(TMath::Abs(GenLeptons->At(0)->Eta()),
1137 TMath::Abs(GenLeptons->At(1)->Eta())),
1138 4.999));
1139 hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(TMath::Abs(GenLeptons->At(0)->Eta()),
1140 TMath::Abs(GenLeptons->At(1)->Eta())),
1141 4.999));
1142 if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
1143 TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5) {
1144 hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
1145 if (GenLeptons->At(0)->Pt() > 20.0) {
1146 hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
1147 if (GenLeptons->At(1)->Pt() > 10.0) {
1148 CompositeParticle *dilepton = new CompositeParticle();
1149 dilepton->AddDaughter(GenLeptons->At(0));
1150 dilepton->AddDaughter(GenLeptons->At(1));
1151 hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
1152 if(dilepton->Mass() > 12.0){
1153 hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(0)->Phi(),
1154 GenLeptons->At(1)->Phi())
1155 * 180./ TMath::Pi());
1156 hDGenLeptons[11]->Fill(MathUtils::DeltaR(*GenLeptons->At(0),
1157 *GenLeptons->At(1)));
1158 }
1159 delete dilepton;
1160 }
1161 }
1162 }
1163 }
1164 // looking at events with three leptons
1165 if (GenLeptons->GetEntries() == 3) {
1166 if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
1167 TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5 &&
1168 TMath::Abs(GenLeptons->At(2)->Eta()) < 2.5) {
1169 hDGenLeptons[12]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
1170 if (GenLeptons->At(0)->Pt() > 20.0) {
1171 hDGenLeptons[13]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
1172 hDGenLeptons[14]->Fill(TMath::Min(GenLeptons->At(2)->Pt(),199.999));
1173 if (GenLeptons->At(1)->Pt() > 10.0 && GenLeptons->At(2)->Pt() > 10.0) {
1174 CompositeParticle *dilepton01 = new CompositeParticle();
1175 dilepton01->AddDaughter(GenLeptons->At(0));
1176 dilepton01->AddDaughter(GenLeptons->At(1));
1177 CompositeParticle *dilepton02 = new CompositeParticle();
1178 dilepton02->AddDaughter(GenLeptons->At(0));
1179 dilepton02->AddDaughter(GenLeptons->At(2));
1180 CompositeParticle *dilepton12 = new CompositeParticle();
1181 dilepton12->AddDaughter(GenLeptons->At(1));
1182 dilepton12->AddDaughter(GenLeptons->At(2));
1183 hDGenLeptons[15]->Fill(TMath::Min(dilepton01->Mass(),999.999));
1184 hDGenLeptons[15]->Fill(TMath::Min(dilepton02->Mass(),999.999));
1185 hDGenLeptons[15]->Fill(TMath::Min(dilepton12->Mass(),999.999));
1186 CompositeParticle *trilepton = new CompositeParticle();
1187 trilepton->AddDaughter(GenLeptons->At(0));
1188 trilepton->AddDaughter(GenLeptons->At(1));
1189 trilepton->AddDaughter(GenLeptons->At(2));
1190 hDGenLeptons[16]->Fill(TMath::Min(trilepton->Mass(),999.999));
1191 Double_t deltaR[3] = {MathUtils::DeltaR(*GenLeptons->At(0),
1192 *GenLeptons->At(1)),
1193 MathUtils::DeltaR(*GenLeptons->At(0),
1194 *GenLeptons->At(2)),
1195 MathUtils::DeltaR(*GenLeptons->At(1),
1196 *GenLeptons->At(2))};
1197 Double_t deltaRMin = deltaR[0];
1198 for(Int_t i=1; i<3; i++)
1199 if(deltaRMin > deltaR[i])
1200 deltaRMin = deltaR[i];
1201 hDGenLeptons[17]->Fill(deltaRMin);
1202
1203 delete dilepton01;
1204 delete dilepton02;
1205 delete dilepton12;
1206 delete trilepton;
1207 }
1208 }
1209 }
1210 }
1211 // looking at events with four leptons
1212 if (GenLeptons->GetEntries() == 4) {
1213 if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
1214 TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5 &&
1215 TMath::Abs(GenLeptons->At(2)->Eta()) < 2.5 &&
1216 TMath::Abs(GenLeptons->At(3)->Eta()) < 2.5) {
1217 hDGenLeptons[18]->Fill(TMath::Min(GenLeptons->At(0)->Pt(),199.999));
1218 if (GenLeptons->At(0)->Pt() > 20.0) {
1219 hDGenLeptons[19]->Fill(TMath::Min(GenLeptons->At(1)->Pt(),199.999));
1220 hDGenLeptons[20]->Fill(TMath::Min(GenLeptons->At(2)->Pt(),199.999));
1221 hDGenLeptons[21]->Fill(TMath::Min(GenLeptons->At(3)->Pt(),199.999));
1222 if (GenLeptons->At(1)->Pt() > 10.0 && GenLeptons->At(2)->Pt() > 10.0 &&
1223 GenLeptons->At(3)->Pt() > 10.0) {
1224 CompositeParticle *dilepton01 = new CompositeParticle();
1225 dilepton01->AddDaughter(GenLeptons->At(0));
1226 dilepton01->AddDaughter(GenLeptons->At(1));
1227 CompositeParticle *dilepton02 = new CompositeParticle();
1228 dilepton02->AddDaughter(GenLeptons->At(0));
1229 dilepton02->AddDaughter(GenLeptons->At(2));
1230 CompositeParticle *dilepton03 = new CompositeParticle();
1231 dilepton03->AddDaughter(GenLeptons->At(0));
1232 dilepton03->AddDaughter(GenLeptons->At(3));
1233 CompositeParticle *dilepton12 = new CompositeParticle();
1234 dilepton12->AddDaughter(GenLeptons->At(1));
1235 dilepton12->AddDaughter(GenLeptons->At(2));
1236 CompositeParticle *dilepton13 = new CompositeParticle();
1237 dilepton13->AddDaughter(GenLeptons->At(1));
1238 dilepton13->AddDaughter(GenLeptons->At(3));
1239 CompositeParticle *dilepton23 = new CompositeParticle();
1240 dilepton23->AddDaughter(GenLeptons->At(2));
1241 dilepton23->AddDaughter(GenLeptons->At(3));
1242 hDGenLeptons[22]->Fill(TMath::Min(dilepton01->Mass(),999.999));
1243 hDGenLeptons[22]->Fill(TMath::Min(dilepton02->Mass(),999.999));
1244 hDGenLeptons[22]->Fill(TMath::Min(dilepton03->Mass(),999.999));
1245 hDGenLeptons[22]->Fill(TMath::Min(dilepton12->Mass(),999.999));
1246 hDGenLeptons[22]->Fill(TMath::Min(dilepton13->Mass(),999.999));
1247 hDGenLeptons[22]->Fill(TMath::Min(dilepton23->Mass(),999.999));
1248 CompositeParticle *fourlepton = new CompositeParticle();
1249 fourlepton->AddDaughter(GenLeptons->At(0));
1250 fourlepton->AddDaughter(GenLeptons->At(1));
1251 fourlepton->AddDaughter(GenLeptons->At(2));
1252 fourlepton->AddDaughter(GenLeptons->At(3));
1253 hDGenLeptons[23]->Fill(TMath::Min(fourlepton->Mass(),999.999));
1254 Double_t deltaR[6] = {MathUtils::DeltaR(*GenLeptons->At(0),
1255 *GenLeptons->At(1)),
1256 MathUtils::DeltaR(*GenLeptons->At(0),
1257 *GenLeptons->At(2)),
1258 MathUtils::DeltaR(*GenLeptons->At(0),
1259 *GenLeptons->At(3)),
1260 MathUtils::DeltaR(*GenLeptons->At(1),
1261 *GenLeptons->At(2)),
1262 MathUtils::DeltaR(*GenLeptons->At(1),
1263 *GenLeptons->At(3)),
1264 MathUtils::DeltaR(*GenLeptons->At(2),
1265 *GenLeptons->At(3))};
1266 Double_t deltaRMin = deltaR[0];
1267 for(Int_t i=1; i<6; i++)
1268 if(deltaRMin > deltaR[i])
1269 deltaRMin = deltaR[i];
1270 hDGenLeptons[24]->Fill(deltaRMin);
1271
1272 delete dilepton01;
1273 delete dilepton02;
1274 delete dilepton03;
1275 delete dilepton12;
1276 delete dilepton13;
1277 delete dilepton23;
1278 delete fourlepton;
1279 }
1280 }
1281 }
1282 }
1283
1284 // all leptons
1285 hDGenAllLeptons[0]->Fill(GenAllLeptons->GetEntries());
1286 for(UInt_t i=0; i<GenAllLeptons->GetEntries(); i++) {
1287 hDGenAllLeptons[1]->Fill(GenAllLeptons->At(i)->Pt());
1288 hDGenAllLeptons[2]->Fill(GenAllLeptons->At(i)->AbsEta());
1289 hDGenAllLeptons[3]->Fill(GenAllLeptons->At(i)->PhiDeg());
1290 }
1291 if(GenAllLeptons->GetEntries() >= 2) hDGenAllLeptons[4]->Fill(GenAllLeptons->At(1)->Pt());
1292
1293 // taus
1294 hDGenTaus[0]->Fill(GenTaus->GetEntries());
1295 for(UInt_t i=0; i<GenTaus->GetEntries(); i++) {
1296 hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
1297 hDGenTaus[2]->Fill(GenTaus->At(i)->AbsEta());
1298 hDGenTaus[3]->Fill(GenTaus->At(i)->PhiDeg());
1299 }
1300
1301 // neutrinos
1302 hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
1303 CompositeParticle *neutrinoTotal = new CompositeParticle();
1304 for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++) {
1305 if (GenNeutrinos->At(i)->HasMother())
1306 neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
1307 }
1308 if (GenNeutrinos->GetEntries() > 0) {
1309 hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
1310 hDGenNeutrinos[2]->Fill(neutrinoTotal->AbsEta());
1311 hDGenNeutrinos[3]->Fill(neutrinoTotal->PhiDeg());
1312 }
1313 delete neutrinoTotal;
1314
1315 // quarks
1316 hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
1317 for(UInt_t i=0; i<GenQuarks->GetEntries(); i++) {
1318 for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++) {
1319 CompositeParticle *dijet = new CompositeParticle();
1320 dijet->AddDaughter(GenQuarks->At(i));
1321 dijet->AddDaughter(GenQuarks->At(j));
1322 hDGenQuarks[1]->Fill(dijet->Pt());
1323 hDGenQuarks[2]->Fill(dijet->Mass());
1324 if (TMath::Abs(GenQuarks->At(i)->Eta()) < 2.5 &&
1325 TMath::Abs(GenQuarks->At(j)->Eta()) < 2.5) {
1326 hDGenQuarks[3]->Fill(dijet->Pt());
1327 hDGenQuarks[4]->Fill(dijet->Mass());
1328 }
1329 delete dijet;
1330 }
1331 // b quark info
1332 if (GenQuarks->At(i)->AbsPdgId() == 5) {
1333 hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
1334 hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());
1335 hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());
1336 if (GenLeptons->GetEntries() >= 2 &&
1337 GenLeptons->At(0)->Pt() > 20 &&
1338 GenLeptons->At(1)->Pt() > 15) {
1339 if (TMath::Abs(GenLeptons->At(0)->Eta()) < 2.5 &&
1340 TMath::Abs(GenLeptons->At(1)->Eta()) < 2.5) {
1341 hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());
1342 hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());
1343 hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());
1344 }
1345 }
1346 }
1347 // t quark info
1348 else if (GenQuarks->At(i)->AbsPdgId() == 6) {
1349 hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
1350 hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());
1351 hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());
1352 }
1353 // light quark info
1354 else {
1355 hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
1356 hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());
1357 hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());
1358 }
1359 }
1360
1361 // wbf
1362 if (GenqqHs->GetEntries() == 2) {
1363 hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
1364 GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
1365 hDGenWBF[1]->Fill(TMath::Abs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
1366 hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
1367 hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
1368 CompositeParticle *diqq = new CompositeParticle();
1369 diqq->AddDaughter(GenqqHs->At(0));
1370 diqq->AddDaughter(GenqqHs->At(1));
1371 hDGenWBF[4]->Fill(diqq->Mass());
1372 delete diqq;
1373 }
1374
1375 // bosons
1376 hDGenBosons[0]->Fill(GenBosons->GetEntries());
1377 for(UInt_t i=0; i<GenBosons->GetEntries(); i++) {
1378 hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
1379 hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
1380 hDGenBosons[3]->Fill(TMath::Min(GenBosons->At(i)->Mass(),1999.999));
1381 hDGenBosons[4]->Fill(TMath::Min(GenBosons->At(i)->Mass(),199.999));
1382 if(GenBosons->At(i)->Is(MCParticle::kW))
1383 hDGenBosons[5]->Fill(TMath::Min(GenBosons->At(i)->Mass(),199.999));
1384 if(GenBosons->At(i)->Is(MCParticle::kZ))
1385 hDGenBosons[6]->Fill(TMath::Min(GenBosons->At(i)->Mass(),199.999));
1386 hDGenBosons[7]->Fill(GenBosons->At(i)->Rapidity());
1387 if(GenBosons->At(i)->Mass() > 60 && GenBosons->At(i)->Mass() < 120){
1388 hDGenBosons[8] ->Fill(GenBosons->At(i)->Pt());
1389 hDGenBosons[9] ->Fill(GenBosons->At(i)->Eta());
1390 hDGenBosons[10]->Fill(GenBosons->At(i)->Rapidity());
1391 }
1392 }
1393 hDGenBosons[11]->Fill(TMath::Min((double)(sumV[0] + 4*sumV[1]),12.4999));
1394
1395 // photons
1396 hDGenPhotons[0]->Fill(GenPhotons->GetEntries());
1397 for(UInt_t i=0; i<GenPhotons->GetEntries(); i++) {
1398 hDGenPhotons[1]->Fill(GenPhotons->At(i)->Pt());
1399 hDGenPhotons[2]->Fill(GenPhotons->At(i)->Eta());
1400 }
1401
1402 // Rad photons
1403 hDGenRadPhotons[0]->Fill(GenRadPhotons->GetEntries());
1404 for(UInt_t i=0; i<GenRadPhotons->GetEntries(); i++) {
1405 hDGenRadPhotons[1]->Fill(TMath::Min(GenRadPhotons->At(i)->Pt(),199.999));
1406 hDGenRadPhotons[2]->Fill(TMath::Min(GenRadPhotons->At(i)->AbsEta(),4.999));
1407 hDGenRadPhotons[3]->Fill(TMath::Min((double)GenRadPhotons->At(i)->DistinctMother()->Status(),
1408 19.499));
1409 hDGenRadPhotons[4]->Fill(GenRadPhotons->At(i)->IsGenerated() +
1410 2*GenRadPhotons->At(i)->IsSimulated());
1411 if(GenRadPhotons->At(i)->DistinctMother()){
1412 hDGenRadPhotons[5]->Fill(TMath::Min(
1413 MathUtils::DeltaR(*GenRadPhotons->At(i),
1414 *GenRadPhotons->At(i)->DistinctMother()),
1415 4.999));
1416 CompositeParticle *object = new CompositeParticle();
1417 object->AddDaughter(GenRadPhotons->At(i));
1418 object->AddDaughter(GenRadPhotons->At(i)->DistinctMother());
1419 hDGenRadPhotons[6]->Fill(TMath::Min(object->Mass(),99.999));
1420 delete object;
1421 }
1422 Int_t Mother = 0;
1423 if(GenRadPhotons->At(i)->DistinctMother() &&
1424 GenRadPhotons->At(i)->DistinctMother()->Is(MCParticle::kMu)) Mother = 1;
1425 hDGenRadPhotons[7]->Fill(Mother);
1426 }
1427
1428 // ISR photons
1429 hDGenISRPhotons[0]->Fill(GenISRPhotons->GetEntries());
1430 for(UInt_t i=0; i<GenISRPhotons->GetEntries(); i++) {
1431 hDGenISRPhotons[1]->Fill(TMath::Min(GenISRPhotons->At(i)->Pt(),199.999));
1432 hDGenISRPhotons[2]->Fill(TMath::Min(GenISRPhotons->At(i)->AbsEta(),4.999));
1433 hDGenISRPhotons[3]->Fill(TMath::Min((Double_t)GenISRPhotons->At(i)->DistinctMother()->Status(),
1434 19.499));
1435 hDGenISRPhotons[4]->Fill(GenISRPhotons->At(i)->IsGenerated() +
1436 2*GenISRPhotons->At(i)->IsSimulated());
1437 if(GenISRPhotons->At(i)->DistinctMother()){
1438 hDGenISRPhotons[5]->Fill(TMath::Min(
1439 MathUtils::DeltaR(*GenISRPhotons->At(i),
1440 *GenISRPhotons->At(i)->DistinctMother()),4.999));
1441 CompositeParticle *object = new CompositeParticle();
1442 object->AddDaughter(GenISRPhotons->At(i));
1443 object->AddDaughter(GenISRPhotons->At(i)->DistinctMother());
1444 hDGenISRPhotons[6]->Fill(TMath::Min(object->Mass(),99.999));
1445 delete object;
1446 }
1447 }
1448 } // Fill histograms
1449
1450 // Apply ISR+Rad filter (but filling all histograms)
1451 if (fApplyISRFilter == kTRUE &&
1452 (GenISRPhotons->GetEntries() > 0 || GenRadPhotons->GetEntries() > 0)) {
1453 SkipEvent();
1454 }
1455 }
1456
1457 //--------------------------------------------------------------------------------------------------
1458 void GeneratorMod::SlaveBegin()
1459 {
1460 // Book branch and histograms if wanted.
1461
1462 if(fIsData == kFALSE){
1463 ReqEventObject(fMCPartName, fParticles, kTRUE);
1464 }
1465
1466 // Publish Arrays For the Output Module
1467 PublishObj(fGenLeptons);
1468 PublishObj(fGenAllLeptons);
1469 PublishObj(fGenTaus);
1470 PublishObj(fGenNeutrinos);
1471 PublishObj(fGenQuarks);
1472 PublishObj(fGenqqHs);
1473 PublishObj(fGenBosons);
1474 PublishObj(fGenPhotons);
1475 PublishObj(fGenRadPhotons);
1476 PublishObj(fGenISRPhotons);
1477
1478 // fill histograms
1479 if (GetFillHist()) {
1480 // MET
1481 AddTH1(hDGenMet[0],"hDGenMet_0","Gen MET Pt;p_{t} [GeV];#",200,0,200);
1482 AddTH1(hDGenMet[1],"hDGenMet_1","Gen MET Px;p_{x} [GeV];#",400,-200,200);
1483 AddTH1(hDGenMet[2],"hDGenMet_2","Gen MET Py;p_{y} [GeV];#",400,-200,200);
1484 AddTH1(hDGenMet[3],"hDGenMet_3","Gen MET Pz;p_{z} [GeV];#",400,-1000,1000);
1485
1486 // pt min for leptons from W/Z
1487 AddTH1(hDGenPtMin, "hDGenPtMin","Pt min leptons from W/Z;p_{t} [GeV];#",200,0.0,200.0);
1488
1489 // leptons from W
1490 AddTH1(hDGenLeptons[0], "hDGenLeptons_0",
1491 "Number of leptons from W/Z;N_{leptons};#",10,-0.5,9.5);
1492 AddTH1(hDGenLeptons[1], "hDGenLeptons_1","Pt leptons from W/Z;p_{t} [GeV];#",100,0.0,200.0);
1493 AddTH1(hDGenLeptons[2], "hDGenLeptons_2","Eta leptons from W/Z;#eta;#",50,0.0,5.0);
1494 AddTH1(hDGenLeptons[3], "hDGenLeptons_3","Phi leptons from W/Z;#phi;#",90,0.0,180.0);
1495 AddTH1(hDGenLeptons[4], "hDGenLeptons_4","Dilepton mass from W/Z;m_{ll};#",1000,0.0,1000.0);
1496 AddTH1(hDGenLeptons[5], "hDGenLeptons_5","Eta Max for 2 lepton case;#eta;#",50,0.0,5.0);
1497 AddTH1(hDGenLeptons[6], "hDGenLeptons_6","Eta Min for 2 lepton case;#eta;#",50,0.0,5.0);
1498 AddTH1(hDGenLeptons[7], "hDGenLeptons_7",
1499 "Pt Max for 2 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1500 AddTH1(hDGenLeptons[8], "hDGenLeptons_8",
1501 "Pt Min for 2 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1502 AddTH1(hDGenLeptons[9], "hDGenLeptons_9",
1503 "DiLepton mass for 2 lepton case;p_{t} [GeV];#",1000,0.0,1000.0);
1504 AddTH1(hDGenLeptons[10],"hDGenLeptons_10",
1505 "Delta Phi ll for 2 lepton case;#Delta#phi_{ll};#",90,0.0,180.0);
1506 AddTH1(hDGenLeptons[11],"hDGenLeptons_11","Delta R ll;#Delta R_{ll};#",100,0.0,5.0);
1507 AddTH1(hDGenLeptons[12],"hDGenLeptons_12",
1508 "Pt Max for 3 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1509 AddTH1(hDGenLeptons[13],"hDGenLeptons_13",
1510 "Pt 2nd for 3 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1511 AddTH1(hDGenLeptons[14],"hDGenLeptons_14",
1512 "Pt Min for 3 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1513 AddTH1(hDGenLeptons[15],"hDGenLeptons_15",
1514 "Dilepton mass for 3 lepton case;m_{ll};#",1000,0.0,1000.0);
1515 AddTH1(hDGenLeptons[16],"hDGenLeptons_16",
1516 "Trilepton mass for 3 lepton case;m_{lll};#",1000,0.0,1000.0);
1517 AddTH1(hDGenLeptons[17],"hDGenLeptons_17",
1518 "Delta R Minimum between leptons for 3 lepton case;#Delta R_{ll};#",100,0.0,5.0);
1519 AddTH1(hDGenLeptons[18],"hDGenLeptons_18",
1520 "Pt Max for 4 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1521 AddTH1(hDGenLeptons[19],"hDGenLeptons_19",
1522 "Pt 2nd for 4 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1523 AddTH1(hDGenLeptons[20],"hDGenLeptons_20",
1524 "Pt 3rd for 4 lepton case;p_{t} [GeV];#",100,0.0,200.0);
1525 AddTH1(hDGenLeptons[21],"hDGenLeptons_21",
1526 "Pt 4th for 4 lepton case;#",100,0.0,200.0);
1527 AddTH1(hDGenLeptons[22],"hDGenLeptons_22",
1528 "Dilepton mass for 4 lepton case;m_{ll};#",1000,0.0,1000.0);
1529 AddTH1(hDGenLeptons[23],"hDGenLeptons_23",
1530 "Fourlepton mass for 3 lepton case;m_{llll};#",1000,0.0,1000.0);
1531 AddTH1(hDGenLeptons[24],"hDGenLeptons_24",
1532 "Delta R Minimum between leptons for 4 lepton case;#Delta R_{ll};#",100,0.0,5.0);
1533
1534 // all leptons
1535 AddTH1(hDGenAllLeptons[0], "hDGenAllLeptons_0",
1536 "Number of all leptons;N_{leptons};#",10,-0.5,9.5);
1537 AddTH1(hDGenAllLeptons[1], "hDGenAllLeptons_1","Pt all leptons;p_{t} [GeV];#",400,0.0,200.0);
1538 AddTH1(hDGenAllLeptons[2], "hDGenAllLeptons_2","Eta all leptons;#eta;#",50,0.0,5.0);
1539 AddTH1(hDGenAllLeptons[3], "hDGenAllLeptons_3","Phi all leptons;#phi;#",90,0.0,180.0);
1540 AddTH1(hDGenAllLeptons[4], "hDGenAllLeptons_4","Pt second lepton;p_{t} [GeV];#",400,0.0,200.0);
1541
1542 // taus
1543 AddTH1(hDGenTaus[0], "hDGenTaus_0","Number of taus;N_{tau};#",10,-0.5,9.5);
1544 AddTH1(hDGenTaus[1], "hDGenTaus_1","Pt taus;p_{t} [GeV];#",100,0.0,200.0);
1545 AddTH1(hDGenTaus[2], "hDGenTaus_2","Eta taus;#eta;#",50,0.0,5.0);
1546 AddTH1(hDGenTaus[3], "hDGenTaus_3","Phi taus;#phi;#",90,0.0,180.0);
1547
1548 // neutrinos
1549 AddTH1(hDGenNeutrinos[0], "hDGenNeutrinos_0","Number of neutrinos;N_{#nu};#",10,-0.5,9.5);
1550 AddTH1(hDGenNeutrinos[1], "hDGenNeutrinos_1","Pt neutrinos;p_{t} [GeV];#",100,0.0,200.0);
1551 AddTH1(hDGenNeutrinos[2], "hDGenNeutrinos_2","Eta neutrinos;#eta;#",50,0.0,5.0);
1552 AddTH1(hDGenNeutrinos[3], "hDGenNeutrinos_3","Phi neutrinos;#phi;#",90,0.0,180.0);
1553
1554 // quarks
1555 AddTH1(hDGenQuarks[0], "hDGenQuarks_0", "Number of quarks;N_{quarks};#",10,-0.5,9.5);
1556 AddTH1(hDGenQuarks[1], "hDGenQuarks_1", "dijet pt for quarks;p_{t} [GeV];#",200,0.0,400.);
1557 AddTH1(hDGenQuarks[2], "hDGenQuarks_2", "dijet mass for quarks;m_{jj};#",2000,0.0,2000.);
1558 AddTH1(hDGenQuarks[3], "hDGenQuarks_3",
1559 "dijet pt for quarks with |eta|<2.5;p_{t} [GeV];#",200,0.0,400.);
1560 AddTH1(hDGenQuarks[4], "hDGenQuarks_4",
1561 "dijet mass for quarks with |eta|<2.5;m_{jj};#",2000,0.0,2000.);
1562 AddTH1(hDGenQuarks[5], "hDGenQuarks_5", "Pt for b quarks;p_{t} [GeV];#",200,0.0,400.);
1563 AddTH1(hDGenQuarks[6], "hDGenQuarks_6", "Eta for b quarks;#eta;#",200,-10.0,10.);
1564 AddTH1(hDGenQuarks[7], "hDGenQuarks_7",
1565 "Phi for b quarks;#phi;#",200,-TMath::Pi(),TMath::Pi());
1566 AddTH1(hDGenQuarks[8], "hDGenQuarks_8",
1567 "Pt for b quarks with |eta|<2.5;p_{t} [GeV];#",200,0.0,400.);
1568 AddTH1(hDGenQuarks[9], "hDGenQuarks_9",
1569 "Eta for b quarks with |eta|<2.5;#eta;#",200,-10.0,10.);
1570 AddTH1(hDGenQuarks[10],"hDGenQuarks_10",
1571 "Phi for b quarks with |eta|<2.5;#phi;#",200,-TMath::Pi(),TMath::Pi());
1572 AddTH1(hDGenQuarks[11],"hDGenQuarks_11","Pt for t quarks;p_{t} [GeV];#",200,0.0,400.);
1573 AddTH1(hDGenQuarks[12],"hDGenQuarks_12","Eta for t quarks;#eta;#",200,-10.0,10.);
1574 AddTH1(hDGenQuarks[13],"hDGenQuarks_13",
1575 "Phi for t quarks;#phi;#",200,-TMath::Pi(),TMath::Pi());
1576 AddTH1(hDGenQuarks[14],"hDGenQuarks_14","Pt for light quarks;p_{t} [GeV];#",200,0.0,400.);
1577 AddTH1(hDGenQuarks[15],"hDGenQuarks_15","Eta for light quarks;#eta;#",200,-10.0,10.);
1578 AddTH1(hDGenQuarks[16],"hDGenQuarks_16",
1579 "Phi for light quarks;#phi;#",200,-TMath::Pi(),TMath::Pi());
1580
1581 // qqH
1582 AddTH1(hDGenWBF[0], "hDGenWBF_0",
1583 "Delta Phi jj for WBF quarks;#Delta Phi_{jj};#",90,0.0,180.);
1584 AddTH1(hDGenWBF[1], "hDGenWBF_1",
1585 "Delta Eta jj for WBF quarks;#Delta #eta_{jj};#",100,0.0,10.);
1586 AddTH1(hDGenWBF[2], "hDGenWBF_2",
1587 "Pt max for WBF quarks;p_{t} [GeV];#",200,0.0,400.);
1588 AddTH1(hDGenWBF[3], "hDGenWBF_3",
1589 "Pt min for WBF quarks;p_{t} [GeV];#",200,0.0,400.);
1590 AddTH1(hDGenWBF[4], "hDGenWBF_4",
1591 "dijet mass for WBF quarks;m_{jj};#",200,0.0,4000.);
1592
1593 // bosons
1594 AddTH1(hDGenBosons[0], "hDGenBosons_0", "Number of bosons;N_{bosons};#",10,-0.5,9.5);
1595 AddTH1(hDGenBosons[1], "hDGenBosons_1", "Pt of bosons;p_{t} [GeV];#",200,0.0,400.0);
1596 AddTH1(hDGenBosons[2], "hDGenBosons_2", "Eta of bosons;#eta;#",100,-10.0,10.0);
1597 AddTH1(hDGenBosons[3], "hDGenBosons_3", "Mass of bosons;Mass;#",2000,0.0,2000.0);
1598 AddTH1(hDGenBosons[4], "hDGenBosons_4", "Mass of bosons;m_{V};#",200,0.0,200.0);
1599 AddTH1(hDGenBosons[5], "hDGenBosons_5", "Mass of W bosons;m_{W};#",200,0.0,200.0);
1600 AddTH1(hDGenBosons[6], "hDGenBosons_6", "Mass of Z bosons;m_{Z};#",200,0.0,200.0);
1601 AddTH1(hDGenBosons[7], "hDGenBosons_7", "Rapidity of bosons;rapidity;#",100,-10.0,10.0);
1602 AddTH1(hDGenBosons[8], "hDGenBosons_8", "Pt of bosons;p_{t} [GeV];#",200,0.0,400.0);
1603 AddTH1(hDGenBosons[9], "hDGenBosons_9", "Eta of bosons;#eta;#",100,-10.0,10.0);
1604 AddTH1(hDGenBosons[10],"hDGenBosons_10","Rapidity of bosons;rapidity;#",100,-10.0,10.0);
1605 AddTH1(hDGenBosons[11],"hDGenBosons_11","Number of W bosons + 4 * Z bosons;Number;#",13,-0.5,12.5);
1606
1607 // photons
1608 AddTH1(hDGenPhotons[0], "hDGenPhotons_0", "Number of photons;N_{photons};#",10,-0.5,9.5);
1609 AddTH1(hDGenPhotons[1], "hDGenPhotons_1", "Pt of photons;p_{t} [GeV];#",200,0.0,400.0);
1610 AddTH1(hDGenPhotons[2], "hDGenPhotons_2", "Eta of photons;#eta;#",100,-5.0,5.0);
1611
1612 // rad photons
1613 AddTH1(hDGenRadPhotons[0], "hDGenRadPhotons_0",
1614 "Number of radiative photons;N_{photons};#",10,-0.5,9.5);
1615 AddTH1(hDGenRadPhotons[1], "hDGenRadPhotons_1",
1616 "Pt of radiative photons;p_{t} [GeV];#",400,0.0,200.0);
1617 AddTH1(hDGenRadPhotons[2], "hDGenRadPhotons_2",
1618 "Eta of radiative photons;#eta;#",100,0.0,5.0);
1619 AddTH1(hDGenRadPhotons[3], "hDGenRadPhotons_3",
1620 "Status of mother of radiative photons;Status;#",20,-0.5,19.5);
1621 AddTH1(hDGenRadPhotons[4], "hDGenRadPhotons_4",
1622 "IsGenerated+2*IsSimulated of radiative photons;IsGenerated+2*IsSimulated;#",
1623 4,-0.5,3.5);
1624 AddTH1(hDGenRadPhotons[5], "hDGenRadPhotons_5",
1625 "Delta R between photon and mother of radiative photons;#Delta R;#",500,0.0,5.0);
1626 AddTH1(hDGenRadPhotons[6], "hDGenRadPhotons_6",
1627 "Mass photon-photon mother;Mass;#",500,0.0,100.0);
1628 AddTH1(hDGenRadPhotons[7], "hDGenRadPhotons_7",
1629 "Number of radiative photon with muon as a mother;Status;#",2,-0.5,1.5);
1630
1631 // ISR photons
1632 AddTH1(hDGenISRPhotons[0], "hDGenISRPhotons_0",
1633 "Number of ISR photons;N_{photons};#",10,-0.5,9.5);
1634 AddTH1(hDGenISRPhotons[1], "hDGenISRPhotons_1",
1635 "Pt of ISR photons;p_{t} [GeV];#",400,0.0,200.0);
1636 AddTH1(hDGenISRPhotons[2], "hDGenISRPhotons_2",
1637 "Eta of ISR photons;#eta;#",100,0.0,5.0);
1638 AddTH1(hDGenISRPhotons[3], "hDGenISRPhotons_3",
1639 "Status of mother of radiative photons;#eta;#",20,-0.5,19.5);
1640 AddTH1(hDGenISRPhotons[4], "hDGenISRPhotons_4",
1641 "IsGenerated+2*IsSimulated of radiative photons;IsGenerated+2*IsSimulated;#",
1642 4,-0.5,3.5);
1643 AddTH1(hDGenISRPhotons[5], "hDGenISRPhotons_5",
1644 "Delta R between photon and mother of ISR photons;#Delta R;#",500,0.0,5.0);
1645 AddTH1(hDGenISRPhotons[6], "hDGenISRPhotons_6",
1646 "Mass photon-photon mother;Mass;#",500,0.0,100.0);
1647
1648 // auxiliar for MG studies
1649 AddTH1(hDVMass[0], "hDVMass_0", "Mass of munu candidates ;Mass [GeV];#",200,0.,200.);
1650 AddTH1(hDVMass[1], "hDVMass_1", "Mass of elnu candidates ;Mass [GeV];#",200,0.,200.);
1651 AddTH1(hDVMass[2], "hDVMass_2", "Mass of taunu candidates ;Mass [GeV];#",200,0.,200.);
1652 AddTH1(hDVMass[3], "hDVMass_3", "Mass of mumu candidates ;Mass [GeV];#",200,0.,200.);
1653 AddTH1(hDVMass[4], "hDVMass_4", "Mass of ee candidates ;Mass [GeV];#",200,0.,200.);
1654 AddTH1(hDVMass[5], "hDVMass_5", "Mass of tautau candidates;Mass [GeV];#",200,0.,200.);
1655 AddTH1(hDVMass[6], "hDVMass_6", "Mass of numunumu candidates;Mass [GeV];#",200,0.,200.);
1656 AddTH1(hDVMass[7], "hDVMass_7", "Mass of nuenue candidates;Mass [GeV];#",200,0.,200.);
1657 AddTH1(hDVMass[8], "hDVMass_8", "Mass of nutaunutau candidates;Mass [GeV];#",200,0.,200.);
1658 AddTH1(hDVMass[9], "hDVMass_9",
1659 "Mass of munu candidates for t events ;Mass [GeV];#",200,0.,200.);
1660 AddTH1(hDVMass[10],"hDVMass_10",
1661 "Mass of elnu candidates for t events ;Mass [GeV];#",200,0.,200.);
1662 AddTH1(hDVMass[11],"hDVMass_11",
1663 "Mass of taunu candidates for t events;Mass [GeV];#",200,0.,200.);
1664 AddTH1(hDVMass[12],"hDVMass_12",
1665 "Mass of qq candidates for t events;Mass [GeV];#",200,0.,200.);
1666
1667 // Special study about VVjets
1668 AddTH1(hDVVMass[0], "hDVVMass_0", "Mass of munu for WW events;Mass [GeV];#",200,0.,200.);
1669 AddTH1(hDVVMass[1], "hDVVMass_1", "Mass of munu WZ events;Mass [GeV];#",200,0.,200.);
1670 AddTH1(hDVVMass[2], "hDVVMass_2", "Mass of elnu WW events;Mass [GeV];#",200,0.,200.);
1671 AddTH1(hDVVMass[3], "hDVVMass_3", "Mass of elnu WZ events;Mass [GeV];#",200,0.,200.);
1672 AddTH1(hDVVMass[4], "hDVVMass_4", "Mass of taunu WW events;Mass [GeV];#",200,0.,200.);
1673 AddTH1(hDVVMass[5], "hDVVMass_5", "Mass of taunu WZ events;Mass [GeV];#",200,0.,200.);
1674 AddTH1(hDVVMass[6], "hDVVMass_6", "Mass of mumu WZ events;Mass [GeV];#",200,0.,200.);
1675 AddTH1(hDVVMass[7], "hDVVMass_7", "Mass of mumu ZZ events;Mass [GeV];#",200,0.,200.);
1676 AddTH1(hDVVMass[8], "hDVVMass_8", "Mass of ee WZ events;Mass [GeV];#",200,0.,200.);
1677 AddTH1(hDVVMass[9], "hDVVMass_9", "Mass of ee ZZ events;Mass [GeV];#",200,0.,200.);
1678 AddTH1(hDVVMass[10],"hDVVMass_10","Mass of tautau WZ events;Mass [GeV];#",200,0.,200.);
1679 AddTH1(hDVVMass[11],"hDVVMass_11","Mass of tautau ZZ events;Mass [GeV];#",200,0.,200.);
1680 AddTH1(hDVVMass[12],"hDVVMass_12","Mass of numunumu WZ events;Mass [GeV];#",200,0.,200.);
1681 AddTH1(hDVVMass[13],"hDVVMass_13","Mass of numunumu ZZ events;Mass [GeV];#",200,0.,200.);
1682 AddTH1(hDVVMass[14],"hDVVMass_14","Mass of nuenue WZ events;Mass [GeV];#",200,0.,200.);
1683 AddTH1(hDVVMass[15],"hDVVMass_15","Mass of nuenue ZZ events;Mass [GeV];#",200,0.,200.);
1684 AddTH1(hDVVMass[16],"hDVVMass_16","Mass of nutaunutau WZ events;Mass [GeV];#",200,0.,200.);
1685 AddTH1(hDVVMass[17],"hDVVMass_17","Mass of nutaunutau ZZ events;Mass [GeV];#",200,0.,200.);
1686 AddTH1(hDVVMass[18],"hDVVMass_18","Ratios for WW events;Type;#",7,-0.5,6.5);
1687 AddTH1(hDVVMass[19],"hDVVMass_19","Ratios for WZ events;Type;#",10,-0.5,9.5);
1688 AddTH1(hDVVMass[20],"hDVVMass_20","Ratios for ZZ2l events;Type;#",7,-0.5,6.5);
1689 AddTH1(hDVVMass[21],"hDVVMass_21","Ratios for ZZ4l events;Type;#",16,-0.5,15.5);
1690 AddTH1(hDVVMass[22],"hDVVMass_22","Maximum mass for WW events;Mass [GeV];#",200,0.,200.);
1691 AddTH1(hDVVMass[23],"hDVVMass_23","Minimum mass for WW events;Mass [GeV];#",200,0.,200.);
1692 AddTH1(hDVVMass[24],"hDVVMass_24","Maximum mass for WZ events;Mass [GeV];#",200,0.,200.);
1693 AddTH1(hDVVMass[25],"hDVVMass_25","Minimum mass for WZ events;Mass [GeV];#",200,0.,200.);
1694 AddTH1(hDVVMass[26],"hDVVMass_26","Maximum mass for ZZ events;Mass [GeV];#",200,0.,200.);
1695 AddTH1(hDVVMass[27],"hDVVMass_27","Minimum mass for ZZ events;Mass [GeV];#",200,0.,200.);
1696 }
1697 }