EcalTrivialCondModules/src/EcalTrivialObjectAnalyzer.cc

0001 //
0002 // Created: 2 Mar 2006
0003 //          Shahram Rahatlou, University of Rome & INFN
0004 //
0005 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialObjectAnalyzer.h"
0006
0007 #include <stdexcept>
0008 #include <string>
0009 #include <iostream>
0010 #include <iomanip>
0011 #include <map>
0012 #include <vector>
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "CondFormats/EcalObjects/interface/EcalXtalGroupId.h"
0016 #include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
0017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0018 #include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
0019
0020 #include "CLHEP/Matrix/Matrix.h"
0021
0022 using namespace std;
0023
0024 EcalTrivialObjectAnalyzer::EcalTrivialObjectAnalyzer(edm::ParameterSet const& p)
0025     : pedestalsToken_(esConsumes()),
0026       adcToGevConstantToken_(esConsumes()),
0027       weightXtalGroupsToken_(esConsumes()),
0028       gainRatiosToken_(esConsumes()),
0029       intercalibConstantsToken_(esConsumes()),
0030       intercalibErrorsToken_(esConsumes()),
0031       timeCalibConstantsToken_(esConsumes()),
0032       timeCalibErrorsToken_(esConsumes()),
0033       timeOffsetConstantToken_(esConsumes()),
0034       tbWeightsToken_(esConsumes()),
0035       clusterLocalContCorrToken_(esConsumes()),
0036       clusterCrackCorrToken_(esConsumes()),
0037       clusterEnergyCorrectionToken_(esConsumes()),
0038       clusterEnergyUncertaintyToken_(esConsumes()),
0039       clusterEnergyCorrectionObjectSpecificToken_(esConsumes()),
0040       laserAlphasToken_(esConsumes()),
0041       laserAPDPNRatiosRefToken_(esConsumes()),
0042       laserAPDPNRatiosToken_(esConsumes()),
0043       channelStatusToken_(esConsumes()),
0044       sampleMaskToken_(esConsumes()) {}
0045
0046 void EcalTrivialObjectAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& context) {
0047   using namespace edm::eventsetup;
0048   // Context is not used.
0049   edm::LogInfo(">>> EcalTrivialObjectAnalyzer: processing run ") << e.id().run() << " event: " << e.id().event();
0050
0051   // ADC -> GeV Scale
0052   const auto& agc = context.getData(adcToGevConstantToken_);
0053   edm::LogInfo("Global ADC->GeV scale: EB ") << std::setprecision(6) << agc.getEBValue()
0054                                              << " GeV/ADC count"
0055                                                 " EE "
0056                                              << std::setprecision(6) << agc.getEEValue() << " GeV/ADC count";
0057
0058   // use a channel to fetch values from DB
0059   double r1 = (double)std::rand() / (double(RAND_MAX) + double(1));
0060   int ieta = int(1 + r1 * 85);
0061   r1 = (double)std::rand() / (double(RAND_MAX) + double(1));
0062   int iphi = int(1 + r1 * 20);
0063
0064   EBDetId ebid(ieta, iphi);  //eta,phi
0065   edm::LogInfo("EcalTrivialObjectAnalyzer: using EBDetId: ") << ebid;
0066
0067   const auto& myped = context.getData(pedestalsToken_);
0068   EcalPedestals::const_iterator it = myped.find(ebid.rawId());
0069   if (it != myped.end()) {
0070     edm::LogInfo("EcalPedestal: ") << "  mean_x1:  " << std::setprecision(8) << (*it).mean_x1
0071                                    << " rms_x1: " << (*it).rms_x1 << "  mean_x6:  " << (*it).mean_x6
0072                                    << " rms_x6: " << (*it).rms_x6 << "  mean_x12: " << (*it).mean_x12
0073                                    << " rms_x12: " << (*it).rms_x12;
0074   } else {
0075     edm::LogInfo("No pedestal found for this xtal! something wrong with EcalPedestals in your DB? ");
0076   }
0077
0078   // fetch map of groups of xtals
0079   const auto& grp = context.getData(weightXtalGroupsToken_);
0080
0081   EcalXtalGroupsMap::const_iterator git = grp.getMap().find(ebid.rawId());
0082   EcalXtalGroupId gid;
0083   if (git != grp.getMap().end()) {
0084     edm::LogInfo("XtalGroupId.id() = ") << std::setprecision(3) << (*git).id();
0085     gid = (*git);
0086   } else {
0087     edm::LogInfo("No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?");
0088   }
0089
0090   // Gain Ratios
0091   const auto& gr = context.getData(gainRatiosToken_);
0092
0093   EcalGainRatioMap::const_iterator grit = gr.getMap().find(ebid.rawId());
0094   EcalMGPAGainRatio mgpa;
0095   if (grit != gr.getMap().end()) {
0096     mgpa = (*grit);
0097
0098     edm::LogInfo("EcalMGPAGainRatio: ") << "gain 12/6 :  " << std::setprecision(4) << mgpa.gain12Over6()
0099                                         << " gain 6/1: " << mgpa.gain6Over1();
0100   } else {
0101     edm::LogInfo("No MGPA Gain Ratio found for this xtal! something wrong with EcalGainRatios in your DB? ");
0102   }
0103
0104   // Intercalib constants
0105   const auto& ical = context.getData(intercalibConstantsToken_);
0106
0107   EcalIntercalibConstantMap::const_iterator icalit = ical.getMap().find(ebid.rawId());
0108   EcalIntercalibConstant icalconst;
0109   if (icalit != ical.getMap().end()) {
0110     icalconst = (*icalit);
0111
0112     edm::LogInfo("EcalIntercalibConstant: ") << std::setprecision(6) << icalconst;
0113   } else {
0114     edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalIntercalibConstants in your DB? ");
0115   }
0116
0117   // Intercalib errors
0118   const auto& icalErr = context.getData(intercalibErrorsToken_);
0119
0120   EcalIntercalibErrorMap::const_iterator icalitErr = icalErr.getMap().find(ebid.rawId());
0121   EcalIntercalibError icalconstErr;
0122   if (icalitErr != icalErr.getMap().end()) {
0123     icalconstErr = (*icalitErr);
0124
0125     edm::LogInfo("EcalIntercalibError: ") << std::setprecision(6) << icalconstErr;
0126   } else {
0127     edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalIntercalibErrors in your DB? ");
0128   }
0129
0130   {  // quick and dirty for cut and paste ;) it is a test program...
0131     // TimeCalib constants
0132     const auto& ical = context.getData(timeCalibConstantsToken_);
0133
0134     EcalTimeCalibConstantMap::const_iterator icalit = ical.getMap().find(ebid.rawId());
0135     EcalTimeCalibConstant icalconst;
0136     if (icalit != ical.getMap().end()) {
0137       icalconst = (*icalit);
0138
0139       edm::LogInfo("EcalTimeCalibConstant: ") << std::setprecision(6) << icalconst;
0140     } else {
0141       edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalTimeCalibConstants in your DB? ");
0142     }
0143
0144     // TimeCalib errors
0145     const auto& icalErr = context.getData(timeCalibErrorsToken_);
0146
0147     EcalTimeCalibErrorMap::const_iterator icalitErr = icalErr.getMap().find(ebid.rawId());
0148     EcalTimeCalibError icalconstErr;
0149     if (icalitErr != icalErr.getMap().end()) {
0150       icalconstErr = (*icalitErr);
0151
0152       edm::LogInfo("EcalTimeCalibError: ") << std::setprecision(6) << icalconstErr;
0153     } else {
0154       edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalTimeCalibErrors in your DB? ");
0155     }
0156   }
0157
0158   // fetch Time Offset
0159   //std::cout <<"Fetching TimeOffsetConstant from DB " << std::endl;
0160
0161   // Time Offset constants
0162   const auto& TOff = context.getData(timeOffsetConstantToken_);
0163
0164   edm::LogInfo("EcalTimeOffsetConstant: ") << " EB " << TOff.getEBValue() << " EE " << TOff.getEEValue();
0165
0166   // fetch TB weights
0167   edm::LogInfo("Fetching EcalTBWeights from DB ");
0168   const auto& wgts = context.getData(tbWeightsToken_);
0169   edm::LogInfo("EcalTBWeightMap.size(): ") << std::setprecision(3) << wgts.getMap().size();
0170
0171   // look up the correct weights for this  xtal
0172   //EcalXtalGroupId gid( (*git) );
0173   EcalTBWeights::EcalTDCId tdcid(1);
0174
0175   edm::LogInfo("Lookup EcalWeightSet for groupid: ") << std::setprecision(3) << gid.id() << " and TDC id " << tdcid;
0176   EcalTBWeights::EcalTBWeightMap::const_iterator wit = wgts.getMap().find(std::make_pair(gid, tdcid));
0177   EcalWeightSet wset;
0178   if (wit != wgts.getMap().end()) {
0179     wset = wit->second;
0180     edm::LogInfo("check size of data members in EcalWeightSet");
0181     //wit->second.print(std::cout);
0182
0183     //typedef std::vector< std::vector<EcalWeight> > EcalWeightMatrix;
0184     const EcalWeightSet::EcalWeightMatrix& mat1 = wit->second.getWeightsBeforeGainSwitch();
0185     const EcalWeightSet::EcalWeightMatrix& mat2 = wit->second.getWeightsAfterGainSwitch();
0186
0187     //       std::cout << "WeightsBeforeGainSwitch.size: " << mat1.size()
0188     //                 << ", WeightsAfterGainSwitch.size: " << mat2.size() << std::endl;
0189
0190     CLHEP::HepMatrix clmat1(3, 10, 0);
0191     CLHEP::HepMatrix clmat2(3, 10, 0);
0192     for (int irow = 0; irow < 3; irow++) {
0193       for (int icol = 0; icol < 10; icol++) {
0194         clmat1[irow][icol] = mat1(irow, icol);
0195         clmat2[irow][icol] = mat2(irow, icol);
0196       }
0197     }
0198     edm::LogInfo("weight matrix before gain switch:") << clmat1;
0199     edm::LogInfo("weight matrix after gain switch:") << clmat2;
0200
0201   } else {
0202     edm::LogInfo("No weights found for EcalGroupId: ") << gid.id() << " and  EcalTDCId: " << tdcid;
0203   }
0204
0205   // cluster functions/corrections
0206   const auto& paramLocalCont = context.getData(clusterLocalContCorrToken_);
0207   edm::LogInfo("LocalContCorrParameters:");
0208   for (EcalFunctionParameters::const_iterator it = paramLocalCont.params().begin(); it != paramLocalCont.params().end();
0209        ++it) {
0210     // edm::LogInfo(" ") << *it;
0211   }
0212   // std::cout << "\n";
0213   const auto& paramCrack = context.getData(clusterCrackCorrToken_);
0214   edm::LogInfo("CrackCorrParameters:");
0215   for (EcalFunctionParameters::const_iterator it = paramCrack.params().begin(); it != paramCrack.params().end(); ++it) {
0216     //  edm::LogInfo(" ") << *it;
0217   }
0218   // std::cout << "\n";
0219   const auto& paramEnergyCorrection = context.getData(clusterEnergyCorrectionToken_);
0220   edm::LogInfo("EnergyCorrectionParameters:");
0221   for (EcalFunctionParameters::const_iterator it = paramEnergyCorrection.params().begin();
0222        it != paramEnergyCorrection.params().end();
0223        ++it) {
0224     // edm::LogInfo(" ") << *it;
0225   }
0226   // std::cout << "\n";
0227   const auto& paramEnergyUncertainty = context.getData(clusterEnergyUncertaintyToken_);
0228   edm::LogInfo("EnergyCorrectionParameters:");
0229   for (EcalFunctionParameters::const_iterator it = paramEnergyUncertainty.params().begin();
0230        it != paramEnergyUncertainty.params().end();
0231        ++it) {
0232     // edm::LogInfo(" ") << *it;
0233   }
0234   // std::cout << "\n";
0235   const auto& paramEnergyCorrectionObjectSpecific = context.getData(clusterEnergyCorrectionObjectSpecificToken_);
0236   edm::LogInfo("EnergyCorrectionObjectSpecificParameters:");
0237   for (EcalFunctionParameters::const_iterator it = paramEnergyCorrectionObjectSpecific.params().begin();
0238        it != paramEnergyCorrectionObjectSpecific.params().end();
0239        ++it) {
0240     //  edm::LogInfo(" ") << *it;
0241   }
0242   // std::cout << "\n";
0243
0244   // laser correction
0245
0246   // laser alphas
0247   const auto& lalpha = context.getData(laserAlphasToken_);
0248
0249   EcalLaserAlphaMap::const_iterator lalphait;
0250   lalphait = lalpha.getMap().find(ebid.rawId());
0251   if (lalphait != lalpha.getMap().end()) {
0252     edm::LogInfo("EcalLaserAlpha: ") << std::setprecision(6) << (*lalphait);
0253   } else {
0254     edm::LogInfo("No laser alpha found for this xtal! something wrong with EcalLaserAlphas in your DB? ");
0255   }
0256
0257   // laser apdpnref
0258   const auto& lref = context.getData(laserAPDPNRatiosRefToken_);
0259
0260   EcalLaserAPDPNRatiosRef::const_iterator lrefit;
0261   lrefit = lref.getMap().find(ebid.rawId());
0262   if (lrefit != lref.getMap().end()) {
0263     //  edm::LogInfo("EcalLaserAPDPNRatiosRef: ")
0264     //         <<std::setprecision(6)
0265     //         << (*lrefit)
0266     ;
0267   } else {
0268     edm::LogInfo("No laser apd/pn ref found for this xtal! something wrong with EcalLaserAPDPNRatiosRef in your DB? ");
0269   }
0270
0271   // laser apdpn ratios
0272   const auto& lratio = context.getData(laserAPDPNRatiosToken_);
0273
0274   EcalLaserAPDPNRatios::EcalLaserAPDPNRatiosMap::const_iterator lratioit;
0275   lratioit = lratio.getLaserMap().find(ebid.rawId());
0276   EcalLaserAPDPNRatios::EcalLaserAPDPNpair lratioconst;
0277
0278   if (lratioit != lratio.getLaserMap().end()) {
0279     lratioconst = (*lratioit);
0280
0281     //  edm::LogInfo("EcalLaserAPDPNRatios: ")
0282     //         <<e.id().run() << " " << e.id().event() << " "
0283     //         << std::setprecision(6)
0284     //         << lratioconst.p1 << " " << lratioconst.p2
0285     ;
0286   } else {
0287     edm::LogInfo("No laser apd/pn ratio found for this xtal! something wrong with EcalLaserAPDPNRatios in your DB? ");
0288   }
0289
0290   // laser timestamps
0291   EcalLaserAPDPNRatios::EcalLaserTimeStamp ltimestamp;
0292   // EcalLaserAPDPNRatios::EcalLaserTimeStampMap::const_iterator ltimeit;
0293   for (int i = 1; i <= 92; i++) {
0294     ltimestamp = lratio.getTimeMap()[i];
0295     //  edm::LogInfo("i = ") << std::setprecision(6) << i
0296     //         << ltimestamp.t1.value() << " " << ltimestamp.t2.value() << " : " ;
0297   }
0298   // edm::LogInfo("Tests finished.");
0299
0300   // channel status
0301   const auto& ch_status = context.getData(channelStatusToken_);
0302
0303   EcalChannelStatusMap::const_iterator chit;
0304   chit = ch_status.getMap().find(ebid.rawId());
0305   if (chit != ch_status.getMap().end()) {
0306     EcalChannelStatusCode ch_code = (*chit);
0307     edm::LogInfo("EcalChannelStatus: ") << std::setprecision(6) << ch_code.getStatusCode();
0308   } else {
0309     edm::LogInfo("No channel status found for this xtal! something wrong with EcalChannelStatus in your DB? ");
0310   }
0311
0312   // laser transparency correction
0313
0314   // Mask to ignore sample
0315   const auto& smask = context.getData(sampleMaskToken_);
0316   edm::LogInfo("Sample Mask EB ") << std::hex << smask.getEcalSampleMaskRecordEB() << " EE " << std::hex
0317                                   << smask.getEcalSampleMaskRecordEE();
0318
0319   /*
0320     std::cout << "make CLHEP matrices from vector<vector<Ecalweight>>" << std::endl;
0321     CLHEP::HepMatrix clmat1(3,8,0);
0322     CLHEP::HepMatrix clmat2(3,8,0);
0323     for(int irow=0; irow<3; irow++) {
0324      for(int icol=0; icol<8; icol++) {
0325        clmat1[irow][icol] = (mat1[irow])[icol]();
0326        clmat2[irow][icol] = (mat2[irow])[icol]();
0327      }
0328    }
0329    std::cout << clmat1 << std::endl;
0330    std::cout << clmat2 << std::endl;
0331 */
0332
0333 }  //end of ::Analyze()
0334 //DEFINE_FWK_MODULE(EcalTrivialObjectAnalyzer);