bbbarAnalysis.cc

Go to the documentation of this file.

// STL 
#include <iostream>
#include <sstream>
#include <stdio.h>      // ???? TODO

#include "HepMC/GenEvent.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/SimpleVector.h"
#include "HepPDT/ParticleData.hh"
#include "CLHEP/Vector/LorentzVector.h"

// ROOT 
#include "TH1.h"
#include "TMath.h"
#include "TLorentzVector.h"
#include "TFile.h"

// fastjet 
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/SISConePlugin.hh"

// Analysis class
#include "../include/bbbarAnalysis.h"

using namespace std;

// ---------------------------------------------------------------------- 
bbbarAnalysis::bbbarAnalysis( char *fileName ) : 
  OutputFileName( fileName )  
{
  clearCounters();
  //initPlots();
}

// ---------------------------------------------------------------------- 
bbbarAnalysis::bbbarAnalysis() : 
  OutputFileName( "plots_lhcb.root" ) 
{
  clearCounters();
  //initPlots();
}

// ---------------------------------------------------------------------- 
bbbarAnalysis::~bbbarAnalysis() 
{
  // write output file with histograms
  // hOutputFile->Write();
  // hOutputFile->Close();
}

// ---------------------------------------------------------------------- 
void bbbarAnalysis::initPlots() 
{
  // output file for histograms
  //hOutputFile = new TFile(OutputFileName, "RECREATE");
  // ROOT histograms definition
  hPtmb = new TH1D("hPtmb", "Pt of stable charged particles in MB events",
           100, 0., 20.);
  hPtbbbar = new TH1D("hPtbbbar", "Pt of stable charged particles in bbbar events",
              100, 0., 20.);
  hPtbhadr = new TH1D("hPtbhadr", "Pt of B hadrons in bbbar events",
              100, 0., 20.);
  hEtabhadr = new TH1D("hEtabhadr", "Eta of B hadrons in bbbar events",
               100, -10., 10.);
}

// ---------------------------------------------------------------------- 
void bbbarAnalysis::clearCounters() 
{
  // reset the number of 'b events'
  Bcount = 0;
  BcountTotal = 0;
  // reset the number of Bd/Bu/Bs/b-baryons
  Bdcount = 0;
  Bucount = 0;
  Bscount = 0;
  Barions = 0;
  // reset the numbers of stable charged particles in MB and b-bbar events
  ChargedPartMB = 0;
  ChargedPartBBbar = 0;
}

// ---------------------------------------------------------------------- 
int bbbarAnalysis::Init( double tr_max_eta, double tr_min_pt )
{
  // specify default eta cut
  m_max_eta=tr_max_eta;

  // specify default pt cut
  m_min_pt=tr_min_pt;

  // default Output file name
  m_outputFileName  = "bbarAnalysis.root";
  m_outputRootDir   = "bbar";

  //m_h_njets   = initHist("njets", "N_jets", "Number of Jets", 15, 0., 15.);

  // ROOT histograms definition
  hPtmb = initHist("hPtmb", "Ptmb", "Pt of stable charged particles in MB events",
           100, 0., 20.);
  hPtbbbar = initHist("hPtbbbar", "Ptbbbar", "Pt of stable charged particles in bbbar events",
              100, 0., 20.);
  hPtbhadr = initHist("hPtbhadr", "Ptbhadr", "Pt of B hadrons in bbbar events",
              100, 0., 20.);
  hEtabhadr = initHist("hEtabhadr", "Etabhadr", "Eta of B hadrons in bbbar events",
               100, -10., 10.);

  //initPlots();
  return true;
}

// ---------------------------------------------------------------------- 
int bbbarAnalysis::Process( HepMC::GenEvent* hepmcevt ) 
{
  getBbbar( hepmcevt );
  return true;
}

// ---------------------------------------------------------------------- 
void bbbarAnalysis::getBbbar(HepMC::GenEvent* evt) 
{
//  std::cout << "getBbbar, in the beginning" << std::endl;

  HepMC::GenEvent::particle_iterator part;
  bool bbbar = false;
  bool bhadr = false;
//  IsFinalState isStable;

  int BcountInEvent = 0;

//  std::cout << "getBbbar, before the loop" << std::endl;

  // Check if b-bbar pair is present and count it
  for (part = evt->particles_begin(); part != evt->particles_end(); part++)
    {
//      std::cout << "getBbbar, first the loop" << std::endl;
      int pdgID = (*part)->pdg_id();
      int apdgID = abs(pdgID);

      // check if b-quark is present
      if (apdgID == 5) {
    bbbar = true;
 //   std::cout << " pdgID = " << pdgID << "  status = " << (*part)->status() << std::endl;
    if ((*part)->status() == 2) BcountInEvent++;
      }

      // measure fractions of Bd/Bu/Bs in b-bbar events
      if (apdgID == 511) { Bdcount++; bhadr = true; }
      if (apdgID == 521) { Bucount++; bhadr = true; }
      if (apdgID == 531) { Bscount++; bhadr = true; }

      // measure fractions of b-baryons in b-bbar events
      if (apdgID == 5122 || apdgID == 5232 || apdgID == 5132) {
    Barions++; bhadr = true; }

      double Pt = (*part)->momentum().perp();
      double Eta = (*part)->momentum().pseudoRapidity();

      if (bhadr) {             // B hadron
    hPtbhadr->Fill(Pt);
    hEtabhadr->Fill(Eta);
      }


      if ( !(*part)->end_vertex() && (*part)->status()==1 )
//      if (isStable(*part))     // stable particle
    {
//      if(pdgID > 3000 && pdgID < 4000) std::cout << "charm-beauty stable particle" << std::endl;

      if ( chargedParticle(pdgID) )    // charged particle
        {
          ChargedPartMB++;
          hPtmb->Fill(Pt); // MB stable charged particle

          if (bbbar) {     // b-bbar stable charged particle
        ChargedPartBBbar++;
        hPtbbbar->Fill(Pt);
          }
        }

    }
    }
  if (bbbar) { Bcount++; BcountTotal += BcountInEvent/2; }
//  std::cout << "getBbbar, at the end" << std::endl;
}

---