File indexing completed on 2024-04-06 12:19:38
0001
0002
0003
0004
0005 #include "L1Trigger/CSCTriggerPrimitives/interface/CSCPatternBank.h"
0006
0007 #include <TCanvas.h>
0008 #include <TGraph.h>
0009 #include <TGraphErrors.h>
0010 #include <TF1.h>
0011 #include <TSystem.h>
0012 #include <TH1F.h>
0013 #include <TFile.h>
0014 #include <TString.h>
0015
0016
0017 #include <fstream>
0018 #include <iostream>
0019 #include <iomanip>
0020 #include <math.h>
0021 #include <memory>
0022 #include <bitset>
0023
0024 using namespace std;
0025
0026
0027 const float HALF_STRIP_ERROR = 0.288675;
0028 const unsigned halfpatternwidth = (CSCConstants::CLCT_PATTERN_WIDTH - 1) / 2;
0029
0030
0031 const int DEBUG = 0;
0032
0033
0034 class CSCPattern {
0035 public:
0036 CSCPattern(unsigned int id, const CSCPatternBank::LCTPattern& pat);
0037 ~CSCPattern() {}
0038
0039 void printCode(const unsigned code) const;
0040 void getLayerPattern(const unsigned code, unsigned layerPattern[CSCConstants::NUM_LAYERS]) const;
0041 int recoverPatternCCCombination(const unsigned code,
0042 int code_hits[CSCConstants::NUM_LAYERS][CSCConstants::CLCT_PATTERN_WIDTH]) const;
0043 string getName() const { return name_; }
0044 unsigned getId() const { return id_; }
0045
0046 private:
0047 const unsigned int id_;
0048 string name_;
0049 CSCPatternBank::LCTPattern pat_;
0050 };
0051
0052 CSCPattern::CSCPattern(unsigned int id, const CSCPatternBank::LCTPattern& pat)
0053 : id_(id), name_(std::to_string(id)), pat_(pat) {}
0054
0055
0056 void CSCPattern::printCode(const unsigned code) const {
0057
0058 unsigned layerPattern[CSCConstants::NUM_LAYERS];
0059 getLayerPattern(code, layerPattern);
0060
0061 std::cout << "Pattern " << id_ << ", Code " << code << " " << std::bitset<12>(code) << std::endl;
0062 for (unsigned int j = 0; j < CSCConstants::NUM_LAYERS; j++) {
0063 unsigned trueCounter = 0;
0064 std::cout << "L" << j + 1 << ": ";
0065 for (unsigned int i = 0; i < CSCConstants::CLCT_PATTERN_WIDTH; i++) {
0066 if (!pat_[j][i]) {
0067 printf("-");
0068 } else {
0069 trueCounter++;
0070 if (trueCounter == layerPattern[j])
0071 printf("X");
0072 else
0073 printf("_");
0074 }
0075 }
0076 std::cout << " -> " << std::bitset<2>(layerPattern[j]) << std::endl;
0077 }
0078
0079 printf(" ");
0080 for (unsigned int i = 0; i < CSCConstants::CLCT_PATTERN_WIDTH; i++) {
0081 printf("%1x", i);
0082 }
0083 printf("\n");
0084 }
0085
0086 void CSCPattern::getLayerPattern(const unsigned code, unsigned layerPattern[CSCConstants::NUM_LAYERS]) const {
0087
0088 std::bitset<12> cc(code);
0089
0090 for (unsigned ilayer = 0; ilayer < CSCConstants::NUM_LAYERS; ilayer++) {
0091 std::bitset<2> cclayer;
0092 cclayer[1] = cc[2 * ilayer + 1];
0093 cclayer[0] = cc[2 * ilayer];
0094 layerPattern[ilayer] = cclayer.to_ulong();
0095 }
0096 }
0097
0098
0099 int CSCPattern::recoverPatternCCCombination(
0100 const unsigned code, int code_hits[CSCConstants::NUM_LAYERS][CSCConstants::CLCT_PATTERN_WIDTH]) const {
0101
0102 unsigned layerPattern[CSCConstants::NUM_LAYERS];
0103 getLayerPattern(code, layerPattern);
0104
0105
0106 for (unsigned int j = 0; j < CSCConstants::NUM_LAYERS; j++) {
0107 unsigned trueCounter = 0;
0108 for (unsigned int i = 0; i < CSCConstants::CLCT_PATTERN_WIDTH; i++) {
0109
0110 if (!pat_[j][i]) {
0111 code_hits[j][i] = 0;
0112 }
0113
0114 else {
0115 trueCounter++;
0116 if (trueCounter == layerPattern[j])
0117 code_hits[j][i] = 1;
0118 else
0119 code_hits[j][i] = 0;
0120 }
0121 }
0122 }
0123 return 0;
0124 }
0125
0126 int convertToLegacyPattern(const int code_hits[CSCConstants::NUM_LAYERS][CSCConstants::CLCT_PATTERN_WIDTH],
0127 unsigned N_LAYER_REQUIREMENT);
0128 void getErrors(const vector<float>& x, const vector<float>& y, float& sigmaM, float& sigmaB);
0129 void writeHeaderPosOffsetLUT(ofstream& file);
0130 void writeHeaderSlopeLUT(ofstream& file);
0131 unsigned firmwareWord(const unsigned quality, const unsigned slope, const unsigned offset);
0132 void setDataWord(unsigned& word, const unsigned newWord, const unsigned shift, const unsigned mask);
0133 unsigned assignPosition(const float fvalue, const float fmin, const float fmax, const unsigned nbits);
0134 unsigned assignBending(const float fvalue);
0135
0136 int CCLUTLinearFitWriter(unsigned N_LAYER_REQUIREMENT = 3) {
0137
0138 std::unique_ptr<std::vector<CSCPattern>> newPatterns(new std::vector<CSCPattern>());
0139 for (unsigned ipat = 0; ipat < 5; ipat++) {
0140 newPatterns->emplace_back(ipat, CSCPatternBank::clct_pattern_run3_[ipat]);
0141 }
0142
0143
0144 const std::string outdir("output_" + std::to_string(N_LAYER_REQUIREMENT) + "layers/");
0145
0146
0147 TFile* file = TFile::Open(TString::Format("figures_%dlayers/fitresults.root", N_LAYER_REQUIREMENT), "RECREATE");
0148
0149 TH1D* all_offsets = new TH1D("all_offsets", "All Half-strip offsets", 200, -2.5, 2.5);
0150 all_offsets->GetYaxis()->SetTitle("Entries");
0151 all_offsets->GetXaxis()->SetTitle("Half-strip offset");
0152
0153 TH1D* all_slopes = new TH1D("all_slopes", "All slopes", 200, -3, 3);
0154 all_slopes->GetYaxis()->SetTitle("Entries");
0155 all_slopes->GetXaxis()->SetTitle("Slope [half-strips/layer]");
0156
0157 std::map<unsigned, TH1D*> offsets;
0158 std::map<unsigned, TH1D*> slopes;
0159 std::map<unsigned, TH1D*> offsetuncs;
0160 std::map<unsigned, TH1D*> slopeuncs;
0161 std::map<unsigned, TH1D*> chi2s;
0162 std::map<unsigned, TH1D*> chi2ndfs;
0163 std::map<unsigned, TH1D*> layers;
0164 std::map<unsigned, TH1D*> legacypatterns;
0165
0166 for (unsigned i = 0; i < 5; ++i) {
0167 std::string title("Half-strip offset, PID " + std::to_string(i));
0168 std::string name("offset_pat" + std::to_string(i));
0169
0170 offsets[i] = new TH1D(name.c_str(), title.c_str(), 100, -2.5, 2.5);
0171 offsets[i]->GetYaxis()->SetTitle("Entries");
0172 offsets[i]->GetXaxis()->SetTitle("Half-strip offset");
0173
0174 title = "Slope, PID " + std::to_string(i);
0175 name = "slope_pat" + std::to_string(i);
0176 slopes[i] = new TH1D(name.c_str(), title.c_str(), 100, -3, 3);
0177 slopes[i]->GetYaxis()->SetTitle("Entries");
0178 slopes[i]->GetXaxis()->SetTitle("Slope [half-strips/layer]");
0179
0180 title = "Half-strip offset uncertainty, PID " + std::to_string(i);
0181 name = "offsetunc_pat" + std::to_string(i);
0182 offsetuncs[i] = new TH1D(name.c_str(), title.c_str(), 100, -2, 2);
0183 offsetuncs[i]->GetYaxis()->SetTitle("Entries");
0184 offsetuncs[i]->GetXaxis()->SetTitle("Half-strip offset uncertainty");
0185
0186 title = "Slope uncertainty, PID " + std::to_string(i);
0187 name = "slopeunc_pat" + std::to_string(i);
0188 slopeuncs[i] = new TH1D(name.c_str(), title.c_str(), 100, -2, 2);
0189 slopeuncs[i]->GetYaxis()->SetTitle("Entries");
0190 slopeuncs[i]->GetXaxis()->SetTitle("Slope uncertainty [half-strips/layer]");
0191
0192 title = "Chi2, PID " + std::to_string(i);
0193 name = "chi2_pat" + std::to_string(i);
0194 chi2s[i] = new TH1D(name.c_str(), title.c_str(), 100, 0, 100);
0195 chi2s[i]->GetYaxis()->SetTitle("Entries");
0196 chi2s[i]->GetXaxis()->SetTitle("#chi^{2}");
0197
0198 title = "Chi2/ndf, PID " + std::to_string(i);
0199 name = "chi2ndf_pat" + std::to_string(i);
0200 chi2ndfs[i] = new TH1D(name.c_str(), title.c_str(), 40, 0, 20);
0201 chi2ndfs[i]->GetYaxis()->SetTitle("Entries");
0202 chi2ndfs[i]->GetXaxis()->SetTitle("#chi^{2}/NDF");
0203
0204 title = "Number of layers, PID " + std::to_string(i);
0205 name = "layers_pat" + std::to_string(i);
0206 layers[i] = new TH1D(name.c_str(), title.c_str(), 7, 0, 7);
0207 layers[i]->GetYaxis()->SetTitle("Entries");
0208 layers[i]->GetXaxis()->SetTitle("Number of layers in pattern");
0209
0210 title = "Run-1/2 pattern, PID " + std::to_string(i);
0211 name = "legacypatterns_pat" + std::to_string(i);
0212 legacypatterns[i] = new TH1D(name.c_str(), title.c_str(), 11, 0, 11);
0213 legacypatterns[i]->GetYaxis()->SetTitle("Entries");
0214 legacypatterns[i]->GetXaxis()->SetTitle("Equivalent Run-1/2 pattern number");
0215 }
0216
0217 for (auto& p : offsets) {
0218 (p.second)->GetYaxis()->SetTitle("Entries");
0219 }
0220
0221
0222 float maxOffset = -1;
0223 float maxPatt = -1;
0224 float maxCode = -1;
0225 float minOffset = -2;
0226 float minPatt = 0;
0227 float minCode = 1;
0228
0229
0230 for (auto patt = newPatterns->begin(); patt != newPatterns->end(); ++patt) {
0231 std::cout << "Processing pattern " << patt - newPatterns->begin() << std::endl;
0232
0233
0234 std::cout << "Create output files for pattern " << patt->getName() << std::endl;
0235
0236
0237 ofstream outoffset_sw;
0238 outoffset_sw.open(outdir + "CSCComparatorCodePosOffsetLUT_pat" + patt->getName() + "_float_v1.txt");
0239 writeHeaderPosOffsetLUT(outoffset_sw);
0240
0241 ofstream outslope_sw;
0242 outslope_sw.open(outdir + "CSCComparatorCodeSlopeLUT_pat" + patt->getName() + "_float_v1.txt");
0243 writeHeaderSlopeLUT(outslope_sw);
0244
0245
0246 ofstream outoffset_sw_bin;
0247 outoffset_sw_bin.open(outdir + "CSCComparatorCodePosOffsetLUT_pat" + patt->getName() + "_v1.txt");
0248 writeHeaderPosOffsetLUT(outoffset_sw_bin);
0249
0250 ofstream outslope_sw_bin;
0251 outslope_sw_bin.open(outdir + "CSCComparatorCodeSlopeLUT_pat" + patt->getName() + "_v1.txt");
0252 writeHeaderSlopeLUT(outslope_sw_bin);
0253
0254
0255 ofstream outfile_fw;
0256 outfile_fw.open(outdir + "rom_pat" + patt->getName() + ".mem");
0257
0258
0259 ofstream outpatternconv;
0260 outpatternconv.open(outdir + "CSCComparatorCodePatternConversionLUT_pat" + patt->getName() + "_v1.txt");
0261 outpatternconv << "#<header> v1.0 12 32 </header>\n";
0262
0263
0264 for (unsigned code = 0; code < CSCConstants::NUM_COMPARATOR_CODES; code++) {
0265 if (DEBUG > 0) {
0266 cout << "Evaluating..." << endl;
0267 patt->printCode(code);
0268 }
0269 int hits[CSCConstants::NUM_LAYERS][CSCConstants::CLCT_PATTERN_WIDTH];
0270
0271 if (patt->recoverPatternCCCombination(code, hits)) {
0272 cout << "Error: CC evaluation has failed" << endl;
0273 }
0274
0275
0276
0277 vector<float> x;
0278 vector<float> y;
0279 vector<float> xe;
0280 vector<float> ye;
0281
0282 for (int i = 0; i < CSCConstants::NUM_LAYERS; i++) {
0283 for (unsigned j = 0; j < CSCConstants::CLCT_PATTERN_WIDTH; j++) {
0284 if (hits[i][j]) {
0285
0286 x.push_back(i - 2);
0287 y.push_back(float(j));
0288 xe.push_back(float(0));
0289 ye.push_back(HALF_STRIP_ERROR);
0290 if (DEBUG > 0)
0291 cout << "L " << x.back() << " HS " << y.back() << endl;
0292 }
0293 }
0294 if (DEBUG > 0)
0295 cout << endl;
0296 }
0297
0298
0299 float offset = 0;
0300 float slope = 0;
0301 float offsetunc = 0;
0302 float slopeunc = 0;
0303 float chi2 = -1;
0304 unsigned ndf;
0305 unsigned layer = 0;
0306 unsigned legacypattern = 0;
0307
0308
0309 if (x.size() >= N_LAYER_REQUIREMENT) {
0310
0311 std::unique_ptr<TGraphErrors> gr(new TGraphErrors(x.size(), &x[0], &y[0], &xe[0], &ye[0]));
0312 std::unique_ptr<TF1> fit(new TF1("fit", "pol1", -3, 4));
0313 gr->Fit("fit", "EMQ");
0314
0315
0316
0317 offset = fit->GetParameter(0) - halfpatternwidth;
0318 slope = fit->GetParameter(1);
0319 chi2 = fit->GetChisquare();
0320 ndf = fit->GetNDF();
0321 layer = ndf + 2;
0322 legacypattern = convertToLegacyPattern(hits, N_LAYER_REQUIREMENT);
0323
0324 getErrors(x, y, offsetunc, slopeunc);
0325
0326
0327 const unsigned ipat(patt - newPatterns->begin());
0328 offsets[ipat]->Fill(offset);
0329 slopes[ipat]->Fill(slope);
0330 offsetuncs[ipat]->Fill(offsetunc);
0331 slopeuncs[ipat]->Fill(slopeunc);
0332 chi2s[ipat]->Fill(chi2);
0333 chi2ndfs[ipat]->Fill(chi2 / ndf);
0334 layers[ipat]->Fill(layer);
0335 legacypatterns[ipat]->Fill(legacypattern);
0336
0337 all_offsets->Fill(offset);
0338 all_slopes->Fill(slope);
0339 }
0340
0341
0342 const float fmaxOffset = 2;
0343 const float fminOffset = -1.75;
0344
0345
0346
0347 const bool slope_sign(slope >= 0);
0348
0349 const unsigned offset_bin = assignPosition(offset, fminOffset, fmaxOffset, 4);
0350 unsigned slope_bin = assignBending(std::abs(slope));
0351 if (slope_sign)
0352 slope_bin += 16;
0353 const unsigned fwword = firmwareWord(0, slope_bin, offset_bin);
0354
0355
0356 outoffset_sw << code << " " << offset << "\n";
0357 outslope_sw << code << " " << slope << "\n";
0358 outoffset_sw_bin << code << " " << offset_bin << "\n";
0359 outslope_sw_bin << code << " " << slope_bin << "\n";
0360 outpatternconv << code << " " << legacypattern << "\n";
0361 outfile_fw << setfill('0');
0362 outfile_fw << setw(5) << std::hex << fwword << "\n";
0363
0364
0365 if (layer >= N_LAYER_REQUIREMENT) {
0366 if (offset < minOffset) {
0367 minOffset = offset;
0368 minPatt = patt->getId();
0369 minCode = code;
0370 }
0371 if (offset > maxOffset) {
0372 maxOffset = offset;
0373 maxPatt = patt->getId();
0374 maxCode = code;
0375 }
0376 }
0377 }
0378
0379
0380 outoffset_sw.close();
0381 outslope_sw.close();
0382 outoffset_sw_bin.close();
0383 outslope_sw_bin.close();
0384 outfile_fw.close();
0385 }
0386
0387 cout << "minOffset = " << minOffset << endl;
0388 for (auto patt = newPatterns->begin(); patt != newPatterns->end(); ++patt) {
0389 if (patt->getId() == minPatt) {
0390 patt->printCode(minCode);
0391 }
0392 }
0393 cout << "maxOffset = " << maxOffset << endl;
0394 for (auto patt = newPatterns->begin(); patt != newPatterns->end(); ++patt) {
0395 if (patt->getId() == maxPatt) {
0396 patt->printCode(maxCode);
0397 }
0398 }
0399
0400
0401 for (auto& q : {offsets, slopes, offsetuncs, slopeuncs, chi2ndfs, layers, legacypatterns}) {
0402 for (auto& p : q) {
0403 TCanvas* c = new TCanvas((p.second)->GetName(), (p.second)->GetTitle(), 800, 600);
0404 c->cd();
0405 gPad->SetLogy(1);
0406 (p.second)->Draw();
0407 c->SaveAs(TString::Format("figures_%dlayers/", N_LAYER_REQUIREMENT) + TString((p.second)->GetName()) + ".pdf");
0408 }
0409 }
0410
0411 file->Write();
0412 file->Close();
0413
0414 return 0;
0415 }
0416
0417
0418
0419 int searchForHits(const int code_hits[CSCConstants::NUM_LAYERS][CSCConstants::CLCT_PATTERN_WIDTH],
0420 int delta,
0421 unsigned N_LAYER_REQUIREMENT) {
0422 unsigned returnValue = 0;
0423 unsigned maxlayers = 0;
0424 for (unsigned iPat = 2; iPat < CSCPatternBank::clct_pattern_legacy_.size(); iPat++) {
0425 unsigned nlayers = 0;
0426 for (int i = 0; i < CSCConstants::NUM_LAYERS; i++) {
0427 bool layerhit = false;
0428 for (unsigned j = 0; j < CSCConstants::CLCT_PATTERN_WIDTH; j++) {
0429
0430 int jdelta = j + delta;
0431 if (jdelta < 0)
0432 jdelta = 0;
0433 if (jdelta >= CSCConstants::CLCT_PATTERN_WIDTH)
0434 jdelta = CSCConstants::CLCT_PATTERN_WIDTH - 1;
0435
0436
0437 if (CSCPatternBank::clct_pattern_legacy_[iPat][i][jdelta]) {
0438
0439 if (code_hits[i][j]) {
0440 layerhit = true;
0441 }
0442 }
0443 }
0444
0445 if (layerhit) {
0446 nlayers++;
0447 }
0448 }
0449 if (nlayers > maxlayers and nlayers >= N_LAYER_REQUIREMENT) {
0450 maxlayers = nlayers;
0451 returnValue = iPat;
0452 }
0453 }
0454 return returnValue;
0455 }
0456
0457
0458 int convertToLegacyPattern(const int code_hits[CSCConstants::NUM_LAYERS][CSCConstants::CLCT_PATTERN_WIDTH],
0459 unsigned N_LAYER_REQUIREMENT) {
0460 int returnValue = searchForHits(code_hits, 0, N_LAYER_REQUIREMENT);
0461
0462 if (!returnValue)
0463 returnValue = searchForHits(code_hits, -1, N_LAYER_REQUIREMENT);
0464
0465 if (!returnValue)
0466 returnValue = searchForHits(code_hits, 1, N_LAYER_REQUIREMENT);
0467 return returnValue;
0468 }
0469
0470 void getErrors(const vector<float>& x, const vector<float>& y, float& sigmaM, float& sigmaB) {
0471 int N = x.size();
0472 if (!N) {
0473 sigmaM = -9;
0474 sigmaB = -9;
0475 return;
0476 }
0477
0478 float sumx = 0;
0479 float sumx2 = 0;
0480 for (int i = 0; i < N; i++) {
0481 sumx += x[i];
0482 sumx2 += x[i] * x[i];
0483 }
0484
0485 float delta = N * sumx2 - sumx * sumx;
0486
0487 sigmaM = HALF_STRIP_ERROR * sqrt(1. * N / delta);
0488 sigmaB = HALF_STRIP_ERROR * sqrt(sumx2 / delta);
0489 }
0490
0491 void writeHeaderPosOffsetLUT(ofstream& file) {
0492 file << "#CSC Position Offset LUT\n"
0493 << "#Structure is comparator code (iCC) -> Stub position offset\n"
0494 << "#iCC bits: 12\n"
0495 << "#<header> v1.0 12 32 </header>\n";
0496 }
0497
0498 void writeHeaderSlopeLUT(ofstream& file) {
0499 file << "#CSC Slope LUT\n"
0500 << "#Structure is comparator code (iCC) -> Stub slope\n"
0501 << "#iCC bits: 12\n"
0502 << "#<header> v1.0 12 32 </header>\n";
0503 }
0504
0505 unsigned assignPosition(const float fvalue, const float fmin, const float fmax, const unsigned nbits) {
0506 bool debug = false;
0507 unsigned value = 0;
0508 const unsigned range = pow(2, nbits);
0509 const unsigned minValue = 0;
0510 const unsigned maxValue = range - 1;
0511 const double fdelta = (fmax - fmin) / range;
0512
0513 if (fvalue >= fmax) {
0514 value = maxValue;
0515 } else if (fvalue <= fmin) {
0516 value = minValue;
0517 } else {
0518 value = std::min(unsigned(std::ceil((fvalue - fmin) / fdelta)), maxValue);
0519 }
0520 if (debug)
0521 std::cout << "fvalue " << fvalue << " " << fmin << " " << fmax << " " << nbits << " " << value << std::endl;
0522
0523 return value;
0524 }
0525
0526 unsigned assignBending(const float fvalue) {
0527 bool debug = false;
0528 unsigned value = 0;
0529
0530
0531 float slopes[17] = {
0532 0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 2.0, 2.5};
0533
0534 for (unsigned i = 0; i < 16; i++) {
0535 if (fvalue >= slopes[i] and fvalue < slopes[i + 1]) {
0536 value = i;
0537 }
0538 }
0539
0540 if (fvalue >= slopes[16]) {
0541 value = 15;
0542 }
0543 return value;
0544 }
0545
0546 unsigned firmwareWord(const unsigned quality, const unsigned slope, const unsigned offset) {
0547
0548
0549
0550
0551
0552
0553 enum Masks { OffsetMask = 0xf, SlopeMask = 0x1f, QualityMask = 0x1ff };
0554 enum Shifts { OffsetShift = 14, SlopeShift = 9, QualityShift = 0 };
0555
0556 unsigned fwword = 0;
0557 setDataWord(fwword, quality, QualityShift, QualityMask);
0558 setDataWord(fwword, slope, SlopeShift, SlopeMask);
0559 setDataWord(fwword, offset, OffsetShift, OffsetMask);
0560
0561 return fwword;
0562 }
0563
0564 void setDataWord(unsigned& word, const unsigned newWord, const unsigned shift, const unsigned mask) {
0565
0566 word &= ~(mask << shift);
0567
0568
0569 word |= newWord << shift;
0570 }