PatCandidates/src/CovarianceParameterization.cc

0001 //#include <iostream>
0002
0003 #include <fmt/printf.h>
0004
0005 #include <TParameter.h>
0006 #include <TVector.h>
0007 #include <TFolder.h>
0008
0009 #include "DataFormats/Math/interface/liblogintpack.h"
0010 #include "DataFormats/Math/interface/libminifloat.h"
0011 #include "DataFormats/PatCandidates/interface/CovarianceParameterization.h"
0012 #include "FWCore/Utilities/interface/FileInPath.h"
0013
0014 uint16_t CompressionElement::pack(float value, float ref) const {
0015   float toCompress = 0;
0016   switch (target) {
0017     case (realValue):
0018       toCompress = value;
0019       break;
0020     case (ratioToRef):
0021       toCompress = value / ref;
0022       break;
0023     case (differenceToRef):
0024       toCompress = value - ref;
0025       break;
0026   }
0027   switch (method) {
0028     case (float16):
0029       return MiniFloatConverter::float32to16(toCompress * params[0]);
0030       break;
0031     case (reduceMantissa):
0032       return MiniFloatConverter::reduceMantissaToNbits(toCompress, params[0]);
0033       break;
0034     case (zero):
0035       return 0;
0036       break;
0037     case (one):
0038       return 1.0;
0039       break;
0040     case (tanLogPack):
0041       return 0;  //FIXME: should be implemented
0042       break;
0043     case (logPack):
0044       int16_t r = logintpack::pack16log(toCompress, params[0], params[1], bits);
0045       return *reinterpret_cast<uint16_t *>(&r);
0046       break;
0047   }
0048   return 0;
0049 }
0050 float CompressionElement::unpack(uint16_t packed, float ref) const {
0051   float unpacked = 0;
0052   switch (method) {
0053     case (float16):
0054       unpacked = MiniFloatConverter::float16to32(packed) / params[0];
0055       break;
0056     case (reduceMantissa):
0057       unpacked = packed;
0058       break;
0059     case (logPack):
0060       unpacked = logintpack::unpack16log(*reinterpret_cast<int16_t *>(&packed), params[0], params[1], bits);
0061       break;
0062     case (zero):
0063       unpacked = 0;
0064       break;
0065     case (one):
0066     case (tanLogPack):
0067       unpacked = 1;  //FIXME: should be implemented
0068   }
0069   switch (target) {
0070     case (realValue):
0071       return unpacked;
0072     case (ratioToRef):
0073       return unpacked * ref;
0074     case (differenceToRef):
0075       return unpacked + ref;
0076   }
0077
0078   return ref;
0079 }
0080
0081 void CovarianceParameterization::load(int version) {
0082   edm::FileInPath fip(
0083       fmt::sprintf("DataFormats/PatCandidates/data/CovarianceParameterization_version%d.root", version));
0084   fileToRead_ = TFile::Open(fip.fullPath().c_str());
0085   TFile &fileToRead = *fileToRead_;
0086   //Read files from here fip.fullPath().c_str();
0087   if (fileToRead.IsOpen()) {
0088     readFile(fileToRead);
0089
0090     TIter next(((TDirectoryFile *)fileToRead.Get("schemas"))->GetListOfKeys());
0091     TKey *key;
0092     while ((key = (TKey *)next())) {
0093       TClass *cl = gROOT->GetClass(key->GetClassName());
0094       if (!cl->InheritsFrom("TDirectoryFile"))
0095         continue;
0096       std::string schemaNumber = key->ReadObj()->GetName();
0097       uint16_t schemaN = std::stoi(schemaNumber);
0098       //for (int folderNumber = 0; folderNumber < 6 ; folderNumber++) {
0099       CompressionSchema schema;
0100       for (int i = 0; i < 5; i++) {
0101         for (int j = i; j < 5; j++) {  //FILLING ONLY THE SCHEMA OF SOME ELEMENTS
0102           std::string folder = "schemas/" + schemaNumber + "/" + char(48 + i) + char(48 + j);
0103           std::string methodString = folder + "/method";
0104           std::string targetString = folder + "/target";
0105           std::string bitString = folder + "/bit";
0106           std::vector<float> vParams;
0107           TVector *p = (TVector *)fileToRead.Get((folder + "/param").c_str());
0108           vParams.reserve(p->GetNoElements());
0109           for (int k = 0; k < p->GetNoElements(); k++) {
0110             vParams.push_back((*p)[k]);
0111           }
0112
0113           schema(i, j) = CompressionElement(
0114               (CompressionElement::Method)((TParameter<int> *)fileToRead.Get(methodString.c_str()))->GetVal(),
0115               (CompressionElement::Target)((TParameter<int> *)fileToRead.Get(targetString.c_str()))->GetVal(),
0116               (int)((TParameter<int> *)fileToRead.Get(bitString.c_str()))->GetVal(),
0117               vParams);
0118         }
0119       }
0120       schemas[schemaN] = schema;
0121     }
0122
0123     loadedVersion_ = version;
0124   } else {
0125     loadedVersion_ = -1;
0126   }
0127 }
0128
0129 void CovarianceParameterization::readFile(TFile &f) {
0130   for (int i = 0; i < 5; i++) {
0131     for (int j = i; j < 5; j++) {
0132       std::string String_first_positive = "_pixel_";
0133       std::string String_second_positive = "_noPixel_";
0134
0135       addTheHistogram(&cov_elements_pixelHit, String_first_positive, i, j, f);
0136       addTheHistogram(&cov_elements_noPixelHit, String_second_positive, i, j, f);
0137     }
0138   }
0139 }
0140
0141 void CovarianceParameterization::addTheHistogram(
0142     std::vector<TH3D *> *HistoVector, std::string StringToAddInTheName, int i, int j, TFile &fileToRead) {
0143   std::string List_covName[5] = {"qoverp", "lambda", "phi", "dxy", "dsz"};
0144
0145   std::string histoNameString = "covariance_" + List_covName[i] + "_" + List_covName[j] + StringToAddInTheName +
0146                                 "parametrization";  // + "_entries";
0147   TH3D *matrixElememtHistogramm = (TH3D *)fileToRead.Get(histoNameString.c_str());
0148   HistoVector->push_back(matrixElememtHistogramm);
0149 }
0150
0151 float CovarianceParameterization::meanValue(
0152     int i, int j, int sign, float pt, float eta, int nHits, int pixelHits, float cii, float cjj) const {
0153   int hitNumberToUse = nHits;
0154   if (hitNumberToUse < 2)
0155     hitNumberToUse = 2;
0156   if (hitNumberToUse > 32)
0157     hitNumberToUse = 32;
0158   int ptBin = cov_elements_pixelHit[0]->GetXaxis()->FindBin(pt);
0159   int etaBin = cov_elements_pixelHit[0]->GetYaxis()->FindBin(std::abs(eta));
0160   int hitBin = cov_elements_pixelHit[0]->GetZaxis()->FindBin(hitNumberToUse);
0161   int min_idx = i;
0162   int max_idx = j;
0163
0164   if (i > j) {
0165     min_idx = j;
0166     max_idx = i;
0167   }
0168
0169   int indexOfTheHitogramInTheList = ((9 - min_idx) * min_idx) / 2 + max_idx;
0170
0171   double meanValue = 0.;
0172   if (pixelHits > 0) {
0173     meanValue = sign * cov_elements_pixelHit[indexOfTheHitogramInTheList]->GetBinContent(ptBin, etaBin, hitBin);
0174   } else {
0175     meanValue = sign * cov_elements_noPixelHit[indexOfTheHitogramInTheList]->GetBinContent(ptBin, etaBin, hitBin);
0176   }
0177   return meanValue;
0178 }
0179
0180 float CovarianceParameterization::pack(
0181     float value, int schema, int i, int j, float pt, float eta, int nHits, int pixelHits, float cii, float cjj) const {
0182   if (i > j)
0183     std::swap(i, j);
0184   float ref = meanValue(i, j, 1., pt, eta, nHits, pixelHits, cii, cjj);
0185   if (ref == 0) {
0186     schema = 0;
0187   }
0188   if (schema == 0 && i == j && (i == 2 || i == 0))
0189     ref = 1. / (pt * pt);
0190   /*  //Used for debugging, to be later removed
0191     uint16_t p=(*schemas.find(schema)).second(i,j).pack(value,ref);
0192     float up=(*schemas.find(schema)).second(i,j).unpack(p,ref);
0193     std::cout << "check " << i << " " << j << " " << value << " " << up << " " << p << " " << ref << " " << schema<< std::endl;*/
0194   return (*schemas.find(schema)).second(i, j).pack(value, ref);
0195 }
0196 float CovarianceParameterization::unpack(
0197     uint16_t packed, int schema, int i, int j, float pt, float eta, int nHits, int pixelHits, float cii, float cjj)
0198     const {
0199   if (i > j)
0200     std::swap(i, j);
0201   float ref = meanValue(i, j, 1., pt, eta, nHits, pixelHits, cii, cjj);
0202   if (ref == 0) {
0203     schema = 0;
0204   }
0205   if (schema == 0 && i == j && (i == 2 || i == 0))
0206     ref = 1. / (pt * pt);
0207   if (i == j && (*schemas.find(schema)).second(i, j).unpack(packed, ref) == 0)
0208     return 1e-9;
0209   else
0210     return (*schemas.find(schema)).second(i, j).unpack(packed, ref);
0211 }