//-------------------------------------------------------------------------
// Implementation of the AliSplineFit class
-// The class performs a spline fit on an incoming TGraph. The graph is
-// divided into several parts (identified by knots between each part).
+// The class performs a spline fit on an incoming TGraph. The graph is
+// divided into several parts (identified by knots between each part).
// Spline fits are performed on each part. According to user parameters,
-// the function, first and second derivative are requested to be continuous
+// the function, first and second derivative are requested to be continuous
// at each knot.
// Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
// Adjustments by Haavard Helstrup, Haavard.Helstrup@cern.ch
//-------------------------------------------------------------------------
-#include <TMath.h>
#include "AliSplineFit.h"
-
+
ClassImp(AliSplineFit)
-TLinearFitter*
-AliSplineFit::fitterStatic()
+TLinearFitter* AliSplineFit::fitterStatic()
{
static TLinearFitter* fit = new TLinearFitter(4,"pol3","");
return fit;
fBDump(kFALSE),
fGraph (0),
fNmin (0),
+ fMinPoints(0),
fSigma (0),
fMaxDelta (0),
fN0 (0),
fBDump (source.fBDump),
fGraph (source.fGraph),
fNmin (source.fNmin),
+ fMinPoints (source.fMinPoints),
fSigma (source.fSigma),
fMaxDelta (source.fMaxDelta),
fN0 (source.fN0),
+ fParams (0),
+ fCovars (0),
+ fIndex (0),
fN (source.fN),
- fChi2 (source.fChi2)
+ fChi2 (0.0),
+ fX (0),
+ fY0 (0),
+ fY1 (0),
+ fChi2I (source.fChi2I)
{
//
// Copy constructor
fParams = (TClonesArray*)source.fParams->Clone();
fCovars = (TClonesArray*)source.fCovars->Clone();
for (Int_t i=0; i<fN0; i++) fIndex[i] = source.fIndex[i];
-
+
fX = new Double_t[fN];
fY0 = new Double_t[fN];
fY1 = new Double_t[fN];
AliSplineFit& AliSplineFit::operator=(const AliSplineFit& source){
//
// assignment operator
-//
+//
if (&source == this) return *this;
//
// reassign memory as previous fit could have a different size
//
- if ( fN0 != source.fN0) {
+ if ( fN0 != source.fN0) {
delete fParams;
delete fCovars;
fParams = new TClonesArray("TVectorD",fN0);
fCovars = new TClonesArray("TMatrixD",fN0);
}
- if ( fN != source.fN) {
+ if ( fN != source.fN) {
delete []fX;
delete []fY0;
delete []fY1;
delete []fChi2I;
- fN = source.fN;
+ fN = source.fN;
fX = new Double_t[fN];
fY0 = new Double_t[fN];
fY1 = new Double_t[fN];
}
// use copy constructor (without reassigning memory) to copy values
-
+
new (this) AliSplineFit(source);
-
+
return *this;
}
-
+
AliSplineFit::~AliSplineFit(){
//
// destructor. Don't delete fGraph, as this normally comes as input parameter
// = 3: 3rd derivative
//
// a2 = -(3*a0 -3*b0 + 2*a1*dx +b1*dx)/(dx*dx)
- // a3 = -(-2*a0+2*b0 - a1*dx - b1*dx)/(dx*dx*dx)
+ // a3 = -(-2*a0+2*b0 - a1*dx - b1*dx)/(dx*dx*dx)
Int_t index = TMath::BinarySearch(fN,fX,x);
if (index<0) index =0;
//
Double_t dx = x-fX[index];
Double_t dxc = fX[index+1]-fX[index];
+ if (dxc<=0) return fY0[index];
Double_t y0 = fY0[index];
Double_t y1 = fY1[index];
Double_t y01 = fY0[index+1];
}
TGraph * graph = new TGraph(npoints,time,value);
-
+
delete [] value;
delete [] time;
return graph;
delete [] time;
return graph;
}
-
+
TGraph * AliSplineFit::MakeGraph(Double_t xmin, Double_t xmax, Int_t npoints, Int_t deriv) const {
//
TGraph * AliSplineFit::MakeDiff(TGraph * graph0) const {
//
- // Make graph of difference to reference graph
+ // Make graph of difference to reference graph
//
Int_t npoints=graph0->GetN();
TH1F * AliSplineFit::MakeDiffHisto(TGraph * graph0) const {
//
- // Make histogram of difference to reference graph
+ // Make histogram of difference to reference graph
//
Int_t npoints=graph0->GetN();
- Float_t min=1e+33,max=-1e+33;
+ Float_t min=1e+39,max=-1e+39;
for (Int_t ip=0; ip<npoints; ip++){
Double_t x = graph0->GetX()[ip];
Double_t y = Eval(x,0)-graph0->GetY()[ip];
}else{
if (y<min) min=y;
if (y>max) max=y;
- }
+ }
}
TH1F *his = new TH1F("hdiff","hdiff", 100, min, max);
void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t maxDelta){
//
- // initialize knots + estimate sigma of noise + make initial parameters
+ // initialize knots + estimate sigma of noise + make initial parameters
//
//
fNmin = min;
fMaxDelta = maxDelta;
Int_t npoints = fGraph->GetN();
+
+ // Make simple spline if too few points in graph
+
+ if (npoints < fMinPoints ) {
+ CopyGraph();
+ return;
+ }
+
fN0 = (npoints/fNmin)+1;
Float_t delta = Double_t(npoints)/Double_t(fN0-1);
for (Int_t iKnot=0; iKnot<fN0; iKnot++){
TMatrixD & cov = *((TMatrixD*)fCovars->At(iKnot));
cov*=fSigma*fSigma;
- }
+ }
OptimizeKnots(iter);
fN = 0;
fChi2I = new Double_t[fN];
Int_t iKnot=0;
for (Int_t i=0; i<fN0; i++){
- if (fIndex[i]<0) continue;
+ if (fIndex[i]<0) continue;
if (iKnot>=fN) {
printf("AliSplineFit::InitKnots: Knot number > Max knot number\n");
break;
TVectorD * param = (TVectorD*) fParams->At(i);
fX[iKnot] = fGraph->GetX()[fIndex[i]];
fY0[iKnot] = (*param)(0);
- fY1[iKnot] = (*param)(1);
+ fY1[iKnot] = (*param)(1);
fChi2I[iKnot] = 0;
iKnot++;
}
TVectorD sPrevious = tPrevious*pPrevious+tNext*pNext;
TVectorD sKnot = tpKnot*sPrevious;
- TVectorD sNext = tpNext*sPrevious;
+ TVectorD sNext = tpNext*sPrevious;
TMatrixD csPrevious00(tPrevious, TMatrixD::kMult,cPrevious);
csPrevious00 *= tPrevious.T();
}
-
+
void AliSplineFit::SplineFit(Int_t nder){
//
// Cubic spline fit of graph
- //
+ //
// nder
// nder<0 - no continuity requirement
// =0 - continous 0 derivative
// =1 - continous 1 derivative
- // >1 - continous 2 derivative
+ // >1 - continous 2 derivative
//
if (!fGraph) return;
TGraph * graph = fGraph;
if (nder>1) nder=2;
Int_t nknots = fN;
+ if (nknots < 2 ) return;
Int_t npoints = graph->GetN();
+ if (npoints < fMinPoints ) return;
//
//
// spline fit
// each knot 4 parameters
- //
+ //
TMatrixD *pmatrix = 0;
TVectorD *pvalues = 0;
if (nder>1){
-void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, char * type){
+void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, Option_t * type){
//
// Interface to GraphSmooth
//
delete h2;
delete hS;
}
- delete pcstream;
+ delete pcstream;
}
void AliSplineFit::Cleanup(){
// deletes extra information to reduce amount of information stored on the data
// base
- delete fGraph; fGraph=0;
+ // delete fGraph; fGraph=0; // Don't delete fGraph -- input parameter
delete fParams; fParams=0;
delete fCovars; fCovars=0;
delete [] fIndex; fIndex=0;
delete [] fChi2I; fChi2I=0;
}
+
+
+void AliSplineFit::CopyGraph() {
+ //
+ // enter graph points directly to fit parameters
+ // (to bee used when too few points are available)
+ //
+ fN = fGraph->GetN();
+ fX = new Double_t[fN];
+ fY0 = new Double_t[fN];
+ fY1 = new Double_t[fN];
+ for (Int_t i=0; i<fN; i++ ) {
+ fX[i] = fGraph->GetX()[i];
+ fY0[i] = fGraph->GetY()[i];
+ fY1[i] = 0;
+ }
+}