Index: fact/tools/rootmacros/zerosearch.C =================================================================== --- fact/tools/rootmacros/zerosearch.C (revision 15367) +++ fact/tools/rootmacros/zerosearch.C (revision 15368) @@ -19,7 +19,7 @@ int VerbosityLevel ){ - vector * ZeroPositions = new vector; -Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0}; +Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0, 0}; + if (step < 1){ if (VerbosityLevel > 0) @@ -67,5 +67,4 @@ } } - } else // search for zero-x-ings on both esges @@ -81,5 +80,4 @@ } } - return ZeroPositions; @@ -172,4 +170,11 @@ << reg->maxPos << " with Value:" << reg->maxVal << endl; } + + if (reg->maxPos <= 0) + { + reg = regions.erase( reg ) ; + --reg; + continue; + } } return regions.size(); @@ -411,34 +416,79 @@ */ +//size_t calcAmplitudeWeightedTime( +// vector ®ions, +// vector &data, + +// ){ + +//} + +size_t calcCFDPositions( + vector ®ions, + vector &cfd_data + ) +{ + vector::iterator reg; + for (reg = regions.begin() ; reg < regions.end() ; reg++) + { + int neg_zero_crossing = -1; + int min_in_cfd = -1; + int pos_zero_crossing = reg->maxPos; + float minVal = 10.0; + + for (int i = pos_zero_crossing; i > pos_zero_crossing - 40; i--) + { + if (i <= 1 || i >= 1439) break; + + if ( cfd_data[i] < minVal) + { + minVal = cfd_data[i]; + min_in_cfd = i; + } + + if (cfd_data[i] > 0 && cfd_data[i+1] <= 0 && cfd_data[i+2] < 0 ) + { + neg_zero_crossing = i+1; + } + } + + reg->cfdMinPos = min_in_cfd; + reg->cfdNegZerocrossing = neg_zero_crossing; + } + return regions.size(); +} + + size_t findTimeOfHalfMaxLeft( - vector ®ions, - vector &data, - float baseline, - int beginRisingEdge, - int endRisingEdge, - int VerbosityLevel) -{ - float pulse_size = 0; // - float thr = 0; - + vector ®ions, + vector &data, + float baseline, + int beginRisingEdge, + int endRisingEdge, + int VerbosityLevel + ) +{ +// float pulse_size = 0; // + float thr = 0; //int begin = beginRisingEdge; - int end = endRisingEdge; - int counter = 0; + int end = 20; + int counter = 0; // local vars for line fitting - float xmean = 0.0; - float ymean = 0.0; - float cov = 0.0; // empiric covarianve between x and y - float var = 0.0; // empiric variance of x + float xmean = 0.0; + float ymean = 0.0; + float cov = 0.0; // empiric covariance between x and y + float var = 0.0; // empiric variance of x // result of line fitting - float slope = 0.0; - float intercept = 0.0; + float slope = 0.0; + float intercept = 0.0; // result of this function -- the position of the half rising edge - float pos_of_thr_Xing = 0.0; - int result = 0; + float pos_of_thr_Xing = 0.0; + int result = 0; + vector::iterator reg; - for (reg = regions.begin() ; reg < regions.end() ; reg++){ - counter = 0; + for (reg = regions.begin() ; reg < regions.end() ; reg++) + { // check if linear rising edge is completely in pipeline // if not --> delete @@ -450,13 +500,4 @@ continue; } - - // The maximum was probably found using smoothed data, - // so I read the value at the position of the maximum from the data - // again. So I rely on a well defined maxPos. - if (VerbosityLevel > 1) cout << "## baseline: " << baseline << endl; - pulse_size = data[reg->maxPos] - baseline; - if (VerbosityLevel > 1) cout << "## pulse_size: " << pulse_size << endl; - thr = pulse_size / 2. + baseline; - if (VerbosityLevel > 1) cout << "## thr: " << thr << endl; // I think 5 slices left of the max there begins the rising edge @@ -469,34 +510,70 @@ // ALARM check for out of range values... + float average_of_max = -1; + int stepwidth = 2; + + counter = 0; + for ( int i = reg->maxPos -stepwidth; i <= reg->maxPos + stepwidth; i++) + { + average_of_max += data[i]; + counter++; + } + average_of_max /= counter; + + thr = average_of_max / 2.0; + + float max_slope = 0; + int max_slope_pos = -1; + int half_height_pos = -1; + bool passed = false; - for (int slice=reg->maxPos; - slice > 1; --slice) - { - if ( data[slice] < 0.9*data[reg->maxPos] && !passed) - { - beginRisingEdge = reg->maxPos - slice; - passed = true; - } - - if ( data[slice] < 0.1*data[reg->maxPos]) - { - endRisingEdge = reg->maxPos - slice; - passed = false; - reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge; - break; - } - } - - + for (int slice=reg->maxPos; slice > reg->maxPos - end; --slice) + { + if (slice <= 0) break; +// if ( data[slice] < 0.9*data[reg->maxPos] && !passed) +// { +// endRisingEdge = slice; +// beginRisingEdge = endRisingEdge -3; +// passed = true; +// } +// else if ( data[slice] < 0) //0.4*data[reg->maxPos] ) +// { +// beginRisingEdge = slice; +// reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge; +// break; +// } + + //calculate max slope on leading edge + if (data[slice] - data[slice-4] > max_slope && slice-4 > 0) + { + max_slope = data[slice] - data[slice-4]; + max_slope_pos = slice-2; + } + + if ( data[slice] <= thr && !passed){ + half_height_pos = slice; + passed = true; + } + } + reg->maxSlope = max_slope/4; + reg->maxSlopePos = max_slope_pos; + reg->halfRisingEdgePos = half_height_pos; + reg->halfHeight = data[half_height_pos]; + + result = half_height_pos; + + endRisingEdge = result +2; + beginRisingEdge = result -2; + + counter = 0; //calculate mean of x and y coordinate - for (int slice=reg->maxPos - beginRisingEdge; - slice > reg->maxPos - endRisingEdge; --slice) - { + for (int slice= beginRisingEdge; slice < endRisingEdge; ++slice) + { + if (slice <= 0) break; + xmean += slice; ymean += data[slice]; counter++; } -// xmean /= beginRisingEdge - endRisingEdge; -// ymean /= beginRisingEdge - endRisingEdge; xmean /= counter; ymean /= counter; @@ -508,7 +585,9 @@ var = 0.0; - for (int slice=reg->maxPos - beginRisingEdge; - slice > reg->maxPos - endRisingEdge; --slice) - { + //calculate covariance and variance + for (int slice=beginRisingEdge; slice < endRisingEdge; ++slice) + { + if (slice <= 0) break; + cov += (data[slice] - ymean) * (slice - xmean); var += (slice - xmean) * (slice - xmean); @@ -517,10 +596,19 @@ if (VerbosityLevel > 2) cout << "## var: " << var << endl; - slope = cov / var; - intercept = ymean - slope * xmean; + slope = cov / var; + intercept = ymean - slope * xmean; + + // The maximum was probably found using smoothed data, + // so I read the value at the position of the maximum from the data + // again. So I rely on a well defined maxPos. +// if (VerbosityLevel > 4) cout << "## baseline: " << baseline << endl; +// pulse_size = data[reg->maxPos] - baseline; +// if (VerbosityLevel > 4) cout << "## pulse_size: " << pulse_size << endl; +// thr = pulse_size / 2. + baseline; +// if (VerbosityLevel > 4) cout << "## thr: " << thr << endl; // now calculate, where the fittet line crosses the threshold - pos_of_thr_Xing = (thr - intercept) / slope; - result = (int)(pos_of_thr_Xing + 0.5); +// pos_of_thr_Xing = (thr - intercept) / slope; +// result = (int)(pos_of_thr_Xing + 0.5); @@ -537,9 +625,14 @@ } - reg->halfRisingEdgePos = result; - - reg->distanceEdgeToMax = reg->maxPos - result; - reg->interceptRisingEdge = intercept; - reg->slopeOfRisingEdge = slope; + reg->distanceEdgeToMax = reg->maxPos - result; + reg->interceptRisingEdge = reg->maxPos - intercept; + reg->slopeOfRisingEdge = slope; + + if (reg->halfRisingEdgePos <= 0 || reg->halfRisingEdgePos >= 1024) + { + reg = regions.erase( reg ) ; + --reg; + continue; + } } return regions.size(); @@ -570,12 +663,12 @@ size_t GetMaxPositions( vector ®ions, - vector &maxPositions, - int VerbosityLevel) -{ - maxPositions.clear(); + vector &positions, + int VerbosityLevel) +{ + positions.clear(); vector::iterator it = regions.begin(); while( it != regions.end() ) { - maxPositions.push_back(it->maxPos); + positions.push_back(it->maxPos); if (VerbosityLevel > 3){ cout << "getting maxPos@ " << it->maxPos << "\t"; @@ -590,12 +683,12 @@ size_t GetEdgePositions( vector ®ions, - vector &edgePos, - int VerbosityLevel) -{ - maxPositions.clear(); + vector &positions, + int VerbosityLevel) +{ + positions.clear(); vector::iterator it = regions.begin(); while( it != regions.end() ) { - maxPositions.push_back(it->halfRisingEdgePos); + positions.push_back(it->halfRisingEdgePos); if (VerbosityLevel > 3){ cout << "getting maxPos@ " << it->halfRisingEdgePos << "\t";