Warning, /CalibCalorimetry/HcalTPGAlgos/test/FixedPed3_UlutCoder_ is written in an unsupported language. File is not indexed.
0001 #include "CalibCalorimetry/HcalTPGAlgos/interface/HcaluLUTTPGCoder.h"
0002 #include "Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryData.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include <iostream>
0007 #include <fstream>
0008 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0009 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0010 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0011 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0012 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0013 #include <cmath>
0014
0015 const float HcaluLUTTPGCoder::nominal_gain = 0.177;
0016
0017 HcaluLUTTPGCoder::HcaluLUTTPGCoder(const char* filename) {
0018 AllocateLUTs();
0019 getRecHitCalib(filename);
0020 }
0021
0022 HcaluLUTTPGCoder::HcaluLUTTPGCoder(const char* filename, const char* fname2) {
0023 throw cms::Exception("PROBLEM: This constructor should never be invoked!");
0024 }
0025
0026 void HcaluLUTTPGCoder::compress(const IntegerCaloSamples& ics, const std::vector<bool>& featureBits, HcalTriggerPrimitiveDigi& tp) const {
0027 throw cms::Exception("PROBLEM: This method should never be invoked!");
0028 }
0029
0030 HcaluLUTTPGCoder::~HcaluLUTTPGCoder() {
0031 for (int i = 0; i < nluts; i++) {
0032 if (inputLUT[i] != 0) delete [] inputLUT[i];
0033 }
0034 delete [] _gain;
0035 delete [] _ped;
0036 }
0037
0038 void HcaluLUTTPGCoder::AllocateLUTs() {
0039 HcalTopology theTopo;
0040 HcalDetId did;
0041
0042 _ped = new float[nluts];
0043 _gain = new float[nluts];
0044 for (int i = 0; i < nluts; i++) inputLUT[i] = 0;
0045 int maxid = 0, minid = 0x7FFFFFFF, rawid = 0;
0046 for (int ieta=-41; ieta <= 41; ieta++) {
0047 for (int iphi = 1; iphi <= 72; iphi++) {
0048 for (int depth = 1; depth <= 3; depth++) {
0049 did=HcalDetId(HcalBarrel,ieta,iphi,depth);
0050 if (theTopo.valid(did)) {
0051 rawid = GetLUTID(HcalBarrel, ieta, iphi, depth);
0052 if (inputLUT[rawid] != 0) std::cout << "Error: LUT with (ieta,iphi,depth) = (" << ieta << "," << iphi << "," << depth << ") has been previously allocated!" << std::endl;
0053 else inputLUT[rawid] = new LUT[INPUT_LUT_SIZE];
0054 if (rawid < minid) minid = rawid;
0055 if (rawid > maxid) maxid = rawid;
0056 }
0057
0058 did=HcalDetId(HcalEndcap,ieta,iphi,depth);
0059 if (theTopo.valid(did)) {
0060 rawid = GetLUTID(HcalEndcap, ieta, iphi, depth);
0061 if (inputLUT[rawid] != 0) std::cout << "Error: LUT with (ieta,iphi,depth) = (" << ieta << "," << iphi << "," << depth << ") has been previously allocated!" << std::endl;
0062 else inputLUT[rawid] = new LUT[INPUT_LUT_SIZE];
0063 if (rawid < minid) minid = rawid;
0064 if (rawid > maxid) maxid = rawid;
0065 }
0066 did=HcalDetId(HcalForward,ieta,iphi,depth);
0067 if (theTopo.valid(did)) {
0068 rawid = GetLUTID(HcalForward, ieta, iphi, depth);
0069 if (inputLUT[rawid] != 0) std::cout << "Error: LUT with (ieta,iphi,depth) = (" << ieta << "," << iphi << "," << depth << ") has been previously allocated!" << std::endl;
0070 else inputLUT[rawid] = new LUT[INPUT_LUT_SIZE];
0071 if (rawid < minid) minid = rawid;
0072 if (rawid > maxid) maxid = rawid;
0073 }
0074 }
0075 }
0076 }
0077
0078 }
0079
0080 int HcaluLUTTPGCoder::GetLUTID(HcalSubdetector id, int ieta, int iphi, int depth) const {
0081 int detid = 0;
0082 if (id == HcalEndcap) detid = 1;
0083 else if (id == HcalForward) detid = 2;
0084 return iphi + 72 * ((ieta + 41) + 83 * (depth + 3 * detid)) - 7777;
0085 }
0086
0087 void HcaluLUTTPGCoder::getRecHitCalib(const char* filename) {
0088
0089 std::ifstream userfile;
0090 userfile.open(filename);
0091 int tool;
0092 float Rec_calib_[87];
0093
0094 if (userfile) {
0095 userfile >> tool;
0096
0097 if (tool != 86) {
0098 std::cout << "Wrong RecHit calibration filesize: " << tool << " (expect 86)" << std::endl;
0099 }
0100 for (int j=1; j<87; j++) {
0101 userfile >> Rec_calib_[j]; // Read the Calib factors
0102 Rcalib[j] = Rec_calib_[j] ;
0103 }
0104
0105 userfile.close();
0106 }
0107 else std::cout << "File " << filename << " with RecHit calibration factors not found" << std::endl;
0108 }
0109
0110 void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {
0111 const HcalQIEShape* shape = conditions.getHcalShape();
0112 HcalCalibrations calibrations;
0113 int id;
0114 float divide;
0115 HcalTopology theTopo;
0116
0117 float cosheta_[41], lsb_ = 1./16.;
0118 for (int i = 0; i < 13; i++) cosheta_[i+29] = cosh((theHFEtaBounds[i+1] + theHFEtaBounds[i])/2.);
0119
0120 for (int depth = 1; depth <= 3; depth++) {
0121 for (int iphi = 1; iphi <= 72; iphi++) {
0122 divide = 1.*nominal_gain;
0123 for (int ieta=-16; ieta <= 16; ieta++) {
0124 HcalDetId cell(HcalBarrel,ieta,iphi,depth);
0125 if (theTopo.valid(cell)) {
0126 id = GetLUTID(HcalBarrel,ieta,iphi,depth);
0127 if (inputLUT[id] == 0) throw cms::Exception("PROBLEM: inputLUT has not been initialized for HB, ieta, iphi, depth, id = ") << ieta << "," << iphi << "," << depth << "," << id << std::endl;
0128 conditions.makeHcalCalibration (cell, &calibrations);
0129 const HcalQIECoder* channelCoder = conditions.getHcalCoder (cell);
0130 HcalCoderDb coder (*channelCoder, *shape);
0131 //float ped_ = (calibrations.pedestal(0)+calibrations.pedestal(1)+calibrations.pedestal(2)+calibrations.pedestal(3))/4;
0132
0133 float gain_= (calibrations.gain(0)+calibrations.gain(1)+calibrations.gain(2)+calibrations.gain(3))/4;
0134 float ped_ = 0;
0135 HBHEDataFrame frame(cell);
0136 frame.setSize(1);
0137 CaloSamples samples(cell, 1);
0138 _ped[id] = ped_;
0139 _gain[id] = gain_;
0140 for (int j = 0; j <= 0x7F; j++) {
0141 HcalQIESample adc(j);
0142 frame.setSample(0,adc);
0143 coder.adc2fC(frame,samples);
0144 float adc2fC_ = samples[0];
0145 if (ieta <0 )inputLUT[id][j] = (LUT) std::min(std::max(0,int(((adc2fC_*gain_)-ped_)*Rcalib[abs(ieta)]/divide)-3), 0x3FF);
0146 else inputLUT[id][j] = (LUT) std::min(std::max(0,int(((adc2fC_*gain_)-ped_)*Rcalib[abs(ieta)+43]/divide)-3), 0x3FF);
0147 //if (ieta <0 )inputLUT[id][j] = (LUT) std::min(std::max(0,int((adc2fC_ - ped_)*gain_*Rcalib[abs(ieta)]/divide)), 0x3FF);
0148 //else inputLUT[id][j] = (LUT) std::min(std::max(0,int((adc2fC_ - ped_)*gain_*Rcalib[abs(ieta)+43]/divide)), 0x3FF);
0149 }
0150 }
0151 }
0152 for (int ieta=-29; ieta <= 29; ieta++) {
0153 HcalDetId cell(HcalEndcap,ieta,iphi,depth);
0154 if (theTopo.valid(cell)) {
0155 if (abs(ieta) < 21) divide = 1.*nominal_gain;
0156 else if (abs(ieta) < 27) divide = 2.*nominal_gain;
0157 else divide = 5.*nominal_gain;
0158 id = GetLUTID(HcalEndcap,ieta,iphi,depth);
0159 if (inputLUT[id] == 0) throw cms::Exception("PROBLEM: inputLUT has not been initialized for HE, ieta, iphi, depth, id = ") << ieta << "," << iphi << "," << depth << "," << id << std::endl;
0160 conditions.makeHcalCalibration (cell, &calibrations);
0161 const HcalQIECoder* channelCoder = conditions.getHcalCoder (cell);
0162 HcalCoderDb coder (*channelCoder, *shape);
0163 float ped_ = 0;//(calibrations.pedestal(0)+calibrations.pedestal(1)+calibrations.pedestal(2)+calibrations.pedestal(3))/4;
0164 float gain_= (calibrations.gain(0)+calibrations.gain(1)+calibrations.gain(2)+calibrations.gain(3))/4;
0165 HBHEDataFrame frame(cell);
0166 frame.setSize(1);
0167 CaloSamples samples(cell, 1);
0168 _ped[id] = ped_;
0169 _gain[id] = gain_;
0170 for (int j = 0; j <= 0x7F; j++) {
0171 HcalQIESample adc(j);
0172 frame.setSample(0,adc);
0173 coder.adc2fC(frame,samples);
0174 float adc2fC_ = samples[0];
0175 if ( ieta < 0 ) inputLUT[id][j] = (LUT) std::min(std::max(0,int(((adc2fC_*gain_)-ped_)*Rcalib[abs(ieta)+1]/divide)-3), 0x3FF);
0176 else inputLUT[id][j] = (LUT) std::min(std::max(0,int(((adc2fC_*gain_)-ped_)*Rcalib[abs(ieta)+44]/divide)-3), 0x3FF);
0177 }
0178 }
0179 }
0180 for (int ieta=-41; ieta <= 41; ieta++) {
0181 HcalDetId cell(HcalForward,ieta,iphi,depth);
0182 if (theTopo.valid(cell)) {
0183 id = GetLUTID(HcalForward,ieta,iphi,depth);
0184 if (inputLUT[id] == 0) throw cms::Exception("PROBLEM: inputLUT has not been initialized for HF, ieta, iphi, depth, id = ") << ieta << "," << iphi << "," << depth << "," << id << std::endl;
0185 conditions.makeHcalCalibration (cell, &calibrations);
0186 const HcalQIECoder* channelCoder = conditions.getHcalCoder (cell);
0187 HcalCoderDb coder (*channelCoder, *shape);
0188 float ped_ = 0;//(calibrations.pedestal(0)+calibrations.pedestal(1)+calibrations.pedestal(2)+calibrations.pedestal(3))/4;
0189 float gain_= (calibrations.gain(0)+calibrations.gain(1)+calibrations.gain(2)+calibrations.gain(3))/4;
0190
0191
0192 HFDataFrame frame(cell);
0193 frame.setSize(1);
0194 CaloSamples samples(cell, 1);
0195 _ped[id] = ped_;
0196 _gain[id] = gain_;
0197 for (int j = 0; j <= 0x7F; j++) {
0198 HcalQIESample adc(j);
0199 frame.setSample(0,adc);
0200 coder.adc2fC(frame,samples);
0201 float adc2fC_ = samples[0];
0202 if (ieta < 0 ) inputLUT[id][j] = (LUT) std::min(std::max(0,int(((adc2fC_*gain_)-ped_)*Rcalib[abs(ieta)+2]*gain_/lsb_/cosheta_[abs(ieta)])-3), 0x3FF);
0203 else inputLUT[id][j] = (LUT) std::min(std::max(0,int(((adc2fC_*gain_)-ped_)*Rcalib[abs(ieta)+45]*gain_/lsb_/cosheta_[abs(ieta)])-3), 0x3FF);
0204 //if (ieta < 0 ) inputLUT[id][j] = (LUT) std::min(std::max(0,int((adc2fC_ - ped_)*Rcalib[abs(ieta)+2]*gain_/lsb_/cosheta_[abs(ieta)])), 0x3FF);
0205 //else inputLUT[id][j] = (LUT) std::min(std::max(0,int((adc2fC_ - ped_)*Rcalib[abs(ieta)+45]*gain_/lsb_/cosheta_[abs(ieta)])), 0x3FF);
0206 }
0207 }
0208 }
0209 }
0210 }
0211 }
0212
0213 void HcaluLUTTPGCoder::adc2Linear(const HBHEDataFrame& df, IntegerCaloSamples& ics) const {
0214 int id = GetLUTID(df.id().subdet(), df.id().ieta(), df.id().iphi(), df.id().depth());
0215 if (inputLUT[id]==0) {
0216 throw cms::Exception("Missing Data") << "No LUT for " << df.id();
0217 } else {
0218 for (int i=0; i<df.size(); i++){
0219 if (df[i].adc() >= INPUT_LUT_SIZE || df[i].adc() < 0) throw cms::Exception("ADC overflow for HBHE tower: ") << i << " adc= " << df[i].adc();
0220 ics[i]=inputLUT[id][df[i].adc()];
0221 }
0222 }
0223 }
0224
0225 void HcaluLUTTPGCoder::adc2Linear(const HFDataFrame& df, IntegerCaloSamples& ics) const{
0226 int id = GetLUTID(df.id().subdet(), df.id().ieta(), df.id().iphi(), df.id().depth());
0227 if (inputLUT[id]==0) {
0228 throw cms::Exception("Missing Data") << "No LUT for " << df.id();
0229 } else {
0230 for (int i=0; i<df.size(); i++){
0231 if (df[i].adc() >= INPUT_LUT_SIZE || df[i].adc() < 0)
0232 throw cms::Exception("ADC overflow for HF tower: ") << i << " adc= " << df[i].adc();
0233 ics[i]=inputLUT[id][df[i].adc()];
0234 }
0235 }
0236 }
0237
0238 unsigned short HcaluLUTTPGCoder::adc2Linear(HcalQIESample sample, HcalDetId id) const {
0239 int ref = GetLUTID(id.subdet(), id.ieta(), id.iphi(), id.depth());
0240 return inputLUT[ref][sample.adc()];
0241 }
0242
0243 float HcaluLUTTPGCoder::getLUTPedestal(HcalDetId id) const {
0244 int ref = GetLUTID(id.subdet(), id.ieta(), id.iphi(), id.depth());
0245 return _ped[ref];
0246 }
0247
0248 float HcaluLUTTPGCoder::getLUTGain(HcalDetId id) const {
0249 int ref = GetLUTID(id.subdet(), id.ieta(), id.iphi(), id.depth());
0250 return _gain[ref];
0251 }
