MET.cc

Go to the documentation of this file.

//
//

#include "DataFormats/PatCandidates/interface/MET.h"


using namespace pat;


MET::MET() {
  initCorMap();
}


MET::MET(const reco::MET & aMET) : PATObject<reco::MET>(aMET) {
    const reco::CaloMET * calo = dynamic_cast<const reco::CaloMET *>(&aMET);
    if (calo != nullptr) caloMET_.push_back(calo->getSpecific());
    const reco::PFMET * pf = dynamic_cast<const reco::PFMET *>(&aMET);
    if (pf != nullptr) pfMET_.push_back(pf->getSpecific());
    const pat::MET * pm = dynamic_cast<const pat::MET *>(&aMET);
    if (pm != nullptr) this->operator=(*pm);

    metSig_ =0.;
    sumPtUnclustered_ = 0.;
    initCorMap();
}


MET::MET(const edm::RefToBase<reco::MET> & aMETRef) : PATObject<reco::MET>(aMETRef) {
    const reco::CaloMET * calo = dynamic_cast<const reco::CaloMET *>(aMETRef.get());
    if (calo != nullptr) caloMET_.push_back(calo->getSpecific());
    const reco::PFMET * pf = dynamic_cast<const reco::PFMET *>(aMETRef.get());
    if (pf != nullptr) pfMET_.push_back(pf->getSpecific());
    const pat::MET * pm = dynamic_cast<const pat::MET *>(aMETRef.get());
    if (pm != nullptr) this->operator=(*pm);

    metSig_ =0.;
    sumPtUnclustered_ = 0.;
    initCorMap();
}

MET::MET(const edm::Ptr<reco::MET> & aMETRef) : PATObject<reco::MET>(aMETRef) {
    const reco::CaloMET * calo = dynamic_cast<const reco::CaloMET *>(aMETRef.get());
    if (calo != nullptr) caloMET_.push_back(calo->getSpecific());
    const reco::PFMET * pf = dynamic_cast<const reco::PFMET *>(aMETRef.get());
    if (pf != nullptr) pfMET_.push_back(pf->getSpecific());
    const pat::MET * pm = dynamic_cast<const pat::MET *>(aMETRef.get());
    if (pm != nullptr) this->operator=(*pm);

    metSig_ =0.;
    sumPtUnclustered_ = 0.;
    initCorMap();
}

MET::MET(MET const& iOther):
PATObject<reco::MET>(iOther),
genMET_(iOther.genMET_),
caloMET_(iOther.caloMET_),
pfMET_(iOther.pfMET_),
metSig_(iOther.metSig_),
sumPtUnclustered_(iOther.sumPtUnclustered_),
uncertaintiesRaw_(iOther.uncertaintiesRaw_), //74X reading compatibility
uncertaintiesType1_(iOther.uncertaintiesType1_), //74X compatibility
uncertaintiesType1p2_(iOther.uncertaintiesType1p2_), //74X compatibility
uncertainties_(iOther.uncertainties_),
corrections_(iOther.corrections_),
caloPackedMet_(iOther.caloPackedMet_) {

  initCorMap();
}

// old uncertainties discarded on purpose to avoid confusion
MET::MET(const reco::MET & corMET, const MET& srcMET ):
PATObject<reco::MET>(corMET),
genMET_(srcMET.genMET_),
caloMET_(srcMET.caloMET_),
pfMET_(srcMET.pfMET_),
metSig_(srcMET.metSig_),
sumPtUnclustered_(srcMET.sumPtUnclustered_),
caloPackedMet_(srcMET.caloPackedMet_) {

  setSignificanceMatrix(srcMET.getSignificanceMatrix());

  initCorMap();
}

MET::~MET() {

}

MET& MET::operator=(MET const& iOther) {
   PATObject<reco::MET>::operator=(iOther);
   genMET_ = iOther.genMET_;
   caloMET_ =iOther.caloMET_;
   pfMET_ =iOther.pfMET_;
   uncertaintiesRaw_ = iOther.uncertaintiesRaw_; //74X compatibility
   uncertaintiesType1_ = iOther.uncertaintiesType1_;
   uncertaintiesType1p2_ = iOther.uncertaintiesType1p2_;
   uncertainties_ = iOther.uncertainties_;
   corrections_ = iOther.corrections_;
   metSig_ = iOther.metSig_;
   sumPtUnclustered_ = iOther.sumPtUnclustered_;
   caloPackedMet_ = iOther.caloPackedMet_;

   return *this;
}

const reco::GenMET * MET::genMET() const {
  return (!genMET_.empty() ? &genMET_.front() : nullptr );
}

void MET::setGenMET(const reco::GenMET & gm) {
  genMET_.clear();
  genMET_.push_back(gm);
}


//Method to set the MET significance
void MET::setMETSignificance(const double& metSig) {
  metSig_ = metSig;
}

double MET::metSignificance() const {
  return metSig_;
}

void MET::setMETSumPtUnclustered(const double &sumPtUnclustered) { sumPtUnclustered_ = sumPtUnclustered; }

double MET::metSumPtUnclustered() const { return sumPtUnclustered_; }

void
MET::initCorMap() {

  std::vector<MET::METCorrectionType> tmpRaw;
  std::vector<MET::METCorrectionType> tmpType1;
  std::vector<MET::METCorrectionType> tmpType01;
  std::vector<MET::METCorrectionType> tmpTypeXY;
  std::vector<MET::METCorrectionType> tmpType1XY;
  std::vector<MET::METCorrectionType> tmpType01XY;
  std::vector<MET::METCorrectionType> tmpType1Smear;
  std::vector<MET::METCorrectionType> tmpType01Smear;
  std::vector<MET::METCorrectionType> tmpType1SmearXY;
  std::vector<MET::METCorrectionType> tmpType01SmearXY;
  
  tmpRaw.push_back(MET::None);
  
  tmpType1.push_back(MET::T1);
  tmpType01.push_back(MET::T1);
  tmpType1XY.push_back(MET::T1);
  tmpType01XY.push_back(MET::T1);
  tmpType1Smear.push_back(MET::T1);
  tmpType01Smear.push_back(MET::T1);
  tmpType1SmearXY.push_back(MET::T1);
  tmpType01SmearXY.push_back(MET::T1);

  tmpType01.push_back(MET::T0);
  tmpType01XY.push_back(MET::T0);
  tmpType01Smear.push_back(MET::T0);
  tmpType01SmearXY.push_back(MET::T0);
  
  tmpType1Smear.push_back(MET::Smear);
  tmpType01Smear.push_back(MET::Smear);
  tmpType1SmearXY.push_back(MET::Smear);
  tmpType01SmearXY.push_back(MET::Smear);

  tmpTypeXY.push_back(MET::TXYForRaw);
  tmpType1XY.push_back(MET::TXY);
  tmpType01XY.push_back(MET::TXYForT01);
  tmpType1SmearXY.push_back(MET::TXYForT1Smear);
  tmpType01SmearXY.push_back(MET::TXYForT01Smear);

  corMap_[MET::Raw] = tmpRaw;
  corMap_[MET::Type1] = tmpType1;
  corMap_[MET::Type01] = tmpType01;
  corMap_[MET::TypeXY] = tmpTypeXY;
  corMap_[MET::Type1XY] = tmpType1XY;
  corMap_[MET::Type01XY] = tmpType01XY;
  corMap_[MET::Type1Smear] = tmpType1Smear;
  corMap_[MET::Type01Smear] = tmpType01Smear;
  corMap_[MET::Type1SmearXY] = tmpType1SmearXY;
  corMap_[MET::Type01SmearXY] = tmpType01SmearXY;
  
  //specific calo case
  std::vector<MET::METCorrectionType> tmpRawCalo;
  tmpRawCalo.push_back(MET::Calo);
  corMap_[MET::RawCalo] = tmpRawCalo;

  //specific chs case
  std::vector<MET::METCorrectionType> tmpRawChs;
  tmpRawChs.push_back(MET::Chs);
  corMap_[MET::RawChs] = tmpRawChs;

  //specific trk case
  std::vector<MET::METCorrectionType> tmpRawTrk;
  tmpRawTrk.push_back(MET::Trk);
  corMap_[MET::RawTrk] = tmpRawTrk;
}

const MET::PackedMETUncertainty
MET::findMETTotalShift(MET::METCorrectionLevel cor, MET::METUncertainty shift) const {

  //find corrections shifts =============================
  std::map<MET::METCorrectionLevel, std::vector<MET::METCorrectionType> >::const_iterator itCor_ = corMap_.find(cor);
  if(itCor_==corMap_.end() ) throw cms::Exception("Unsupported", "Specified MET correction scheme does not exist");

  bool isSmeared=false;
  MET::PackedMETUncertainty totShift;
  unsigned int scor=itCor_->second.size();
  for(unsigned int i=0; i<scor;i++) {
    totShift.add( corrections_[ itCor_->second[i] ].dpx(),
          corrections_[ itCor_->second[i] ].dpy(),
          corrections_[ itCor_->second[i] ].dsumEt() );

    if(itCor_->first>=MET::Type1Smear)
      isSmeared=true;
  }

  

  //find uncertainty shift =============================

  if (uncertainties_.empty())
      return totShift;

  if(shift>=MET::METUncertaintySize) throw cms::Exception("Unsupported", "MET uncertainty does not exist");
  if(isSmeared && shift<=MET::JetResDown) shift = (MET::METUncertainty)(MET::METUncertaintySize+shift+1);
                              
  totShift.add( uncertainties_[ shift ].dpx(),
        uncertainties_[ shift ].dpy(),
        uncertainties_[ shift ].dsumEt() );

  return totShift;
}


MET::Vector2 MET::shiftedP2(MET::METUncertainty shift, MET::METCorrectionLevel cor)  const {
  
  Vector2 vo;

  //backward compatibility with 74X samples -> the only one
  // with uncertaintiesType1_/uncertaintiesRaw_ not empty
  //will be removed once 74X is not used anymore
  if(!uncertaintiesType1_.empty() || !uncertaintiesRaw_.empty()) {
    if(cor!=MET::METCorrectionLevel::RawCalo) {
      vo = shiftedP2_74x(shift, cor);
    } else {
      Vector2 ret{ caloPackedMet_.dpx(), caloPackedMet_.dpy() };
      vo = ret;
    }
  }
  else {
    const MET::PackedMETUncertainty& v = findMETTotalShift(cor,shift);
    Vector2 ret{ (px() + v.dpx()), (py() + v.dpy()) };
    //return ret;
    vo = ret;
  }
  return vo;
}
MET::Vector MET::shiftedP3(MET::METUncertainty shift, MET::METCorrectionLevel cor)  const {

  Vector vo;

  //backward compatibility with 74X samples -> the only one
  // with uncertaintiesType1_/uncertaintiesRaw_ not empty
  //will be removed once 74X is not used anymore
  if(!uncertaintiesType1_.empty() || !uncertaintiesRaw_.empty()) {
    if(cor!=MET::METCorrectionLevel::RawCalo) {
      vo = shiftedP3_74x(shift, cor);
    } else {
      Vector tmp(caloPackedMet_.dpx(), caloPackedMet_.dpy(), 0);
      vo = tmp;
    }
  }
  else {
    const MET::PackedMETUncertainty& v = findMETTotalShift(cor,shift);
    //return Vector(px() + v.dpx(), py() + v.dpy(), 0);
    Vector tmp(px() + v.dpx(), py() + v.dpy(), 0);
    vo = tmp;
  }
  return vo;
}
MET::LorentzVector MET::shiftedP4(METUncertainty shift, MET::METCorrectionLevel cor)  const {

  LorentzVector vo;

  //backward compatibility with 74X samples -> the only one
  // with uncertaintiesType1_/uncertaintiesRaw_ not empty
  //will be removed once 74X is not used anymore
  if(!uncertaintiesType1_.empty() || !uncertaintiesRaw_.empty()) {
    if(cor!=MET::METCorrectionLevel::RawCalo) {
      vo = shiftedP4_74x(shift, cor);
    } else {
      double x = caloPackedMet_.dpx(), y = caloPackedMet_.dpy();
      LorentzVector tmp(x, y, 0, std::hypot(x,y));
      vo =  tmp;
    }
  }
  else {
    const MET::PackedMETUncertainty& v = findMETTotalShift(cor,shift);
    double x = px() + v.dpx(), y = py() + v.dpy();
    //return LorentzVector(x, y, 0, std::hypot(x,y));
    LorentzVector tmp(x, y, 0, std::hypot(x,y));
    vo =  tmp;
  }
  return vo;
}
double MET::shiftedSumEt(MET::METUncertainty shift, MET::METCorrectionLevel cor) const {

  double sumEto;

  //backward compatibility with 74X samples -> the only one
  // with uncertaintiesType1_/uncertaintiesRaw_ not empty
  //will be removed once 74X is not used anymore
  if(!uncertaintiesType1_.empty() || !uncertaintiesRaw_.empty()) {
    if(cor!=MET::METCorrectionLevel::RawCalo) {
      sumEto = shiftedSumEt_74x(shift, cor);
    } else {
      sumEto = caloPackedMet_.dsumEt();
    }
  }
  else {
    const MET::PackedMETUncertainty& v = findMETTotalShift(cor,shift);
    //return sumEt() + v.dsumEt();
    sumEto = sumEt() + v.dsumEt();
  }
  return sumEto;
}

MET::Vector2 MET::corP2(MET::METCorrectionLevel cor)  const {
  return shiftedP2(MET::NoShift, cor );
}
MET::Vector MET::corP3(MET::METCorrectionLevel cor)  const {
  return shiftedP3(MET::NoShift, cor );
}
MET::LorentzVector MET::corP4(MET::METCorrectionLevel cor)  const {
  return shiftedP4(MET::NoShift, cor );
}
double MET::corSumEt(MET::METCorrectionLevel cor) const {
  return shiftedSumEt(MET::NoShift, cor );
}

MET::Vector2 MET::uncorP2()  const {
  return shiftedP2(MET::NoShift, MET::Raw );
}
MET::Vector MET::uncorP3()  const {
  return shiftedP3(MET::NoShift, MET::Raw );
}
MET::LorentzVector MET::uncorP4()  const {
  return shiftedP4(MET::NoShift, MET::Raw );
}
double MET::uncorSumEt() const {
  return shiftedSumEt(MET::NoShift, MET::Raw );
}


void MET::setUncShift(double px, double py, double sumEt, METUncertainty shift, bool isSmeared) {
  if (uncertainties_.empty()) { 
    uncertainties_.resize(METUncertainty::METFullUncertaintySize);
  }

  if(isSmeared && shift<=MET::JetResDown) {
    //changing reference to only get the uncertainty shift and not the smeared one
    // which is performed independently
    shift = (MET::METUncertainty)(METUncertainty::METUncertaintySize+shift+1);
    const PackedMETUncertainty& ref = uncertainties_[METUncertainty::NoShift];
    uncertainties_[shift].set(px + ref.dpx() - this->px(), py + ref.dpy() - this->py(), sumEt + ref.dsumEt() - this->sumEt() );
  }
  else
    uncertainties_[shift].set(px - this->px(), py - this->py(), sumEt - this->sumEt());
  
}

void MET::setCorShift(double px, double py, double sumEt, MET::METCorrectionType level) {
  if (corrections_.empty()) {
    corrections_.resize(MET::METCorrectionType::METCorrectionTypeSize);
  }
  
  corrections_[level].set(px - this->px(), py - this->py(), sumEt - this->sumEt());
}


MET::Vector2 MET::caloMETP2() const {
  return shiftedP2(MET::METUncertainty::NoShift, MET::METCorrectionLevel::RawCalo );
}

double MET::caloMETPt() const {
  return caloMETP2().pt();
}

double MET::caloMETPhi() const {
  return caloMETP2().phi();
}

double MET::caloMETSumEt() const {
  return shiftedSumEt(MET::NoShift, MET::RawCalo );
}

// functions to access to 74X samples ========================================================
MET::Vector2 MET::shiftedP2_74x(MET::METUncertainty shift, MET::METCorrectionLevel level)  const {
  if (level != Type1 && level != Raw) throw cms::Exception("Unsupported", "MET uncertainties only supported for Raw and Type1 in 74X samples \n");
    const std::vector<PackedMETUncertainty> &v = (level == Type1 ? uncertaintiesType1_ :  uncertaintiesRaw_);
    if (v.empty()) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type\n");
    if (v.size() == 1) {
        if (shift != MET::METUncertainty::NoShift) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type (only central value available)\n");
        return Vector2{ (px() + v.front().dpx()), (py() + v.front().dpy()) };
    }
    Vector2 ret{ (px() + v[shift].dpx()), (py() + v[shift].dpy()) };
    return ret;
}

MET::Vector MET::shiftedP3_74x(MET::METUncertainty shift, MET::METCorrectionLevel level)  const {
  if (level != Type1 && level != Raw) throw cms::Exception("Unsupported", "MET uncertainties only supported for Raw and Type1 in 74X samples \n");
    const std::vector<PackedMETUncertainty> &v = (level == Type1 ? uncertaintiesType1_ :  uncertaintiesRaw_);
    if (v.empty()) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type\n");
    if (v.size() == 1) {
        if (shift != MET::METUncertainty::NoShift) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type (only central value available)\n");
        return Vector(px() + v.front().dpx(), py() + v.front().dpy(), 0);
    }
    return Vector(px() + v[shift].dpx(), py() + v[shift].dpy(), 0);
}

MET::LorentzVector MET::shiftedP4_74x(METUncertainty shift, MET::METCorrectionLevel level)  const {
  if (level != Type1 && level != Raw) throw cms::Exception("Unsupported", "MET uncertainties only supported for Raw and Type1 in 74X samples\n");
    const std::vector<PackedMETUncertainty> &v = (level == Type1 ? uncertaintiesType1_ : uncertaintiesRaw_);
    if (v.empty()) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type\n");
    if (v.size() == 1) {
        if (shift != MET::METUncertainty::NoShift) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type (only central value available)\n");
        double x = px() + v.front().dpx(), y = py() + v.front().dpy();
        return LorentzVector(x, y, 0, std::hypot(x,y));
    }
    double x = px() + v[shift].dpx(), y = py() + v[shift].dpy();
    return LorentzVector(x, y, 0, std::hypot(x,y));
}

double MET::shiftedSumEt_74x(MET::METUncertainty shift, MET::METCorrectionLevel level) const {
  if (level != Type1 && level != Raw) throw cms::Exception("Unsupported", "MET uncertainties only supported for Raw and Type1 in 74X samples\n");
    const std::vector<PackedMETUncertainty> &v = (level == Type1 ? uncertaintiesType1_ : uncertaintiesRaw_);
    if (v.empty()) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type\n");
    if (v.size() == 1) {
        if (shift != MET::METUncertainty::NoShift) throw cms::Exception("Unsupported", "MET uncertainties not available for the specified correction type (only central value available)\n");
        return sumEt() + v.front().dsumEt();
    }
    return sumEt() + v[shift].dsumEt();
}




#include "DataFormats/Math/interface/libminifloat.h"

void MET::PackedMETUncertainty::pack() {
  packedDpx_  =  MiniFloatConverter::float32to16(dpx_);
  packedDpy_  =  MiniFloatConverter::float32to16(dpy_);
  packedDSumEt_  =  MiniFloatConverter::float32to16(dsumEt_);
}
void  MET::PackedMETUncertainty::unpack() const {
  unpacked_=true;
  dpx_=MiniFloatConverter::float16to32(packedDpx_);
  dpy_=MiniFloatConverter::float16to32(packedDpy_);
  dsumEt_=MiniFloatConverter::float16to32(packedDSumEt_);
}