Ignore:
Timestamp:
01/29/05 11:35:42 (20 years ago)
Author:
mazin
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/msignal
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc

    r5697 r6107  
    7272//
    7373MExtractFixedWindowSpline::MExtractFixedWindowSpline(const char *name, const char *title)
    74     : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL),
    75       fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL)
    7674{
    7775
    … …  
    8078
    8179  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
    82 }
    83 
    84 MExtractFixedWindowSpline::~MExtractFixedWindowSpline()
    85 {
    86  
    87   if (fHiGainFirstDeriv)
    88     delete [] fHiGainFirstDeriv;
    89   if (fLoGainFirstDeriv)
    90     delete [] fLoGainFirstDeriv;
    91   if (fHiGainSecondDeriv)
    92     delete [] fHiGainSecondDeriv;
    93   if (fLoGainSecondDeriv)
    94     delete [] fLoGainSecondDeriv;
    95  
    9680}
    9781
    … …  
    124108    {
    125109      *fLog << warn << GetDescriptor()
    126             << Form("%s%2i%s%2i",": Hi Gain window size has to be even, set last slice from "
     110            << Form("%s%2i%s%2i",": Hi Gain window size has to be uneven, set last slice from "
    127111                    ,(int)hilast," to ",(int)(hilast-1)) << endl;
    128112      hilast -= 1;
    … …  
    133117      *fLog << warn << GetDescriptor()
    134118            << Form("%s%2i%s%2i",": Hi Gain window is smaller than 2 FADC sampes, set last slice from"
    135                     ,(int)hilast," to ",(int)(hifirst+1)) << endl;
    136       hilast = hifirst+1;
     119                    ,(int)hilast," to ",(int)(hifirst+2)) << endl;
     120      hilast = hifirst+2;
    137121    }
    138122
    … …  
    145129        {
    146130          *fLog << warn << GetDescriptor()
    147                 << Form("%s%2i%s%2i",": Lo Gain window size has to be even, set last slice from "
     131                << Form("%s%2i%s%2i",": Lo Gain window size has to be uneven, set last slice from "
    148132                        ,(int)lolast," to ",(int)(lolast-1)) << endl;
    149133          lolast -= 1;
    … …  
    155139          *fLog << warn << GetDescriptor()
    156140                << Form("%s%2i%s%2i",": Lo Gain window is smaller than 2 FADC sampes, set last slice from"
    157                         ,(int)lolast," to ",(int)(lofirst+1)) << endl;
    158           lolast = lofirst+1;       
     141                        ,(int)lolast," to ",(int)(lofirst+2)) << endl;
     142          lolast = lofirst+2;       
    159143        }
    160144    }
    … …  
    164148
    165149  //
    166   // Very important: Because the spline interpolates between the slice,
     150  // Very important: Because the spline interpolates between the slices,
    167151  //                 the number of samples for the pedestal subtraction
    168152  //                 is now 1 less than with e.g. MExtractFixedWindow
    … …  
    200184                              fLoGainFirst, fLoGainLast, fNumLoGainSamples);
    201185
    202   if (fHiGainFirstDeriv)
    203     delete [] fHiGainFirstDeriv;
    204   if (fLoGainFirstDeriv)
    205     delete [] fLoGainFirstDeriv;
    206   if (fHiGainSecondDeriv)
    207     delete [] fHiGainSecondDeriv;
    208   if (fLoGainSecondDeriv)
    209     delete [] fLoGainSecondDeriv;
    210  
    211186  Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
    212187
    213   fHiGainFirstDeriv  = new Float_t[range];
    214   memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
    215   fHiGainSecondDeriv = new Float_t[range];
    216   memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
     188  fHiGainFirstDeriv.Set(range);
     189  fHiGainSecondDeriv.Set(range);
    217190
    218191  range = fLoGainLast - fLoGainFirst + 1;
    219192
    220   fLoGainFirstDeriv  = new Float_t[range];
    221   memset(fLoGainFirstDeriv,0,range*sizeof(Float_t));
    222   fLoGainSecondDeriv = new Float_t[range];
    223   memset(fLoGainSecondDeriv,0,range*sizeof(Float_t));
     193  fLoGainFirstDeriv.Set(range);
     194  fLoGainSecondDeriv.Set(range);
    224195
    225196  return kTRUE;
    … …  
    243214 
    244215  const Byte_t *end = ptr + fHiGainLast - fHiGainFirst;
    245 
    246216  Int_t range = fHiGainLast - fHiGainFirst + fHiLoLast + 1;
    247217 
    248218  Float_t pp;
    249   Int_t   i = 0;
     219  //  Int_t   i = 0;
    250220
    251221  Int_t summ = 0;
    252   sum = (Float_t)*ptr++/2.;
    253   //
     222  //
     223  // Take half of the first slice content
     224  //
     225  Float_t *firstderiv = fHiGainFirstDeriv.GetArray();
     226  Float_t *secondderiv = fHiGainSecondDeriv.GetArray();
     227  sum = (Float_t)*ptr/2.;
     228  //
     229  // The first slice has already been treated now!
     230  //
     231  ptr++; // i++;
     232  secondderiv++;
     233  firstderiv++;
     234 //
    254235  // Check for saturation in all other slices
    255236  //
    … …  
    258239
    259240      summ += *ptr;
    260       i++;
    261 
    262       pp = fHiGainSecondDeriv[i-1] + 4.;
    263       fHiGainSecondDeriv[i] = -1.0/pp;
    264       fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
    265       fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
     241
     242      // pp = fHiGainSecondDeriv[i-1] + 4.;
     243      // fHiGainSecondDeriv[i] = -1.0/pp;
     244      // fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
     245      // fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
     246
     247      pp = *(secondderiv-1) + 4.;
     248      *secondderiv = -1.0/pp;
     249      *firstderiv  = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
     250      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
    266251
    267252      if (*ptr++ >= fSaturationLimit)
    268253        sat++;
    269     }
    270  
    271   if (*ptr++ >= fSaturationLimit)
    272     sat++;
    273 
    274   if (fHiLoLast == 0)
    275     {
     254
     255      secondderiv++;
     256      firstderiv++;
     257
     258      //      i++;
     259    }
     260 
     261  switch (fHiLoLast)
     262    {
     263    case 0:
     264      // Treat the last slice of the high-gain as half slice:
    276265      sum += (Float_t)*ptr/2.;
    277       fHiGainSecondDeriv[++i] = 0.;     
    278     }
    279   else
    280     {
    281 
    282       //
    283       // There are two overlap slices which we have to treat sepatately
    284       //
    285       // First, the last high-gain slice as center
    286       //
     266      break;
     267    case 1:
     268      // Treat the last slice of the high-gain as full slice:
    287269      summ += *ptr;
    288       i++;
    289 
     270      pp    = *(secondderiv-1) + 4.;
     271      *secondderiv = -1.0/pp;
     272      *firstderiv  = *(logain) - 2.* *(ptr) + *(ptr-1);
     273      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
     274      secondderiv++;
     275      firstderiv++;
    290276      if (*logain >= fSaturationLimit)
    291277        sat++;
    292 
    293       if (fHiLoLast == 1)
    294         sum += (Float_t)*logain/2;
    295       else
    296         {
    297           pp = fHiGainSecondDeriv[i-1] + 4.;
    298           fHiGainSecondDeriv[i] = -1.0/pp;
    299           fHiGainFirstDeriv [i] = *(logain) - 2.* *(ptr) + *(ptr-1);
    300           fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
    301      
    302           //
    303           // Second, the first low-gain slice as center
    304           //
    305           summ += *logain;
    306           i++;
    307          
    308           pp = fHiGainSecondDeriv[i-1] + 4.;
    309           fHiGainSecondDeriv[i] = -1.0/pp;
    310           fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(ptr);
    311           fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
    312          
    313           if (*logain++ >= fSaturationLimit)
    314             sat++;
    315          
    316           end  = logain + fHiLoLast - 2;
    317          
    318           while (logain<end)
    319             {
    320              
    321               summ += *logain;
    322               i++;
    323              
    324               pp = fHiGainSecondDeriv[i-1] + 4.;
    325               fHiGainSecondDeriv[i] = -1.0/pp;
    326               fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(logain-1);
    327               fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
    328              
    329               if (*logain++ >= fSaturationLimit)
    330                 sat++;
    331              
    332             }
    333         }
    334     }
    335  
    336   fHiGainSecondDeriv[range-1] = 0.;
     278      // Treat the first slice of the low-gain as half slice:
     279      sum  += (Float_t)*logain/2;
     280      break;
     281    case 2:
     282      // Treat the last slice of the high-gain as full slice:
     283      summ += *ptr;
     284      pp    = *(secondderiv-1) + 4.;
     285      *secondderiv = -1.0/pp;
     286      *firstderiv  = *(logain) - 2.* *(ptr) + *(ptr-1);
     287      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
     288      secondderiv++;
     289      firstderiv++;
     290      // Treat the last first slice of the low-gain as full slice:
     291      summ += *logain;
     292      pp    = *(secondderiv-1) + 4.;
     293      *secondderiv = -1.0/pp;
     294      *firstderiv  = *(logain+1) - 2.* *(logain) + *(ptr);
     295      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
     296      secondderiv++;
     297      firstderiv++;
     298      if (*logain++ >= fSaturationLimit)
     299        sat++;
     300      // Treat the second slice of the low-gain as half slice:
     301      sum  += (Float_t)*logain/2;
     302      if (*logain >= fSaturationLimit)
     303        sat++;
     304      break;
     305    default:
     306      // Treat the last slice of the high-gain as full slice:
     307      summ += *ptr;
     308      pp    = *(secondderiv-1) + 4.;
     309      *secondderiv = -1.0/pp;
     310      *firstderiv  = *(logain) - 2.* *(ptr) + *(ptr-1);
     311      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
     312      secondderiv++;
     313      firstderiv++;
     314      // Treat the last first slice of the low-gain as full slice:
     315      summ += *logain;
     316      pp    = *(secondderiv-1) + 4.;
     317      *secondderiv = -1.0/pp;
     318      *firstderiv  = *(logain+1) - 2.* *(logain) + *(ptr);
     319      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
     320      secondderiv++;
     321      firstderiv++;
     322      if (*logain++ >= fSaturationLimit)
     323        sat++;
     324      // Treat the rest of the slices:
     325      const Byte_t *end = logain+fHiLoLast;
     326      while (logain<end)
     327        {
     328          summ += *logain;
     329          pp    = *(secondderiv-1) + 4.;
     330          *secondderiv = -1.0/pp;
     331          *firstderiv  = *(logain+1) - 2.* *(logain) + *(logain-1);
     332          *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
     333          //      pp = fHiGainSecondDeriv[i-1] + 4.;
     334          //      fHiGainSecondDeriv[i] = -1.0/pp;
     335          //      fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(logain-1);
     336          //      fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
     337          secondderiv++;
     338          firstderiv++;
     339          if (*logain++ >= fSaturationLimit)
     340            sat++;
     341        }
     342      break;
     343    }
     344 
     345  //
     346  // Go back to last but one element:
     347  //
     348  secondderiv--;
     349  firstderiv--;
    337350
    338351  for (Int_t k=range-2;k>0;k--)
    339352    {
    340       fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
    341       sum += 0.25*fHiGainSecondDeriv[k];
    342     }
    343  
     353      *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv;
     354      sum += 0.25* *secondderiv;
     355      secondderiv--;
     356      //      fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
     357      //      sum += 0.25*fHiGainSecondDeriv[k];
     358    }
     359
    344360  sum += (Float_t)summ;
    345361}
    … …  
    360376 
    361377  Float_t pp;
    362   Int_t   i = 0;
    363 
    364   fLoGainSecondDeriv[0] = 0.;
    365   fLoGainFirstDeriv[0]  = 0.;
     378  //  Int_t   i = 0;
    366379
    367380  Int_t summ = 0;
    368   sum = (Float_t)*ptr++/2.;
     381  //
     382  // Take half of the first slice content
     383  //
     384  Float_t *firstderiv = fLoGainFirstDeriv.GetArray();
     385  Float_t *secondderiv = fLoGainSecondDeriv.GetArray();
     386  sum = (Float_t)*ptr/2.;
     387  //
     388  // The first slice has already been treated now!
     389  //
     390  ptr++; // i++;
     391  secondderiv++;
     392  firstderiv++;
    369393  //
    370394  // Check for saturation in all other slices
    … …  
    374398
    375399      summ += *ptr;
    376       i++;
    377 
    378       pp = fLoGainSecondDeriv[i-1] + 4.;
    379       fLoGainSecondDeriv[i] = -1.0/pp;
    380       fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
    381       fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
     400      //      i++;
     401
     402      // pp = fLoGainSecondDeriv[i-1] + 4.;
     403      // fLoGainSecondDeriv[i] = -1.0/pp;
     404      // fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
     405      // fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
     406
     407      pp = *(secondderiv-1) + 4.;
     408      *secondderiv = -1.0/pp;
     409      *firstderiv  = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
     410      *firstderiv  = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
    382411
    383412      if (*ptr++ >= fSaturationLimit)
    384413        sat++;
    385     }
    386  
    387   if (*ptr++ >= fSaturationLimit)
    388     sat++;
    389 
     414
     415      secondderiv++;
     416      firstderiv++;
     417    }
     418 
    390419  sum += (Float_t)*ptr/2.;
    391   fLoGainSecondDeriv[++i] = 0.;     
    392  
     420 
     421  //
     422  // Go back to last but one element:
     423  //
     424  secondderiv--;
     425  firstderiv--;
     426
    393427  for (Int_t k=range-2;k>0;k--)
    394428    {
    395       fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
    396       sum += 0.25*fLoGainSecondDeriv[k];
     429      *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv;
     430      sum += 0.25* *secondderiv;
     431      secondderiv--;
     432      //      fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
     433      //      sum += 0.25*fLoGainSecondDeriv[k];
    397434    }
    398435 

  • trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h

    r4340 r6107  
    44#ifndef MARS_MExtractor
    55#include "MExtractor.h"
     6#endif
     7
     8#ifndef MARS_MArrayF   
     9#include "MArrayF.h"
    610#endif
    711
    … …  
    1620  static const Byte_t  fgLoGainLast;     // Default for fLoGainLast   (now set to: 14)
    1721
    18   Float_t *fHiGainFirstDeriv;
    19   Float_t *fLoGainFirstDeriv; 
    20   Float_t *fHiGainSecondDeriv;
    21   Float_t *fLoGainSecondDeriv; 
     22  MArrayF fHiGainFirstDeriv;             //! Temporary storage
     23  MArrayF fLoGainFirstDeriv;             //! Temporary storage
     24  MArrayF fHiGainSecondDeriv;            //! Temporary storage
     25  MArrayF fLoGainSecondDeriv;            //! Temporary storage
    2226
    2327  Bool_t ReInit    (MParList *pList);
    … …  
    2933
    3034  MExtractFixedWindowSpline(const char *name=NULL, const char *title=NULL);
    31   ~MExtractFixedWindowSpline();
     35  ~MExtractFixedWindowSpline() {}
    3236 
    3337  void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0);
Note: See TracChangeset for help on using the changeset viewer.