EBPhase2TPGTools/plugins/EcalEBPhase2TPParamProducer.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ParameterSet/interface/FileInPath.h"
0006 //
0007 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0009 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0010 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0011 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0012 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
0013 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
0014 #include "CondFormats/EcalObjects/interface/EcalTPGCrystalStatus.h"
0015 #include "CondFormats/EcalObjects/interface/EcalLiteDTUPedestals.h"
0016 #include "CondFormats/DataRecord/interface/EcalLiteDTUPedestalsRcd.h"
0017 #include "DataFormats/EcalDigi/interface/EcalConstants.h"
0018 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0019 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0020
0021 #include <iostream>
0022 #include <string>
0023 #include <sstream>
0024 #include <vector>
0025 #include <ctime>
0026 #include <fstream>
0027 #include <iomanip>
0028 #include <TGraph.h>
0029 #include <TFile.h>
0030 #include <TMatrix.h>
0031 #include <zlib.h>
0032 /**
0033 \class EcalEBPhase2TPParamProducer
0034 \author L. Lutton, N. Marinelli - Univ. of Notre Dame
0035 \brief TPG Param Builder for Phase2
0036 */
0037
0038 class EcalEBPhase2TPParamProducer : public edm::one::EDAnalyzer<> {
0039 public:
0040   explicit EcalEBPhase2TPParamProducer(edm::ParameterSet const& pSet);
0041   ~EcalEBPhase2TPParamProducer() override;
0042   void analyze(const edm::Event& evt, const edm::EventSetup& evtSetup) override;
0043   void beginJob() override;
0044   static void fillDescriptions(edm::ConfigurationDescriptions&);
0045
0046 private:
0047   std::vector<int> computeWeights(int type);
0048
0049   void getNumericalDeriv(TGraph graph, TGraph& deriv);
0050   void fillFMat(std::vector<unsigned int> clockSampleSet,
0051                 bool useThirdPulse,
0052                 std::vector<float> sampleSet,
0053                 std::vector<float> sampleDotSet,
0054                 TMatrix& FMat,
0055                 unsigned int binOfMaximum);
0056   void getGMatrix(TMatrix FMat, float scaleMatrixBy, TMatrix& GMat);
0057   void getPulseSampleSet(TGraph pulseGraph, float phaseShift, std::vector<float>& sampleSet);
0058   bool computeLinearizerParam(double theta, double gainRatio, double calibCoeff, int& shift, int& mult);
0059
0060   const edm::ESGetToken<CaloSubdetectorGeometry, EcalBarrelGeometryRecord> theBarrelGeometryToken_;
0061   const edm::FileInPath inFile_;
0062   const std::string outFile_;
0063   const int nSamplesToUse_;
0064   const bool useBXPlusOne_;
0065   const double phaseShift_;
0066   const unsigned int nWeightGroups_;
0067   const edm::ESGetToken<EcalLiteDTUPedestalsMap, EcalLiteDTUPedestalsRcd> theEcalTPGPedestals_Token_;
0068
0069   gzFile out_file_;
0070   TGraph* thePulse_;
0071   TGraph* pulseDot_;
0072
0073   const UInt_t NPoints_ = 1599;  //With the CMSSW pulse
0074
0075   static constexpr float norm_ = 1 / 503.109;  // with the CMSSW pulse shape
0076   static constexpr float offset_ = 0.;         // with the CMSSW pulse shape
0077   int multToInt_ = 0x1000;
0078
0079   int i2cSub_[2] = {0, 0};
0080
0081   const double et_sat_;
0082   const double xtal_LSB_;
0083   const unsigned int binOfMaximum_;
0084   static const int linTopRange_;
0085 };
0086
0087 EcalEBPhase2TPParamProducer::EcalEBPhase2TPParamProducer(edm::ParameterSet const& pSet)
0088     : theBarrelGeometryToken_(esConsumes(edm::ESInputTag("", "EcalBarrel"))),
0089       inFile_(pSet.getParameter<edm::FileInPath>("inputFile")),
0090       outFile_(pSet.getUntrackedParameter<std::string>("outputFile")),
0091       nSamplesToUse_(pSet.getParameter<unsigned int>("nSamplesToUse")),
0092       useBXPlusOne_(pSet.getParameter<bool>("useBXPlusOne")),
0093       phaseShift_(pSet.getParameter<double>("phaseShift")),
0094       nWeightGroups_(pSet.getParameter<unsigned int>("nWeightGroups")),
0095       theEcalTPGPedestals_Token_(esConsumes(edm::ESInputTag("EcalLiteDTUPedestals", ""))),
0096       et_sat_(pSet.getParameter<double>("Et_sat")),
0097       xtal_LSB_(pSet.getParameter<double>("xtal_LSB")),
0098       binOfMaximum_(pSet.getParameter<unsigned int>("binOfMaximum"))
0099
0100 {
0101   out_file_ = gzopen(outFile_.c_str(), "wb");
0102
0103   std::string filename = inFile_.fullPath();
0104   TFile* inFile = new TFile(filename.c_str(), "READ");
0105
0106   inFile->GetObject("average-pulse", thePulse_);
0107   delete inFile;
0108
0109   if (binOfMaximum_ != 6 && binOfMaximum_ != 8)
0110     edm::LogError("EcalEBPhase2TPParamProducer")
0111         << " Value for binOfMaximum " << binOfMaximum_ << " is wrong, The default binOfMaximum=6  will be used";
0112
0113   if (nSamplesToUse_ != 6 && nSamplesToUse_ != 8 && nSamplesToUse_ != 12)
0114     edm::LogError("EcalEBPhase2TPParamProducer")
0115         << " Value for nSamplesToUse " << nSamplesToUse_ << " is wrong, The default nSamplesToUse=8 will be used";
0116 }
0117
0118 EcalEBPhase2TPParamProducer::~EcalEBPhase2TPParamProducer() { gzclose(out_file_); }
0119
0120 void EcalEBPhase2TPParamProducer::beginJob() {}
0121
0122 void EcalEBPhase2TPParamProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0123   edm::ParameterSetDescription desc;
0124   desc.add<edm::FileInPath>("inputFile");
0125   desc.addUntracked<std::string>("outputFile");
0126   desc.add<unsigned int>("nSamplesToUse", 8);
0127   desc.add<bool>("useBXPlusOne", false);
0128   desc.add<double>("phaseShift", 2.581);
0129   desc.add<unsigned int>("nWeightGroups", 61200);
0130   desc.add<double>("Et_sat", 1998.36);
0131   desc.add<double>("xtal_LSB", 0.0488);
0132   desc.add<unsigned int>("binOfMaximum", 6);
0133   descriptions.add("ecalEBPhase2TPParamProducerDefault", desc);
0134 }
0135
0136 void EcalEBPhase2TPParamProducer::analyze(const edm::Event& evt, const edm::EventSetup& evtSetup) {
0137   using namespace edm;
0138   using namespace std;
0139
0140   const EcalLiteDTUPedestals* peds = nullptr;
0141   const auto* theBarrelGeometry = &evtSetup.getData(theBarrelGeometryToken_);
0142   const auto* theEcalTPPedestals = &evtSetup.getData(theEcalTPGPedestals_Token_);
0143
0144   std::string tmpStringConv;
0145   const char* tmpStringOut;
0146
0147   // Compute weights  //
0148   std::vector<int> ampWeights[nWeightGroups_];
0149   std::vector<int> timeWeights[nWeightGroups_];
0150
0151   for (unsigned int iGr = 0; iGr < nWeightGroups_; iGr++) {
0152     ampWeights[iGr] = computeWeights(1);
0153     timeWeights[iGr] = computeWeights(2);
0154   }
0155
0156   /* write to compressed file  */
0157   std::stringstream toCompressStream("");
0158   for (unsigned int iGr = 0; iGr < nWeightGroups_; iGr++) {
0159     toCompressStream << "  WEIGHTAMP  " << dec << iGr << std::endl;
0160     for (long unsigned int i = 0; i < ampWeights[iGr].size(); i++) {
0161       if (ampWeights[iGr][i] < 0)
0162         toCompressStream << "-0x" << std::hex << abs(ampWeights[iGr][i]) << " ";
0163       else
0164         toCompressStream << "0x" << std::hex << ampWeights[iGr][i] << " ";
0165     }
0166     toCompressStream << "\n";
0167   }
0168   toCompressStream << "\n";
0169   tmpStringConv = toCompressStream.str();
0170   tmpStringOut = tmpStringConv.c_str();
0171   gzwrite(out_file_, tmpStringOut, std::strlen(tmpStringOut));
0172   toCompressStream.str(std::string());
0173
0174   for (unsigned int iGr = 0; iGr < nWeightGroups_; iGr++) {
0175     toCompressStream << "WEIGHTTIME " << dec << iGr << std::endl;
0176     for (long unsigned int i = 0; i < timeWeights[iGr].size(); i++) {
0177       if (timeWeights[iGr][i] < 0)
0178         toCompressStream << "-0x" << std::hex << abs(timeWeights[iGr][i]) << " ";
0179       else
0180         toCompressStream << "0x" << std::hex << timeWeights[iGr][i] << " ";
0181     }
0182     toCompressStream << "\n";
0183   }
0184
0185   toCompressStream << "\n";
0186   tmpStringConv = toCompressStream.str();
0187   tmpStringOut = tmpStringConv.c_str();
0188   gzwrite(out_file_, tmpStringOut, std::strlen(tmpStringOut));
0189   toCompressStream.str(std::string());
0190
0191   //  fill  map between xTals and groups. If each xTal is a group there is a one-to-one map
0192   const std::vector<DetId>& ebCells = theBarrelGeometry->getValidDetIds(DetId::Ecal, EcalBarrel);
0193   std::map<int, int> mapXtalToGroup;
0194
0195   int iGroup = 0;
0196   for (const auto& it : ebCells) {
0197     EBDetId id(it);
0198     std::pair<int, int> xTalToGroup(id.rawId(), iGroup);
0199     mapXtalToGroup.insert(xTalToGroup);
0200     iGroup++;
0201   }
0202
0203   //write to file
0204
0205   for (std::map<int, int>::const_iterator it = mapXtalToGroup.begin(); it != mapXtalToGroup.end(); it++) {
0206     toCompressStream << "CRYSTAL " << dec << it->first << std::endl;
0207     toCompressStream << it->second << std::endl;
0208   }
0209   tmpStringConv = toCompressStream.str();
0210   tmpStringOut = tmpStringConv.c_str();
0211   gzwrite(out_file_, tmpStringOut, std::strlen(tmpStringOut));
0212   toCompressStream.str(std::string());
0213
0214   /////////////////////////////////////
0215
0216   for (const auto& it : ebCells) {
0217     EBDetId id(it);
0218     toCompressStream << "LINCONST " << dec << id.rawId() << std::endl;
0219     double theta = theBarrelGeometry->getGeometry(id)->getPosition().theta();
0220     EcalLiteDTUPedestalsMap::const_iterator itped = theEcalTPPedestals->getMap().find(id);
0221
0222     if (itped != theEcalTPPedestals->end()) {
0223       peds = &(*itped);
0224
0225     } else {
0226       edm::LogError("EcalEBPhase2TPParamProducer") << " could not find EcalLiteDTUPedestal entry for " << id;
0227       throw cms::Exception("could not find pedestals");
0228     }
0229
0230     int shift, mult;
0231     double calibCoeff = 1.;
0232     bool ok;
0233     for (unsigned int i = 0; i < ecalPh2::NGAINS; ++i) {
0234       ok = computeLinearizerParam(theta, ecalph2::gains[ecalPh2::NGAINS - 1 - i], calibCoeff, shift, mult);
0235       if (!ok) {
0236         edm::LogError("EcalEBPhase2TPParamProducer")
0237             << "unable to compute the parameters for SM=" << id.ism() << " xt=" << id.ic() << " " << id.rawId();
0238         throw cms::Exception("unable to compute the parameters");
0239
0240       } else {
0241         int tmpPedByGain = (int)(peds->mean(i) + 0.5);
0242         toCompressStream << std::hex << " 0x" << tmpPedByGain << " 0x" << mult << " 0x" << shift << " " << i2cSub_[i]
0243                          << std::endl;
0244       }
0245     }
0246   }
0247   tmpStringConv = toCompressStream.str();
0248   tmpStringOut = tmpStringConv.c_str();
0249   gzwrite(out_file_, tmpStringOut, std::strlen(tmpStringOut));
0250   toCompressStream.str(std::string());
0251 }
0252
0253 std::vector<int> EcalEBPhase2TPParamProducer::computeWeights(int type) {
0254   std::vector<float> sampleSet;
0255   std::vector<float> sampleDotSet;
0256   std::vector<unsigned int> clockSampleSet;
0257   double scaleMatrixBy = 1.;
0258   int lbinOfMaximum = binOfMaximum_;
0259
0260   switch (binOfMaximum_) {
0261     case 6:
0262       break;
0263     case 8:
0264       break;
0265     default:
0266       lbinOfMaximum = 6;
0267       break;
0268   }
0269
0270   switch (nSamplesToUse_) {
0271     case 12:
0272       clockSampleSet = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
0273       break;
0274     case 8:
0275       switch (lbinOfMaximum) {
0276         case 8:
0277           clockSampleSet = {2, 3, 4, 5, 6, 7, 8, 9};
0278           break;
0279         case 6:
0280           clockSampleSet = {0, 1, 2, 3, 4, 5, 6, 7};
0281           break;
0282       }
0283       break;
0284
0285     case 6:
0286       switch (lbinOfMaximum) {
0287         case 8:
0288           clockSampleSet = {3, 4, 6, 7, 8, 9};
0289           break;
0290         case 6:
0291           clockSampleSet = {1, 2, 4, 5, 6, 7};
0292           break;
0293       }
0294       break;
0295
0296     default:
0297       clockSampleSet = {0, 1, 2, 3, 4, 5, 6, 7};
0298       break;
0299   }
0300
0301   getPulseSampleSet(*thePulse_, phaseShift_, sampleSet);
0302   pulseDot_ = new TGraph();
0303   getNumericalDeriv(*thePulse_, *pulseDot_);
0304   getPulseSampleSet(*pulseDot_, phaseShift_, sampleDotSet);
0305
0306   unsigned int fMatColumns = useBXPlusOne_ ? 6 : 4;
0307
0308   TMatrix fMat(clockSampleSet.size(), fMatColumns);
0309   fillFMat(clockSampleSet, useBXPlusOne_, sampleSet, sampleDotSet, fMat, lbinOfMaximum);
0310   TMatrix gMat(fMatColumns, clockSampleSet.size());
0311
0312   getGMatrix(fMat, scaleMatrixBy, gMat);
0313
0314   std::vector<int> tmpWeightVec;
0315   std::vector<int> tmpTimeWeightVec;
0316   unsigned int iClock = 0;
0317   for (unsigned int iSample = 0; iSample < 12; iSample++) {
0318     bool inSampleSet = false;
0319     for (unsigned int clockSample = 0; clockSample < clockSampleSet.size(); clockSample++) {
0320       if (iSample == clockSampleSet[clockSample]) {
0321         inSampleSet = true;
0322         iClock = clockSample;
0323         break;
0324       }
0325     }
0326     if (inSampleSet) {
0327       if (type == 1)
0328         tmpWeightVec.push_back(round(gMat(2, iClock) * multToInt_));  // amp weights
0329       if (type == 2)
0330         tmpWeightVec.push_back(round(gMat(3, iClock) * multToInt_));  // time weights
0331     } else {
0332       if (type == 1)
0333         tmpWeightVec.push_back(0);  // amp weights
0334       if (type == 2)
0335         tmpWeightVec.push_back(0);  // time weights
0336     }
0337   }
0338
0339   return tmpWeightVec;
0340 }
0341
0342 void EcalEBPhase2TPParamProducer::getNumericalDeriv(TGraph graph, TGraph& deriv) {
0343   UInt_t numPoints = graph.GetN();
0344   if (numPoints != NPoints_) {
0345     edm::LogWarning("EcalEBPhase2TPParamProducer") << "Error! Wrong amount of points in pulse graph! ";
0346   }
0347   Double_t xval;
0348   Double_t yval;
0349   Double_t xvalPOne;
0350   Double_t yvalPOne;
0351
0352   for (UInt_t p = 0; p < NPoints_ - 1; p++) {
0353     graph.GetPoint(p, xval, yval);
0354     graph.GetPoint(p + 1, xvalPOne, yvalPOne);
0355     float midpoint = (xvalPOne + xval) / 2;
0356     float rise = yvalPOne - yval;
0357     float run = xvalPOne - xval;
0358     deriv.SetPoint(deriv.GetN(), midpoint, rise / run);
0359   }
0360   deriv.SetName("pulse_prime");
0361 }
0362
0363 void EcalEBPhase2TPParamProducer::fillFMat(std::vector<UInt_t> clockSampleSet,
0364                                            bool useThirdPulse,
0365                                            std::vector<float> sampleSet,
0366                                            std::vector<float> sampleDotSet,
0367                                            TMatrix& fMat,
0368                                            uint binOfMaximum) {
0369   Int_t iShift = 8 - binOfMaximum;
0370   for (UInt_t i = 0; i < clockSampleSet.size(); i++) {
0371     Int_t tmpClockToSample = clockSampleSet[i] + iShift;
0372     fMat(i, 0) = sampleSet[tmpClockToSample];
0373     fMat(i, 1) = sampleDotSet[tmpClockToSample];
0374     if (tmpClockToSample > 4) {
0375       fMat(i, 2) = sampleSet[tmpClockToSample - 4];
0376       fMat(i, 3) = sampleDotSet[tmpClockToSample - 4];
0377     }
0378     if (clockSampleSet[i] > 8 && useThirdPulse) {
0379       fMat(i, 4) = sampleSet[tmpClockToSample - 8];
0380       fMat(i, 5) = sampleDotSet[tmpClockToSample - 8];
0381     }
0382   }
0383 }
0384
0385 void EcalEBPhase2TPParamProducer::getGMatrix(TMatrix fMat, float scaleMatrixBy, TMatrix& gMat) {
0386   TMatrix FT = fMat;
0387   FT.T();
0388   TMatrix tmpFT = FT;
0389   TMatrix FTDotF = TMatrix(tmpFT, TMatrix::kMult, fMat);
0390   TMatrix InvFTDotF = FTDotF;
0391
0392   //Possible for this bit to fail depending on the sample set and phase shift
0393   InvFTDotF.Invert();
0394
0395   TMatrix tmpMat(InvFTDotF, TMatrix::kMult, FT);
0396   gMat = tmpMat;
0397   gMat *= scaleMatrixBy;
0398 }
0399
0400 void EcalEBPhase2TPParamProducer::getPulseSampleSet(TGraph pulseGraph,
0401                                                     float phaseShift,
0402                                                     std::vector<float>& sampleSet) {
0403   for (UInt_t i = 0; i < ecalPh2::sampleSize; i++) {
0404     float t = (ecalPh2::Samp_Period * i) + phaseShift;
0405     float y = pulseGraph.Eval(t + offset_) * norm_;
0406     sampleSet.push_back(y);
0407   }
0408 }
0409
0410 bool EcalEBPhase2TPParamProducer::computeLinearizerParam(
0411     double theta, double gainRatio, double calibCoeff, int& shift, int& mult) {
0412   bool result = false;
0413
0414   static constexpr double linTopRange_ = 16383.;
0415   //  linTopRange_ 16383 = (2**14)-1  is setting the top of the range for the linearizer output
0416   double factor = (linTopRange_ * (xtal_LSB_ * gainRatio * calibCoeff * sin(theta))) / et_sat_;
0417   //first with shift_ = 0
0418   //add 0.5 (for rounding) and set to int
0419   //Here we are getting mult with a max bit length of 8
0420   //and shift_ with a max bit length of 4
0421   mult = (int)(factor + 0.5);
0422   for (shift = 0; shift < 15; shift++) {
0423     if (mult >= 128 && mult < 256) {
0424       result = true;
0425       break;
0426     }
0427     factor *= 2;
0428     mult = (int)(factor + 0.5);
0429   }
0430
0431   return result;
0432 }
0433
0434 // DEfine this module as a plug-in
0435 DEFINE_FWK_MODULE(EcalEBPhase2TPParamProducer);