SiStripCommissioningSources/src/PedsFullNoiseTask.cc

0001
0002 #include "DQM/SiStripCommissioningSources/interface/PedsFullNoiseTask.h"
0003
0004 #include "DataFormats/SiStripCommon/interface/SiStripConstants.h"
0005 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0006 #include "DQM/SiStripCommon/interface/ExtractTObject.h"
0007 #include "DQM/SiStripCommon/interface/UpdateTProfile.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 // -----------------------------------------------------------------------------
0013 //
0014 PedsFullNoiseTask::PedsFullNoiseTask(DQMStore* dqm, const FedChannelConnection& conn, const edm::ParameterSet& pset)
0015     : CommissioningTask(dqm, conn, "PedsFullNoiseTask"), nstrips_(256) {
0016   LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
0017                                   << " Constructing object...";
0018   edm::ParameterSet params = pset.getParameter<edm::ParameterSet>("PedsFullNoiseParameters");
0019   nskip_ = params.getParameter<int>("NrEvToSkipAtStart");
0020   skipped_ = false;
0021   nevpeds_ = params.getParameter<int>("NrEvForPeds");
0022   pedsdone_ = false;
0023   nadcnoise_ = params.getParameter<int>("NrPosBinsNoiseHist");
0024   fillnoiseprofile_ = params.getParameter<bool>("FillNoiseProfile");
0025   useavgcm_ = params.getParameter<bool>("UseAverageCommonMode");
0026   usefloatpeds_ = params.getParameter<bool>("UseFloatPedestals");
0027 }
0028
0029 // -----------------------------------------------------------------------------
0030 //
0031 PedsFullNoiseTask::~PedsFullNoiseTask() {
0032   LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
0033                                   << " Destructing object...";
0034 }
0035
0036 // -----------------------------------------------------------------------------
0037 //
0038 void PedsFullNoiseTask::book() {
0039   LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]";
0040
0041   // pedestal profile histo
0042   pedhist_.isProfile_ = true;
0043   pedhist_.explicitFill_ = false;
0044   if (!pedhist_.explicitFill_) {
0045     pedhist_.vNumOfEntries_.resize(nstrips_, 0);
0046     pedhist_.vSumOfContents_.resize(nstrips_, 0);
0047     pedhist_.vSumOfSquares_.resize(nstrips_, 0);
0048   }
0049   std::string titleped = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0050                                            sistrip::PEDS_FULL_NOISE,
0051                                            sistrip::FED_KEY,
0052                                            fedKey(),
0053                                            sistrip::LLD_CHAN,
0054                                            connection().lldChannel(),
0055                                            sistrip::extrainfo::pedestals_)
0056                              .title();
0057   pedhist_.histo(dqm()->bookProfile(titleped, titleped, nstrips_, -0.5, nstrips_ * 1. - 0.5, 1025, 0., 1025.));
0058
0059   // Noise profile
0060   noiseprof_.isProfile_ = true;
0061   noiseprof_.explicitFill_ = false;
0062   if (!noiseprof_.explicitFill_) {
0063     noiseprof_.vNumOfEntries_.resize(nstrips_, 0);
0064     noiseprof_.vSumOfContents_.resize(nstrips_, 0);
0065     noiseprof_.vSumOfSquares_.resize(nstrips_, 0);
0066   }
0067   std::string titlenoise = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0068                                              sistrip::PEDS_FULL_NOISE,
0069                                              sistrip::FED_KEY,
0070                                              fedKey(),
0071                                              sistrip::LLD_CHAN,
0072                                              connection().lldChannel(),
0073                                              sistrip::extrainfo::noiseProfile_)
0074                                .title();
0075   noiseprof_.histo(dqm()->bookProfile(titlenoise, titlenoise, nstrips_, -0.5, nstrips_ * 1. - 0.5, 1025, 0., 1025.));
0076
0077   // noise 2D compact histo
0078   noisehist_.explicitFill_ = false;
0079   if (!noisehist_.explicitFill_) {
0080     noisehist_.vNumOfEntries_.resize((nstrips_ + 2) * 2 * (nadcnoise_ + 2), 0);
0081   }
0082   std::string titlenoise2d = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0083                                                sistrip::PEDS_FULL_NOISE,
0084                                                sistrip::FED_KEY,
0085                                                fedKey(),
0086                                                sistrip::LLD_CHAN,
0087                                                connection().lldChannel(),
0088                                                sistrip::extrainfo::noise2D_)
0089                                  .title();
0090   noisehist_.histo(dqm()->book2S(
0091       titlenoise2d, titlenoise2d, 2 * nadcnoise_, -nadcnoise_, nadcnoise_, nstrips_, -0.5, nstrips_ * 1. - 0.5));
0092   hist2d_ = (TH2S*)noisehist_.histo()->getTH2S();
0093 }
0094
0095 // -----------------------------------------------------------------------------
0096 //
0097 void PedsFullNoiseTask::fill(const SiStripEventSummary& summary, const edm::DetSet<SiStripRawDigi>& digis) {
0098   // Check number of digis
0099   uint16_t nbins = digis.data.size();
0100   if (nbins != nstrips_) {
0101     edm::LogWarning(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
0102                                            << " " << nstrips_ << " digis expected, but got " << nbins << ". Skipping.";
0103     return;
0104   }
0105
0106   // get the event number of the first event, not necessarily 1 (parallel processing on FUs)
0107   static int32_t firstev = summary.event();
0108
0109   // skipping events
0110   if (!skipped_) {
0111     if (static_cast<int32_t>(summary.event()) - firstev < nskip_) {
0112       return;
0113     } else {  // when all events are skipped
0114       skipped_ = true;
0115       if (nskip_ > 0)
0116         LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
0117                                         << " Done skipping events. Now starting pedestals.";
0118     }
0119   }
0120
0121   // determine pedestals - decoupled from noise determination
0122   if (!pedsdone_) {
0123     if (static_cast<int32_t>(summary.event()) - firstev < nskip_ + nevpeds_) {
0124       // estimate the pedestals
0125       for (uint16_t istrip = 0; istrip < nstrips_; ++istrip) {
0126         updateHistoSet(pedhist_, istrip, digis.data[istrip].adc());
0127       }
0128       return;
0129     } else {  // when pedestals are done
0130       pedsdone_ = true;
0131       // cache the pedestal values for use in the 2D noise estimation
0132       peds_.clear();
0133       pedsfl_.clear();
0134       for (uint16_t iapv = 0; iapv < 2; ++iapv) {
0135         for (uint16_t ibin = 0; ibin < 128; ++ibin) {
0136           uint16_t istrip = (iapv * 128) + ibin;
0137           if (usefloatpeds_) {
0138             pedsfl_.push_back(1. * pedhist_.vSumOfContents_.at(istrip) / pedhist_.vNumOfEntries_.at(istrip));
0139           } else {
0140             peds_.push_back(
0141                 static_cast<int16_t>(1. * pedhist_.vSumOfContents_.at(istrip) / pedhist_.vNumOfEntries_.at(istrip)));
0142           }
0143         }
0144       }
0145       LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
0146                                       << " Rough pedestals done. Now starting noise measurements.";
0147     }
0148   }
0149
0150   // fill (or not) the old-style niose profile
0151   if (fillnoiseprofile_) {
0152     // Calc common mode for both APVs
0153     std::vector<int32_t> cm;
0154     cm.resize(2, 0);
0155     std::vector<uint16_t> adc;
0156     for (uint16_t iapv = 0; iapv < 2; iapv++) {
0157       adc.clear();
0158       adc.reserve(128);
0159       for (uint16_t ibin = 0; ibin < 128; ibin++) {
0160         if ((iapv * 128) + ibin < nbins) {
0161           adc.push_back(digis.data.at((iapv * 128) + ibin).adc());
0162         }
0163       }
0164       sort(adc.begin(), adc.end());
0165       // take median as common mode
0166       uint16_t index = adc.size() % 2 ? adc.size() / 2 : adc.size() / 2 - 1;
0167       cm[iapv] = static_cast<int16_t>(adc[index]);
0168     }
0169     // 1D noise profile - see also further processing in the update() method
0170     for (uint16_t istrip = 0; istrip < nstrips_; ++istrip) {
0171       // calculate the noise in the old way, by subtracting the common mode, but without pedestal subtraction
0172       int16_t noiseval = static_cast<int16_t>(digis.data.at(istrip).adc()) - cm[istrip / 128];
0173       updateHistoSet(noiseprof_, istrip, noiseval);
0174     }
0175   }
0176
0177   // 2D noise histogram
0178   std::vector<int16_t> noisevals, noisevalssorted;
0179   std::vector<float> noisevalsfl, noisevalssortedfl;
0180   for (uint16_t iapv = 0; iapv < 2; ++iapv) {
0181     float totadc = 0;
0182     noisevals.clear();
0183     noisevalsfl.clear();
0184     noisevalssorted.clear();
0185     noisevalssortedfl.clear();
0186     for (uint16_t ibin = 0; ibin < 128; ++ibin) {
0187       uint16_t istrip = (iapv * 128) + ibin;
0188       // calculate the noise after subtracting the pedestal
0189       if (usefloatpeds_) {  // if float pedestals -> before FED processing
0190         noisevalsfl.push_back(static_cast<float>(digis.data.at(istrip).adc()) - pedsfl_.at(istrip));
0191         // now we still have a possible constant shift of the adc values with respect to 0, so we prepare to calculate the median of this shift
0192         if (useavgcm_) {  // if average CM -> before FED processing
0193           totadc += noisevalsfl[ibin];
0194         } else {  // if median CM -> after FED processing
0195           noisevalssortedfl.push_back(noisevalsfl[ibin]);
0196         }
0197       } else {  // if integer pedestals -> after FED processing
0198         noisevals.push_back(static_cast<int16_t>(digis.data.at(istrip).adc()) - peds_.at(istrip));
0199         // now we still have a possible constant shift of the adc values with respect to 0, so we prepare to calculate the median of this shift
0200         if (useavgcm_) {  // if average CM -> before FED processing
0201           totadc += noisevals[ibin];
0202         } else {  // if median CM -> after FED processing
0203           noisevalssorted.push_back(noisevals[ibin]);
0204         }
0205       }
0206     }
0207     // calculate the common mode shift to apply
0208     float cmshift = 0;
0209     if (useavgcm_) {        // if average CM -> before FED processing
0210       if (usefloatpeds_) {  // if float pedestals -> before FED processing
0211         cmshift = totadc / 128;
0212       } else {  // if integer pedestals -> after FED processing
0213         cmshift = static_cast<int16_t>(totadc / 128);
0214       }
0215     } else {                // if median CM -> after FED processing
0216       if (usefloatpeds_) {  // if float pedestals -> before FED processing
0217         // get the median common mode
0218         sort(noisevalssortedfl.begin(), noisevalssortedfl.end());
0219         uint16_t index = noisevalssortedfl.size() % 2 ? noisevalssortedfl.size() / 2 : noisevalssortedfl.size() / 2 - 1;
0220         cmshift = noisevalssortedfl[index];
0221       } else {  // if integer pedestals -> after FED processing
0222         // get the median common mode
0223         sort(noisevalssorted.begin(), noisevalssorted.end());
0224         uint16_t index = noisevalssorted.size() % 2 ? noisevalssorted.size() / 2 : noisevalssorted.size() / 2 - 1;
0225         cmshift = noisevalssorted[index];
0226       }
0227     }
0228     // now loop again to calculate the CM+pedestal subtracted noise values
0229     for (uint16_t ibin = 0; ibin < 128; ++ibin) {
0230       uint16_t istrip = (iapv * 128) + ibin;
0231       // subtract the remaining common mode after subtraction of the rough pedestals
0232       float noiseval = (usefloatpeds_ ? noisevalsfl[ibin] : noisevals[ibin]) - cmshift;
0233       // retrieve the linear binnr through the histogram
0234       uint32_t binnr = hist2d_->GetBin(static_cast<int>(noiseval + nadcnoise_), istrip + 1);
0235       // store the noise value in the 2D histo
0236       updateHistoSet(noisehist_, binnr);  // no value, so weight 1
0237     }
0238   }
0239 }
0240
0241 // -----------------------------------------------------------------------------
0242 //
0243 void PedsFullNoiseTask::update() {
0244   // pedestals
0245   updateHistoSet(pedhist_);
0246
0247   if (fillnoiseprofile_) {
0248     // noise profile (does not use HistoSet directly, as want to plot noise as "contents", not "error")
0249     TProfile* histo = ExtractTObject<TProfile>().extract(noiseprof_.histo());
0250     for (uint16_t ii = 0; ii < noiseprof_.vNumOfEntries_.size(); ++ii) {
0251       float mean = 0.;
0252       float spread = 0.;
0253       float entries = noiseprof_.vNumOfEntries_[ii];
0254       if (entries > 0.) {
0255         mean = noiseprof_.vSumOfContents_[ii] / entries;
0256         spread = sqrt(noiseprof_.vSumOfSquares_[ii] / entries -
0257                       mean * mean);  // nice way to calculate std dev: Sum (x-<x>)^2 / N
0258       }
0259       float error = spread / sqrt(entries);  // uncertainty on std.dev. when no uncertainty on mean
0260       UpdateTProfile::setBinContent(histo, ii + 1, entries, spread, error);
0261     }
0262   }
0263
0264   // noise 2D histo
0265   updateHistoSet(noisehist_);
0266 }
0267 // -----------------------------------------------------------------------------