Index: trunk/MagicSoft/Mars/macros/spline.C =================================================================== --- trunk/MagicSoft/Mars/macros/spline.C (revision 3022) +++ trunk/MagicSoft/Mars/macros/spline.C (revision 3022) @@ -0,0 +1,109 @@ +/* This macro is defined as a class for debugging (with CInt) reasons + +To use it at the root prompt type: + +root [0] .L spline.C +root [1] TestSpline::sp() +*/ +/* Example of a spline. You can use it as Test. If you think there are some + bugs in the MCubicSpline class please mail to: raducci@fisica.uniud.it */ + +class TestSpline +{ +public: + void sp(); +}; + +void TestSpline::sp() +{ +gROOT -> Reset(); + +//Here are defined the points. X goes from 0 to 14 (as the fadc slices...) +//Y are arbitrary values + +/* User Change */ +Byte_t y[]={0x0F,0x10,0x2F,0x7F,0xAA,0x6C,0x14,0x13,0x15,0x18,0x21,0x12,0x11,0x14,0x13}; +/* End user Change */ +Byte_t x[]={0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,0x0C,0x0D,0x0E}; + +/*This cast is needed only to show graphically the output. Don' t needed if you + use the spline to calc the arrival times */ + +Int_t *newX = new Int_t[15]; +Int_t *newY = new Int_t[15]; + +for (Int_t i = 0; i < 15; i++) +{ + newX[i] = (Int_t) x[i]; + newY[i] = (Int_t) y[i]; +} + +//Canvas to display output +TCanvas *c = new TCanvas ("c1","c1",800,600); + +//Graph containting only the points (displayed as stars) +TGraph *g1 = new TGraph(15,newX,newY); + +g1 -> Draw("A*"); + +//Spline constructor(specialized for 15 slices using Bytes as values. There exist another constructor. +MCubicSpline *s = new MCubicSpline(y); + +//*spline and *ab are two arrays containing some values evaluated from the spline +Double_t *spline = new Double_t[139]; +Double_t *ab = new Double_t[139]; +Double_t step = 0.0; + +for (Int_t i = 0; i < 139; i++) +{ + spline[i] = s->Eval(step); + ab[i] = step; + step += 0.1; +} + +//Graph of the sline. The points calculated are joined with a red line. If the stars lie +//on the red line, then the Spline class is working properly +TGraph *g2 = new TGraph(139,ab,spline); + +g2 -> SetLineColor(2); +g2 -> Draw("C"); + +//Maximum and minimum evaluation +Double_t *mm = new Double_t[2]; +Double_t *abmm = new Double_t[2]; + +mm[0] = s->EvalMin(); +mm[1] = s->EvalMax(); +abmm[0] = s->EvalAbMin(); +abmm[1] = s->EvalAbMax(); + +//Display the max and the min using two black squares. If they lie on the max and min +//of the red line, then the Spline class is working properly + +TGraph *g3 = new TGraph(2,abmm,mm); + +g3 -> SetMarkerStyle(21); +g3 -> Draw("P"); + +//Test of the Cardan formula. Find the point(abval) where the Spline value is val +Double_t val = 82.5; +Double_t abval = s->FindVal(val, abmm[1], 'l'); + +//Display this point. Again, if it lies on the red line, then we are right. +//It's a black triangle + +TGraph *g4 = new TGraph(1,&abval,&val); + +g4 -> SetMarkerStyle(22); +g4 -> Draw("P"); + +//Free memory +s->~MCubicSpline(); +delete [] newX; +delete [] newY; +delete [] spline; +delete [] ab; +delete [] mm; +delete [] abmm; + +} Index: trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc (revision 3022) +++ trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc (revision 3022) @@ -0,0 +1,226 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Sebastian Raducci 01/2004 +! +! Copyright: MAGIC Software Development, 2001-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// Cubic Spline Interpolation // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MCubicCoeff.h" + +#include "TMath.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCubicCoeff); + +using namespace std; + +//---------------------------------------------------------------------------- +// +// Constructor (The spline is: fA(x-fX)3+fB(x-fX)2+fC(x-fX)+fY +// 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) : + fX(x), fXNext(xNext), fA(a), fB(b), fC(c), fY(y), fYNext(yNext) +{ + fH = fXNext - fX; + if(!EvalMinMax()) + { + gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl; + fMin = 0; + fMax = 0; + } +} + +//---------------------------------------------------------------------------- +// +// Evaluate the spline at a given point +// + +Double_t MCubicCoeff::Eval(Double_t x) +{ + Double_t dx = x - fX; + return (fY+dx*(fC+dx*(fB+dx*fA))); +} + +//---------------------------------------------------------------------------- +// +// Find min and max using derivatives. The min and max could be at the begin +// or at the end of the interval or somewhere inside the interval (in this case +// a comparison between the first derivative and zero is made) +// The first derivative coefficients are obviously: 3*fA, 2*fB, fC +// +Bool_t MCubicCoeff::EvalMinMax() +{ + fMin = fMax = fY; + fAbMin = fAbMax = fX; + if (fYNext < fMin) + { + fMin = fYNext; + fAbMin = fXNext; + } + if (fYNext > fMax) + { + 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; + } + } + } + return kTRUE; +} +//------------------------------------------------------------------------- +// +// Given y finds x using the cubic (cardan) formula. +// +// we consider the following form: x3 + ax2 + bx + c = 0 where +// a = fB/fA, b = fC/fA, c = (fY - y)/fA +// +// 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; + + 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; + for (Int_t i = 0; i < 3; i++) + if (x[i] >= 0.0 && x[i] <= fH) + { + x[i] = x[i] + fX; + whichRoot = i; + break; + } + } + 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 true if x is in this interval +// + +Bool_t MCubicCoeff :: IsIn(Double_t x) +{ + if (x >= fX && x <= fXNext) + return kTRUE; + else + return kFALSE; +} Index: trunk/MagicSoft/Mars/mtools/MCubicCoeff.h =================================================================== --- trunk/MagicSoft/Mars/mtools/MCubicCoeff.h (revision 3022) +++ trunk/MagicSoft/Mars/mtools/MCubicCoeff.h (revision 3022) @@ -0,0 +1,42 @@ +#ifndef MARS_MCubicCoeff +#define MARS_MCubicCoeff + +#ifndef MARS_MAGIC +#include "MAGIC.h" +#endif + +class MCubicCoeff : public TObject +{ + private: + Double_t fX; // abscissa + Double_t fXNext; // abscissa of the next point + Double_t fA; // 3rd order coeff + Double_t fB; // 2nd order coeff + Double_t fC; // 1st order coeff + Double_t fY; // constant term + Double_t fYNext; // value in the next point + Double_t fH; // interval width + Double_t fMin; // minimum value + Double_t fMax; // maximum value + Double_t fAbMin; // abscissa of the min + Double_t fAbMax; // abscissa of the max + public: + MCubicCoeff(){} + MCubicCoeff(Double_t x, Double_t xNext, Double_t y, Double_t yNext, + Double_t a, Double_t b, Double_t c); + Double_t GetA() { return fA; } + Double_t GetB() { return fB; } + Double_t GetC() { return fC; } + Double_t GetMin() { return fMin; } + Double_t GetMax() { return fMax; } + Double_t GetAbMin() { return fAbMin; } + Double_t GetAbMax() { return fAbMax; } + Double_t Eval(Double_t x); //Evaluate the spline at a point x + Bool_t EvalMinMax(); //Finds min & max + Short_t FindCardanRoot(Double_t y, Double_t *x); //Evaluate the abscissa of the spline given y + Bool_t IsIn(Double_t x); + + ClassDef(MCubicCoeff, 0) //Class to contain spline coefficients +}; + +#endif Index: trunk/MagicSoft/Mars/mtools/MCubicSpline.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MCubicSpline.cc (revision 3022) +++ trunk/MagicSoft/Mars/mtools/MCubicSpline.cc (revision 3022) @@ -0,0 +1,273 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Sebastian Raducci 01/2004 +! +! Copyright: MAGIC Software Development, 2001-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// Cubic Spline Interpolation // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MCubicSpline.h" + +#include "MCubicCoeff.h" + +#include "TMath.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCubicSpline); + +using namespace std; + +//--------------------------------------------------------------------------- +// +// Contructor +// +// +MCubicSpline::MCubicSpline(Byte_t *y, Byte_t *x, Bool_t areAllEq, + Int_t n, Double_t begSD, Double_t endSD): + fN(n) +{ + Init(y,x,areAllEq,n,begSD,endSD); +} + +//--------------------------------------------------------------------------- +// +// Constructor for FADC slice (only the FADC counts are needed) +// +// +MCubicSpline::MCubicSpline(Byte_t *y) +{ + 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); +} + +//--------------------------------------------------------------------------- +// +// Constructors common part +// +// +void MCubicSpline::Init(Byte_t *y, Byte_t *x, Bool_t areAllEq, + Int_t n, Double_t begSD, Double_t endSD) + +{ + 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); + + 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; + } + 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; + 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; + } + 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; + } +} + +MCubicSpline::~MCubicSpline() +{ + fCoeff->Delete(); +} + +//--------------------------------------------------------------------------- +// +// Evaluate the spline at a given point +// +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); + gLog << warn << "Cannot evaluate Spline at " << x << "; returning 0"; + return 0.0; +} + +//---------------------------------------------------------------------------- +// +// Search for max +// +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; + } + return max; +} + +//---------------------------------------------------------------------------- +// +// Search for min +// +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; + } + return min; +} + +//---------------------------------------------------------------------------- +// +// Search for abscissa of the max +// +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; +} + +//---------------------------------------------------------------------------- +// +// Search for abscissa of the min +// +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; +} + +//---------------------------------------------------------------------------- +// +// Finds the abscissa where the spline reaches y starting from x0 going in +// direction direction +// You have to give as input a starting point and a direction ("l" or "r") +// +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; + } + } + } + } + } + gLog << warn << "Nothing found calling MCubicSpline :: FindVal(), returning 0" << endl; + return 0.0; +} Index: trunk/MagicSoft/Mars/mtools/MCubicSpline.h =================================================================== --- trunk/MagicSoft/Mars/mtools/MCubicSpline.h (revision 3022) +++ trunk/MagicSoft/Mars/mtools/MCubicSpline.h (revision 3022) @@ -0,0 +1,35 @@ +#ifndef MARS_MCubicSpline +#define MARS_MCubicSpline + +#ifndef MARS_MAGIC +#include "MAGIC.h" +#endif + +#ifndef ROOT_TObjArray +#include "TObjArray.h" +#endif + +class MCubicCoeff; + +class MCubicSpline : public TObject +{ + 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); + + public: + MCubicSpline(Byte_t *y, Byte_t *x, Bool_t areAllEq, Int_t n, Double_t begSD=0.0, Double_t endSD=0.0); + MCubicSpline(Byte_t *y); + ~MCubicSpline(); + Double_t Eval(Double_t x); //Eval the spline at a point x + Double_t EvalMax(); //Eval the max + Double_t EvalMin(); //Eval the min + Double_t EvalAbMax(); //Eval the abscissa of the max + Double_t EvalAbMin(); //Eval the abscissa of the min + Double_t FindVal(Double_t y, Double_t x0, Char_t direction); //Finds the abscissa where the spline reaches y + + ClassDef(MCubicSpline, 0) //Class to contain spline coefficients +}; + +#endif Index: trunk/MagicSoft/Mars/mtools/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtools/Makefile (revision 3021) +++ trunk/MagicSoft/Mars/mtools/Makefile (revision 3022) @@ -36,4 +36,6 @@ MSimulatedAnnealing.cc \ MFFT.cc \ + MCubicCoeff.cc \ + MCubicSpline.cc \ MagicReversi.cc \ MagicSnake.cc \ @@ -42,4 +44,5 @@ MagicCivilization.cc \ MineSweeper.cc + SRCS = $(SRCFILES) Index: trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h (revision 3021) +++ trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h (revision 3022) @@ -11,4 +11,6 @@ #pragma link C++ class MHSimulatedAnnealing+; #pragma link C++ class MSimulatedAnnealing+; +#pragma link C++ class MCubicCoeff+; +#pragma link C++ class MCubicSpline+; #pragma link C++ class MineSweeper+;