CovarianceParameterization.cc

Go to the documentation of this file.

#include "DataFormats/PatCandidates/interface/CovarianceParameterization.h"
#include "DataFormats/PatCandidates/interface/liblogintpack.h"
#include "DataFormats/PatCandidates/interface/libminifloat.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include <boost/format.hpp>
#include <iostream>
#include <TParameter.h>
#include <TVector.h>
#include <TFolder.h>

uint16_t CompressionElement::pack(float value, float ref) const
{
    float toCompress=0;
    switch(target) {
        case(realValue):
          toCompress=value;
          break;
        case(ratioToRef):
          toCompress=value/ref;
          break;
        case(differenceToRef):
          toCompress=value-ref;
          break;
    }
    switch(method) {
        case(float16):
          return MiniFloatConverter::float32to16(toCompress*params[0]);
          break;
        case(reduceMantissa):
          return MiniFloatConverter::reduceMantissaToNbits(toCompress,params[0]);
          break;
        case(zero):
          return 0;
          break;
        case(one):
          return 1.0;
          break;
        case(tanLogPack):
          return 0; //FIXME: should be implemented
          break;
        case(logPack):
          int16_t r=logintpack::pack16log(toCompress,params[0],params[1],bits);
          return * reinterpret_cast<uint16_t *>(&r); 
          break;
      
    }
  return 0;
}
float CompressionElement::unpack(uint16_t packed, float ref) const
{
    float unpacked=0;
    switch(method) {
        case(float16):
          unpacked= MiniFloatConverter::float16to32(packed)/params[0];
          break;
        case(reduceMantissa):
          unpacked=packed;
          break;
        case(logPack):
          unpacked=logintpack::unpack16log(* reinterpret_cast<int16_t *>(&packed),params[0],params[1],bits);
          break;
        case(zero):
          unpacked=0;
          break;
        case(one):
        case(tanLogPack):
          unpacked=1; //FIXME: should be implemented
    }
    switch(target) {
        case(realValue):
          return unpacked;
        case(ratioToRef):
          return unpacked*ref;
        case(differenceToRef):
          return unpacked+ref;
    }

    return ref;
  
}






void CovarianceParameterization::load(int version)
{
 edm::FileInPath fip((boost::format("DataFormats/PatCandidates/data/CovarianceParameterization_version%d.root") % version).str());
 fileToRead_ = TFile::Open(fip.fullPath().c_str()); 
 TFile & fileToRead=*fileToRead_;
//Read files from here fip.fullPath().c_str();
 if(fileToRead.IsOpen())  {
     readFile(fileToRead);

     TIter next(((TDirectoryFile*)fileToRead.Get("schemas"))->GetListOfKeys());
     TKey *key;
     while ((key = (TKey*)next())) {
      TClass *cl = gROOT->GetClass(key->GetClassName());
      if (!cl->InheritsFrom("TDirectoryFile")) continue;
      std::string schemaNumber = key->ReadObj()->GetName();
      uint16_t schemaN = std::stoi(schemaNumber);
     //for (int folderNumber = 0; folderNumber < 6 ; folderNumber++) {
      CompressionSchema schema; 
      for (int i = 0; i < 5; i++) {
            for (int j = i; j < 5; j++) {        //FILLING ONLY THE SCHEMA OF SOME ELEMENTS
                  std::string folder = "schemas/" + schemaNumber + "/"  + char(48+i) + char(48+j);
                std::string methodString = folder + "/method";
                std::string targetString = folder + "/target";
                std::string bitString = folder + "/bit";
                std::vector<float> vParams ;
                TVector *p=(TVector*) fileToRead.Get((folder+"/param").c_str());
                for(int k = 0 ; k < p->GetNoElements() ; k++){
                    vParams.push_back((*p)[k]);
                }

                schema(i,j)=CompressionElement((CompressionElement::Method) ((TParameter<int>*) fileToRead.Get(methodString.c_str()))->GetVal(),
        (CompressionElement::Target) ((TParameter<int>*) fileToRead.Get(targetString.c_str()))->GetVal(),
         (int) ((TParameter<int>*) fileToRead.Get(bitString.c_str()))->GetVal(), vParams);
        
                
            }
        }
     schemas[schemaN]=schema; 
     }

    loadedVersion_=version; 
 } else {loadedVersion_=-1;}

}



void CovarianceParameterization::readFile( TFile & f) {

    for (int i = 0; i < 5; i++) {
        for (int j = i; j < 5; j++) {

            std::string String_first_positive = "_pixel_";
            std::string String_second_positive = "_noPixel_";

            addTheHistogram(&cov_elements_pixelHit, String_first_positive, i, j,f);
            addTheHistogram(&cov_elements_noPixelHit, String_second_positive, i, j,f);


        }
    }
}


void CovarianceParameterization::addTheHistogram(std::vector<TH3D *> * HistoVector, std::string StringToAddInTheName, int i, int j,  TFile & fileToRead) {

    std::string  List_covName[5] = {"qoverp", "lambda", "phi", "dxy", "dsz"};

    std::string histoNameString = "covariance_" + List_covName[i] + "_" + List_covName[j] + StringToAddInTheName+"parametrization" ;// + "_entries";
    TH3D * matrixElememtHistogramm = (TH3D*) fileToRead.Get(histoNameString.c_str());
    HistoVector->push_back(matrixElememtHistogramm);
}



float CovarianceParameterization::meanValue(int i,int j,int sign,float pt, float eta, int nHits,int pixelHits,  float cii,float cjj) const {
    int hitNumberToUse = nHits;
    if (hitNumberToUse < 2 ) hitNumberToUse = 2;
    if (hitNumberToUse > 32 ) hitNumberToUse = 32;
    int  ptBin = cov_elements_pixelHit[0]->GetXaxis()->FindBin(pt);
    int etaBin = cov_elements_pixelHit[0]->GetYaxis()->FindBin(std::abs(eta));
    int hitBin = cov_elements_pixelHit[0]->GetZaxis()->FindBin(hitNumberToUse);
    int min_idx = i;
    int max_idx = j;

    if (i>j) {
        min_idx = j;
        max_idx = i;
    }

    int indexOfTheHitogramInTheList = ((9 - min_idx)*min_idx)/2 + max_idx;

    double meanValue = 0.;
    if (pixelHits > 0) {
            meanValue =sign* cov_elements_pixelHit[indexOfTheHitogramInTheList]->GetBinContent(ptBin, etaBin, hitBin);
        }
        else {
            meanValue = sign*cov_elements_noPixelHit[indexOfTheHitogramInTheList]->GetBinContent(ptBin, etaBin, hitBin);
        }
    return meanValue;

}

float CovarianceParameterization::pack(float value, int schema, int i,int j,float pt, float eta, int nHits,int pixelHits,  float cii,float cjj) const {
    if(i>j) std::swap(i,j);
    float ref=meanValue(i,j,1.,pt,eta,nHits,pixelHits,cii,cjj);
    if(ref==0) {
      schema=0;
    }
    if(schema==0 && i==j && (i==2 || i==0) ) ref=1./(pt*pt);
/*  //Used for debugging, to be later removed  
    uint16_t p=(*schemas.find(schema)).second(i,j).pack(value,ref);
    float up=(*schemas.find(schema)).second(i,j).unpack(p,ref);
    std::cout << "check " << i << " " << j << " " << value << " " << up << " " << p << " " << ref << " " << schema<< std::endl;*/
    return (*schemas.find(schema)).second(i,j).pack(value,ref);
}
float CovarianceParameterization::unpack(uint16_t packed, int schema, int i,int j,float pt, float eta, int nHits,int pixelHits,  float cii,float cjj) const {
    if(i>j) std::swap(i,j);
    float ref=meanValue(i,j,1.,pt,eta,nHits,pixelHits,cii,cjj);
    if(ref==0) {
      schema=0;
    }
    if(schema==0 && i==j && (i==2 || i==0) ) ref=1./(pt*pt);
    if(i==j && (*schemas.find(schema)).second(i,j).unpack(packed,ref)==0) return 1e-9;
    else
    return (*schemas.find(schema)).second(i,j).unpack(packed,ref);
 
}