TtFullHadSignalSel.cc

Go to the documentation of this file.

#include "TopQuarkAnalysis/TopEventSelection/interface/TtFullHadSignalSel.h"
#include "PhysicsTools/CandUtils/interface/EventShapeVariables.h"
#include "PhysicsTools/CandUtils/interface/Thrust.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "TLorentzVector.h"
#include "TVector3.h"

TtFullHadSignalSel::TtFullHadSignalSel()
{
}

std::vector<math::XYZVector> makeVecForEventShape(std::vector<pat::Jet> jets, bool only6Jets = true, ROOT::Math::Boost boost = ROOT::Math::Boost(0.,0.,0.)) {
  std::vector<math::XYZVector> p;
  bool doBoost = (boost == ROOT::Math::Boost(0.,0.,0.)) ? false : true;
  for(std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet){
    math::XYZVector Vjet;
    if(doBoost){
      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > Ljet(jet->px(), jet->py(), jet->pz(), jet->energy());
      Vjet = math::XYZVector(boost(Ljet).Px(), boost(Ljet).Py(), boost(Ljet).Pz());
    }
    else
      Vjet = math::XYZVector(jet->px(), jet->py(), jet->pz());
    p.push_back(Vjet);
    if(only6Jets && jet-jets.begin()==5) break;
  }
  return p;
}

TtFullHadSignalSel::TtFullHadSignalSel(const std::vector<pat::Jet>& jets)
{

  H_      = -1.;
  Ht_     = -1.;
  Ht123_  = -1.;
  Ht3jet_ = -1.;
  Et56_   = -1.;
  sqrt_s_ = -1.;
  M3_     = -1.;
  
  TCHE_Bjets_   = 0.;
  TCHP_Bjets_   = 0.;
  SSVHE_Bjets_  = 0.;
  SSVHP_Bjets_  = 0.;
  CSV_Bjets_    = 0.;
  CSVMVA_Bjets_ = 0.;
  SM_Bjets_     = 0.;

  jets_etaetaMoment_ = 0.;
  jets_etaphiMoment_ = 0.;
  jets_phiphiMoment_ = 0.;

  jets_etaetaMomentLogEt_ = 0.;
  jets_etaphiMomentLogEt_ = 0.;
  jets_phiphiMomentLogEt_ = 0.;

  jets_etaetaMomentNoB_ = 0.;
  jets_etaphiMomentNoB_ = 0.;
  jets_phiphiMomentNoB_ = 0.;

  aplanarity_  = -1.;
  sphericity_  = -1.;
  circularity_ = -1.;
  isotropy_ = -1.;
  C_ = -1.;
  D_ = -1.;

  aplanarityAll_  = -1.;
  sphericityAll_  = -1.;
  circularityAll_ = -1.;
  isotropyAll_ = -1.;
  CAll_ = -1.;
  DAll_ = -1.;

  aplanarityAllCMS_  = -1.;
  sphericityAllCMS_  = -1.;
  circularityAllCMS_ = -1.;
  isotropyAllCMS_ = -1.;
  CAllCMS_ = -1.;
  DAllCMS_ = -1.;

  thrust_ = -1.;
  thrustCMS_ = -1.;

  TCHE_BJet_Discs_   = std::vector<double>(0);
  TCHP_BJet_Discs_   = std::vector<double>(0);
  SSVHE_BJet_Discs_  = std::vector<double>(0);
  SSVHP_BJet_Discs_  = std::vector<double>(0);
  CSV_BJet_Discs_    = std::vector<double>(0);
  CSVMVA_BJet_Discs_ = std::vector<double>(0);
  SM_BJet_Discs_     = std::vector<double>(0);

  pts_               = std::vector<double>(0);
  EtSin2Thetas_      = std::vector<double>(0);
  thetas_            = std::vector<double>(0);
  thetaStars_        = std::vector<double>(0);
  EtStars_           = std::vector<double>(0);

  EtSin2Theta3jet_  = 0.;
  theta3jet_        = 0.;
  thetaStar3jet_    = 0.;
  sinTheta3jet_     = 0.;
  sinThetaStar3jet_ = 0.;
  EtStar3jet_       = 0.;

  etaetaMoments_ = std::vector<double>(0);
  etaphiMoments_ = std::vector<double>(0);
  phiphiMoments_ = std::vector<double>(0);
 
  etaetaMomentsLogEt_ = std::vector<double>(0);
  etaphiMomentsLogEt_ = std::vector<double>(0);
  phiphiMomentsLogEt_ = std::vector<double>(0);
 
  etaetaMomentsMoment_ = std::vector<double>(0);
  etaphiMomentsMoment_ = std::vector<double>(0);
  phiphiMomentsMoment_ = std::vector<double>(0);
 
  etaetaMomentsMomentLogEt_ = std::vector<double>(0);
  etaphiMomentsMomentLogEt_ = std::vector<double>(0);
  phiphiMomentsMomentLogEt_ = std::vector<double>(0);
 
  etaetaMomentsNoB_ = std::vector<double>(0);
  etaphiMomentsNoB_ = std::vector<double>(0);
  phiphiMomentsNoB_ = std::vector<double>(0);
 
  dR_      = std::vector<double>(0);
  dRMass_  = std::vector<double>(0);
  dRAngle_ = std::vector<double>(0);

  dR3Jets_     = std::vector<double>(0);
  dR3JetsMass_ = std::vector<double>(0);

  massDiffMWCands_ = std::vector<double>(0);

  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > totalSystem(0.,0.,0.,0.);

  std::vector< std::pair< double, std::vector<unsigned short> > > dRs(0);
  std::vector< std::pair< double, std::vector<unsigned short> > > dRs3Jets(0);

  std::vector< std::pair< double, std::vector<unsigned short> > > M3s(0);

  std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > fourVectors(0);

  unsigned short nonBJets = 0;
  for(std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet){

    H_      += jet->energy();
    Ht_     += jet->et();

    if(jet - jets.begin() < 3){
      Ht123_  += jet->et();
    }
    if(jet - jets.begin() == 4 || jet - jets.begin() == 5)
      Et56_ += jet->et();

    if(jet - jets.begin() > 1){
      Ht3jet_ += jet->et();
      EtSin2Theta3jet_ += jet->et()*pow(sin(jet->theta()),2);
      theta3jet_ += (jet->theta() > M_PI/2.)? (M_PI - jet->theta()) : jet->theta();
      sinTheta3jet_ += sin(jet->theta());
    }

    TCHE_BJet_Discs_  .push_back( jet->bDiscriminator("trackCountingHighEffBJetTags")         );
    TCHP_BJet_Discs_  .push_back( jet->bDiscriminator("trackCountingHighPurBJetTags")         );
    SSVHE_BJet_Discs_ .push_back( jet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags") );
    SSVHP_BJet_Discs_ .push_back( jet->bDiscriminator("simpleSecondaryVertexHighPurBJetTags") );
    CSV_BJet_Discs_   .push_back( jet->bDiscriminator("combinedSecondaryVertexBJetTags")      );
    CSVMVA_BJet_Discs_.push_back( jet->bDiscriminator("combinedSecondaryVertexMVABJetTags")   );
    SM_BJet_Discs_    .push_back( jet->bDiscriminator("softMuonBJetTags")                     );
    
    pts_         .push_back(jet->pt());
    EtSin2Thetas_.push_back(jet->et()*pow(sin(jet->theta()),2));
    thetas_      .push_back( (jet->theta() > M_PI/2.)? (M_PI - jet->theta()) : jet->theta() );

    fourVectors.push_back(jet->p4());

    if(jet->bDiscriminator("trackCountingHighEffBJetTags")         > 3.3  ) ++TCHE_Bjets_;
    if(jet->bDiscriminator("trackCountingHighPurBJetTags")         > 3.41 ) ++TCHP_Bjets_;
    if(jet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags") > 1.74 ) ++SSVHE_Bjets_;
    if(jet->bDiscriminator("simpleSecondaryVertexHighPurBJetTags") > 2.0  ) ++SSVHP_Bjets_;
    if(jet->bDiscriminator("combinedSecondaryVertexBJetTags")      > 0.75 ) ++CSV_Bjets_;
    if(jet->bDiscriminator("combinedSecondaryVertexMVABJetTags")   > 0.75 ) ++CSVMVA_Bjets_;
    if(jet->bDiscriminator("softMuonBJetTags")                     > 0.3  ) ++SM_Bjets_;
    
    if(jet->nConstituents() > 0){
      //if( jet->daughterPtr(0).productGetter()->getIt(jet->daughterPtr(0).id()) != 0 ){
    etaetaMoments_.push_back(         jet->etaetaMoment() );
    etaphiMoments_.push_back(std::abs(jet->etaphiMoment()));
    phiphiMoments_.push_back(         jet->phiphiMoment() );
    
    jets_etaetaMoment_ +=          jet->etaetaMoment();
    jets_etaphiMoment_ += std::abs(jet->etaphiMoment());
    jets_phiphiMoment_ +=          jet->phiphiMoment();

    etaetaMomentsLogEt_.push_back(         jet->etaetaMoment() *log(jet->et()));
    etaphiMomentsLogEt_.push_back(std::abs(jet->etaphiMoment())*log(jet->et()));
    phiphiMomentsLogEt_.push_back(         jet->phiphiMoment() *log(jet->et()));

    jets_etaetaMomentLogEt_ +=          jet->etaetaMoment() *log(jet->et());
    jets_etaphiMomentLogEt_ += std::abs(jet->etaphiMoment())*log(jet->et());
    jets_phiphiMomentLogEt_ +=          jet->phiphiMoment() *log(jet->et());
 
    if(jet->bDiscriminator("trackCountingHighEffBJetTags")         < 3.3  && jet->bDiscriminator("trackCountingHighPurBJetTags")         < 1.93 &&
       jet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags") < 1.74 && jet->bDiscriminator("simpleSecondaryVertexHighPurBJetTags") < 2.0  ){

      ++nonBJets;

      etaetaMomentsNoB_.push_back(         jet->etaetaMoment() );
      etaphiMomentsNoB_.push_back(std::abs(jet->etaphiMoment()));
      phiphiMomentsNoB_.push_back(         jet->phiphiMoment() );
      
      jets_etaetaMomentNoB_ +=          jet->etaetaMoment();
      jets_etaphiMomentNoB_ += std::abs(jet->etaphiMoment());
      jets_phiphiMomentNoB_ +=          jet->phiphiMoment();
    }
    //}
    }

    for(std::vector<pat::Jet>::const_iterator jet2 = jet+1; jet2 != jets.end(); ++jet2){
      unsigned short comb2A[2] = { (unsigned short)(jet-jets.begin()) , (unsigned short)(jet2-jets.begin()) };
      std::vector<unsigned short> comb2(comb2A, comb2A + sizeof(comb2A) / sizeof(unsigned short));
      dRs.push_back( std::make_pair( deltaR( jet->phi(), jet->eta(), jet2->phi(), jet2->eta() ), comb2 ) );

      for(std::vector<pat::Jet>::const_iterator jet3 = jet2+1; jet3 != jets.end(); ++jet3){
    unsigned short comb3A[3] = { (unsigned short)(jet-jets.begin()) , (unsigned short)(jet2-jets.begin()) , (unsigned short)(jet3-jets.begin()) };
    std::vector<unsigned short> comb3(comb3A, comb3A + sizeof(comb3A) / sizeof(unsigned short));
    double dR1 = deltaR( jet ->eta(), jet ->phi(), jet2->eta(), jet2->phi() );
    double dR2 = deltaR( jet ->eta(), jet ->phi(), jet3->eta(), jet3->phi() );
    double dR3 = deltaR( jet2->eta(), jet2->phi(), jet3->eta(), jet3->phi() );
    dRs3Jets.push_back( std::make_pair( dR1+dR2+dR3, comb3 ) );
    M3s.push_back( std::make_pair( ( jet->p4() + jet2->p4() + jet3->p4() ).pt(), comb3 ) );
      }
    }

    totalSystem += jet->p4();
  }

  ROOT::Math::Boost CoMBoostTotal(totalSystem.BoostToCM());
  std::vector<reco::LeafCandidate> boostedJets;

  for(std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet){
    boostedJets.push_back(reco::LeafCandidate(jet->charge(), CoMBoostTotal(jet->p4()), jet->vertex(), jet->pdgId(), jet->status(), true));
  }

  EtSin2Theta3jet_ /= ((double)(jets.size()-3));
  theta3jet_       /= ((double)(jets.size()-3));
  sinTheta3jet_    /= ((double)(jets.size()-3));

  jets_etaetaMoment_ /= (double)jets.size();
  jets_etaphiMoment_ /= (double)jets.size();
  jets_phiphiMoment_ /= (double)jets.size();

  jets_etaetaMomentLogEt_ /= (double)jets.size();
  jets_etaphiMomentLogEt_ /= (double)jets.size();
  jets_phiphiMomentLogEt_ /= (double)jets.size();

  if(nonBJets){
    jets_etaetaMomentNoB_ /= (double)nonBJets;
    jets_etaphiMomentNoB_ /= (double)nonBJets;
    jets_phiphiMomentNoB_ /= (double)nonBJets;
  }

  for(unsigned short i = 0 ; i < etaetaMoments_.size() ; ++i){
    etaetaMomentsMoment_.push_back(etaetaMoments_.at(i)/jets_etaetaMoment_);
    etaphiMomentsMoment_.push_back(etaphiMoments_.at(i)/jets_etaphiMoment_);
    phiphiMomentsMoment_.push_back(phiphiMoments_.at(i)/jets_phiphiMoment_);

    etaetaMomentsMomentLogEt_.push_back(etaetaMomentsLogEt_.at(i)/jets_etaetaMomentLogEt_);
    etaphiMomentsMomentLogEt_.push_back(etaphiMomentsLogEt_.at(i)/jets_etaphiMomentLogEt_);
    phiphiMomentsMomentLogEt_.push_back(phiphiMomentsLogEt_.at(i)/jets_phiphiMomentLogEt_);
  }

  std::sort(dRs     .begin(), dRs     .end());
  std::sort(dRs3Jets.begin(), dRs3Jets.end());
  
  for(std::vector< std::pair< double, std::vector<unsigned short> > >::const_iterator dR = dRs.begin(); dR != dRs.end(); ++dR){
    dR_.push_back(dR->first);
    dRMass_.push_back((jets.at(dR->second.at(0)).p4()+jets.at(dR->second.at(1)).p4()).mass());

    ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > wHypo = jets.at(dR->second.at(0)).p4()+jets.at(dR->second.at(1)).p4();
    TLorentzVector wHypoHelper(wHypo.Px(), wHypo.Py(), wHypo.Pz(), wHypo.E());
    wHypoHelper.SetVectM(TVector3(wHypo.Px(), wHypo.Py(), wHypo.Pz()), 80.4);
    wHypo.SetPxPyPzE(wHypoHelper.Px(), wHypoHelper.Py(), wHypoHelper.Pz(), wHypoHelper.E());
    ROOT::Math::Boost CoMBoostWHypo(wHypo.BoostToCM());
    dRAngle_.push_back(ROOT::Math::VectorUtil::Angle(CoMBoostWHypo(jets.at(dR->second.at(0)).p4()), CoMBoostWHypo(jets.at(dR->second.at(1)).p4())));
  }

  for(std::vector< std::pair< double, std::vector<unsigned short> > >::const_iterator dR = dRs3Jets.begin(); dR != dRs3Jets.end(); ++dR){
    dR3Jets_.push_back(dR->first);
    dR3JetsMass_.push_back((jets.at(dR->second.at(0)).p4()+jets.at(dR->second.at(1)).p4()+jets.at(dR->second.at(2)).p4()).mass());
  }

  std::vector< std::pair< double, unsigned short > > massDiff2W;

  for(std::vector< double >::const_iterator mass = dRMass_.begin(); mass != dRMass_.end(); ++mass){
    massDiff2W.push_back(std::make_pair(std::abs((*mass)-80.4), mass - dRMass_.begin()));
  }

  std::sort(massDiff2W.begin(), massDiff2W.end());
  
  //std::vector<std::pair< double, std::vector<unsigned short> > > massDiff;

  for(std::vector< std::pair< double, unsigned short > >::const_iterator i = massDiff2W.begin(); i != massDiff2W.end(); ++i){
    unsigned int mass1 = i->second;
    for(std::vector< std::pair< double, unsigned short > >::const_iterator j = i + 1; j != massDiff2W.end(); ++j){
      unsigned int mass2 = j->second;
      if(dRs.at(mass1).second.at(0) != dRs.at(mass2).second.at(0) && dRs.at(mass1).second.at(0) != dRs.at(mass2).second.at(1) &&
     dRs.at(mass1).second.at(1) != dRs.at(mass2).second.at(0) && dRs.at(mass1).second.at(1) != dRs.at(mass2).second.at(1)){
    //unsigned short combA[2] = { mass1 , mass2 };
    //std::vector<unsigned short> comb(combA, combA + sizeof(combA) / sizeof(unsigned short));
    //massDiff.push_back(std::make_pair(std::abs(dRMass_.at(mass1)-dRMass_.at(mass2)), comb));
    massDiffMWCands_.push_back(std::abs(dRMass_.at(mass1)-dRMass_.at(mass2)));
      }
    }
    if(massDiffMWCands_.size() > 20) break;
  }

  //std::sort(massDiff.begin(), massDiff.end());
  /*
  for(std::vector<std::pair< double, std::vector<unsigned short> > >::const_iterator diff = massDiff.begin(); diff != massDiff.end() ; ++diff){
    std::cout << "| "   << dRMass_.at(diff->second.at(0)) << "(" << diff->second.at(0)
          << ") - " << dRMass_.at(diff->second.at(1)) << "(" << diff->second.at(1)
          << ") | = " << diff->first << std::endl;
  }
  std::cout << "---------------------------------------------" << std::endl;
  */

  std::sort(TCHE_BJet_Discs_  .begin(), TCHE_BJet_Discs_  .end());
  std::sort(TCHP_BJet_Discs_  .begin(), TCHP_BJet_Discs_  .end());
  std::sort(SSVHE_BJet_Discs_ .begin(), SSVHE_BJet_Discs_ .end());
  std::sort(SSVHP_BJet_Discs_ .begin(), SSVHP_BJet_Discs_ .end());
  std::sort(CSV_BJet_Discs_   .begin(), CSV_BJet_Discs_   .end());
  std::sort(CSVMVA_BJet_Discs_.begin(), CSVMVA_BJet_Discs_.end());
  std::sort(SM_BJet_Discs_    .begin(), SM_BJet_Discs_    .end());

  std::sort(etaetaMoments_.begin(), etaetaMoments_.end());
  std::sort(etaphiMoments_.begin(), etaphiMoments_.end());
  std::sort(phiphiMoments_.begin(), phiphiMoments_.end());

  std::sort(etaetaMomentsLogEt_.begin(), etaetaMomentsLogEt_.end());
  std::sort(etaphiMomentsLogEt_.begin(), etaphiMomentsLogEt_.end());
  std::sort(phiphiMomentsLogEt_.begin(), phiphiMomentsLogEt_.end());

  std::sort(etaetaMomentsMoment_.begin(), etaetaMomentsMoment_.end());
  std::sort(etaphiMomentsMoment_.begin(), etaphiMomentsMoment_.end());
  std::sort(phiphiMomentsMoment_.begin(), phiphiMomentsMoment_.end());

  std::sort(etaetaMomentsMomentLogEt_.begin(), etaetaMomentsMomentLogEt_.end());
  std::sort(etaphiMomentsMomentLogEt_.begin(), etaphiMomentsMomentLogEt_.end());
  std::sort(phiphiMomentsMomentLogEt_.begin(), phiphiMomentsMomentLogEt_.end());

  std::sort(etaetaMomentsNoB_.begin(), etaetaMomentsNoB_.end());
  std::sort(etaphiMomentsNoB_.begin(), etaphiMomentsNoB_.end());
  std::sort(phiphiMomentsNoB_.begin(), phiphiMomentsNoB_.end());

  std::sort(M3s.begin(), M3s.end());
  M3_ = ( jets.at((M3s.back().second.at(0))).p4() +  jets.at((M3s.back().second.at(1))).p4() + jets.at((M3s.back().second.at(2))).p4() ).mass();

  sqrt_s_ = totalSystem.mass();

  for(std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > >::const_iterator jet = fourVectors.begin(); jet != fourVectors.end(); ++jet){
    ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > boostedJet = CoMBoostTotal( *jet );
    if(jet - fourVectors.begin() > 1){
      thetaStar3jet_ += (boostedJet.Theta() > M_PI/2.)? (M_PI - boostedJet.Theta()) : boostedJet.Theta();
      sinThetaStar3jet_ += sin(boostedJet.Theta());
      EtStar3jet_ += jet->Et() * pow( sin( boostedJet.Theta() ), 2);
    }
    thetaStars_.push_back((boostedJet.Theta() > M_PI/2.)? (M_PI - boostedJet.Theta()) : boostedJet.Theta());
    EtStars_.push_back( jet->Et() * pow( sin( boostedJet.Theta() ), 2) );
  }

  theta3jet_        /= (double)fourVectors.size() - 2.;
  sinTheta3jet_     /= (double)fourVectors.size() - 2.;
  thetaStar3jet_    /= (double)fourVectors.size() - 2.;
  sinThetaStar3jet_ /= (double)fourVectors.size() - 2.;
  
  EventShapeVariables eventshape(makeVecForEventShape(jets));

  aplanarity_  = eventshape.aplanarity();
  sphericity_  = eventshape.sphericity();
  circularity_ = eventshape.circularity();
  isotropy_    = eventshape.isotropy();
  C_           = eventshape.C();
  D_           = eventshape.D();

  EventShapeVariables eventshapeAll(makeVecForEventShape(jets,false));

  aplanarityAll_  = eventshapeAll.aplanarity();
  sphericityAll_  = eventshapeAll.sphericity();
  circularityAll_ = eventshapeAll.circularity();
  isotropyAll_    = eventshapeAll.isotropy();
  CAll_           = eventshapeAll.C();
  DAll_           = eventshapeAll.D();

  EventShapeVariables eventshapeAllCMS(makeVecForEventShape(jets,false,CoMBoostTotal));

  aplanarityAllCMS_  = eventshapeAllCMS.aplanarity();
  sphericityAllCMS_  = eventshapeAllCMS.sphericity();
  circularityAllCMS_ = eventshapeAllCMS.circularity();
  isotropyAllCMS_    = eventshapeAllCMS.isotropy();
  CAllCMS_           = eventshapeAllCMS.C();
  DAllCMS_           = eventshapeAllCMS.D();

  Thrust thrustAlgo(jets.begin(), jets.end());
  thrust_ = thrustAlgo.thrust();

  Thrust thrustAlgoCMS(boostedJets.begin(), boostedJets.end());
  thrustCMS_ = thrustAlgoCMS.thrust();

}

TtFullHadSignalSel::~TtFullHadSignalSel() 
{
}