X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliSplineFit.cxx;h=ae375245c013d1951d4405a0e958ea9a7310b3ef;hb=ce812f4fc89a4aa5e40ec7ad182063b1dabb2deb;hp=52c096ffacc77a3ae46c35f7ed9e7154550395af;hpb=0dd3a2ac0e41f5a75f1f2112e9b3f151cff308a7;p=u%2Fmrichter%2FAliRoot.git

diff --git a/STEER/AliSplineFit.cxx b/STEER/AliSplineFit.cxx
index 52c096ffacc..ae375245c01 100644
--- a/STEER/AliSplineFit.cxx
+++ b/STEER/AliSplineFit.cxx
@@ -13,17 +13,33 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+//-------------------------------------------------------------------------
+//                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).
+//   Spline fits are performed on each part. According to user parameters,
+//   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 "AliSplineFit.h"
- 
-ClassImp(AliSplineFit);
 
-TLinearFitter AliSplineFit::fitterStatic = TLinearFitter(4,"pol3","");
+ClassImp(AliSplineFit)
+
+TLinearFitter* AliSplineFit::fitterStatic()
+{
+  static TLinearFitter* fit = new TLinearFitter(4,"pol3","");
+  return fit;
+}
 
 AliSplineFit::AliSplineFit() :
   fBDump(kFALSE),
   fGraph (0),
   fNmin (0),
+  fMinPoints(0),
   fSigma (0),
   fMaxDelta (0),
   fN0 (0),
@@ -48,11 +64,19 @@ AliSplineFit::AliSplineFit(const AliSplineFit& source) :
   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
@@ -63,7 +87,7 @@ AliSplineFit::AliSplineFit(const AliSplineFit& source) :
   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];
@@ -72,19 +96,20 @@ AliSplineFit::AliSplineFit(const AliSplineFit& source) :
     fX[i]  = source.fX[i];
     fY0[i] = source.fY0[i];
     fY1[i] = source.fY1[i];
+    fChi2I[i] = source.fChi2I[i];
   }
 }
 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;
@@ -95,13 +120,13 @@ AliSplineFit& AliSplineFit::operator=(const AliSplineFit& source){
     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];
@@ -109,13 +134,13 @@ AliSplineFit& AliSplineFit::operator=(const AliSplineFit& source){
   }
 
 // 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
@@ -138,7 +163,7 @@ Double_t   AliSplineFit::Eval(Double_t x, Int_t deriv) const{
   //         = 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;
@@ -146,6 +171,7 @@ Double_t   AliSplineFit::Eval(Double_t x, Int_t deriv) const{
   //
   Double_t dx   = x-fX[index];
   Double_t dxc  = fX[index+1]-fX[index];
+  if (dxc<=0)    return 0;
   Double_t y0   = fY0[index];
   Double_t y1   = fY1[index];
   Double_t y01  = fY0[index+1];
@@ -205,7 +231,7 @@ TGraph * AliSplineFit::GenerGraph(Int_t npoints, Double_t fraction, Double_t s1,
   }
 
   TGraph * graph = new TGraph(npoints,time,value);
- 
+
   delete [] value;
   delete [] time; 
   return graph;  
@@ -230,7 +256,7 @@ TGraph * AliSplineFit::GenerNoise(TGraph * graph0, Double_t sigma0){
   delete [] time;
   return graph;  
 }
- 
+
 
 TGraph * AliSplineFit::MakeGraph(Double_t xmin, Double_t xmax, Int_t npoints, Int_t deriv) const {
   //
@@ -258,7 +284,7 @@ TGraph * AliSplineFit::MakeGraph(Double_t xmin, Double_t xmax, Int_t npoints, In
 
 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();
@@ -278,7 +304,7 @@ TGraph * AliSplineFit::MakeDiff(TGraph * graph0) const {
 
 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();
@@ -292,7 +318,7 @@ TH1F * AliSplineFit::MakeDiffHisto(TGraph * graph0) const {
     }else{
       if (y<min) min=y;
       if (y>max) max=y;
-    }    
+    }
   }
 
   TH1F *his = new TH1F("hdiff","hdiff", 100, min, max);
@@ -309,7 +335,7 @@ TH1F * AliSplineFit::MakeDiffHisto(TGraph * graph0) const {
 
 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
   //
   //
 
@@ -318,6 +344,14 @@ void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t max
   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);
 
@@ -361,7 +395,7 @@ void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t max
   for (Int_t iKnot=0; iKnot<fN0; iKnot++){
     TMatrixD & cov = *((TMatrixD*)fCovars->At(iKnot));
     cov*=fSigma*fSigma;
-  }  
+  }
   OptimizeKnots(iter);
 
   fN = 0;
@@ -372,7 +406,7 @@ void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t max
   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;
@@ -380,7 +414,7 @@ void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t max
     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++;
   }
@@ -436,7 +470,6 @@ Bool_t   AliSplineFit::RefitKnot(Int_t iKnot){
   //
   //
   //
-  const Double_t kEpsilon = 1.e-7;
 
   Int_t iPrevious=(iKnot>0)  ?iKnot-1: 0;
   Int_t iNext    =(iKnot<fN0)?iKnot+1: fN0-1;
@@ -446,7 +479,7 @@ Bool_t   AliSplineFit::RefitKnot(Int_t iKnot){
   if (iNext>=fN0) iNext=fN0-1;
   
   Double_t startX = fGraph->GetX()[fIndex[iKnot]]; 
-  AliSplineFit::fitterStatic.ClearPoints();
+  AliSplineFit::fitterStatic()->ClearPoints();
   Int_t indPrev = fIndex[iPrevious];
   Int_t indNext = fIndex[iNext];
   Double_t *graphX = fGraph->GetX();
@@ -468,8 +501,8 @@ Bool_t   AliSplineFit::RefitKnot(Int_t iKnot){
 //    ePoint[indVec] =  fSigma+TMath::Abs(y)*kEpsilon;
 //    AliSplineFit::fitterStatic.AddPoint(&dxl,y,fSigma+TMath::Abs(y)*kEpsilon);
   }
-  AliSplineFit::fitterStatic.AssignData(nPoints,1,xPoint,yPoint,ePoint);
-  AliSplineFit::fitterStatic.Eval();
+  AliSplineFit::fitterStatic()->AssignData(nPoints,1,xPoint,yPoint,ePoint);
+  AliSplineFit::fitterStatic()->Eval();
 
 //  delete temporary arrays
 
@@ -477,8 +510,8 @@ Bool_t   AliSplineFit::RefitKnot(Int_t iKnot){
   
   TMatrixD   * covar = (TMatrixD*)fCovars->At(iKnot);
   TVectorD   * param = (TVectorD*)fParams->At(iKnot);
-  AliSplineFit::fitterStatic.GetParameters(*param);
-  AliSplineFit::fitterStatic.GetCovarianceMatrix(*covar);
+  AliSplineFit::fitterStatic()->GetParameters(*param);
+  AliSplineFit::fitterStatic()->GetCovarianceMatrix(*covar);
   return 0;
 }
 
@@ -536,7 +569,7 @@ Float_t AliSplineFit::CheckKnot(Int_t 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();
@@ -583,27 +616,29 @@ Float_t AliSplineFit::CheckKnot(Int_t iKnot){
 
 
 }
-                  
+
 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){
@@ -857,7 +892,7 @@ void AliSplineFit::MakeKnots0(TGraph * graph, Double_t maxdelta, Int_t minpoints
 
 
 
-void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, char * type){
+void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, Option_t * type){
   //
   // Interface to GraphSmooth
   //
@@ -981,5 +1016,34 @@ void AliSplineFit::Test(Int_t npoints, Int_t ntracks, Float_t snoise){
     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;   // 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;
+  }
 }