JetCombinatorics.h

Go to the documentation of this file.

#ifndef JetCombinatorics_h
#define JetCombinatorics_h

#include "TLorentzVector.h"
#include "TString.h"
#include "TH1F.h"
#include "TFile.h"
#include "TMath.h"
#include <map>
#include <vector>
#include <iostream>

class Combo {


  public:

    Combo() {
        
      MW = 84.2;//79.8;
      Mtop_h = 180.7;//175.;
      Mtop_l = 174.9;
      sigmaHadW = 10.5;//2.*7.6;
      sigmaHadt = 19.2;//2.*12.5;
      sigmaLept = 24.2;//2.*15.6;

      SumEt_ = 0.;
      usebtag_ = false;
      useMtop_ = true;

      useFlv_ = false;
      Wp_flv_ = Wq_flv_ = Hadb_flv_ = Lepb_flv_ = 1.;
    }
    ~Combo(){};

    void SetWp(const TLorentzVector& Wp) { Wp_ = Wp; }
    void SetWq(const TLorentzVector& Wq) { Wq_ = Wq; }
    void SetHadb(const TLorentzVector& Hadb) { Hadb_ = Hadb; }
    void SetLepW(const TLorentzVector& LepW) { LepW_ = LepW; }
    void SetLepb(const TLorentzVector& Lepb) { Lepb_ = Lepb; }
    // flavor corrections
    void ApplyFlavorCorrections(bool option=true){ useFlv_ = option;}
    void SetFlvCorrWp( double corr ) { Wp_flv_ = corr; }
    void SetFlvCorrWq( double corr ) { Wq_flv_ = corr; }
    void SetFlvCorrHadb( double corr ) { Hadb_flv_ = corr; }
    void SetFlvCorrLepb( double corr ) { Lepb_flv_ = corr; }
    // b tagging
    void SetWp_disc(double disc) { Wp_disc_ = disc;}
    void SetWq_disc(double disc) { Wq_disc_= disc;}
    void SetHadb_disc(double disc) { Hadb_disc_= disc;}
    void SetLepb_disc(double disc) { Lepb_disc_= disc;}
    void SetbDiscPdf(const TString& filename) { 
      pdffile_ = TFile::Open(filename);
      hdisc_b_ = (TH1F*) gDirectory->Get("hdiscNorm_b");
      hdisc_cl_ = (TH1F*) gDirectory->Get("hdiscNorm_cl");
    }
    void SetSigmas(int type=0) {

      // type == 0 take defaults
      if (type==1) {
        // JES +10%
        MW = 87.2;
        Mtop_h = 193.2;
        Mtop_l = 179.0;
        sigmaHadW = 13.0;
        sigmaHadt = 22.8;
        sigmaLept = 26.3;
      }
      if (type==-1) {
            // JES -10%
            MW = 81.6;
            Mtop_h = 169.3;
            Mtop_l = 171.4;
            sigmaHadW =8.9;
            sigmaHadt =17.9;
            sigmaLept =22.6;
      }

    }
    void Usebtagging(bool option = true) { usebtag_ = option;}
    void SetMinMassLepW( double mass ) { minMassLepW_ = mass; }
    void SetMaxMassLepW( double mass ) { maxMassLepW_ = mass; }
    void SetMinMassHadW( double mass ) { minMassHadW_ = mass; }
    void SetMaxMassHadW( double mass ) { maxMassHadW_ = mass; }
    void SetMinMassLepTop( double mass ) { minMassLepTop_ = mass; }
    void SetMaxMassLepTop( double mass ) { maxMassLepTop_ = mass; }
    void UseMtopConstraint(bool option=true) { useMtop_ = option; }
        
    void analyze() {

        if ( useFlv_ ) {
            Wp_ = Wp_flv_ * Wp_;
            Wq_ = Wq_flv_ * Wq_;
            Hadb_ = Hadb_flv_ * Hadb_;
            Lepb_ = Lepb_flv_ * Lepb_;
        }

        HadW_ = Wp_ + Wq_;
        HadTop_ = HadW_ + Hadb_;
        LepTop_ = LepW_ + Lepb_;
        TopPair_ = HadTop_ + LepTop_;

        //double sigmaHadW = 10.5;//2.*7.6;
        //double sigmaHadt = 19.2;//2.*12.5;
        //double sigmaLept = 24.2;//2.*15.6;
        
        double chiHadW = (HadW_.M() - MW)/sigmaHadW;
        double chiHadt = (HadTop_.M() - Mtop_h)/sigmaHadt;
        double chiLept = (LepTop_.M() - Mtop_l)/sigmaLept;

        if ( useMtop_ ) {
            chi2_ = chiHadW*chiHadW + chiHadt*chiHadt + chiLept*chiLept;
            Ndof_ = 3;
        } else {
            chi2_ = chiHadW*chiHadW + (HadTop_.M() - LepTop_.M())*(HadTop_.M() - LepTop_.M())/(sigmaHadt*sigmaHadt+sigmaLept*sigmaLept);
            Ndof_ = 2;
        }
        
        SumEt_ = HadTop_.Pt();

        if ( usebtag_ ) {

            double gauss_norm = (2.)*TMath::Log(sigmaHadW*TMath::Sqrt(2*TMath::Pi())) +
                (2.)*TMath::Log(sigmaHadt*TMath::Sqrt(2*TMath::Pi())) + (2.)*TMath::Log(sigmaLept*TMath::Sqrt(2*TMath::Pi()));

            double LR_Wp; double LR_Wq;
            double LR_Hadb; double LR_Lepb;

            double LR_den = 0;
            LR_den = ( getPdfValue("cl", Wp_disc_) + getPdfValue("b", Wp_disc_));
            if (LR_den == 0 ) LR_Wp = 1e-5;
            else LR_Wp = getPdfValue( "cl", Wp_disc_ )/ LR_den;

            LR_den = ( getPdfValue("cl", Wq_disc_) + getPdfValue("b", Wq_disc_));
            if (LR_den == 0 ) LR_Wq = 1e-5;
            else LR_Wq = getPdfValue( "cl", Wq_disc_ )/ LR_den;

            LR_den = ( getPdfValue("cl", Hadb_disc_) + getPdfValue("b", Hadb_disc_));
            if (LR_den == 0 ) LR_Hadb = 1e-5;
            else LR_Hadb = getPdfValue( "b", Hadb_disc_ )/ LR_den;

            LR_den = ( getPdfValue("cl", Lepb_disc_) + getPdfValue("b", Lepb_disc_));
            if (LR_den == 0 ) LR_Lepb = 1e-5;
            else LR_Lepb = getPdfValue( "b", Lepb_disc_ )/ LR_den;

            double btag_norm = (-0.25-TMath::Log(4)/2);
            double btag_N2LL = btag_norm*4.*( LR_Wp * TMath::Log(LR_Wp/4) + LR_Wq*TMath::Log(LR_Wq/4) + LR_Hadb*TMath::Log(LR_Hadb/4) + LR_Lepb*TMath::Log(LR_Lepb/4) );
          
            chi2_ += btag_N2LL + gauss_norm;
            Ndof_ += 3;
            pdffile_->Close();
        }
    }

    TLorentzVector GetWp() { return Wp_; }
    TLorentzVector GetWq() { return Wq_; }
    TLorentzVector GetHadW() { return HadW_; }
    TLorentzVector GetLepW() { return LepW_; }
    TLorentzVector GetHadb() { return Hadb_; }
    TLorentzVector GetLepb() { return Lepb_; }
    TLorentzVector GetHadTop() { return HadTop_; }
    TLorentzVector GetLepTop() { return LepTop_; }
    TLorentzVector GetTopPair() { return TopPair_; }
    double GetChi2() const { return chi2_; }
    double GetNdof() { return Ndof_; }
    double GetSumEt() const { return SumEt_; }
    int GetIdHadb() { return IdHadb_;}
    int GetIdWp() { return IdWp_; }
    int GetIdWq() { return IdWq_; }
    int GetIdLepb() { return IdLepb_;}
    void SetIdHadb(int id) { IdHadb_ = id;}
    void SetIdWp(int id) { IdWp_ = id; }
    void SetIdWq(int id) { IdWq_ = id; }
    void SetIdLepb(int id) { IdLepb_ = id;}
    void Print() {
      std::cout << " jet Wp  : px = " << Wp_.Px() << " py = " <<  Wp_.Py() << " pz = " << Wp_.Pz() << " e = " << Wp_.E() << std::endl;
      std::cout << " jet Wq  : px = " << Wq_.Px() << " py = " <<  Wq_.Py() << " pz = " << Wq_.Pz() << " e = "<< Wq_.E() << std::endl;
      std::cout << " jet Hadb: px = " << Hadb_.Px() << " py = " <<  Hadb_.Py() <<" pz = " << Hadb_.Pz() <<" e = "<< Hadb_.E() << std::endl;
      std::cout << " jet Lepb: px = " << Lepb_.Px() << " py = " <<  Lepb_.Py() <<" pz = " << Lepb_.Pz() <<" e = "<< Lepb_.E() << std::endl;
      std::cout << " chi-squared = " << chi2_ << " sumEt = " << SumEt_ << std::endl;
    }
    double getPdfValue(std::string flavor, double disc) {
      double pdf= 0;
      TH1F *hpdf;
      if ( flavor == "b" ) hpdf = hdisc_b_;
      else hpdf = hdisc_cl_;
      int bin = hpdf->GetXaxis()->FindBin( disc );
      pdf = hpdf->GetBinContent( bin );
      if ( disc < -10 || disc >50 ) return 0;
      //if ( pdf == 0 ) return 1.e-7;
      return pdf;
    }
    
  private:
    
    TLorentzVector Wp_;
    TLorentzVector Wq_;
    TLorentzVector HadW_;
    TLorentzVector Hadb_;
    TLorentzVector HadTop_;
    TLorentzVector LepW_;
    TLorentzVector Lepb_;   
    TLorentzVector LepTop_;
    TLorentzVector TopPair_;
    
    bool usebtag_;
    bool useMtop_;
    double Wp_disc_;
    double Wq_disc_;
    double Hadb_disc_;
    double Lepb_disc_;
    TFile *pdffile_;
    TH1F *hdisc_b_;
    TH1F *hdisc_cl_;

    double Wp_flv_, Wq_flv_, Hadb_flv_, Lepb_flv_;
    bool useFlv_;
    double chi2_;
    double Ndof_;
    double SumEt_;
    double minMassLepW_;
    double maxMassLepW_;
    double minMassHadW_;
    double maxMassHadW_;
    
    double minMassLepTop_;
    double maxMassLepTop_;

    double MW;
    double Mtop_h;
    double Mtop_l;
    double sigmaHadW;
    double sigmaHadt;
    double sigmaLept;


    int IdHadb_;
    int IdWp_;
    int IdWq_;
    int IdLepb_;
    
};

struct minChi2
{
  bool operator()(const Combo& s1, const Combo& s2) const
  {
    return s1.GetChi2() <= s2.GetChi2();
  }
};

struct maxSumEt
{
  bool operator()(const Combo& s1,const Combo& s2) const
  {
    return s1.GetSumEt() >= s2.GetSumEt();
  }
};


class JetCombinatorics {

  public:

    JetCombinatorics();
    ~JetCombinatorics();

    void Verbose() {
      verbosef = true;
    }

    std::map< int, std::string > Combinatorics(int k, int max = 6);
    std::map< int, std::string > NestedCombinatorics();

    void FourJetsCombinations(const std::vector<TLorentzVector>& jets, const std::vector<double>& bdiscriminators );
    void SetFlavorCorrections(const std::vector<double >& vector ) { flavorCorrections_ = vector; }
    void SetMaxNJets(int n) { maxNJets_ = n; }
    Combo GetCombination(int n=0);
    Combo GetCombinationSumEt(int n=0);
    int GetNumberOfCombos() { return ( (int)allCombos_.size() ); } 
    //void SetCandidate( std::vector< TLorentzVector > JetCandidates );

    void SetSigmas(int type = 0) {
      SigmasTypef = type;
    }
    void SetLeptonicW( const TLorentzVector& LepW ) { theLepW_ = LepW; }

    void SetMinMassLepW( double mass ) { minMassLepW_ = mass; }
    void SetMaxMassLepW( double mass ) { maxMassLepW_ = mass; }
    void SetMinMassHadW( double mass ) { minMassHadW_ = mass; }
    void SetMaxMassHadW( double mass ) { maxMassHadW_ = mass; }
    void SetMinMassLepTop( double mass ) { minMassLepTop_ = mass; }
    void SetMaxMassLepTop( double mass ) { maxMassLepTop_ = mass; }

    void UsebTagging( bool option = true ) { UsebTagging_ = option; }
    void ApplyFlavorCorrection( bool option = true ) { UseFlv_ = option; }
    void UseMtopConstraint( bool option = true) { UseMtop_ = option; }
    void SetbTagPdf(const TString& name ) { bTagPdffilename_ = name; }
    void Clear();

    std::vector< TLorentzVector > TwoCombos();
    std::vector< TLorentzVector > ThreeCombos();

    void RemoveDuplicates( bool option) { removeDuplicates_ = option; }

    std::vector< TLorentzVector > GetComposites();
    void AnalyzeCombos();


  private:

    //int kcombos_;
    //int maxcombos_;
    int SigmasTypef;
    bool verbosef;
    std::map< int, std::string > Template4jCombos_;
    std::map< int, std::string > Template5jCombos_;
    std::map< int, std::string > Template6jCombos_;
    std::map< int, std::string > Template7jCombos_;

    std::vector< double > flavorCorrections_;
    int maxNJets_;
    bool UsebTagging_;
    bool UseMtop_;
    TString bTagPdffilename_;
    bool UseFlv_;
    
    TLorentzVector theLepW_;

    double minMassLepW_;
    double maxMassLepW_;
    double minMassHadW_;
    double maxMassHadW_;
    double minMassLepTop_;
    double maxMassLepTop_;
    
    std::map< Combo, int, minChi2 > allCombos_;
    std::map< Combo, int, maxSumEt > allCombosSumEt_;

    Double_t minPhi_;
    double chi2_;
    int ndf_;
    bool removeDuplicates_;
    
    std::vector< TLorentzVector > cand1_;
    std::vector< TLorentzVector > cand2_;
    std::vector< TLorentzVector > cand3_;

    //int nLists_;
    
    //std::vector< TLorentzVector > composites_;
    
};

#endif