Changeset 8546 for trunk


Ignore:
Timestamp:
06/11/07 16:49:42 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/NEWS

    r8518 r8546  
    77     Also the old values seemed not exactly the PulsePos used for
    88     teh check.
     9
     10   - callisto: improved calculation of spline coefficients a lot. This
     11     leads to a further improvement of the event rate calibrating MUX
     12     data of about 15% (175evt/s instead of 150evt/s)
    913
    1014

  • trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

    r8545 r8546  
    193193Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
    194194{
    195 /*
    196     // The number of steps is calculated directly from the integration
    197     // window. This is the only way to ensure we are not dealing with
    198     // numerical rounding uncertanties, because we always get the same
    199     // value under the same conditions -- it might still be different on
    200     // other machines!
    201     const Float_t start = pos-fRiseTime;
    202     const Float_t step  = 0.2;
    203     const Float_t width = fRiseTime+fFallTime;
    204     const Float_t max   = fNum-1 - (width+step);
    205     const Int_t   num   = TMath::Nint(width/step);
    206 
    207     // The order is important. In some cases (limlo-/limup-check) it can
    208     // happen that max<0. In this case we start at 0
    209     if (start > max)
    210         start = max;
    211     if (start < 0)
    212         start = 0;
    213 
    214     start += step/2;
    215 
    216     Double_t sum = 0.;
    217     for (Int_t i=0; i<num; i++)
    218     {
    219         // Note: if x is close to one integer number (= a FADC sample)
    220         // we get the same result by using that sample as klo, and the
    221         // next one as khi, or using the sample as khi and the previous
    222         // one as klo (the spline is of course continuous). So we do not
    223         // expect problems from rounding issues in the argument of
    224         // Floor() above (we have noticed differences in roundings
    225         // depending on the compilation options).
    226 
    227         sum += EvalAt(start + i*step);
    228 
    229         // FIXME? Perhaps the integral should be done analitically
    230         // between every two FADC slices, instead of numerically
    231     }
    232     sum *= step; // Transform sum in integral
    233 
    234     return sum;
    235     */
    236 
    237195    // In the future we will calculate the intgeral analytically.
    238196    // It has been tested that it gives identical results within
    … …  
    280238    if (fNum<2)
    281239        return;
    282 /*
    283     //
    284     // Allow no saturated slice and
    285     // Don't start if the maxpos is too close to the limits.
    286     //
    287 
    288     const Bool_t limlo = maxbin <      TMath::Ceil(fRiseTime);
    289     const Bool_t limup = maxbin > fNum-TMath::Ceil(fFallTime)-1;
    290     if (sat || limlo || limup)
    291     {
    292         fTimeDev = 1.0;
    293         if (fExtractionType == kAmplitude)
    294         {
    295             fSignal    = fVal[maxbin];
    296             fTime      = maxbin;
    297             fSignalDev = 0;  // means: is valid
    298             return;
    299         }
    300 
    301         fSignal    = CalcIntegral(limlo ? 0 : fNum);
    302         fTime      = maxbin - 1;
    303         fSignalDev = 0;  // means: is valid
    304         return;
    305     }
    306 
    307     //
    308     // Now find the maximum
    309     //
    310 
    311     Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
    312 
    313     Int_t klo = maxbin-1;
    314 
    315     Float_t maxpos = maxbin;//! Current position of the maximum of the spline
    316     Float_t max = fVal[maxbin];//! Current maximum of the spline
    317 
    318     //
    319     // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
    320     // If no maximum is found, go to interval maxpos+1.
    321     //
    322     for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++)
    323     {
    324         const Float_t x = klo + step*(i+1);
    325         //const Float_t y = Eval(klo, x);
    326         const Float_t y = Eval(klo, x-klo);
    327         if (y > max)
    328         {
    329             max    = y;
    330             maxpos = x;
    331         }
    332     }
    333 
    334     //
    335     // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2
    336     //
    337     if (maxpos > maxbin - 0.1)
    338     {
    339         klo = maxbin;
    340 
    341         for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++)
    342         {
    343             const Float_t x = klo + step*(i+1);
    344             //const Float_t y = Eval(klo, x);
    345             const Float_t y = Eval(klo, x-klo);
    346             if (y > max)
    347             {
    348                 max    = y;
    349                 maxpos = x;
    350             }
    351         }
    352     }
    353 
    354     //
    355     // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
    356     // Try a better precision.
    357     //
    358     const Float_t up = maxpos+step - 3.0*fResolution;
    359     const Float_t lo = maxpos-step + 3.0*fResolution;
    360     const Float_t abmaxpos = maxpos;
    361 
    362     step  = 2.*fResolution; // step size of 0.1 FADC slices
    363 
    364     for (int i=0; i<TMath::Nint(TMath::Ceil((up-abmaxpos)/step)); i++)
    365     {
    366         const Float_t x = abmaxpos + (i+1)*step;
    367         //const Float_t y = Eval(klo, x);
    368         const Float_t y = Eval(klo, x-klo);
    369         if (y > max)
    370         {
    371             max    = y;
    372             maxpos = x;
    373         }
    374     }
    375 
    376     //
    377     // Second, try from time down to time-0.2 in steps of fResolution.
    378     //
    379 
    380     //
    381     // Test the possibility that the absolute maximum has not been found between
    382     // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
    383     // which requires new setting of klocont and khicont
    384     //
    385     if (abmaxpos < klo + fResolution)
    386         klo--;
    387 
    388     for (int i=TMath::Nint(TMath::Ceil((abmaxpos-lo)/step))-1; i>=0; i--)
    389     {
    390         const Float_t x = abmaxpos - (i+1)*step;
    391         //const Float_t y = Eval(klo, x);
    392         const Float_t y = Eval(klo, x-klo);
    393         if (y > max)
    394         {
    395             max    = y;
    396             maxpos = x;
    397         }
    398     }
    399 
    400     fTime      = maxpos;
    401     fTimeDev   = fResolution;
    402     fSignal    = CalcIntegral(maxpos);
    403     fSignalDev = 0;  // means: is valid
    404 
    405     return;
    406 */
    407     // --- Start NEW ---
    408 
    409     // This block extracts values very similar to the old algorithm...
    410     // for max>10
    411     /* Most accurate (old identical) version:
    412 
    413     Float_t xmax=maxpos, ymax=Eval(maxpos-1, 1);
    414     Int_t rc = GetMaxPos(maxpos-1, xmax, ymax);
    415     if (xmax==maxpos)
    416     {
    417         GetMaxPos(maxpos, xmax, ymax);
    418 
    419         Float_t y = Eval(maxpos, 1);
    420         if (y>ymax)
    421         {
    422             ymax = y;
    423             xmax = maxpos+1;
    424         }
    425     }*/
    426240
    427241    Float_t maxpos;
    … …  
    461275        fWidthDev = 0;
    462276    }
    463 
    464     //
    465     // Loop from the beginning of the slice upwards to reach the maxhalf:
    466     // With means of bisection:
    467     //
    468     /*
    469     static const Float_t sqrt2 = TMath::Sqrt(2.);
    470 
    471     step = sqrt2*3*0.061;//fRiseTime;
    472     Float_t x = maxpos-0.86-3*0.061;//fRiseTime*1.25;
    473 
    474 //    step = sqrt2*0.5;//fRiseTime;
    475 //    Float_t x = maxpos-1.25;//fRiseTime*1.25;
    476 
    477     Int_t  cnt  =0;
    478     while (cnt++<30)
    479     {
    480         const Float_t y=EvalAt(x);
    481 
    482         if (TMath::Abs(y-maxval/2)<fResolution)
    483             break;
    484 
    485         step /= sqrt2; // /2
    486         x += y>maxval/2 ? -step : +step;
    487     }
    488     */
    489 
    490     //
    491     // Now integrate the whole thing!
    492     //
    493     //   fTime      = cnt==31 ? -1 : x;
    494     //   fTimeDev   = fResolution;
    495     //   fSignal    = cnt==31 ? CalcIntegral(x) : CalcIntegral(maxpos);
    496     //   fSignalDev = 0;  // means: is valid
    497 }
     277}

  • trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

    r8545 r8546  
    7777    Double_t SearchY(Float_t maxpos, Float_t y) const;
    7878    Double_t SearchYup(Float_t maxpos, Float_t y) const;
    79 /*
    80     // Evaluate first solution for a possible maximum (x|first deriv==0)
    81     inline Double_t EvalDerivEq0S1(const Int_t i) const
    82     {
    83         // return the x value [0;1[ at which the derivative is zero (solution1)
    84 
    85         Double_t sumder = fDer2[i]+fDer2[i+1];
    86         Double_t difder = fDer2[i]-fDer2[i+1];
    87 
    88         Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
    89         Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
    90 
    91         Double_t x = 3*fDer2[i] - sqrt(3*sqt1 + 3*sqt2);
    92 
    93         Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
    94 
    95         return -x/denom;
    96     }
    97 
    98     // Evaluate second solution for a possible maximum (x|first deriv==0)
    99     inline Double_t EvalDerivEq0S2(const Int_t i) const
    100     {
    101         // return the x value [0;1[ at which the derivative is zero (solution2)
    102 
    103         Double_t sumder = fDer2[i]+fDer2[i+1];
    104         Double_t difder = fDer2[i]-fDer2[i+1];
    105 
    106         Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
    107         Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);
    108 
    109         Double_t x = 3*fDer2[i] + sqrt(3*sqt1 + 3*sqt2);
    110 
    111         Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
    112 
    113         return -x/denom;
    114     }
    115     */
    11679
    11780    Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
    … …  
    289252            xmax = xmax2;
    290253        }
    291   /*
    292         // Search real maximum in [i-0.5, i+1.5]
    293         Float_t xmax1, xmax2, xmax3;
    294         Float_t ymax1, ymax2, ymax3;
    295 
    296         Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1, 0.5, 1.0);
    297         Bool_t rc2 = GetMax(i,   xmax2, ymax2, 0.0, 1.0);
    298         Bool_t rc3 = i<fNum-1 && GetMax(i+1, xmax3, ymax3, 0.0, 0.5);
    299 
    300         // In case the medium bin is the first or last bin
    301         // take the lower or upper edge of the region into account.
    302         if (i==0)
    303         {
    304             xmax1 = 0;
    305             ymax1 = Eval(0, 0);
    306             rc1 = kTRUE;
    307         }
    308         if (i==fNum-1)
    309         {
    310             xmax3 = fNum-1e-5;
    311             ymax3 = Eval(fNum-1, 1);
    312             rc3 = kTRUE;
    313         }
    314 
    315         // Take a real default in case no maximum is found
    316         xmax=i+0.5;
    317         ymax=Eval(i, 0.5);
    318 
    319         //if (!rc1 && !rc2 && !rc3)
    320         //    cout << "!!!!!!!!!!!!!!!" << endl;
    321 
    322         if (rc1)
    323         {
    324             ymax = ymax1;
    325             xmax = xmax1;
    326         }
    327         else
    328             if (rc2)
    329             {
    330                 ymax = ymax2;
    331                 xmax = xmax2;
    332             }
    333             else
    334                 if (rc3)
    335                 {
    336                     ymax = ymax3;
    337                     xmax = xmax3;
    338                 }
    339 
    340         if (rc2 && ymax2>ymax)
    341         {
    342             ymax = ymax2;
    343             xmax = xmax2;
    344         }
    345         if (rc3 && ymax3>ymax)
    346         {
    347             ymax = ymax3;
    348             xmax = xmax3;
    349         }
    350         */
    351254    }
    352255
    … …  
    359262            return kFALSE;
    360263
    361         // const Float_t x1 = EvalDerivEq0S1(i);
    362         // const Float_t x2 = EvalDerivEq0S2(i);
    363 
    364264        const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(x1, i)<0;
    365265        const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(x2, i)<0;
    … …  
    384284        // Somehting must be wrong...
    385285        return kFALSE;
    386         /*
    387         std::cout << "?????????????" << std::endl;
    388 
    389         const Double_t y1 = Eval(i, x1);
    390         const Double_t y2 = Eval(i, x2);
    391 
    392         if (y1>y2)
    393         {
    394             xmax = i+x1;
    395             ymax = Eval(i, x1);
    396             return kTRUE;
    397         }
    398         else
    399         {
    400             xmax = i+x2;
    401             ymax = Eval(i, x2);
    402             return kTRUE;
    403         }
    404 
    405         return kFALSE;*/
    406     }
    407 /*
    408     inline Int_t GetMaxPos(Int_t i, Float_t &xmax, Float_t &ymax) const
    409     {
    410         Double_t x[3];
    411 
    412         x[0] = 0;
    413         // x[1] = 1; // This means we miss a possible maximum at the
    414                             // upper edge of the last interval...
    415 
    416         x[1] = EvalDerivEq0S1(i);
    417         x[2] = EvalDerivEq0S2(i);
    418 
    419         //y[0] = Eval(i, x[0]);
    420         //y[1] = Eval(i, x[1]);
    421         //y[1] = Eval(i, x[1]);
    422         //y[2] = Eval(i, x[2]);
    423 
    424         Int_t rc = 0;
    425         Double_t max = Eval(i, x[0]);
    426 
    427         for (Int_t j=1; j<3; j++)
    428         {
    429             if (x[j]<=0 || x[j]>=1)
    430                 continue;
    431 
    432             const Float_t y = Eval(i, x[j]);
    433             if (y>max)
    434             {
    435                 max = y;
    436                 rc = j;
    437             }
    438         }
    439 
    440         if (max>ymax)
    441         {
    442             xmax = x[rc]+i;
    443             ymax = max;
    444         }
    445 
    446         return rc;
    447     }
    448 
    449     inline void GetMaxPos(Int_t min, Int_t max, Float_t &xmax, Float_t &ymax) const
    450     {
    451         Float_t xmax=-1;
    452         Float_t ymax=-FLT_MAX;
    453 
    454         for (int i=min; i<max; i++)
    455             GetMaxPos(i, xmax, ymax);
    456 
    457         for (int i=min+1; i<max; i++)
    458         {
    459             Float_t y = Eval(i, 0);
    460             if (y>ymax)
    461             {
    462                 ymax = y;
    463                 xmax = i;
    464             }
    465         }
    466 
    467     }*/
    468 
     286    }
    469287
    470288    void InitDerivatives() const;
Note: See TracChangeset for help on using the changeset viewer.