Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4261) +++ /trunk/MagicSoft/Mars/Changelog (revision 4262) @@ -52,4 +52,6 @@ documentation + * msignal/MExtractTimeFastSpline.cc + - fixed some smaller bugs affecting a small part of the signals Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc (revision 4261) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc (revision 4262) @@ -203,5 +203,5 @@ for (Int_t k=range-2;k>=0;k--) - fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k]; + fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.; // @@ -211,5 +211,254 @@ Float_t lower = (Float_t)maxpos-1.; Float_t upper = (Float_t)maxpos; - Float_t abmax = 0.; + Float_t x = lower; + Float_t y = 0.; + Float_t a = 1.; + Float_t b = 0.; + Int_t klo = maxpos-1; + Int_t khi = maxpos; + Float_t klocont = (Float_t)*(first+klo); + Float_t khicont = (Float_t)*(first+khi); + time = lower; + Float_t abmax = klocont; + + // + // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to + // interval maxpos+1. + // + while (x abmax) + { + abmax = y; + time = x; + } + + } + + if (time < lower+0.1) + { + + upper = (Float_t)maxpos-1.; + lower = (Float_t)maxpos-2.; + x = lower; + a = 1.; + b = 0.; + khi = maxpos-1; + klo = maxpos-2; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+khi); + + while (x abmax) + { + abmax = y; + time = x; + } + } + } + + const Float_t up = time+step-0.055; + const Float_t lo = time-step+0.055; + const Float_t maxpossave = time; + + x = time; + a = upper - x; + b = x - lower; + + step = 0.04; // step size of 83 ps + + while (x abmax) + { + abmax = y; + time = x; + } + + } + + x = maxpossave; + a = upper - x; + b = x - lower; + + while (x>lo) + { + + x -= step; + a += step; + b -= step; + + y = a* klocont + + b* khicont + + (a*a*a-a)*fHiGainSecondDeriv[klo] + + (b*b*b-b)*fHiGainSecondDeriv[khi]; + + if (y > abmax) + { + abmax = y; + time = x; + } + + } + + const Float_t pedes = ped.GetPedestal(); + const Float_t halfmax = pedes + (abmax - pedes)/2.; + + // + // Now, loop from the maximum bin leftward down in order to find the position of the half maximum. + // First, find the right FADC slice: + // + klo = maxpos - 1; + while (klo > maxpos-4) + { + if (*(first+klo) < (Byte_t)halfmax) + break; + klo--; + } + + // + // Loop from the beginning of the slice upwards to reach the halfmax: + // With means of bisection: + // + x = (Float_t)klo; + a = 1.; + b = 0.; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+klo+1); + time = x; + + step = 0.5; + Bool_t back = kFALSE; + + while (step > fResolution) + { + + if (back) + { + x -= step; + a += step; + b -= step; + } + else + { + x += step; + a -= step; + b += step; + } + + y = a*klocont + + b*khicont + + (a*a*a-a)*fHiGainSecondDeriv[klo] + + (b*b*b-b)*fHiGainSecondDeriv[khi]; + + if (y >= halfmax) + back = kTRUE; + else + back = kFALSE; + + step /= 2.; + + } + + time = (Float_t)fHiGainFirst + x; + dtime = fResolution; +} + + +// -------------------------------------------------------------------------- +// +// Calculates the arrival time for each pixel +// +void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped) const +{ + + const Int_t range = fLoGainLast - fLoGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + Byte_t max = 0; + Byte_t maxpos = 0; + + // + // Check for saturation in all other slices + // + while (++p max) + { + max = *p; + maxpos = p-first; + } + + if (*p >= fSaturationLimit) + { + sat++; + break; + } + } + + if (sat) + return; + + if (maxpos < 2) + return; + + Float_t pp; + + p = first; + fLoGainSecondDeriv[0] = 0.; + fLoGainFirstDeriv[0] = 0.; + + for (Int_t i=1;i=0;k--) + fLoGainSecondDeriv[k] = (fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k])/6.; + + // + // Now find the maximum + // + Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one + Float_t lower = (Float_t)maxpos-1.; + Float_t upper = (Float_t)maxpos; Float_t x = lower; Float_t y = 0.; @@ -221,4 +470,5 @@ Float_t khicont = (Float_t)*(first+khi); time = lower; + Float_t abmax = klocont; // @@ -226,5 +476,5 @@ // interval maxpos+1. // - while (x abmax) @@ -246,18 +496,18 @@ } - if (time > upper-0.1) - { - - upper = (Float_t)maxpos+1.; - lower = (Float_t)maxpos; + if (time < lower+0.1) + { + + upper = (Float_t)maxpos-1.; + lower = (Float_t)maxpos-2.; x = lower; a = 1.; b = 0.; - khi = maxpos+1; - klo = maxpos; + khi = maxpos-1; + klo = maxpos-2; klocont = (Float_t)*(first+klo); khicont = (Float_t)*(first+khi); - while (x abmax) @@ -279,5 +529,7 @@ } - const Float_t up = time+step+0.015; + const Float_t up = time+step-0.055; + const Float_t lo = time-step+0.055; + const Float_t maxpossave = time; x = time; @@ -285,5 +537,5 @@ b = x - lower; - step = 0.04; // step size of 83 ps + step = 0.025; // step size of 165 ps while (x abmax) @@ -307,13 +559,7 @@ } - Float_t lo = time-0.225; - - x = time; + x = maxpossave; a = upper - x; b = x - lower; - khi--; - klo--; - klocont = (Float_t)*(first+klo); - khicont = (Float_t)*(first+khi); while (x>lo) @@ -326,6 +572,6 @@ y = a* klocont + b* khicont - + ( (a*a*a-a)*fHiGainSecondDeriv[klo] - +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0; + + (a*a*a-a)*fLoGainSecondDeriv[klo] + + (b*b*b-b)*fLoGainSecondDeriv[khi]; if (y > abmax) @@ -384,6 +630,6 @@ y = a*klocont + b*khicont - + ( (a*a*a-a)*fHiGainSecondDeriv[klo] - +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0; + + (a*a*a-a)*fLoGainSecondDeriv[klo] + + (b*b*b-b)*fLoGainSecondDeriv[khi]; if (y >= halfmax) @@ -393,257 +639,5 @@ step /= 2.; - } - - time = (Float_t)fHiGainFirst + x; - dtime = fResolution; -} - - -// -------------------------------------------------------------------------- -// -// Calculates the arrival time for each pixel -// -void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped) const -{ - - const Int_t range = fLoGainLast - fLoGainFirst + 1; - const Byte_t *end = first + range; - Byte_t *p = first; - Byte_t max = 0; - Byte_t maxpos = 0; - - // - // Check for saturation in all other slices - // - while (++p max) - { - max = *p; - maxpos = p-first; - } - - if (*p >= fSaturationLimit) - { - sat++; - break; - } - } - - if (sat) - return; - - if (maxpos < 2) - return; - - Float_t pp; - - p = first; - fLoGainSecondDeriv[0] = 0.; - fLoGainFirstDeriv[0] = 0.; - - for (Int_t i=1;i=0;k--) - fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k]; - - // - // Now find the maximum - // - Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one - Float_t lower = (Float_t)maxpos-1.; - Float_t upper = (Float_t)maxpos; - Float_t abmax = 0.; - Float_t x = lower; - Float_t y = 0.; - Float_t a = 1.; - Float_t b = 0.; - Int_t klo = maxpos-1; - Int_t khi = maxpos; - Float_t klocont = (Float_t)*(first+klo); - Float_t khicont = (Float_t)*(first+khi); - time = lower; - - // - // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to - // interval maxpos+1. - // - while (x abmax) - { - abmax = y; - time = x; - } - - } - - if (time > upper-0.1) - { - - upper = (Float_t)maxpos+1.; - lower = (Float_t)maxpos; - x = lower; - a = 1.; - b = 0.; - khi = maxpos+1; - klo = maxpos; - klocont = (Float_t)*(first+klo); - khicont = (Float_t)*(first+khi); - - while (x abmax) - { - abmax = y; - time = x; - } - } - } - - const Float_t up = time+step+0.015; - - x = time; - a = upper - x; - b = x - lower; - - step = 0.025; // step size of 165 ps - - while (x abmax) - { - abmax = y; - time = x; - } - - } - - Float_t lo = time-0.225; - - x = time; - a = upper - x; - b = x - lower; - khi--; - klo--; - klocont = (Float_t)*(first+klo); - khicont = (Float_t)*(first+khi); - - while (x>lo) - { - - x -= step; - a += step; - b -= step; - - y = a* klocont - + b* khicont - + ( (a*a*a-a)*fLoGainSecondDeriv[klo] - +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0; - - if (y > abmax) - { - abmax = y; - time = x; - } - - } - - const Float_t pedes = ped.GetPedestal(); - const Float_t halfmax = pedes + (abmax - pedes)/2.; - - // - // Now, loop from the maximum bin leftward down in order to find the position of the half maximum. - // First, find the right FADC slice: - // - klo = maxpos; - while (klo > maxpos-4) - { - if (*(first+klo) < (Byte_t)halfmax) - break; - klo--; - } - - // - // Loop from the beginning of the slice upwards to reach the halfmax: - // With means of bisection: - // - x = (Float_t)klo; - a = 1.; - b = 0.; - klocont = (Float_t)*(first+klo); - khicont = (Float_t)*(first+klo+1); - time = x; - - step = 0.5; - Bool_t back = kFALSE; - - while (step > fResolution) - { - - if (back) - { - x -= step; - a += step; - b -= step; - } - else - { - x += step; - a -= step; - b += step; - } - - y = a*klocont - + b*khicont - + ( (a*a*a-a)*fLoGainSecondDeriv[klo] - +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0; - - if (y >= halfmax) - back = kTRUE; - else - back = kFALSE; - - step /= 2.; + }