EcalTrivialCondModules/src/ESTrivialConditionRetriever.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include <vector>
0004
0005 #include "CalibCalorimetry/EcalTrivialCondModules/interface/ESTrivialConditionRetriever.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008
0009 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0010
0011 //#include "DataFormats/Provenance/interface/Timestamp.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 using namespace edm;
0016
0017 ESTrivialConditionRetriever::ESTrivialConditionRetriever(const edm::ParameterSet& ps) {
0018   // initilize parameters used to produce cond DB objects
0019   adcToGeVLowConstant_ = ps.getUntrackedParameter<double>("adcToGeVLowConstant", 1.0);
0020   adcToGeVHighConstant_ = ps.getUntrackedParameter<double>("adcToGeVHighConstant", 1.0);
0021
0022   intercalibConstantMean_ = ps.getUntrackedParameter<double>("intercalibConstantMean", 1.0);
0023   intercalibConstantSigma_ = ps.getUntrackedParameter<double>("intercalibConstantSigma", 0.0);
0024
0025   //  intercalibErrorMean_ = ps.getUntrackedParameter<double>("IntercalibErrorMean",0.0);
0026
0027   ESpedMean_ = ps.getUntrackedParameter<double>("ESpedMean", 200.);
0028   ESpedRMS_ = ps.getUntrackedParameter<double>("ESpedRMS", 1.00);
0029
0030   getWeightsFromFile_ = ps.getUntrackedParameter<bool>("getWeightsFromFile", false);
0031
0032   std::string path = "CalibCalorimetry/EcalTrivialCondModules/data/";
0033   std::string weightType;
0034   std::ostringstream str;
0035
0036   weightType = str.str();
0037
0038   amplWeightsFile_ = ps.getUntrackedParameter<std::string>("amplWeightsFile", path + "ampWeightsES" + weightType);
0039
0040   // default weights for MGPA shape after pedestal subtraction
0041   getWeightsFromConfiguration(ps);
0042
0043   producedESPedestals_ = ps.getUntrackedParameter<bool>("producedESPedestals", true);
0044   producedESWeights_ = ps.getUntrackedParameter<bool>("producedESWeights", true);
0045
0046   producedESADCToGeVConstant_ = ps.getUntrackedParameter<bool>("producedESADCToGeVConstant", true);
0047
0048   verbose_ = ps.getUntrackedParameter<int>("verbose", 0);
0049
0050   //Tell Producer what we produce
0051   //setWhatproduce(this);
0052   if (producedESPedestals_)
0053     setWhatProduced(this, &ESTrivialConditionRetriever::produceESPedestals);
0054
0055   if (producedESWeights_) {
0056     setWhatProduced(this, &ESTrivialConditionRetriever::produceESWeightStripGroups);
0057     setWhatProduced(this, &ESTrivialConditionRetriever::produceESTBWeights);
0058   }
0059
0060   if (producedESADCToGeVConstant_)
0061     setWhatProduced(this, &ESTrivialConditionRetriever::produceESADCToGeVConstant);
0062
0063   // intercalibration constants
0064   producedESIntercalibConstants_ = ps.getUntrackedParameter<bool>("producedESIntercalibConstants", true);
0065   intercalibConstantsFile_ = ps.getUntrackedParameter<std::string>("intercalibConstantsFile", "");
0066
0067   if (producedESIntercalibConstants_) {      // user asks to produce constants
0068     if (intercalibConstantsFile_.empty()) {  // if file provided read constants
0069       //    setWhatProduced (this, &ESTrivialConditionRetriever::getIntercalibConstantsFromConfiguration ) ;
0070       //  } else { // set all constants to 1. or smear as specified by user
0071       setWhatProduced(this, &ESTrivialConditionRetriever::produceESIntercalibConstants);
0072     }
0073     findingRecord<ESIntercalibConstantsRcd>();
0074   }
0075
0076   // intercalibration constants errors
0077   /*  producedESIntercalibErrors_ = ps.getUntrackedParameter<bool>("producedESIntercalibErrors",true);
0078       intercalibErrorsFile_ = ps.getUntrackedParameter<std::string>("intercalibErrorsFile","") ;
0079
0080       if (producedESIntercalibErrors_) { // user asks to produce constants
0081       if(intercalibErrorsFile_ != "") {  // if file provided read constants
0082       setWhatProduced (this, &ESTrivialConditionRetriever::getIntercalibErrorsFromConfiguration ) ;
0083       } else { // set all constants to 1. or smear as specified by user
0084       setWhatProduced (this, &ESTrivialConditionRetriever::produceESIntercalibErrors ) ;
0085       }
0086       findingRecord<ESIntercalibErrorsRcd> () ;
0087       }
0088   */
0089
0090   // channel status
0091   producedESChannelStatus_ = ps.getUntrackedParameter<bool>("producedESChannelStatus", true);
0092   channelStatusFile_ = ps.getUntrackedParameter<std::string>("channelStatusFile", "");
0093
0094   if (producedESChannelStatus_) {
0095     if (!channelStatusFile_.empty()) {  // if file provided read channel map
0096       setWhatProduced(this, &ESTrivialConditionRetriever::getChannelStatusFromConfiguration);
0097     } else {  // set all channels to working -- FIXME might be changed
0098       setWhatProduced(this, &ESTrivialConditionRetriever::produceESChannelStatus);
0099     }
0100     findingRecord<ESChannelStatusRcd>();
0101   }
0102
0103   //Tell Finder what records we find
0104   if (producedESPedestals_)
0105     findingRecord<ESPedestalsRcd>();
0106
0107   if (producedESWeights_) {
0108     findingRecord<ESWeightStripGroupsRcd>();
0109     findingRecord<ESTBWeightsRcd>();
0110   }
0111
0112   if (producedESADCToGeVConstant_)
0113     findingRecord<ESADCToGeVConstantRcd>();
0114 }
0115
0116 ESTrivialConditionRetriever::~ESTrivialConditionRetriever() {}
0117
0118 //
0119 // member functions
0120 //
0121 void ESTrivialConditionRetriever::setIntervalFor(const edm::eventsetup::EventSetupRecordKey& rk,
0122                                                  const edm::IOVSyncValue& iTime,
0123                                                  edm::ValidityInterval& oValidity) {
0124   if (verbose_ >= 1)
0125     std::cout << "ESTrivialConditionRetriever::setIntervalFor(): record key = " << rk.name()
0126               << "\ttime: " << iTime.time().value() << std::endl;
0127   //For right now, we will just use an infinite interval of validity
0128   oValidity = edm::ValidityInterval(edm::IOVSyncValue::beginOfTime(), edm::IOVSyncValue::endOfTime());
0129 }
0130
0131 //produce methods
0132 std::unique_ptr<ESPedestals> ESTrivialConditionRetriever::produceESPedestals(const ESPedestalsRcd&) {
0133   std::cout << " producing pedestals" << std::endl;
0134   auto peds = std::make_unique<ESPedestals>();
0135   ESPedestals::Item ESitem;
0136   ESitem.mean = ESpedMean_;
0137   ESitem.rms = ESpedRMS_;
0138
0139   for (int istrip = ESDetId::ISTRIP_MIN; istrip <= ESDetId::ISTRIP_MAX; istrip++) {
0140     for (int ix = ESDetId::IX_MIN; ix <= ESDetId::IX_MAX; ix++) {
0141       for (int iy = ESDetId::IY_MIN; iy <= ESDetId::IY_MAX; iy++) {
0142         for (int iplane = 1; iplane <= 2; iplane++) {
0143           for (int izeta = -1; izeta <= 1; ++izeta) {
0144             if (izeta == 0)
0145               continue;
0146             try {
0147               //ESDetId Plane iplane Zside izeta
0148               ESDetId aPositiveId(istrip, ix, iy, iplane, izeta);
0149               peds->insert(std::make_pair(aPositiveId.rawId(), ESitem));
0150             } catch (cms::Exception& e) {
0151             }
0152           }
0153         }
0154       }
0155     }
0156   }
0157   //return std::unique_ptr<ESPedestals>( peds );
0158   std::cout << " produced pedestals" << std::endl;
0159   return peds;
0160 }
0161
0162 std::unique_ptr<ESWeightStripGroups> ESTrivialConditionRetriever::produceESWeightStripGroups(
0163     const ESWeightStripGroupsRcd&) {
0164   auto xtalGroups = std::make_unique<ESWeightStripGroups>();
0165   ESStripGroupId defaultGroupId(1);
0166   std::cout << "entering produce weight groups" << std::endl;
0167   for (int istrip = ESDetId::ISTRIP_MIN; istrip <= ESDetId::ISTRIP_MAX; istrip++) {
0168     for (int ix = ESDetId::IX_MIN; ix <= ESDetId::IX_MAX; ix++) {
0169       for (int iy = ESDetId::IY_MIN; iy <= ESDetId::IY_MAX; iy++) {
0170         for (int iplane = 1; iplane <= 2; iplane++) {
0171           for (int izeta = -1; izeta <= 1; ++izeta) {
0172             if (izeta == 0)
0173               continue;
0174             try {
0175               //ESDetId Plane iplane Zside izeta
0176               ESDetId anESId(istrip, ix, iy, iplane, izeta);
0177               //      xtalGroups->setValue(ebid.rawId(), ESStripGroupId(ieta) ); // define rings in eta
0178               xtalGroups->setValue(anESId.rawId(), defaultGroupId);  // define rings in eta
0179             } catch (cms::Exception& e) {
0180             }
0181           }
0182         }
0183       }
0184     }
0185   }
0186   std::cout << "done with produce weight groups" << std::endl;
0187
0188   return xtalGroups;
0189 }
0190
0191 std::unique_ptr<ESIntercalibConstants> ESTrivialConditionRetriever::produceESIntercalibConstants(
0192     const ESIntercalibConstantsRcd&) {
0193   auto ical = std::make_unique<ESIntercalibConstants>();
0194   std::cout << "entring produce intercalib " << std::endl;
0195
0196   for (int istrip = ESDetId::ISTRIP_MIN; istrip <= ESDetId::ISTRIP_MAX; istrip++) {
0197     for (int ix = ESDetId::IX_MIN; ix <= ESDetId::IX_MAX; ix++) {
0198       for (int iy = ESDetId::IY_MIN; iy <= ESDetId::IY_MAX; iy++) {
0199         for (int iplane = 1; iplane <= 2; iplane++) {
0200           for (int izeta = -1; izeta <= 1; ++izeta) {
0201             if (izeta == 0)
0202               continue;
0203             try {
0204               //ESDetId Plane iplane Zside izeta
0205               ESDetId anESId(istrip, ix, iy, iplane, izeta);
0206               double r = (double)std::rand() / (double(RAND_MAX) + double(1));
0207               ical->setValue(anESId.rawId(), intercalibConstantMean_ + r * intercalibConstantSigma_);
0208             } catch (cms::Exception& e) {
0209             }
0210           }
0211         }
0212       }
0213     }
0214   }
0215   std::cout << "done produce intercalib" << std::endl;
0216
0217   return ical;
0218 }
0219
0220 std::unique_ptr<ESADCToGeVConstant> ESTrivialConditionRetriever::produceESADCToGeVConstant(
0221     const ESADCToGeVConstantRcd&) {
0222   return std::make_unique<ESADCToGeVConstant>(adcToGeVLowConstant_, adcToGeVHighConstant_);
0223 }
0224
0225 std::unique_ptr<ESTBWeights> ESTrivialConditionRetriever::produceESTBWeights(const ESTBWeightsRcd&) {
0226   // create weights for the test-beam
0227   auto tbwgt = std::make_unique<ESTBWeights>();
0228
0229   int igrp = 1;
0230   ESWeightSet wgt = ESWeightSet(amplWeights_);
0231   //  ESWeightSet::ESWeightMatrix& mat1 = wgt.getWeights();
0232
0233   tbwgt->setValue(igrp, wgt);
0234
0235   return tbwgt;
0236 }
0237
0238 void ESTrivialConditionRetriever::getWeightsFromConfiguration(const edm::ParameterSet& ps) {
0239   ESWeightSet::ESWeightMatrix vampl;
0240
0241   if (!getWeightsFromFile_) {
0242     //      vampl.set(1.);
0243
0244     // amplwgtv[0]= ps.getUntrackedParameter< std::vector<double> >("amplWeights", vampl);
0245   } else if (getWeightsFromFile_) {
0246     edm::LogInfo("ESTrivialConditionRetriever")
0247         << "Reading amplitude weights from file " << edm::FileInPath(amplWeightsFile_).fullPath().c_str();
0248     std::ifstream amplFile(edm::FileInPath(amplWeightsFile_).fullPath().c_str());
0249     while (!amplFile.eof()) {
0250       for (int j = 0; j < 2; ++j) {
0251         std::vector<float> vec(3);
0252         for (int k = 0; k < 3; ++k) {
0253           float ww;
0254           amplFile >> ww;
0255           vec[k] = ww;
0256         }
0257         // vampl.putRow(vec);
0258       }
0259     }
0260   } else {
0261     //Not supported
0262     edm::LogError("ESTrivialConditionRetriever") << "Configuration not supported. Exception is raised ";
0263     throw cms::Exception("WrongConfig");
0264   }
0265
0266   amplWeights_ = ESWeightSet(vampl);
0267 }
0268
0269 // --------------------------------------------------------------------------------
0270
0271 std::unique_ptr<ESChannelStatus> ESTrivialConditionRetriever::getChannelStatusFromConfiguration(
0272     const ESChannelStatusRcd&) {
0273   auto ecalStatus = std::make_unique<ESChannelStatus>();
0274
0275   // start by setting all statuses to 0
0276
0277   for (int istrip = ESDetId::ISTRIP_MIN; istrip <= ESDetId::ISTRIP_MAX; istrip++) {
0278     for (int ix = ESDetId::IX_MIN; ix <= ESDetId::IX_MAX; ix++) {
0279       for (int iy = ESDetId::IY_MIN; iy <= ESDetId::IY_MAX; iy++) {
0280         for (int iplane = 1; iplane <= 2; iplane++) {
0281           for (int izeta = -1; izeta <= 1; ++izeta) {
0282             if (izeta == 0)
0283               continue;
0284             try {
0285               //ESDetId Plane iplane Zside izeta
0286               ESDetId anESId(istrip, ix, iy, iplane, izeta);
0287               //      xtalGroups->setValue(ebid.rawId(), ESStripGroupId(ieta) ); // define rings in eta
0288               ecalStatus->setValue(anESId, 0);
0289             } catch (cms::Exception& e) {
0290             }
0291           }
0292         }
0293       }
0294     }
0295   }
0296
0297   // overwrite the statuses which are in the file
0298
0299   edm::LogInfo("ESTrivialConditionRetriever")
0300       << "Reading channel statuses from file " << edm::FileInPath(channelStatusFile_).fullPath().c_str();
0301   std::ifstream statusFile(edm::FileInPath(channelStatusFile_).fullPath().c_str());
0302   if (!statusFile.good()) {
0303     edm::LogError("ESTrivialConditionRetriever") << "*** Problems opening file: " << channelStatusFile_;
0304     throw cms::Exception("Cannot open ECAL channel status file");
0305   }
0306
0307   std::string ESSubDet;
0308   std::string str;
0309   int hashIndex(0);
0310   int status(0);
0311
0312   while (!statusFile.eof()) {
0313     statusFile >> ESSubDet;
0314     if (ESSubDet != std::string("ES")) {
0315       std::getline(statusFile, str);
0316       continue;
0317     } else {
0318       statusFile >> hashIndex >> status;
0319     }
0320     // std::cout << ESSubDet << " " << hashIndex << " " << status;
0321
0322     if (ESSubDet == std::string("ES")) {
0323       ESDetId esid = ESDetId::unhashIndex(hashIndex);
0324       ecalStatus->setValue(esid, status);
0325     } else {
0326       edm::LogError("ESTrivialConditionRetriever") << " *** " << ESSubDet << " is not ES ";
0327     }
0328   }
0329   // the file is supposed to be in the form
0330   // ES hashed_index status
0331   // ES 132332 1  --> higher than 0  means bad
0332
0333   statusFile.close();
0334   return ecalStatus;
0335 }
0336
0337 std::unique_ptr<ESChannelStatus> ESTrivialConditionRetriever::produceESChannelStatus(const ESChannelStatusRcd&) {
0338   auto ical = std::make_unique<ESChannelStatus>();
0339   for (int istrip = ESDetId::ISTRIP_MIN; istrip <= ESDetId::ISTRIP_MAX; istrip++) {
0340     for (int ix = ESDetId::IX_MIN; ix <= ESDetId::IX_MAX; ix++) {
0341       for (int iy = ESDetId::IY_MIN; iy <= ESDetId::IY_MAX; iy++) {
0342         for (int iplane = 1; iplane <= 2; iplane++) {
0343           for (int izeta = -1; izeta <= 1; ++izeta) {
0344             if (izeta == 0)
0345               continue;
0346             try {
0347               //ESDetId Plane iplane Zside izeta
0348               ESDetId anESId(istrip, ix, iy, iplane, izeta);
0349               //      xtalGroups->setValue(ebid.rawId(), ESStripGroupId(ieta) ); // define rings in eta
0350               ical->setValue(anESId, 0);
0351             } catch (cms::Exception& e) {
0352             }
0353           }
0354         }
0355       }
0356     }
0357   }
0358
0359   return ical;
0360 }
0361
0362 // --------------------------------------------------------------------------------