-
Notifications
You must be signed in to change notification settings - Fork 4.3k
/
CaloTPGTranscoderULUT.cc
320 lines (277 loc) · 12.9 KB
/
CaloTPGTranscoderULUT.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
#include "CalibCalorimetry/CaloTPG/src/CaloTPGTranscoderULUT.h"
#include "DataFormats/HcalDetId/interface/HcalGenericDetId.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Framework/interface/ESTransientHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include <iostream>
#include <fstream>
#include <cmath>
//#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "CondFormats/DataRecord/interface/HcalLutMetadataRcd.h"
#include "CalibCalorimetry/HcalTPGAlgos/interface/HcaluLUTTPGCoder.h"
using namespace std;
CaloTPGTranscoderULUT::CaloTPGTranscoderULUT(const std::string& compressionFile,
const std::string& decompressionFile)
: theTopology(nullptr),
nominal_gain_(0.), lsb_factor_(0.), rct_factor_(1.), nct_factor_(1.),
lin8_factor_(1.), lin11_factor_(1.),
compressionFile_(compressionFile),
decompressionFile_(decompressionFile)
{
}
CaloTPGTranscoderULUT::~CaloTPGTranscoderULUT() {
}
void CaloTPGTranscoderULUT::loadHCALCompress(HcalLutMetadata const& lutMetadata,
HcalTrigTowerGeometry const& theTrigTowerGeometry) {
if (!theTopology) {
throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
}
std::array<unsigned int, OUTPUT_LUT_SIZE> analyticalLUT;
std::array<unsigned int, OUTPUT_LUT_SIZE> linearQIE8LUT;
std::array<unsigned int, OUTPUT_LUT_SIZE> linearQIE11LUT;
std::array<unsigned int, OUTPUT_LUT_SIZE> linearRctLUT;
std::array<unsigned int, OUTPUT_LUT_SIZE> linearNctLUT;
// Compute compression LUT
for (unsigned int i=0; i < OUTPUT_LUT_SIZE; i++) {
analyticalLUT[i] = min(static_cast<unsigned int>(sqrt(14.94 * log(1. + i / 14.94) * i) + 0.5), TPGMAX - 1);
linearQIE8LUT[i] = min(static_cast<unsigned int>((i + .5) / lin8_factor_), TPGMAX - 1);
linearQIE11LUT[i] = min(static_cast<unsigned int>((i + .5) / lin11_factor_), TPGMAX - 1);
linearRctLUT[i] = min(static_cast<unsigned int>((i + .5) / rct_factor_), TPGMAX - 1);
linearNctLUT[i] = min(static_cast<unsigned int>((i + .5) / nct_factor_), TPGMAX - 1);
}
std::vector<DetId> allChannels = lutMetadata.getAllChannels();
for(std::vector<DetId>::iterator i=allChannels.begin(); i!=allChannels.end(); ++i){
if (not HcalGenericDetId(*i).isHcalTrigTowerDetId()) {
if ((not HcalGenericDetId(*i).isHcalDetId()) and
(not HcalGenericDetId(*i).isHcalZDCDetId()) and
(not HcalGenericDetId(*i).isHcalCastorDetId()))
edm::LogWarning("CaloTPGTranscoderULUT") << "Encountered invalid HcalDetId " << HcalGenericDetId(*i);
continue;
}
HcalTrigTowerDetId id(*i);
if(!theTopology->validHT(id)) continue;
unsigned int index = getOutputLUTId(id);
const HcalLutMetadatum *meta = lutMetadata.getValues(id);
unsigned int threshold = meta->getOutputLutThreshold();
int ieta=id.ieta();
int version=id.version();
bool isHBHE = (abs(ieta) < theTrigTowerGeometry.firstHFTower(version));
unsigned int lutsize = getOutputLUTSize(id);
outputLUT_[index].resize(lutsize);
for (unsigned int i = 0; i < threshold; ++i)
outputLUT_[index][i] = 0;
if (isHBHE) {
for (unsigned int i = threshold; i < lutsize; ++i)
if (allLinear_) {
outputLUT_[index][i] = isOnlyQIE11(id) ? linearQIE11LUT[i] : linearQIE8LUT[i];
} else {
outputLUT_[index][i] = analyticalLUT[i];
}
} else {
for (unsigned int i = threshold; i < lutsize; ++i)
outputLUT_[index][i] = version == 0 ? linearRctLUT[i] : linearNctLUT[i];
}
double eta_low = 0., eta_high = 0.;
theTrigTowerGeometry.towerEtaBounds(ieta,version,eta_low,eta_high);
double cosh_ieta = fabs(cosh((eta_low + eta_high)/2.));
double granularity = meta->getLutGranularity();
if(isHBHE){
if (allLinear_) {
LUT tpg = outputLUT_[index][0];
hcaluncomp_[index][tpg] = 0;
for (unsigned int i = 0; i < lutsize; ++i){
if (outputLUT_[index][i] != tpg){
tpg = outputLUT_[index][i];
hcaluncomp_[index][tpg] = lsb_factor_ * i / (isOnlyQIE11(id) ? lin11_factor_ : lin8_factor_);
}
}
} else {
double factor = nominal_gain_ / cosh_ieta * granularity;
LUT tpg = outputLUT_[index][0];
int low = 0;
for (unsigned int i = 0; i < lutsize; ++i){
if (outputLUT_[index][i] != tpg){
unsigned int mid = (low + i)/2;
hcaluncomp_[index][tpg] = (tpg == 0 ? low : factor * mid);
low = i;
tpg = outputLUT_[index][i];
}
}
hcaluncomp_[index][tpg] = factor * low;
}
}
else{
LUT tpg = outputLUT_[index][0];
hcaluncomp_[index][tpg]=0;
for (unsigned int i = 0; i < lutsize; ++i){
if (outputLUT_[index][i] != tpg){
tpg = outputLUT_[index][i];
hcaluncomp_[index][tpg] = lsb_factor_ * i / (version==0?rct_factor_:nct_factor_);
}
}
}
}
}
HcalTriggerPrimitiveSample CaloTPGTranscoderULUT::hcalCompress(const HcalTrigTowerDetId& id, unsigned int sample, int fineGrain) const {
unsigned int itower = getOutputLUTId(id);
if (sample >= getOutputLUTSize(id))
throw cms::Exception("Out of Range")
<< "LUT has " << getOutputLUTSize(id) << " entries for " << id << " but " << sample << " was requested.";
if(itower >= outputLUT_.size())
throw cms::Exception("Out of Range") << "No decompression LUT found for " << id;
return HcalTriggerPrimitiveSample(outputLUT_[itower][sample], fineGrain);
}
double CaloTPGTranscoderULUT::hcaletValue(const int& ieta, const int& iphi, const int& version, const int& compET) const {
double etvalue = 0.;
int itower = getOutputLUTId(ieta,iphi, version);
if (itower < 0) {
edm::LogError("CaloTPGTranscoderULUT") << "No decompression LUT found for ieta, iphi = " << ieta << ", " << iphi;
} else if (compET < 0 || compET >= (int) TPGMAX) {
edm::LogError("CaloTPGTranscoderULUT") << "Compressed value out of range: eta, phi, cET = " << ieta << ", " << iphi << ", " << compET;
} else {
etvalue = hcaluncomp_[itower][compET];
}
return(etvalue);
}
double CaloTPGTranscoderULUT::hcaletValue(const HcalTrigTowerDetId& hid, const HcalTriggerPrimitiveSample& hc) const {
int compET = hc.compressedEt(); // to be within the range by the class
int itower = getOutputLUTId(hid);
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 version) const {
HcalTrigTowerDetId id(ieta, iphiin);
id.setVersion(version);
if (!theTopology) {
throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
}
return theTopology->validHT(id);
}
int CaloTPGTranscoderULUT::getOutputLUTId(const HcalTrigTowerDetId& id) const {
if (!theTopology) {
throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
}
return theTopology->detId2denseIdHT(id);
}
int CaloTPGTranscoderULUT::getOutputLUTId(const int ieta, const int iphiin, const int version) const {
if (!theTopology) {
throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
}
HcalTrigTowerDetId id(ieta, iphiin);
id.setVersion(version);
return theTopology->detId2denseIdHT(id);
}
unsigned int
CaloTPGTranscoderULUT::getOutputLUTSize(const HcalTrigTowerDetId& id) const
{
if (!theTopology)
throw cms::Exception("CaloTPGTranscoderULUT")
<< "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
switch (theTopology->triggerMode()) {
case HcalTopologyMode::TriggerMode_2009:
case HcalTopologyMode::TriggerMode_2016:
return QIE8_OUTPUT_LUT_SIZE;
case HcalTopologyMode::TriggerMode_2017:
if (id.ietaAbs() <= theTopology->lastHERing())
return QIE8_OUTPUT_LUT_SIZE;
else
return QIE10_OUTPUT_LUT_SIZE;
case HcalTopologyMode::TriggerMode_2017plan1:
if (plan1_towers_.find(id) != plan1_towers_.end())
return QIE11_OUTPUT_LUT_SIZE;
else if (id.ietaAbs() <= theTopology->lastHERing())
return QIE8_OUTPUT_LUT_SIZE;
else
return QIE10_OUTPUT_LUT_SIZE;
case HcalTopologyMode::TriggerMode_2018legacy:
case HcalTopologyMode::TriggerMode_2018:
if (id.ietaAbs() <= theTopology->lastHBRing())
return QIE8_OUTPUT_LUT_SIZE;
else if (id.ietaAbs() <= theTopology->lastHERing())
return QIE11_OUTPUT_LUT_SIZE;
else
return QIE10_OUTPUT_LUT_SIZE;
case HcalTopologyMode::TriggerMode_2019:
if (id.ietaAbs() <= theTopology->lastHERing())
return QIE11_OUTPUT_LUT_SIZE;
else
return QIE10_OUTPUT_LUT_SIZE;
default:
throw cms::Exception("CaloTPGTranscoderULUT")
<< "Unknown trigger mode used by the topology!";
}
}
bool
CaloTPGTranscoderULUT::isOnlyQIE11(const HcalTrigTowerDetId& id) const
{
if (!theTopology)
throw cms::Exception("CaloTPGTranscoderULUT")
<< "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
switch (theTopology->triggerMode()) {
case HcalTopologyMode::TriggerMode_2017plan1:
if (plan1_towers_.find(id) != plan1_towers_.end() and id.ietaAbs() != theTopology->lastHBRing())
return true;
return false;
case HcalTopologyMode::TriggerMode_2018legacy:
case HcalTopologyMode::TriggerMode_2018:
if (id.ietaAbs() <= theTopology->lastHBRing())
return false;
return true;
case HcalTopologyMode::TriggerMode_2019:
return true;
default:
throw cms::Exception("CaloTPGTranscoderULUT")
<< "Unknown trigger mode used by the topology!";
}
}
const std::vector<unsigned int> CaloTPGTranscoderULUT::getCompressionLUT(const HcalTrigTowerDetId& id) const {
int itower = getOutputLUTId(id);
auto lut = outputLUT_[itower];
std::vector<unsigned int> result(lut.begin(), lut.end());
return result;
}
void CaloTPGTranscoderULUT::setup(HcalLutMetadata const& lutMetadata, HcalTrigTowerGeometry const& theTrigTowerGeometry, int nctScaleShift, int rctScaleShift, double lsbQIE8, double lsbQIE11, bool allLinear)
{
theTopology = lutMetadata.topo();
nominal_gain_ = lutMetadata.getNominalGain();
lsb_factor_ = lutMetadata.getRctLsb();
allLinear_ = allLinear;
rct_factor_ = lsb_factor_/(HcaluLUTTPGCoder::lsb_*(1<<rctScaleShift));
nct_factor_ = lsb_factor_/(HcaluLUTTPGCoder::lsb_*(1<<nctScaleShift));
lin8_factor_ = lsb_factor_/lsbQIE8;
lin11_factor_ = lsb_factor_/lsbQIE11;
outputLUT_.resize(theTopology->getHTSize());
hcaluncomp_.resize(theTopology->getHTSize());
plan1_towers_.clear();
for (const auto& id: lutMetadata.getAllChannels()) {
if (not (id.det() == DetId::Hcal and theTopology->valid(id)))
continue;
HcalDetId cell(id);
if (not theTopology->dddConstants()->isPlan1(cell))
continue;
for (const auto& tower: theTrigTowerGeometry.towerIds(cell))
plan1_towers_.emplace(tower);
}
if (compressionFile_.empty() && decompressionFile_.empty()) {
loadHCALCompress(lutMetadata,theTrigTowerGeometry);
}
else {
throw cms::Exception("Not Implemented") << "setup of CaloTPGTranscoderULUT from text files";
}
}