00001 #include "CalibTracker/SiStripAPVAnalysis/interface/TT6PedestalCalculator.h" 00002 00003 #include <cmath> 00004 #include <numeric> 00005 #include <algorithm> 00006 using namespace std; 00007 TT6PedestalCalculator::TT6PedestalCalculator(int evnt_ini, 00008 int evnt_iter, float sig_cut) : 00009 numberOfEvents(0), 00010 alreadyUsedEvent(false) 00011 { 00012 if (0) cout << "Constructing TT6PedestalCalculator " << endl; 00013 eventsRequiredToCalibrate = evnt_ini; 00014 eventsRequiredToUpdate = evnt_iter; 00015 cutToAvoidSignal = sig_cut; 00016 init(); 00017 } 00018 // 00019 // Initialization. 00020 // 00021 void TT6PedestalCalculator::init() { 00022 theRawNoise.clear(); 00023 thePedestal.clear(); 00024 thePedSum.clear(); 00025 thePedSqSum.clear(); 00026 theEventPerStrip.clear(); 00027 theStatus.setCalibrating(); 00028 } 00029 // 00030 // -- Destructor 00031 // 00032 TT6PedestalCalculator::~TT6PedestalCalculator() { 00033 if (0) cout << "Destructing TT6PedestalCalculator " << endl; 00034 } 00035 // 00036 // -- Set Pedestal Update Status 00037 // 00038 void TT6PedestalCalculator::updateStatus(){ 00039 00040 if (theStatus.isCalibrating() && 00041 numberOfEvents >= eventsRequiredToCalibrate) { 00042 theStatus.setUpdating(); 00043 } 00044 } 00045 // 00046 // -- Initialize or Update (when needed) Pedestal Values 00047 // 00048 void TT6PedestalCalculator::updatePedestal(ApvAnalysis::RawSignalType& in) { 00049 if (alreadyUsedEvent == false) { 00050 alreadyUsedEvent = true; 00051 numberOfEvents++; 00052 if (theStatus.isCalibrating()) { 00053 initializePedestal(in); 00054 } else if (theStatus.isUpdating()) { 00055 refinePedestal(in); 00056 } 00057 updateStatus(); 00058 } 00059 } 00060 // 00061 // -- Initialize Pedestal Values using a set of events (eventsRequiredToCalibrate) 00062 // 00063 void TT6PedestalCalculator::initializePedestal(ApvAnalysis::RawSignalType& in) { 00064 if (numberOfEvents == 1) { 00065 00066 thePedSum.clear(); 00067 thePedSqSum.clear(); 00068 theEventPerStrip.clear(); 00069 00070 thePedSum.reserve(128); 00071 thePedSqSum.reserve(128); 00072 theEventPerStrip.reserve(128); 00073 00074 thePedSum.resize(in.data.size(), 0.0); 00075 thePedSqSum.resize(in.data.size(), 0.0); 00076 theEventPerStrip.resize(in.data.size(), 0); 00077 } 00078 if (numberOfEvents <= eventsRequiredToCalibrate) { 00079 edm::DetSet<SiStripRawDigi>::const_iterator i = in.data.begin(); 00080 int ii=0; 00081 for (;i!=in.data.end() ; i++) { 00082 thePedSum[ii] += (*i).adc(); 00083 thePedSqSum[ii] += ((*i).adc())*((*i).adc()); 00084 theEventPerStrip[ii]++; 00085 ii++; 00086 } 00087 } 00088 if (numberOfEvents == eventsRequiredToCalibrate) { 00089 thePedestal.clear(); 00090 theRawNoise.clear(); 00091 edm::DetSet<SiStripRawDigi>::const_iterator i = in.data.begin(); 00092 int ii=0; 00093 for (;i!=in.data.end() ; i++) { 00094 double avVal = (theEventPerStrip[ii]) ? thePedSum[ii]/theEventPerStrip[ii]:0.0; 00095 double sqAvVal = (theEventPerStrip[ii]) ? thePedSqSum[ii]/theEventPerStrip[ii]:0.0; 00096 double corr_fac = (theEventPerStrip[ii] > 1) ? (theEventPerStrip[ii]/(theEventPerStrip[ii]-1)) : 1.0; 00097 double rmsVal = (sqAvVal - avVal*avVal > 0.0) ? sqrt(corr_fac * (sqAvVal - avVal*avVal)) : 0.0; 00098 00099 thePedestal.push_back(static_cast<float>(avVal)); 00100 theRawNoise.push_back(static_cast<float>(rmsVal)); 00101 ii++; 00102 } 00103 } 00104 } 00105 // 00106 // -- Update Pedestal Values when needed. 00107 // 00108 void TT6PedestalCalculator::refinePedestal(ApvAnalysis::RawSignalType& in) { 00109 if (((numberOfEvents-eventsRequiredToCalibrate)%eventsRequiredToUpdate) == 1) { 00110 00111 thePedSum.clear(); 00112 thePedSqSum.clear(); 00113 theEventPerStrip.clear(); 00114 00115 thePedSum.reserve(128); 00116 thePedSqSum.reserve(128); 00117 theEventPerStrip.reserve(128); 00118 00119 thePedSum.resize(in.data.size(), 0.0); 00120 thePedSqSum.resize(in.data.size(), 0.0); 00121 theEventPerStrip.resize(in.data.size(), 0); 00122 } 00123 unsigned int ii=0; 00124 ApvAnalysis::RawSignalType::const_iterator i= in.data.begin(); 00125 for (; i < in.data.end(); i++) { 00126 if (fabs((*i).adc()-thePedestal[ii]) < cutToAvoidSignal*theRawNoise[ii]) { 00127 thePedSum[ii] += (*i).adc(); 00128 thePedSqSum[ii] += ((*i).adc())*((*i).adc()); 00129 theEventPerStrip[ii]++; 00130 } 00131 ii++; 00132 } 00133 if (((numberOfEvents-eventsRequiredToCalibrate) % eventsRequiredToUpdate) == 0) { 00134 00135 for (unsigned int iii = 0; iii < in.data.size(); iii++) { 00136 if (theEventPerStrip[iii] > 10 ) { 00137 double avVal = (theEventPerStrip[iii]) ? thePedSum[iii]/theEventPerStrip[iii]:0.0; 00138 double sqAvVal = (theEventPerStrip[iii]) ? thePedSqSum[iii]/theEventPerStrip[iii]:0.0; 00139 double rmsVal = (sqAvVal - avVal*avVal > 0.0) ? sqrt(sqAvVal - avVal*avVal) : 0.0; 00140 00141 00142 if (avVal != 0 ) { 00143 thePedestal[iii] = static_cast<float>(avVal); 00144 theRawNoise[iii] = static_cast<float>(rmsVal); 00145 } 00146 } 00147 } 00148 thePedSum.clear(); 00149 thePedSqSum.clear(); 00150 theEventPerStrip.clear(); 00151 } 00152 } 00153 // 00154 // Define New Event 00155 // 00156 void TT6PedestalCalculator::newEvent(){ 00157 alreadyUsedEvent = false; 00158 }