#include "TLinearFitter.h"
#include "TDecompChol.h"
#include "TGraph.h"
#include "TGraph2D.h"
#include "TMultiGraph.h"
#include "TRandom.h"
#include "TObjString.h"
#include "TF1.h"
ClassImp(TLinearFitter)
TLinearFitter::TLinearFitter()
{
fChisquare=0;
fNpoints=0;
fY2=0;
fNfixed=0;
fIsSet=kFALSE;
fFormula=0;
fFixedParams=0;
fSpecial=0;
fInputFunction=0;
fStoreData=kTRUE;
fRobust=kFALSE;
}
TLinearFitter::TLinearFitter(Int_t ndim)
{
fNdim=ndim;
fNpoints=0;
fY2=0;
fNfixed=0;
fFixedParams=0;
fFormula=0;
fIsSet=kFALSE;
fChisquare=0;
fSpecial=0;
fInputFunction=0;
fStoreData=kTRUE;
fRobust = kFALSE;
}
TLinearFitter::TLinearFitter(Int_t ndim, const char *formula, Option_t *opt)
{
fNdim=ndim;
fNpoints=0;
fChisquare=0;
fY2=0;
fNfixed=0;
fFixedParams=0;
fSpecial=0;
fInputFunction=0;
TString option=opt;
option.ToUpper();
if (option.Contains("D"))
fStoreData=kTRUE;
else
fStoreData=kFALSE;
fRobust=kFALSE;
SetFormula(formula);
}
TLinearFitter::TLinearFitter(TFormula *function, Option_t *opt)
{
fNdim=function->GetNdim();
if (!function->IsLinear()){
Int_t number=function->GetNumber();
if (number<299 || number>310){
Error("TLinearFitter", "Trying to fit with a nonlinear function");
return;
}
}
fNpoints=0;
fChisquare=0;
fY2=0;
fNfixed=0;
fFixedParams=0;
fSpecial=0;
fFormula = 0;
TString option=opt;
option.ToUpper();
if (option.Contains("D"))
fStoreData=kTRUE;
else
fStoreData=kFALSE;
fIsSet=kTRUE;
fRobust=kFALSE;
SetFormula(function);
}
TLinearFitter::TLinearFitter(const TLinearFitter& tlf) :
TVirtualFitter(tlf),
fParams(tlf.fParams),
fParCovar(tlf.fParCovar),
fTValues(tlf.fTValues),
fParSign(tlf.fParSign),
fDesign(tlf.fDesign),
fDesignTemp(tlf.fDesignTemp),
fDesignTemp2(tlf.fDesignTemp2),
fDesignTemp3(tlf.fDesignTemp3),
fAtb(tlf.fAtb),
fAtbTemp(tlf.fAtbTemp),
fAtbTemp2(tlf.fAtbTemp2),
fAtbTemp3(tlf.fAtbTemp3),
fFunctions(tlf.fFunctions),
fY(tlf.fY),
fY2(tlf.fY2),
fY2Temp(tlf.fY2Temp),
fX(tlf.fX),
fE(tlf.fE),
fInputFunction((TFormula*)tlf.fInputFunction->Clone()),
fNpoints(tlf.fNpoints),
fNfunctions(tlf.fNfunctions),
fFormulaSize(tlf.fFormulaSize),
fNdim(tlf.fNdim),
fNfixed(tlf.fNfixed),
fSpecial(tlf.fSpecial),
fIsSet(tlf.fIsSet),
fStoreData(tlf.fStoreData),
fChisquare(tlf.fChisquare),
fH(tlf.fH),
fRobust(tlf.fRobust),
fFitsample(tlf.fFitsample)
{
fFixedParams=new Bool_t[fNfixed];
for(Int_t i=0; i<fNfixed; ++i)
fFixedParams[i]=tlf.fFixedParams[i];
strcpy(fFormula,tlf.fFormula);
}
TLinearFitter::~TLinearFitter()
{
if (fFormula)
delete [] fFormula;
fFormula = 0;
if (fFixedParams) delete [] fFixedParams;
fFixedParams = 0;
fInputFunction = 0;
fFunctions.Delete();
}
TLinearFitter& TLinearFitter::operator=(const TLinearFitter& tlf)
{
if(this!=&tlf) {
TVirtualFitter::operator=(tlf);
fParams=tlf.fParams;
fParCovar=tlf.fParCovar;
fTValues=tlf.fTValues;
fParSign=tlf.fParSign;
fDesign=tlf.fDesign;
fDesignTemp=tlf.fDesignTemp;
fDesignTemp2=tlf.fDesignTemp2;
fDesignTemp3=tlf.fDesignTemp3;
fAtb=tlf.fAtb;
fAtbTemp=tlf.fAtbTemp;
fAtbTemp2=tlf.fAtbTemp2;
fAtbTemp3=tlf.fAtbTemp3;
fFixedParams=new Bool_t[tlf.fNfixed];
for(Int_t i=0; i<tlf.fNfixed; ++i)
fFixedParams[i]=tlf.fFixedParams[i];
fFunctions=tlf.fFunctions;
fY=tlf.fY;
fY2=tlf.fY2;
fY2Temp=tlf.fY2Temp;
fX=tlf.fX;
fE=tlf.fE;
fInputFunction=(TFormula*)tlf.fInputFunction->Clone();
fNpoints=tlf.fNpoints;
fNfunctions=tlf.fNfunctions;
fFormulaSize=tlf.fFormulaSize;
fNdim=tlf.fNdim;
fNfixed=tlf.fNfixed;
fSpecial=tlf.fSpecial;
strcpy(fFormula,tlf.fFormula);
fIsSet=tlf.fIsSet;
fStoreData=tlf.fStoreData;
fChisquare=tlf.fChisquare;
fH=tlf.fH;
fRobust=tlf.fRobust;
fFitsample=tlf.fFitsample;
}
return *this;
}
void TLinearFitter::AddPoint(Double_t *x, Double_t y, Double_t e)
{
Int_t size;
fNpoints++;
if (fStoreData){
size=fY.GetNoElements();
if (size<fNpoints){
fY.ResizeTo(fNpoints+fNpoints/2);
fE.ResizeTo(fNpoints+fNpoints/2);
fX.ResizeTo(fNpoints+fNpoints/2, fNdim);
}
Int_t j=fNpoints-1;
fY(j)=y;
fE(j)=e;
for (Int_t i=0; i<fNdim; i++)
fX(j,i)=x[i];
}
if (!fFunctions.IsEmpty() || fInputFunction || fSpecial>199 || !fRobust)
AddToDesign(x, y, e);
else if (!fStoreData)
Error("AddPoint", "Point can't be added, because the formula hasn't been set and data is not stored");
}
void TLinearFitter::AssignData(Int_t npoints, Int_t xncols, Double_t *x, Double_t *y, Double_t *e)
{
if (npoints<fNpoints){
Error("AddData", "Those points are already added");
return;
}
Bool_t same=kFALSE;
if (fX.GetMatrixArray()==x && fY.GetMatrixArray()==y){
if (e){
if (fE.GetMatrixArray()==e)
same=kTRUE;
}
}
fX.Use(npoints, xncols, x);
fY.Use(npoints, y);
if (e)
fE.Use(npoints, e);
else {
fE.ResizeTo(npoints);
fE=1;
}
Int_t xfirst;
if (!fFunctions.IsEmpty() || fInputFunction || fSpecial>199) {
if (same)
xfirst=fNpoints;
else
xfirst=0;
for (Int_t i=xfirst; i<npoints; i++)
AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
}
fNpoints=npoints;
}
void TLinearFitter::AddToDesign(Double_t *x, Double_t y, Double_t e)
{
Int_t i, j, ii;
y/=e;
Double_t val[100];
if ((fSpecial>100)&&(fSpecial<200)){
Int_t npar=fSpecial-100;
val[0]=1;
for (i=1; i<npar; i++)
val[i]=val[i-1]*x[0];
for (i=0; i<npar; i++)
val[i]/=e;
} else {
if (fSpecial>200){
Int_t npar=fSpecial-201;
val[0]=1./e;
for (i=0; i<npar; i++)
val[i+1]=x[i]/e;
} else {
for (ii=0; ii<fNfunctions; ii++){
if (!fFunctions.IsEmpty()){
TFormula *f1 = (TFormula*)(fFunctions.UncheckedAt(ii));
val[ii]=f1->EvalPar(x)/e;
} else {
TFormula *f=(TFormula*)fInputFunction->GetLinearPart(ii);
val[ii]=f->EvalPar(x)/e;
}
}
}
}
for (i=0; i<fNfunctions; i++){
for (j=0; j<i; j++)
fDesignTemp3(j, i)+=val[i]*val[j];
fDesignTemp3(i, i)+=val[i]*val[i];
fAtbTemp3(i)+=val[i]*y;
}
fY2Temp+=y*y;
fIsSet=kTRUE;
if (fNpoints % 100 == 0 && fNpoints>100){
fDesignTemp2+=fDesignTemp3;
fDesignTemp3.Zero();
fAtbTemp2+=fAtbTemp3;
fAtbTemp3.Zero();
if (fNpoints % 10000 == 0 && fNpoints>10000){
fDesignTemp+=fDesignTemp2;
fDesignTemp2.Zero();
fAtbTemp+=fAtbTemp2;
fAtbTemp2.Zero();
fY2+=fY2Temp;
fY2Temp=0;
if (fNpoints % 1000000 == 0 && fNpoints>1000000){
fDesign+=fDesignTemp;
fDesignTemp.Zero();
fAtb+=fAtbTemp;
fAtbTemp.Zero();
}
}
}
}
void TLinearFitter::Clear(Option_t * )
{
fParams.Clear();
fParCovar.Clear();
fTValues.Clear();
fParSign.Clear();
fDesign.Clear();
fDesignTemp.Clear();
fDesignTemp2.Clear();
fDesignTemp3.Clear();
fAtb.Clear();
fAtbTemp.Clear();
fAtbTemp2.Clear();
fAtbTemp3.Clear();
fFunctions.Clear();
fInputFunction=0;
fY.Clear();
fX.Clear();
fE.Clear();
fNpoints=0;
fNfunctions=0;
fFormulaSize=0;
fNdim=0;
delete [] fFormula;
fFormula=0;
fIsSet=0;
delete [] fFixedParams;
fFixedParams=0;
fChisquare=0;
fY2=0;
fSpecial=0;
fRobust=kFALSE;
fFitsample.Clear();
}
void TLinearFitter::ClearPoints()
{
fDesign.Zero();
fAtb.Zero();
fDesignTemp.Zero();
fDesignTemp2.Zero();
fDesignTemp3.Zero();
fAtbTemp.Zero();
fAtbTemp2.Zero();
fAtbTemp3.Zero();
fParams.Zero();
fParCovar.Zero();
fTValues.Zero();
fParSign.Zero();
for (Int_t i=0; i<fNfunctions; i++)
fFixedParams[i]=0;
fChisquare=0;
fNpoints=0;
}
void TLinearFitter::Chisquare()
{
Int_t i, j;
Double_t sumtotal2;
Double_t temp, temp2;
if (!fStoreData){
sumtotal2 = 0;
for (i=0; i<fNfunctions; i++){
for (j=0; j<i; j++){
sumtotal2 += 2*fParams(i)*fParams(j)*fDesign(j, i);
}
sumtotal2 += fParams(i)*fParams(i)*fDesign(i, i);
sumtotal2 -= 2*fParams(i)*fAtb(i);
}
sumtotal2 += fY2;
} else {
sumtotal2 = 0;
if (fInputFunction){
for (i=0; i<fNpoints; i++){
temp = fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
temp2 = (fY(i)-temp)*(fY(i)-temp);
temp2 /= fE(i)*fE(i);
sumtotal2 += temp2;
}
} else {
sumtotal2 = 0;
Double_t val[100];
for (Int_t point=0; point<fNpoints; point++){
temp = 0;
if ((fSpecial>100)&&(fSpecial<200)){
Int_t npar = fSpecial-100;
val[0] = 1;
for (i=1; i<npar; i++)
val[i] = val[i-1]*fX(point, 0);
for (i=0; i<npar; i++)
temp += fParams(i)*val[i];
} else {
if (fSpecial>200) {
Int_t npar = fSpecial-201;
temp+=fParams(0);
for (i=0; i<npar; i++)
temp += fParams(i+1)*fX(point, i);
} else {
for (j=0; j<fNfunctions; j++) {
TFormula *f1 = (TFormula*)(fFunctions.UncheckedAt(j));
val[j] = f1->EvalPar(TMatrixDRow(fX, point).GetPtr());
temp += fParams(j)*val[j];
}
}
}
temp2 = (fY(point)-temp)*(fY(point)-temp);
temp2 /= fE(point)*fE(point);
sumtotal2 += temp2;
}
}
}
fChisquare = sumtotal2;
}
void TLinearFitter::ComputeTValues()
{
for (Int_t i=0; i<fNfunctions; i++){
fTValues(i) = fParams(i)/(TMath::Sqrt(fParCovar(i, i)));
fParSign(i) = 2*(1-TMath::StudentI(TMath::Abs(fTValues(i)),fNpoints-fNfunctions+fNfixed));
}
}
Int_t TLinearFitter::Eval()
{
Double_t e;
if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<200)){
Error("TLinearFitter::Eval", "The formula hasn't been set");
return 1;
}
fParams.ResizeTo(fNfunctions);
fTValues.ResizeTo(fNfunctions);
fParSign.ResizeTo(fNfunctions);
fParCovar.ResizeTo(fNfunctions,fNfunctions);
fChisquare=0;
if (!fIsSet){
Bool_t update = UpdateMatrix();
if (!update){
fParams.Zero();
fParCovar.Zero();
fTValues.Zero();
fParSign.Zero();
fChisquare=0;
if (fInputFunction){
((TF1*)fInputFunction)->SetParameters(fParams.GetMatrixArray());
for (Int_t i=0; i<fNfunctions; i++)
((TF1*)fInputFunction)->SetParError(i, 0);
((TF1*)fInputFunction)->SetChisquare(0);
((TF1*)fInputFunction)->SetNDF(0);
((TF1*)fInputFunction)->SetNumberFitPoints(0);
}
return 1;
}
}
fDesignTemp2+=fDesignTemp3;
fDesignTemp+=fDesignTemp2;
fDesign+=fDesignTemp;
fDesignTemp3.Zero();
fDesignTemp2.Zero();
fDesignTemp.Zero();
fAtbTemp2+=fAtbTemp3;
fAtbTemp+=fAtbTemp2;
fAtb+=fAtbTemp;
fAtbTemp3.Zero();
fAtbTemp2.Zero();
fAtbTemp.Zero();
fY2+=fY2Temp;
fY2Temp=0;
Int_t i, ii, j=0;
if (fNfixed>0){
for (ii=0; ii<fNfunctions; ii++)
fDesignTemp(ii, fNfixed) = fAtb(ii);
for (i=0; i<fNfunctions; i++){
if (fFixedParams[i]){
for (ii=0; ii<i; ii++)
fDesignTemp(ii, j) = fDesign(ii, i);
for (ii=i; ii<fNfunctions; ii++)
fDesignTemp(ii, j) = fDesign(i, ii);
j++;
for (ii=0; ii<fNfunctions; ii++){
fAtb(ii)-=fParams(i)*(fDesignTemp(ii, j-1));
}
}
}
for (i=0; i<fNfunctions; i++){
if (fFixedParams[i]){
for (ii=0; ii<fNfunctions; ii++){
fDesign(ii, i) = 0;
fDesign(i, ii) = 0;
}
fDesign (i, i) = 1;
fAtb(i) = fParams(i);
}
}
}
TDecompChol chol(fDesign);
Bool_t ok;
TVectorD coef(fNfunctions);
coef=chol.Solve(fAtb, ok);
if (!ok){
Error("Eval","Matrix inversion failed");
fParams.Zero();
fParCovar.Zero();
fTValues.Zero();
fParSign.Zero();
return 1;
}
fParams=coef;
fParCovar=chol.Invert();
if (fInputFunction){
fInputFunction->SetParameters(fParams.GetMatrixArray());
for (i=0; i<fNfunctions; i++){
e = TMath::Sqrt(fParCovar(i, i));
((TF1*)fInputFunction)->SetParError(i, e);
}
if (!fObjectFit)
((TF1*)fInputFunction)->SetChisquare(GetChisquare());
((TF1*)fInputFunction)->SetNDF(fNpoints-fNfunctions+fNfixed);
((TF1*)fInputFunction)->SetNumberFitPoints(fNpoints);
}
j = 0;
if (fNfixed>0){
for (i=0; i<fNfunctions; i++){
if (fFixedParams[i]){
for (ii=0; ii<i; ii++){
fDesign(ii, i) = fDesignTemp(ii, j);
fAtb(ii) = fDesignTemp(ii, fNfixed);
}
for (ii=i; ii<fNfunctions; ii++){
fDesign(i, ii) = fDesignTemp(ii, j);
fAtb(ii) = fDesignTemp(ii, fNfixed);
}
j++;
}
}
}
return 0;
}
void TLinearFitter::FixParameter(Int_t ipar)
{
if (fParams.NonZeros()<1){
Error("FixParameter", "no value available to fix the parameter");
return;
}
if (ipar>fNfunctions || ipar<0){
Error("FixParameter", "illegal parameter value");
return;
}
if (fNfixed==fNfunctions) {
Error("FixParameter", "no free parameters left");
return;
}
if (!fFixedParams)
fFixedParams = new Bool_t[fNfunctions];
fFixedParams[ipar] = 1;
fNfixed++;
}
void TLinearFitter::FixParameter(Int_t ipar, Double_t parvalue)
{
if (ipar>fNfunctions || ipar<0){
Error("FixParameter", "illegal parameter value");
return;
}
if (fNfixed==fNfunctions) {
Error("FixParameter", "no free parameters left");
return;
}
if(!fFixedParams)
fFixedParams = new Bool_t[fNfunctions];
fFixedParams[ipar] = 1;
if (fParams.GetNoElements()<fNfunctions)
fParams.ResizeTo(fNfunctions);
fParams(ipar) = parvalue;
fNfixed++;
}
void TLinearFitter::ReleaseParameter(Int_t ipar)
{
if (ipar>fNfunctions || ipar<0){
Error("ReleaseParameter", "illegal parameter value");
return;
}
if (!fFixedParams[ipar]){
Warning("ReleaseParameter","This parameter is not fixed\n");
return;
} else {
fFixedParams[ipar] = 0;
fNfixed--;
}
}
Double_t TLinearFitter::GetChisquare()
{
if (fChisquare > 1e-16)
return fChisquare;
else {
Chisquare();
return fChisquare;
}
}
void TLinearFitter::GetConfidenceIntervals(Int_t n, Int_t ndim, const Double_t *x, Double_t *ci, Double_t cl)
{
if (fInputFunction){
Double_t *grad = new Double_t[fNfunctions];
Double_t *sum_vector = new Double_t[fNfunctions];
Double_t c=0;
Int_t df = fNpoints-fNfunctions+fNfixed;
Double_t t = TMath::StudentQuantile(cl, df);
Double_t chidf = TMath::Sqrt(fChisquare/df);
for (Int_t ipoint=0; ipoint<n; ipoint++){
c=0;
((TF1*)(fInputFunction))->GradientPar(x+ndim*ipoint, grad);
for (Int_t irow=0; irow<fNfunctions; irow++){
sum_vector[irow]=0;
for (Int_t icol=0; icol<fNfunctions; icol++){
sum_vector[irow]+=fParCovar(irow,icol)*grad[icol];
}
}
for (Int_t i=0; i<fNfunctions; i++)
c+=grad[i]*sum_vector[i];
c=TMath::Sqrt(c);
ci[ipoint]=c*t*chidf;
}
delete [] grad;
delete [] sum_vector;
}
}
void TLinearFitter::GetConfidenceIntervals(TObject *obj, Double_t cl)
{
if (!fInputFunction) {
Error("GetConfidenceIntervals", "The case of fitting not with a TFormula is not yet implemented");
return;
}
if (obj->InheritsFrom(TGraph::Class())) {
TGraph *gr = (TGraph*)obj;
if (!gr->GetEY()){
Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a graph");
return;
}
if (fObjectFit->InheritsFrom(TGraph2D::Class())){
Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a graph");
return;
}
if (fObjectFit->InheritsFrom(TH1::Class())){
if (((TH1*)(fObjectFit))->GetDimension()>1){
Error("GetConfidenceIntervals", "A TGraph2DErrors or a TH23 should be passed instead of a graph");
return;
}
}
GetConfidenceIntervals(gr->GetN(),1,gr->GetX(), gr->GetEY(), cl);
for (Int_t i=0; i<gr->GetN(); i++)
gr->SetPoint(i, gr->GetX()[i], fInputFunction->Eval(gr->GetX()[i]));
}
else if (obj->InheritsFrom(TGraph2D::Class())) {
TGraph2D *gr2 = (TGraph2D*)obj;
if (!gr2->GetEZ()){
Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a TGraph2D");
return;
}
if (fObjectFit->InheritsFrom(TGraph::Class())){
Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a TGraph2D");
return;
}
if (fObjectFit->InheritsFrom(TH1::Class())){
if (((TH1*)(fObjectFit))->GetDimension()==1){
Error("GetConfidenceIntervals", "A TGraphErrors or a TH1 should be passed instead of a graph");
return;
}
}
Double_t xy[2];
Int_t np = gr2->GetN();
Double_t *grad = new Double_t[fNfunctions];
Double_t *sum_vector = new Double_t[fNfunctions];
Double_t *x = gr2->GetX();
Double_t *y = gr2->GetY();
Double_t t = TMath::StudentQuantile(cl, ((TF1*)(fInputFunction))->GetNDF());
Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
Double_t c = 0;
for (Int_t ipoint=0; ipoint<np; ipoint++){
c=0;
xy[0]=x[ipoint];
xy[1]=y[ipoint];
((TF1*)(fInputFunction))->GradientPar(xy, grad);
for (Int_t irow=0; irow<fNfunctions; irow++){
sum_vector[irow]=0;
for (Int_t icol=0; icol<fNfunctions; icol++)
sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
}
for (Int_t i=0; i<fNfunctions; i++)
c+=grad[i]*sum_vector[i];
c=TMath::Sqrt(c);
gr2->SetPoint(ipoint, xy[0], xy[1], fInputFunction->EvalPar(xy));
gr2->GetEZ()[ipoint]=c*t*chidf;
}
delete [] grad;
delete [] sum_vector;
}
else if (obj->InheritsFrom(TH1::Class())) {
if (fObjectFit->InheritsFrom(TGraph::Class())){
if (((TH1*)obj)->GetDimension()>1){
Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
return;
}
}
if (fObjectFit->InheritsFrom(TGraph2D::Class())){
if (((TH1*)obj)->GetDimension()!=2){
Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
return;
}
}
if (fObjectFit->InheritsFrom(TH1::Class())){
if (((TH1*)(fObjectFit))->GetDimension()!=((TH1*)(obj))->GetDimension()){
Error("GetConfidenceIntervals", "Fitted and passed histograms have different number of dimensions");
return;
}
}
TH1 *hfit = (TH1*)obj;
Double_t *grad = new Double_t[fNfunctions];
Double_t *sum_vector = new Double_t[fNfunctions];
Double_t x[3];
Int_t hxfirst = hfit->GetXaxis()->GetFirst();
Int_t hxlast = hfit->GetXaxis()->GetLast();
Int_t hyfirst = hfit->GetYaxis()->GetFirst();
Int_t hylast = hfit->GetYaxis()->GetLast();
Int_t hzfirst = hfit->GetZaxis()->GetFirst();
Int_t hzlast = hfit->GetZaxis()->GetLast();
TAxis *xaxis = hfit->GetXaxis();
TAxis *yaxis = hfit->GetYaxis();
TAxis *zaxis = hfit->GetZaxis();
Double_t t = TMath::StudentQuantile(cl, ((TF1*)(fInputFunction))->GetNDF());
Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
Double_t c=0;
for (Int_t binz=hzfirst; binz<=hzlast; binz++){
x[2]=zaxis->GetBinCenter(binz);
for (Int_t biny=hyfirst; biny<=hylast; biny++) {
x[1]=yaxis->GetBinCenter(biny);
for (Int_t binx=hxfirst; binx<=hxlast; binx++) {
x[0]=xaxis->GetBinCenter(binx);
((TF1*)(fInputFunction))->GradientPar(x, grad);
c=0;
for (Int_t irow=0; irow<fNfunctions; irow++){
sum_vector[irow]=0;
for (Int_t icol=0; icol<fNfunctions; icol++)
sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
}
for (Int_t i=0; i<fNfunctions; i++)
c+=grad[i]*sum_vector[i];
c=TMath::Sqrt(c);
hfit->SetBinContent(binx, biny, binz, fInputFunction->EvalPar(x));
hfit->SetBinError(binx, biny, binz, c*t*chidf);
}
}
}
delete [] grad;
delete [] sum_vector;
}
else {
Error("GetConfidenceIntervals", "This object type is not supported");
return;
}
}
Double_t* TLinearFitter::GetCovarianceMatrix() const
{
Double_t *p = const_cast<Double_t*>(fParCovar.GetMatrixArray());
return p;
}
void TLinearFitter::GetCovarianceMatrix(TMatrixD &matr)
{
if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
matr.ResizeTo(fNfunctions, fNfunctions);
}
matr = fParCovar;
}
void TLinearFitter::GetErrors(TVectorD &vpar)
{
if (vpar.GetNoElements()!=fNfunctions) {
vpar.ResizeTo(fNfunctions);
}
for (Int_t i=0; i<fNfunctions; i++)
vpar(i) = TMath::Sqrt(fParCovar(i, i));
}
void TLinearFitter::GetParameters(TVectorD &vpar)
{
if (vpar.GetNoElements()!=fNfunctions) {
vpar.ResizeTo(fNfunctions);
}
vpar=fParams;
}
Int_t TLinearFitter::GetParameter(Int_t ipar,char* name,Double_t& value,Double_t& ,Double_t& , Double_t& ) const
{
if (ipar<0 || ipar>fNfunctions) {
Error("GetParError", "illegal value of parameter");
return 0;
}
value = fParams(ipar);
if (fInputFunction)
strcpy(name, fInputFunction->GetParName(ipar));
else
name = "";
return 1;
}
Double_t TLinearFitter::GetParError(Int_t ipar) const
{
if (ipar<0 || ipar>fNfunctions) {
Error("GetParError", "illegal value of parameter");
return 0;
}
return TMath::Sqrt(fParCovar(ipar, ipar));
}
const char *TLinearFitter::GetParName(Int_t ipar) const
{
if (ipar<0 || ipar>fNfunctions) {
Error("GetParError", "illegal value of parameter");
return 0;
}
if (fInputFunction)
return fInputFunction->GetParName(ipar);
return "";
}
Double_t TLinearFitter::GetParTValue(Int_t ipar)
{
if (ipar<0 || ipar>fNfunctions) {
Error("GetParTValue", "illegal value of parameter");
return 0;
}
if (!fTValues.NonZeros())
ComputeTValues();
return fTValues(ipar);
}
Double_t TLinearFitter::GetParSignificance(Int_t ipar)
{
if (ipar<0 || ipar>fNfunctions) {
Error("GetParSignificance", "illegal value of parameter");
return 0;
}
if (!fParSign.NonZeros())
ComputeTValues();
return fParSign(ipar);
}
void TLinearFitter::GetFitSample(TBits &bits)
{
if (!fRobust){
Error("GetFitSample", "there is no fit sample in ordinary least-squares fit");
return;
}
for (Int_t i=0; i<fNpoints; i++)
bits.SetBitNumber(i, fFitsample.TestBitNumber(i));
}
void TLinearFitter::SetDim(Int_t ndim)
{
fNdim=ndim;
fY.ResizeTo(ndim+1);
fX.ResizeTo(ndim+1, ndim);
fE.ResizeTo(ndim+1);
fNpoints=0;
fIsSet=kFALSE;
}
void TLinearFitter::SetFormula(const char *formula)
{
Int_t size, special = 0;
Int_t i;
Bool_t isHyper = kFALSE;
if (fInputFunction)
fInputFunction = 0;
fFormulaSize = strlen(formula);
fFormula = new char[fFormulaSize+1];
strcpy(fFormula, formula);
fSpecial = 0;
char *fstring;
fstring = (char *)strstr(fFormula, "hyp");
if (fstring!=0){
isHyper = kTRUE;
fstring+=3;
sscanf(fstring, "%d", &size);
size++;
fSpecial=200+size;
}
if (fSpecial==0) {
TString sstring(fFormula);
sstring = sstring.ReplaceAll("++", 2, "|", 1);
TString replaceformula;
TObjArray *oa = sstring.Tokenize("|");
if (!fFunctions.IsEmpty())
fFunctions.Clear();
fNfunctions = oa->GetEntriesFast();
fFunctions.Expand(fNfunctions);
char pattern[5];
char replacement[6];
for (i=0; i<fNdim; i++){
sprintf(pattern, "x%d", i);
sprintf(replacement, "x[%d]", i);
sstring = sstring.ReplaceAll(pattern, Int_t(i/10)+2, replacement, Int_t(i/10)+4);
}
oa = sstring.Tokenize("|");
for (Int_t i=0; i<fNfunctions; i++) {
replaceformula = ((TObjString *)oa->UncheckedAt(i))->GetString();
TFormula *f = new TFormula("f", replaceformula.Data());
if (!f) {
Error("TLinearFitter", "f_linear not allocated");
return;
}
special=f->GetNumber();
fFunctions.Add(f);
}
if ((fNfunctions==1)&&(special>299)&&(special<310)){
size=special-299;
fSpecial=100+size;
} else
size=fNfunctions;
oa->Delete();
delete oa;
}
fNfunctions=size;
fDesign.ResizeTo(size, size);
fAtb.ResizeTo(size);
fDesignTemp.ResizeTo(size, size);
fDesignTemp2.ResizeTo(size, size);
fDesignTemp3.ResizeTo(size, size);
fAtbTemp.ResizeTo(size);
fAtbTemp2.ResizeTo(size);
fAtbTemp3.ResizeTo(size);
if (fFixedParams)
delete [] fFixedParams;
fFixedParams=new Bool_t[size];
fDesign.Zero();
fAtb.Zero();
fDesignTemp.Zero();
fDesignTemp2.Zero();
fDesignTemp3.Zero();
fAtbTemp.Zero();
fAtbTemp2.Zero();
fAtbTemp3.Zero();
fY2Temp=0;
fY2=0;
for (i=0; i<size; i++)
fFixedParams[i]=0;
fIsSet=kFALSE;
fChisquare=0;
}
void TLinearFitter::SetFormula(TFormula *function)
{
Int_t special, size;
fInputFunction=function;
fNfunctions=fInputFunction->GetNpar();
fSpecial = 0;
special=fInputFunction->GetNumber();
if (!fFunctions.IsEmpty())
fFunctions.Delete();
if ((special>299)&&(special<310)){
size=special-299;
fSpecial=100+size;
} else
size=fNfunctions;
fNfunctions=size;
fDesign.ResizeTo(size, size);
fAtb.ResizeTo(size);
fDesignTemp.ResizeTo(size, size);
fAtbTemp.ResizeTo(size);
fDesignTemp2.ResizeTo(size, size);
fDesignTemp3.ResizeTo(size, size);
fAtbTemp2.ResizeTo(size);
fAtbTemp3.ResizeTo(size);
if (fFixedParams)
delete [] fFixedParams;
fFixedParams=new Bool_t[size];
fDesign.Zero();
fAtb.Zero();
fDesignTemp.Zero();
fAtbTemp.Zero();
fDesignTemp2.Zero();
fDesignTemp3.Zero();
fAtbTemp2.Zero();
fAtbTemp3.Zero();
fY2Temp=0;
fY2=0;
for (Int_t i=0; i<size; i++)
fFixedParams[i]=0;
Double_t al,bl;
for (Int_t i=0;i<fNfunctions;i++) {
((TF1*)function)->GetParLimits(i,al,bl);
if (al*bl != 0 && al >= bl) {
FixParameter(i, function->GetParameter(i));
}
}
fIsSet=kFALSE;
fChisquare=0;
}
Bool_t TLinearFitter::UpdateMatrix()
{
if (fStoreData) {
for (Int_t i=0; i<fNpoints; i++) {
AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
}
return 1;
} else
return 0;
}
Int_t TLinearFitter::ExecuteCommand(const char *command, Double_t *args, Int_t )
{
if (!strcmp(command, "FitGraph")){
if (args) return GraphLinearFitter(args[0]);
else return GraphLinearFitter(0);
}
if (!strcmp(command, "FitGraph2D")){
if (args) return Graph2DLinearFitter(args[0]);
else return Graph2DLinearFitter(0);
}
if (!strcmp(command, "FitMultiGraph")){
if (args) return MultiGraphLinearFitter(args[0]);
else return MultiGraphLinearFitter(0);
}
if (!strcmp(command, "FitHist")) return HistLinearFitter();
return 0;
}
void TLinearFitter::PrintResults(Int_t level, Double_t ) const
{
if (level==3){
if (!fRobust){
printf("Fitting results:\nParameters:\nNO.\t\tVALUE\t\tERROR\n");
for (Int_t i=0; i<fNfunctions; i++){
printf("%d\t%e\t%e\n", i, fParams(i), TMath::Sqrt(fParCovar(i, i)));
}
} else {
printf("Fitting results:\nParameters:\nNO.\t\tVALUE\n");
for (Int_t i=0; i<fNfunctions; i++){
printf("%d\t%e\n", i, fParams(i));
}
}
}
}
Int_t TLinearFitter::GraphLinearFitter(Double_t h)
{
StoreData(kFALSE);
TGraph *grr=(TGraph*)GetObjectFit();
TF1 *f1=(TF1*)GetUserFunc();
Foption_t fitOption=GetFitOption();
Double_t *x=grr->GetX();
Double_t *y=grr->GetY();
Double_t e;
Int_t fitResult = 0;
SetDim(1);
SetFormula(f1);
if (fitOption.Robust){
fRobust=kTRUE;
StoreData(kTRUE);
}
Int_t n=grr->GetN();
for (Int_t i=0; i<n; i++){
if (!f1->IsInside(&x[i])) continue;
e=grr->GetErrorY(i);
if (e<0 || fitOption.W1)
e=1;
AddPoint(&x[i], y[i], e);
}
if (fitOption.Robust){
return EvalRobust(h);
}
fitResult = Eval();
if (!fitOption.Nochisq){
Double_t temp, temp2, sumtotal=0;
for (Int_t i=0; i<n; i++){
if (!f1->IsInside(&x[i])) continue;
temp=f1->Eval(x[i]);
temp2=(y[i]-temp)*(y[i]-temp);
e=grr->GetErrorY(i);
if (e<0 || fitOption.W1)
e=1;
temp2/=(e*e);
sumtotal+=temp2;
}
fChisquare=sumtotal;
f1->SetChisquare(fChisquare);
}
return fitResult;
}
Int_t TLinearFitter::Graph2DLinearFitter(Double_t h)
{
StoreData(kFALSE);
TGraph2D *gr=(TGraph2D*)GetObjectFit();
TF2 *f2=(TF2*)GetUserFunc();
Foption_t fitOption=GetFitOption();
Int_t n = gr->GetN();
Double_t *gx = gr->GetX();
Double_t *gy = gr->GetY();
Double_t *gz = gr->GetZ();
Double_t x[2];
Double_t z, e;
Int_t fitResult=0;
SetDim(2);
SetFormula(f2);
if (fitOption.Robust){
fRobust=kTRUE;
StoreData(kTRUE);
}
for (Int_t bin=0;bin<n;bin++) {
x[0] = gx[bin];
x[1] = gy[bin];
if (!f2->IsInside(x)) {
continue;
}
z = gz[bin];
e=gr->GetErrorZ(bin);
if (e<0 || fitOption.W1)
e=1;
AddPoint(x, z, e);
}
if (fitOption.Robust){
return EvalRobust(h);
}
fitResult = Eval();
if (!fitOption.Nochisq){
Double_t temp, temp2, sumtotal=0;
for (Int_t bin=0; bin<n; bin++){
x[0] = gx[bin];
x[1] = gy[bin];
if (!f2->IsInside(x)) {
continue;
}
z = gz[bin];
temp=f2->Eval(x[0], x[1]);
temp2=(z-temp)*(z-temp);
e=gr->GetErrorZ(bin);
if (e<0 || fitOption.W1)
e=1;
temp2/=(e*e);
sumtotal+=temp2;
}
fChisquare=sumtotal;
f2->SetChisquare(fChisquare);
}
return fitResult;
}
Int_t TLinearFitter::MultiGraphLinearFitter(Double_t h)
{
Int_t n, i;
Double_t *gx, *gy;
Double_t e;
TVirtualFitter *grFitter = TVirtualFitter::GetFitter();
TMultiGraph *mg = (TMultiGraph*)grFitter->GetObjectFit();
TF1 *f1 = (TF1*)grFitter->GetUserFunc();
Foption_t fitOption = grFitter->GetFitOption();
Int_t fitResult=0;
SetDim(1);
if (fitOption.Robust){
fRobust=kTRUE;
StoreData(kTRUE);
}
SetFormula(f1);
TGraph *gr;
TIter next(mg->GetListOfGraphs());
while ((gr = (TGraph*) next())) {
n = gr->GetN();
gx = gr->GetX();
gy = gr->GetY();
for (i=0; i<n; i++){
if (!f1->IsInside(&gx[i])) continue;
e=gr->GetErrorY(i);
if (e<0 || fitOption.W1)
e=1;
AddPoint(&gx[i], gy[i], e);
}
}
if (fitOption.Robust){
return EvalRobust(h);
}
fitResult = Eval();
if (!fitOption.Nochisq){
Double_t temp, temp2, sumtotal=0;
next.Reset();
while((gr = (TGraph*)next())) {
n = gr->GetN();
gx = gr->GetX();
gy = gr->GetY();
for (i=0; i<n; i++){
if (!f1->IsInside(&gx[i])) continue;
temp=f1->Eval(gx[i]);
temp2=(gy[i]-temp)*(gy[i]-temp);
e=gr->GetErrorY(i);
if (e<0 || fitOption.W1)
e=1;
temp2/=(e*e);
sumtotal+=temp2;
}
}
fChisquare=sumtotal;
f1->SetChisquare(fChisquare);
}
return fitResult;
}
Int_t TLinearFitter::HistLinearFitter()
{
StoreData(kFALSE);
Double_t cu,eu;
Double_t x[3];
Int_t bin,binx,biny,binz;
TH1 *hfit = (TH1*)GetObjectFit();
TF1 *f1 = (TF1*)GetUserFunc();
Int_t fitResult;
Foption_t fitOption = GetFitOption();
SetDim(hfit->GetDimension());
SetFormula(f1);
Int_t hxfirst = GetXfirst();
Int_t hxlast = GetXlast();
Int_t hyfirst = GetYfirst();
Int_t hylast = GetYlast();
Int_t hzfirst = GetZfirst();
Int_t hzlast = GetZlast();
TAxis *xaxis = hfit->GetXaxis();
TAxis *yaxis = hfit->GetYaxis();
TAxis *zaxis = hfit->GetZaxis();
for (binz=hzfirst;binz<=hzlast;binz++) {
x[2] = zaxis->GetBinCenter(binz);
for (biny=hyfirst;biny<=hylast;biny++) {
x[1] = yaxis->GetBinCenter(biny);
for (binx=hxfirst;binx<=hxlast;binx++) {
x[0] = xaxis->GetBinCenter(binx);
if (!f1->IsInside(x)) continue;
bin = hfit->GetBin(binx,biny,binz);
cu = hfit->GetBinContent(bin);
if (fitOption.W1) {
if (fitOption.W1==1 && cu == 0) continue;
eu = 1;
} else {
eu = hfit->GetBinError(bin);
if (eu <= 0) continue;
}
AddPoint(x, cu, eu);
}
}
}
fitResult = Eval();
if (!fitOption.Nochisq){
Double_t temp, temp2, sumtotal=0;
for (binz=hzfirst;binz<=hzlast;binz++) {
x[2] = zaxis->GetBinCenter(binz);
for (biny=hyfirst;biny<=hylast;biny++) {
x[1] = yaxis->GetBinCenter(biny);
for (binx=hxfirst;binx<=hxlast;binx++) {
x[0] = xaxis->GetBinCenter(binx);
if (!f1->IsInside(x)) continue;
bin = hfit->GetBin(binx,biny,binz);
cu = hfit->GetBinContent(bin);
if (fitOption.W1) {
if (fitOption.W1==1 && cu == 0) continue;
eu = 1;
} else {
eu = hfit->GetBinError(bin);
if (eu <= 0) continue;
}
temp=f1->EvalPar(x);
temp2=(cu-temp)*(cu-temp);
temp2/=(eu*eu);
sumtotal+=temp2;
}
}
}
fChisquare=sumtotal;
f1->SetChisquare(fChisquare);
}
return fitResult;
}
Int_t TLinearFitter::EvalRobust(Double_t h)
{
fRobust = kTRUE;
Double_t kEps = 1e-13;
Int_t nmini = 300;
Int_t i, j, maxind=0, k, k1 = 500;
Int_t nbest = 10;
Double_t chi2;
Double_t *bestchi2 = new Double_t[nbest];
for (i=0; i<nbest; i++)
bestchi2[i]=1e30;
Int_t hdef=Int_t((fNpoints+fNfunctions+1)/2);
if (h<0.000001) fH = hdef;
else if (h>0 && h<1 && fNpoints*h > hdef)
fH = Int_t(fNpoints*h);
else {
Warning("Fitting:", "illegal value of H, default is taken");
fH=hdef;
}
fDesign.ResizeTo(fNfunctions, fNfunctions);
fAtb.ResizeTo(fNfunctions);
fParams.ResizeTo(fNfunctions);
Int_t *index = new Int_t[fNpoints];
Double_t *residuals = new Double_t[fNpoints];
if (fNpoints < 2*nmini) {
TMatrixD cstock(fNfunctions, nbest);
for (k = 0; k < k1; k++) {
CreateSubset(fNpoints, fH, index);
chi2 = CStep(1, fH, residuals,index, index, -1, -1);
chi2 = CStep(2, fH, residuals,index, index, -1, -1);
maxind = TMath::LocMax(nbest, bestchi2);
if (chi2 < bestchi2[maxind]) {
bestchi2[maxind] = chi2;
for (i=0; i<fNfunctions; i++)
cstock(i, maxind) = fParams(i);
}
}
Int_t *bestindex = new Int_t[fH];
Double_t currentbest;
for (i=0; i<nbest; i++) {
for (j=0; j<fNfunctions; j++)
fParams(j) = cstock(j, i);
chi2 = 1;
while (chi2 > kEps) {
chi2 = CStep(2, fH, residuals,index, index, -1, -1);
if (TMath::Abs(chi2 - bestchi2[i]) < kEps)
break;
else
bestchi2[i] = chi2;
}
currentbest = TMath::MinElement(nbest, bestchi2);
if (chi2 <= currentbest + kEps) {
for (j=0; j<fH; j++){
bestindex[j]=index[j];
}
maxind = i;
}
for (j=0; j<fNfunctions; j++)
cstock(j, i) = fParams(j);
}
for (j=0; j<fNfunctions; j++)
fParams(j) = cstock(j, maxind);
fFitsample.SetBitNumber(fNpoints, kFALSE);
for (j=0; j<fH; j++){
fFitsample.SetBitNumber(bestindex[j]);
}
if (fInputFunction){
((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
((TF1*)fInputFunction)->SetNumberFitPoints(fH);
((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
}
delete [] index;
delete [] bestindex;
delete [] residuals;
return 0;
}
Int_t indsubdat[5];
for (i=0; i<5; i++)
indsubdat[i] = 0;
Int_t nsub = Partition(nmini, indsubdat);
Int_t hsub;
Int_t sum = TMath::Min(nmini*5, fNpoints);
Int_t *subdat = new Int_t[sum];
RDraw(subdat, indsubdat);
TMatrixD cstockbig(fNfunctions, nbest*5);
Int_t *beststock = new Int_t[nbest];
Int_t i_start = 0;
Int_t i_end = indsubdat[0];
Int_t k2 = Int_t(k1/nsub);
for (Int_t kgroup = 0; kgroup < nsub; kgroup++) {
hsub = Int_t(fH * indsubdat[kgroup]/fNpoints);
for (i=0; i<nbest; i++)
bestchi2[i] = 1e16;
for (k=0; k<k2; k++) {
CreateSubset(indsubdat[kgroup], hsub, index);
chi2 = CStep(1, hsub, residuals, index, subdat, i_start, i_end);
chi2 = CStep(2, hsub, residuals, index, subdat, i_start, i_end);
maxind = TMath::LocMax(nbest, bestchi2);
if (chi2 < bestchi2[maxind]){
for (i=0; i<fNfunctions; i++)
cstockbig(i, nbest*kgroup + maxind) = fParams(i);
bestchi2[maxind] = chi2;
}
}
if (kgroup != nsub - 1){
i_start += indsubdat[kgroup];
i_end += indsubdat[kgroup+1];
}
}
for (i=0; i<nbest; i++)
bestchi2[i] = 1e30;
Int_t hsub2 = Int_t(fH*sum/fNpoints);
for (k=0; k<nbest*5; k++) {
for (i=0; i<fNfunctions; i++)
fParams(i)=cstockbig(i, k);
chi2 = CStep(1, hsub2, residuals, index, subdat, 0, sum);
chi2 = CStep(2, hsub2, residuals, index, subdat, 0, sum);
maxind = TMath::LocMax(nbest, bestchi2);
if (chi2 < bestchi2[maxind]){
beststock[maxind] = k;
bestchi2[maxind] = chi2;
}
}
for (k=0; k<nbest; k++) {
for (i=0; i<fNfunctions; i++)
fParams(i) = cstockbig(i, beststock[k]);
chi2 = CStep(1, fH, residuals, index, index, -1, -1);
chi2 = CStep(2, fH, residuals, index, index, -1, -1);
bestchi2[k] = chi2;
}
maxind = TMath::LocMin(nbest, bestchi2);
for (i=0; i<fNfunctions; i++)
fParams(i)=cstockbig(i, beststock[maxind]);
chi2 = 1;
while (chi2 > kEps) {
chi2 = CStep(2, fH, residuals, index, index, -1, -1);
if (TMath::Abs(chi2 - bestchi2[maxind]) < kEps)
break;
else
bestchi2[maxind] = chi2;
}
fFitsample.SetBitNumber(fNpoints, kFALSE);
for (j=0; j<fH; j++)
fFitsample.SetBitNumber(index[j]);
if (fInputFunction){
((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
((TF1*)fInputFunction)->SetNumberFitPoints(fH);
((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
}
delete [] subdat;
delete [] beststock;
delete [] bestchi2;
delete [] residuals;
delete [] index;
return 0;
}
void TLinearFitter::CreateSubset(Int_t ntotal, Int_t h, Int_t *index)
{
Int_t i, j;
Bool_t repeat=kFALSE;
Int_t nindex=0;
Int_t num;
for(i=0; i<ntotal; i++)
index[i] = ntotal+1;
TRandom r;
for (i=0; i<fNfunctions; i++) {
num=Int_t(r.Uniform(0, 1)*(ntotal-1));
if (i>0){
for(j=0; j<=i-1; j++) {
if(index[j]==num)
repeat = kTRUE;
}
}
if(repeat==kTRUE) {
i--;
repeat = kFALSE;
} else {
index[i] = num;
nindex++;
}
}
fDesign.Zero();
fAtb.Zero();
for (i=0; i<fNfunctions; i++){
AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
}
Bool_t ok;
ok = Linf();
while (!ok && (nindex < h)) {
repeat=kFALSE;
do{
num=Int_t(r.Uniform(0,1)*(ntotal-1));
repeat=kFALSE;
for(i=0; i<nindex; i++) {
if(index[i]==num) {
repeat=kTRUE;
break;
}
}
} while(repeat==kTRUE);
index[nindex] = num;
nindex++;
AddToDesign(TMatrixDRow(fX, index[nindex-1]).GetPtr(), fY(index[nindex-1]), fE(index[nindex-1]));
ok = Linf();
}
}
Double_t TLinearFitter::CStep(Int_t step, Int_t h, Double_t *residuals, Int_t *index, Int_t *subdat, Int_t start, Int_t end)
{
Int_t i, j, itemp, n;
Double_t func;
Double_t val[100];
Int_t npar;
if (start > -1) {
n = end - start;
for (i=0; i<n; i++) {
func = 0;
itemp = subdat[start+i];
if (fInputFunction){
fInputFunction->SetParameters(fParams.GetMatrixArray());
func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
} else {
func=0;
if ((fSpecial>100)&&(fSpecial<200)){
npar = fSpecial-100;
val[0] = 1;
for (j=1; j<npar; j++)
val[j] = val[j-1]*fX(itemp, 0);
for (j=0; j<npar; j++)
func += fParams(j)*val[j];
} else {
if (fSpecial>200) {
npar = fSpecial-201;
func+=fParams(0);
for (j=0; j<npar; j++)
func += fParams(j+1)*fX(itemp, j);
} else {
for (j=0; j<fNfunctions; j++) {
TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
func += fParams(j)*val[j];
}
}
}
}
residuals[i] = (fY(itemp) - func)*(fY(itemp) - func)/(fE(i)*fE(i));
}
} else {
n=fNpoints;
for (i=0; i<fNpoints; i++) {
func = 0;
if (fInputFunction){
fInputFunction->SetParameters(fParams.GetMatrixArray());
func=fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
} else {
func=0;
if ((fSpecial>100)&&(fSpecial<200)){
Int_t npar = fSpecial-100;
val[0] = 1;
for (j=1; j<npar; j++)
val[j] = val[j-1]*fX(i, 0);
for (j=0; j<npar; j++)
func += fParams(j)*val[j];
} else {
if (fSpecial>200) {
Int_t npar = fSpecial-201;
func+=fParams(0);
for (j=0; j<npar; j++)
func += fParams(j+1)*fX(i, j);
} else {
for (j=0; j<fNfunctions; j++) {
TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
val[j] = f1->EvalPar(TMatrixDRow(fX, i).GetPtr());
func += fParams(j)*val[j];
}
}
}
}
residuals[i] = (fY(i) - func)*(fY(i) - func)/(fE(i)*fE(i));
}
}
TMath::KOrdStat(n, residuals, h-1, index);
fDesign.Zero();
fAtb.Zero();
for (i=0; i<h; i++)
AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
Linf();
if (step==1) return 0;
Double_t sum=0;
if (start > -1) {
for (i=0; i<h; i++) {
itemp = subdat[start+index[i]];
if (fInputFunction){
fInputFunction->SetParameters(fParams.GetMatrixArray());
func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
} else {
func=0;
if ((fSpecial>100)&&(fSpecial<200)){
npar = fSpecial-100;
val[0] = 1;
for (j=1; j<npar; j++)
val[j] = val[j-1]*fX(itemp, 0);
for (j=0; j<npar; j++)
func += fParams(j)*val[j];
} else {
if (fSpecial>200) {
npar = fSpecial-201;
func+=fParams(0);
for (j=0; j<npar; j++)
func += fParams(j+1)*fX(itemp, j);
} else {
for (j=0; j<fNfunctions; j++) {
TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
func += fParams(j)*val[j];
}
}
}
}
sum+=(fY(itemp)-func)*(fY(itemp)-func)/(fE(itemp)*fE(itemp));
}
} else {
for (i=0; i<h; i++) {
if (fInputFunction){
fInputFunction->SetParameters(fParams.GetMatrixArray());
func=fInputFunction->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
} else {
func=0;
if ((fSpecial>100)&&(fSpecial<200)){
Int_t npar = fSpecial-100;
val[0] = 1;
for (j=1; j<npar; j++)
val[j] = val[j-1]*fX(index[i], 0);
for (j=0; j<npar; j++)
func += fParams(j)*val[j];
} else {
if (fSpecial>200) {
Int_t npar = fSpecial-201;
func+=fParams(0);
for (j=0; j<npar; j++)
func += fParams(j+1)*fX(index[i], j);
} else {
for (j=0; j<fNfunctions; j++) {
TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
val[j] = f1->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
func += fParams(j)*val[j];
}
}
}
}
sum+=(fY(index[i])-func)*(fY(index[i])-func)/(fE(index[i])*fE(index[i]));
}
}
return sum;
}
Bool_t TLinearFitter::Linf()
{
fDesignTemp2+=fDesignTemp3;
fDesignTemp+=fDesignTemp2;
fDesign+=fDesignTemp;
fDesignTemp3.Zero();
fDesignTemp2.Zero();
fDesignTemp.Zero();
fAtbTemp2+=fAtbTemp3;
fAtbTemp+=fAtbTemp2;
fAtb+=fAtbTemp;
fAtbTemp3.Zero();
fAtbTemp2.Zero();
fAtbTemp.Zero();
fY2+=fY2Temp;
fY2Temp=0;
TDecompChol chol(fDesign);
TVectorD temp(fNfunctions);
Bool_t ok;
temp = chol.Solve(fAtb, ok);
if (!ok){
Error("Linf","Matrix inversion failed");
fParams.Zero();
return kFALSE;
}
fParams = temp;
return ok;
}
Int_t TLinearFitter::Partition(Int_t nmini, Int_t *indsubdat)
{
Int_t nsub;
if ((fNpoints>=2*nmini) && (fNpoints<=(3*nmini-1))) {
if (fNpoints%2==1){
indsubdat[0]=Int_t(fNpoints*0.5);
indsubdat[1]=Int_t(fNpoints*0.5)+1;
} else
indsubdat[0]=indsubdat[1]=Int_t(fNpoints/2);
nsub=2;
}
else{
if((fNpoints>=3*nmini) && (fNpoints<(4*nmini -1))) {
if(fNpoints%3==0){
indsubdat[0]=indsubdat[1]=indsubdat[2]=Int_t(fNpoints/3);
} else {
indsubdat[0]=Int_t(fNpoints/3);
indsubdat[1]=Int_t(fNpoints/3)+1;
if (fNpoints%3==1) indsubdat[2]=Int_t(fNpoints/3);
else indsubdat[2]=Int_t(fNpoints/3)+1;
}
nsub=3;
}
else{
if((fNpoints>=4*nmini)&&(fNpoints<=(5*nmini-1))){
if (fNpoints%4==0) indsubdat[0]=indsubdat[1]=indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
else {
indsubdat[0]=Int_t(fNpoints/4);
indsubdat[1]=Int_t(fNpoints/4)+1;
if(fNpoints%4==1) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
if(fNpoints%4==2) {
indsubdat[2]=Int_t(fNpoints/4)+1;
indsubdat[3]=Int_t(fNpoints/4);
}
if(fNpoints%4==3) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4)+1;
}
nsub=4;
} else {
for(Int_t i=0; i<5; i++)
indsubdat[i]=nmini;
nsub=5;
}
}
}
return nsub;
}
void TLinearFitter::RDraw(Int_t *subdat, Int_t *indsubdat)
{
Int_t jndex = 0;
Int_t nrand;
Int_t i, k, m, j;
Int_t ngroup=0;
for (i=0; i<5; i++) {
if (indsubdat[i]!=0)
ngroup++;
}
TRandom r;
for (k=1; k<=ngroup; k++) {
for (m=1; m<=indsubdat[k-1]; m++) {
nrand = Int_t(r.Uniform(0, 1) * (fNpoints-jndex)) + 1;
jndex++;
if (jndex==1) {
subdat[0] = nrand;
} else {
subdat[jndex-1] = nrand + jndex - 2;
for (i=1; i<=jndex-1; i++) {
if(subdat[i-1] > nrand+i-2) {
for(j=jndex; j>=i+1; j--) {
subdat[j-1] = subdat[j-2];
}
subdat[i-1] = nrand+i-2;
break;
}
}
}
}
}
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.