Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3147) +++ /trunk/MagicSoft/Mars/Changelog (revision 3148) @@ -4,4 +4,29 @@ -*-*- END OF LINE -*-*- + + 2004/02/14: Thomas Bretz + + * manalysis/MCerPhotEvt.[h,cc]: + - added 'Iterator' facility, this will replace some for-loops + in the near future + + * mbase/MTime.[h,cc]: + - added a more powerfull interface to get and interprete the + MTime contents as string + - added a new constructor + + * mreport/MReportTrigger.h: + - fixed GetPixContent + + * mtools/MCubicCoeff.cc, mtools/MCubicSpline.[h,cc]: + - many small changes to simple details (like order of includes) + - some speed improvements + - many small simplifications + - changed parts of the code to be more C++ like (eg Iterators + instead of for-loops) + - disentangles some if-cases + - replaced some math.h function by TMath:: + - removed data-member fN (obsolete with iterators) + Index: /trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.cc (revision 3147) +++ /trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.cc (revision 3148) @@ -54,4 +54,5 @@ ClassImp(MCerPhotEvt); +ClassImp(MCerPhotEvtIter); using namespace std; @@ -428,2 +429,11 @@ *fLog << warn << "MCerPhotEvt::DrawPixelContent - not available." << endl; } + +TObject *MCerPhotEvtIter::Next() +{ + MCerPhotPix *pix; + while ((pix = (MCerPhotPix*)TObjArrayIter::Next())) + if (pix->IsPixelUsed()) + return pix; + return pix; +} Index: /trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h (revision 3147) +++ /trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h (revision 3148) @@ -17,7 +17,9 @@ class MGeomCam; class MCerPhotPix; +class MCerPhotEvtIter; class MCerPhotEvt : public MParContainer, public MCamEvent { + friend class MCerPhotEvtIter; private: UInt_t fNumPixels; @@ -78,9 +80,24 @@ void Clear(Option_t *opt=NULL) { Reset(); } + // class MCamEvent Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; void DrawPixelContent(Int_t num) const; + + operator TIterator*() const; ClassDef(MCerPhotEvt, 2) // class for an event containing cerenkov photons }; +class MCerPhotEvtIter : public TObjArrayIter +{ +private: + Bool_t fUsedOnly; +public: + MCerPhotEvtIter(const MCerPhotEvt *evt, Bool_t usedonly=kTRUE, Bool_t dir=kIterForward) : TObjArrayIter(evt->fPixels, dir), fUsedOnly(usedonly) { } + TObject *Next(); + ClassDef(MCerPhotEvtIter, 0) +}; + +inline MCerPhotEvt::operator TIterator*() const { return new MCerPhotEvtIter(this); } + #endif Index: /trunk/MagicSoft/Mars/mbase/MTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTime.cc (revision 3147) +++ /trunk/MagicSoft/Mars/mbase/MTime.cc (revision 3148) @@ -29,4 +29,21 @@ // A generalized MARS time stamp. // +// +// We do not use floating point values here, because of several reasons: +// - having the times stored in integers only is more accurate and +// more reliable in comparison conditions +// - storing only integers gives similar bit-pattern for similar times +// which makes compression (eg gzip algorithm in TFile) more +// successfull +// +// Note, that there are many conversion function converting the day time +// into a readable string. Also a direct interface to SQL time strings +// is available. +// +// If you are using MTime containers as axis lables in root histograms +// use GetAxisTime(). Make sure that you use the correct TimeFormat +// on your TAxis (see GetAxisTime()) +// +// // WARNING: Be carefull changing this class. It is also used in the // MAGIC drive software cosy! @@ -34,4 +51,5 @@ // Remarke: If you encounter strange behaviour, check the casting. // Note, that on Linux machines ULong_t and UInt_t is the same. +// // // Version 1: @@ -52,4 +70,6 @@ #include + +#include // struct tm #include // struct timeval @@ -65,4 +85,18 @@ const UInt_t MTime::kHour = 3600000; // [ms] one hour const UInt_t MTime::kDay = MTime::kHour*24; // [ms] one day + +// -------------------------------------------------------------------------- +// +// Constructor. Calls SetMjd(d) for d>0 in all other cases the time +// is set to the current UTC time. +// +MTime::MTime(Double_t d) +{ + Init(0, 0); + if (d<=0) + Now(); + else + SetMjd(d); +} // -------------------------------------------------------------------------- @@ -192,4 +226,55 @@ Set(mjd, tm, ms*1000); +} + +// -------------------------------------------------------------------------- +// +// Return contents as a TString of the form: +// "dd.mm.yyyy hh:mm:ss.fff" +// +Bool_t MTime::SetString(const char *str) +{ + if (!str) + return kFALSE; + + UInt_t y, mon, d, h, m, s, ms; + const Int_t n = sscanf(str, "%02u.%02u.%04u %02u:%02u:%02u.%03u", + &d, &mon, &y, &h, &m, &s, &ms); + + return n==7 ? Set(y, mon, d, h, m, s, ms) : kFALSE; +} + +// -------------------------------------------------------------------------- +// +// Return contents as a TString of the form: +// "yyyy-mm-dd hh:mm:ss" +// +Bool_t MTime::SetSqlDateTime(const char *str) +{ + if (!str) + return kFALSE; + + UInt_t y, mon, d, h, m, s; + const Int_t n = sscanf(str, "%04u-%02u-%02u %02u:%02u:%02u", + &y, &mon, &d, &h, &m, &s); + + return n==6 ? Set(y, mon, d, h, m, s) : kFALSE; +} + +// -------------------------------------------------------------------------- +// +// Return contents as a TString of the form: +// "yyyymmddhhmmss" +// +Bool_t MTime::SetSqlTimeStamp(const char *str) +{ + if (!str) + return kFALSE; + + UInt_t y, mon, d, h, m, s; + const Int_t n = sscanf(str, "%04u%02u%02u%02u%02u%02u", + &y, &mon, &d, &h, &m, &s); + + return n==6 ? Set(y, mon, d, h, m, s) : kFALSE; } @@ -255,5 +340,5 @@ // // Return contents as a TString of the form: -// "yyyy/mm/dd hh:mm:ss.fff" +// "dd.mm.yyyy hh:mm:ss.fff" // TString MTime::GetString() const @@ -265,5 +350,88 @@ GetTime(h, m, s, ms); - return TString(Form("%4d/%02d/%02d %02d:%02d:%02d.%03d", y, mon, d, h, m, s, ms)); + return TString(Form("%02d.%02d.%04d %02d:%02d:%02d.%03d", d, mon, y, h, m, s, ms)); +} + +// -------------------------------------------------------------------------- +// +// Return contents as a string format'd with strftime: +// Here is a short summary of the most important formats. For more +// information see the man page (or any other description) of +// strftime... +// +// %a The abbreviated weekday name according to the current locale. +// %A The full weekday name according to the current locale. +// %b The abbreviated month name according to the current locale. +// %B The full month name according to the current locale. +// %c The preferred date and time representation for the current locale. +// %d The day of the month as a decimal number (range 01 to 31). +// %e Like %d, the day of the month as a decimal number, +// but a leading zero is replaced by a space. +// %H The hour as a decimal number using a 24-hour clock (range 00 to 23) +// %k The hour (24-hour clock) as a decimal number (range 0 to 23); +// single digits are preceded by a blank. +// %m The month as a decimal number (range 01 to 12). +// %M The minute as a decimal number (range 00 to 59). +// %R The time in 24-hour notation (%H:%M). For a +// version including the seconds, see %T below. +// %S The second as a decimal number (range 00 to 61). +// %T The time in 24-hour notation (%H:%M:%S). +// %x The preferred date representation for the current +// locale without the time. +// %X The preferred time representation for the current +// locale without the date. +// %y The year as a decimal number without a century (range 00 to 99). +// %Y The year as a decimal number including the century. +// %+ The date and time in date(1) format. +// +// The default is: Tuesday 16.February 2004 12:17:22 +// +// The maximum size of the return string is 128 (incl. NULL) +// +TString MTime::GetStringFmt(const char *fmt) const +{ + if (!fmt) + fmt = "%A %e.%B %Y %H:%M:%S"; + + UShort_t y, ms; + Byte_t mon, d, h, m, s; + + GetDate(y, mon, d); + GetTime(h, m, s, ms); + + struct tm time; + time.tm_sec = s; + time.tm_min = m; + time.tm_hour = h; + time.tm_mday = d; + time.tm_mon = mon-1; + time.tm_year = y-1900; + time.tm_isdst = 0; + + // recalculate tm_yday and tm_wday + mktime(&time); + + char ret[128]; + return TString(strftime(ret, 127, fmt, &time) ? ret : ""); +} + +// -------------------------------------------------------------------------- +// +// Return contents as a TString of the form: +// "yyyy-mm-dd hh:mm:ss" +// +TString MTime::GetSqlDateTime() const +{ + return GetStringFmt("%Y-%m-%d %H:%M:%S"); +} + +// -------------------------------------------------------------------------- +// +// Return contents as a TString of the form: +// "yyyymmddhhmmss" +// +TString MTime::GetSqlTimeStamp() const +{ + return GetStringFmt("%Y%m%d%H%M%S"); } @@ -275,11 +443,5 @@ TString MTime::GetFileName() const { - UShort_t y; - Byte_t mon, d, h, m, s; - - GetDate(y, mon, d); - GetTime(h, m, s); - - return TString(Form("%04d%02d%02d_%02d%02d%02d", y, mon, d, h, m, s)); + return GetStringFmt("%Y%m%d_%H%M%S"); } Index: /trunk/MagicSoft/Mars/mbase/MTime.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTime.h (revision 3147) +++ /trunk/MagicSoft/Mars/mbase/MTime.h (revision 3148) @@ -56,4 +56,5 @@ Set(tm); } + MTime(Double_t d); MTime(const MTime& t) : fMjd(t.fMjd), fTime(t.fTime), fNanoSec(t.fNanoSec) { @@ -77,4 +78,7 @@ 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); void Set(const struct timeval &tv); + Bool_t SetString(const char *str); + Bool_t SetSqlDateTime(const char *str); + Bool_t SetSqlTimeStamp(const char *str); void SetCT1Time(UInt_t mjd, UInt_t t1, UInt_t t0); Bool_t UpdMagicTime(Byte_t h, Byte_t m, Byte_t s, UInt_t ns); @@ -82,4 +86,7 @@ Double_t GetMjd() const; TString GetString() const; + TString GetStringFmt(const char *fmt=0) const; + TString GetSqlDateTime() const; + TString GetSqlTimeStamp() const; TString GetFileName() const; void GetDate(UShort_t &y, Byte_t &m, Byte_t &d) const; Index: /trunk/MagicSoft/Mars/mreport/MReportTrigger.h =================================================================== --- /trunk/MagicSoft/Mars/mreport/MReportTrigger.h (revision 3147) +++ /trunk/MagicSoft/Mars/mreport/MReportTrigger.h (revision 3148) @@ -20,5 +20,5 @@ TArrayF fPrescalerRates; //[Hz] L2 prescaler rates - //TArrayF fRates; //[Hz] curently undefined + //TArrayF fRates; //[Hz] currently undefined Int_t InterpreteBody(TString &str); @@ -31,16 +31,15 @@ Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const - { - if(idx<19) - { - val = fPrescalerRates[idx]; - return val>0; - } - else - val = kFALSE; - } + { + if(idx>18) + return kFALSE; + + val = fPrescalerRates[idx]; + return kTRUE; + } + void DrawPixelContent(Int_t num) const - { - } + { + } ClassDef(MReportTrigger, 1) // Class for TRIGGER-REPORT information Index: /trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc =================================================================== --- /trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc (revision 3147) +++ /trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc (revision 3148) @@ -24,12 +24,11 @@ ////////////////////////////////////////////////////////////////////////////// -// // -// Cubic Spline Interpolation // -// // +// +// Cubic Spline Interpolation +// ////////////////////////////////////////////////////////////////////////////// - #include "MCubicCoeff.h" -#include "TMath.h" +#include #include "MLog.h" @@ -45,5 +44,4 @@ // where x is the independent variable, 2 and 3 are exponents // - MCubicCoeff::MCubicCoeff(Double_t x, Double_t xNext, Double_t y, Double_t yNext, Double_t a, Double_t b, Double_t c) : @@ -51,10 +49,10 @@ { fH = fXNext - fX; - if(!EvalMinMax()) - { - gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl; - fMin = 0; - fMax = 0; - } + if (EvalMinMax()) + return; + + gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl; + fMin = 0; + fMax = 0; } @@ -66,6 +64,6 @@ Double_t MCubicCoeff::Eval(Double_t x) { - Double_t dx = x - fX; - return (fY+dx*(fC+dx*(fB+dx*fA))); + const Double_t dx = x - fX; + return fY + dx*(fC + dx*(fB + dx*fA)); } @@ -79,71 +77,82 @@ Bool_t MCubicCoeff::EvalMinMax() { - fMin = fMax = fY; - fAbMin = fAbMax = fX; + fMin = fY; + fMax = fY; + + fAbMin = fX; + fAbMax = fX; + if (fYNext < fMin) { - fMin = fYNext; + fMin = fYNext; fAbMin = fXNext; } if (fYNext > fMax) { - fMax = fYNext; + fMax = fYNext; fAbMax = fXNext; } - Double_t tempMinMax; - Double_t delta = 4.0*fB*fB - 12.0*fA*fC; - if (delta >= 0.0 && fA != 0.0) - { - Double_t sqrtDelta = sqrt(delta); - Double_t xPlus = (-2.0*fB + sqrtDelta)/(6.0*fA); - Double_t xMinus = (-2.0*fB - sqrtDelta)/(6.0*fA); - if (xPlus >= 0.0 && xPlus <= fH) - { - tempMinMax = this->Eval(fX+xPlus); - if (tempMinMax < fMin) - { - fMin = tempMinMax; - fAbMin = fX + xPlus; - } - if (tempMinMax > fMax) - { - fMax = tempMinMax; - fAbMax = fX + xPlus; - } - } - if (xMinus >= 0.0 && xMinus <= fH) - { - tempMinMax = this->Eval(fX+xMinus); - if (tempMinMax < fMin) - { - fMin = tempMinMax; - fAbMin = fX + xMinus; - } - if (tempMinMax > fMax) - { - fMax = tempMinMax; - fAbMax = fX + xMinus; - } - } - } -/* if fA is zero then we have only one possible solution */ - else if (fA == 0.0 && fB != 0.0) - { - Double_t xSolo = - (fC/(2.0*fB)); - if (xSolo >= 0.0 && xSolo <= fH) - { - tempMinMax = this->Eval(fX+xSolo); - if (tempMinMax < fMin) - { - fMin = tempMinMax; - fAbMin = fX + xSolo; - } - if (tempMinMax > fMax) - { - fMax = tempMinMax; - fAbMax = fX + xSolo; - } - } - } + + const Double_t delta = fB*fB*4 - fA*fC*12; + if (delta >= 0 && fA != 0) + { + const Double_t sqrtDelta = TMath::Sqrt(delta); + + const Double_t xPlus = (-fB*2 + sqrtDelta)/(fA*6); + const Double_t xMinus = (-fB*2 - sqrtDelta)/(fA*6); + + if (xPlus >= 0 && xPlus <= fH) + { + const Double_t tempMinMax = Eval(fX+xPlus); + if (tempMinMax < fMin) + { + fMin = tempMinMax; + fAbMin = fX + xPlus; + } + if (tempMinMax > fMax) + { + fMax = tempMinMax; + fAbMax = fX + xPlus; + } + } + if (xMinus >= 0 && xMinus <= fH) + { + const Double_t tempMinMax = Eval(fX+xMinus); + if (tempMinMax < fMin) + { + fMin = tempMinMax; + fAbMin = fX + xMinus; + } + if (tempMinMax > fMax) + { + fMax = tempMinMax; + fAbMax = fX + xMinus; + } + } + return kTRUE; + } + + /* if fA is zero then we have only one possible solution */ + if (fA == 0 && fB != 0) + { + const Double_t xSolo = -fC/(fB*2); + + if (xSolo < 0 || xSolo > fH) + return kTRUE; + + const Double_t tempMinMax = Eval(fX+xSolo); + if (tempMinMax < fMin) + { + fMin = tempMinMax; + fAbMin = fX + xSolo; + } + if (tempMinMax > fMax) + { + fMax = tempMinMax; + fAbMax = fX + xSolo; + } + return kTRUE; + } + return kTRUE; } @@ -157,58 +166,45 @@ // There could be three or one real solutions // - Short_t MCubicCoeff::FindCardanRoot(Double_t y, Double_t *x) { - - Short_t whichRoot = -1; - - Double_t a = fB/fA; - Double_t b = fC/fA; - Double_t c = (fY - y)/fA; - - Double_t q = (a*a-3.0*b)/9.0; - Double_t r = (2.0*a*a*a-9.0*a*b+27.0*c)/54.0; - - Double_t aOver3 = a/3.0; - Double_t r2 = r*r; - Double_t q3 = q*q*q; + const Double_t a = fB/fA; + const Double_t b = fC/fA; + const Double_t c = (fY - y)/fA; + + const Double_t q = (a*a - b*3)/9; + const Double_t r = (a*a*a*2 - a*b*9 + c*27)/54; + + const Double_t aOver3 = a/3; + const Double_t r2 = r*r; + const Double_t q3 = q*q*q; if (r2 < q3) //3 real sol { - Double_t sqrtQ = TMath::Sqrt(q); - Double_t min2SqQ = -2.0*sqrtQ; - Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ)); - - x[0] = min2SqQ * TMath::Cos(theta/3.0) - aOver3; - x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3.0) - aOver3; - x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3.0) - aOver3; + const Double_t sqrtQ = TMath::Sqrt(q); + const Double_t min2SqQ = -sqrtQ*2; + const Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ)); + + x[0] = min2SqQ * TMath::Cos(theta/3) - aOver3; + x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3) - aOver3; + x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3) - aOver3; + for (Int_t i = 0; i < 3; i++) - if (x[i] >= 0.0 && x[i] <= fH) + if (x[i] >= 0 && x[i] <= fH) { - x[i] = x[i] + fX; - whichRoot = i; - break; + x[i] += fX; + return i; } - } - else //only 1 real sol - { - Double_t s; - if (r == 0.0) - s = 0.0; - else if (r < 0.0) - s = TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0); - else // r > 0.0 - s = - TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0); - if (s == 0.0) - x[0] = - aOver3; - else - x[0] = (s + q/s) - aOver3; - if (x[0] >= 0.0 && x[0] <= fH) - { - x[0] = x[0] + fX; - whichRoot = 0; - } - } - return whichRoot; + return -1; + } + + const Double_t s = r==0 ? 0 : -TMath::Sign(TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3), 1./3), r); + + x[0] = s==0 ? - aOver3 : (s + q/s) - aOver3; + + if (x[0] < 0 || x[0] > fH) + return -1; + + x[0] += fX; + return 0; } @@ -220,7 +216,4 @@ Bool_t MCubicCoeff :: IsIn(Double_t x) { - if (x >= fX && x <= fXNext) - return kTRUE; - else - return kFALSE; -} + return x >= fX && x <= fXNext; +} Index: /trunk/MagicSoft/Mars/mtools/MCubicSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/mtools/MCubicSpline.cc (revision 3147) +++ /trunk/MagicSoft/Mars/mtools/MCubicSpline.cc (revision 3148) @@ -24,18 +24,17 @@ ////////////////////////////////////////////////////////////////////////////// -// // -// Cubic Spline Interpolation // -// // +// +// Cubic Spline Interpolation +// ////////////////////////////////////////////////////////////////////////////// - #include "MCubicSpline.h" -#include "MCubicCoeff.h" - -#include "TMath.h" +#include #include "MLog.h" #include "MLogManip.h" +#include "MCubicCoeff.h" + ClassImp(MCubicSpline); @@ -48,6 +47,5 @@ // MCubicSpline::MCubicSpline(Byte_t *y, Byte_t *x, Bool_t areAllEq, - Int_t n, Double_t begSD, Double_t endSD): - fN(n) + Int_t n, Double_t begSD, Double_t endSD) { Init(y,x,areAllEq,n,begSD,endSD); @@ -62,5 +60,4 @@ { Byte_t x[]={0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,0x0C,0x0D,0x0E}; - fN = 15; Init(y,x,kTRUE,15,0.0,0.0); } @@ -75,58 +72,55 @@ { - Double_t *temp = new Double_t[fN-1]; - Double_t *ysd = new Double_t[fN-1]; - Double_t p,h; - MCubicCoeff *tempCoeff; - fCoeff = new TObjArray(fN-1,0); + Double_t *temp = new Double_t[n-1]; + Double_t *ysd = new Double_t[n-1]; + + fCoeff = new TObjArray(n-1,0); + + ysd[0] =begSD; + temp[0] =begSD; + ysd[n-1]=endSD; + Double_t h = x[1]-x[0]; + if (areAllEq) { - h = x[1]-x[0]; - ysd[0]=temp[0]=begSD; - ysd[n-1]=endSD; for(Int_t i = 1; i < n-1; i++) { - p = 0.5*ysd[i-1]+2.0; - ysd[i] = (-0.5)/p; - temp[i] = (y[i+1]-2*y[i]+y[i-1])/h; - temp[i] = (6.0*temp[i]/h-0.5*temp[i-1])/p; - } - for(Int_t i = n-2; i > 0; i--) - ysd[i]=ysd[i]*ysd[i+1]+temp[i]; - for(Int_t i = 0; i < n-1; i++) - { - tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h), - ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6); - fCoeff->AddAt(tempCoeff,i); - } - delete [] temp; - delete [] ysd; + const Double_t p = ysd[i-1]/2+2; + + ysd[i] = -0.5/p; + temp[i] = (y[i+1] - y[i]*2 + y[i-1])/h; + temp[i] = (temp[i]*6/h-temp[i-1]/2)/p; + } } else { - Double_t sig; - ysd[0]=temp[0]=begSD; - ysd[n-1]=endSD; for(Int_t i = 1; i < n-1; i++) { - sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]); - p = sig*ysd[i-1]+2.0; - ysd[i] = (sig-1.0)/p; + const Double_t sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]); + + const Double_t p = sig*ysd[i-1]+2; + + ysd[i] = (sig-1.0)/p; temp[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]); - temp[i] = (6.0*temp[i]/(x[i+1]-x[i-1])-sig*temp[i-1])/p; + temp[i] = (temp[i]*6/(x[i+1]-x[i-1])-sig*temp[i-1])/p; } - for(Int_t i = n-2; i > 0; i--) - ysd[i]=ysd[i]*ysd[i+1]+temp[i]; - for(Int_t i = 0; i < n-1; i++) - { - h = x[i+1]-x[i]; - tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h), - ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6); - fCoeff->AddAt(tempCoeff,i); - } - delete [] temp; - delete [] ysd; - } + } + + for(Int_t i = n-2; i > 0; i--) + ysd[i] = ysd[i]*ysd[i+1] + temp[i]; + + for(Int_t i = 0; i < n-1; i++) + { + if (!areAllEq) + h = x[i+1]-x[i]; + + MCubicCoeff *c = new MCubicCoeff(x[i], x[i+1], y[i], y[i+1], (ysd[i+1]-ysd[i])/(h*6), + ysd[i]/2, (y[i+1]-y[i])/h-(h*(ysd[i+1]+ysd[i]*2))/6); + fCoeff->AddAt(c, i); + } + + delete [] temp; + delete [] ysd; } @@ -134,4 +128,5 @@ { fCoeff->Delete(); + delete fCoeff; } @@ -142,9 +137,15 @@ Double_t MCubicSpline :: Eval(Double_t x) { - for (Int_t i = 0; i < fN-1; i++) - if (((MCubicCoeff*)fCoeff->At(i))->IsIn(x)) - return ((MCubicCoeff*)fCoeff->At(i))->Eval(x); + const Int_t n = fCoeff->GetSize()-1; + for (Int_t i = 0; i < n; i++) + { + MCubicCoeff *c = (MCubicCoeff*)fCoeff->UncheckedAt(i); + if (c->IsIn(x)) + return c->Eval(x); + } + gLog << warn << "Cannot evaluate Spline at " << x << "; returning 0"; - return 0.0; + + return 0; } @@ -155,12 +156,11 @@ Double_t MCubicSpline :: EvalMax() { - Double_t temp_max; - Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax(); - for (Int_t i = 1; i < fN-1; i++) - { - temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax(); - if (temp_max > max) - max = temp_max; - } + Double_t max = -FLT_MAX; + + TIter Next(fCoeff); + MCubicCoeff *c; + while ((c=(MCubicCoeff*)Next())) + max = TMath::Max(max, c->GetMax()); + return max; } @@ -172,12 +172,11 @@ Double_t MCubicSpline :: EvalMin() { - Double_t temp_min; - Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin(); - for (Int_t i = 1; i < fN-1; i++) - { - temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin(); - if (temp_min < min) - min = temp_min; - } + Double_t min = FLT_MAX; + + TIter Next(fCoeff); + MCubicCoeff *c; + while ((c=(MCubicCoeff*)Next())) + min = TMath::Min(min, c->GetMax()); + return min; } @@ -189,17 +188,22 @@ Double_t MCubicSpline :: EvalAbMax() { - Double_t temp_max; - Double_t abMax = ((MCubicCoeff*)fCoeff->At(0))->GetAbMax(); - Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax(); - for (Int_t i = 1; i < fN-1; i++) - { - temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax(); - if (temp_max > max) - { - max = temp_max; - abMax = ((MCubicCoeff*)fCoeff->At(i))->GetAbMax(); - } - } - return abMax; + Double_t max = -FLT_MAX; + + TIter Next(fCoeff); + + MCubicCoeff *c; + MCubicCoeff *cmax=0; + + while ((c=(MCubicCoeff*)Next())) + { + const Double_t temp = c->GetMax(); + if (temp <= max) + continue; + + max = temp; + cmax = c; + } + + return cmax ? cmax->GetAbMax() : -FLT_MAX; } @@ -210,17 +214,22 @@ Double_t MCubicSpline :: EvalAbMin() { - Double_t temp_min; - Double_t abMin = ((MCubicCoeff*)fCoeff->At(0))->GetAbMin(); - Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin(); - for (Int_t i = 1; i < fN-1; i++) - { - temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin(); - if (temp_min < min) - { - min = temp_min; - abMin = ((MCubicCoeff*)fCoeff->At(i))->GetAbMin(); - } - } - return abMin; + Double_t min = FLT_MAX; + + TIter Next(fCoeff); + + MCubicCoeff *c; + MCubicCoeff *cmin=0; + + while ((c=(MCubicCoeff*)Next())) + { + const Double_t temp = c->GetMax(); + if (temp >= min) + continue; + + min = temp; + cmin = c; + } + + return cmin ? cmin->GetAbMax() : FLT_MAX; } @@ -233,41 +242,37 @@ Double_t MCubicSpline :: FindVal(Double_t y, Double_t x0, Char_t direction = 'l') { - Short_t whichRoot; - Double_t tempRoot; - Double_t *roots = new Double_t[3]; - - for (Int_t i = 0; i < fN-1; i++) - { - if(((MCubicCoeff*)fCoeff->At(i))->IsIn(x0)) - { - if(direction == 'l') - { - for (Int_t j = i; j >= 0; j--) - { - whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots); - if (whichRoot >= 0 ) - { - tempRoot = roots[whichRoot]; - delete [] roots; - return tempRoot; - } - } - } - if(direction == 'r') - { - for (Int_t j = i; j < fN-1; j++) - { - whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots); - if (whichRoot >= 0) - { - tempRoot = roots[whichRoot]; - delete [] roots; - return tempRoot; - } - } - } - } - } + Double_t roots[3] = { 0, 0, 0 }; + + const Int_t n = fCoeff->GetSize()-1; + + for (Int_t i = 0; i < n; i++) + { + if (!((MCubicCoeff*)fCoeff->At(i))->IsIn(x0)) + continue; + + switch (direction) + { + case 'l': + for (Int_t j = i; j >= 0; j--) + { + const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots); + if (whichRoot >= 0 ) + return roots[whichRoot]; + } + break; + + case 'r': + for (Int_t j = i; j < n; j++) + { + const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots); + if (whichRoot >= 0) + return roots[whichRoot]; + } + break; + } + } + gLog << warn << "Nothing found calling MCubicSpline :: FindVal(), returning 0" << endl; - return 0.0; -} + + return 0; +} Index: /trunk/MagicSoft/Mars/mtools/MCubicSpline.h =================================================================== --- /trunk/MagicSoft/Mars/mtools/MCubicSpline.h (revision 3147) +++ /trunk/MagicSoft/Mars/mtools/MCubicSpline.h (revision 3148) @@ -16,5 +16,5 @@ private: TObjArray *fCoeff; //array of the coefficients - Int_t fN; //number of points + void Init(Byte_t *y, Byte_t *x, Bool_t areAllEq, Int_t n, Double_t begSD, Double_t endSD);