Changeset 3148 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
02/14/04 11:48:43 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r3147 r3148  
    44
    55                                                 -*-*- END OF LINE -*-*-
     6
     7 2004/02/14: Thomas Bretz
     8
     9   * manalysis/MCerPhotEvt.[h,cc]:
     10     - added 'Iterator' facility, this will replace some for-loops
     11       in the near future
     12     
     13   * mbase/MTime.[h,cc]:
     14     - added a more powerfull interface to get and interprete the
     15       MTime contents as string
     16     - added a new constructor
     17
     18   * mreport/MReportTrigger.h:
     19     - fixed GetPixContent
     20
     21   * mtools/MCubicCoeff.cc, mtools/MCubicSpline.[h,cc]:
     22     - many small changes to simple details (like order of includes)
     23     - some speed improvements
     24     - many small simplifications
     25     - changed parts of the code to be more C++ like (eg Iterators
     26       instead of for-loops)
     27     - disentangles some if-cases
     28     - replaced some math.h function by TMath::
     29     - removed data-member fN (obsolete with iterators)
     30
    631
    732

  • trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.cc

    r2859 r3148  
    5454
    5555ClassImp(MCerPhotEvt);
     56ClassImp(MCerPhotEvtIter);
    5657
    5758using namespace std;
    … …  
    428429    *fLog << warn << "MCerPhotEvt::DrawPixelContent - not available." << endl;
    429430}
     431
     432TObject *MCerPhotEvtIter::Next()
     433{
     434    MCerPhotPix *pix;
     435    while ((pix = (MCerPhotPix*)TObjArrayIter::Next()))
     436        if (pix->IsPixelUsed())
     437            return pix;
     438    return pix;
     439}

  • trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h

    r2958 r3148  
    1717class MGeomCam;
    1818class MCerPhotPix;
     19class MCerPhotEvtIter;
    1920
    2021class MCerPhotEvt : public MParContainer, public MCamEvent
    2122{
     23    friend class MCerPhotEvtIter;
    2224private:
    2325    UInt_t        fNumPixels;
    … …  
    7880    void Clear(Option_t *opt=NULL) { Reset(); }
    7981
     82    // class MCamEvent
    8083    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
    8184    void DrawPixelContent(Int_t num) const;
     85
     86    operator TIterator*() const;
    8287
    8388    ClassDef(MCerPhotEvt, 2)    // class for an event containing cerenkov photons
    8489};
    8590
     91class MCerPhotEvtIter : public TObjArrayIter
     92{
     93private:
     94    Bool_t fUsedOnly;
     95public:
     96    MCerPhotEvtIter(const MCerPhotEvt *evt, Bool_t usedonly=kTRUE, Bool_t dir=kIterForward) : TObjArrayIter(evt->fPixels, dir), fUsedOnly(usedonly) { }
     97    TObject *Next();
     98    ClassDef(MCerPhotEvtIter, 0)
     99};
     100
     101inline MCerPhotEvt::operator TIterator*() const { return new MCerPhotEvtIter(this); }
     102
    86103#endif

  • trunk/MagicSoft/Mars/mbase/MTime.cc

    r3138 r3148  
    2929// A generalized MARS time stamp.
    3030//
     31//
     32// We do not use floating point values here, because of several reasons:
     33//  - having the times stored in integers only is more accurate and
     34//    more reliable in comparison conditions
     35//  - storing only integers gives similar bit-pattern for similar times
     36//    which makes compression (eg gzip algorithm in TFile) more
     37//    successfull
     38//
     39// Note, that there are many conversion function converting the day time
     40// into a readable string. Also a direct interface to SQL time strings
     41// is available.
     42//
     43// If you are using MTime containers as axis lables in root histograms
     44// use GetAxisTime(). Make sure that you use the correct TimeFormat
     45// on your TAxis (see GetAxisTime())
     46//
     47//
    3148// WARNING: Be carefull changing this class. It is also used in the
    3249//          MAGIC drive software cosy!
    … …  
    3451// Remarke: If you encounter strange behaviour, check the casting.
    3552//          Note, that on Linux machines ULong_t and UInt_t is the same.
     53//
    3654//
    3755// Version 1:
    … …  
    5270
    5371#include <iomanip>
     72
     73#include <time.h>     // struct tm
    5474#include <sys/time.h> // struct timeval
    5575
    … …  
    6585const UInt_t MTime::kHour = 3600000;         // [ms] one hour
    6686const UInt_t MTime::kDay  = MTime::kHour*24; // [ms] one day
     87
     88// --------------------------------------------------------------------------
     89//
     90// Constructor. Calls SetMjd(d) for d>0 in all other cases the time
     91// is set to the current UTC time.
     92//
     93MTime::MTime(Double_t d)
     94{
     95    Init(0, 0);
     96    if (d<=0)
     97        Now();
     98    else
     99        SetMjd(d);
     100}
    67101
    68102// --------------------------------------------------------------------------
    … …  
    192226
    193227    Set(mjd, tm, ms*1000);
     228}
     229
     230// --------------------------------------------------------------------------
     231//
     232// Return contents as a TString of the form:
     233//   "dd.mm.yyyy hh:mm:ss.fff"
     234//
     235Bool_t MTime::SetString(const char *str)
     236{
     237    if (!str)
     238        return kFALSE;
     239
     240    UInt_t y, mon, d, h, m, s, ms;
     241    const Int_t n = sscanf(str, "%02u.%02u.%04u %02u:%02u:%02u.%03u",
     242                           &d, &mon, &y, &h, &m, &s, &ms);
     243
     244    return n==7 ? Set(y, mon, d, h, m, s, ms) : kFALSE;
     245}
     246
     247// --------------------------------------------------------------------------
     248//
     249// Return contents as a TString of the form:
     250//   "yyyy-mm-dd hh:mm:ss"
     251//
     252Bool_t MTime::SetSqlDateTime(const char *str)
     253{
     254    if (!str)
     255        return kFALSE;
     256
     257    UInt_t y, mon, d, h, m, s;
     258    const Int_t n = sscanf(str, "%04u-%02u-%02u %02u:%02u:%02u",
     259                           &y, &mon, &d, &h, &m, &s);
     260
     261    return n==6 ? Set(y, mon, d, h, m, s) : kFALSE;
     262}
     263
     264// --------------------------------------------------------------------------
     265//
     266// Return contents as a TString of the form:
     267//   "yyyymmddhhmmss"
     268//
     269Bool_t MTime::SetSqlTimeStamp(const char *str)
     270{
     271    if (!str)
     272        return kFALSE;
     273
     274    UInt_t y, mon, d, h, m, s;
     275    const Int_t n = sscanf(str, "%04u%02u%02u%02u%02u%02u",
     276                           &y, &mon, &d, &h, &m, &s);
     277
     278    return n==6 ? Set(y, mon, d, h, m, s) : kFALSE;
    194279}
    195280
    … …  
    255340//
    256341// Return contents as a TString of the form:
    257 //   "yyyy/mm/dd hh:mm:ss.fff"
     342//   "dd.mm.yyyy hh:mm:ss.fff"
    258343//
    259344TString MTime::GetString() const
    … …  
    265350    GetTime(h, m, s, ms);
    266351
    267     return TString(Form("%4d/%02d/%02d %02d:%02d:%02d.%03d", y, mon, d, h, m, s, ms));
     352    return TString(Form("%02d.%02d.%04d %02d:%02d:%02d.%03d", d, mon, y, h, m, s, ms));
     353}
     354
     355// --------------------------------------------------------------------------
     356//
     357// Return contents as a string format'd with strftime:
     358// Here is a short summary of the most important formats. For more
     359// information see the man page (or any other description) of
     360// strftime...
     361//
     362//  %a  The abbreviated weekday name according to the current locale.
     363//  %A  The full weekday name according to the current locale.
     364//  %b  The abbreviated month name according to the current locale.
     365//  %B  The full month name according to the current locale.
     366//  %c  The preferred date and time representation for the current locale.
     367//  %d  The day of the month as a decimal number (range  01 to 31).
     368//  %e  Like %d, the day of the month as a decimal number,
     369//      but a leading zero is replaced by a space.
     370//  %H  The hour as a decimal number using a 24-hour clock (range 00 to 23)
     371//  %k  The hour (24-hour clock) as a decimal number (range 0 to 23);
     372//      single digits are preceded by a blank.
     373//  %m  The month as a decimal number (range 01 to 12).
     374//  %M  The minute as a decimal number (range 00 to 59).
     375//  %R  The time in 24-hour notation (%H:%M).  For a
     376//      version including the seconds, see %T below.
     377//  %S  The second as a decimal number (range 00 to 61).
     378//  %T  The time in 24-hour notation (%H:%M:%S).
     379//  %x  The preferred date representation for the current
     380//      locale without the time.
     381//  %X  The preferred time representation for the current
     382//      locale without the date.
     383//  %y  The year as a decimal number without a century (range 00 to 99).
     384//  %Y  The year as a decimal number including the century.
     385//  %+  The date and time in date(1) format.
     386//
     387// The default is: Tuesday 16.February 2004 12:17:22
     388//
     389// The maximum size of the return string is 128 (incl. NULL)
     390//
     391TString MTime::GetStringFmt(const char *fmt) const
     392{
     393    if (!fmt)
     394        fmt = "%A %e.%B %Y %H:%M:%S";
     395
     396    UShort_t y, ms;
     397    Byte_t mon, d, h, m, s;
     398
     399    GetDate(y, mon, d);
     400    GetTime(h, m, s, ms);
     401
     402    struct tm time;
     403    time.tm_sec   = s;
     404    time.tm_min   = m;
     405    time.tm_hour  = h;
     406    time.tm_mday  = d;
     407    time.tm_mon   = mon-1;
     408    time.tm_year  = y-1900;
     409    time.tm_isdst = 0;
     410
     411    // recalculate tm_yday and tm_wday
     412    mktime(&time);
     413
     414    char ret[128];
     415    return TString(strftime(ret, 127, fmt, &time) ? ret : "");
     416}
     417
     418// --------------------------------------------------------------------------
     419//
     420// Return contents as a TString of the form:
     421//   "yyyy-mm-dd hh:mm:ss"
     422//
     423TString MTime::GetSqlDateTime() const
     424{
     425    return GetStringFmt("%Y-%m-%d %H:%M:%S");
     426}
     427
     428// --------------------------------------------------------------------------
     429//
     430// Return contents as a TString of the form:
     431//   "yyyymmddhhmmss"
     432//
     433TString MTime::GetSqlTimeStamp() const
     434{
     435    return GetStringFmt("%Y%m%d%H%M%S");
    268436}
    269437
    … …  
    275443TString MTime::GetFileName() const
    276444{
    277     UShort_t y;
    278     Byte_t mon, d, h, m, s;
    279 
    280     GetDate(y, mon, d);
    281     GetTime(h, m, s);
    282 
    283     return TString(Form("%04d%02d%02d_%02d%02d%02d", y, mon, d, h, m, s));
     445    return GetStringFmt("%Y%m%d_%H%M%S");
    284446}
    285447

  • trunk/MagicSoft/Mars/mbase/MTime.h

    r3138 r3148  
    5656        Set(tm);
    5757    }
     58    MTime(Double_t d);
    5859    MTime(const MTime& t) : fMjd(t.fMjd), fTime(t.fTime), fNanoSec(t.fNanoSec)
    5960    {
    … …  
    7778    Bool_t   Set(UShort_t y, Byte_t m, Byte_t d, Byte_t h=13, Byte_t min=0, Byte_t s=0, UShort_t ms=0, UInt_t ns=0);
    7879    void     Set(const struct timeval &tv);
     80    Bool_t   SetString(const char *str);
     81    Bool_t   SetSqlDateTime(const char *str);
     82    Bool_t   SetSqlTimeStamp(const char *str);
    7983    void     SetCT1Time(UInt_t mjd, UInt_t t1, UInt_t t0);
    8084    Bool_t   UpdMagicTime(Byte_t h, Byte_t m, Byte_t s, UInt_t ns);
    … …  
    8286    Double_t GetMjd() const;
    8387    TString  GetString() const;
     88    TString  GetStringFmt(const char *fmt=0) const;
     89    TString  GetSqlDateTime() const;
     90    TString  GetSqlTimeStamp() const;
    8491    TString  GetFileName() const;
    8592    void     GetDate(UShort_t &y, Byte_t &m, Byte_t &d) const;

  • trunk/MagicSoft/Mars/mreport/MReportTrigger.h

    r3082 r3148  
    2020
    2121    TArrayF fPrescalerRates;    //[Hz] L2 prescaler rates
    22     //TArrayF fRates;           //[Hz] curently undefined
     22    //TArrayF fRates;           //[Hz] currently undefined
    2323
    2424    Int_t InterpreteBody(TString &str);
    … …  
    3131
    3232    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
    33       {
    34         if(idx<19)
    35           {
    36             val = fPrescalerRates[idx];
    37             return val>0;
    38           }
    39         else
    40           val = kFALSE;
    41       }
     33    {
     34        if(idx>18)
     35            return kFALSE;
     36
     37        val = fPrescalerRates[idx];
     38        return kTRUE;
     39    }
     40
    4241    void DrawPixelContent(Int_t num) const
    43       {
    44       }
     42    {
     43    }
    4544
    4645    ClassDef(MReportTrigger, 1) // Class for TRIGGER-REPORT information

  • trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc

    r3022 r3148  
    2424
    2525//////////////////////////////////////////////////////////////////////////////
    26 //                                                                          //
    27 //  Cubic Spline Interpolation                                              //
    28 //                                                                          //
     26//
     27//  Cubic Spline Interpolation
     28//
    2929//////////////////////////////////////////////////////////////////////////////
    30 
    3130#include "MCubicCoeff.h"
    3231
    33 #include "TMath.h"
     32#include <TMath.h>
    3433
    3534#include "MLog.h"
    … …  
    4544// where x is the independent variable, 2 and 3 are exponents
    4645//
    47 
    4846MCubicCoeff::MCubicCoeff(Double_t x, Double_t xNext, Double_t y, Double_t yNext,
    4947                         Double_t a, Double_t b, Double_t c) :
    … …  
    5149{
    5250    fH = fXNext - fX;
    53     if(!EvalMinMax())
    54     {
    55         gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
    56         fMin = 0;
    57         fMax = 0;
    58     }
     51    if (EvalMinMax())
     52        return;
     53
     54    gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
     55    fMin = 0;
     56    fMax = 0;
    5957}
    6058
    … …  
    6664Double_t MCubicCoeff::Eval(Double_t x)
    6765{
    68     Double_t dx = x - fX;
    69     return (fY+dx*(fC+dx*(fB+dx*fA)));
     66    const Double_t dx = x - fX;
     67    return fY + dx*(fC + dx*(fB + dx*fA));
    7068}
    7169
    … …  
    7977Bool_t MCubicCoeff::EvalMinMax()
    8078{
    81     fMin = fMax = fY;
    82     fAbMin = fAbMax = fX;
     79    fMin = fY;
     80    fMax = fY;
     81
     82    fAbMin = fX;
     83    fAbMax = fX;
     84
    8385    if (fYNext < fMin)
    8486    {
    85         fMin = fYNext;
     87        fMin   = fYNext;
    8688        fAbMin = fXNext;
    8789    }
    8890    if (fYNext > fMax)
    8991    {
    90         fMax = fYNext;
     92        fMax   = fYNext;
    9193        fAbMax = fXNext;
    9294    }
    93     Double_t tempMinMax;
    94     Double_t delta = 4.0*fB*fB - 12.0*fA*fC;
    95         if (delta >= 0.0 && fA != 0.0)
    96         {
    97             Double_t sqrtDelta = sqrt(delta);
    98             Double_t xPlus  = (-2.0*fB + sqrtDelta)/(6.0*fA);
    99             Double_t xMinus = (-2.0*fB - sqrtDelta)/(6.0*fA);
    100             if (xPlus >= 0.0 && xPlus <= fH)
    101             {
    102                 tempMinMax = this->Eval(fX+xPlus);
    103                 if (tempMinMax < fMin)
    104                 {
    105                     fMin = tempMinMax;
    106                     fAbMin = fX + xPlus;
    107                 }
    108                 if (tempMinMax > fMax)
    109                 {
    110                     fMax = tempMinMax;
    111                     fAbMax = fX + xPlus;
    112                 }
    113             }
    114             if (xMinus >= 0.0 && xMinus <= fH)
    115             {
    116                 tempMinMax = this->Eval(fX+xMinus);
    117                 if (tempMinMax < fMin)
    118                 {
    119                     fMin = tempMinMax;
    120                     fAbMin = fX + xMinus;
    121                 }
    122                 if (tempMinMax > fMax)
    123                 {
    124                     fMax = tempMinMax;
    125                     fAbMax = fX + xMinus;
    126                 }
    127             }
    128         }
    129 /* if fA is zero then we have only one possible solution */
    130         else if (fA == 0.0 && fB != 0.0)
    131         {
    132             Double_t xSolo = - (fC/(2.0*fB));
    133             if (xSolo >= 0.0 && xSolo <= fH)
    134             {
    135                 tempMinMax = this->Eval(fX+xSolo);
    136                 if (tempMinMax < fMin)
    137                 {
    138                     fMin = tempMinMax;
    139                     fAbMin = fX + xSolo;
    140                 }
    141                 if (tempMinMax > fMax)
    142                 {
    143                     fMax = tempMinMax;
    144                     fAbMax = fX + xSolo;
    145                 }
    146             }
    147         }
     95
     96    const Double_t delta = fB*fB*4 - fA*fC*12;
     97    if (delta >= 0 && fA != 0)
     98    {
     99        const Double_t sqrtDelta = TMath::Sqrt(delta);
     100
     101        const Double_t xPlus  = (-fB*2 + sqrtDelta)/(fA*6);
     102        const Double_t xMinus = (-fB*2 - sqrtDelta)/(fA*6);
     103
     104        if (xPlus >= 0 && xPlus <= fH)
     105        {
     106            const Double_t tempMinMax = Eval(fX+xPlus);
     107            if (tempMinMax < fMin)
     108            {
     109                fMin = tempMinMax;
     110                fAbMin = fX + xPlus;
     111            }
     112            if (tempMinMax > fMax)
     113            {
     114                fMax = tempMinMax;
     115                fAbMax = fX + xPlus;
     116            }
     117        }
     118        if (xMinus >= 0 && xMinus <= fH)
     119        {
     120            const Double_t tempMinMax = Eval(fX+xMinus);
     121            if (tempMinMax < fMin)
     122            {
     123                fMin = tempMinMax;
     124                fAbMin = fX + xMinus;
     125            }
     126            if (tempMinMax > fMax)
     127            {
     128                fMax = tempMinMax;
     129                fAbMax = fX + xMinus;
     130            }
     131        }
     132        return kTRUE;
     133    }
     134
     135    /* if fA is zero then we have only one possible solution */
     136    if (fA == 0 && fB != 0)
     137    {
     138        const Double_t xSolo = -fC/(fB*2);
     139
     140        if (xSolo < 0 || xSolo > fH)
     141            return kTRUE;
     142
     143        const Double_t tempMinMax = Eval(fX+xSolo);
     144        if (tempMinMax < fMin)
     145        {
     146            fMin = tempMinMax;
     147            fAbMin = fX + xSolo;
     148        }
     149        if (tempMinMax > fMax)
     150        {
     151            fMax = tempMinMax;
     152            fAbMax = fX + xSolo;
     153        }
     154        return kTRUE;
     155    }
     156
    148157    return kTRUE;
    149158}
    … …  
    157166// There could be three or one real solutions
    158167//
    159 
    160168Short_t MCubicCoeff::FindCardanRoot(Double_t y, Double_t *x)
    161169{
    162 
    163     Short_t whichRoot = -1;
    164 
    165     Double_t a = fB/fA;
    166     Double_t b = fC/fA;
    167     Double_t c = (fY - y)/fA;
    168 
    169     Double_t q = (a*a-3.0*b)/9.0;
    170     Double_t r = (2.0*a*a*a-9.0*a*b+27.0*c)/54.0;
    171 
    172     Double_t aOver3 = a/3.0;
    173     Double_t r2 = r*r;
    174     Double_t q3 = q*q*q;
     170    const Double_t a = fB/fA;
     171    const Double_t b = fC/fA;
     172    const Double_t c = (fY - y)/fA;
     173
     174    const Double_t q = (a*a - b*3)/9;
     175    const Double_t r = (a*a*a*2 - a*b*9 + c*27)/54;
     176
     177    const Double_t aOver3 = a/3;
     178    const Double_t r2 = r*r;
     179    const Double_t q3 = q*q*q;
    175180   
    176181    if (r2 < q3) //3 real sol
    177182    {
    178         Double_t sqrtQ = TMath::Sqrt(q);
    179         Double_t min2SqQ = -2.0*sqrtQ;
    180         Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));
    181 
    182         x[0] = min2SqQ * TMath::Cos(theta/3.0) - aOver3;
    183         x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3.0) - aOver3;
    184         x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3.0) - aOver3;
     183        const Double_t sqrtQ = TMath::Sqrt(q);
     184        const Double_t min2SqQ = -sqrtQ*2;
     185        const Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));
     186
     187        x[0] = min2SqQ * TMath::Cos(theta/3) - aOver3;
     188        x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3) - aOver3;
     189        x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3) - aOver3;
     190
    185191        for (Int_t i = 0; i < 3; i++)
    186             if (x[i] >= 0.0 && x[i] <= fH)
     192            if (x[i] >= 0 && x[i] <= fH)
    187193            {
    188                 x[i] = x[i] + fX;
    189                 whichRoot = i;
    190                 break;
     194                x[i] += fX;
     195                return i;
    191196            }
    192     }
    193     else //only 1 real sol
    194     {
    195         Double_t s;
    196         if (r == 0.0)
    197             s = 0.0;
    198         else if (r < 0.0)
    199             s = TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
    200         else // r > 0.0
    201             s = - TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
    202         if (s == 0.0)
    203             x[0] = - aOver3;
    204         else
    205             x[0] = (s + q/s) - aOver3;
    206         if (x[0] >= 0.0 && x[0] <= fH)
    207         {
    208             x[0] = x[0] + fX;
    209             whichRoot = 0;
    210         }
    211     }
    212     return whichRoot;
     197        return -1;
     198    }
     199
     200    const Double_t s = r==0 ? 0 : -TMath::Sign(TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3), 1./3), r);
     201
     202    x[0] = s==0 ? - aOver3 : (s + q/s) - aOver3;
     203
     204    if (x[0] < 0 || x[0] > fH)
     205        return -1;
     206
     207    x[0] += fX;
     208    return 0;
    213209}
    214210
    … …  
    220216Bool_t MCubicCoeff :: IsIn(Double_t x)
    221217{
    222     if (x >= fX && x <= fXNext)
    223         return kTRUE;
    224     else
    225         return kFALSE;
    226 }
     218    return x >= fX && x <= fXNext;
     219}

  • trunk/MagicSoft/Mars/mtools/MCubicSpline.cc

    r3022 r3148  
    2424
    2525//////////////////////////////////////////////////////////////////////////////
    26 //                                                                          //
    27 //  Cubic Spline Interpolation                                              //
    28 //                                                                          //
     26//
     27//  Cubic Spline Interpolation
     28//
    2929//////////////////////////////////////////////////////////////////////////////
    30 
    3130#include "MCubicSpline.h"
    3231
    33 #include "MCubicCoeff.h"
    34 
    35 #include "TMath.h"
     32#include <TMath.h>
    3633
    3734#include "MLog.h"
    3835#include "MLogManip.h"
    3936
     37#include "MCubicCoeff.h"
     38
    4039ClassImp(MCubicSpline);
    4140
    … …  
    4847//
    4948MCubicSpline::MCubicSpline(Byte_t *y, Byte_t *x, Bool_t areAllEq,
    50                            Int_t n, Double_t begSD, Double_t endSD):
    51     fN(n)
     49                           Int_t n, Double_t begSD, Double_t endSD)
    5250{
    5351    Init(y,x,areAllEq,n,begSD,endSD);
    … …  
    6260{
    6361    Byte_t x[]={0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,0x0C,0x0D,0x0E};
    64     fN = 15;
    6562    Init(y,x,kTRUE,15,0.0,0.0);
    6663}
    … …  
    7572
    7673{
    77     Double_t *temp = new Double_t[fN-1];
    78     Double_t *ysd = new Double_t[fN-1];
    79     Double_t p,h;
    80     MCubicCoeff *tempCoeff;
    81     fCoeff = new TObjArray(fN-1,0);
     74    Double_t *temp = new Double_t[n-1];
     75    Double_t *ysd  = new Double_t[n-1];
     76
     77    fCoeff = new TObjArray(n-1,0);
     78
     79    ysd[0]  =begSD;
     80    temp[0] =begSD;
     81    ysd[n-1]=endSD;
    8282   
     83    Double_t h = x[1]-x[0];
     84
    8385    if (areAllEq)
    8486    {
    85         h = x[1]-x[0];
    86         ysd[0]=temp[0]=begSD;
    87         ysd[n-1]=endSD;
    8887        for(Int_t i = 1; i < n-1; i++)
    8988        {
    90             p = 0.5*ysd[i-1]+2.0;
    91             ysd[i] = (-0.5)/p;
    92             temp[i] = (y[i+1]-2*y[i]+y[i-1])/h;
    93             temp[i] = (6.0*temp[i]/h-0.5*temp[i-1])/p;
    94         }
    95         for(Int_t i = n-2; i > 0; i--)
    96             ysd[i]=ysd[i]*ysd[i+1]+temp[i];
    97         for(Int_t i = 0; i < n-1; i++)
    98         {
    99             tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h),
    100                                         ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6);
    101             fCoeff->AddAt(tempCoeff,i);
    102         }
    103         delete [] temp;
    104         delete [] ysd;
     89            const Double_t p = ysd[i-1]/2+2;
     90
     91            ysd[i]  = -0.5/p;
     92            temp[i] = (y[i+1] - y[i]*2 + y[i-1])/h;
     93            temp[i] = (temp[i]*6/h-temp[i-1]/2)/p;
     94        }
    10595    }
    10696    else
    10797    {
    108         Double_t sig;
    109         ysd[0]=temp[0]=begSD;
    110         ysd[n-1]=endSD;
    11198        for(Int_t i = 1; i < n-1; i++)
    11299        {
    113             sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
    114             p = sig*ysd[i-1]+2.0;
    115             ysd[i] = (sig-1.0)/p;
     100            const Double_t sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
     101
     102            const Double_t p = sig*ysd[i-1]+2;
     103
     104            ysd[i]  = (sig-1.0)/p;
    116105            temp[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
    117             temp[i] = (6.0*temp[i]/(x[i+1]-x[i-1])-sig*temp[i-1])/p;
     106            temp[i] = (temp[i]*6/(x[i+1]-x[i-1])-sig*temp[i-1])/p;
    118107        }
    119         for(Int_t i = n-2; i > 0; i--)
    120             ysd[i]=ysd[i]*ysd[i+1]+temp[i];
    121         for(Int_t i = 0; i < n-1; i++)
    122         {
    123             h = x[i+1]-x[i];
    124             tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h),
    125                                         ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6);
    126             fCoeff->AddAt(tempCoeff,i);
    127         }
    128         delete [] temp;
    129         delete [] ysd;
    130     }
     108    }
     109
     110    for(Int_t i = n-2; i > 0; i--)
     111        ysd[i] = ysd[i]*ysd[i+1] + temp[i];
     112
     113    for(Int_t i = 0; i < n-1; i++)
     114    {
     115        if (!areAllEq)
     116            h = x[i+1]-x[i];
     117
     118        MCubicCoeff *c = new MCubicCoeff(x[i], x[i+1], y[i], y[i+1], (ysd[i+1]-ysd[i])/(h*6),
     119                                         ysd[i]/2, (y[i+1]-y[i])/h-(h*(ysd[i+1]+ysd[i]*2))/6);
     120        fCoeff->AddAt(c, i);
     121    }
     122
     123    delete [] temp;
     124    delete [] ysd;
    131125}
    132126
    … …  
    134128{
    135129    fCoeff->Delete();
     130    delete fCoeff;
    136131}
    137132
    … …  
    142137Double_t MCubicSpline :: Eval(Double_t x)
    143138{
    144     for (Int_t i = 0; i < fN-1; i++)
    145         if (((MCubicCoeff*)fCoeff->At(i))->IsIn(x))
    146             return ((MCubicCoeff*)fCoeff->At(i))->Eval(x);
     139    const Int_t n = fCoeff->GetSize()-1;
     140    for (Int_t i = 0; i < n; i++)
     141    {
     142        MCubicCoeff *c = (MCubicCoeff*)fCoeff->UncheckedAt(i);
     143        if (c->IsIn(x))
     144            return c->Eval(x);
     145    }
     146
    147147    gLog << warn << "Cannot evaluate Spline at " << x << "; returning 0";
    148     return 0.0;
     148
     149    return 0;
    149150}
    150151
    … …  
    155156Double_t MCubicSpline :: EvalMax()
    156157{
    157     Double_t temp_max;
    158     Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax();
    159     for (Int_t i = 1; i < fN-1; i++)
    160     {
    161         temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax();
    162         if (temp_max > max)
    163             max = temp_max;
    164     }
     158    Double_t max = -FLT_MAX;
     159
     160    TIter Next(fCoeff);
     161    MCubicCoeff *c;
     162    while ((c=(MCubicCoeff*)Next()))
     163        max = TMath::Max(max, c->GetMax());
     164
    165165    return max;
    166166}
    … …  
    172172Double_t MCubicSpline :: EvalMin()
    173173{
    174     Double_t temp_min;
    175     Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin();
    176     for (Int_t i = 1; i < fN-1; i++)
    177     {
    178         temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin();
    179         if (temp_min < min)
    180             min = temp_min;
    181     }
     174    Double_t min = FLT_MAX;
     175
     176    TIter Next(fCoeff);
     177    MCubicCoeff *c;
     178    while ((c=(MCubicCoeff*)Next()))
     179        min = TMath::Min(min, c->GetMax());
     180
    182181    return min;
    183182}
    … …  
    189188Double_t MCubicSpline :: EvalAbMax()
    190189{
    191     Double_t temp_max;
    192     Double_t abMax = ((MCubicCoeff*)fCoeff->At(0))->GetAbMax();
    193     Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax();
    194     for (Int_t i = 1; i < fN-1; i++)
    195     {
    196         temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax();
    197         if (temp_max > max)
    198         {
    199             max = temp_max;
    200             abMax = ((MCubicCoeff*)fCoeff->At(i))->GetAbMax();
    201         }
    202     }
    203     return abMax;
     190    Double_t max = -FLT_MAX;
     191
     192    TIter Next(fCoeff);
     193
     194    MCubicCoeff *c;
     195    MCubicCoeff *cmax=0;
     196
     197    while ((c=(MCubicCoeff*)Next()))
     198    {
     199        const Double_t temp = c->GetMax();
     200        if (temp <= max)
     201            continue;
     202
     203        max = temp;
     204        cmax = c;
     205    }
     206
     207    return cmax ? cmax->GetAbMax() : -FLT_MAX;
    204208}
    205209
    … …  
    210214Double_t MCubicSpline :: EvalAbMin()
    211215{
    212     Double_t temp_min;
    213     Double_t abMin = ((MCubicCoeff*)fCoeff->At(0))->GetAbMin();
    214     Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin();
    215     for (Int_t i = 1; i < fN-1; i++)
    216     {
    217         temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin();
    218         if (temp_min < min)
    219         {
    220             min = temp_min;
    221             abMin = ((MCubicCoeff*)fCoeff->At(i))->GetAbMin();
    222         }
    223     }
    224     return abMin;
     216    Double_t min = FLT_MAX;
     217
     218    TIter Next(fCoeff);
     219
     220    MCubicCoeff *c;
     221    MCubicCoeff *cmin=0;
     222
     223    while ((c=(MCubicCoeff*)Next()))
     224    {
     225        const Double_t temp = c->GetMax();
     226        if (temp >= min)
     227            continue;
     228
     229        min = temp;
     230        cmin = c;
     231    }
     232
     233    return cmin ? cmin->GetAbMax() : FLT_MAX;
    225234}
    226235
    … …  
    233242Double_t MCubicSpline :: FindVal(Double_t y, Double_t x0, Char_t direction = 'l')
    234243{
    235     Short_t whichRoot;
    236     Double_t tempRoot;
    237     Double_t *roots = new Double_t[3];
    238 
    239     for (Int_t i = 0; i < fN-1; i++)
    240     {
    241         if(((MCubicCoeff*)fCoeff->At(i))->IsIn(x0))
    242         {
    243             if(direction == 'l')
    244             {
    245                 for (Int_t j = i; j >= 0; j--)
    246                 {
    247                     whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
    248                     if (whichRoot >= 0 )
    249                     {
    250                         tempRoot = roots[whichRoot];
    251                         delete [] roots;
    252                         return tempRoot;
    253                     }
    254                 }
    255             }
    256             if(direction == 'r')
    257             {
    258                 for (Int_t j = i; j < fN-1; j++)
    259                 {
    260                     whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
    261                     if (whichRoot >= 0)
    262                     {
    263                         tempRoot = roots[whichRoot];
    264                         delete [] roots;
    265                         return tempRoot;
    266                     }
    267                 }
    268             }
    269         }
    270     }
     244    Double_t roots[3] = { 0, 0, 0 };
     245
     246    const Int_t n = fCoeff->GetSize()-1;
     247
     248    for (Int_t i = 0; i < n; i++)
     249    {
     250        if (!((MCubicCoeff*)fCoeff->At(i))->IsIn(x0))
     251            continue;
     252
     253        switch (direction)
     254        {
     255        case 'l':
     256            for (Int_t j = i; j >= 0; j--)
     257            {
     258                const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
     259                if (whichRoot >= 0 )
     260                    return roots[whichRoot];
     261            }
     262            break;
     263
     264        case 'r':
     265            for (Int_t j = i; j < n; j++)
     266            {
     267                const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
     268                if (whichRoot >= 0)
     269                    return roots[whichRoot];
     270            }
     271            break;
     272        }
     273    }
     274
    271275    gLog << warn << "Nothing found calling MCubicSpline :: FindVal(), returning 0" << endl;
    272     return 0.0;
    273 }
     276
     277    return 0;
     278}

  • trunk/MagicSoft/Mars/mtools/MCubicSpline.h

    r3022 r3148  
    1616 private:
    1717    TObjArray *fCoeff; //array of the coefficients
    18     Int_t fN; //number of points
     18
    1919    void Init(Byte_t *y, Byte_t *x, Bool_t areAllEq, Int_t n, Double_t begSD, Double_t endSD);
    2020
Note: See TracChangeset for help on using the changeset viewer.