Ignore:
Timestamp:
06/02/04 02:34:43 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3996 r4262  
    203203
    204204  for (Int_t k=range-2;k>=0;k--)
    205     fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
     205    fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
    206206 
    207207  //
    … …  
    211211  Float_t lower = (Float_t)maxpos-1.;
    212212  Float_t upper = (Float_t)maxpos;
    213   Float_t abmax = 0.;
     213  Float_t x     = lower;
     214  Float_t y     = 0.;
     215  Float_t a     = 1.;
     216  Float_t b     = 0.;
     217  Int_t   klo = maxpos-1;
     218  Int_t   khi = maxpos;
     219  Float_t klocont = (Float_t)*(first+klo);
     220  Float_t khicont = (Float_t)*(first+khi);
     221  time          = lower;
     222  Float_t abmax = klocont;
     223
     224  //
     225  // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
     226  // interval maxpos+1.
     227  //
     228  while (x<upper-0.1)
     229    {
     230
     231      x += step;
     232      a -= step;
     233      b += step;
     234
     235      y = a*klocont
     236        + b*khicont
     237        + (a*a*a-a)*fHiGainSecondDeriv[klo]
     238        + (b*b*b-b)*fHiGainSecondDeriv[khi];
     239
     240      if (y > abmax)
     241        {
     242          abmax  = y;
     243          time   = x;
     244        }
     245
     246    }
     247
     248 if (time < lower+0.1)
     249    {
     250
     251      upper = (Float_t)maxpos-1.;
     252      lower = (Float_t)maxpos-2.;
     253      x     = lower;
     254      a     = 1.;
     255      b     = 0.;
     256      khi   = maxpos-1;
     257      klo   = maxpos-2;
     258      klocont = (Float_t)*(first+klo);
     259      khicont = (Float_t)*(first+khi);
     260
     261      while (x<upper-0.1)
     262        {
     263
     264          x += step;
     265          a -= step;
     266          b += step;
     267         
     268          y = a* klocont
     269            + b* khicont
     270            + (a*a*a-a)*fHiGainSecondDeriv[klo]
     271            + (b*b*b-b)*fHiGainSecondDeriv[khi];
     272         
     273          if (y > abmax)
     274            {
     275              abmax  = y;
     276              time   = x;
     277            }
     278        }
     279  }
     280 
     281  const Float_t up = time+step-0.055;
     282  const Float_t lo = time-step+0.055;
     283  const Float_t maxpossave = time;
     284
     285  x     = time;
     286  a     = upper - x;
     287  b     = x - lower;
     288
     289  step  = 0.04; // step size of 83 ps
     290
     291  while (x<up)
     292    {
     293
     294      x += step;
     295      a -= step;
     296      b += step;
     297     
     298      y = a* klocont
     299        + b* khicont
     300        + (a*a*a-a)*fHiGainSecondDeriv[klo]
     301        + (b*b*b-b)*fHiGainSecondDeriv[khi];
     302     
     303      if (y > abmax)
     304        {
     305          abmax  = y;
     306          time   = x;
     307        }
     308     
     309    }
     310
     311  x     = maxpossave;
     312  a     = upper - x;
     313  b     = x - lower;
     314
     315  while (x>lo)
     316    {
     317
     318      x -= step;
     319      a += step;
     320      b -= step;
     321     
     322      y = a* klocont
     323        + b* khicont
     324        + (a*a*a-a)*fHiGainSecondDeriv[klo]
     325        + (b*b*b-b)*fHiGainSecondDeriv[khi];
     326     
     327      if (y > abmax)
     328        {
     329          abmax  = y;
     330          time   = x;
     331        }
     332     
     333    }
     334
     335  const Float_t pedes   = ped.GetPedestal();
     336  const Float_t halfmax = pedes + (abmax - pedes)/2.;
     337
     338  //
     339  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
     340  // First, find the right FADC slice:
     341  //
     342  klo   = maxpos - 1;
     343  while (klo > maxpos-4)
     344    {
     345      if (*(first+klo) < (Byte_t)halfmax)
     346        break;
     347      klo--;
     348    }
     349
     350  //
     351  // Loop from the beginning of the slice upwards to reach the halfmax:
     352  // With means of bisection:
     353  //
     354  x     = (Float_t)klo;
     355  a     = 1.;
     356  b     = 0.;
     357  klocont = (Float_t)*(first+klo);
     358  khicont = (Float_t)*(first+klo+1);
     359  time  = x;
     360
     361  step = 0.5;
     362  Bool_t back = kFALSE;
     363
     364  while (step > fResolution)
     365    {
     366     
     367      if (back)
     368        {
     369          x -= step;
     370          a += step;
     371          b -= step;
     372        }
     373      else
     374        {
     375          x += step;
     376          a -= step;
     377          b += step;
     378        }
     379     
     380      y = a*klocont
     381        + b*khicont
     382        + (a*a*a-a)*fHiGainSecondDeriv[klo]
     383        + (b*b*b-b)*fHiGainSecondDeriv[khi];
     384
     385      if (y >= halfmax)
     386        back = kTRUE;
     387      else
     388        back = kFALSE;
     389
     390      step /= 2.;
     391     
     392    }
     393
     394  time  = (Float_t)fHiGainFirst + x;
     395  dtime = fResolution;
     396}
     397
     398
     399// --------------------------------------------------------------------------
     400//
     401// Calculates the arrival time for each pixel
     402//
     403void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime,
     404                                        Byte_t &sat, const MPedestalPix &ped) const
     405{
     406 
     407  const Int_t range = fLoGainLast - fLoGainFirst + 1;
     408  const Byte_t *end = first + range;
     409  Byte_t *p = first;
     410  Byte_t max    = 0;
     411  Byte_t maxpos = 0;
     412
     413  //
     414  // Check for saturation in all other slices
     415  //
     416  while (++p<end)
     417    {
     418      if (*p > max)
     419        {
     420          max    = *p;
     421          maxpos =  p-first;
     422        }
     423
     424      if (*p >= fSaturationLimit)
     425        {
     426          sat++;
     427          break;
     428        }
     429    }
     430 
     431  if (sat)
     432    return;
     433
     434  if (maxpos < 2)
     435    return;
     436 
     437  Float_t pp;
     438
     439  p = first;
     440  fLoGainSecondDeriv[0] = 0.;
     441  fLoGainFirstDeriv[0]  = 0.;
     442
     443  for (Int_t i=1;i<range-1;i++)
     444    {
     445      pp = fLoGainSecondDeriv[i-1] + 4.;
     446      fLoGainSecondDeriv[i] = -1.0/pp;
     447      fLoGainFirstDeriv [i] = *(p+2) - 2.* *(p+1) + *(p);
     448      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
     449      p++;
     450    }
     451
     452  fLoGainSecondDeriv[range-1] = 0.;
     453
     454  for (Int_t k=range-2;k>=0;k--)
     455    fLoGainSecondDeriv[k] = (fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k])/6.;
     456 
     457  //
     458  // Now find the maximum 
     459  //
     460  Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
     461  Float_t lower = (Float_t)maxpos-1.;
     462  Float_t upper = (Float_t)maxpos;
    214463  Float_t x     = lower;
    215464  Float_t y     = 0.;
    … …  
    221470  Float_t khicont = (Float_t)*(first+khi);
    222471  time  = lower;
     472  Float_t abmax = klocont;
    223473
    224474  //
    … …  
    226476  // interval maxpos+1.
    227477  //
    228   while (x<upper+0.1)
     478  while (x<upper-0.1)
    229479    {
    230480
    … …  
    235485      y = a*klocont
    236486        + b*khicont
    237         + (  (a*a*a-a)*fHiGainSecondDeriv[klo]
    238              +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
     487        + (a*a*a-a)*fLoGainSecondDeriv[klo]
     488        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    239489
    240490      if (y > abmax)
    … …  
    246496    }
    247497
    248  if (time > upper-0.1)
    249     {
    250 
    251       upper = (Float_t)maxpos+1.;
    252       lower = (Float_t)maxpos;
     498 if (time < lower+0.1)
     499    {
     500
     501      upper = (Float_t)maxpos-1.;
     502      lower = (Float_t)maxpos-2.;
    253503      x     = lower;
    254504      a     = 1.;
    255505      b     = 0.;
    256       khi   = maxpos+1;
    257       klo   = maxpos;
     506      khi   = maxpos-1;
     507      klo   = maxpos-2;
    258508      klocont = (Float_t)*(first+klo);
    259509      khicont = (Float_t)*(first+khi);
    260510
    261       while (x<upper+0.1)
     511      while (x<upper-0.1)
    262512        {
    263513
    … …  
    268518          y = a* klocont
    269519            + b* khicont
    270             + (  (a*a*a-a)*fHiGainSecondDeriv[klo]
    271                  +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
     520            + (a*a*a-a)*fLoGainSecondDeriv[klo]
     521            + (b*b*b-b)*fLoGainSecondDeriv[khi];
    272522         
    273523          if (y > abmax)
    … …  
    279529  }
    280530 
    281   const Float_t up = time+step+0.015;
     531  const Float_t up = time+step-0.055;
     532  const Float_t lo = time-step+0.055;
     533  const Float_t maxpossave = time;
    282534
    283535  x     = time;
    … …  
    285537  b     = x - lower;
    286538
    287   step  = 0.04; // step size of 83 ps
     539  step  = 0.025; // step size of 165 ps
    288540
    289541  while (x<up)
    … …  
    296548      y = a* klocont
    297549        + b* khicont
    298         + (  (a*a*a-a)*fHiGainSecondDeriv[klo]
    299              +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
     550        + (a*a*a-a)*fLoGainSecondDeriv[klo]
     551        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    300552     
    301553      if (y > abmax)
    … …  
    307559    }
    308560
    309   Float_t lo = time-0.225;
    310 
    311   x     = time;
     561  x     = maxpossave;
    312562  a     = upper - x;
    313563  b     = x - lower;
    314   khi--;
    315   klo--;
    316   klocont = (Float_t)*(first+klo);
    317   khicont = (Float_t)*(first+khi);
    318564
    319565  while (x>lo)
    … …  
    326572      y = a* klocont
    327573        + b* khicont
    328         + (  (a*a*a-a)*fHiGainSecondDeriv[klo]
    329              +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
     574        + (a*a*a-a)*fLoGainSecondDeriv[klo]
     575        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    330576     
    331577      if (y > abmax)
    … …  
    384630      y = a*klocont
    385631        + b*khicont
    386         + (  (a*a*a-a)*fHiGainSecondDeriv[klo]
    387              +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
     632        + (a*a*a-a)*fLoGainSecondDeriv[klo]
     633        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    388634
    389635      if (y >= halfmax)
    … …  
    393639
    394640      step /= 2.;
    395     }
    396 
    397   time  = (Float_t)fHiGainFirst + x;
    398   dtime = fResolution;
    399 }
    400 
    401 
    402 // --------------------------------------------------------------------------
    403 //
    404 // Calculates the arrival time for each pixel
    405 //
    406 void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime,
    407                                         Byte_t &sat, const MPedestalPix &ped) const
    408 {
    409  
    410   const Int_t range = fLoGainLast - fLoGainFirst + 1;
    411   const Byte_t *end = first + range;
    412   Byte_t *p = first;
    413   Byte_t max    = 0;
    414   Byte_t maxpos = 0;
    415 
    416   //
    417   // Check for saturation in all other slices
    418   //
    419   while (++p<end)
    420     {
    421       if (*p > max)
    422         {
    423           max    = *p;
    424           maxpos =  p-first;
    425         }
    426 
    427       if (*p >= fSaturationLimit)
    428         {
    429           sat++;
    430           break;
    431         }
    432     }
    433  
    434   if (sat)
    435     return;
    436 
    437   if (maxpos < 2)
    438     return;
    439  
    440   Float_t pp;
    441 
    442   p = first;
    443   fLoGainSecondDeriv[0] = 0.;
    444   fLoGainFirstDeriv[0]  = 0.;
    445 
    446   for (Int_t i=1;i<range-1;i++)
    447     {
    448       pp = fLoGainSecondDeriv[i-1] + 4.;
    449       fLoGainSecondDeriv[i] = -1.0/pp;
    450       fLoGainFirstDeriv [i] = *(p+2) - 2.* *(p+1) + *(p);
    451       fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
    452       p++;
    453     }
    454 
    455   fLoGainSecondDeriv[range-1] = 0.;
    456 
    457   for (Int_t k=range-2;k>=0;k--)
    458     fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
    459  
    460   //
    461   // Now find the maximum 
    462   //
    463   Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
    464   Float_t lower = (Float_t)maxpos-1.;
    465   Float_t upper = (Float_t)maxpos;
    466   Float_t abmax = 0.;
    467   Float_t x     = lower;
    468   Float_t y     = 0.;
    469   Float_t a     = 1.;
    470   Float_t b     = 0.;
    471   Int_t   klo = maxpos-1;
    472   Int_t   khi = maxpos;
    473   Float_t klocont = (Float_t)*(first+klo);
    474   Float_t khicont = (Float_t)*(first+khi);
    475   time  = lower;
    476 
    477   //
    478   // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
    479   // interval maxpos+1.
    480   //
    481   while (x<upper+0.1)
    482     {
    483 
    484       x += step;
    485       a -= step;
    486       b += step;
    487 
    488       y = a*klocont
    489         + b*khicont
    490         + (  (a*a*a-a)*fLoGainSecondDeriv[klo]
    491              +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
    492 
    493       if (y > abmax)
    494         {
    495           abmax  = y;
    496           time   = x;
    497         }
    498 
    499     }
    500 
    501  if (time > upper-0.1)
    502     {
    503 
    504       upper = (Float_t)maxpos+1.;
    505       lower = (Float_t)maxpos;
    506       x     = lower;
    507       a     = 1.;
    508       b     = 0.;
    509       khi   = maxpos+1;
    510       klo   = maxpos;
    511       klocont = (Float_t)*(first+klo);
    512       khicont = (Float_t)*(first+khi);
    513 
    514       while (x<upper+0.1)
    515         {
    516 
    517           x += step;
    518           a -= step;
    519           b += step;
    520          
    521           y = a* klocont
    522             + b* khicont
    523             + (  (a*a*a-a)*fLoGainSecondDeriv[klo]
    524                  +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
    525          
    526           if (y > abmax)
    527             {
    528               abmax  = y;
    529               time   = x;
    530             }
    531         }
    532   }
    533  
    534   const Float_t up = time+step+0.015;
    535 
    536   x     = time;
    537   a     = upper - x;
    538   b     = x - lower;
    539 
    540   step  = 0.025; // step size of 165 ps
    541 
    542   while (x<up)
    543     {
    544 
    545       x += step;
    546       a -= step;
    547       b += step;
    548      
    549       y = a* klocont
    550         + b* khicont
    551         + (  (a*a*a-a)*fLoGainSecondDeriv[klo]
    552              +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
    553      
    554       if (y > abmax)
    555         {
    556           abmax  = y;
    557           time   = x;
    558         }
    559      
    560     }
    561 
    562   Float_t lo = time-0.225;
    563 
    564   x     = time;
    565   a     = upper - x;
    566   b     = x - lower;
    567   khi--;
    568   klo--;
    569   klocont = (Float_t)*(first+klo);
    570   khicont = (Float_t)*(first+khi);
    571 
    572   while (x>lo)
    573     {
    574 
    575       x -= step;
    576       a += step;
    577       b -= step;
    578      
    579       y = a* klocont
    580         + b* khicont
    581         + (  (a*a*a-a)*fLoGainSecondDeriv[klo]
    582              +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
    583      
    584       if (y > abmax)
    585         {
    586           abmax  = y;
    587           time   = x;
    588         }
    589      
    590     }
    591 
    592   const Float_t pedes   = ped.GetPedestal();
    593   const Float_t halfmax = pedes + (abmax - pedes)/2.;
    594 
    595   //
    596   // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
    597   // First, find the right FADC slice:
    598   //
    599   klo   = maxpos;
    600   while (klo > maxpos-4)
    601     {
    602       if (*(first+klo) < (Byte_t)halfmax)
    603         break;
    604       klo--;
    605     }
    606 
    607   //
    608   // Loop from the beginning of the slice upwards to reach the halfmax:
    609   // With means of bisection:
    610   //
    611   x     = (Float_t)klo;
    612   a     = 1.;
    613   b     = 0.;
    614   klocont = (Float_t)*(first+klo);
    615   khicont = (Float_t)*(first+klo+1);
    616   time  = x;
    617 
    618   step = 0.5;
    619   Bool_t back = kFALSE;
    620 
    621   while (step > fResolution)
    622     {
    623      
    624       if (back)
    625         {
    626           x -= step;
    627           a += step;
    628           b -= step;
    629         }
    630       else
    631         {
    632           x += step;
    633           a -= step;
    634           b += step;
    635         }
    636      
    637       y = a*klocont
    638         + b*khicont
    639         + (  (a*a*a-a)*fLoGainSecondDeriv[klo]
    640              +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
    641 
    642       if (y >= halfmax)
    643         back = kTRUE;
    644       else
    645         back = kFALSE;
    646 
    647       step /= 2.;
     641
    648642    }
    649643
Note: See TracChangeset for help on using the changeset viewer.