Ignore:
Timestamp:
01/08/05 20:17:34 (20 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/msignal
Files:
4 edited

Legend:

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

    r5745 r5747  
    252252
    253253      if (*p++ >= fSaturationLimit)
    254         sat++;
     254        if (!sat)
     255          sat = ids-3;
    255256    }
    256257
     
    272273         
    273274          if (*logain++ >= fSaturationLimit)
    274             sat++;
     275            if (!sat)
     276              sat = ids-3;
    275277
    276278          range++;
  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc

    r5746 r5747  
    231231      if (*p++ >= fSaturationLimit)
    232232        if (!sat)
    233             sat = p-first+fHiGainFirst-1;
     233            sat = ids-2;
    234234     
    235235      count++;
     
    242242    if (*p++ >= fSaturationLimit)
    243243      if (!sat)
    244         sat = p-first+fHiGainFirst-1;
     244        sat = ids-2;
    245245 
    246246  if (IsNoiseCalculation())
     
    287287          if (*l++ >= fSaturationLimit)
    288288            if (!sat)
    289               sat = l-logain+fHiGainLast-1;
     289              sat = ids-2;
    290290         
    291291          if (sum>max)
     
    308308              if (*l++ >= fSaturationLimit)
    309309                if (!sat)
    310                   sat = l-logain+fHiGainLast-1;
     310                  sat = ids-2;
    311311
    312312              if (sum>max)
  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

    r5613 r5747  
    6666//   which simplifies the Numerical Recipes algorithm.
    6767//   (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.)
    68 //
    6968//
    7069//   The algorithm to search the time proceeds as follows:
     
    102101//     (Int_t)(fAbMaxPos + fFallTime)
    103102//     (default: fRiseTime: 1.5, fFallTime: 4.5)
    104 //  3) Sum only half the values of the edge slices
    105 //  4) Sum 1.5*fHiGainSecondDeriv of the not-edge slices using the "natural cubic
    106 //     spline with second derivatives set to 0. at the edges.
    107 //     (Remember that fHiGainSecondDeriv had been divided by 6.)
    108103//   
    109104//  The values: fNumHiGainSamples and fNumLoGainSamples are set to:
    110105//  1) If Charge Type: kAmplitude was chosen: 1.
    111 //  2) If Charge Type: kIntegral was chosen: TMath::Floor(fRiseTime + fFallTime)
    112 //                 or: TMath::Floor(fRiseTime + fFallTime + 1.) in the case of the low-gain
     106//  2) If Charge Type: kIntegral was chosen: fRiseTime + fFallTime
     107//                 or: fRiseTime + fFallTime + 1. in the case of the low-gain
    113108//
    114109//  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast)
     
    127122//        for the case, the extraction type kIntegral has been chosen.
    128123//
    129 //  Call: - SetTimeType(MExtractTimeAndChargeSpline::kMaximum) for extraction
     124//  Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the
     125//          computation of the amplitude at the maximum (default) and extraction
    130126//          the position of the maximum (default)
    131 //          --> needs 22 function evaluations
    132 //        - SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum) for extraction
    133 //          the position of the half maximum at the rising edge.
    134 //          --> needs max. 34 function evaluations
    135 //        - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the
    136 //          computation of the amplitude at the maximum (default)
    137127//          --> no further function evaluation needed
    138128//        - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the
    139129//          computation of the integral beneith the spline between fRiseTime
    140130//          from the position of the maximum to fFallTime after the position of
    141 //          the maximum. The Low Gain is computed with one more slice at the falling
    142 //          edge.
    143 //          --> needs one more simple summation loop over 7 slices.
     131//          the maximum. The Low Gain is computed with half a slice more at the rising
     132//          edge and half a slice more at the falling edge.
     133//          The time of the half maximum is returned.
     134//          --> needs one function evaluations but is more precise
    144135//       
    145136//////////////////////////////////////////////////////////////////////////////
     
    159150const Byte_t  MExtractTimeAndChargeSpline::fgLoGainFirst  = 2;
    160151const Byte_t  MExtractTimeAndChargeSpline::fgLoGainLast   = 14;
    161 const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.025;
     152const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.05;
    162153const Float_t MExtractTimeAndChargeSpline::fgRiseTime     = 1.5;
    163154const Float_t MExtractTimeAndChargeSpline::fgFallTime     = 4.5;
     
    177168//
    178169MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title)
    179     : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), fRandomIter(0)
     170    : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), fHalfMaxPos(0.), fRandomIter(0)
    180171{
    181172
     
    189180  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
    190181
    191   SetTimeType();
    192182  SetChargeType();
    193183
     
    212202    SetChargeType(kAmplitude);
    213203 
    214 }
    215 
    216 
    217 //-------------------------------------------------------------------
    218 //
    219 // Set the Time Extraction type. Possible are:
    220 // - kMaximum:     Search the maximum of the spline and return its position
    221 // - kHalfMaximum: Search the half maximum left from the maximum and return
    222 //                 its position
    223 //
    224 void MExtractTimeAndChargeSpline::SetTimeType( ExtractionType_t typ )
    225 {
    226 
    227   CLRBIT(fFlags,kMaximum);
    228   CLRBIT(fFlags,kHalfMaximum);
    229   SETBIT(fFlags,typ);
    230 
    231204}
    232205
     
    316289  const Byte_t *end = first + range;
    317290  Byte_t       *p  = first;
    318   Int_t  count     = 0;
     291
     292  sat = 0;
    319293
    320294  const Float_t pedes  = ped.GetPedestal();
    321295  const Float_t ABoffs = ped.GetPedestalABoffset();
    322296
    323   Float_t pedmean[2];
    324   pedmean[0] = pedes + ABoffs;
    325   pedmean[1] = pedes - ABoffs;
     297  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    326298
    327299  fAbMax        = 0.;
    328300  fAbMaxPos     = 0.;
     301  fHalfMax      = 0.;
     302  fHalfMaxPos   = 0.;
    329303  Int_t  maxpos = 0;
     304  Int_t  max    = 0;
    330305
    331306  //
    332307  // Check for saturation in all other slices
    333308  //
     309  Int_t ids = fHiGainFirst;
     310  Float_t *sample = fHiGainSignal.GetArray();
    334311  while (p<end)
    335312    {
    336313
    337       const Int_t ids      = fHiGainFirst + count ;
    338       const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
    339       fHiGainSignal[count] = signal;
    340 
    341       if (signal > fAbMax + 0.1) /* the 0.1 is necessary for the ultra-high enery events saturating many slices */
    342         {
    343           fAbMax = signal;
    344           maxpos = count;
     314      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
     315
     316      if (*p > max)
     317        {
     318          maxpos = ids-fHiGainFirst-1;
     319          max    = *p;
    345320        }
    346321
    347322      if (*p++ >= fSaturationLimit)
    348             sat++;
    349      
    350       count++;                                                   
     323        if (!sat)
     324          sat = ids-3;
     325     
    351326    }
    352327 
     
    359334        {
    360335
    361           const Int_t ids      = fHiGainFirst + range ;
    362           const Float_t signal = (Float_t)*logain - pedmean[(ids+abflag) & 0x1];
    363           fHiGainSignal[range] = signal;
    364          
    365           if (signal > fAbMax)
     336          *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
     337
     338          if (*logain > max)
    366339            {
    367               fAbMax = signal;
    368               maxpos = range;
     340              maxpos = ids-fHiGainFirst-1;
     341              max    = *logain;
    369342            }
    370343         
    371           if (*logain >= fSaturationLimit)
    372             sat++;
     344          if (*logain++ >= fSaturationLimit)
     345            if (!sat)
     346              sat = ids-3;
    373347
    374348          range++;
    375           logain++;
    376         }
    377     }
    378  
     349        }
     350    }
     351 
     352  fAbMax = fHiGainSignal[maxpos];
     353
    379354  Float_t pp;
    380355
     
    434409
    435410  //
    436   // Allow one saturated slice
     411  // Allow no saturated slice
    437412  // and
    438413  // Don't start if the maxpos is too close to the limits.
    439414  //
    440   if (sat > 1 || maxpos < 1 || maxpos > range-2)
    441     {
    442       time =  IsExtractionType(kMaximum)
    443         ? (Float_t)(fHiGainFirst + maxpos)
    444         : (Float_t)(fHiGainFirst + maxpos - 1);
    445 
     415  if (sat || maxpos < 1 || maxpos > range-2)
     416    {
     417      dtime = 0.5;
    446418      if (IsExtractionType(kAmplitude))
    447419        {
    448           sum = fAbMax;
     420          sum  = fAbMax;
     421          time = (Float_t)(fHiGainFirst + maxpos);
    449422          return;
    450423        }
     
    455428        CalcIntegralHiGain(sum, 0.001, fRiseTime + fFallTime);
    456429
     430      time =  (Float_t)(fHiGainFirst + maxpos - 1);
    457431      return;
    458432    }
    459433     
     434  dtime = fResolution;
     435
    460436  //
    461437  // Now find the maximum 
     
    533509  // Try a better precision.
    534510  //
    535   const Float_t up = fAbMaxPos+step - 1.5*fResolution;
    536   const Float_t lo = fAbMaxPos-step + 1.5*fResolution;
     511  const Float_t up = fAbMaxPos+step - 3.0*fResolution;
     512  const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
    537513  const Float_t maxpossave = fAbMaxPos;
    538514 
     
    541517  b     = x - lower;
    542518 
    543   step  = fResolution; // step size of 83 ps
     519  step  = 2.*fResolution; // step size of 0.1 FADC slices
    544520 
    545521  while (x<up)
     
    569545  //
    570546  // Test the possibility that the absolute maximum has not been found between
    571   // maxpos and maxpos+0.025, then we have to look between maxpos-0.025 and maxpos
     547  // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
    572548  // which requires new setting of klocont and khicont
    573549  //
    574   if (x < lower + fResolution/2.)
     550  if (x < lower + fResolution)
    575551    {
    576552      klo--;
     
    602578    }
    603579
    604   if (IsExtractionType(kMaximum))
    605     {
    606       time = (Float_t)fHiGainFirst + fAbMaxPos;
    607       dtime = fResolution;
    608     }
    609   else
    610     {
    611       fHalfMax = fAbMax/2.;
    612      
    613       //
    614       // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
    615       // First, find the right FADC slice:
    616       //
    617       klo  = maxpos;
    618       while (klo >= 0)
    619         {
    620           klo--;
    621           if (fHiGainSignal[klo] < fHalfMax)
    622             break;
    623         }
    624      
    625       khi = klo+1;
    626       //
    627       // Loop from the beginning of the slice upwards to reach the fHalfMax:
    628       // With means of bisection:
    629       //
    630       x     = (Float_t)klo;
    631       a     = 1.;
    632       b     = 0.;
    633      
    634       step = 0.5;
    635       Bool_t back = kFALSE;
    636      
    637       Int_t maxcnt = 20;
    638       Int_t cnt    = 0;
    639 
    640       while (TMath::Abs(y-fHalfMax) > fResolution)
    641         {
    642          
    643           if (back)
    644             {
    645               x -= step;
    646               a += step;
    647               b -= step;
    648             }
    649           else
    650             {
    651               x += step;
    652               a -= step;
    653               b += step;
    654             }
    655          
    656           y = a*fHiGainSignal[klo]
    657             + b*fHiGainSignal[khi]
    658             + (a*a*a-a)*fHiGainSecondDeriv[klo]
    659             + (b*b*b-b)*fHiGainSecondDeriv[khi];
    660          
    661           if (y > fHalfMax)
    662             back = kTRUE;
    663           else
    664             back = kFALSE;
    665 
    666           if (++cnt > maxcnt)
    667             break;
    668          
    669           step /= 2.;
    670         }
    671 
    672       time  = (Float_t)fHiGainFirst + x;
    673       dtime = fResolution;
    674     }
    675  
    676580  if (IsExtractionType(kAmplitude))
    677581    {
    678       sum = fAbMax;
     582      time  = fAbMaxPos + (Int_t)fHiGainFirst;
     583      sum   = fAbMax;
    679584      return;
    680585    }
    681586
    682   if (IsExtractionType(kIntegral))
    683     {
    684       //
    685       // Now integrate the whole thing!
    686       //
    687 
    688       Float_t start = fAbMaxPos - fRiseTime;
    689       Float_t last  = fAbMaxPos + fFallTime;
    690      
    691       const Int_t diff = int(last) - range;
    692 
    693       if (diff > 0)
    694         {
    695           last  -= diff;
    696           start -= diff;
    697         }
    698 
    699       CalcIntegralHiGain(sum, start, last);
    700     }
     587  fHalfMax = fAbMax/2.;
     588 
     589  //
     590  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
     591  // First, find the right FADC slice:
     592  //
     593  klo  = maxpos;
     594  while (klo >= 0)
     595    {
     596      klo--;
     597      if (fHiGainSignal[klo] < fHalfMax)
     598        break;
     599    }
     600 
     601  khi = klo+1;
     602  //
     603  // Loop from the beginning of the slice upwards to reach the fHalfMax:
     604  // With means of bisection:
     605  //
     606  x     = (Float_t)klo;
     607  a     = 1.;
     608  b     = 0.;
     609 
     610  step = 0.5;
     611  Bool_t back = kFALSE;
     612 
     613  Int_t maxcnt = 20;
     614  Int_t cnt    = 0;
     615 
     616  while (TMath::Abs(y-fHalfMax) > fResolution)
     617    {
     618     
     619      if (back)
     620        {
     621          x -= step;
     622          a += step;
     623          b -= step;
     624        }
     625      else
     626        {
     627          x += step;
     628          a -= step;
     629          b += step;
     630        }
     631     
     632      y = a*fHiGainSignal[klo]
     633        + b*fHiGainSignal[khi]
     634        + (a*a*a-a)*fHiGainSecondDeriv[klo]
     635        + (b*b*b-b)*fHiGainSecondDeriv[khi];
     636     
     637      if (y > fHalfMax)
     638        back = kTRUE;
     639      else
     640        back = kFALSE;
     641     
     642      if (++cnt > maxcnt)
     643        break;
     644     
     645      step /= 2.;
     646    }
     647 
     648  time  = (Float_t)fHiGainFirst + x;
     649  //
     650  // Now integrate the whole thing!
     651  //
     652 
     653  Float_t start = fAbMaxPos - fRiseTime;
     654  Float_t last  = fAbMaxPos + fFallTime;
     655 
     656  const Int_t diff = int(last) - range;
     657 
     658  if (diff > 0)
     659    {
     660      last  -= diff;
     661      start -= diff;
     662    }
     663 
     664  CalcIntegralHiGain(sum, start, last);
    701665}
    702666
     
    718682  const Float_t ABoffs = ped.GetPedestalABoffset();
    719683
    720   Float_t pedmean[2];
    721   pedmean[0] = pedes + ABoffs;
    722   pedmean[1] = pedes - ABoffs;
     684  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    723685
    724686  fAbMax        = 0.;
    725687  fAbMaxPos     = 0.;
    726688  Int_t  maxpos = 0;
    727   Int_t  count  = 0;
     689  Int_t  max    = 0;
    728690
    729691  //
    730692  // Check for saturation in all other slices
    731693  //
     694  Int_t    ids    = fLoGainFirst;
     695  Float_t *sample = fLoGainSignal.GetArray();
    732696  while (p<end)
    733697    {
    734698
    735       const Int_t ids      =  count + fLoGainFirst;
    736       const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
    737       fLoGainSignal[count] = signal;
    738 
    739       if (signal > fAbMax + 0.1)
    740         {
    741           fAbMax = signal;
    742           maxpos = count;
     699      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
     700
     701      if (*p > max)
     702        {
     703          maxpos = ids-fLoGainFirst-1;
     704          max    = *p;
    743705        }
    744706
    745707      if (*p++ >= fSaturationLimit)
    746708        sat++;
    747      
    748       count++;                                                   
    749     }
     709    }
     710 
     711  fAbMax = fLoGainSignal[maxpos];
    750712 
    751713  Float_t pp;
     
    808770  // Don't start if the maxpos is too close to the limits.
    809771  //
    810   if (sat || maxpos < 2 || maxpos > range-2)
    811     {
    812       time =  IsExtractionType(kMaximum)
    813         ? (Float_t)(fLoGainFirst + maxpos)
    814         : (Float_t)(fLoGainFirst + maxpos - 1);
    815      
     772  if (sat || maxpos < TMath::Ceil(fRiseTime+0.45) || maxpos > range-2)
     773    {
     774      dtime = 0.5;
    816775      if (IsExtractionType(kAmplitude))
    817776        {
     777          time = (Float_t)(fLoGainFirst + maxpos);
    818778          sum = fAbMax;
    819779          return;
     
    825785        CalcIntegralLoGain(sum, 0.001, fRiseTime + fFallTime + 1.);
    826786
     787      time = (Float_t)(fLoGainFirst + maxpos - 1);
    827788      return;
    828789    }
    829790
    830   if (maxpos < (Int_t)(fRiseTime+2.))
    831     {
    832       time =  IsExtractionType(kMaximum)
    833         ? (Float_t)(fLoGainFirst + maxpos)
    834         : (Float_t)(fLoGainFirst + maxpos - 1);
    835 
    836       if (maxpos > range-2)
    837         CalcIntegralLoGain(sum, (Float_t)range - fRiseTime - fFallTime-1., (Float_t)range - 0.001);
    838       else
    839         CalcIntegralLoGain(sum, 0.001, fRiseTime + fFallTime + 1.);
    840 
    841       return;
    842     }
     791  dtime = fResolution;
    843792
    844793  //
     
    920869  // Try a better precision.
    921870  //
    922   const Float_t up = fAbMaxPos+step - 1.5*fResolution;
    923   const Float_t lo = fAbMaxPos-step + 1.5*fResolution;
     871  const Float_t up = fAbMaxPos+step - 3.0*fResolution;
     872  const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
    924873  const Float_t maxpossave = fAbMaxPos;
    925874 
     
    928877  b     = x - lower;
    929878 
    930   step  = fResolution; // step size of fResolution (33 ps )
     879  step  = 2.*fResolution; // step size of 0.1 FADC slice
    931880 
    932881  while (x<up)
     
    956905  //
    957906  // Test the possibility that the absolute maximum has not been found between
    958   // maxpos and maxpos+0.02, then we have to look between maxpos-0.02 and maxpos
     907  // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
    959908  // which requires new setting of klocont and khicont
    960909  //
    961   if (x < lower + fResolution/2.)
     910  if (x < lower + fResolution)
    962911    {
    963912      klo--;
     
    989938    }
    990939
    991   if (IsExtractionType(kMaximum))
     940  if (IsExtractionType(kAmplitude))
    992941    {
    993942      time = fAbMaxPos + (Int_t)fLoGainFirst;
    994       dtime = fResolution;
    995     }
    996   else
    997     {
    998       fHalfMax = fAbMax/2.;
    999      
    1000       //
    1001       // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
    1002       // First, find the right FADC slice:
    1003       //
    1004       klo  = maxpos;
    1005       while (klo > 0)
    1006         {
    1007           klo--;
    1008           if (fLoGainSignal[klo] < fHalfMax)
    1009             break;
    1010         }
    1011      
    1012       khi = klo+1;
    1013       //
    1014       // Loop from the beginning of the slice upwards to reach the fHalfMax:
    1015       // With means of bisection:
    1016       //
    1017       x     = (Float_t)klo;
    1018       a     = 1.;
    1019       b     = 0.;
    1020      
    1021       step = 0.5;
    1022       Bool_t back = kFALSE;
    1023      
    1024       Int_t maxcnt = 20;
    1025       Int_t cnt    = 0;
    1026 
    1027       while (TMath::Abs(y-fHalfMax) > fResolution)
    1028         {
    1029          
    1030           if (back)
    1031             {
    1032               x -= step;
    1033               a += step;
    1034               b -= step;
    1035             }
    1036           else
    1037             {
    1038               x += step;
    1039               a -= step;
    1040               b += step;
    1041             }
    1042          
    1043           y = a*fLoGainSignal[klo]
    1044             + b*fLoGainSignal[khi]
    1045             + (a*a*a-a)*fLoGainSecondDeriv[klo]
    1046             + (b*b*b-b)*fLoGainSecondDeriv[khi];
    1047          
    1048           if (y > fHalfMax)
    1049             back = kTRUE;
    1050           else
    1051             back = kFALSE;
    1052 
    1053           if (++cnt > maxcnt)
    1054             break;
    1055          
    1056           step /= 2.;
    1057         }
    1058 
    1059       time  = x + (Int_t)fLoGainFirst;
    1060       dtime = fResolution;
    1061     }
    1062  
    1063   if (IsExtractionType(kAmplitude))
    1064     {
    1065       sum = fAbMax;
     943      sum  = fAbMax;
    1066944      return;
    1067945    }
    1068 
    1069   if (IsExtractionType(kIntegral))
    1070     {
    1071       //
    1072       // Now integrate the whole thing!
    1073       //
    1074       Float_t start = fAbMaxPos - fRiseTime - 0.5;
    1075       Float_t last  = fAbMaxPos + fFallTime + 0.5;
    1076      
    1077       const Int_t diff = int(last) - range;
    1078 
    1079       if (diff > 0)
    1080         {
    1081           last  -= diff;
    1082           start -= diff;
    1083         }
    1084       CalcIntegralLoGain(sum, start, last);
    1085       //      *fLog << inf << time << "  " << sum << "  " << start << "  " << last << endl;
    1086     }
    1087 
     946 
     947  fHalfMax = fAbMax/2.;
     948 
     949  //
     950  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
     951  // First, find the right FADC slice:
     952  //
     953  klo  = maxpos;
     954  while (klo > 0)
     955    {
     956      klo--;
     957      if (fLoGainSignal[klo] < fHalfMax)
     958        break;
     959    }
     960 
     961  khi = klo+1;
     962  //
     963  // Loop from the beginning of the slice upwards to reach the fHalfMax:
     964  // With means of bisection:
     965  //
     966  x     = (Float_t)klo;
     967  a     = 1.;
     968  b     = 0.;
     969 
     970  step = 0.5;
     971  Bool_t back = kFALSE;
     972 
     973  Int_t maxcnt = 20;
     974  Int_t cnt    = 0;
     975 
     976  while (TMath::Abs(y-fHalfMax) > fResolution)
     977    {
     978     
     979      if (back)
     980        {
     981          x -= step;
     982          a += step;
     983          b -= step;
     984        }
     985      else
     986        {
     987          x += step;
     988          a -= step;
     989          b += step;
     990        }
     991     
     992      y = a*fLoGainSignal[klo]
     993        + b*fLoGainSignal[khi]
     994        + (a*a*a-a)*fLoGainSecondDeriv[klo]
     995        + (b*b*b-b)*fLoGainSecondDeriv[khi];
     996     
     997      if (y > fHalfMax)
     998        back = kTRUE;
     999      else
     1000        back = kFALSE;
     1001     
     1002      if (++cnt > maxcnt)
     1003        break;
     1004     
     1005      step /= 2.;
     1006    }
     1007 
     1008  time  = x + (Int_t)fLoGainFirst;
     1009
     1010  //
     1011  // Now integrate the whole thing!
     1012  //
     1013  Float_t start = fAbMaxPos - fRiseTime - 0.5;
     1014  Float_t last  = fAbMaxPos + fFallTime + 0.5;
     1015 
     1016  const Int_t diff = int(last) - range;
     1017 
     1018  if (diff > 0)
     1019    {
     1020      last  -= diff;
     1021      start -= diff;
     1022    }
     1023  CalcIntegralLoGain(sum, start, last);
    10881024}
    10891025
     
    12481184      rc = kTRUE;
    12491185    }
    1250   if (IsEnvDefined(env, prefix, "Maximum", print))
    1251     {
    1252       b = GetEnvValue(env, prefix, "Maximum", IsExtractionType(kMaximum));
    1253       if (b)
    1254         SetTimeType(kMaximum);
    1255       rc = kTRUE;
    1256     }
    1257   if (IsEnvDefined(env, prefix, "HalfMaximum", print))
    1258     {
    1259       b = GetEnvValue(env, prefix, "HalfMaximum", IsExtractionType(kHalfMaximum));
    1260       if (b)
    1261         SetTimeType(kHalfMaximum);
    1262       rc = kTRUE;
    1263     }
    12641186
    12651187  return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

    r5684 r5747  
    7272  void SetFallTime   ( const Float_t f=fgFallTime    )  { fFallTime    = f;  }
    7373
    74   void SetTimeType   ( const ExtractionType_t typ=kHalfMaximum );
    7574  void SetChargeType ( const ExtractionType_t typ=kAmplitude);
    7675 
Note: See TracChangeset for help on using the changeset viewer.