PackedCandidate.cc

Go to the documentation of this file.

#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/PatCandidates/interface/libminifloat.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include "DataFormats/PatCandidates/interface/liblogintpack.h"
using namespace logintpack;

CovarianceParameterization pat::PackedCandidate::covarianceParameterization_;
std::once_flag pat::PackedCandidate::covariance_load_flag;

void pat::PackedCandidate::pack(bool unpackAfterwards) {
    packedPt_  =  MiniFloatConverter::float32to16(p4_.load()->Pt());
    packedEta_ =  int16_t(std::round(p4_.load()->Eta()/6.0f*std::numeric_limits<int16_t>::max()));
    packedPhi_ =  int16_t(std::round(p4_.load()->Phi()/3.2f*std::numeric_limits<int16_t>::max()));
    packedM_   =  MiniFloatConverter::float32to16(p4_.load()->M());
    if (unpackAfterwards) {
      delete p4_.exchange(nullptr);
      delete p4c_.exchange(nullptr);
      unpack(); // force the values to match with the packed ones
    }
}

void pat::PackedCandidate::packVtx(bool unpackAfterwards) {
    reco::VertexRef pvRef = vertexRef();
    Point pv = pvRef.isNonnull() ? pvRef->position() : Point();
    float dxPV = vertex_.load()->X() - pv.X(), dyPV = vertex_.load()->Y() - pv.Y(); //, rPV = std::hypot(dxPV, dyPV);
    float s = std::sin(float(p4_.load()->Phi())+dphi_), c = std::cos(float(p4_.load()->Phi()+dphi_)); // not the fastest option, but we're in reduced precision already, so let's avoid more roundoffs
    dxy_  = - dxPV * s + dyPV * c;    
    // if we want to go back to the full x,y,z we need to store also
    // float dl = dxPV * c + dyPV * s; 
    // float xRec = - dxy_ * s + dl * c, yRec = dxy_ * c + dl * s;
    float pzpt = p4_.load()->Pz()/p4_.load()->Pt();
    dz_ = vertex_.load()->Z() - pv.Z() - (dxPV*c + dyPV*s) * pzpt;
    packedDxy_ = MiniFloatConverter::float32to16(dxy_*100);
    packedDz_   = pvRef.isNonnull() ? MiniFloatConverter::float32to16(dz_*100) : int16_t(std::round(dz_/40.f*std::numeric_limits<int16_t>::max()));
    packedDPhi_ =  int16_t(std::round(dphi_/3.2f*std::numeric_limits<int16_t>::max()));
    packedDEta_ =  MiniFloatConverter::float32to16(deta_);
    packedDTrkPt_ = MiniFloatConverter::float32to16(dtrkpt_);
 
    if (unpackAfterwards) {
      delete vertex_.exchange(nullptr);
      unpackVtx();
    }
}

void pat::PackedCandidate::unpack() const {
    float pt = MiniFloatConverter::float16to32(packedPt_);
    double shift = (pt<1. ? 0.1*pt : 0.1/pt); // shift particle phi to break degeneracies in angular separations
    double sign = ( ( int(pt*10) % 2 == 0 ) ? 1 : -1 ); // introduce a pseudo-random sign of the shift
    double phi = int16_t(packedPhi_)*3.2f/std::numeric_limits<int16_t>::max() + sign*shift*3.2/std::numeric_limits<int16_t>::max();
    auto p4 = std::make_unique<PolarLorentzVector>(pt,
                             int16_t(packedEta_)*6.0f/std::numeric_limits<int16_t>::max(),
                             phi,
                             MiniFloatConverter::float16to32(packedM_));
    auto p4c = std::make_unique<LorentzVector>( *p4 );
    PolarLorentzVector* expectp4= nullptr;
    if( p4_.compare_exchange_strong(expectp4,p4.get()) ) {
      p4.release();
    }

    //p4c_ works as the guard for unpacking so it
    // must be set last
    LorentzVector* expectp4c = nullptr;
    if(p4c_.compare_exchange_strong(expectp4c, p4c.get()) ) {
      p4c.release();
    }
}

void pat::PackedCandidate::packCovariance(const reco::TrackBase::CovarianceMatrix &m, bool unpackAfterwards){
    packedCovariance_.dptdpt = packCovarianceElement(m,0,0);
    packedCovariance_.detadeta = packCovarianceElement(m,1,1);
    packedCovariance_.dphidphi = packCovarianceElement(m,2,2); 
    packedCovariance_.dxydxy =packCovarianceElement(m,3,3);
    packedCovariance_.dzdz = packCovarianceElement(m,4,4);
    packedCovariance_.dxydz = packCovarianceElement(m,3,4);
    packedCovariance_.dlambdadz = packCovarianceElement(m,1,4);
    packedCovariance_.dphidxy = packCovarianceElement(m,2,3);
   //unpack afterwards
   if(unpackAfterwards) unpackCovariance();
}

void pat::PackedCandidate::unpackCovariance() const {
    const CovarianceParameterization & p=covarianceParameterization();
    if(p.isValid()) 
    {
      auto m = std::make_unique<reco::TrackBase::CovarianceMatrix>() ;
      for(int i=0;i<5;i++)
        for(int j=0;j<5;j++){
          (*m)(i,j)=0;
      }
      unpackCovarianceElement(*m,packedCovariance_.dptdpt,0,0);
      unpackCovarianceElement(*m,packedCovariance_.detadeta,1,1);
      unpackCovarianceElement(*m,packedCovariance_.dphidphi,2,2);
      unpackCovarianceElement(*m,packedCovariance_.dxydxy,3,3);
      unpackCovarianceElement(*m,packedCovariance_.dzdz,4,4);
      unpackCovarianceElement(*m,packedCovariance_.dxydz,3,4);
      unpackCovarianceElement(*m,packedCovariance_.dlambdadz,1,4);
      unpackCovarianceElement(*m,packedCovariance_.dphidxy,2,3);
      reco::TrackBase::CovarianceMatrix* expected = nullptr;
      if( m_.compare_exchange_strong(expected,m.get()) ) {
         m.release();
      }

    } else {
     throw edm::Exception(edm::errors::UnimplementedFeature)
     << "You do not have a valid track parameters file loaded. "
     << "Please check that the release version is compatible with your input data"
     <<"or avoid accessing track parameter uncertainties. ";
    }
}

void pat::PackedCandidate::unpackVtx() const {
    reco::VertexRef pvRef = vertexRef();
    dphi_ = int16_t(packedDPhi_)*3.2f/std::numeric_limits<int16_t>::max(),
    deta_ = MiniFloatConverter::float16to32(packedDEta_);
    dtrkpt_ = MiniFloatConverter::float16to32(packedDTrkPt_);
    dxy_ = MiniFloatConverter::float16to32(packedDxy_)/100.;
    dz_   = pvRef.isNonnull() ? MiniFloatConverter::float16to32(packedDz_)/100. : int16_t(packedDz_)*40.f/std::numeric_limits<int16_t>::max();
    Point pv = pvRef.isNonnull() ? pvRef->position() : Point();
    float phi = p4_.load()->Phi()+dphi_, s = std::sin(phi), c = std::cos(phi);
    auto vertex = std::make_unique<Point>(pv.X() - dxy_ * s,
                    pv.Y() + dxy_ * c,
                    pv.Z() + dz_ ); // for our choice of using the PCA to the PV, by definition the remaining term -(dx*cos(phi) + dy*sin(phi))*(pz/pt) is zero

     

 
    Point* expected = nullptr;
    if( vertex_.compare_exchange_strong(expected,vertex.get()) ) {
      vertex.release();
    }
}

pat::PackedCandidate::~PackedCandidate() { 
  delete p4_.load();
  delete p4c_.load();
  delete vertex_.load();
  delete track_.load();
  delete m_.load();
}


float pat::PackedCandidate::dxy(const Point &p) const {
    maybeUnpackBoth();
    const float phi = float(p4_.load()->Phi())+dphi_;
    return -(vertex_.load()->X()-p.X()) * std::sin(phi) + (vertex_.load()->Y()-p.Y()) * std::cos(phi);
}
float pat::PackedCandidate::dz(const Point &p) const {
    maybeUnpackBoth();
    const float phi = float(p4_.load()->Phi())+dphi_;
    return (vertex_.load()->Z()-p.Z())  - ((vertex_.load()->X()-p.X()) * std::cos(phi) + (vertex_.load()->Y()-p.Y()) * std::sin(phi)) * p4_.load()->Pz()/p4_.load()->Pt();
}

void pat::PackedCandidate::unpackTrk() const {
    maybeUnpackBoth();
    math::RhoEtaPhiVector p3(ptTrk(),etaAtVtx(),phiAtVtx());
    maybeUnpackCovariance();
    int numberOfStripLayers = stripLayersWithMeasurement(), numberOfPixelLayers = pixelLayersWithMeasurement();
    int numberOfPixelHits = this->numberOfPixelHits();
    int numberOfHits = this->numberOfHits();

    int ndof = numberOfHits+numberOfPixelHits-5;
    LostInnerHits innerLost = lostInnerHits();
    
    auto track = std::make_unique<reco::Track>(normalizedChi2_*ndof,ndof,*vertex_,math::XYZVector(p3.x(),p3.y(),p3.z()),charge(),*(m_.load()),reco::TrackBase::undefAlgorithm,reco::TrackBase::loose);
    int i=0;
    if ( firstHit_ == 0) { //Backward compatible 
       if(innerLost == validHitInFirstPixelBarrelLayer){
          track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 1, 0, TrackingRecHit::valid);       
          i=1;
       } 
    } else {
       track->appendHitPattern(firstHit_,TrackingRecHit::valid);
    }

    if(firstHit_!=0 && reco::HitPattern::pixelHitFilter(firstHit_)) i=1;

    // add hits to match the number of laters and validHitInFirstPixelBarrelLayer
    if(innerLost == validHitInFirstPixelBarrelLayer){
        // then to encode the number of layers, we add more hits on distinct layers (B2, B3, B4, F1, ...)
        for(; i<numberOfPixelLayers; i++) {
            if (i <= 3) { 
                track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, i+1, 0, TrackingRecHit::valid); 
            } else {    
                track->appendTrackerHitPattern(PixelSubdetector::PixelEndcap, i-3, 0, TrackingRecHit::valid); 
            }
        }
    } else {
        // to encode the information on the layers, we add one valid hits per layer but skipping PXB1
        int iOffset=0;  
        if(firstHit_!=0 && reco::HitPattern::pixelHitFilter(firstHit_)) {
         iOffset=reco::HitPattern::getLayer(firstHit_);
         if(reco::HitPattern::getSubStructure(firstHit_)==PixelSubdetector::PixelEndcap) iOffset+=3;
    } else {iOffset=1; }
        for(;i<numberOfPixelLayers; i++) {
            if (i+iOffset <= 2 ) { track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, i+iOffset+1, 0, TrackingRecHit::valid);   }
        else {track->appendTrackerHitPattern(PixelSubdetector::PixelEndcap, i+iOffset-3+1, 0, TrackingRecHit::valid);  }

        }
    }
    // add extra hits (overlaps, etc), all on the first layer with a hit - to avoid increasing the layer count
    for(;i<numberOfPixelHits; i++) { 
       if(firstHit_ !=0 && reco::HitPattern::pixelHitFilter(firstHit_)) { 
          track->appendTrackerHitPattern(reco::HitPattern::getSubStructure(firstHit_), reco::HitPattern::getLayer(firstHit_), 0, TrackingRecHit::valid); 
    } else {
          track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, (innerLost == validHitInFirstPixelBarrelLayer ? 1 : 2), 0, TrackingRecHit::valid); 
    }
    }
    // now start adding strip layers, putting one hit on each layer so that the hitPattern.stripLayersWithMeasurement works.
    // we don't know what the layers where, so we just start with TIB (4 layers), then TOB (6 layers), then TEC (9)
    // and then TID(3), so that we can get a number of valid strip layers up to 4+6+9+3
    if(firstHit_!=0 && reco::HitPattern::stripHitFilter(firstHit_)) i+=1;
    int slOffset=0; 
    if(firstHit_!=0 && reco::HitPattern::stripHitFilter(firstHit_)) { 
         slOffset=reco::HitPattern::getLayer(firstHit_)-1;
         if(reco::HitPattern::getSubStructure(firstHit_)==StripSubdetector::TID) slOffset+=4;
         if(reco::HitPattern::getSubStructure(firstHit_)==StripSubdetector::TOB) slOffset+=7;
         if(reco::HitPattern::getSubStructure(firstHit_)==StripSubdetector::TEC) slOffset+=13;
    }
    for(int sl=slOffset; sl < numberOfStripLayers+slOffset; ++sl, ++i) {
        if      (sl < 4)    track->appendTrackerHitPattern(StripSubdetector::TIB,   sl   +1, 1, TrackingRecHit::valid);
        else if (sl < 4+3)  track->appendTrackerHitPattern(StripSubdetector::TID, (sl- 4)+1, 1, TrackingRecHit::valid);
        else if (sl < 7+6) track->appendTrackerHitPattern(StripSubdetector::TOB, (sl-7)+1, 1, TrackingRecHit::valid);
        else if (sl < 13+9) track->appendTrackerHitPattern(StripSubdetector::TEC, (sl-13)+1, 1, TrackingRecHit::valid);
        else break; // wtf?
    }
    // finally we account for extra strip hits beyond the one-per-layer added above. we put them all on TIB1,
    // to avoid incrementing the number of layersWithMeasurement.
    for(;i<numberOfHits;i++) {
      if(reco::HitPattern::stripHitFilter(firstHit_)) { 
              track->appendTrackerHitPattern(reco::HitPattern::getSubStructure(firstHit_),  reco::HitPattern::getLayer(firstHit_), 1, TrackingRecHit::valid);
      } else {
              track->appendTrackerHitPattern(StripSubdetector::TIB, 1, 1, TrackingRecHit::valid);
      }
    }


    switch (innerLost) {
        case validHitInFirstPixelBarrelLayer:
            break;
        case noLostInnerHits:
            break;
        case oneLostInnerHit:
            track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 1, 0, TrackingRecHit::missing_inner);
            break;
        case moreLostInnerHits:
            track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 1, 0, TrackingRecHit::missing_inner);
            track->appendTrackerHitPattern(PixelSubdetector::PixelBarrel, 2, 0, TrackingRecHit::missing_inner);
            break;
    };

    if (trackHighPurity()) track->setQuality(reco::TrackBase::highPurity);
    
    reco::Track* expected = nullptr;
    if( track_.compare_exchange_strong(expected,track.get()) ) {
      track.release();
    }


}


const reco::CandidateBaseRef & pat::PackedCandidate::masterClone() const {
  throw cms::Exception("Invalid Reference")
    << "this Candidate has no master clone reference."
    << "Can't call masterClone() method.\n";
}

bool pat::PackedCandidate::hasMasterClone() const {
  return false;
}

bool pat::PackedCandidate::hasMasterClonePtr() const {
  return false;
}


const reco::CandidatePtr & pat::PackedCandidate::masterClonePtr() const {
  throw cms::Exception("Invalid Reference")
    << "this Candidate has no master clone ptr."
    << "Can't call masterClonePtr() method.\n";
}

size_t pat::PackedCandidate::numberOfDaughters() const { 
  return 0; 
}

size_t pat::PackedCandidate::numberOfMothers() const { 
  return 0; 
}

bool pat::PackedCandidate::overlap( const reco::Candidate & o ) const { 
  return  p4() == o.p4() && vertex() == o.vertex() && charge() == o.charge();
//  return  p4() == o.p4() && charge() == o.charge();
}

const reco::Candidate * pat::PackedCandidate::daughter( size_type ) const {
  return nullptr;
}

const reco::Candidate * pat::PackedCandidate::mother( size_type ) const {
  return nullptr;
}

const reco::Candidate * pat::PackedCandidate::daughter(const std::string&) const {
  throw edm::Exception(edm::errors::UnimplementedFeature)
    << "This Candidate type does not implement daughter(std::string). "
    << "Please use CompositeCandidate or NamedCompositeCandidate.\n";
}

reco::Candidate * pat::PackedCandidate::daughter(const std::string&) {
  throw edm::Exception(edm::errors::UnimplementedFeature)
    << "This Candidate type does not implement daughter(std::string). "
    << "Please use CompositeCandidate or NamedCompositeCandidate.\n";
}



reco::Candidate * pat::PackedCandidate::daughter( size_type ) {
  return nullptr;
}

double pat::PackedCandidate::vertexChi2() const {
  return 0;
}

double pat::PackedCandidate::vertexNdof() const {
  return 0;
}

double pat::PackedCandidate::vertexNormalizedChi2() const {
  return 0;
}

double pat::PackedCandidate::vertexCovariance(int i, int j) const {
  throw edm::Exception(edm::errors::UnimplementedFeature)
    << "reco::ConcreteCandidate does not implement vertex covariant matrix.\n";
}

void pat::PackedCandidate::fillVertexCovariance(CovarianceMatrix & err) const {
  throw edm::Exception(edm::errors::UnimplementedFeature)
    << "reco::ConcreteCandidate does not implement vertex covariant matrix.\n";
}


bool pat::PackedCandidate::longLived() const {return false;}

bool pat::PackedCandidate::massConstraint() const {return false;}

// puppiweight
void pat::PackedCandidate::setPuppiWeight(float p, float p_nolep) {
  // Set both weights at once to avoid misconfigured weights if called in the wrong order
  packedPuppiweight_ = pack8logClosed((p-0.5)*2,-2,0,64);
  packedPuppiweightNoLepDiff_ = pack8logClosed((p_nolep-0.5)*2,-2,0,64) - packedPuppiweight_;
}

float pat::PackedCandidate::puppiWeight() const { return unpack8logClosed(packedPuppiweight_,-2,0,64)/2. + 0.5;}

float pat::PackedCandidate::puppiWeightNoLep() const { return unpack8logClosed(packedPuppiweightNoLepDiff_+packedPuppiweight_,-2,0,64)/2. + 0.5;}

void pat::PackedCandidate::setRawCaloFraction(float p) {
  if(100*p>std::numeric_limits<uint8_t>::max())
    rawCaloFraction_ = std::numeric_limits<uint8_t>::max(); // Set to overflow value
  else
    rawCaloFraction_ = 100*p;
}

void pat::PackedCandidate::setHcalFraction(float p) {
  hcalFraction_ = 100*p;
}

void pat::PackedCandidate::setIsIsolatedChargedHadron(bool p) {
  isIsolatedChargedHadron_ = p;
}

void pat::PackedCandidate::setDTimeAssociatedPV(float aTime, float aTimeError) {
    if (aTime == 0 && aTimeError == 0) {
        packedTime_ = 0; packedTimeError_ = 0;
    } else if (aTimeError == 0) {
        packedTimeError_ = 0;
        packedTime_ = packTimeNoError(aTime);
    } else {
        packedTimeError_ = packTimeError(aTimeError);
        aTimeError = unpackTimeError(packedTimeError_); // for reproducibility
        packedTime_ = packTimeWithError(aTime, aTimeError);
    }
}

uint8_t pat::PackedCandidate::packTimeError(float timeError) {
    if (timeError <= 0) return 0;
    // log-scale packing.
    // for MIN_TIMEERROR = 0.002, EXPO_TIMEERROR = 5:
    //      minimum value 0.002 = 2ps (packed as 1)
    //      maximum value 0.5 ns      (packed as 255)
    //      constant *relative* precision of about 2%
    return std::max<uint8_t>( std::min(std::round(std::ldexp(std::log2(timeError/MIN_TIMEERROR), +EXPO_TIMEERROR)), 255.f), 1);
}
float pat::PackedCandidate::unpackTimeError(uint8_t timeError) {
    return timeError > 0 ? MIN_TIMEERROR * std::exp2(std::ldexp(float(timeError),-EXPO_TIMEERROR)) : -1.0f;
}
float pat::PackedCandidate::unpackTimeNoError(int16_t time) {
    if (time == 0) return 0.f;
    return (time > 0 ? MIN_TIME_NOERROR : -MIN_TIME_NOERROR) * std::exp2(std::ldexp(float(std::abs(time)),-EXPO_TIME_NOERROR));
}
int16_t pat::PackedCandidate::packTimeNoError(float time) {
    // encoding in log scale to store times in a large range with few bits.
    // for MIN_TIME_NOERROR = 0.0002 and EXPO_TIME_NOERROR = 6:
    //    smallest non-zero time = 0.2 ps (encoded as +/-1)
    //    one BX, +/- 12.5 ns, is fully covered with 11 bits (+/- 1023)
    //    12 bits cover by far any plausible value (+/-2047 corresponds to about +/- 0.8 ms!)
    //    constant *relative* ~1% precision
    if (std::abs(time) < MIN_TIME_NOERROR) return 0; // prevent underflows
    float fpacked = std::ldexp(std::log2(std::abs(time/MIN_TIME_NOERROR)),+EXPO_TIME_NOERROR);
    return (time > 0 ? +1 : -1)*std::min(std::round(fpacked), 2047.f);
}
float pat::PackedCandidate::unpackTimeWithError(int16_t time, uint8_t timeError) {
    if (time % 2 == 0) {
        // no overflow: drop rightmost bit and unpack in units of timeError
        return std::ldexp(unpackTimeError(timeError), EXPO_TIME_WITHERROR) * float(time/2);
    } else {
        // overflow: drop rightmost bit, unpack using the noError encoding
        return pat::PackedCandidate::unpackTimeNoError(time/2);
    }
}
int16_t pat::PackedCandidate::packTimeWithError(float   time, float   timeError) {
    // Encode in units of timeError * 2^EXPO_TIME_WITHERROR (~1.6% if EXPO_TIME_WITHERROR = -6)
    // the largest value that can be stored in 14 bits + sign bit + overflow bit is about 260 sigmas
    // values larger than that will be stored using the no-timeError packing (with less precision).
    // overflows of these kinds should happen only for particles that are late arriving, out-of-time,
    // or mis-reconstructed, as timeError is O(20ps) and the beam spot witdth is O(200ps)
    float fpacked = std::round(time/std::ldexp(timeError, EXPO_TIME_WITHERROR));
    if (std::abs(fpacked) < 16383.f) { // 16383 = (2^14 - 1) = largest absolute value for a signed 15 bit integer
        return int16_t(fpacked) * 2; // make it even, and fit in a signed 16 bit int
    } else {
        int16_t packed = packTimeNoError(time); // encode
        return packed * 2 + (time > 0 ? +1 : -1); // make it odd, to signal that there was an overlow
    }
}