/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/CalibCalorimetry/CaloTPG/src/CaloTPGTranscoderULUT.cc

Go to the documentation of this file.

#include "CalibCalorimetry/CaloTPG/src/CaloTPGTranscoderULUT.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <iostream>
#include <fstream>
#include <math.h>

//#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "CondFormats/DataRecord/interface/HcalLutMetadataRcd.h"

using namespace std;

HcalTrigTowerGeometry theTrigTowerGeometry;

CaloTPGTranscoderULUT::CaloTPGTranscoderULUT(const std::string& compressionFile,
                                             const std::string& decompressionFile)
                                                : isLoaded_(false), nominal_gain_(0.), rctlsb_factor_(0.),
                                                  compressionFile_(compressionFile),
                                                  decompressionFile_(decompressionFile)
{
  for (int i = 0; i < NOUTLUTS; i++) outputLUT_[i] = 0;
}

CaloTPGTranscoderULUT::~CaloTPGTranscoderULUT() {
  for (int i = 0; i < NOUTLUTS; i++) {
    if (outputLUT_[i] != 0) delete [] outputLUT_[i];
  }
}

void CaloTPGTranscoderULUT::loadHCALCompress() const{
// Initialize analytical compression LUT's here
   // TODO cms::Error log
  if (OUTPUT_LUT_SIZE != (unsigned int) 0x400) std::cout << "Error: Analytic compression expects 10-bit LUT; found LUT with " << OUTPUT_LUT_SIZE << " entries instead" << std::endl;

  std::vector<unsigned int> analyticalLUT(OUTPUT_LUT_SIZE, 0);
  std::vector<unsigned int> identityLUT(OUTPUT_LUT_SIZE, 0);

  // Compute compression LUT
  for (unsigned int i=0; i < OUTPUT_LUT_SIZE; i++) {
        analyticalLUT[i] = (unsigned int)(sqrt(14.94*log(1.+i/14.94)*i) + 0.5);
        identityLUT[i] = min(i,0xffu);
  }
 
  for (int ieta=-32; ieta <= 32; ieta++){
     for (int iphi = 1; iphi <= 72; iphi++){
        if (!HTvalid(ieta,iphi)) continue;
        int lutId = getOutputLUTId(ieta,iphi);
        // TODO cms::Error log
        if (outputLUT_[lutId] != 0){
           std::cout << "Error: LUT with (ieta,iphi) = (" << ieta << "," << iphi << ") has been previously allocated!" << std::endl;
           continue;
        }

        outputLUT_[lutId] = new LUT[OUTPUT_LUT_SIZE];

        HcalTrigTowerDetId id(ieta, iphi);
        const HcalLutMetadatum *meta = lutMetadata_->getValues(id);
        int threshold = meta->getOutputLutThreshold();

        for (int i = 0; i < threshold; ++i)
           outputLUT_[lutId][i] = 0;

        for (unsigned int i = threshold; i < OUTPUT_LUT_SIZE; ++i)
           outputLUT_[lutId][i] = (abs(ieta) < theTrigTowerGeometry.firstHFTower()) ? analyticalLUT[i] : identityLUT[i];
     } //for iphi
  } //for ieta
}

void CaloTPGTranscoderULUT::loadHCALCompress(const std::string& filename) const{
  int tool;
  std::ifstream userfile;
  std::vector< std::vector<LUT> > outputluts;

  std::cout << "Initializing compression LUT's from " << (char *)filename.data() << std::endl;
  for (int i = 0; i < NOUTLUTS; i++) outputLUT_[i] = 0;
  int maxid = 0, minid = 0x7FFFFFFF, rawid = 0;
  for (int ieta=-32; ieta <= 32; ieta++) {
    for (int iphi = 1; iphi <= 72; iphi++) {
                if (HTvalid(ieta,iphi)) {
          rawid = getOutputLUTId(ieta, iphi);
          if (outputLUT_[rawid] != 0) std::cout << "Error: LUT with (ieta,iphi) = (" << ieta << "," << iphi << ") has been previously allocated!" << std::endl;
          else outputLUT_[rawid] = new LUT[OUTPUT_LUT_SIZE];
          if (rawid < minid) minid = rawid;
          if (rawid > maxid) maxid = rawid;
        }
        }
  }

  userfile.open((char *)filename.data());

  if( userfile ) {
    int nluts = 0;
        std::string s;
    std::vector<int> loieta,hiieta;
    std::vector<int> loiphi,hiiphi;
    getline(userfile,s);

    getline(userfile,s);

        unsigned int index = 0;
        while (index < s.length()) {
          while (isspace(s[index])) index++;
          if (index < s.length()) nluts++;
          while (!isspace(s[index])) index++;
        }
        for (unsigned int i=0; i<=s.length(); i++) userfile.unget(); //rewind last line
    outputluts.resize(nluts);
    for (int i=0; i<nluts; i++) outputluts[i].resize(OUTPUT_LUT_SIZE);
    
    
    for (int i=0; i<nluts; i++) {
      userfile >> tool;
      loieta.push_back(tool);
    }
    
    for (int i=0; i<nluts; i++) {
      userfile >> tool;
      hiieta.push_back(tool);    
    }
   
    for (int i=0; i<nluts; i++) {
      userfile >> tool;
      loiphi.push_back(tool);
        
    }
    
    for (int i=0; i<nluts; i++) {
      userfile >> tool;
      hiiphi.push_back(tool);
    
    }
        
    for (unsigned int j=0; j<OUTPUT_LUT_SIZE; j++) { 
      for(int i=0; i <nluts; i++) {
                userfile >> tool;
                if (tool < 0) {
                        std::cout << "Error: LUT can't have negative numbers; 0 used instead: " << i << ", " << j << " : = " << tool << std::endl;
                        tool = 0;
                } else if (tool > 0xff) {
                        std::cout << "Error: LUT can't have >8-bit numbers; 0xff used instead: " << i << ", " << j << " : = " << tool << std::endl;
                        tool = 0xff;
                }
                outputluts[i][j] = tool;
                if (userfile.eof()) std::cout << "Error: LUT file is truncated or has a wrong format: " << i << "," << j << std::endl;
          }
    }
    userfile.close();
        
        HcalDetId cell;
        int id, ntot = 0;
        for (int i=0; i < nluts; i++) {
                int nini = 0;
                for (int iphi = loiphi[i]; iphi <= hiiphi[i]; iphi++) {      
                        for (int ieta=loieta[i]; ieta <= hiieta[i]; ieta++) {
                                        if (HTvalid(ieta,iphi)) {  
                                                id = getOutputLUTId(ieta,iphi);
                                                if (outputLUT_[id] == 0) throw cms::Exception("PROBLEM: outputLUT_ has not been initialized for ieta, iphi, id = ") << ieta << ", " << iphi << ", " << id << std::endl;
                                        for (int j = 0; j <= 0x3ff; j++) outputLUT_[id][j] = outputluts[i][j];
                                                nini++;
                                                ntot++;
                                        }
                        }
       }
        
    }
        
  } else {
    
        loadHCALCompress();
  }
}

void CaloTPGTranscoderULUT::loadHCALUncompress() const {
   hcaluncomp_.clear();
   for (int i = 0; i < NOUTLUTS; i++){
      RCTdecompression decompressionTable(TPGMAX,0);
      hcaluncomp_.push_back(decompressionTable);
   }

   for (int ieta = -32; ieta <= 32; ++ieta){

      double eta_low = 0., eta_high = 0.;
                theTrigTowerGeometry.towerEtaBounds(ieta,eta_low,eta_high); 
      double cosh_ieta = fabs(cosh((eta_low + eta_high)/2.));

                for (int iphi = 1; iphi <= 72; iphi++) {
         if (!HTvalid(ieta, iphi)) continue;

                        int lutId = getOutputLUTId(ieta,iphi);
         HcalTrigTowerDetId id(ieta, iphi);

         const HcalLutMetadatum *meta = lutMetadata_->getValues(id);
         double factor = 0.;

         // HF
         if (abs(ieta) >= theTrigTowerGeometry.firstHFTower())
            factor = rctlsb_factor_;
         // HBHE
         else 
            factor = nominal_gain_ / cosh_ieta * meta->getLutGranularity();

         // tpg - compressed value
         unsigned int tpg = outputLUT_[lutId][0];
         int low = 0;
         for (unsigned int i = 0; i < OUTPUT_LUT_SIZE; ++i){
            if (outputLUT_[lutId][i] != tpg){
               unsigned int mid = (low + i)/2;
               hcaluncomp_[lutId][tpg] = (tpg == 0 ? low : factor * mid);
               low = i;
               tpg = outputLUT_[lutId][i];
            }
         }
         hcaluncomp_[lutId][tpg] = factor * low;
                } // for phi
        } // for ieta
}

void CaloTPGTranscoderULUT::loadHCALUncompress(const std::string& filename) const {
  std::ifstream userfile;
  userfile.open(filename.c_str());
  
  hcaluncomp_.resize(NOUTLUTS);
  for (int i = 0; i < NOUTLUTS; i++) hcaluncomp_[i].resize(TPGMAX);
  
  static const int etabound = 32;
  if( userfile ) {
         double et;
     for (int i=0; i<TPGMAX; i++) { 
      for(int j = 1; j <=etabound; j++) {
            userfile >> et;
                for (int iphi = 1; iphi <= 72; iphi++) {
                  int itower = getOutputLUTId(j,iphi);
                  if (itower >= 0) hcaluncomp_[itower][i] = et;
                  itower = getOutputLUTId(-j,iphi);
                  if (itower >= 0) hcaluncomp_[itower][i] = et;           
                }
          }
         }
         userfile.close();
  }
  else {

        loadHCALUncompress();
  }
}

HcalTriggerPrimitiveSample CaloTPGTranscoderULUT::hcalCompress(const HcalTrigTowerDetId& id, unsigned int sample, bool fineGrain) const {
  int ieta = id.ieta();
  int iphi = id.iphi();
//  if (abs(ieta) > 28) iphi = iphi/4 + 1; // Changing iphi index from 1, 5, ..., 69 to 1, 2, ..., 18
  int itower = getOutputLUTId(ieta,iphi);

  if (itower < 0) cms::Exception("Invalid Data") << "No trigger tower found for ieta, iphi = " << ieta << ", " << iphi;
  if (sample >= OUTPUT_LUT_SIZE) {

    throw cms::Exception("Out of Range") << "LUT has 1024 entries for " << itower << " but " << sample << " was requested.";
    sample=OUTPUT_LUT_SIZE - 1;
  }

  return HcalTriggerPrimitiveSample(outputLUT_[itower][sample],fineGrain,0,0);
}

double CaloTPGTranscoderULUT::hcaletValue(const int& ieta, const int& iphi, const int& compET) const {
  if (hcaluncomp_.empty()) {

        CaloTPGTranscoderULUT::loadHCALUncompress(decompressionFile_);
  }
  double etvalue = 0.;
  int itower = getOutputLUTId(ieta,iphi);
  if (itower < 0) std::cout << "hcaletValue error: no decompression LUT found for ieta, iphi = " << ieta << ", " << iphi << std::endl;
  else if (compET < 0 || compET > 0xff) std::cout << "hcaletValue error: compressed value out of range: eta, phi, cET = " << ieta << ", " << iphi << ", " << compET << std::endl;
  else etvalue = hcaluncomp_[itower][compET];
  return(etvalue);
}

double CaloTPGTranscoderULUT::hcaletValue(const int& ieta, const int& compET) const {
// This is now an obsolete method; we return the AVERAGE over all the allowed iphi channels if it's invoked
// The user is encouraged to use hcaletValue(const int& ieta, const int& iphi, const int& compET) instead

  if (hcaluncomp_.empty()) {
        std::cout << "Initializing the RCT decompression table from the file: " << decompressionFile_ << std::endl;
        CaloTPGTranscoderULUT::loadHCALUncompress(decompressionFile_);
  }

  double etvalue = 0.;
  if (compET < 0 || compET > 0xff) std::cout << "hcaletValue error: compressed value out of range: eta, cET = " << ieta << ", " << compET << std::endl;
  else {
        int nphi = 0;
        for (int iphi=1; iphi <= 72; iphi++) {
                if (HTvalid(ieta,iphi)) {
                        nphi++;
                        int itower = getOutputLUTId(ieta,iphi);
                        etvalue += hcaluncomp_[itower][compET];
                }
        }
        if (nphi > 0) etvalue /= nphi;
        else std::cout << "hcaletValue error: no decompression LUTs found for any iphi for ieta = " << ieta << std::endl;
  }
  return(etvalue);
}

double CaloTPGTranscoderULUT::hcaletValue(const HcalTrigTowerDetId& hid, const HcalTriggerPrimitiveSample& hc) const {
  if (hcaluncomp_.empty()) loadHCALUncompress(decompressionFile_);

  int ieta = hid.ieta();                        // No need to check the validity,
  int iphi = hid.iphi();                        // as the values are guaranteed
  int compET = hc.compressedEt();       // to be within the range by the class
  int itower = getOutputLUTId(ieta,iphi);
  double etvalue = hcaluncomp_[itower][compET];
  return(etvalue);
}

EcalTriggerPrimitiveSample CaloTPGTranscoderULUT::ecalCompress(const EcalTrigTowerDetId& id, unsigned int sample, bool fineGrain) const {
  throw cms::Exception("Not Implemented") << "CaloTPGTranscoderULUT::ecalCompress";
}

void CaloTPGTranscoderULUT::rctEGammaUncompress(const HcalTrigTowerDetId& hid, const HcalTriggerPrimitiveSample& hc,
                                                const EcalTrigTowerDetId& eid, const EcalTriggerPrimitiveSample& ec, 
                                                unsigned int& et, bool& egVecto, bool& activity) const {
  throw cms::Exception("Not Implemented") << "CaloTPGTranscoderULUT::rctEGammaUncompress";
}
void CaloTPGTranscoderULUT::rctJetUncompress(const HcalTrigTowerDetId& hid, const HcalTriggerPrimitiveSample& hc,
                                             const EcalTrigTowerDetId& eid, const EcalTriggerPrimitiveSample& ec, 
                                             unsigned int& et) const {
  throw cms::Exception("Not Implemented") << "CaloTPGTranscoderULUT::rctJetUncompress";
}

bool CaloTPGTranscoderULUT::HTvalid(const int ieta, const int iphiin) const {
        int iphi = iphiin;
        if (iphi <= 0 || iphi > 72 || ieta == 0 || abs(ieta) > 32) return false;
        if (abs(ieta) > 28) {
          if (newHFphi) {
            if ((iphi/4)*4 + 1 != iphi) return false;
            iphi = iphi/4 + 1;
          }
          if (iphi > 18) return false;
        }
        return true;
}

int CaloTPGTranscoderULUT::getOutputLUTId(const int ieta, const int iphiin) const {
        int iphi = iphiin;
        if (HTvalid(ieta, iphi)) {
                int offset = 0, ietaabs;
                ietaabs = abs(ieta);
                if (ieta < 0) offset = NOUTLUTS/2;
                if (ietaabs < 29) return 72*(ietaabs - 1) + (iphi - 1) + offset;
                else {
                  if (newHFphi) iphi = iphi/4 + 1;
                  return 18*(ietaabs - 29) + iphi + 2015 + offset;
                }
        } else return -1;       
}

std::vector<unsigned char> CaloTPGTranscoderULUT::getCompressionLUT(HcalTrigTowerDetId id) const {
   std::vector<unsigned char> lut;
   int itower = getOutputLUTId(id.ieta(),id.iphi());
   if (itower >= 0) {
         lut.resize(OUTPUT_LUT_SIZE);
         for (unsigned int i = 0; i < OUTPUT_LUT_SIZE; i++) lut[i]=outputLUT_[itower][i];
   } 
   return lut;
}

void CaloTPGTranscoderULUT::setup(const edm::EventSetup& es, Mode mode=All) const{
   if (isLoaded_) return;
   // TODO Try/except
   es.get<HcalLutMetadataRcd>().get(lutMetadata_);

   nominal_gain_ = lutMetadata_->getNominalGain();
   float rctlsb =lutMetadata_->getRctLsb();
   if (rctlsb != 0.25 && rctlsb != 0.5)
      throw cms::Exception("RCTLSB") << " value=" << rctlsb << " (should be 0.25 or 0.5)" << std::endl;
   rctlsb_factor_ = rctlsb;

   if (compressionFile_.empty() && decompressionFile_.empty()) {
      loadHCALCompress();
   }
   else {
      // TODO Message to discourage using txt.
      std::cout << "From Text File:" << std::endl;
      loadHCALCompress(compressionFile_);
   }
   isLoaded_ = true;
}

void CaloTPGTranscoderULUT::printDecompression() const{
   std::cout << "RCT Decompression table" << std::endl;                          
   for (int i=0; i < 256; i++) {                                                 
      for (int j=1; j <= theTrigTowerGeometry.nTowers(); j++)
         std::cout << int(hcaletValue(j,i)*100. + 0.5)/100. << " ";
      std::cout << std::endl;                                               
      for (int j=1; j <= theTrigTowerGeometry.nTowers(); j++)               
         if (hcaletValue(j,i) != hcaletValue(-j,i))                    
            cout << "Error: decompression table for ieta = +/- " << j << " disagree! " << hcaletValue(-j,i) << ", " << hcaletValue(j,i) << endl;        
   }
}

//int CaloTPGTranscoderULUT::getLutGranularity(const DetId& id) const{
//   int ieta = id.ietaAbs();
//   if (ieta < 18) return 1;
//   else if (ieta < 27) return 2;
//   else if (ieta < 29) return 5;
//   return 0;
//}
//
//int CaloTPGTranscoderULUT::getLutThreshold(const DetId& id) const{
//   int ieta = id.ietaAbs();
//   if (ieta < 18) return 4;
//   else if (ieta < 27) return 2;
//   else if (ieta < 29) return 1;
//   return 0;
//}
//