Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 5212) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 5213) @@ -60,6 +60,7 @@ const Byte_t MExtractTimeAndChargeSpline::fgLoGainFirst = 3; const Byte_t MExtractTimeAndChargeSpline::fgLoGainLast = 14; -const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.001; -const Float_t MExtractTimeAndChargeSpline::fgRatioMax2Fall = 0.05; +const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.001; +const Float_t MExtractTimeAndChargeSpline::fgRiseTime = 2.0; +const Float_t MExtractTimeAndChargeSpline::fgFallTime = 4.0; // -------------------------------------------------------------------------- // @@ -85,7 +86,11 @@ SetResolution(); - SetRatioMax2Fall(); + SetRiseTime(); + SetFallTime(); + + SetTimeType(); + SetChargeType(); + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); - fNumHiGainSamples = 1.; @@ -201,6 +206,4 @@ const Byte_t *end = first + range; Byte_t *p = first; - Byte_t max = 0; - Byte_t maxpos = 0; Int_t count = 0; @@ -212,4 +215,8 @@ PedMean[1] = pedes - ABoffs; + fAbMax = 0.; + fAbMaxPos = 0.; + Byte_t maxpos = 0; + // // Check for saturation in all other slices @@ -219,15 +226,17 @@ const Int_t ids = fHiGainFirst + count ; - fHiGainSignal[count] = (Float_t)*p - PedMean[(ids+abflag) & 0x1]; - - if (*p > max) - { - max = *p; + const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1]; + fHiGainSignal[count] = signal; + + if (signal > fAbMax) + { + fAbMax = signal; maxpos = p-first; } - if (*p++ >= fSaturationLimit) + if (*p >= fSaturationLimit) sat++; - + + p++; count++; } @@ -242,31 +251,34 @@ const Int_t ids = fHiGainFirst + range ; - fHiGainSignal[range] = (Float_t)*logain - PedMean[(ids+abflag) & 0x1]; + const Float_t signal = (Float_t)*logain - PedMean[(ids+abflag) & 0x1]; + fHiGainSignal[range] = signal; range++; - if (*logain > max) - { - max = *logain; + if (signal > fAbMax) + { + fAbMax = signal; maxpos = logain-first; } - if (*logain++ >= fSaturationLimit) + if (*logain >= fSaturationLimit) sat++; - } - } - - // - // allow no saturated slice - // - if (sat) - return; - - // + logain++; + } + } + + // + // Allow no saturated slice + // and // Don't start if the maxpos is too close to the left limit. // - if (maxpos < 1) - return; - + if (sat || maxpos < 2) + { + time = IsTimeType(kMaximum) + ? (Float_t)(fHiGainFirst + maxpos) + : (Float_t)(fHiGainFirst + maxpos - 1); + return; + } + Float_t pp; @@ -292,4 +304,5 @@ Float_t lower = (Float_t)maxpos-1.; Float_t upper = (Float_t)maxpos; + fAbMaxPos = upper; Float_t x = lower; Float_t y = 0.; @@ -298,6 +311,4 @@ Int_t klo = maxpos-1; Int_t khi = maxpos; - fAbMax = fHiGainSignal[khi]; - fAbMaxPos = upper; // @@ -439,105 +450,107 @@ } - fHalfMax = fAbMax/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-3) - { - if (fHiGainSignal[klo] < fHalfMax) - break; - klo--; - } - - // - // Loop from the beginning of the slice upwards to reach the fHalfMax: - // With means of bisection: - // - x = (Float_t)klo; - a = 1.; - b = 0.; - - 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*fHiGainSignal[klo] - + b*fHiGainSignal[khi] - + (a*a*a-a)*fHiGainSecondDeriv[klo] - + (b*b*b-b)*fHiGainSecondDeriv[khi]; - - if (y > fHalfMax) - back = kTRUE; - else - back = kFALSE; - - // *fLog << inf << x << " " << y << " " << fHalfMax << endl; - - step /= 2.; - } - - - time = (Float_t)fHiGainFirst + fAbMaxPos; - // *fLog << inf << time << endl; - // time = (Float_t)fHiGainFirst + x; - dtime = fResolution; - sum = fAbMax; - -#if 0 - // - // Now integrate the whole thing! - // - Int_t startslice = (Int_t)(x - 0.75); - Int_t lastslice = (Int_t)(fAbMaxPos + fRatioMax2Fall*fAbMax); - Int_t coverage = lastslice-startslice; - - // - // make them even - // - if (coverage != (coverage & ~1)) - lastslice++; - - if (startslice < 0) - { + if (IsTimeType(kMaximum)) + { + time = (Float_t)fHiGainFirst + fAbMaxPos; + dtime = 0.02; + } + else + { + fHalfMax = fAbMax/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 >= 0) + { + if (fHiGainSignal[klo] < fHalfMax) + break; + klo--; + } + + // + // Loop from the beginning of the slice upwards to reach the fHalfMax: + // With means of bisection: + // + x = (Float_t)klo; + a = 1.; + b = 0.; + + step = 0.5; + Bool_t back = kFALSE; + + Int_t maxcnt = 1000; + Int_t cnt = 0; + + while (TMath::Abs(y-fHalfMax) > fResolution) + { + + if (back) + { + x -= step; + a += step; + b -= step; + } + else + { + x += step; + a -= step; + b += step; + } + + y = a*fHiGainSignal[klo] + + b*fHiGainSignal[khi] + + (a*a*a-a)*fHiGainSecondDeriv[klo] + + (b*b*b-b)*fHiGainSecondDeriv[khi]; + + if (y > fHalfMax) + back = kTRUE; + else + back = kFALSE; + + if (++cnt > maxcnt) + { + *fLog << inf << x << " " << y << " " << fHalfMax << endl; + break; + } + + step /= 2.; + } + + time = (Float_t)fHiGainFirst + x; + dtime = fResolution; + } + + if (IsChargeType(kAmplitude)) + { + sum = fAbMax; return; - // gLog << warn << GetDescriptor() - // << ": High Gain Bounce against left border, your range is too limited! " << endl; - startslice = 0; - } - - if (lastslice > range-1) - { - return; - // gLog << warn << GetDescriptor() - // << ": High Gain Bounce against right border, your range is too limited! " << endl; - lastslice = range-1; - } - - Int_t i = startslice; - sum = 0.5*fHiGainSignal[i] + 0.75*fHiGainSecondDeriv[i]; - - for (i=startslice+1; i max) - { - max = *p; + while (p fAbMax) + { + fAbMax = signal; maxpos = p-first; } if (*p >= fSaturationLimit) - { - sat++; - break; - } - count++; - } - - if (sat) - return; - - if (maxpos < 2) - return; - + sat++; + + p++; + count++; + } + + // + // Allow no saturated slice + // and + // Don't start if the maxpos is too close to the left limit. + // + if (sat || maxpos < 2) + { + time = IsTimeType(kMaximum) + ? (Float_t)(fLoGainFirst + maxpos) + : (Float_t)(fLoGainFirst + maxpos - 1); + return; + } + Float_t pp; - p = first; fLoGainSecondDeriv[0] = 0.; fLoGainFirstDeriv[0] = 0.; @@ -606,5 +628,4 @@ fLoGainFirstDeriv [i] = fLoGainSignal[i+1] - fLoGainSignal[i] - fLoGainSignal[i] + fLoGainSignal[i-1]; fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp; - p++; } @@ -616,24 +637,20 @@ // 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.; - Float_t a = 1.; - Float_t b = 0.; - Int_t klo = maxpos-1; - Int_t khi = maxpos; - Float_t klocont = fLoGainSignal[klo]; - Float_t khicont = fLoGainSignal[khi]; - fAbMax = klocont; - fAbMaxPos = lower; - - - // - // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to - // interval maxpos+1. - // - while (x fAbMax) { - fAbMax = y; + fAbMax = y; fAbMaxPos = x; } - } - - if (fAbMaxPos < 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 = fLoGainSignal[klo]; - khicont = fLoGainSignal[khi]; - - while (x upper-0.1) + { + + upper = (Float_t)maxpos+1.; + lower = (Float_t)maxpos; + x = lower; + a = 1.; + b = 0.; + khi = maxpos+1; + klo = maxpos; + + while (x fAbMax) { @@ -685,17 +705,24 @@ fAbMaxPos = x; } + // *fLog << inf << x << " " << y << " " << fAbMaxPos << endl; + } } - - const Float_t up = fAbMaxPos+step-0.035; - const Float_t lo = fAbMaxPos-step+0.035; - const Float_t maxpossave = fAbMaxPos; - + + + // + // Now, the time, abmax and khicont and klocont are set correctly within the previous precision. + // Try a better precision. + // + const Float_t up = fAbMaxPos+step-0.055; + const Float_t lo = fAbMaxPos-step+0.055; + const Float_t maxpossave = fAbMaxPos; + x = fAbMaxPos; a = upper - x; b = x - lower; - - step = 0.025; // step size of 165 ps - + + step = 0.02; // step size of 42 ps + while (xlo) { @@ -729,6 +773,6 @@ b -= step; - y = a* klocont - + b* khicont + y = a*fLoGainSignal[klo] + + b*fLoGainSignal[khi] + (a*a*a-a)*fLoGainSecondDeriv[klo] + (b*b*b-b)*fLoGainSecondDeriv[khi]; @@ -739,112 +783,110 @@ fAbMaxPos = x; } - - } - - fHalfMax = fAbMax/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)fHalfMax) - break; - klo--; - } - - // - // Loop from the beginning of the slice upwards to reach the fHalfMax: - // With means of bisection: - // - x = (Float_t)klo; - a = 1.; - b = 0.; - klocont = fLoGainSignal[klo]; - khicont = fLoGainSignal[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]; - - if (y >= fHalfMax) - back = kTRUE; - else - back = kFALSE; - - step /= 2.; - } - - time = (Float_t)fLoGainFirst + x; - dtime = fResolution; - - // - // Now integrate the whole thing! - // - Int_t startslice = (Int_t)(time - 1.5); - Int_t lastslice = (Int_t)(fAbMaxPos + 2.*fRatioMax2Fall*fAbMax-1.); - Int_t coverage = lastslice-startslice; - - // - // make them even - // - if (coverage != (coverage & ~1)) - lastslice++; - - if (startslice < 0) - { - // gLog << warn << GetDescriptor() - // << ": Low Gain Bounce against left border, your range is too limited! " << endl; - sum = 0.; + // *fLog << warn << x << " " << y << " " << fAbMaxPos << endl; + } + + if (IsTimeType(kMaximum)) + { + time = (Float_t)fLoGainFirst + fAbMaxPos; + dtime = 0.02; + } + else + { + fHalfMax = fAbMax/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 >= 0) + { + if (fLoGainSignal[klo] < fHalfMax) + break; + klo--; + } + + // + // Loop from the beginning of the slice upwards to reach the fHalfMax: + // With means of bisection: + // + x = (Float_t)klo; + a = 1.; + b = 0.; + + step = 0.5; + Bool_t back = kFALSE; + + Int_t maxcnt = 1000; + Int_t cnt = 0; + + while (TMath::Abs(y-fHalfMax) > fResolution) + { + + if (back) + { + x -= step; + a += step; + b -= step; + } + else + { + x += step; + a -= step; + b += step; + } + + y = a*fLoGainSignal[klo] + + b*fLoGainSignal[khi] + + (a*a*a-a)*fLoGainSecondDeriv[klo] + + (b*b*b-b)*fLoGainSecondDeriv[khi]; + + if (y > fHalfMax) + back = kTRUE; + else + back = kFALSE; + + if (++cnt > maxcnt) + { + *fLog << inf << x << " " << y << " " << fHalfMax << endl; + break; + } + + step /= 2.; + } + + time = (Float_t)fLoGainFirst + x; + dtime = fResolution; + } + + if (IsChargeType(kAmplitude)) + { + sum = fAbMax; return; - startslice = 0; - } - - if (lastslice > range-1) - { - sum = 0.; - return; - // gLog << warn << GetDescriptor() - // << ": Low Gain Bounce against right border, your range is too limited! " << endl; - lastslice = range-1; - } - - Int_t i = startslice; - sum = 0.5*fLoGainSignal[i] + 0.75*fLoGainSecondDeriv[i]; - - for (i=startslice+1; i