Changeset 5152


Ignore:
Timestamp:
09/30/04 15:08:58 (20 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/msignal
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

    r5146 r5152  
    8080// - SetBinningResolution();
    8181//
     82// Sets all weights to 1.
     83//
    8284MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title)
    8385{
     
    9092  SetBinningResolution();
    9193
    92   for (int i=0; i<60; i++){
    93     fw_amp[i]  = 1;
    94     fw_time[i] = 1;
    95   }
     94  SetWeightsFile("");
    9695}
    9796
     
    161160// This is mandatory for the extraction
    162161//
     162// If filenname is empty, then all weights will be set to 1.
     163//
    163164void MExtractTimeAndChargeDigitalFilter::SetWeightsFile(TString filename)
    164165{
     166
     167  fAmpWeightsHiGain .Set(fBinningResolution*fWindowSizeHiGain);
     168  fAmpWeightsLoGain .Set(fBinningResolution*fWindowSizeLoGain);
     169  fTimeWeightsHiGain.Set(fBinningResolution*fWindowSizeHiGain);
     170  fTimeWeightsLoGain.Set(fBinningResolution*fWindowSizeLoGain);
     171 
     172  if (filename.IsNull())
     173    {
     174      for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++)
     175        {
     176          fAmpWeightsHiGain [i] = 1.;
     177          fTimeWeightsHiGain[i] = 1.;
     178        }
     179      for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++)
     180        {
     181          fAmpWeightsLoGain [i] = 1.;
     182          fTimeWeightsLoGain[i] = 1.;
     183        }
     184      return;
     185    }
    165186
    166187  TFile *f = new TFile(filename);
     
    182203    }
    183204 
    184   for (int i=0; i<fBinningResolution*fWindowSizeHiGain;i++)
    185     {
    186     fw_amp [i] = hw_amp->GetBinContent(i+1);
    187     fw_time[i] = hw_time->GetBinContent(i+1);
    188     //    y[i]=hshape->GetBinContent(i+1);
    189     //    derivative[i]=hderivative->GetBinContent(i+1);
    190  
    191     //    cout << "w_amp[" << i << "] = " << fw_amp[i] << endl;
    192   }
     205  for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++)
     206    {
     207      fAmpWeightsHiGain [i] = hw_amp->GetBinContent(i+1);
     208      fTimeWeightsHiGain[i] = hw_time->GetBinContent(i+1);
     209    }
     210  for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++)
     211    {
     212      fAmpWeightsLoGain [i] = hw_amp->GetBinContent(i+1);
     213      fTimeWeightsLoGain[i] = hw_time->GetBinContent(i+1);
     214    }
     215  //    cout << "w_amp[" << i << "] = " << fw_amp[i] << endl;
    193216  delete f;
     217}
     218
     219// --------------------------------------------------------------------------
     220//
     221// ReInit
     222//
     223// Calls:
     224// - MExtractor::ReInit(pList);
     225// - Creates new arrays according to the extraction range
     226//
     227Bool_t MExtractTimeAndChargeDigitalFilter::ReInit(MParList *pList)
     228{
     229
     230  if (!MExtractTimeAndCharge::ReInit(pList))
     231    return kFALSE;
     232
     233  Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
     234
     235  fHiGainSignal.Set(range);
     236
     237  range = fLoGainLast - fLoGainFirst + 1;
     238
     239  fLoGainSignal.Set(range);
     240
     241  return kTRUE;
    194242}
    195243
     
    199247{
    200248
    201   const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
     249  Int_t range = fHiGainLast - fHiGainFirst + 1;
     250  const Byte_t *end = ptr + range;
     251  Byte_t *p     = ptr;
     252  Byte_t maxpos = 0;
     253  Byte_t max    = 0;
     254  Int_t count   = 0;
     255
     256  //
     257  // Check for saturation in all other slices
     258  //
     259  while (p<end)
     260    {
     261
     262      fHiGainSignal[count++] = (Float_t)*p;
     263
     264      if (*p > max)
     265        {
     266          max    = *p;
     267          maxpos =  p-ptr;
     268        }
     269
     270      if (*p++ >= fSaturationLimit)
     271        sat++;
     272    }
     273 
     274  if (fHiLoLast != 0)
     275    {
     276
     277      end = logain + fHiLoLast;
     278
     279      while (logain<end)
     280        {
     281
     282          fHiGainSignal[count++] = (Float_t)*logain;
     283          range++;
     284
     285          if (*logain > max)
     286            {
     287              max    = *logain;
     288              maxpos =  range;
     289            }
     290         
     291          if (*logain++ >= fSaturationLimit)
     292            sat++;
     293        }
     294    }
     295 
     296  //
     297  // allow one saturated slice
     298  //
     299  if (sat > 0)
     300    return;
     301
     302#if 0
     303  if (maxpos < 2)
     304    {
     305      time = (Float_t)fHiGainFirst;
     306      return;
     307    }
     308#endif
    202309 
    203310  Float_t time_sum  = 0.;
    204311  Float_t fmax      = 0.;
    205312  Float_t ftime_max = 0.;
    206  
    207   Float_t pedes = ped.GetPedestal();
    208 
     313  Int_t   max_p     = 0;
     314 
     315  const Float_t pedes  = ped.GetPedestal();
    209316  const Float_t ABoffs = ped.GetPedestalABoffset();
    210317       
     
    216323  // Calculate the sum of the first fWindowSize slices
    217324  //
    218   Byte_t *p     = ptr;
    219   Byte_t *max_p = ptr;
    220 
    221   while (p < end-fWindowSizeHiGain)
     325  for (Int_t i=0;i<range-fWindowSizeHiGain;i++)
    222326    {
    223327      sum      = 0.;
    224328      time_sum = 0.;
    225 
     329     
    226330      //
    227331      // Slide with a window of size fWindowSizeHiGain over the sample
    228332      // and multiply the entries with the corresponding weights
    229333      //
    230       Byte_t *p1 = p;
    231334      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
    232335      {
    233336        const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
    234         const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-ptr)+abflag) & 0x1];
    235         sum      += fw_amp [idx]*pex;
    236         time_sum += fw_time[idx]*pex;
     337        const Float_t pex = fHiGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1];
     338        sum              += fAmpWeightsHiGain [idx]*pex;
     339        time_sum         += fTimeWeightsHiGain[idx]*pex;
     340      }
     341
     342      if (sum>fmax)
     343      {
     344        fmax      = sum;
     345        ftime_max = time_sum;
     346        max_p     = i;
     347      }
     348    } /*   for (Int_t i=0;i<range-fWindowSizeHiGain;i++) */
     349 
     350  if (fmax!=0)
     351    {
     352      ftime_max        /= fmax;
     353      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
     354      Int_t sample_iter = 0;
     355
     356      while ( t_iter > fBinningResolutionHalf-1 || t_iter < -fBinningResolutionHalf )
     357        {
     358          if (t_iter > fBinningResolutionHalf-1)
     359            {
     360              t_iter -= fBinningResolution;
     361              max_p--;
     362              sample_iter--;
     363            }
     364          if (t_iter < -fBinningResolutionHalf)
     365            {
     366              t_iter += fBinningResolution;
     367              max_p++;
     368              sample_iter++;
     369            }
     370        }
     371     
     372      sum = 0.;
     373      //
     374      // Slide with a window of size fWindowSizeHiGain over the sample
     375      // and multiply the entries with the corresponding weights
     376      //
     377      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
     378      {
     379        const Int_t   idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter;
     380        const Int_t   ids = max_p + sample;
     381        const Float_t pex = ids < 0 ? 0. :
     382          ( ids > range ? 0. : fHiGainSignal[ids]-PedMean[(ids+abflag) & 0x1]);
     383        sum              += fAmpWeightsHiGain [idx]*pex;
     384        time_sum         += fTimeWeightsHiGain[idx]*pex;
     385        //        cout << idx << " " << fw_amp[idx] << "  " << pex << "  " << t_iter << "  " << sum << "  " << ((Int_t(p2-ptr)+abflag) & 0x1) << "  " << PedMean[(Int_t(p2-ptr)+abflag) & 0x1] << endl;
     386      }
     387
     388      if (sum != 0.)
     389        time = max_p + 1. + sample_iter
     390             - ((Float_t)t_iter)/fBinningResolution - 0.4 - time_sum/sum
     391             + (Float_t)fHiGainFirst;
     392      else
     393        time = 0.;
     394    } /* if (max!=0) */
     395    else
     396      time=0.;
     397   
     398  return;
     399}
     400
     401void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
     402                                                                 Float_t &time, Float_t &dtime,
     403                                                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
     404{
     405
     406  const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
     407 
     408  Float_t time_sum  = 0;
     409  Float_t fmax      = 0;
     410  Float_t ftime_max = 0;
     411 
     412  Float_t pedes = ped.GetPedestal();
     413
     414  const Float_t ABoffs = ped.GetPedestalABoffset();
     415       
     416  Float_t PedMean[2];
     417  PedMean[0] = pedes + ABoffs;
     418  PedMean[1] = pedes - ABoffs;
     419
     420  //
     421  // Calculate the sum of the first fWindowSize slices
     422  //
     423  sat           = 0;
     424  Byte_t *p     = ptr;
     425  Byte_t *max_p = ptr;
     426
     427  while (p < end-fWindowSizeLoGain)
     428    {
     429      sum      = 0;
     430      time_sum = 0;
     431
     432      //
     433      // Slide with a window of size fWindowSizeLoGain over the sample
     434      // and multiply the entries with the corresponding weights
     435      //
     436      Byte_t *p1 = p++;
     437      for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
     438      {
     439        const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf;
     440        const Int_t adx = (Int_t(p1-ptr)+abflag) & 0x1;
     441        sum      += fAmpWeightsLoGain [idx]*(*p1-PedMean[adx]);
     442        time_sum += fTimeWeightsLoGain[idx]*(*p1-PedMean[adx]);
    237443        p1++;
    238444      }
     
    244450        max_p     = p;
    245451      }
    246       p++;
    247     } /* while (p<ptr+fWindowSizeHiGain-fHiLoLast) */
    248  
    249 #if 0
    250   //
    251   // If part of the "low-Gain" slices are used,
    252   // repeat steps one and two for the logain part until fHiLoLast
    253   //
    254   Byte_t *l = logain;
    255   while (l<logain+fHiLoLast)
    256     {
    257       sum      = 0;
    258       time_sum = 0;
    259 
    260       //
    261       // Slide with a window of size fWindowSizeHiGain over the sample
    262       // and multiply the entries with the corresponding weights
    263       //
    264       Byte_t *p1 = l++;
    265       for (Int_t sample=0; sample < fHiLoLast; sample++)
    266       {
    267         const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
    268         const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-logain)+abflag) & 0x1];
    269         sum      += fw_amp [idx]*pex;
    270         time_sum += fw_time[idx]*pex;
    271         p1++;
    272       }
    273 
    274       if (sum>fmax)
    275       {
    276         fmax      = sum;
    277         ftime_max = time_sum;
    278         max_p     = p;
    279       }
    280 
    281     }
    282 #endif 
     452    } /* while (p<end-fWindowSizeLoGain) */
    283453 
    284454    if (fmax!=0)
     
    305475        }
    306476     
    307       sum        = 0.;
     477      sum       = 0;
    308478      Byte_t *p2 = p;
    309479     
    310480      //
    311       // Slide with a window of size fWindowSizeHiGain over the sample
     481      // Slide with a window of size fWindowSizeLoGain over the sample
    312482      // and multiply the entries with the corresponding weights
    313483      //
    314       for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
     484      for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
    315485      {
    316         const Int_t   idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter;
    317         const Float_t pex = (Float_t)*p2-PedMean[(Int_t(p2-ptr)+abflag) & 0x1];
    318         sum      += fw_amp [idx]*pex;
    319         //        cout << idx << " " << fw_amp[idx] << "  " << pex << "  " << t_iter << "  " << sum << "  " << ((Int_t(p2-ptr)+abflag) & 0x1) << "  " << PedMean[(Int_t(p2-ptr)+abflag) & 0x1] << endl;
    320         time_sum += fw_time[idx]*pex;
     486        const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf+t_iter;
     487        const Int_t adx = (Int_t(p2-ptr)+abflag) & 0x1;
     488        sum      += fAmpWeightsLoGain [idx]*(*p2-PedMean[adx]);
     489        time_sum += fTimeWeightsLoGain[idx]*(*p2-PedMean[adx]);
    321490        p2++;
    322491      }
     
    324493      if (sum != 0)
    325494        time = (Float_t)(max_p - ptr) + 1. + sample_iter
    326              - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ;
     495             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum
     496             + (Float_t)fLoGainFirst;
    327497      else
    328498        time = 0.;
     
    334504}
    335505
    336 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
    337                                                                  Float_t &time, Float_t &dtime,
    338                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    339 {
    340 
    341   const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
    342  
    343   Float_t time_sum  = 0;
    344   Float_t fmax      = 0;
    345   Float_t ftime_max = 0;
    346  
    347   Float_t pedes = ped.GetPedestal();
    348 
    349   const Float_t ABoffs = ped.GetPedestalABoffset();
    350        
    351   Float_t PedMean[2];
    352   PedMean[0] = pedes + ABoffs;
    353   PedMean[1] = pedes - ABoffs;
    354 
    355   //
    356   // Calculate the sum of the first fWindowSize slices
    357   //
    358   sat           = 0;
    359   Byte_t *p     = ptr;
    360   Byte_t *max_p = ptr;
    361 
    362   while (p < end-fWindowSizeLoGain)
    363     {
    364       sum      = 0;
    365       time_sum = 0;
    366 
    367       //
    368       // Slide with a window of size fWindowSizeLoGain over the sample
    369       // and multiply the entries with the corresponding weights
    370       //
    371       Byte_t *p1 = p++;
    372       for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
    373       {
    374         const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf;
    375         const Int_t adx = (Int_t(p1-ptr)+abflag) & 0x1;
    376         sum      += fw_amp [idx]*(*p1-PedMean[adx]);
    377         time_sum += fw_time[idx]*(*p1-PedMean[adx]);
    378         p1++;
    379       }
    380 
    381       if (sum>fmax)
    382       {
    383         fmax      = sum;
    384         ftime_max = time_sum;
    385         max_p     = p;
    386       }
    387     } /* while (p<end-fWindowSizeLoGain) */
    388  
    389     if (fmax!=0)
    390     {
    391       ftime_max        /= fmax;
    392       Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
    393       Int_t sample_iter = 0;
    394       p                 = max_p;
    395 
    396       while((t_iter)>fBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf)
    397         {
    398           if (t_iter>fBinningResolutionHalf-1)
    399             {
    400               t_iter -= fBinningResolution;
    401               p--;
    402               sample_iter--;
    403             }
    404           if (t_iter < -fBinningResolutionHalf)
    405             {
    406               t_iter += fBinningResolution;
    407               p++;
    408               sample_iter++;
    409             }
    410         }
    411      
    412       sum       = 0;
    413       Byte_t *p2 = p;
    414      
    415       //
    416       // Slide with a window of size fWindowSizeLoGain over the sample
    417       // and multiply the entries with the corresponding weights
    418       //
    419       for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
    420       {
    421         const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf+t_iter;
    422         const Int_t adx = (Int_t(p2-ptr)+abflag) & 0x1;
    423         sum      += fw_amp [idx]*(*p2-PedMean[adx]);
    424         time_sum += fw_time[idx]*(*p2-PedMean[adx]);
    425         p2++;
    426       }
    427 
    428       if (sum != 0)
    429         time = (Float_t)(max_p - ptr) + 1. + sample_iter
    430              - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ;
    431       else
    432         time = 0.;
    433     } /* if (max!=0) */
    434     else
    435       time=0.;
    436    
    437   return;
    438 }
    439 
    440506Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
    441507{
  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h

    r5146 r5152  
    44#ifndef MARS_MExtractTimeAndCharge
    55#include "MExtractTimeAndCharge.h"
     6#endif
     7
     8#ifndef MARS_MArrayF
     9#include "MArrayF.h"
    610#endif
    711
     
    1923  static const Int_t  fgBinningResolution;
    2024 
     25  MArrayF fHiGainSignal;               //! Need fast access to the signals in a float way
     26  MArrayF fLoGainSignal;               //! Store them in separate arrays
     27
    2128  Int_t   fWindowSizeHiGain;           
    2229  Int_t   fWindowSizeLoGain;           
     
    2532  Int_t   fBinningResolutionHalf;
    2633 
    27   Float_t fw_amp [60];
    28   Float_t fw_time[60];
     34  MArrayF fAmpWeightsHiGain;
     35  MArrayF fTimeWeightsHiGain;
     36  MArrayF fAmpWeightsLoGain;
     37  MArrayF fTimeWeightsLoGain;
     38
     39  Bool_t ReInit( MParList *pList );
    2940
    3041  Int_t  ReadEnv(const TEnv &env, TString prefix, Bool_t print);
Note: See TracChangeset for help on using the changeset viewer.