Mods/src/GeneratorMod.cc

// $Id: GeneratorMod.cc,v 1.13 2008/12/02 09:35:54 loizides Exp $

#include "MitPhysics/Mods/interface/GeneratorMod.h"
#include "MitCommon/MathTools/interface/MathUtils.h"
#include "MitPhysics/Init/interface/ModNames.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(ModNames::gkMCLeptonsName),
  fMCAllLeptonsName(ModNames::gkMCAllLeptonsName),
  fMCTausName(ModNames::gkMCTausName),
  fMCNeutrinosName(ModNames::gkMCNeutrinosName),
  fMCQuarksName(ModNames::gkMCQuarksName),
  fMCqqHsName(ModNames::gkMCqqHsName),
  fMCBosonsName(ModNames::gkMCBosonsName),
  fMCPhotonsName(ModNames::gkMCPhotonsName),
  fParticles(0),
  fPtLeptonMin(0.0),
  fEtaLeptonMax(5.0),
  fPtPhotonMin(0.0),
  fEtaPhotonMax(5.0)
{
  // Constructor.
}

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

  // these arrays will be filled in the loop of particles
  MCParticleOArr *GenLeptons    = new MCParticleOArr;
  GenLeptons->SetName(fMCLeptonsName);
  MCParticleOArr *GenAllLeptons = new MCParticleOArr;
  GenAllLeptons->SetName(fMCAllLeptonsName);
  MCParticleOArr *GenTaus       = new MCParticleOArr; 
  GenTaus->SetName(fMCTausName);
  GenTaus->SetOwner(kTRUE);
  MCParticleOArr *GenNeutrinos  = new MCParticleOArr;
  GenNeutrinos->SetName(fMCNeutrinosName);
  MCParticleOArr *GenQuarks     = new MCParticleOArr;
  GenQuarks->SetName(fMCQuarksName);
  MCParticleOArr *GenqqHs       = new MCParticleOArr;
  GenqqHs->SetName(fMCqqHsName);
  MCParticleOArr *GenBosons     = new MCParticleOArr;
  GenBosons->SetName(fMCBosonsName);
  MCParticleOArr *GenPhotons    = new MCParticleOArr;
  GenPhotons->SetName(fMCPhotonsName);

  // load MCParticle branch
  LoadBranch(fMCPartName);

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

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

    // muons/electrons from W/Z decays
    if ((p->Is(MCParticle::kEl) || p->Is(MCParticle::kMu)) && p->Status() == 1) {
      if (p->Pt() > fPtLeptonMin && p->AbsEta() < fEtaLeptonMax) {
        GenAllLeptons->Add(p);
      }
      Bool_t isGoodLepton = kFALSE;
      const MCParticle *pm = p;
      while (pm->HasMother() && isGoodLepton == kFALSE) {
        if (pm->Mother()->Is(MCParticle::kZ) || pm->Mother()->Is(MCParticle::kW)) {
          GenLeptons->Add(p);
          isGoodLepton = kTRUE;
          break;
        } else if (pm->Mother()->Is(MCParticle::kPi0) || pm->Mother()->Is(MCParticle::kEta)) {
          // this is fake, but it is a trick to get rid of these cases and abort the loop
          isGoodLepton = kTRUE;
          break;
        } 
        pm = pm->Mother();
      }
    }

    // hadronic taus
    else if (p->Is(MCParticle::kTau) && p->Status() == 2) {
      if (!p->HasDaughter(MCParticle::kEl) && !p->HasDaughter(MCParticle::kMu)) {
        const MCParticle *tv = p->FindDaughter(MCParticle::kTauNu);
        if (tv) {
          MCParticle *pm_f = new MCParticle(*p);
          pm_f->SetMom(p->Px()-tv->Px(), p->Py()-tv->Py(),
                       p->Pz()-tv->Pz(), p->E()-tv->E());
          GenTaus->AddOwned(pm_f);
        } else {
          SendError(kWarning, "Process", "Could not find tau neutrino!");
        }
      }
    }

    // neutrinos
    else if (p->Status() == 1 && p->IsNeutrino()) {
      GenNeutrinos->Add(p);
    }

    // quarks from W/Z decays or top particles
    else if (p->IsQuark() && p->HasMother()) {
      if (p->Mother()->Is(MCParticle::kZ) || p->Mother()->Is(MCParticle::kW) ||
          p->Is(MCParticle::kTop)         || p->Mother()->Is(MCParticle::kTop)) {
        GenQuarks->Add(p);
      }
    }

    // qqH, information about the forward jets
    else if (isqqH == kFALSE && p->Is(MCParticle::kH)) {
      isqqH = kTRUE;
      MCParticle *pq1 = fParticles->At(i-1);
      MCParticle *pq2 = fParticles->At(i-2);
      if (!pq1 || !pq2) {
          SendError(kWarning, "Process", "Could not find quark pair!");
      } else if (pq1->IsQuark()   && pq2->IsQuark()   && 
                 pq1->HasMother() && pq2->HasMother() &&
                 pq1->Mother()->PdgId() == p->Mother()->PdgId() &&
                 pq2->Mother()->PdgId() == p->Mother()->PdgId()) {
        GenqqHs->Add(pq1);
        GenqqHs->Add(pq2);
      }
      if (p->Status() == 3)  
        GenBosons->Add(p); // take higgs boson in  account here rather in next else if 
    }

    // information about bosons: W, Z, h, Z', W', H0, A0, H+
    else if (p->Status() == 3 &&
             (p->Is(MCParticle::kZ)  || p->Is(MCParticle::kW)   || p->Is(MCParticle::kH) ||
              p->Is(MCParticle::kZp) || p->Is(MCParticle::kZpp) ||
              p->Is(MCParticle::kH0) || p->Is(MCParticle::kA0)  || p->Is(MCParticle::kHp))) {
      GenBosons->Add(p);
    }

    // photons
    else if (p->Status() == 1 && p->Is(MCParticle::kGamma) &&
             p->Pt() > fPtPhotonMin && p->AbsEta() < fEtaPhotonMax) {
      GenPhotons->Add(p);
    }

  }

  // add objects to this event for other modules to use
  AddObjThisEvt(GenLeptons);  
  AddObjThisEvt(GenAllLeptons);
  AddObjThisEvt(GenTaus);
  AddObjThisEvt(GenNeutrinos);
  AddObjThisEvt(GenQuarks);
  AddObjThisEvt(GenqqHs);
  AddObjThisEvt(GenBosons);
  AddObjThisEvt(GenPhotons);

  // 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());
    } 

    // photons
    hDGenPhotons[0]->Fill(GenPhotons->GetEntries());
    for(UInt_t i=0; i<GenPhotons->GetEntries(); i++) {
      hDGenPhotons[1]->Fill(GenPhotons->At(i)->Pt());
      hDGenPhotons[2]->Fill(GenPhotons->At(i)->Eta());
    } 
  }
}

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

  ReqBranch(fMCPartName, fParticles);

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

    // leptons from W
    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]);

    // all leptons
    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,100,-5.0,5.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]);

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