00001 #include "DQM/SiStripCommissioningAnalysis/interface/DaqScopeModeAlgorithm.h" 00002 #include "CondFormats/SiStripObjects/interface/DaqScopeModeAnalysis.h" 00003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h" 00004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h" 00005 #include "FWCore/MessageLogger/interface/MessageLogger.h" 00006 #include "TH1.h" 00007 #include <iostream> 00008 #include <iomanip> 00009 #include <cmath> 00010 00011 using namespace sistrip; 00012 00013 // ---------------------------------------------------------------------------- 00014 // 00015 DaqScopeModeAlgorithm::DaqScopeModeAlgorithm( const edm::ParameterSet & pset, DaqScopeModeAnalysis* const anal ) 00016 : CommissioningAlgorithm(anal), 00017 histo_(0,"") 00018 {;} 00019 00020 // ---------------------------------------------------------------------------- 00021 // 00022 void DaqScopeModeAlgorithm::extract( const std::vector<TH1*>& histos ) { 00023 00024 if ( !anal() ) { 00025 edm::LogWarning(mlCommissioning_) 00026 << "[DaqScopeModeAlgorithm::" << __func__ << "]" 00027 << " NULL pointer to base Analysis object!"; 00028 return; 00029 } 00030 00031 // Check 00032 if ( histos.size() != 1 ) { 00033 anal()->addErrorCode(sistrip::numberOfHistos_); 00034 } 00035 00036 // Extract FED key from histo title 00037 if ( !histos.empty() ) { anal()->fedKey( extractFedKey( histos.front() ) ); } 00038 00039 // Extract histograms 00040 std::vector<TH1*>::const_iterator ihis = histos.begin(); 00041 for ( ; ihis != histos.end(); ihis++ ) { 00042 00043 // Check for NULL pointer 00044 if ( !(*ihis) ) { continue; } 00045 00046 // Check name 00047 SiStripHistoTitle title( (*ihis)->GetName() ); 00048 if ( title.runType() != sistrip::DAQ_SCOPE_MODE ) { 00049 anal()->addErrorCode(sistrip::unexpectedTask_); 00050 continue; 00051 } 00052 00053 // Extract timing histo 00054 histo_.first = *ihis; 00055 histo_.second = (*ihis)->GetName(); 00056 00057 } 00058 00059 } 00060 00061 // ---------------------------------------------------------------------------- 00062 // 00063 void DaqScopeModeAlgorithm::analyse() { 00064 00065 if ( !anal() ) { 00066 edm::LogWarning(mlCommissioning_) 00067 << "[DaqScopeModeAlgorithm::" << __func__ << "]" 00068 << " NULL pointer to base Analysis object!"; 00069 return; 00070 } 00071 00072 CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() ); 00073 DaqScopeModeAnalysis* anal = dynamic_cast<DaqScopeModeAnalysis*>( tmp ); 00074 if ( !anal ) { 00075 edm::LogWarning(mlCommissioning_) 00076 << "[DaqScopeModeAlgorithm::" << __func__ << "]" 00077 << " NULL pointer to derived Analysis object!"; 00078 return; 00079 } 00080 00081 if ( !histo_.first ) { 00082 anal->addErrorCode(sistrip::nullPtr_); 00083 return; 00084 } 00085 00086 // Some initialization 00087 std::vector<float> median; 00088 float max_value = -1. * sistrip::invalid_; 00089 float max_contents = -1. * sistrip::invalid_; 00090 float sum = 0.; 00091 float sum2 = 0.; 00092 00093 // Entries, min, mode 00094 uint16_t nbins = static_cast<uint16_t>( histo_.first->GetNbinsX() ); 00095 for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) { 00096 float value = histo_.first->GetBinLowEdge(ibin+1) + histo_.first->GetBinWidth(ibin+1) / 2.; 00097 float contents = histo_.first->GetBinContent(ibin+1); 00098 //float errors = histo_.first->GetBinError(ibin+1); 00099 if ( contents ) { 00100 if ( contents > max_contents ) { anal->mode_ = contents; max_contents = contents; } 00101 if ( value > max_value ) { max_value = value; } 00102 if ( value < anal->min_ ) { anal->min_ = value; } 00103 sum += value * contents; 00104 sum2 += value * contents* value * contents; 00105 median.insert( median.end(), static_cast<uint32_t>(contents), value ); 00106 } 00107 anal->entries_ += contents; 00108 } 00109 00110 // Max 00111 if ( max_value > -1. * sistrip::maximum_ ) { anal->max_ = max_value; } 00112 00113 // Median 00114 sort( median.begin(), median.end() ); 00115 if ( !median.empty() ) { anal->median_ = median[ median.size()%2 ? median.size()/2 : median.size()/2 ]; } 00116 00117 // Mean, rms 00118 if ( anal->entries_ ) { 00119 sum /= static_cast<float>(anal->entries_); 00120 sum2 /= static_cast<float>(anal->entries_); 00121 anal->mean_ = sum; 00122 if ( sum2 > sum*sum ) { anal->rms_ = sqrt( sum2 - sum*sum ); } 00123 } 00124 00125 }