Methods to get derivatives added

sevdokim · Jan 28, 2016 · fa32159 · fa32159
1 parent dcb7787
commit fa32159
Show file tree

Hide file tree

Showing 2 changed files with 81 additions and 2 deletions.
diff --git a/STAT/AliNDLocalRegression.cxx b/STAT/AliNDLocalRegression.cxx
@@ -512,12 +512,12 @@ Double_t AliNDLocalRegression::Eval(Double_t *point ){
   TVectorD &vecParam = *((TVectorD*)fLocalFitParam->At(ibin));
   Double_t value=vecParam[0];
   if (!rangeOK) return value;
-  if (fUseBinNorm){
+  if (fUseBinNorm) {
     for (Int_t ipar=0; ipar<fNParameters; ipar++){
       Double_t delta=(point[ipar]-fBinCenter[ipar])/fBinWidth[ipar];
       value+=(vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
     }
-  }else{
+  } else {
     for (Int_t ipar=0; ipar<fNParameters; ipar++){
       Double_t delta=(point[ipar]-fBinCenter[ipar]);
       value+=(vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
@@ -547,6 +547,83 @@ Double_t AliNDLocalRegression::EvalError(Double_t *point ){
   return value;
 }
 
+Bool_t AliNDLocalRegression::Derivative(Double_t *point, Double_t *d)
+{
+  // fill d by partial derivatives
+  //
+  const Double_t almost0=0.00000001;
+  for (Int_t iDim=0; iDim<fNParameters; iDim++){
+    const TAxis* ax = fHistPoints->GetAxis(iDim);
+    if      (point[iDim]<=ax->GetXmin()) point[iDim] = ax->GetXmin()+almost0*ax->GetBinWidth(0);
+    else if (point[iDim]>=ax->GetXmax()) point[iDim] = ax->GetXmax()-almost0*ax->GetBinWidth(0);
+  }
+
+  Int_t ibin = fHistPoints->GetBin(point);
+  if (ibin>=fLocalFitParam->GetEntriesFast() ||
+      !fLocalFitParam->UncheckedAt(ibin) ) return kFALSE;
+  //
+  fHistPoints->GetBinContent(ibin,fBinIndex); 
+  for (Int_t idim=0; idim<fNParameters; idim++){
+    const TAxis* ax = fHistPoints->GetAxis(idim);
+    fBinCenter[idim] = ax->GetBinCenter(fBinIndex[idim]);
+    fBinWidth[idim] = ax->GetBinWidth(fBinIndex[idim]);
+  } 
+  TVectorD &vecParam = *((TVectorD*)fLocalFitParam->At(ibin));
+
+  if (fUseBinNorm){
+    for (Int_t ipar=0; ipar<fNParameters; ipar++){
+      Double_t delta=(point[ipar]-fBinCenter[ipar])/fBinWidth[ipar];
+      d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
+    }
+  }else{
+    for (Int_t ipar=0; ipar<fNParameters; ipar++){
+      Double_t delta=(point[ipar]-fBinCenter[ipar]);
+      d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
+    }
+  }
+  return kTRUE;
+}
+
+Bool_t AliNDLocalRegression::EvalAndDerivative(Double_t *point, Double_t &val, Double_t *d)
+{
+  // fill d by partial derivatives and calculate value in one go
+  //
+  const Double_t almost0=0.00000001;
+  for (Int_t iDim=0; iDim<fNParameters; iDim++){
+    const TAxis* ax = fHistPoints->GetAxis(iDim);
+    if      (point[iDim]<=ax->GetXmin()) point[iDim] = ax->GetXmin()+almost0*ax->GetBinWidth(0);
+    else if (point[iDim]>=ax->GetXmax()) point[iDim] = ax->GetXmax()-almost0*ax->GetBinWidth(0);
+  }
+  val = 0;
+  Int_t ibin = fHistPoints->GetBin(point);
+  if (ibin>=fLocalFitParam->GetEntriesFast() ||
+      !fLocalFitParam->UncheckedAt(ibin) ) return kFALSE;
+  //
+  fHistPoints->GetBinContent(ibin,fBinIndex); 
+  for (Int_t idim=0; idim<fNParameters; idim++){
+    const TAxis* ax = fHistPoints->GetAxis(idim);
+    fBinCenter[idim] = ax->GetBinCenter(fBinIndex[idim]);
+    fBinWidth[idim] = ax->GetBinWidth(fBinIndex[idim]);
+  } 
+  TVectorD &vecParam = *((TVectorD*)fLocalFitParam->At(ibin));
+
+  val = vecParam[0];
+  if (fUseBinNorm){
+    for (Int_t ipar=0; ipar<fNParameters; ipar++){
+      Double_t delta=(point[ipar]-fBinCenter[ipar])/fBinWidth[ipar];
+      d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
+      val += (vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
+    }
+  } else {
+    for (Int_t ipar=0; ipar<fNParameters; ipar++){
+      Double_t delta=(point[ipar]-fBinCenter[ipar]);
+      d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
+      val += (vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
+    }
+  }
+  return kTRUE;
+}
+
 
 Int_t  AliNDLocalRegression::GetVisualCorrectionIndex(const char *corName){
   //

diff --git a/STAT/AliNDLocalRegression.h b/STAT/AliNDLocalRegression.h
@@ -33,6 +33,8 @@ class AliNDLocalRegression : public TNamed {
 
   Double_t Eval(Double_t *point);
   Double_t EvalError(Double_t *point);
+  Bool_t   Derivative(Double_t *point, Double_t *d);
+  Bool_t   EvalAndDerivative(Double_t *point, Double_t &val, Double_t *d);
   const THn *GetHistogram() {return fHistPoints;}
   const TObjArray *   GetFitParam(){ return fLocalFitParam;}
   void SetCuts(Double_t nSigma=6, Double_t robustFraction=0.95, Int_t estimator=1);