Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6106) +++ trunk/MagicSoft/Mars/Changelog (revision 6107) @@ -20,4 +20,10 @@ -*-*- END OF LINE -*-*- + 2005/01/29 Markus Gaug + + * msignal/MExtractFixedWindowSpline.[h,cc] + - fixed a bug occurring with overlapping windows between high-gain + and low-gain + 2005/01/28 Abelardo Moralejo Index: trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc (revision 6106) +++ trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc (revision 6107) @@ -72,6 +72,4 @@ // MExtractFixedWindowSpline::MExtractFixedWindowSpline(const char *name, const char *title) - : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL), - fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL) { @@ -80,18 +78,4 @@ SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); -} - -MExtractFixedWindowSpline::~MExtractFixedWindowSpline() -{ - - if (fHiGainFirstDeriv) - delete [] fHiGainFirstDeriv; - if (fLoGainFirstDeriv) - delete [] fLoGainFirstDeriv; - if (fHiGainSecondDeriv) - delete [] fHiGainSecondDeriv; - if (fLoGainSecondDeriv) - delete [] fLoGainSecondDeriv; - } @@ -124,5 +108,5 @@ { *fLog << warn << GetDescriptor() - << Form("%s%2i%s%2i",": Hi Gain window size has to be even, set last slice from " + << Form("%s%2i%s%2i",": Hi Gain window size has to be uneven, set last slice from " ,(int)hilast," to ",(int)(hilast-1)) << endl; hilast -= 1; @@ -133,6 +117,6 @@ *fLog << warn << GetDescriptor() << Form("%s%2i%s%2i",": Hi Gain window is smaller than 2 FADC sampes, set last slice from" - ,(int)hilast," to ",(int)(hifirst+1)) << endl; - hilast = hifirst+1; + ,(int)hilast," to ",(int)(hifirst+2)) << endl; + hilast = hifirst+2; } @@ -145,5 +129,5 @@ { *fLog << warn << GetDescriptor() - << Form("%s%2i%s%2i",": Lo Gain window size has to be even, set last slice from " + << Form("%s%2i%s%2i",": Lo Gain window size has to be uneven, set last slice from " ,(int)lolast," to ",(int)(lolast-1)) << endl; lolast -= 1; @@ -155,6 +139,6 @@ *fLog << warn << GetDescriptor() << Form("%s%2i%s%2i",": Lo Gain window is smaller than 2 FADC sampes, set last slice from" - ,(int)lolast," to ",(int)(lofirst+1)) << endl; - lolast = lofirst+1; + ,(int)lolast," to ",(int)(lofirst+2)) << endl; + lolast = lofirst+2; } } @@ -164,5 +148,5 @@ // - // Very important: Because the spline interpolates between the slice, + // Very important: Because the spline interpolates between the slices, // the number of samples for the pedestal subtraction // is now 1 less than with e.g. MExtractFixedWindow @@ -200,26 +184,13 @@ fLoGainFirst, fLoGainLast, fNumLoGainSamples); - if (fHiGainFirstDeriv) - delete [] fHiGainFirstDeriv; - if (fLoGainFirstDeriv) - delete [] fLoGainFirstDeriv; - if (fHiGainSecondDeriv) - delete [] fHiGainSecondDeriv; - if (fLoGainSecondDeriv) - delete [] fLoGainSecondDeriv; - Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; - fHiGainFirstDeriv = new Float_t[range]; - memset(fHiGainFirstDeriv,0,range*sizeof(Float_t)); - fHiGainSecondDeriv = new Float_t[range]; - memset(fHiGainSecondDeriv,0,range*sizeof(Float_t)); + fHiGainFirstDeriv.Set(range); + fHiGainSecondDeriv.Set(range); range = fLoGainLast - fLoGainFirst + 1; - fLoGainFirstDeriv = new Float_t[range]; - memset(fLoGainFirstDeriv,0,range*sizeof(Float_t)); - fLoGainSecondDeriv = new Float_t[range]; - memset(fLoGainSecondDeriv,0,range*sizeof(Float_t)); + fLoGainFirstDeriv.Set(range); + fLoGainSecondDeriv.Set(range); return kTRUE; @@ -243,13 +214,23 @@ const Byte_t *end = ptr + fHiGainLast - fHiGainFirst; - Int_t range = fHiGainLast - fHiGainFirst + fHiLoLast + 1; Float_t pp; - Int_t i = 0; + // Int_t i = 0; Int_t summ = 0; - sum = (Float_t)*ptr++/2.; - // + // + // Take half of the first slice content + // + Float_t *firstderiv = fHiGainFirstDeriv.GetArray(); + Float_t *secondderiv = fHiGainSecondDeriv.GetArray(); + sum = (Float_t)*ptr/2.; + // + // The first slice has already been treated now! + // + ptr++; // i++; + secondderiv++; + firstderiv++; + // // Check for saturation in all other slices // @@ -258,88 +239,123 @@ summ += *ptr; - i++; - - pp = fHiGainSecondDeriv[i-1] + 4.; - fHiGainSecondDeriv[i] = -1.0/pp; - fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1); - fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; + + // pp = fHiGainSecondDeriv[i-1] + 4.; + // fHiGainSecondDeriv[i] = -1.0/pp; + // fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1); + // fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; + + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(ptr+1) - 2.* *(ptr) + *(ptr-1); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; if (*ptr++ >= fSaturationLimit) sat++; - } - - if (*ptr++ >= fSaturationLimit) - sat++; - - if (fHiLoLast == 0) - { + + secondderiv++; + firstderiv++; + + // i++; + } + + switch (fHiLoLast) + { + case 0: + // Treat the last slice of the high-gain as half slice: sum += (Float_t)*ptr/2.; - fHiGainSecondDeriv[++i] = 0.; - } - else - { - - // - // There are two overlap slices which we have to treat sepatately - // - // First, the last high-gain slice as center - // + break; + case 1: + // Treat the last slice of the high-gain as full slice: summ += *ptr; - i++; - + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(logain) - 2.* *(ptr) + *(ptr-1); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; + secondderiv++; + firstderiv++; if (*logain >= fSaturationLimit) sat++; - - if (fHiLoLast == 1) - sum += (Float_t)*logain/2; - else - { - pp = fHiGainSecondDeriv[i-1] + 4.; - fHiGainSecondDeriv[i] = -1.0/pp; - fHiGainFirstDeriv [i] = *(logain) - 2.* *(ptr) + *(ptr-1); - fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; - - // - // Second, the first low-gain slice as center - // - summ += *logain; - i++; - - pp = fHiGainSecondDeriv[i-1] + 4.; - fHiGainSecondDeriv[i] = -1.0/pp; - fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(ptr); - fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; - - if (*logain++ >= fSaturationLimit) - sat++; - - end = logain + fHiLoLast - 2; - - while (logain= fSaturationLimit) - sat++; - - } - } - } - - fHiGainSecondDeriv[range-1] = 0.; + // Treat the first slice of the low-gain as half slice: + sum += (Float_t)*logain/2; + break; + case 2: + // Treat the last slice of the high-gain as full slice: + summ += *ptr; + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(logain) - 2.* *(ptr) + *(ptr-1); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; + secondderiv++; + firstderiv++; + // Treat the last first slice of the low-gain as full slice: + summ += *logain; + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(logain+1) - 2.* *(logain) + *(ptr); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; + secondderiv++; + firstderiv++; + if (*logain++ >= fSaturationLimit) + sat++; + // Treat the second slice of the low-gain as half slice: + sum += (Float_t)*logain/2; + if (*logain >= fSaturationLimit) + sat++; + break; + default: + // Treat the last slice of the high-gain as full slice: + summ += *ptr; + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(logain) - 2.* *(ptr) + *(ptr-1); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; + secondderiv++; + firstderiv++; + // Treat the last first slice of the low-gain as full slice: + summ += *logain; + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(logain+1) - 2.* *(logain) + *(ptr); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; + secondderiv++; + firstderiv++; + if (*logain++ >= fSaturationLimit) + sat++; + // Treat the rest of the slices: + const Byte_t *end = logain+fHiLoLast; + while (logain= fSaturationLimit) + sat++; + } + break; + } + + // + // Go back to last but one element: + // + secondderiv--; + firstderiv--; for (Int_t k=range-2;k>0;k--) { - fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k]; - sum += 0.25*fHiGainSecondDeriv[k]; - } - + *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv; + sum += 0.25* *secondderiv; + secondderiv--; + // fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k]; + // sum += 0.25*fHiGainSecondDeriv[k]; + } + sum += (Float_t)summ; } @@ -360,11 +376,19 @@ Float_t pp; - Int_t i = 0; - - fLoGainSecondDeriv[0] = 0.; - fLoGainFirstDeriv[0] = 0.; + // Int_t i = 0; Int_t summ = 0; - sum = (Float_t)*ptr++/2.; + // + // Take half of the first slice content + // + Float_t *firstderiv = fLoGainFirstDeriv.GetArray(); + Float_t *secondderiv = fLoGainSecondDeriv.GetArray(); + sum = (Float_t)*ptr/2.; + // + // The first slice has already been treated now! + // + ptr++; // i++; + secondderiv++; + firstderiv++; // // Check for saturation in all other slices @@ -374,25 +398,38 @@ summ += *ptr; - i++; - - pp = fLoGainSecondDeriv[i-1] + 4.; - fLoGainSecondDeriv[i] = -1.0/pp; - fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1); - fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp; + // i++; + + // pp = fLoGainSecondDeriv[i-1] + 4.; + // fLoGainSecondDeriv[i] = -1.0/pp; + // fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1); + // fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp; + + pp = *(secondderiv-1) + 4.; + *secondderiv = -1.0/pp; + *firstderiv = *(ptr+1) - 2.* *(ptr) + *(ptr-1); + *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp; if (*ptr++ >= fSaturationLimit) sat++; - } - - if (*ptr++ >= fSaturationLimit) - sat++; - + + secondderiv++; + firstderiv++; + } + sum += (Float_t)*ptr/2.; - fLoGainSecondDeriv[++i] = 0.; - + + // + // Go back to last but one element: + // + secondderiv--; + firstderiv--; + for (Int_t k=range-2;k>0;k--) { - fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k]; - sum += 0.25*fLoGainSecondDeriv[k]; + *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv; + sum += 0.25* *secondderiv; + secondderiv--; + // fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k]; + // sum += 0.25*fLoGainSecondDeriv[k]; } Index: trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h (revision 6106) +++ trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h (revision 6107) @@ -4,4 +4,8 @@ #ifndef MARS_MExtractor #include "MExtractor.h" +#endif + +#ifndef MARS_MArrayF +#include "MArrayF.h" #endif @@ -16,8 +20,8 @@ static const Byte_t fgLoGainLast; // Default for fLoGainLast (now set to: 14) - Float_t *fHiGainFirstDeriv; - Float_t *fLoGainFirstDeriv; - Float_t *fHiGainSecondDeriv; - Float_t *fLoGainSecondDeriv; + MArrayF fHiGainFirstDeriv; //! Temporary storage + MArrayF fLoGainFirstDeriv; //! Temporary storage + MArrayF fHiGainSecondDeriv; //! Temporary storage + MArrayF fLoGainSecondDeriv; //! Temporary storage Bool_t ReInit (MParList *pList); @@ -29,5 +33,5 @@ MExtractFixedWindowSpline(const char *name=NULL, const char *title=NULL); - ~MExtractFixedWindowSpline(); + ~MExtractFixedWindowSpline() {} void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0);