Mods/src/GeneratorMod.cc

// $Id: GeneratorMod.cc,v 1.4 2008/11/11 21:22:54 ceballos Exp $

#include "MitPhysics/Mods/interface/GeneratorMod.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitAna/DataCont/interface/ObjArray.h"
#include "MitAna/DataUtil/interface/Debug.h"
#include "MitCommon/MathTools/interface/MathUtils.h"
#include <TH1D.h>
#include <TH2D.h>

using namespace mithep;

ClassImp(mithep::GeneratorMod)

//--------------------------------------------------------------------------------------------------
  GeneratorMod::GeneratorMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fFillHist(kFALSE),
  fMCPartName(Names::gkMCPartBrn),
  fMCLeptonsName(Names::gkMCLeptonsName),
  fMCAllLeptonsName(Names::gkMCAllLeptonsName),
  fMCTausName(Names::gkMCTausName),
  fMCNeutrinosName(Names::gkMCNeutrinosName),
  fMCQuarksName(Names::gkMCQuarksName),
  fMCqqHsName(Names::gkMCqqHsName),
  fMCBosonsName(Names::gkMCBosonsName),
  fParticles(0),
  fNEventsProcessed(0)
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void GeneratorMod::Process()
{
  // Process entries of the tree.

  // these arrays will be filled in the loop of particles
  ObjArray<MCParticle> *GenLeptons    = new ObjArray<MCParticle>;
  ObjArray<MCParticle> *GenAllLeptons = new ObjArray<MCParticle>;
  ObjArray<MCParticle> *GenTaus       = new ObjArray<MCParticle>; 
  GenTaus->SetOwner(true);
  ObjArray<MCParticle> *GenNeutrinos  = new ObjArray<MCParticle>;
  ObjArray<MCParticle> *GenQuarks     = new ObjArray<MCParticle>;
  ObjArray<MCParticle> *GenqqHs       = new ObjArray<MCParticle>;
  ObjArray<MCParticle> *GenBosons     = new ObjArray<MCParticle>;


  // load MCParticle branch
  LoadBranch(fMCPartName);

  bool isqqH = kFALSE;
  for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
    MCParticle *p = fParticles->At(i);

    if (!p->IsGenerated()) continue;

    // muons/electrons from W/Z decays
    if ((p->AbsPdgId() == 11 || p->AbsPdgId() == 13) && p->Status() == 1) {
      if (p->Pt() > 3.0 && TMath::Abs(p->Eta()) < 3.0) {
        GenAllLeptons->Add(p);
      }
      bool isGoodLepton = kFALSE;
      MCParticle *pm = p;
      while (pm->HasMother() && isGoodLepton == kFALSE) {
        if     (pm->Mother()->AbsPdgId() == 23 || pm->Mother()->AbsPdgId() == 24) {
          GenLeptons->Add(p);
          isGoodLepton = kTRUE;
        }
        else if(pm->Mother()->AbsPdgId() == 111 || pm->Mother()->AbsPdgId() == 221) {
          // this is fake, but it is a trick to get rid of these cases
          isGoodLepton = kTRUE;
        } else {
          pm = (MCParticle *)pm->Mother();
        }
      }
    }

    // taus
    else if(p->AbsPdgId() == 16 && p->Status() == 1) { 

      if (p->DistinctMother()) {
        MCParticle *pm = (mithep::MCParticle*)p->DistinctMother();
        if(pm->AbsPdgId() == 15) {
          MCParticle *pm_f = new MCParticle(*pm);
          pm_f->SetMom(pm->Px()-p->Px(), pm->Py()-p->Py(),
                       pm->Pz()-p->Pz(), pm->E()-p->E());
          GenTaus->Add(pm_f);
        }
      }
    }

    // neutrinos
    else if (p->Status() == 1 &&
            (p->AbsPdgId() == 12 || p->AbsPdgId() == 14 || p->AbsPdgId() == 16)) {
      GenNeutrinos->Add(p);
    }

    // quarks from W/Z decays or top particles
    else if (p->AbsPdgId() >=1 && p->AbsPdgId() <=6 && p->HasMother()) {
      if(p->Mother()->AbsPdgId() == 23 || p->Mother()->AbsPdgId() == 24 ||
         p->AbsPdgId() == 6 || p->Mother()->AbsPdgId() == 6) {
        GenQuarks->Add(p);
      }
    }

    // qqH, information about the forward jets
    else if(isqqH == kFALSE && p->AbsPdgId() == 25) {
      isqqH = kTRUE;
      MCParticle *pq1 = fParticles->At(i-1);
      MCParticle *pq2 = fParticles->At(i-2);

      if(pq1->HasMother() && pq2->HasMother() &&
         pq1->Mother()->PdgId() == p->Mother()->PdgId() &&
         pq2->Mother()->PdgId() == p->Mother()->PdgId() &&
         pq1->AbsPdgId() < 7 &&  pq2->AbsPdgId() < 7 &&
         pq1->AbsPdgId() > 0 &&  pq2->AbsPdgId() > 0) {
        GenqqHs->Add(pq1);
        GenqqHs->Add(pq2);
      }
    }

    // information about bosons: W, Z, h, Z', W', H0, A0, H+
    else if (p->Status() == 2 &&
            (p->AbsPdgId() == 23 || p->AbsPdgId() == 24 || p->AbsPdgId() == 25 ||
             p->AbsPdgId() == 32 || p->AbsPdgId() == 34 ||
             p->AbsPdgId() == 35 || p->AbsPdgId() == 36 || p->AbsPdgId() == 37)) {
      GenBosons->Add(p);
    }
  }

  // save Objects for Other Modules to use
  AddObjThisEvt(GenLeptons,   fMCLeptonsName.Data());  
  AddObjThisEvt(GenAllLeptons,fMCAllLeptonsName.Data());  
  AddObjThisEvt(GenTaus,      fMCTausName.Data());  
  AddObjThisEvt(GenNeutrinos, fMCNeutrinosName.Data());  
  AddObjThisEvt(GenQuarks,    fMCQuarksName.Data());  
  AddObjThisEvt(GenqqHs,      fMCqqHsName.Data());  
  AddObjThisEvt(GenBosons,    fMCBosonsName.Data());
  
  // fill histograms if requested
  if (fFillHist) {

    // Leptons
    hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
    Int_t   idxMaxLep[2] = {-1, -1};
    Double_t ptMaxLep[2] = {-1, -1};
    for(UInt_t i=0; i<GenLeptons->GetEntries(); i++) {
      hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
      hDGenLeptons[2]->Fill(TMath::Abs(GenLeptons->At(i)->Eta()));
      hDGenLeptons[3]->Fill(GenLeptons->At(i)->PhiDeg());
      for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++) {
        CompositeParticle *dilepton = new CompositeParticle();
        dilepton->AddDaughter(GenLeptons->At(i));
        dilepton->AddDaughter(GenLeptons->At(j));
        hDGenLeptons[4]->Fill(dilepton->Mass());
        delete dilepton;
      }
      // Selecting the two highest Pt leptons
      if     (GenLeptons->At(i)->Pt() > ptMaxLep[0]) {
        ptMaxLep[1]  = ptMaxLep[0];
        idxMaxLep[1] = idxMaxLep[0];
        ptMaxLep[0]  = GenLeptons->At(i)->Pt();
        idxMaxLep[0] = i;
      }
      else if (GenLeptons->At(i)->Pt() > ptMaxLep[1]) {
        ptMaxLep[1]  = GenLeptons->At(i)->Pt();
        idxMaxLep[1] = i;
      }
    }
    // Looking at events with at least two leptons
    if (ptMaxLep[0] > 0 && ptMaxLep[1] > 0) {
      hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()),
                                                  TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta())),
                                       4.999));
      hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()),
                                                  TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta())),
                                       4.999));
      if (TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
          TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5) {
        hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[0])->Pt(),199.999));
        if (GenLeptons->At(idxMaxLep[0])->Pt() > 20.0) {
          hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[1])->Pt(),199.999));
          if (GenLeptons->At(idxMaxLep[1])->Pt() > 10.0) {
            CompositeParticle *dilepton = new CompositeParticle();
            dilepton->AddDaughter(GenLeptons->At(idxMaxLep[0]));
            dilepton->AddDaughter(GenLeptons->At(idxMaxLep[1]));
            hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
            hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(idxMaxLep[0])->Phi(),
                                                       GenLeptons->At(idxMaxLep[1])->Phi())
                                                       * 180./ TMath::Pi());
            delete dilepton;
          }
        }
      }
    }

    // All Leptons
    hDGenAllLeptons[0]->Fill(GenAllLeptons->GetEntries());
    for(UInt_t i=0; i<GenAllLeptons->GetEntries(); i++) {
      hDGenAllLeptons[1]->Fill(GenAllLeptons->At(i)->Pt());
      hDGenAllLeptons[2]->Fill(GenAllLeptons->At(i)->Eta());
      hDGenAllLeptons[3]->Fill(GenAllLeptons->At(i)->PhiDeg());
    }

    // Taus
    hDGenTaus[0]->Fill(GenTaus->GetEntries());
    for(UInt_t i=0; i<GenTaus->GetEntries(); i++) {
      hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
      hDGenTaus[2]->Fill(GenTaus->At(i)->Eta());
      hDGenTaus[3]->Fill(GenTaus->At(i)->PhiDeg());
    }

    // Neutrinos
    hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
    CompositeParticle *neutrinoTotal = new CompositeParticle();
    for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++) {
      if (GenNeutrinos->At(i)->HasMother())
        neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
    }
    if (GenNeutrinos->GetEntries() > 0) {
      hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
      hDGenNeutrinos[2]->Fill(neutrinoTotal->Eta());
      hDGenNeutrinos[3]->Fill(neutrinoTotal->PhiDeg());    
    }
    delete neutrinoTotal;
    
    // Quarks
    hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
    for(UInt_t i=0; i<GenQuarks->GetEntries(); i++) {
      for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++) {
        CompositeParticle *dijet = new CompositeParticle();
        dijet->AddDaughter(GenQuarks->At(i));
        dijet->AddDaughter(GenQuarks->At(j));
        hDGenQuarks[1]->Fill(dijet->Pt());
        hDGenQuarks[2]->Fill(dijet->Mass());
        if (TMath::Abs(GenQuarks->At(i)->Eta()) < 2.5 && 
            TMath::Abs(GenQuarks->At(j)->Eta()) < 2.5) {
          hDGenQuarks[3]->Fill(dijet->Pt());
          hDGenQuarks[4]->Fill(dijet->Mass());
        }
        delete dijet;
      }
      // b quark info
      if    (GenQuarks->At(i)->AbsPdgId() == 5) {
        hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
        hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());      
        hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());      
        if (ptMaxLep[0] > 20 && ptMaxLep[1] > 15) {
          if (TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
              TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5) {
            hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());       
            hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());      
            hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());     
          }
        }
      }
      // t quark info
      else if (GenQuarks->At(i)->AbsPdgId() == 6) {
        hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
        hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());      
        hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());      
      }
      // light quark info
      else {
        hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
        hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());      
        hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());      
      }
    }

    // WBF
    if (GenqqHs->GetEntries() == 2) {
      hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
                                            GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
      hDGenWBF[1]->Fill(TMath::Abs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
      hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
      hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
      CompositeParticle *diqq = new CompositeParticle();
      diqq->AddDaughter(GenqqHs->At(0));
      diqq->AddDaughter(GenqqHs->At(1));
      hDGenWBF[4]->Fill(diqq->Mass());
      delete diqq;
    }

    // Bosons
    hDGenBosons[0]->Fill(GenBosons->GetEntries());
    for(UInt_t i=0; i<GenBosons->GetEntries(); i++) {
      hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
      hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
      hDGenBosons[3]->Fill(GenBosons->At(i)->Mass());
    } 
  }
}

//--------------------------------------------------------------------------------------------------
void GeneratorMod::SlaveBegin()
{
  // Book branch and histograms if wanted.

  ReqBranch(fMCPartName, fParticles);

  // fill histograms
  if (fFillHist == kTRUE) {
    char sb[1024];

    // Leptons
    sprintf(sb,"hDGenLeptons_%d", 0);  hDGenLeptons[0]  = new TH1D(sb,sb,10,-0.5,9.5); 
    sprintf(sb,"hDGenLeptons_%d", 1);  hDGenLeptons[1]  = new TH1D(sb,sb,100,0.0,200.0); 
    sprintf(sb,"hDGenLeptons_%d", 2);  hDGenLeptons[2]  = new TH1D(sb,sb,50,0.0,5.0); 
    sprintf(sb,"hDGenLeptons_%d", 3);  hDGenLeptons[3]  = new TH1D(sb,sb,90,0.0,180.0); 
    sprintf(sb,"hDGenLeptons_%d", 4);  hDGenLeptons[4]  = new TH1D(sb,sb,1000,0.0,1000.0); 
    sprintf(sb,"hDGenLeptons_%d", 5);  hDGenLeptons[5]  = new TH1D(sb,sb,50,0.0,5.0); 
    sprintf(sb,"hDGenLeptons_%d", 6);  hDGenLeptons[6]  = new TH1D(sb,sb,50,0.0,5.0); 
    sprintf(sb,"hDGenLeptons_%d", 7);  hDGenLeptons[7]  = new TH1D(sb,sb,100,0.0,200.0); 
    sprintf(sb,"hDGenLeptons_%d", 8);  hDGenLeptons[8]  = new TH1D(sb,sb,100,0.0,200.0); 
    sprintf(sb,"hDGenLeptons_%d", 9);  hDGenLeptons[9]  = new TH1D(sb,sb,1000,0.0,1000.0); 
    sprintf(sb,"hDGenLeptons_%d",10);  hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0); 
    for(Int_t i=0; i<11; i++) AddOutput(hDGenLeptons[i]);

    // AllLeptons
    sprintf(sb,"hDGenAllLeptons_%d", 0);  hDGenAllLeptons[0]  = new TH1D(sb,sb,10,-0.5,9.5); 
    sprintf(sb,"hDGenAllLeptons_%d", 1);  hDGenAllLeptons[1]  = new TH1D(sb,sb,100,0.0,200.0); 
    sprintf(sb,"hDGenAllLeptons_%d", 2);  hDGenAllLeptons[2]  = new TH1D(sb,sb,60,-3.0,3.0); 
    sprintf(sb,"hDGenAllLeptons_%d", 3);  hDGenAllLeptons[3]  = new TH1D(sb,sb,90,0.0,180.0); 
    for(Int_t i=0; i<4; i++) AddOutput(hDGenAllLeptons[i]);

    // Taus
    sprintf(sb,"hDGenTaus_%d", 0);  hDGenTaus[0]  = new TH1D(sb,sb,10,-0.5,9.5); 
    sprintf(sb,"hDGenTaus_%d", 1);  hDGenTaus[1]  = new TH1D(sb,sb,100,0.0,200.0); 
    sprintf(sb,"hDGenTaus_%d", 2);  hDGenTaus[2]  = new TH1D(sb,sb,100,-5.0,5.0); 
    sprintf(sb,"hDGenTaus_%d", 3);  hDGenTaus[3]  = new TH1D(sb,sb,90,0.0,180.0); 
    for(Int_t i=0; i<4; i++) AddOutput(hDGenTaus[i]);

    // Neutrinos
    sprintf(sb,"hDGenNeutrinos_%d", 0);  hDGenNeutrinos[0]  = new TH1D(sb,sb,10,-0.5,9.5); 
    sprintf(sb,"hDGenNeutrinos_%d", 1);  hDGenNeutrinos[1]  = new TH1D(sb,sb,100,0.0,200.0); 
    sprintf(sb,"hDGenNeutrinos_%d", 2);  hDGenNeutrinos[2]  = new TH1D(sb,sb,100,-5.0,5.0); 
    sprintf(sb,"hDGenNeutrinos_%d", 3);  hDGenNeutrinos[3]  = new TH1D(sb,sb,90,0.0,180.0); 
    for(Int_t i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);

    // Quarks
    sprintf(sb,"hDGenQuarks_%d", 0);  hDGenQuarks[0]  = new TH1D(sb,sb,10,-0.5,9.5); 
    sprintf(sb,"hDGenQuarks_%d", 1);  hDGenQuarks[1]  = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenQuarks_%d", 2);  hDGenQuarks[2]  = new TH1D(sb,sb,2000,0.0,2000.);
    sprintf(sb,"hDGenQuarks_%d", 3);  hDGenQuarks[3]  = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenQuarks_%d", 4);  hDGenQuarks[4]  = new TH1D(sb,sb,2000,0.0,2000.);
    sprintf(sb,"hDGenQuarks_%d", 5);  hDGenQuarks[5]  = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenQuarks_%d", 6);  hDGenQuarks[6]  = new TH1D(sb,sb,200,-10.0,10.); 
    sprintf(sb,"hDGenQuarks_%d", 7);  hDGenQuarks[7]  = new TH1D(sb,sb,200,-TMath::Pi(),
                                                                 TMath::Pi()); 
    sprintf(sb,"hDGenQuarks_%d", 8);  hDGenQuarks[8]  = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenQuarks_%d", 9);  hDGenQuarks[9]  = new TH1D(sb,sb,200,-10.0,10.); 
    sprintf(sb,"hDGenQuarks_%d",10);  hDGenQuarks[10] = new TH1D(sb,sb,200,-TMath::Pi(),
                                                                 TMath::Pi()); 
    sprintf(sb,"hDGenQuarks_%d",11);  hDGenQuarks[11] = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenQuarks_%d",12);  hDGenQuarks[12] = new TH1D(sb,sb,200,-10.0,10.); 
    sprintf(sb,"hDGenQuarks_%d",13);  hDGenQuarks[13] = new TH1D(sb,sb,200,-TMath::Pi(),
                                                                 TMath::Pi()); 
    sprintf(sb,"hDGenQuarks_%d",14);  hDGenQuarks[14] = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenQuarks_%d",15);  hDGenQuarks[15] = new TH1D(sb,sb,200,-10.0,10.); 
    sprintf(sb,"hDGenQuarks_%d",16);  hDGenQuarks[16] = new TH1D(sb,sb,200,-TMath::Pi(),
                                                                 TMath::Pi()); 
    for(Int_t i=0; i<17; i++) AddOutput(hDGenQuarks[i]);

    // qqH
    sprintf(sb,"hDGenWBF_%d", 0);  hDGenWBF[0]  = new TH1D(sb,sb,90,0.0,180.);
    sprintf(sb,"hDGenWBF_%d", 1);  hDGenWBF[1]  = new TH1D(sb,sb,100,0.0,10.);   
    sprintf(sb,"hDGenWBF_%d", 2);  hDGenWBF[2]  = new TH1D(sb,sb,200,0.0,400.); 
    sprintf(sb,"hDGenWBF_%d", 3);  hDGenWBF[3]  = new TH1D(sb,sb,200,0.0,400.);
    sprintf(sb,"hDGenWBF_%d", 4);  hDGenWBF[4]  = new TH1D(sb,sb,200,0.0,4000.);
    for(Int_t i=0; i<5; i++) AddOutput(hDGenWBF[i]);

    // Bosons
    sprintf(sb,"hDGenBosons_%d", 0);  hDGenBosons[0]  = new TH1D(sb,sb,10,-0.5,9.5); 
    sprintf(sb,"hDGenBosons_%d", 1);  hDGenBosons[1]  = new TH1D(sb,sb,200,0.0,400.0); 
    sprintf(sb,"hDGenBosons_%d", 2);  hDGenBosons[2]  = new TH1D(sb,sb,100,-5.0,5.0); 
    sprintf(sb,"hDGenBosons_%d", 3);  hDGenBosons[3]  = new TH1D(sb,sb,2000,0.0,2000.0); 
    for(Int_t i=0; i<4; i++) AddOutput(hDGenBosons[i]);
  }
}
Revision:	1.5
Committed:	Wed Nov 19 15:44:35 2008 UTC (16 years, 5 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+143 -175 lines
Log Message:	Cleaned plus coding conventions
#	User	Rev	Content
1	loizides	1.5	// $Id: GeneratorMod.cc,v 1.4 2008/11/11 21:22:54 ceballos Exp $
2	loizides	1.1
3			#include "MitPhysics/Mods/interface/GeneratorMod.h"
4			#include "MitAna/DataTree/interface/Names.h"
5			#include "MitAna/DataCont/interface/ObjArray.h"
6	loizides	1.5	#include "MitAna/DataUtil/interface/Debug.h"
7	loizides	1.1	#include "MitCommon/MathTools/interface/MathUtils.h"
8			#include <TH1D.h>
9			#include <TH2D.h>
10
11			using namespace mithep;
12
13			ClassImp(mithep::GeneratorMod)
14
15			//--------------------------------------------------------------------------------------------------
16			GeneratorMod::GeneratorMod(const char name, const char title) :
17			BaseMod(name,title),
18	loizides	1.5	fFillHist(kFALSE),
19	loizides	1.1	fMCPartName(Names::gkMCPartBrn),
20			fMCLeptonsName(Names::gkMCLeptonsName),
21	ceballos	1.3	fMCAllLeptonsName(Names::gkMCAllLeptonsName),
22	loizides	1.1	fMCTausName(Names::gkMCTausName),
23			fMCNeutrinosName(Names::gkMCNeutrinosName),
24			fMCQuarksName(Names::gkMCQuarksName),
25			fMCqqHsName(Names::gkMCqqHsName),
26			fMCBosonsName(Names::gkMCBosonsName),
27			fParticles(0),
28			fNEventsProcessed(0)
29			{
30			// Constructor.
31			}
32
33			//--------------------------------------------------------------------------------------------------
34			void GeneratorMod::Process()
35			{
36	loizides	1.5	// Process entries of the tree.
37	loizides	1.1
38	loizides	1.5	// these arrays will be filled in the loop of particles
39	ceballos	1.3	ObjArray<MCParticle> *GenLeptons = new ObjArray<MCParticle>;
40			ObjArray<MCParticle> *GenAllLeptons = new ObjArray<MCParticle>;
41	loizides	1.5	ObjArray<MCParticle> *GenTaus = new ObjArray<MCParticle>;
42			GenTaus->SetOwner(true);
43	ceballos	1.3	ObjArray<MCParticle> *GenNeutrinos = new ObjArray<MCParticle>;
44			ObjArray<MCParticle> *GenQuarks = new ObjArray<MCParticle>;
45			ObjArray<MCParticle> *GenqqHs = new ObjArray<MCParticle>;
46			ObjArray<MCParticle> *GenBosons = new ObjArray<MCParticle>;
47	loizides	1.1
48	ceballos	1.4
49	loizides	1.5	// load MCParticle branch
50			LoadBranch(fMCPartName);
51
52			bool isqqH = kFALSE;
53			for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
54			MCParticle *p = fParticles->At(i);
55
56			if (!p->IsGenerated()) continue;
57
58			// muons/electrons from W/Z decays
59			if ((p->AbsPdgId() == 11 \|\| p->AbsPdgId() == 13) && p->Status() == 1) {
60			if (p->Pt() > 3.0 && TMath::Abs(p->Eta()) < 3.0) {
61			GenAllLeptons->Add(p);
62			}
63			bool isGoodLepton = kFALSE;
64			MCParticle *pm = p;
65			while (pm->HasMother() && isGoodLepton == kFALSE) {
66			if (pm->Mother()->AbsPdgId() == 23 \|\| pm->Mother()->AbsPdgId() == 24) {
67			GenLeptons->Add(p);
68			isGoodLepton = kTRUE;
69			}
70			else if(pm->Mother()->AbsPdgId() == 111 \|\| pm->Mother()->AbsPdgId() == 221) {
71			// this is fake, but it is a trick to get rid of these cases
72			isGoodLepton = kTRUE;
73			} else {
74			pm = (MCParticle *)pm->Mother();
75			}
76	loizides	1.1	}
77	loizides	1.5	}
78
79			// taus
80			else if(p->AbsPdgId() == 16 && p->Status() == 1) {
81	loizides	1.1
82	loizides	1.5	if (p->DistinctMother()) {
83			MCParticle pm = (mithep::MCParticle)p->DistinctMother();
84			if(pm->AbsPdgId() == 15) {
85			MCParticle pm_f = new MCParticle(pm);
86			pm_f->SetMom(pm->Px()-p->Px(), pm->Py()-p->Py(),
87			pm->Pz()-p->Pz(), pm->E()-p->E());
88			GenTaus->Add(pm_f);
89			}
90	loizides	1.1	}
91	loizides	1.5	}
92	loizides	1.1
93	loizides	1.5	// neutrinos
94			else if (p->Status() == 1 &&
95			(p->AbsPdgId() == 12 \|\| p->AbsPdgId() == 14 \|\| p->AbsPdgId() == 16)) {
96			GenNeutrinos->Add(p);
97			}
98	loizides	1.1
99	loizides	1.5	// quarks from W/Z decays or top particles
100			else if (p->AbsPdgId() >=1 && p->AbsPdgId() <=6 && p->HasMother()) {
101			if(p->Mother()->AbsPdgId() == 23 \|\| p->Mother()->AbsPdgId() == 24 \|\|
102			p->AbsPdgId() == 6 \|\| p->Mother()->AbsPdgId() == 6) {
103			GenQuarks->Add(p);
104	loizides	1.1	}
105	loizides	1.5	}
106	loizides	1.1
107	loizides	1.5	// qqH, information about the forward jets
108			else if(isqqH == kFALSE && p->AbsPdgId() == 25) {
109			isqqH = kTRUE;
110			MCParticle *pq1 = fParticles->At(i-1);
111			MCParticle *pq2 = fParticles->At(i-2);
112
113			if(pq1->HasMother() && pq2->HasMother() &&
114			pq1->Mother()->PdgId() == p->Mother()->PdgId() &&
115			pq2->Mother()->PdgId() == p->Mother()->PdgId() &&
116			pq1->AbsPdgId() < 7 && pq2->AbsPdgId() < 7 &&
117			pq1->AbsPdgId() > 0 && pq2->AbsPdgId() > 0) {
118			GenqqHs->Add(pq1);
119			GenqqHs->Add(pq2);
120	loizides	1.1	}
121	loizides	1.5	}
122	loizides	1.1
123	loizides	1.5	// information about bosons: W, Z, h, Z', W', H0, A0, H+
124			else if (p->Status() == 2 &&
125			(p->AbsPdgId() == 23 \|\| p->AbsPdgId() == 24 \|\| p->AbsPdgId() == 25 \|\|
126			p->AbsPdgId() == 32 \|\| p->AbsPdgId() == 34 \|\|
127			p->AbsPdgId() == 35 \|\| p->AbsPdgId() == 36 \|\| p->AbsPdgId() == 37)) {
128			GenBosons->Add(p);
129	loizides	1.1	}
130	loizides	1.5	}
131	loizides	1.1
132	loizides	1.5	// save Objects for Other Modules to use
133	ceballos	1.3	AddObjThisEvt(GenLeptons, fMCLeptonsName.Data());
134			AddObjThisEvt(GenAllLeptons,fMCAllLeptonsName.Data());
135			AddObjThisEvt(GenTaus, fMCTausName.Data());
136			AddObjThisEvt(GenNeutrinos, fMCNeutrinosName.Data());
137			AddObjThisEvt(GenQuarks, fMCQuarksName.Data());
138			AddObjThisEvt(GenqqHs, fMCqqHsName.Data());
139			AddObjThisEvt(GenBosons, fMCBosonsName.Data());
140	loizides	1.1
141	loizides	1.5	// fill histograms if requested
142			if (fFillHist) {
143
144	loizides	1.1	// Leptons
145			hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
146	loizides	1.5	Int_t idxMaxLep[2] = {-1, -1};
147			Double_t ptMaxLep[2] = {-1, -1};
148			for(UInt_t i=0; i<GenLeptons->GetEntries(); i++) {
149	loizides	1.1	hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
150	loizides	1.5	hDGenLeptons[2]->Fill(TMath::Abs(GenLeptons->At(i)->Eta()));
151			hDGenLeptons[3]->Fill(GenLeptons->At(i)->PhiDeg());
152			for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++) {
153	loizides	1.1	CompositeParticle *dilepton = new CompositeParticle();
154			dilepton->AddDaughter(GenLeptons->At(i));
155			dilepton->AddDaughter(GenLeptons->At(j));
156			hDGenLeptons[4]->Fill(dilepton->Mass());
157			delete dilepton;
158			}
159			// Selecting the two highest Pt leptons
160	loizides	1.5	if (GenLeptons->At(i)->Pt() > ptMaxLep[0]) {
161	loizides	1.1	ptMaxLep[1] = ptMaxLep[0];
162			idxMaxLep[1] = idxMaxLep[0];
163			ptMaxLep[0] = GenLeptons->At(i)->Pt();
164			idxMaxLep[0] = i;
165			}
166	loizides	1.5	else if (GenLeptons->At(i)->Pt() > ptMaxLep[1]) {
167	loizides	1.1	ptMaxLep[1] = GenLeptons->At(i)->Pt();
168			idxMaxLep[1] = i;
169			}
170			}
171			// Looking at events with at least two leptons
172	loizides	1.5	if (ptMaxLep[0] > 0 && ptMaxLep[1] > 0) {
173			hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()),
174			TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta())),
175			4.999));
176			hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()),
177			TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta())),
178			4.999));
179			if (TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
180			TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5) {
181	loizides	1.1	hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[0])->Pt(),199.999));
182	loizides	1.5	if (GenLeptons->At(idxMaxLep[0])->Pt() > 20.0) {
183	loizides	1.1	hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[1])->Pt(),199.999));
184	loizides	1.5	if (GenLeptons->At(idxMaxLep[1])->Pt() > 10.0) {
185	loizides	1.1	CompositeParticle *dilepton = new CompositeParticle();
186			dilepton->AddDaughter(GenLeptons->At(idxMaxLep[0]));
187			dilepton->AddDaughter(GenLeptons->At(idxMaxLep[1]));
188			hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
189			hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(idxMaxLep[0])->Phi(),
190			GenLeptons->At(idxMaxLep[1])->Phi())
191			* 180./ TMath::Pi());
192			delete dilepton;
193			}
194			}
195			}
196			}
197
198	ceballos	1.3	// All Leptons
199			hDGenAllLeptons[0]->Fill(GenAllLeptons->GetEntries());
200	loizides	1.5	for(UInt_t i=0; i<GenAllLeptons->GetEntries(); i++) {
201	ceballos	1.3	hDGenAllLeptons[1]->Fill(GenAllLeptons->At(i)->Pt());
202			hDGenAllLeptons[2]->Fill(GenAllLeptons->At(i)->Eta());
203	loizides	1.5	hDGenAllLeptons[3]->Fill(GenAllLeptons->At(i)->PhiDeg());
204	ceballos	1.3	}
205
206	loizides	1.1	// Taus
207			hDGenTaus[0]->Fill(GenTaus->GetEntries());
208	loizides	1.5	for(UInt_t i=0; i<GenTaus->GetEntries(); i++) {
209	loizides	1.1	hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
210			hDGenTaus[2]->Fill(GenTaus->At(i)->Eta());
211	loizides	1.5	hDGenTaus[3]->Fill(GenTaus->At(i)->PhiDeg());
212	loizides	1.1	}
213
214			// Neutrinos
215			hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
216			CompositeParticle *neutrinoTotal = new CompositeParticle();
217	loizides	1.5	for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++) {
218			if (GenNeutrinos->At(i)->HasMother())
219	loizides	1.1	neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
220			}
221	loizides	1.5	if (GenNeutrinos->GetEntries() > 0) {
222	loizides	1.1	hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
223			hDGenNeutrinos[2]->Fill(neutrinoTotal->Eta());
224	loizides	1.5	hDGenNeutrinos[3]->Fill(neutrinoTotal->PhiDeg());
225	loizides	1.1	}
226			delete neutrinoTotal;
227
228			// Quarks
229			hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
230	loizides	1.5	for(UInt_t i=0; i<GenQuarks->GetEntries(); i++) {
231			for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++) {
232	loizides	1.1	CompositeParticle *dijet = new CompositeParticle();
233			dijet->AddDaughter(GenQuarks->At(i));
234			dijet->AddDaughter(GenQuarks->At(j));
235			hDGenQuarks[1]->Fill(dijet->Pt());
236			hDGenQuarks[2]->Fill(dijet->Mass());
237	loizides	1.5	if (TMath::Abs(GenQuarks->At(i)->Eta()) < 2.5 &&
238			TMath::Abs(GenQuarks->At(j)->Eta()) < 2.5) {
239	loizides	1.1	hDGenQuarks[3]->Fill(dijet->Pt());
240			hDGenQuarks[4]->Fill(dijet->Mass());
241			}
242			delete dijet;
243			}
244	ceballos	1.2	// b quark info
245	loizides	1.5	if (GenQuarks->At(i)->AbsPdgId() == 5) {
246	ceballos	1.2	hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
247			hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());
248			hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());
249	loizides	1.5	if (ptMaxLep[0] > 20 && ptMaxLep[1] > 15) {
250			if (TMath::Abs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
251			TMath::Abs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5) {
252	ceballos	1.2	hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());
253			hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());
254			hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());
255			}
256			}
257			}
258			// t quark info
259	loizides	1.5	else if (GenQuarks->At(i)->AbsPdgId() == 6) {
260	ceballos	1.2	hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
261			hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());
262			hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());
263			}
264			// light quark info
265			else {
266			hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
267			hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());
268			hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());
269			}
270	loizides	1.1	}
271
272			// WBF
273	loizides	1.5	if (GenqqHs->GetEntries() == 2) {
274	loizides	1.1	hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
275			GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
276	loizides	1.5	hDGenWBF[1]->Fill(TMath::Abs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
277	loizides	1.1	hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
278			hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
279			CompositeParticle *diqq = new CompositeParticle();
280			diqq->AddDaughter(GenqqHs->At(0));
281			diqq->AddDaughter(GenqqHs->At(1));
282			hDGenWBF[4]->Fill(diqq->Mass());
283			delete diqq;
284			}
285
286			// Bosons
287			hDGenBosons[0]->Fill(GenBosons->GetEntries());
288	loizides	1.5	for(UInt_t i=0; i<GenBosons->GetEntries(); i++) {
289	loizides	1.1	hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
290			hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
291			hDGenBosons[3]->Fill(GenBosons->At(i)->Mass());
292			}
293			}
294			}
295
296			//--------------------------------------------------------------------------------------------------
297			void GeneratorMod::SlaveBegin()
298			{
299	loizides	1.5	// Book branch and histograms if wanted.
300
301	loizides	1.1	ReqBranch(fMCPartName, fParticles);
302
303	loizides	1.5	// fill histograms
304			if (fFillHist == kTRUE) {
305			char sb[1024];
306
307	ceballos	1.3	// Leptons
308	loizides	1.1	sprintf(sb,"hDGenLeptons_%d", 0); hDGenLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
309			sprintf(sb,"hDGenLeptons_%d", 1); hDGenLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
310			sprintf(sb,"hDGenLeptons_%d", 2); hDGenLeptons[2] = new TH1D(sb,sb,50,0.0,5.0);
311			sprintf(sb,"hDGenLeptons_%d", 3); hDGenLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
312			sprintf(sb,"hDGenLeptons_%d", 4); hDGenLeptons[4] = new TH1D(sb,sb,1000,0.0,1000.0);
313			sprintf(sb,"hDGenLeptons_%d", 5); hDGenLeptons[5] = new TH1D(sb,sb,50,0.0,5.0);
314			sprintf(sb,"hDGenLeptons_%d", 6); hDGenLeptons[6] = new TH1D(sb,sb,50,0.0,5.0);
315			sprintf(sb,"hDGenLeptons_%d", 7); hDGenLeptons[7] = new TH1D(sb,sb,100,0.0,200.0);
316			sprintf(sb,"hDGenLeptons_%d", 8); hDGenLeptons[8] = new TH1D(sb,sb,100,0.0,200.0);
317			sprintf(sb,"hDGenLeptons_%d", 9); hDGenLeptons[9] = new TH1D(sb,sb,1000,0.0,1000.0);
318			sprintf(sb,"hDGenLeptons_%d",10); hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0);
319	loizides	1.5	for(Int_t i=0; i<11; i++) AddOutput(hDGenLeptons[i]);
320	loizides	1.1
321	ceballos	1.3	// AllLeptons
322			sprintf(sb,"hDGenAllLeptons_%d", 0); hDGenAllLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
323			sprintf(sb,"hDGenAllLeptons_%d", 1); hDGenAllLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
324			sprintf(sb,"hDGenAllLeptons_%d", 2); hDGenAllLeptons[2] = new TH1D(sb,sb,60,-3.0,3.0);
325			sprintf(sb,"hDGenAllLeptons_%d", 3); hDGenAllLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
326	loizides	1.5	for(Int_t i=0; i<4; i++) AddOutput(hDGenAllLeptons[i]);
327	ceballos	1.3
328	loizides	1.1	// Taus
329			sprintf(sb,"hDGenTaus_%d", 0); hDGenTaus[0] = new TH1D(sb,sb,10,-0.5,9.5);
330			sprintf(sb,"hDGenTaus_%d", 1); hDGenTaus[1] = new TH1D(sb,sb,100,0.0,200.0);
331			sprintf(sb,"hDGenTaus_%d", 2); hDGenTaus[2] = new TH1D(sb,sb,100,-5.0,5.0);
332			sprintf(sb,"hDGenTaus_%d", 3); hDGenTaus[3] = new TH1D(sb,sb,90,0.0,180.0);
333	loizides	1.5	for(Int_t i=0; i<4; i++) AddOutput(hDGenTaus[i]);
334	loizides	1.1
335			// Neutrinos
336			sprintf(sb,"hDGenNeutrinos_%d", 0); hDGenNeutrinos[0] = new TH1D(sb,sb,10,-0.5,9.5);
337			sprintf(sb,"hDGenNeutrinos_%d", 1); hDGenNeutrinos[1] = new TH1D(sb,sb,100,0.0,200.0);
338			sprintf(sb,"hDGenNeutrinos_%d", 2); hDGenNeutrinos[2] = new TH1D(sb,sb,100,-5.0,5.0);
339			sprintf(sb,"hDGenNeutrinos_%d", 3); hDGenNeutrinos[3] = new TH1D(sb,sb,90,0.0,180.0);
340	loizides	1.5	for(Int_t i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);
341	loizides	1.1
342			// Quarks
343			sprintf(sb,"hDGenQuarks_%d", 0); hDGenQuarks[0] = new TH1D(sb,sb,10,-0.5,9.5);
344			sprintf(sb,"hDGenQuarks_%d", 1); hDGenQuarks[1] = new TH1D(sb,sb,200,0.0,400.);
345			sprintf(sb,"hDGenQuarks_%d", 2); hDGenQuarks[2] = new TH1D(sb,sb,2000,0.0,2000.);
346			sprintf(sb,"hDGenQuarks_%d", 3); hDGenQuarks[3] = new TH1D(sb,sb,200,0.0,400.);
347			sprintf(sb,"hDGenQuarks_%d", 4); hDGenQuarks[4] = new TH1D(sb,sb,2000,0.0,2000.);
348	ceballos	1.2	sprintf(sb,"hDGenQuarks_%d", 5); hDGenQuarks[5] = new TH1D(sb,sb,200,0.0,400.);
349			sprintf(sb,"hDGenQuarks_%d", 6); hDGenQuarks[6] = new TH1D(sb,sb,200,-10.0,10.);
350	loizides	1.5	sprintf(sb,"hDGenQuarks_%d", 7); hDGenQuarks[7] = new TH1D(sb,sb,200,-TMath::Pi(),
351			TMath::Pi());
352	ceballos	1.2	sprintf(sb,"hDGenQuarks_%d", 8); hDGenQuarks[8] = new TH1D(sb,sb,200,0.0,400.);
353			sprintf(sb,"hDGenQuarks_%d", 9); hDGenQuarks[9] = new TH1D(sb,sb,200,-10.0,10.);
354	loizides	1.5	sprintf(sb,"hDGenQuarks_%d",10); hDGenQuarks[10] = new TH1D(sb,sb,200,-TMath::Pi(),
355			TMath::Pi());
356	ceballos	1.2	sprintf(sb,"hDGenQuarks_%d",11); hDGenQuarks[11] = new TH1D(sb,sb,200,0.0,400.);
357			sprintf(sb,"hDGenQuarks_%d",12); hDGenQuarks[12] = new TH1D(sb,sb,200,-10.0,10.);
358	loizides	1.5	sprintf(sb,"hDGenQuarks_%d",13); hDGenQuarks[13] = new TH1D(sb,sb,200,-TMath::Pi(),
359			TMath::Pi());
360	ceballos	1.2	sprintf(sb,"hDGenQuarks_%d",14); hDGenQuarks[14] = new TH1D(sb,sb,200,0.0,400.);
361			sprintf(sb,"hDGenQuarks_%d",15); hDGenQuarks[15] = new TH1D(sb,sb,200,-10.0,10.);
362	loizides	1.5	sprintf(sb,"hDGenQuarks_%d",16); hDGenQuarks[16] = new TH1D(sb,sb,200,-TMath::Pi(),
363			TMath::Pi());
364			for(Int_t i=0; i<17; i++) AddOutput(hDGenQuarks[i]);
365	loizides	1.1
366			// qqH
367			sprintf(sb,"hDGenWBF_%d", 0); hDGenWBF[0] = new TH1D(sb,sb,90,0.0,180.);
368			sprintf(sb,"hDGenWBF_%d", 1); hDGenWBF[1] = new TH1D(sb,sb,100,0.0,10.);
369			sprintf(sb,"hDGenWBF_%d", 2); hDGenWBF[2] = new TH1D(sb,sb,200,0.0,400.);
370			sprintf(sb,"hDGenWBF_%d", 3); hDGenWBF[3] = new TH1D(sb,sb,200,0.0,400.);
371			sprintf(sb,"hDGenWBF_%d", 4); hDGenWBF[4] = new TH1D(sb,sb,200,0.0,4000.);
372	loizides	1.5	for(Int_t i=0; i<5; i++) AddOutput(hDGenWBF[i]);
373	loizides	1.1
374			// Bosons
375			sprintf(sb,"hDGenBosons_%d", 0); hDGenBosons[0] = new TH1D(sb,sb,10,-0.5,9.5);
376			sprintf(sb,"hDGenBosons_%d", 1); hDGenBosons[1] = new TH1D(sb,sb,200,0.0,400.0);
377			sprintf(sb,"hDGenBosons_%d", 2); hDGenBosons[2] = new TH1D(sb,sb,100,-5.0,5.0);
378			sprintf(sb,"hDGenBosons_%d", 3); hDGenBosons[3] = new TH1D(sb,sb,2000,0.0,2000.0);
379	loizides	1.5	for(Int_t i=0; i<4; i++) AddOutput(hDGenBosons[i]);
380	loizides	1.1	}
381			}