Reference Manual

Go to the documentation of this file.
00001 
00002 #include "DQM/SiStripCommissioningSources/interface/PedsFullNoiseTask.h"
00003 
00004 #include "DataFormats/SiStripCommon/interface/SiStripConstants.h"
00005 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
00006 #include "DQM/SiStripCommon/interface/ExtractTObject.h"
00007 #include "DQM/SiStripCommon/src/UpdateTProfile.cc"
00008 #include "DQMServices/Core/interface/DQMStore.h"
00009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00011 
00012 #include "boost/lexical_cast.hpp"
00013 
00014 
00015 // -----------------------------------------------------------------------------
00016 //
00017 PedsFullNoiseTask::PedsFullNoiseTask(DQMStore * dqm, const FedChannelConnection & conn, const edm::ParameterSet & pset) :
00018   CommissioningTask(dqm, conn, "PedsFullNoiseTask"),
00019   nstrips_(256)
00020 {
00021   LogTrace(sistrip::mlDqmSource_)
00022     << "[PedsFullNoiseTask::" << __func__ << "]"
00023     << " Constructing object...";
00024   edm::ParameterSet params = pset.getParameter<edm::ParameterSet>("PedsFullNoiseParameters");
00025   nskip_            = params.getParameter<int>("NrEvToSkipAtStart");
00026   skipped_          = false;
00027   nevpeds_          = params.getParameter<int>("NrEvForPeds");
00028   pedsdone_         = false;
00029   nadcnoise_        = params.getParameter<int>("NrPosBinsNoiseHist");
00030   fillnoiseprofile_ = params.getParameter<bool>("FillNoiseProfile");
00031   useavgcm_         = params.getParameter<bool>("UseAverageCommonMode");
00032   usefloatpeds_     = params.getParameter<bool>("UseFloatPedestals");
00033 }
00034 
00035 // -----------------------------------------------------------------------------
00036 //
00037 PedsFullNoiseTask::~PedsFullNoiseTask()
00038 {
00039   LogTrace(sistrip::mlDqmSource_)
00040     << "[PedsFullNoiseTask::" << __func__ << "]"
00041     << " Destructing object...";
00042 }
00043 
00044 // -----------------------------------------------------------------------------
00045 //
00046 void PedsFullNoiseTask::book()
00047 {
00048 
00049   LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]";
00050 
00051   // pedestal profile histo
00052   pedhist_.isProfile_ = true;
00053   pedhist_.explicitFill_ = false;
00054   if (!pedhist_.explicitFill_) {
00055     pedhist_.vNumOfEntries_.resize(nstrips_,0);
00056     pedhist_.vSumOfContents_.resize(nstrips_,0);
00057     pedhist_.vSumOfSquares_.resize(nstrips_,0);
00058   }
00059   std::string titleped = SiStripHistoTitle( sistrip::EXPERT_HISTO,
00060                                             sistrip::PEDS_FULL_NOISE,
00061                                             sistrip::FED_KEY,
00062                                             fedKey(),
00063                                             sistrip::LLD_CHAN,
00064                                             connection().lldChannel(),
00065                                             sistrip::extrainfo::pedestals_).title();
00066   pedhist_.histo( dqm()->bookProfile( titleped, titleped,
00067                                       nstrips_, -0.5, nstrips_*1.-0.5,
00068                                       1025, 0., 1025. ) );
00069 
00070   // Noise profile
00071   noiseprof_.isProfile_ = true;
00072   noiseprof_.explicitFill_ = false;
00073   if (!noiseprof_.explicitFill_) {
00074     noiseprof_.vNumOfEntries_.resize(nstrips_,0);
00075     noiseprof_.vSumOfContents_.resize(nstrips_,0);
00076     noiseprof_.vSumOfSquares_.resize(nstrips_,0);
00077   }
00078   std::string titlenoise = SiStripHistoTitle( sistrip::EXPERT_HISTO, 
00079                                               sistrip::PEDS_FULL_NOISE, 
00080                                               sistrip::FED_KEY, 
00081                                               fedKey(),
00082                                               sistrip::LLD_CHAN, 
00083                                               connection().lldChannel(),
00084                                               sistrip::extrainfo::noiseProfile_ ).title();
00085   noiseprof_.histo( dqm()->bookProfile( titlenoise, titlenoise,
00086                                         nstrips_, -0.5, nstrips_*1.-0.5,
00087                                         1025, 0., 1025. ) );
00088 
00089   // noise 2D compact histo
00090   noisehist_.explicitFill_ = false;
00091   if (!noisehist_.explicitFill_) {
00092     noisehist_.vNumOfEntries_.resize((nstrips_+2)*2*(nadcnoise_+2), 0);
00093   }
00094   std::string titlenoise2d = SiStripHistoTitle( sistrip::EXPERT_HISTO, 
00095                                                 sistrip::PEDS_FULL_NOISE, 
00096                                                 sistrip::FED_KEY, 
00097                                                 fedKey(),
00098                                                 sistrip::LLD_CHAN, 
00099                                                 connection().lldChannel(),
00100                                                 sistrip::extrainfo::noise2D_ ).title();
00101   noisehist_.histo( dqm()->book2S( titlenoise2d, titlenoise2d,
00102                                    2*nadcnoise_, -nadcnoise_, nadcnoise_,
00103                                    nstrips_, -0.5, nstrips_*1.-0.5 ) );
00104   hist2d_ = (TH2S *) noisehist_.histo()->getTH2S();
00105 
00106 }
00107 
00108 // -----------------------------------------------------------------------------
00109 //
00110 void PedsFullNoiseTask::fill( const SiStripEventSummary & summary,
00111                          const edm::DetSet<SiStripRawDigi> & digis )
00112 {
00113 
00114   // Check number of digis
00115   uint16_t nbins = digis.data.size();
00116   if (nbins != nstrips_) {
00117     edm::LogWarning(sistrip::mlDqmSource_)
00118       << "[PedsFullNoiseTask::" << __func__ << "]"
00119       << " " << nstrips_ << " digis expected, but got " << nbins << ". Skipping.";
00120     return;
00121   }
00122 
00123   // get the event number of the first event, not necessarily 1 (parallel processing on FUs)
00124   static int32_t firstev = summary.event();
00125 
00126   // skipping events
00127   if (!skipped_) {
00128     if (static_cast<int32_t>(summary.event()) - firstev < nskip_) {
00129       return;
00130     } else { // when all events are skipped
00131       skipped_ = true;
00132       if (nskip_ > 0) LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
00133                                                       << " Done skipping events. Now starting pedestals.";
00134     }
00135   }
00136 
00137   // determine pedestals - decoupled from noise determination
00138   if (!pedsdone_) {
00139     if (static_cast<int32_t>(summary.event()) - firstev < nskip_ + nevpeds_) {
00140       // estimate the pedestals
00141       for ( uint16_t istrip = 0; istrip < nstrips_; ++istrip ) {
00142         updateHistoSet( pedhist_, istrip, digis.data[istrip].adc() );
00143       }
00144       return;
00145     } else { // when pedestals are done
00146       pedsdone_ = true;
00147       // cache the pedestal values for use in the 2D noise estimation
00148       peds_.clear(); pedsfl_.clear();
00149       for ( uint16_t iapv = 0; iapv < 2; ++iapv ) {
00150         for ( uint16_t ibin = 0; ibin < 128; ++ibin ) {
00151           uint16_t istrip = (iapv*128)+ibin;
00152           if (usefloatpeds_) {
00153             pedsfl_.push_back(1.*pedhist_.vSumOfContents_.at(istrip)/pedhist_.vNumOfEntries_.at(istrip));
00154           } else {
00155             peds_.push_back(static_cast<int16_t>(1.*pedhist_.vSumOfContents_.at(istrip)/pedhist_.vNumOfEntries_.at(istrip)));
00156           }
00157         }
00158       }
00159       LogTrace(sistrip::mlDqmSource_) << "[PedsFullNoiseTask::" << __func__ << "]"
00160                                       << " Rough pedestals done. Now starting noise measurements.";
00161     }
00162   }
00163 
00164   // fill (or not) the old-style niose profile
00165   if (fillnoiseprofile_) {
00166     // Calc common mode for both APVs
00167     std::vector<int32_t> cm; cm.resize(2, 0);
00168     std::vector<uint16_t> adc;
00169     for ( uint16_t iapv = 0; iapv < 2; iapv++ ) { 
00170       adc.clear(); adc.reserve(128);
00171       for ( uint16_t ibin = 0; ibin < 128; ibin++ ) { 
00172         if ( (iapv*128)+ibin < nbins ) { 
00173           adc.push_back( digis.data.at((iapv*128)+ibin).adc() );
00174         }
00175       }
00176       sort( adc.begin(), adc.end() ); 
00177       // take median as common mode
00178       uint16_t index = adc.size()%2 ? adc.size()/2 : adc.size()/2-1;
00179       cm[iapv] = static_cast<int16_t>(adc[index]);
00180     }
00181     // 1D noise profile - see also further processing in the update() method
00182     for ( uint16_t istrip = 0; istrip < nstrips_; ++istrip ) {
00183       // calculate the noise in the old way, by subtracting the common mode, but without pedestal subtraction
00184       int16_t noiseval = static_cast<int16_t>(digis.data.at(istrip).adc()) - cm[istrip/128];
00185       updateHistoSet( noiseprof_, istrip, noiseval );
00186     }
00187   }
00188 
00189   // 2D noise histogram
00190   std::vector<int16_t> noisevals, noisevalssorted;
00191   std::vector<float> noisevalsfl, noisevalssortedfl;
00192   for ( uint16_t iapv = 0; iapv < 2; ++iapv ) { 
00193     float totadc = 0;
00194     noisevals.clear();       noisevalsfl.clear();
00195     noisevalssorted.clear(); noisevalssortedfl.clear();
00196     for ( uint16_t ibin = 0; ibin < 128; ++ibin ) {
00197       uint16_t istrip = (iapv*128)+ibin;
00198       // calculate the noise after subtracting the pedestal
00199       if (usefloatpeds_) { // if float pedestals -> before FED processing
00200         noisevalsfl.push_back( static_cast<float>(digis.data.at(istrip).adc()) - pedsfl_.at(istrip) );
00201         // 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
00202         if (useavgcm_) { // if average CM -> before FED processing
00203           totadc += noisevalsfl[ibin];
00204         } else { // if median CM -> after FED processing
00205           noisevalssortedfl.push_back( noisevalsfl[ibin] );
00206         }
00207       } else { // if integer pedestals -> after FED processing
00208         noisevals.push_back( static_cast<int16_t>(digis.data.at(istrip).adc()) - peds_.at(istrip) );
00209         // 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
00210         if (useavgcm_) { // if average CM -> before FED processing
00211           totadc += noisevals[ibin];
00212         } else { // if median CM -> after FED processing
00213           noisevalssorted.push_back( noisevals[ibin] );
00214         }
00215       }
00216     }
00217     // calculate the common mode shift to apply
00218     float cmshift = 0;
00219     if (useavgcm_) { // if average CM -> before FED processing
00220       if (usefloatpeds_) { // if float pedestals -> before FED processing
00221         cmshift = totadc/128;
00222       } else { // if integer pedestals -> after FED processing
00223         cmshift = static_cast<int16_t>(totadc/128);
00224       }
00225     } else { // if median CM -> after FED processing
00226       if (usefloatpeds_) { // if float pedestals -> before FED processing
00227         // get the median common mode
00228         sort( noisevalssortedfl.begin(), noisevalssortedfl.end() );
00229         uint16_t index = noisevalssortedfl.size()%2 ? noisevalssortedfl.size()/2 : noisevalssortedfl.size()/2-1;
00230         cmshift = noisevalssortedfl[index];
00231       } else { // if integer pedestals -> after FED processing
00232         // get the median common mode
00233         sort( noisevalssorted.begin(), noisevalssorted.end() );
00234         uint16_t index = noisevalssorted.size()%2 ? noisevalssorted.size()/2 : noisevalssorted.size()/2-1;
00235         cmshift = noisevalssorted[index];
00236       }
00237     }
00238     // now loop again to calculate the CM+pedestal subtracted noise values
00239     for ( uint16_t ibin = 0; ibin < 128; ++ibin ) {
00240       uint16_t istrip = (iapv*128)+ibin;
00241       // subtract the remaining common mode after subtraction of the rough pedestals
00242       float noiseval = (usefloatpeds_ ? noisevalsfl[ibin] : noisevals[ibin]) - cmshift;
00243       // retrieve the linear binnr through the histogram
00244       uint32_t binnr = hist2d_->GetBin(static_cast<int>(noiseval+nadcnoise_), istrip+1);
00245       // store the noise value in the 2D histo
00246       updateHistoSet( noisehist_, binnr ); // no value, so weight 1
00247     }
00248   }
00249 
00250 }
00251 
00252 // -----------------------------------------------------------------------------
00253 //
00254 void PedsFullNoiseTask::update()
00255 {
00256 
00257   // pedestals 
00258   updateHistoSet( pedhist_ );
00259 
00260   if (fillnoiseprofile_) {
00261   // noise profile (does not use HistoSet directly, as want to plot noise as "contents", not "error")
00262   TProfile* histo = ExtractTObject<TProfile>().extract( noiseprof_.histo() );
00263     for ( uint16_t ii = 0; ii < noiseprof_.vNumOfEntries_.size(); ++ii ) {
00264       float mean    = 0.;
00265       float spread  = 0.;
00266       float entries = noiseprof_.vNumOfEntries_[ii];
00267       if ( entries > 0. ) {
00268         mean = noiseprof_.vSumOfContents_[ii] / entries;
00269         spread = sqrt( noiseprof_.vSumOfSquares_[ii] / entries - mean * mean );  // nice way to calculate std dev: Sum (x-<x>)^2 / N
00270       }
00271       float error = spread / sqrt(entries); // uncertainty on std.dev. when no uncertainty on mean
00272       UpdateTProfile::setBinContent( histo, ii+1, entries, spread, error );
00273     }
00274   }
00275 
00276   // noise 2D histo
00277   updateHistoSet( noisehist_ );
00278 
00279 }
00280 // -----------------------------------------------------------------------------
/data/refman/pasoursint/CMSSW_5_3_3/src/DQM/SiStripCommissioningSources/src/PedsFullNoiseTask.cc
