JetCombinatorics.cc

Go to the documentation of this file.

#include "JetCombinatorics.h"
#include "TMath.h"

#include <iostream>

std::string itoa(int i) {
  char temp[20];
  sprintf(temp,"%d",i);
  return((std::string)temp);
}



//_______________________________________________________________
JetCombinatorics::JetCombinatorics() {

    this->Clear();
        
    minMassLepW_ = -999999.;
    maxMassLepW_ = 999999.;
    minMassHadW_ = -999999.;
    maxMassHadW_ = 999999.;
    minMassLepTop_ = -999999.;
    maxMassLepTop_ = 999999.;
    
    minPhi_ = -1.;
    removeDuplicates_ = true;
    maxNJets_ = 9999;
    verbosef = false;
    UsebTagging_ = false;
    UseMtop_ = true;
    SigmasTypef = 0;
    UseFlv_ = false;
    
    Template4jCombos_ = NestedCombinatorics(); // 12 combinations
    Template5jCombos_ = Combinatorics(4,5); // 5 combinations of 4 combos
    Template6jCombos_ = Combinatorics(4,6); // 15 combinations of 4 combos
    Template7jCombos_ = Combinatorics(4,7); // xx combinations of 4 combos
    
}

//_______________________________________________________________
JetCombinatorics::~JetCombinatorics() {
  this->Clear();
}

//_______________________________________________________________
void JetCombinatorics::Clear() {

    allCombos_.clear();
    allCombosSumEt_.clear();
    Template4jCombos_.clear();
    Template5jCombos_.clear();
    Template6jCombos_.clear();
    Template7jCombos_.clear();
    cand1_.clear();
    
}
    

//_______________________________________________________________
std::map< int, std::string > JetCombinatorics::Combinatorics(int n, int max) {

    // find a combinatorics template
    // This is a simple stupid function to make algebratic combinatorics

    int kcombos = n;
    int maxcombos = max;

    std::string list;

    for ( int m=0; m<maxcombos; m++) { list = list + (itoa(m));}

    std::string seed;
    for ( int m=0; m<kcombos; m++) { seed = seed + (itoa(m));}

    
    std::map< int, std::string > aTemplateCombos;
    aTemplateCombos.clear();
    
    aTemplateCombos[0] = seed;

    int i = 0;
    int totalmatches = seed.size();
    int totalIte = list.size();

    for ( int ite = 0; ite < ((int)totalIte); ite++) {

        //cout << "iteration " << ite << endl;
        //i = 0;
        //for ( Itevec = seed.begin(); Itevec != seed.end(); ++Itevec) {
        for ( i=0; i< (int) totalmatches; i++) {

            std::string newseed = aTemplateCombos[ite];
            std::string newseed2;
            /*
            cout << "newseed size= " << newseed.size() << " : ";
            for (std::vector< std::string>::iterator iite = newseed.begin();
                 iite != newseed.end(); ++iite) {

                cout << *iite << " ";
            }
            cout << endl;
            */
            for ( int itemp=0; itemp<(int)newseed.size(); itemp++) {
                if (itemp!=i) newseed2 = newseed2 + (newseed[itemp]);
            }
            /*
            cout << "newseed2: ";
            for (std::vector< std::string>::iterator iite = newseed2.begin();
                 iite != newseed2.end(); ++iite) {

                cout << *iite << " ";
            }
            cout << endl;
            */
            for ( int j=0; j<(int) list.size(); j++) {
                //cout << " j = " << j << endl;
                bool Isnewelement = true;
                std::string newelement = "0";
                //bool Isnewcombo = true;
                for (int k=0; k< (int)newseed2.size(); k++) {
                    if ( list[j] == newseed2[k] ) Isnewelement = false;
                }
                if (Isnewelement) {

                    newelement = list[j];
                    //cout << "new element: " << newelement << endl;

                    std::string candseed = newseed2;
                    candseed = candseed + newelement;

                    bool IsnewCombo = true;
                    for (int ic=0; ic<(int)aTemplateCombos.size(); ++ic ) {

                        int nmatch = 0;
                        for ( int ij=0; ij<(int)(aTemplateCombos[ic]).size(); ij++) {

                            for (int ik=0; ik<(int)candseed.size(); ik++) {
                                if ( candseed[ik] == aTemplateCombos[ic][ij] ) nmatch++;
                            }
                        }
                        if (nmatch == (int)totalmatches)
                            IsnewCombo = false;

                    }
                    if (IsnewCombo) {
                        //cout << "new combo"<< " before combo size=" << aTemplateCombos.size() << endl;
                        aTemplateCombos[(int)aTemplateCombos.size()] = candseed;
                        //cout << " after size = " << aTemplateCombos.size() << endl;
                    }
                }

            }
        }
    }//close iterations

    // debug info
    
    //std::cout << " build templates for total combos = " << aTemplateCombos.size() << std::endl;
    //std::cout << " template combos: " << std::endl;
    //for (size_t ic=0; ic != aTemplateCombos.size(); ++ic) {

    //std::cout << aTemplateCombos[ic] << std::endl;
    //}
    
    return aTemplateCombos;
    
    
    
}


//______________________________________________________________
std::map< int, std::string > JetCombinatorics::NestedCombinatorics() {

    // build by hand 12 combinations for semileptonic top decays

    std::map< int, std::string > aTemplateCombos;
    aTemplateCombos.clear();
    
    aTemplateCombos[0] = "0123";
    aTemplateCombos[1] = "0132";
    aTemplateCombos[2] = "0213";
    aTemplateCombos[3] = "0231";
    aTemplateCombos[4] = "0312";
    aTemplateCombos[5] = "0321";
    aTemplateCombos[6] = "1203";
    aTemplateCombos[7] = "1230";
    aTemplateCombos[8] = "1302";
    aTemplateCombos[9] = "1320";
    aTemplateCombos[10] = "2301";
    aTemplateCombos[11] = "2310";
        
    return aTemplateCombos;

}

//______________________________________________________________
void JetCombinatorics::FourJetsCombinations(const std::vector<TLorentzVector>& jets, const std::vector<double>& bdiscriminators ) {


    int n = 0; // total number of combos
    std::map< Combo, int, minChi2 > allCombos;
    std::map< Combo, int, maxSumEt > allCombosSumEt;
    
    std::map< int, std::string > aTemplateCombos;
    aTemplateCombos.clear();

    if ( jets.size() == 4 ) aTemplateCombos[0] = std::string("0123");
    if ( jets.size() == 5 ) aTemplateCombos = Template5jCombos_;
    if ( jets.size() == 6 ) aTemplateCombos = Template6jCombos_;
    if ( jets.size() == 7 ) aTemplateCombos = Template7jCombos_;    

    // force to use only 4 jets
    if ( maxNJets_ == 4 ) aTemplateCombos[0] = std::string("0123");
    
    if (verbosef) std::cout << "[JetCombinatorics] size of vector of jets = " << jets.size() << std::endl;
    
    for (size_t ic=0; ic != aTemplateCombos.size(); ++ic) {

        if (verbosef) std::cout << "[JetCombinatorics] get 4 jets from the list, cluster # " << ic << "/"<< aTemplateCombos.size()-1 << std::endl;
        
        // get a template
        std::string aTemplate = aTemplateCombos[ic];

        if (verbosef) std::cout << "[JetCombinatorics] template of 4 jets = " << aTemplate << std::endl;
        
        // make a list of 4 jets
        std::vector< TLorentzVector > the4jets;
        std::vector< int > the4Ids;
        std::vector< double > thebdisc;
        std::vector< double > theFlvCorr;
        //the4jets[0] = jets[0];
        
        for (int ij=0; ij<4; ij++) {
            //std::cout << "ij= " << ij << std::endl;
            //std::cout << "atoi = " << atoi((aTemplate.substr(0,1)).c_str()) << std::endl;
            //std::cout << "jets[].Pt = " << jets[ij].Pt() << std::endl;
            int tmpi = atoi((aTemplate.substr(ij,1)).c_str());
            //std::cout << "tmpi= " << tmpi << std::endl;
            the4jets.push_back(jets[tmpi]);
            the4Ids.push_back(tmpi);
            if ( UsebTagging_ ) thebdisc.push_back( bdiscriminators[tmpi] );
            if ( UseFlv_ ) theFlvCorr.push_back( flavorCorrections_[tmpi] );
        }

        if (verbosef) std::cout<< "[JetCombinatorics] with these 4 jets, make 12 combinations: " <<std::endl;

        //std::cout << " the4jets[ij].size = " << the4jets.size() << std::endl;
            
        for (size_t itemplate=0; itemplate!= Template4jCombos_.size(); ++itemplate) {
            
            std::string a4template = Template4jCombos_[itemplate];

            if (verbosef) std::cout << "[JetCombinatorics] ==> combination: " << a4template << " is # " << itemplate << "/"<<  Template4jCombos_.size()-1 << std::endl;
            
            Combo acombo;
            
            acombo.SetWp( the4jets[atoi((a4template.substr(0,1)).c_str())] );
            acombo.SetWq( the4jets[atoi((a4template.substr(1,1)).c_str())] );
            acombo.SetHadb( the4jets[atoi((a4template.substr(2,1)).c_str())] );
            acombo.SetLepb( the4jets[atoi((a4template.substr(3,1)).c_str())] );
            acombo.SetLepW( theLepW_ );

            acombo.SetIdWp( the4Ids[atoi((a4template.substr(0,1)).c_str())] );
            acombo.SetIdWq( the4Ids[atoi((a4template.substr(1,1)).c_str())] );
            acombo.SetIdHadb( the4Ids[atoi((a4template.substr(2,1)).c_str())] );
            acombo.SetIdLepb( the4Ids[atoi((a4template.substr(3,1)).c_str())] );
            //std::cout << " acombo setup" << std::endl;

            if ( UseFlv_ ) {
                acombo.SetFlvCorrWp( theFlvCorr[atoi((a4template.substr(0,1)).c_str())] );
                acombo.SetFlvCorrWq( theFlvCorr[atoi((a4template.substr(1,1)).c_str())] );
                acombo.SetFlvCorrHadb( theFlvCorr[atoi((a4template.substr(2,1)).c_str())] );
                acombo.SetFlvCorrLepb( theFlvCorr[atoi((a4template.substr(3,1)).c_str())] );
                acombo.ApplyFlavorCorrections();
            }
            if ( UsebTagging_ ) {

                acombo.Usebtagging();
                acombo.SetbDiscPdf(bTagPdffilename_);
                acombo.SetWp_disc( thebdisc[atoi((a4template.substr(0,1)).c_str())] );
                acombo.SetWq_disc( thebdisc[atoi((a4template.substr(1,1)).c_str())] );
                acombo.SetHadb_disc( thebdisc[atoi((a4template.substr(2,1)).c_str())] );
                acombo.SetLepb_disc( thebdisc[atoi((a4template.substr(3,1)).c_str())] );
                
            }

            acombo.UseMtopConstraint(UseMtop_);
            // choose value of sigmas
            acombo.SetSigmas(SigmasTypef);

            acombo.analyze();

            if (verbosef) {

              std::cout << "[JetCombinatorics] ==> combination done:" << std::endl;
              acombo.Print();
            }

            // invariant mass cuts
            TLorentzVector aHadWP4 = acombo.GetHadW();
            TLorentzVector aLepWP4 = acombo.GetLepW();
            TLorentzVector aLepTopP4=acombo.GetLepTop();
            
            if ( ( aHadWP4.M() > minMassHadW_ && aHadWP4.M() < maxMassHadW_ ) &&
                 ( aLepWP4.M() > minMassLepW_ && aLepWP4.M() < maxMassLepW_ ) &&
                 ( aLepTopP4.M() > minMassLepTop_ && aLepTopP4.M() < maxMassLepTop_) ) {
            
                allCombos[acombo] = n;
                allCombosSumEt[acombo] = n;
            
                n++;
            }
        
        }
    }

    allCombos_ = allCombos;
    allCombosSumEt_ = allCombosSumEt;
       
}

Combo JetCombinatorics::GetCombination(int n) {

    int j = 0;
    Combo a;
    for ( std::map<Combo,int,minChi2>::const_iterator ite=allCombos_.begin();
          ite!=allCombos_.end(); ++ite) {
        
        if (j == n) a = ite->first;
        j++;
    }

    return a;

    
}

Combo JetCombinatorics::GetCombinationSumEt(int n) {

    int j = 0;
    Combo a;
    for ( std::map<Combo,int,maxSumEt>::const_iterator ite=allCombosSumEt_.begin();
          ite!=allCombosSumEt_.end(); ++ite) {
        
        if (j == n) a = ite->first;
        j++;
    }

    return a;

    
}