Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5746) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5747) @@ -252,5 +252,6 @@ if (*p++ >= fSaturationLimit) - sat++; + if (!sat) + sat = ids-3; } @@ -272,5 +273,6 @@ if (*logain++ >= fSaturationLimit) - sat++; + if (!sat) + sat = ids-3; range++; Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc (revision 5746) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc (revision 5747) @@ -231,5 +231,5 @@ if (*p++ >= fSaturationLimit) if (!sat) - sat = p-first+fHiGainFirst-1; + sat = ids-2; count++; @@ -242,5 +242,5 @@ if (*p++ >= fSaturationLimit) if (!sat) - sat = p-first+fHiGainFirst-1; + sat = ids-2; if (IsNoiseCalculation()) @@ -287,5 +287,5 @@ if (*l++ >= fSaturationLimit) if (!sat) - sat = l-logain+fHiGainLast-1; + sat = ids-2; if (sum>max) @@ -308,5 +308,5 @@ if (*l++ >= fSaturationLimit) if (!sat) - sat = l-logain+fHiGainLast-1; + sat = ids-2; if (sum>max) Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 5746) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 5747) @@ -66,5 +66,4 @@ // which simplifies the Numerical Recipes algorithm. // (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.) -// // // The algorithm to search the time proceeds as follows: @@ -102,13 +101,9 @@ // (Int_t)(fAbMaxPos + fFallTime) // (default: fRiseTime: 1.5, fFallTime: 4.5) -// 3) Sum only half the values of the edge slices -// 4) Sum 1.5*fHiGainSecondDeriv of the not-edge slices using the "natural cubic -// spline with second derivatives set to 0. at the edges. -// (Remember that fHiGainSecondDeriv had been divided by 6.) // // The values: fNumHiGainSamples and fNumLoGainSamples are set to: // 1) If Charge Type: kAmplitude was chosen: 1. -// 2) If Charge Type: kIntegral was chosen: TMath::Floor(fRiseTime + fFallTime) -// or: TMath::Floor(fRiseTime + fFallTime + 1.) in the case of the low-gain +// 2) If Charge Type: kIntegral was chosen: fRiseTime + fFallTime +// or: fRiseTime + fFallTime + 1. in the case of the low-gain // // Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) @@ -127,19 +122,15 @@ // for the case, the extraction type kIntegral has been chosen. // -// Call: - SetTimeType(MExtractTimeAndChargeSpline::kMaximum) for extraction +// Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the +// computation of the amplitude at the maximum (default) and extraction // the position of the maximum (default) -// --> needs 22 function evaluations -// - SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum) for extraction -// the position of the half maximum at the rising edge. -// --> needs max. 34 function evaluations -// - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the -// computation of the amplitude at the maximum (default) // --> no further function evaluation needed // - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the // computation of the integral beneith the spline between fRiseTime // from the position of the maximum to fFallTime after the position of -// the maximum. The Low Gain is computed with one more slice at the falling -// edge. -// --> needs one more simple summation loop over 7 slices. +// the maximum. The Low Gain is computed with half a slice more at the rising +// edge and half a slice more at the falling edge. +// The time of the half maximum is returned. +// --> needs one function evaluations but is more precise // ////////////////////////////////////////////////////////////////////////////// @@ -159,5 +150,5 @@ const Byte_t MExtractTimeAndChargeSpline::fgLoGainFirst = 2; const Byte_t MExtractTimeAndChargeSpline::fgLoGainLast = 14; -const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.025; +const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.05; const Float_t MExtractTimeAndChargeSpline::fgRiseTime = 1.5; const Float_t MExtractTimeAndChargeSpline::fgFallTime = 4.5; @@ -177,5 +168,5 @@ // MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title) - : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), fRandomIter(0) + : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), fHalfMaxPos(0.), fRandomIter(0) { @@ -189,5 +180,4 @@ SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); - SetTimeType(); SetChargeType(); @@ -212,21 +202,4 @@ SetChargeType(kAmplitude); -} - - -//------------------------------------------------------------------- -// -// Set the Time Extraction type. Possible are: -// - kMaximum: Search the maximum of the spline and return its position -// - kHalfMaximum: Search the half maximum left from the maximum and return -// its position -// -void MExtractTimeAndChargeSpline::SetTimeType( ExtractionType_t typ ) -{ - - CLRBIT(fFlags,kMaximum); - CLRBIT(fFlags,kHalfMaximum); - SETBIT(fFlags,typ); - } @@ -316,37 +289,39 @@ const Byte_t *end = first + range; Byte_t *p = first; - Int_t count = 0; + + sat = 0; const Float_t pedes = ped.GetPedestal(); const Float_t ABoffs = ped.GetPedestalABoffset(); - Float_t pedmean[2]; - pedmean[0] = pedes + ABoffs; - pedmean[1] = pedes - ABoffs; + const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; fAbMax = 0.; fAbMaxPos = 0.; + fHalfMax = 0.; + fHalfMaxPos = 0.; Int_t maxpos = 0; + Int_t max = 0; // // Check for saturation in all other slices // + Int_t ids = fHiGainFirst; + Float_t *sample = fHiGainSignal.GetArray(); while (p fAbMax + 0.1) /* the 0.1 is necessary for the ultra-high enery events saturating many slices */ - { - fAbMax = signal; - maxpos = count; + *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1]; + + if (*p > max) + { + maxpos = ids-fHiGainFirst-1; + max = *p; } if (*p++ >= fSaturationLimit) - sat++; - - count++; + if (!sat) + sat = ids-3; + } @@ -359,22 +334,22 @@ { - const Int_t ids = fHiGainFirst + range ; - const Float_t signal = (Float_t)*logain - pedmean[(ids+abflag) & 0x1]; - fHiGainSignal[range] = signal; - - if (signal > fAbMax) + *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1]; + + if (*logain > max) { - fAbMax = signal; - maxpos = range; + maxpos = ids-fHiGainFirst-1; + max = *logain; } - if (*logain >= fSaturationLimit) - sat++; + if (*logain++ >= fSaturationLimit) + if (!sat) + sat = ids-3; range++; - logain++; - } - } - + } + } + + fAbMax = fHiGainSignal[maxpos]; + Float_t pp; @@ -434,17 +409,15 @@ // - // Allow one saturated slice + // Allow no saturated slice // and // Don't start if the maxpos is too close to the limits. // - if (sat > 1 || maxpos < 1 || maxpos > range-2) - { - time = IsExtractionType(kMaximum) - ? (Float_t)(fHiGainFirst + maxpos) - : (Float_t)(fHiGainFirst + maxpos - 1); - + if (sat || maxpos < 1 || maxpos > range-2) + { + dtime = 0.5; if (IsExtractionType(kAmplitude)) { - sum = fAbMax; + sum = fAbMax; + time = (Float_t)(fHiGainFirst + maxpos); return; } @@ -455,7 +428,10 @@ CalcIntegralHiGain(sum, 0.001, fRiseTime + fFallTime); + time = (Float_t)(fHiGainFirst + maxpos - 1); return; } + dtime = fResolution; + // // Now find the maximum @@ -533,6 +509,6 @@ // Try a better precision. // - const Float_t up = fAbMaxPos+step - 1.5*fResolution; - const Float_t lo = fAbMaxPos-step + 1.5*fResolution; + const Float_t up = fAbMaxPos+step - 3.0*fResolution; + const Float_t lo = fAbMaxPos-step + 3.0*fResolution; const Float_t maxpossave = fAbMaxPos; @@ -541,5 +517,5 @@ b = x - lower; - step = fResolution; // step size of 83 ps + step = 2.*fResolution; // step size of 0.1 FADC slices while (x= 0) - { - klo--; - if (fHiGainSignal[klo] < fHalfMax) - break; - } - - khi = klo+1; - // - // 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 = 20; - 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) - break; - - step /= 2.; - } - - time = (Float_t)fHiGainFirst + x; - dtime = fResolution; - } - if (IsExtractionType(kAmplitude)) { - sum = fAbMax; + time = fAbMaxPos + (Int_t)fHiGainFirst; + sum = fAbMax; return; } - if (IsExtractionType(kIntegral)) - { - // - // Now integrate the whole thing! - // - - Float_t start = fAbMaxPos - fRiseTime; - Float_t last = fAbMaxPos + fFallTime; - - const Int_t diff = int(last) - range; - - if (diff > 0) - { - last -= diff; - start -= diff; - } - - CalcIntegralHiGain(sum, start, last); - } + 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 >= 0) + { + klo--; + if (fHiGainSignal[klo] < fHalfMax) + break; + } + + khi = klo+1; + // + // 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 = 20; + 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) + break; + + step /= 2.; + } + + time = (Float_t)fHiGainFirst + x; + // + // Now integrate the whole thing! + // + + Float_t start = fAbMaxPos - fRiseTime; + Float_t last = fAbMaxPos + fFallTime; + + const Int_t diff = int(last) - range; + + if (diff > 0) + { + last -= diff; + start -= diff; + } + + CalcIntegralHiGain(sum, start, last); } @@ -718,34 +682,32 @@ const Float_t ABoffs = ped.GetPedestalABoffset(); - Float_t pedmean[2]; - pedmean[0] = pedes + ABoffs; - pedmean[1] = pedes - ABoffs; + const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs }; fAbMax = 0.; fAbMaxPos = 0.; Int_t maxpos = 0; - Int_t count = 0; + Int_t max = 0; // // Check for saturation in all other slices // + Int_t ids = fLoGainFirst; + Float_t *sample = fLoGainSignal.GetArray(); while (p fAbMax + 0.1) - { - fAbMax = signal; - maxpos = count; + *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1]; + + if (*p > max) + { + maxpos = ids-fLoGainFirst-1; + max = *p; } if (*p++ >= fSaturationLimit) sat++; - - count++; - } + } + + fAbMax = fLoGainSignal[maxpos]; Float_t pp; @@ -808,12 +770,10 @@ // Don't start if the maxpos is too close to the limits. // - if (sat || maxpos < 2 || maxpos > range-2) - { - time = IsExtractionType(kMaximum) - ? (Float_t)(fLoGainFirst + maxpos) - : (Float_t)(fLoGainFirst + maxpos - 1); - + if (sat || maxpos < TMath::Ceil(fRiseTime+0.45) || maxpos > range-2) + { + dtime = 0.5; if (IsExtractionType(kAmplitude)) { + time = (Float_t)(fLoGainFirst + maxpos); sum = fAbMax; return; @@ -825,20 +785,9 @@ CalcIntegralLoGain(sum, 0.001, fRiseTime + fFallTime + 1.); + time = (Float_t)(fLoGainFirst + maxpos - 1); return; } - if (maxpos < (Int_t)(fRiseTime+2.)) - { - time = IsExtractionType(kMaximum) - ? (Float_t)(fLoGainFirst + maxpos) - : (Float_t)(fLoGainFirst + maxpos - 1); - - if (maxpos > range-2) - CalcIntegralLoGain(sum, (Float_t)range - fRiseTime - fFallTime-1., (Float_t)range - 0.001); - else - CalcIntegralLoGain(sum, 0.001, fRiseTime + fFallTime + 1.); - - return; - } + dtime = fResolution; // @@ -920,6 +869,6 @@ // Try a better precision. // - const Float_t up = fAbMaxPos+step - 1.5*fResolution; - const Float_t lo = fAbMaxPos-step + 1.5*fResolution; + const Float_t up = fAbMaxPos+step - 3.0*fResolution; + const Float_t lo = fAbMaxPos-step + 3.0*fResolution; const Float_t maxpossave = fAbMaxPos; @@ -928,5 +877,5 @@ b = x - lower; - step = fResolution; // step size of fResolution (33 ps ) + step = 2.*fResolution; // step size of 0.1 FADC slice while (x 0) - { - klo--; - if (fLoGainSignal[klo] < fHalfMax) - break; - } - - khi = klo+1; - // - // 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 = 20; - 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) - break; - - step /= 2.; - } - - time = x + (Int_t)fLoGainFirst; - dtime = fResolution; - } - - if (IsExtractionType(kAmplitude)) - { - sum = fAbMax; + sum = fAbMax; return; } - - if (IsExtractionType(kIntegral)) - { - // - // Now integrate the whole thing! - // - Float_t start = fAbMaxPos - fRiseTime - 0.5; - Float_t last = fAbMaxPos + fFallTime + 0.5; - - const Int_t diff = int(last) - range; - - if (diff > 0) - { - last -= diff; - start -= diff; - } - CalcIntegralLoGain(sum, start, last); - // *fLog << inf << time << " " << sum << " " << start << " " << last << endl; - } - + + 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 > 0) + { + klo--; + if (fLoGainSignal[klo] < fHalfMax) + break; + } + + khi = klo+1; + // + // 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 = 20; + 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) + break; + + step /= 2.; + } + + time = x + (Int_t)fLoGainFirst; + + // + // Now integrate the whole thing! + // + Float_t start = fAbMaxPos - fRiseTime - 0.5; + Float_t last = fAbMaxPos + fFallTime + 0.5; + + const Int_t diff = int(last) - range; + + if (diff > 0) + { + last -= diff; + start -= diff; + } + CalcIntegralLoGain(sum, start, last); } @@ -1248,18 +1184,4 @@ rc = kTRUE; } - if (IsEnvDefined(env, prefix, "Maximum", print)) - { - b = GetEnvValue(env, prefix, "Maximum", IsExtractionType(kMaximum)); - if (b) - SetTimeType(kMaximum); - rc = kTRUE; - } - if (IsEnvDefined(env, prefix, "HalfMaximum", print)) - { - b = GetEnvValue(env, prefix, "HalfMaximum", IsExtractionType(kHalfMaximum)); - if (b) - SetTimeType(kHalfMaximum); - rc = kTRUE; - } return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc; Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h (revision 5746) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h (revision 5747) @@ -72,5 +72,4 @@ void SetFallTime ( const Float_t f=fgFallTime ) { fFallTime = f; } - void SetTimeType ( const ExtractionType_t typ=kHalfMaximum ); void SetChargeType ( const ExtractionType_t typ=kAmplitude);