Adding function neccessary for multidimentsional linear fits
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 20 Feb 2011 20:36:20 +0000 (20:36 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 20 Feb 2011 20:36:20 +0000 (20:36 +0000)
with constraints.
Marian

STAT/TStatToolkit.cxx
STAT/TStatToolkit.h

index 727993b..24493ea 100644 (file)
@@ -1021,3 +1021,150 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
    delete[] values;
    return preturnFormula;
 }
+
+
+
+
+
+Int_t TStatToolkit::GetFitIndex(TString fString, TString subString){
+  //
+  // fitString - ++ separated list of fits
+  // substring - ++ separated list of the requiered substrings
+  //
+  // return the last occurance of substring in fit string
+  // 
+  TObjArray *arrFit = fString.Tokenize("++");
+  TObjArray *arrSub = subString.Tokenize("++");
+  Int_t index=-1;
+  for (Int_t i=0; i<arrFit->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str =arrFit->At(i)->GetName();
+    for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
+      if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
+    }
+    if (isOK) index=i;
+  }
+  return index;
+}
+
+
+TString  TStatToolkit::FilterFit(TString &input, TString filter, TVectorD &param, TMatrixD & covar){
+  //
+  // Filter fit expression make sub-fit
+  //
+  TObjArray *array0= input.Tokenize("++");
+  TObjArray *array1= filter.Tokenize("++");
+  //TString *presult=new TString("(0");
+  TString result="(0.0";
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str(array0->At(i)->GetName());
+    for (Int_t j=0; j<array1->GetEntries(); j++){
+      if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;      
+    }
+    if (isOK) {
+      result+="+"+str;
+      result+=Form("*(%f)",param[i+1]);
+      printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());    
+    }
+  }
+  result+="-0.)";
+  return result;
+}
+
+void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
+  //
+  // Update parameters and covariance - with one measurement
+  // Input:
+  // vecXk - input vector - Updated in function 
+  // covXk - covariance matrix - Updated in function
+  // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
+  const Int_t knMeas=1;
+  Int_t knElem=vecXk.GetNrows();
+  TMatrixD mat1(knElem,knElem);            // update covariance matrix
+  TMatrixD matHk(1,knElem);        // vector to mesurement
+  TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
+  TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
+  TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
+  TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
+  TMatrixD covXk2(knElem,knElem);  // helper matrix
+  TMatrixD covXk3(knElem,knElem);  // helper matrix
+  TMatrixD vecZk(1,1);
+  TMatrixD measR(1,1);
+  vecZk(0,0)=delta;
+  measR(0,0)=sigma*sigma;
+  //
+  // reset matHk
+  for (Int_t iel=0;iel<knElem;iel++) 
+    for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
+  //mat1
+  for (Int_t iel=0;iel<knElem;iel++) {
+    for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
+    mat1(iel,iel)=1;
+  }
+  //
+  matHk(0, s1)=1;
+  vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
+  matHkT=matHk.T(); matHk.T();
+  matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
+  matSk.Invert();
+  matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
+  vecXk += matKk*vecYk;                    //  updated vector 
+  covXk2= (mat1-(matKk*matHk));
+  covXk3 =  covXk2*covXk;          
+  covXk = covXk3;  
+  Int_t nrows=covXk3.GetNrows();
+  
+  for (Int_t irow=0; irow<nrows; irow++)
+    for (Int_t icol=0; icol<nrows; icol++){
+      // rounding problems - make matrix again symteric
+      covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; 
+    }
+}
+
+
+
+void   TStatToolkit::Constrain1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
+  //
+  // constrain linear fit
+  // input  - string description of fit function
+  // filter - string filter to select sub fits
+  // param,covar - parameters and covariance matrix of the fit
+  // mean,sigma  - new measurement uning which the fit is updated
+  //
+  TObjArray *array0= input.Tokenize("++");
+  TObjArray *array1= filter.Tokenize("++");
+  TMatrixD paramM(param.GetNrows(),1);
+  for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
+  
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str(array0->At(i)->GetName());
+    for (Int_t j=0; j<array1->GetEntries(); j++){
+      if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;      
+    }
+    if (isOK) {
+      TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
+    }
+  }
+  for (Int_t i=0; i<=array0->GetEntries(); i++){
+    param(i)=paramM(i,0);
+  }
+}
+
+TString  TStatToolkit::MakeFitString(TString &input, TVectorD &param, TMatrixD & covar){
+  //
+  //
+  //
+  TObjArray *array0= input.Tokenize("++");
+  TString result="(0.0";
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    TString str(array0->At(i)->GetName());
+    result+="+"+str;
+    result+=Form("*(%f)",param[i+1]);
+    printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());    
+  }
+  result+="-0.)";
+  return result;
+}
index e940ac0..79aedce 100644 (file)
@@ -44,16 +44,21 @@ class TStatToolkit : public TObject
   //
   static TString* FitPlane(TTree * tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints,  TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac=-1, Int_t start=0, Int_t stop=10000000, Bool_t fix0=kFALSE);
   static TString* FitPlaneFixed(TTree * tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints,  TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac=-1, Int_t start=0, Int_t stop=10000000);
-
-
+  //
+  //Linear fitter helper function
+  //
   static TString* FitPlaneConstrain(TTree * tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints,  TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac=-1, Int_t start=0, Int_t stop=10000000, Double_t constrain=-1);
-
+  static Int_t GetFitIndex(TString fString, TString subString);
+ static TString FilterFit(TString &input, TString filter, TVectorD &vec, TMatrixD &covar);
+ static void Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &param, TMatrixD &covar);
+  static void   Constrain1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma);
+  static TString  MakeFitString(TString &input, TVectorD &param, TMatrixD & covar);
 
   //
   // TestFunctions:
   //
  static  void TestGausFit(Int_t nhistos=5000);
-    
+
  ClassDef(TStatToolkit,0) // Various mathematical tools for physics analysis - which are not included in ROOT TMath
  
 };