/* ======================================================================== *\ ! ! * ! * 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; }