File indexing completed on 2024-04-06 12:27:58
0001 #include "RecoTBCalo/HcalTBObjectUnpacker/interface/HcalTBTDCUnpacker.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include <cmath>
0004
0005
0006 static const int lcTriggerTime = 1;
0007 static const int lcTTCL1ATime = 2;
0008 static const int lcBeamCoincidence = 3;
0009 static const int lcLaserFlash = 4;
0010 static const int lcQIEPhase = 5;
0011
0012 static const int lcMuon1 = 90;
0013 static const int lcMuon2 = 91;
0014 static const int lcMuon3 = 92;
0015 static const int lcScint1 = 93;
0016 static const int lcScint2 = 94;
0017 static const int lcScint3 = 95;
0018 static const int lcScint4 = 96;
0019 static const int lcBeamHalo1 = 50;
0020 static const int lcBeamHalo2 = 51;
0021 static const int lcBeamHalo3 = 55;
0022 static const int lcBeamHalo4 = 63;
0023 static const int lcTOF1S = 129;
0024 static const int lcTOF1J = 130;
0025 static const int lcTOF2S = 131;
0026 static const int lcTOF2J = 132;
0027
0028 namespace hcaltb {
0029
0030 HcalTBTDCUnpacker::HcalTBTDCUnpacker(bool include_unmatched_hits)
0031 : includeUnmatchedHits_(include_unmatched_hits), dumpObs_(nullptr) {
0032
0033
0034 }
0035 void HcalTBTDCUnpacker::setCalib(const std::vector<std::vector<std::string> >& calibLines_) {
0036 for (int i = 0; i < 161; i++) {
0037 tdc_ped[i] = 0.;
0038 tdc_convers[i] = 1.;
0039 }
0040 for (unsigned int ii = 0; ii < calibLines_.size(); ii++) {
0041
0042 if (calibLines_[ii][0] == "TDC") {
0043 if (calibLines_[ii].size() == 4) {
0044 int channel = atoi(calibLines_[ii][1].c_str());
0045 tdc_ped[channel] = atof(calibLines_[ii][2].c_str());
0046 tdc_convers[channel] = atof(calibLines_[ii][3].c_str());
0047
0048 } else {
0049 throw cms::Exception("Incomplete configuration")
0050 << "Wrong TDC configuration format : expected 3 parameters, got " << calibLines_[ii].size() - 1;
0051 }
0052 }
0053
0054
0055 if (calibLines_[ii][0] == "WC") {
0056 if (calibLines_[ii].size() == 6) {
0057 int plane = atoi(calibLines_[ii][1].c_str());
0058 wc_[plane].b0 = atof(calibLines_[ii][2].c_str());
0059 wc_[plane].b1 = atof(calibLines_[ii][3].c_str());
0060 wc_[plane].mean = atof(calibLines_[ii][4].c_str());
0061 wc_[plane].sigma = atof(calibLines_[ii][5].c_str());
0062
0063
0064 } else {
0065 throw cms::Exception("Incomplete configuration")
0066 << "Wrong Wire Chamber configuration format : expected 5 parameters, got " << calibLines_[ii].size() - 1;
0067 }
0068 }
0069
0070 }
0071 }
0072
0073 void HcalTBTDCUnpacker::unpack(const FEDRawData& raw, HcalTBEventPosition& pos, HcalTBTiming& timing) const {
0074 std::vector<Hit> hits;
0075
0076 unpackHits(raw, hits, timing);
0077
0078 reconstructWC(hits, pos);
0079 reconstructTiming(hits, timing);
0080 }
0081
0082 struct ClassicTDCDataFormat {
0083 unsigned int cdfHeader0, cdfHeader1, cdfHeader2, cdfHeader3;
0084 unsigned int n_max_hits;
0085 unsigned int n_hits;
0086 unsigned int hits[2];
0087 };
0088
0089 struct CombinedTDCQDCDataFormat {
0090 unsigned int cdfHeader0, cdfHeader1, cdfHeader2, cdfHeader3;
0091 unsigned int n_qdc_hits;
0092 unsigned int n_tdc_hits;
0093 unsigned short qdc_values[4];
0094 };
0095
0096
0097
0098 void HcalTBTDCUnpacker::unpackHits(const FEDRawData& raw, std::vector<Hit>& hits, HcalTBTiming& timing) const {
0099 const ClassicTDCDataFormat* tdc = (const ClassicTDCDataFormat*)raw.data();
0100
0101 if (raw.size() < 3 * 8) {
0102 throw cms::Exception("Missing Data") << "No data in the TDC block";
0103 }
0104
0105 const unsigned int* hitbase = nullptr;
0106 unsigned int totalhits = 0;
0107
0108
0109 if (tdc->n_max_hits != 192) {
0110 const CombinedTDCQDCDataFormat* qdctdc = (const CombinedTDCQDCDataFormat*)raw.data();
0111 hitbase = (const unsigned int*)(qdctdc);
0112 hitbase += 6;
0113 hitbase += qdctdc->n_qdc_hits / 2;
0114 totalhits = qdctdc->n_tdc_hits & 0xFFFF;
0115
0116
0117
0118
0119 } else {
0120 hitbase = &(tdc->hits[0]);
0121 totalhits = tdc->n_hits;
0122 }
0123
0124 for (unsigned int i = 0; i < totalhits; i++) {
0125 Hit h;
0126 h.channel = (hitbase[i] & 0x7FC00000) >> 22;
0127 h.time = (hitbase[i] & 0xFFFFF) * tdc_convers[h.channel];
0128 hits.push_back(h);
0129
0130 }
0131
0132
0133 int v775[32];
0134 for (int i = 0; i < 32; i++)
0135 v775[i] = -1;
0136 if (tdc->n_max_hits != 192) {
0137 const CombinedTDCQDCDataFormat* qdctdc = (const CombinedTDCQDCDataFormat*)raw.data();
0138 hitbase = (const unsigned int*)(qdctdc);
0139 hitbase += 6;
0140 hitbase += qdctdc->n_qdc_hits / 2;
0141 hitbase += (qdctdc->n_tdc_hits & 0xFFFF);
0142 totalhits = (qdctdc->n_tdc_hits & 0xFFFF0000) >> 16;
0143 for (unsigned int i = 0; i < totalhits; i++) {
0144 Hit h;
0145
0146 h.channel = 129 + ((hitbase[i] & 0x3F0000) >> 16);
0147 h.time = (hitbase[i] & 0xFFF) * tdc_convers[h.channel];
0148 hits.push_back(h);
0149 if ((h.channel - 129) < 32)
0150 v775[(h.channel - 129)] = (hitbase[i] & 0xFFF);
0151
0152 }
0153 }
0154 timing.setV775(v775);
0155 }
0156
0157 void HcalTBTDCUnpacker::reconstructTiming(const std::vector<Hit>& hits, HcalTBTiming& timing) const {
0158 std::vector<Hit>::const_iterator j;
0159 double trigger_time = 0;
0160 double ttc_l1a_time = 0;
0161 double laser_flash = 0;
0162 double qie_phase = 0;
0163 double TOF1S_time = 0;
0164 double TOF1J_time = 0;
0165 double TOF2S_time = 0;
0166 double TOF2J_time = 0;
0167
0168 std::vector<double> m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits, bh1hits, bh2hits, bh3hits, bh4hits,
0169 beam_coinc;
0170
0171 for (j = hits.begin(); j != hits.end(); j++) {
0172 switch (j->channel) {
0173 case lcTriggerTime:
0174 trigger_time = j->time - tdc_ped[lcTriggerTime];
0175 break;
0176 case lcTTCL1ATime:
0177 ttc_l1a_time = j->time - tdc_ped[lcTTCL1ATime];
0178 break;
0179 case lcBeamCoincidence:
0180 beam_coinc.push_back(j->time - tdc_ped[lcBeamCoincidence]);
0181 break;
0182 case lcLaserFlash:
0183 laser_flash = j->time - tdc_ped[lcLaserFlash];
0184 break;
0185 case lcQIEPhase:
0186 qie_phase = j->time - tdc_ped[lcQIEPhase];
0187 break;
0188 case lcMuon1:
0189 m1hits.push_back(j->time - tdc_ped[lcMuon1]);
0190 break;
0191 case lcMuon2:
0192 m2hits.push_back(j->time - tdc_ped[lcMuon2]);
0193 break;
0194 case lcMuon3:
0195 m3hits.push_back(j->time - tdc_ped[lcMuon3]);
0196 break;
0197 case lcScint1:
0198 s1hits.push_back(j->time - tdc_ped[lcScint1]);
0199 break;
0200 case lcScint2:
0201 s2hits.push_back(j->time - tdc_ped[lcScint2]);
0202 break;
0203 case lcScint3:
0204 s3hits.push_back(j->time - tdc_ped[lcScint3]);
0205 break;
0206 case lcScint4:
0207 s4hits.push_back(j->time - tdc_ped[lcScint4]);
0208 break;
0209 case lcTOF1S:
0210 TOF1S_time = j->time - tdc_ped[lcTOF1S];
0211 break;
0212 case lcTOF1J:
0213 TOF1J_time = j->time - tdc_ped[lcTOF1J];
0214 break;
0215 case lcTOF2S:
0216 TOF2S_time = j->time - tdc_ped[lcTOF2S];
0217 break;
0218 case lcTOF2J:
0219 TOF2J_time = j->time - tdc_ped[lcTOF2J];
0220 break;
0221 case lcBeamHalo1:
0222 bh1hits.push_back(j->time - tdc_ped[lcBeamHalo1]);
0223 break;
0224 case lcBeamHalo2:
0225 bh2hits.push_back(j->time - tdc_ped[lcBeamHalo2]);
0226 break;
0227 case lcBeamHalo3:
0228 bh3hits.push_back(j->time - tdc_ped[lcBeamHalo3]);
0229 break;
0230 case lcBeamHalo4:
0231 bh4hits.push_back(j->time - tdc_ped[lcBeamHalo4]);
0232 break;
0233 default:
0234 break;
0235 }
0236 }
0237
0238 timing.setTimes(trigger_time, ttc_l1a_time, laser_flash, qie_phase, TOF1S_time, TOF1J_time, TOF2S_time, TOF2J_time);
0239 timing.setHits(
0240 m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits, bh1hits, bh2hits, bh3hits, bh4hits, beam_coinc);
0241 }
0242
0243 const int HcalTBTDCUnpacker::WC_CHANNELIDS[PLANECOUNT * 3] = {
0244 12, 13, 14,
0245 10, 11, 14,
0246 22, 23, 24,
0247 20, 21, 24,
0248 32, 33, 34,
0249 30, 31, 34,
0250 101, 102, 104,
0251 107, 108, 110,
0252 113, 114, 116,
0253 97, 98, 99,
0254 42, 43, -1,
0255 44, 60, -1,
0256 40, 41, -1,
0257 45, 61, -1,
0258 52, 53, -1,
0259 54, 62, -1
0260 };
0261
0262 static const double TDC_OFFSET_CONSTANT = 12000;
0263 static const double N_SIGMA = 2.5;
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297 void HcalTBTDCUnpacker::reconstructWC(const std::vector<Hit>& hits, HcalTBEventPosition& ep) const {
0298
0299 const int MAX_HITS = 100;
0300 float hits1[MAX_HITS], hits2[MAX_HITS], hitsA[MAX_HITS];
0301 int n1, n2, nA, chan1, chan2, chanA;
0302
0303 std::vector<double> x;
0304
0305 for (int plane = 0; plane < PLANECOUNT; plane++) {
0306 n1 = 0;
0307 n2 = 0;
0308 nA = 0;
0309
0310 std::vector<double> plane_hits;
0311 double hit_time;
0312 hits1[0] = -1000;
0313 hits2[0] = -1000;
0314 hitsA[0] = -1000;
0315 hitsA[1] = -1000;
0316
0317 chan1 = WC_CHANNELIDS[plane * 3];
0318 chan2 = WC_CHANNELIDS[plane * 3 + 1];
0319 chanA = WC_CHANNELIDS[plane * 3 + 2];
0320
0321 for (std::vector<Hit>::const_iterator j = hits.begin(); j != hits.end(); j++) {
0322 if (j->channel == chan1 && n1 < MAX_HITS) {
0323 hits1[n1] = j->time - TDC_OFFSET_CONSTANT;
0324 n1++;
0325 }
0326 if (j->channel == chan2 && n2 < MAX_HITS) {
0327 hits2[n2] = j->time - TDC_OFFSET_CONSTANT;
0328 n2++;
0329 }
0330 if (j->channel == chanA && nA < MAX_HITS) {
0331 hitsA[nA] = j->time - TDC_OFFSET_CONSTANT;
0332 nA++;
0333 }
0334 }
0335
0336 if (n1 != 0 && n2 != 0 && dumpObs_ != nullptr) {
0337 fprintf(dumpObs_, "%d,%f,%f,%f,%f\n", plane, hits1[0], hits2[0], hitsA[0], hitsA[1]);
0338 fflush(dumpObs_);
0339 }
0340
0341
0342 for (int ii = 0; ii < n1; ii++) {
0343 int jmin = -1, lmin = -1;
0344 float dsumMin = 99999;
0345 for (int jj = 0; jj < n2; jj++) {
0346 for (int ll = 0; ll < nA; ll++) {
0347 float dsum = fabs(wc_[plane].mean - hits1[ii] - hits2[jj] + 2.0 * hitsA[ll]);
0348 if (dsum < (N_SIGMA * wc_[plane].sigma) && dsum < dsumMin) {
0349 jmin = jj;
0350 lmin = ll;
0351 dsumMin = dsum;
0352 }
0353 }
0354 }
0355 if (jmin >= 0) {
0356 hit_time = wc_[plane].b0 + wc_[plane].b1 * (hits1[ii] - hits2[jmin]);
0357 if ((plane % 2) == 0) {
0358 plane_hits.push_back(-hit_time);
0359 } else {
0360 plane_hits.push_back(hit_time);
0361 }
0362 hits1[ii] = -99999;
0363 hits2[jmin] = -99999;
0364 hitsA[lmin] = 99999;
0365 }
0366 }
0367
0368 if (includeUnmatchedHits_ || plane > 9)
0369 for (int ii = 0; ii < n1; ii++) {
0370 if (hits1[ii] < -99990)
0371 continue;
0372 for (int jj = 0; jj < n2; jj++) {
0373 if (hits2[jj] < -99990)
0374 continue;
0375 hit_time = wc_[plane].b0 + wc_[plane].b1 * (hits1[ii] - hits2[jj]);
0376 if ((plane % 2) == 0) {
0377 plane_hits.push_back(-hit_time);
0378 } else {
0379 plane_hits.push_back(hit_time);
0380 }
0381 }
0382 }
0383
0384 if ((plane % 2) == 0)
0385 x = plane_hits;
0386 else {
0387 char chamber = 'A' + plane / 2;
0388 ep.setChamberHits(chamber, x, plane_hits);
0389 }
0390 }
0391 }
0392
0393 }