Changeset 5213 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
10/09/04 12:21:26 (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/MExtractTimeAndChargeSpline.cc

    r5207 r5213  
    6060const Byte_t  MExtractTimeAndChargeSpline::fgLoGainFirst  = 3;
    6161const Byte_t  MExtractTimeAndChargeSpline::fgLoGainLast   = 14;
    62 const Float_t MExtractTimeAndChargeSpline::fgResolution    = 0.001;
    63 const Float_t MExtractTimeAndChargeSpline::fgRatioMax2Fall = 0.05; 
     62const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.001;
     63const Float_t MExtractTimeAndChargeSpline::fgRiseTime     = 2.0;
     64const Float_t MExtractTimeAndChargeSpline::fgFallTime     = 4.0;
    6465// --------------------------------------------------------------------------
    6566//
    … …  
    8586
    8687  SetResolution();
    87   SetRatioMax2Fall();
     88  SetRiseTime();
     89  SetFallTime();
     90 
     91  SetTimeType();
     92  SetChargeType();
     93
    8894  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
    89 
    9095
    9196  fNumHiGainSamples = 1.;
    … …  
    201206  const Byte_t *end = first + range;
    202207  Byte_t       *p  = first;
    203   Byte_t max       = 0;
    204   Byte_t maxpos    = 0;
    205208  Int_t  count     = 0;
    206209
    … …  
    212215  PedMean[1] = pedes - ABoffs;
    213216
     217  fAbMax        = 0.;
     218  fAbMaxPos     = 0.;
     219  Byte_t maxpos = 0;
     220
    214221  //
    215222  // Check for saturation in all other slices
    … …  
    219226
    220227      const Int_t ids      = fHiGainFirst + count ;
    221       fHiGainSignal[count] = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
    222 
    223       if (*p > max)
    224         {
    225           max    = *p;
     228      const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
     229      fHiGainSignal[count] = signal;
     230
     231      if (signal > fAbMax)
     232        {
     233          fAbMax = signal;
    226234          maxpos =  p-first;
    227235        }
    228236
    229       if (*p++ >= fSaturationLimit)
     237      if (*p >= fSaturationLimit)
    230238        sat++;
    231                                                    
     239           
     240      p++;
    232241      count++;                                                   
    233242    }
    … …  
    242251
    243252          const Int_t ids      = fHiGainFirst + range ;
    244           fHiGainSignal[range] = (Float_t)*logain - PedMean[(ids+abflag) & 0x1];
     253          const Float_t signal = (Float_t)*logain - PedMean[(ids+abflag) & 0x1];
     254          fHiGainSignal[range] = signal;
    245255          range++;
    246256         
    247           if (*logain > max)
    248             {
    249               max    = *logain;
     257          if (signal > fAbMax)
     258            { 
     259              fAbMax = signal;
    250260              maxpos =  logain-first;
    251261            }
    252262         
    253           if (*logain++ >= fSaturationLimit)
     263          if (*logain >= fSaturationLimit)
    254264            sat++;
    255265
    256         }
    257     }
    258  
    259   //
    260   // allow no saturated slice
    261   //
    262   if (sat)
    263     return;
    264 
    265   //
     266          logain++;
     267        }
     268    }
     269 
     270  //
     271  // Allow no saturated slice
     272  // and
    266273  // Don't start if the maxpos is too close to the left limit.
    267274  //
    268   if (maxpos < 1)
    269     return;
    270 
     275  if (sat || maxpos < 2)
     276    {
     277      time =  IsTimeType(kMaximum)
     278        ? (Float_t)(fHiGainFirst + maxpos)
     279        : (Float_t)(fHiGainFirst + maxpos - 1);
     280      return;
     281    }
     282     
    271283  Float_t pp;
    272284
    … …  
    292304  Float_t lower   = (Float_t)maxpos-1.;
    293305  Float_t upper   = (Float_t)maxpos;
     306  fAbMaxPos       = upper;
    294307  Float_t x       = lower;
    295308  Float_t y       = 0.;
    … …  
    298311  Int_t   klo     = maxpos-1;
    299312  Int_t   khi     = maxpos;
    300   fAbMax          = fHiGainSignal[khi];
    301   fAbMaxPos       = upper;
    302313
    303314  //
    … …  
    439450    }
    440451
    441   fHalfMax = fAbMax/2.;
    442 
    443   //
    444   // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
    445   // First, find the right FADC slice:
    446   //
    447   klo  = maxpos;
    448   while (klo > maxpos-3)
    449     {
    450       if (fHiGainSignal[klo] < fHalfMax)
    451         break;
    452       klo--;
    453     }
    454 
    455   //
    456   // Loop from the beginning of the slice upwards to reach the fHalfMax:
    457   // With means of bisection:
    458   //
    459   x     = (Float_t)klo;
    460   a     = 1.;
    461   b     = 0.;
    462 
    463   step = 0.5;
    464   Bool_t back = kFALSE;
    465 
    466   while (step > fResolution)
    467     {
    468      
    469       if (back)
    470         {
    471           x -= step;
    472           a += step;
    473           b -= step;
    474         }
    475       else
    476         {
    477           x += step;
    478           a -= step;
    479           b += step;
    480         }
    481      
    482       y = a*fHiGainSignal[klo]
    483         + b*fHiGainSignal[khi]
    484         + (a*a*a-a)*fHiGainSecondDeriv[klo]
    485         + (b*b*b-b)*fHiGainSecondDeriv[khi];
    486 
    487       if (y > fHalfMax)
    488         back = kTRUE;
    489       else
    490         back = kFALSE;
    491 
    492       //      *fLog << inf << x << "  " << y << " " << fHalfMax << endl;
    493 
    494       step /= 2.;
    495     }
    496 
    497 
    498   time  = (Float_t)fHiGainFirst + fAbMaxPos;
    499   //  *fLog << inf << time << endl;
    500   //  time  = (Float_t)fHiGainFirst + x;
    501   dtime = fResolution;
    502   sum = fAbMax;
    503 
    504 #if 0
    505   //
    506   // Now integrate the whole thing!
    507   //
    508   Int_t startslice = (Int_t)(x - 0.75);
    509   Int_t lastslice  = (Int_t)(fAbMaxPos + fRatioMax2Fall*fAbMax);
    510   Int_t coverage   = lastslice-startslice;
    511 
    512   //
    513   // make them even
    514   //
    515   if (coverage != (coverage & ~1))
    516     lastslice++;
    517 
    518   if (startslice < 0)
    519     {
     452  if (IsTimeType(kMaximum))
     453    {
     454      time = (Float_t)fHiGainFirst + fAbMaxPos;
     455      dtime = 0.02;
     456    }
     457  else
     458    {
     459      fHalfMax = fAbMax/2.;
     460     
     461      //
     462      // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
     463      // First, find the right FADC slice:
     464      //
     465      klo  = maxpos - 1;
     466      while (klo >= 0)
     467        {
     468          if (fHiGainSignal[klo] < fHalfMax)
     469            break;
     470          klo--;
     471        }
     472     
     473      //
     474      // Loop from the beginning of the slice upwards to reach the fHalfMax:
     475      // With means of bisection:
     476      //
     477      x     = (Float_t)klo;
     478      a     = 1.;
     479      b     = 0.;
     480     
     481      step = 0.5;
     482      Bool_t back = kFALSE;
     483     
     484      Int_t maxcnt = 1000;
     485      Int_t cnt    = 0;
     486
     487      while (TMath::Abs(y-fHalfMax) > fResolution)
     488        {
     489         
     490          if (back)
     491            {
     492              x -= step;
     493              a += step;
     494              b -= step;
     495            }
     496          else
     497            {
     498              x += step;
     499              a -= step;
     500              b += step;
     501            }
     502         
     503          y = a*fHiGainSignal[klo]
     504            + b*fHiGainSignal[khi]
     505            + (a*a*a-a)*fHiGainSecondDeriv[klo]
     506            + (b*b*b-b)*fHiGainSecondDeriv[khi];
     507         
     508          if (y > fHalfMax)
     509            back = kTRUE;
     510          else
     511            back = kFALSE;
     512
     513          if (++cnt > maxcnt)
     514            {
     515              *fLog << inf << x << "  " << y << " " << fHalfMax << endl;
     516              break;
     517            }
     518         
     519          step /= 2.;
     520        }
     521
     522      time  = (Float_t)fHiGainFirst + x;
     523      dtime = fResolution;
     524    }
     525 
     526  if (IsChargeType(kAmplitude))
     527    {
     528      sum = fAbMax;
    520529      return;
    521       //      gLog << warn << GetDescriptor()
    522       //           << ": High Gain Bounce against left border, your range is too limited! " << endl;
    523       startslice = 0;
    524     }
    525 
    526   if (lastslice > range-1)
    527     {
    528       return;
    529       //      gLog << warn << GetDescriptor()
    530       //           << ": High Gain Bounce against right border, your range is too limited! " << endl;
    531       lastslice = range-1;
    532     }
    533 
    534   Int_t i = startslice;
    535   sum = 0.5*fHiGainSignal[i] + 0.75*fHiGainSecondDeriv[i];
    536 
    537   for (i=startslice+1; i<lastslice; i++)
    538     sum += fHiGainSignal[i] + 1.5*fHiGainSecondDeriv[i];
    539  
    540   sum += 0.5*fHiGainSignal[lastslice] + 0.75*fHiGainSecondDeriv[lastslice];
    541 #endif
     530    }
     531
     532  if (IsChargeType(kIntegral))
     533    {
     534      //
     535      // Now integrate the whole thing!
     536      //
     537      Int_t startslice = (Int_t)(fAbMaxPos - fRiseTime);
     538      Int_t lastslice  = (Int_t)(fAbMaxPos + fFallTime);
     539     
     540      if (startslice < 0)
     541        {
     542          lastslice -= startslice;
     543          startslice = 0;
     544        }
     545
     546      Int_t i = startslice;
     547      sum = 0.5*fHiGainSignal[i];
     548     
     549      for (i=startslice+1; i<lastslice; i++)
     550        sum += fHiGainSignal[i] + 1.5*fHiGainSecondDeriv[i];
     551     
     552      sum += 0.5*fHiGainSignal[lastslice];
     553    }
     554 
    542555}
    543556
    … …  
    552565{
    553566 
    554   const Int_t range = fLoGainLast - fLoGainFirst + 1;
     567  Int_t range = fLoGainLast - fLoGainFirst + 1;
    555568  const Byte_t *end = first + range;
    556   Byte_t *p = first;
    557   Byte_t max    = 0;
    558   Byte_t maxpos = 0;
    559   Int_t count   = 0;
     569  Byte_t       *p  = first;
     570  Int_t  count     = 0;
    560571
    561572  Float_t pedes        = ped.GetPedestal();
    562573  const Float_t ABoffs = ped.GetPedestalABoffset();
     574
    563575  Float_t PedMean[2];
    564576  PedMean[0] = pedes + ABoffs;
    565577  PedMean[1] = pedes - ABoffs;
    566578
     579  fAbMax        = 0.;
     580  fAbMaxPos     = 0.;
     581  Byte_t maxpos = 0;
     582
    567583  //
    568584  // Check for saturation in all other slices
    569585  //
    570   while (++p<end)
    571     {
    572 
    573       const Int_t ids       = fLoGainFirst + count ;
    574       fLoGainSignal[count] = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
    575 
    576       if (*p > max)
    577         {
    578           max    = *p;
     586  while (p<end)
     587    {
     588
     589      const Int_t ids      = fLoGainFirst + count ;
     590      const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
     591      fLoGainSignal[count] = signal;
     592
     593      if (signal > fAbMax)
     594        {
     595          fAbMax = signal;
    579596          maxpos =  p-first;
    580597        }
    581598
    582599      if (*p >= fSaturationLimit)
    583         {
    584           sat++;
    585           break;
    586         }
    587       count++;
    588     }
    589  
    590   if (sat)
    591     return;
    592 
    593   if (maxpos < 2)
    594     return;
    595  
     600        sat++;
     601           
     602      p++;
     603      count++;                                                   
     604    }
     605 
     606  //
     607  // Allow no saturated slice
     608  // and
     609  // Don't start if the maxpos is too close to the left limit.
     610  //
     611  if (sat || maxpos < 2)
     612    {
     613      time =  IsTimeType(kMaximum)
     614        ? (Float_t)(fLoGainFirst + maxpos)
     615        : (Float_t)(fLoGainFirst + maxpos - 1);
     616      return;
     617    }
     618     
    596619  Float_t pp;
    597620
    598   p = first;
    599621  fLoGainSecondDeriv[0] = 0.;
    600622  fLoGainFirstDeriv[0]  = 0.;
    … …  
    606628      fLoGainFirstDeriv [i] = fLoGainSignal[i+1] - fLoGainSignal[i] - fLoGainSignal[i] + fLoGainSignal[i-1];
    607629      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
    608       p++;
    609630    }
    610631
    … …  
    616637  // Now find the maximum 
    617638  //
    618   Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
    619   Float_t lower = (Float_t)maxpos-1.;
    620   Float_t upper = (Float_t)maxpos;
    621   Float_t x     = lower;
    622   Float_t y     = 0.;
    623   Float_t a     = 1.;
    624   Float_t b     = 0.;
    625   Int_t   klo = maxpos-1;
    626   Int_t   khi = maxpos;
    627   Float_t klocont = fLoGainSignal[klo];
    628   Float_t khicont = fLoGainSignal[khi];
    629   fAbMax    = klocont;
    630   fAbMaxPos = lower;
    631 
    632 
    633   //
    634   // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
    635   // interval maxpos+1.
    636   //
    637   while (x<upper-0.1)
     639  Float_t step    = 0.2; // start with step size of 1ns and loop again with the smaller one
     640  Float_t lower   = (Float_t)maxpos-1.;
     641  Float_t upper   = (Float_t)maxpos;
     642  fAbMaxPos       = upper;
     643  Float_t x       = lower;
     644  Float_t y       = 0.;
     645  Float_t a       = 1.;
     646  Float_t b       = 0.;
     647  Int_t   klo     = maxpos-1;
     648  Int_t   khi     = maxpos;
     649
     650  //
     651  // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
     652  // If no maximum is found, go to interval maxpos+1.
     653  //
     654  while ( x < upper - 0.3 )
    638655    {
    639656
    … …  
    642659      b += step;
    643660
    644       y = a*klocont
    645         + b*khicont
     661      y = a*fLoGainSignal[klo]
     662        + b*fLoGainSignal[khi]
    646663        + (a*a*a-a)*fLoGainSecondDeriv[klo]
    647664        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    … …  
    649666      if (y > fAbMax)
    650667        {
    651           fAbMax   = y;
     668          fAbMax    = y;
    652669          fAbMaxPos = x;
    653670        }
    654671
    655     }
    656 
    657  if (fAbMaxPos < lower+0.1)
    658     {
    659 
    660       upper = (Float_t)maxpos-1.;
    661       lower = (Float_t)maxpos-2.;
    662       x     = lower;
    663       a     = 1.;
    664       b     = 0.;
    665       khi   = maxpos-1;
    666       klo   = maxpos-2;
    667       klocont = fLoGainSignal[klo];
    668       khicont = fLoGainSignal[khi];
    669 
    670       while (x<upper-0.1)
     672      //      *fLog << err << x << "  " << y << " " << fAbMaxPos<< endl;
     673    }
     674
     675  //
     676  // Test the possibility that the absolute maximum has not been found before the
     677  // maxpos and search from maxpos to maxpos+1 in steps of 0.2
     678  //
     679  if (fAbMaxPos > upper-0.1)
     680    {
     681
     682      upper   = (Float_t)maxpos+1.;
     683      lower   = (Float_t)maxpos;
     684      x       = lower;
     685      a       = 1.;
     686      b       = 0.;
     687      khi     = maxpos+1;
     688      klo     = maxpos;
     689
     690      while (x<upper-0.3)
    671691        {
    672692
    … …  
    675695          b += step;
    676696         
    677           y = a* klocont
    678             + b* khicont
     697          y = a*fLoGainSignal[klo]
     698            + b*fLoGainSignal[khi]
    679699            + (a*a*a-a)*fLoGainSecondDeriv[klo]
    680700            + (b*b*b-b)*fLoGainSecondDeriv[khi];
    681          
     701
    682702          if (y > fAbMax)
    683703            {
    … …  
    685705              fAbMaxPos = x;
    686706            }
     707          //          *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;
     708
    687709        }
    688710  }
    689  
    690  const Float_t up = fAbMaxPos+step-0.035;
    691  const Float_t lo = fAbMaxPos-step+0.035;
    692  const Float_t maxpossave = fAbMaxPos;
    693  
     711
     712
     713  //
     714  // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
     715  // Try a better precision.
     716  //
     717  const Float_t up = fAbMaxPos+step-0.055;
     718  const Float_t lo = fAbMaxPos-step+0.055;
     719  const Float_t maxpossave = fAbMaxPos;
     720 
    694721  x     = fAbMaxPos;
    695722  a     = upper - x;
    696723  b     = x - lower;
    697 
    698   step  = 0.025; // step size of 165 ps
    699 
     724 
     725  step  = 0.02; // step size of 42 ps
     726 
    700727  while (x<up)
    701728    {
    … …  
    705732      b += step;
    706733     
    707       y = a* klocont
    708         + b* khicont
     734      y = a*fLoGainSignal[klo]
     735        + b*fLoGainSignal[khi]
    709736        + (a*a*a-a)*fLoGainSecondDeriv[klo]
    710737        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    … …  
    715742          fAbMaxPos = x;
    716743        }
    717      
    718     }
    719 
     744      //      *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;     
     745    }
     746
     747  //
     748  // Second, try from time down to time-0.2 in steps of 0.04.
     749  //
    720750  x     = maxpossave;
     751
     752  //
     753  // Test the possibility that the absolute maximum has not been found between
     754  // maxpos and maxpos+0.02, then we have to look between maxpos-0.02 and maxpos
     755  // which requires new setting of klocont and khicont
     756  //
     757  if (x < klo + 0.02)
     758    {
     759      klo--;
     760      khi--;
     761      upper--;
     762      lower--;
     763    }
     764
    721765  a     = upper - x;
    722766  b     = x - lower;
    723 
     767 
    724768  while (x>lo)
    725769    {
    … …  
    729773      b -= step;
    730774     
    731       y = a* klocont
    732         + b* khicont
     775      y = a*fLoGainSignal[klo]
     776        + b*fLoGainSignal[khi]
    733777        + (a*a*a-a)*fLoGainSecondDeriv[klo]
    734778        + (b*b*b-b)*fLoGainSecondDeriv[khi];
    … …  
    739783          fAbMaxPos = x;
    740784        }
    741      
    742     }
    743 
    744   fHalfMax = fAbMax/2.;
    745 
    746   //
    747   // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
    748   // First, find the right FADC slice:
    749   //
    750   klo   = maxpos -1;
    751   while (klo > maxpos-4)
    752     {
    753       if (*(first+klo) < (Byte_t)fHalfMax)
    754         break;
    755       klo--;
    756     }
    757 
    758   //
    759   // Loop from the beginning of the slice upwards to reach the fHalfMax:
    760   // With means of bisection:
    761   //
    762   x     = (Float_t)klo;
    763   a     = 1.;
    764   b     = 0.;
    765   klocont = fLoGainSignal[klo];
    766   khicont = fLoGainSignal[klo+1];
    767   time  = x;
    768 
    769   step = 0.5;
    770   Bool_t back = kFALSE;
    771 
    772   while (step > fResolution)
    773     {
    774      
    775       if (back)
    776         {
    777           x -= step;
    778           a += step;
    779           b -= step;
    780         }
    781       else
    782         {
    783           x += step;
    784           a -= step;
    785           b += step;
    786         }
    787      
    788       y = a*klocont
    789         + b*khicont
    790         + (a*a*a-a)*fLoGainSecondDeriv[klo]
    791         + (b*b*b-b)*fLoGainSecondDeriv[khi];
    792 
    793       if (y >= fHalfMax)
    794         back = kTRUE;
    795       else
    796         back = kFALSE;
    797 
    798       step /= 2.;
    799     }
    800 
    801   time  = (Float_t)fLoGainFirst + x;
    802   dtime = fResolution;
    803 
    804   //
    805   // Now integrate the whole thing!
    806   //
    807   Int_t startslice = (Int_t)(time - 1.5);
    808   Int_t lastslice  = (Int_t)(fAbMaxPos + 2.*fRatioMax2Fall*fAbMax-1.);
    809   Int_t coverage   = lastslice-startslice;
    810 
    811   //
    812   // make them even
    813   //
    814   if (coverage != (coverage & ~1))
    815     lastslice++;
    816 
    817   if (startslice < 0)
    818     {
    819       //      gLog << warn << GetDescriptor()
    820       //           << ": Low Gain Bounce against left border, your range is too limited! " << endl;
    821       sum = 0.;
     785      //      *fLog << warn << x << "  " << y << "  " << fAbMaxPos << endl;
     786    }
     787
     788  if (IsTimeType(kMaximum))
     789    {
     790      time = (Float_t)fLoGainFirst + fAbMaxPos;
     791      dtime = 0.02;
     792    }
     793  else
     794    {
     795      fHalfMax = fAbMax/2.;
     796     
     797      //
     798      // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
     799      // First, find the right FADC slice:
     800      //
     801      klo  = maxpos - 1;
     802      while (klo >= 0)
     803        {
     804          if (fLoGainSignal[klo] < fHalfMax)
     805            break;
     806          klo--;
     807        }
     808     
     809      //
     810      // Loop from the beginning of the slice upwards to reach the fHalfMax:
     811      // With means of bisection:
     812      //
     813      x     = (Float_t)klo;
     814      a     = 1.;
     815      b     = 0.;
     816     
     817      step = 0.5;
     818      Bool_t back = kFALSE;
     819     
     820      Int_t maxcnt = 1000;
     821      Int_t cnt    = 0;
     822
     823      while (TMath::Abs(y-fHalfMax) > fResolution)
     824        {
     825         
     826          if (back)
     827            {
     828              x -= step;
     829              a += step;
     830              b -= step;
     831            }
     832          else
     833            {
     834              x += step;
     835              a -= step;
     836              b += step;
     837            }
     838         
     839          y = a*fLoGainSignal[klo]
     840            + b*fLoGainSignal[khi]
     841            + (a*a*a-a)*fLoGainSecondDeriv[klo]
     842            + (b*b*b-b)*fLoGainSecondDeriv[khi];
     843         
     844          if (y > fHalfMax)
     845            back = kTRUE;
     846          else
     847            back = kFALSE;
     848
     849          if (++cnt > maxcnt)
     850            {
     851              *fLog << inf << x << "  " << y << " " << fHalfMax << endl;
     852              break;
     853            }
     854         
     855          step /= 2.;
     856        }
     857
     858      time  = (Float_t)fLoGainFirst + x;
     859      dtime = fResolution;
     860    }
     861 
     862  if (IsChargeType(kAmplitude))
     863    {
     864      sum = fAbMax;
    822865      return;
    823       startslice = 0;
    824     }
    825 
    826   if (lastslice > range-1)
    827     {
    828       sum = 0.;
    829       return;
    830       //      gLog << warn << GetDescriptor()
    831       //           << ": Low Gain Bounce against right border, your range is too limited! " << endl;
    832       lastslice = range-1;
    833     }
    834  
    835   Int_t i = startslice;
    836   sum = 0.5*fLoGainSignal[i] + 0.75*fLoGainSecondDeriv[i];
    837 
    838   for (i=startslice+1; i<lastslice; i++)
    839     sum += fLoGainSignal[i] + 1.5*fLoGainSecondDeriv[i];
    840 
    841  
    842   sum += 0.5*fLoGainSignal[lastslice] + 0.75*fLoGainSecondDeriv[lastslice];
    843   //
    844   // subtract the pedestal
    845   //
    846   //  sum -= pedes*(lastslice-startslice);
    847  
    848   //  sig.SetNumLoGainSlices(fNumLoGainSamples);     
     866    }
     867
     868  if (IsChargeType(kIntegral))
     869    {
     870      //
     871      // Now integrate the whole thing!
     872      //
     873      Int_t startslice = (Int_t)(fAbMaxPos - fRiseTime);
     874      Int_t lastslice  = (Int_t)(fAbMaxPos + fFallTime);
     875     
     876      if (startslice < 0)
     877        {
     878          lastslice -= startslice;
     879          startslice = 0;
     880        }
     881
     882      Int_t i = startslice;
     883      sum = 0.5*fLoGainSignal[i];
     884     
     885      for (i=startslice+1; i<lastslice; i++)
     886        sum += fLoGainSignal[i] + 1.5*fLoGainSecondDeriv[i];
     887     
     888      sum += 0.5*fLoGainSignal[lastslice];
     889    }
     890 
    849891
    850892}

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

    r5205 r5213  
    1212private:
    1313 
    14   static const Byte_t  fgHiGainFirst;    // Default for fHiGainFirst  (now set to: 2)
    15   static const Byte_t  fgHiGainLast;     // Default for fHiGainLast   (now set to: 14)
    16   static const Byte_t  fgLoGainFirst;    // Default for fLOGainFirst  (now set to: 3)
    17   static const Byte_t  fgLoGainLast;     // Default for fLoGainLast   (now set to: 14)
    18   static const Float_t fgResolution;     // Default for fResolution   (now set to: 0.003)
    19   static const Float_t fgRatioMax2Fall;  // Default ratio signal maximum to fall time (now set to: 0.03)
     14  static const Byte_t  fgHiGainFirst;    //! Default for fHiGainFirst  (now set to: 2)
     15  static const Byte_t  fgHiGainLast;     //! Default for fHiGainLast   (now set to: 14)
     16  static const Byte_t  fgLoGainFirst;    //! Default for fLOGainFirst  (now set to: 3)
     17  static const Byte_t  fgLoGainLast;     //! Default for fLoGainLast   (now set to: 14)
     18  static const Float_t fgResolution;     //! Default for fResolution   (now set to: 0.003)
     19  static const Float_t fgRiseTime;       //! Default for fRiseTime     (now set to: 1.5)
     20  static const Float_t fgFallTime;       //! Default for fFallTime     (now set to: 4.5)
    2021
    21   Float_t *fHiGainSignal;               //! Need fast access to the signals in a float way
    22   Float_t *fLoGainSignal;               //! Store them in separate arrays
    23   Float_t *fHiGainFirstDeriv;           //!
    24   Float_t *fLoGainFirstDeriv;           //!
    25   Float_t *fHiGainSecondDeriv;          //!
    26   Float_t *fLoGainSecondDeriv;          //!
     22  Float_t *fHiGainSignal;                //! Need fast access to the signals in a float way
     23  Float_t *fLoGainSignal;                //! Store them in separate arrays
     24  Float_t *fHiGainFirstDeriv;            //!
     25  Float_t *fLoGainFirstDeriv;            //!
     26  Float_t *fHiGainSecondDeriv;           //!
     27  Float_t *fLoGainSecondDeriv;           //!
    2728
    28   Float_t fResolution;                  // The time resolution in FADC units
    29   Float_t fRatioMax2Fall;               // Relation fall time to maximum slice content
     29  Float_t fResolution;                   // The time resolution in FADC units
     30  Float_t fRiseTime;                     // The usual rise time of the pulse
     31  Float_t fFallTime;                     // The usual fall time of the pulse
    3032
    31   Float_t fAbMax;
    32   Float_t fAbMaxPos;
    33   Float_t fHalfMax;
     33  Float_t fAbMax;                        // Current maximum of the spline
     34  Float_t fAbMaxPos;                     // Current position of the maximum of the spline
     35  Float_t fHalfMax;                      // Current half maximum of the spline
     36
     37  Byte_t  fTimeFlags;                    // Flag to hold the time extraction type
     38  Byte_t  fChargeFlags;                  // Flag to hold the charge extraction type
    3439 
    3540  Int_t  PreProcess(MParList *pList);
    … …  
    4550public:
    4651
     52  enum TimeType_t   { kMaximum,   kHalfMaximum };    //! Possible time   extraction types
     53  enum ChargeType_t { kAmplitude, kIntegral    };    //! Possible charge extraction types
     54
     55private:
     56
     57  Bool_t IsTimeType  ( TimeType_t   typ )  { return TESTBIT(fTimeFlags  , typ); }
     58  Bool_t IsChargeType( ChargeType_t typ )  { return TESTBIT(fChargeFlags, typ); } 
     59 
     60public:
     61
    4762  MExtractTimeAndChargeSpline(const char *name=NULL, const char *title=NULL);
    4863  ~MExtractTimeAndChargeSpline(); 
    49    
    50   void SetResolution     ( Float_t f=fgResolution    )  { fResolution    = f;  }
    51   void SetRatioMax2Fall  ( Float_t f=fgRatioMax2Fall )  { fRatioMax2Fall = f;  }
     64
     65  void SetResolution   ( Float_t f=fgResolution  )  { fResolution  = f;  }
     66  void SetRiseTime     ( Float_t f=fgRiseTime    )  { fRiseTime    = f;  }
     67  void SetFallTime     ( Float_t f=fgFallTime    )  { fFallTime    = f;  }
     68
     69  void SetTimeType     ( TimeType_t typ=kMaximum    ) { fTimeFlags = 0;   SETBIT(fTimeFlags,typ);  }
     70  void SetChargeType   ( ChargeType_t typ=kAmplitude) { fChargeFlags = 0; SETBIT(fChargeFlags,typ);}
    5271 
    53   ClassDef(MExtractTimeAndChargeSpline, 0)   // Task to Extract the Arrival Times using a Fast Spline
     72  ClassDef(MExtractTimeAndChargeSpline, 0)   // Task to Extract Arrival Times and Charges using a Fast Cubic Spline
    5473};
    5574
Note: See TracChangeset for help on using the changeset viewer.