[u/mrichter/AliRoot.git] / STEER / AliMathBase.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////
// Class AliMathBase
// 
// Subset of  matheamtical functions  not included in the TMath
//

///////////////////////////////////////////////////////////////////////////
#include "TMath.h"
#include "AliMathBase.h"
#include "Riostream.h"
#include "TH1F.h"
#include "TF1.h"
#include "TLinearFitter.h"
 
ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
 
AliMathBase::AliMathBase() : TObject()
{
// Default constructor
}
///////////////////////////////////////////////////////////////////////////
AliMathBase::~AliMathBase()
{
// Destructor
}


//_____________________________________________________________________________
void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
                           , Double_t &sigma, Int_t hh)
{
  //
  // Robust estimator in 1D case MI version - (faster than ROOT version)
  //
  // For the univariate case
  // estimates of location and scatter are returned in mean and sigma parameters
  // the algorithm works on the same principle as in multivariate case -
  // it finds a subset of size hh with smallest sigma, and then returns mean and
  // sigma of this subset
  //

  if (hh==0)
    hh=(nvectors+2)/2;
  Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
  Int_t *index=new Int_t[nvectors];
  TMath::Sort(nvectors, data, index, kFALSE);
  
  Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
  Double_t factor = faclts[TMath::Max(0,nquant-1)];
  
  Double_t sumx  =0;
  Double_t sumx2 =0;
  Int_t    bestindex = -1;
  Double_t bestmean  = 0; 
  Double_t bestsigma = data[index[nvectors-1]]-data[index[0]];   // maximal possible sigma
  for (Int_t i=0; i<hh; i++){
    sumx  += data[index[i]];
    sumx2 += data[index[i]]*data[index[i]];
  }
  
  Double_t norm = 1./Double_t(hh);
  Double_t norm2 = 1./Double_t(hh-1);
  for (Int_t i=hh; i<nvectors; i++){
    Double_t cmean  = sumx*norm;
    Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
    if (csigma<bestsigma){
      bestmean  = cmean;
      bestsigma = csigma;
      bestindex = i-hh;
    }
    
    sumx  += data[index[i]]-data[index[i-hh]];
    sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
  }
  
  Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
  mean  = bestmean;
  sigma = bstd;
  delete [] index;

}


void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh,  Float_t externalfactor)
{
  // Modified version of ROOT robust EvaluateUni
  // robust estimator in 1D case MI version
  // added external factor to include precision of external measurement
  // 

  if (hh==0)
    hh=(nvectors+2)/2;
  Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
  Int_t *index=new Int_t[nvectors];
  TMath::Sort(nvectors, data, index, kFALSE);
  //
  Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
  Double_t factor = faclts[0];
  if (nquant>0){
    // fix proper normalization - Anja
    factor = faclts[nquant-1];
  }

  //
  //
  Double_t sumx  =0;
  Double_t sumx2 =0;
  Int_t    bestindex = -1;
  Double_t bestmean  = 0; 
  Double_t bestsigma = -1;
  for (Int_t i=0; i<hh; i++){
    sumx  += data[index[i]];
    sumx2 += data[index[i]]*data[index[i]];
  }
  //   
  Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
  Double_t norm = 1./Double_t(hh);
  for (Int_t i=hh; i<nvectors; i++){
    Double_t cmean  = sumx*norm;
    Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
    if (csigma<bestsigma ||  bestsigma<0){
      bestmean  = cmean;
      bestsigma = csigma;
      bestindex = i-hh;
    }
    //
    //
    sumx  += data[index[i]]-data[index[i-hh]];
    sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
  }
  
  Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
  mean  = bestmean;
  sigma = bstd;
  delete [] index;
}


//_____________________________________________________________________________
Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist
                        , Int_t *outlist, Bool_t down)
{    
  //
  //  Sort eleements according occurancy 
  //  The size of output array has is 2*n 
  //

  Int_t * sindexS = new Int_t[n];     // temp array for sorting
  Int_t * sindexF = new Int_t[2*n];   
  for (Int_t i=0;i<n;i++) sindexF[i]=0;
  //
  TMath::Sort(n,inlist, sindexS, down);  
  Int_t last      = inlist[sindexS[0]];
  Int_t val       = last;
  sindexF[0]      = 1;
  sindexF[0+n]    = last;
  Int_t countPos  = 0;
  //
  //  find frequency
  for(Int_t i=1;i<n; i++){
    val = inlist[sindexS[i]];
    if (last == val)   sindexF[countPos]++;
    else{      
      countPos++;
      sindexF[countPos+n] = val;
      sindexF[countPos]++;
      last =val;
    }
  }
  if (last==val) countPos++;
  // sort according frequency
  TMath::Sort(countPos, sindexF, sindexS, kTRUE);
  for (Int_t i=0;i<countPos;i++){
    outlist[2*i  ] = sindexF[sindexS[i]+n];
    outlist[2*i+1] = sindexF[sindexS[i]];
  }
  delete [] sindexS;
  delete [] sindexF;
  
  return countPos;

}

//___AliMathBase__________________________________________________________________________
void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
  //
  //
  //
  Int_t nbins    = his->GetNbinsX();
  Float_t nentries = his->GetEntries();
  Float_t sum      =0;
  Float_t mean   = 0;
  Float_t sigma2 = 0;
  Float_t ncumul=0;  
  for (Int_t ibin=1;ibin<nbins; ibin++){
    ncumul+= his->GetBinContent(ibin);
    Float_t fraction = Float_t(ncumul)/Float_t(nentries);
    if (fraction>down && fraction<up){
      sum+=his->GetBinContent(ibin);
      mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
      sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);      
    }
  }
  mean/=sum;
  sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
  if (param){
    (*param)[0] = his->GetMaximum();
    (*param)[1] = mean;
    (*param)[2] = sigma2;
    
  }
  if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
}

void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction,  Bool_t verbose){
  //
  // LTM
  //
  Int_t nbins    = his->GetNbinsX();
  Int_t nentries = (Int_t)his->GetEntries();
  Double_t *data  = new Double_t[nentries];
  Int_t npoints=0;
  for (Int_t ibin=1;ibin<nbins; ibin++){
    Float_t entriesI = his->GetBinContent(ibin);
    Float_t xcenter= his->GetBinCenter(ibin);
    for (Int_t ic=0; ic<entriesI; ic++){
      if (npoints<nentries){
	data[npoints]= xcenter;
	npoints++;
      }
    }
  }
  Double_t mean, sigma;
  Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
  npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
  AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2);
  delete [] data;
  if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
    (*param)[0] = his->GetMaximum();
    (*param)[1] = mean;
    (*param)[2] = sigma;    
  }
}

Double_t  AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
  //
  //  Fit histogram with gaussian function
  //  
  //  Prameters:
  //       return value- chi2 - if negative ( not enough points)
  //       his        -  input histogram
  //       param      -  vector with parameters 
  //       xmin, xmax -  range to fit - if xmin=xmax=0 - the full histogram range used
  //  Fitting:
  //  1. Step - make logarithm
  //  2. Linear  fit (parabola) - more robust - always converge
  //  3. In case of small statistic bins are averaged
  //  
  static TLinearFitter fitter(3,"pol2");
  TVectorD  par(3);
  TVectorD  sigma(3);
  TMatrixD mat(3,3);
  if (his->GetMaximum()<4) return -1;  
  if (his->GetEntries()<12) return -1;  
  if (his->GetRMS()<mat.GetTol()) return -1;
  Float_t maxEstimate   = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((2.*TMath::Pi()*his->GetRMS()));  
  Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));

  if (maxEstimate<1) return -1;
  Int_t nbins    = his->GetNbinsX();
  Int_t npoints=0;
  //


  if (xmin>=xmax){
    xmin = his->GetXaxis()->GetXmin();
    xmax = his->GetXaxis()->GetXmax();
  }
  for (Int_t iter=0; iter<2; iter++){
    fitter.ClearPoints();
    npoints=0;
    for (Int_t ibin=1;ibin<nbins-1; ibin++){      
      Int_t countB=1;
      Float_t entriesI =  his->GetBinContent(ibin);
      for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
	if (ibin+delta>1 &&ibin+delta<nbins-1){
	  entriesI +=  his->GetBinContent(ibin+delta);
	  countB++;
	}
      }
      entriesI/=countB;
      Double_t xcenter= his->GetBinCenter(ibin);
      if (xcenter<xmin || xcenter>xmax) continue;
      Double_t error=1./TMath::Sqrt(countB);
      Float_t   cont=2;
      if (iter>0){
	if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
	cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
	if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
      }
      if (entriesI>1&&cont>1){
	fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
	npoints++;
      }
    }  
    if (npoints>3){
      fitter.Eval();
      fitter.GetParameters(par);
    }else{
      break;
    }
  }
  if (npoints<=3){
    return -1;
  }
  fitter.GetParameters(par);
  fitter.GetCovarianceMatrix(mat);
  if (TMath::Abs(par[1])<mat.GetTol()) return -1;
  if (TMath::Abs(par[2])<mat.GetTol()) return -1;
  Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
  //fitter.GetParameters();
  if (!param)  param  = new TVectorD(3);
  if (!matrix) matrix = new TMatrixD(3,3);
  (*param)[1] = par[1]/(-2.*par[2]);
  (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
  (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1]);
  if (verbose){
    par.Print();
    mat.Print();
    param->Print();
    printf("Chi2=%f\n",chi2);
    TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
    f1->SetParameter(0, (*param)[0]);
    f1->SetParameter(1, (*param)[1]);
    f1->SetParameter(2, (*param)[2]);    
    f1->Draw("same");
  }
  return chi2;
}
Commit	Line	Data
284050f7	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16
	17	///////////////////////////////////////////////////////////////////////////
	18	// Class AliMathBase
	19	//
	20	// Subset of matheamtical functions not included in the TMath
	21	//
	22
	23	///////////////////////////////////////////////////////////////////////////
	24	#include "TMath.h"
	25	#include "AliMathBase.h"
	26	#include "Riostream.h"
f6659a9d	27	#include "TH1F.h"
	28	#include "TF1.h"
	29	#include "TLinearFitter.h"
284050f7	30
	31	ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
	32
	33	AliMathBase::AliMathBase() : TObject()
	34	{
	35	// Default constructor
	36	}
	37	///////////////////////////////////////////////////////////////////////////
	38	AliMathBase::~AliMathBase()
	39	{
	40	// Destructor
	41	}
	42
	43
	44	//_____________________________________________________________________________
	45	void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
	46	, Double_t &sigma, Int_t hh)
	47	{
	48	//
	49	// Robust estimator in 1D case MI version - (faster than ROOT version)
	50	//
	51	// For the univariate case
	52	// estimates of location and scatter are returned in mean and sigma parameters
	53	// the algorithm works on the same principle as in multivariate case -
	54	// it finds a subset of size hh with smallest sigma, and then returns mean and
	55	// sigma of this subset
	56	//
	57
	58	if (hh==0)
	59	hh=(nvectors+2)/2;
	60	Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
	61	Int_t *index=new Int_t[nvectors];
	62	TMath::Sort(nvectors, data, index, kFALSE);
	63
	64	Int_t nquant = TMath::Min(Int_t(Double_t(((hh1./nvectors)-0.5)40))+1, 11);
d9e9045c	65	Double_t factor = faclts[TMath::Max(0,nquant-1)];
284050f7	66
	67	Double_t sumx =0;
	68	Double_t sumx2 =0;
	69	Int_t bestindex = -1;
	70	Double_t bestmean = 0;
	71	Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
	72	for (Int_t i=0; i<hh; i++){
	73	sumx += data[index[i]];
	74	sumx2 += data[index[i]]*data[index[i]];
	75	}
	76
	77	Double_t norm = 1./Double_t(hh);
	78	Double_t norm2 = 1./Double_t(hh-1);
	79	for (Int_t i=hh; i<nvectors; i++){
	80	Double_t cmean = sumx*norm;
	81	Double_t csigma = (sumx2 - hhcmeancmean)*norm2;
	82	if (csigma<bestsigma){
	83	bestmean = cmean;
	84	bestsigma = csigma;
	85	bestindex = i-hh;
	86	}
	87
	88	sumx += data[index[i]]-data[index[i-hh]];
	89	sumx2 += data[index[i]]data[index[i]]-data[index[i-hh]]data[index[i-hh]];
	90	}
	91
	92	Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
	93	mean = bestmean;
	94	sigma = bstd;
	95	delete [] index;
	96
	97	}
	98
	99
	100
	101	void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor)
	102	{
	103	// Modified version of ROOT robust EvaluateUni
	104	// robust estimator in 1D case MI version
	105	// added external factor to include precision of external measurement
	106	//
	107
	108	if (hh==0)
	109	hh=(nvectors+2)/2;
	110	Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
	111	Int_t *index=new Int_t[nvectors];
	112	TMath::Sort(nvectors, data, index, kFALSE);
	113	//
	114	Int_t nquant = TMath::Min(Int_t(Double_t(((hh1./nvectors)-0.5)40))+1, 11);
	115	Double_t factor = faclts[0];
	116	if (nquant>0){
	117	// fix proper normalization - Anja
	118	factor = faclts[nquant-1];
	119	}
	120
	121	//
	122	//
	123	Double_t sumx =0;
	124	Double_t sumx2 =0;
	125	Int_t bestindex = -1;
	126	Double_t bestmean = 0;
	127	Double_t bestsigma = -1;
	128	for (Int_t i=0; i<hh; i++){
	129	sumx += data[index[i]];
130	sumx2 += data[index[i]]*data[index[i]];
131	}
132	//
133	Double_t kfactor = 2.externalfactor - externalfactorexternalfactor;
134	Double_t norm = 1./Double_t(hh);
135	for (Int_t i=hh; i<nvectors; i++){
136	Double_t cmean = sumx*norm;
137	Double_t csigma = (sumx2norm - cmeancmean*kfactor);
138	if (csigma<bestsigma \|\| bestsigma<0){
139	bestmean = cmean;
140	bestsigma = csigma;
141	bestindex = i-hh;
142	}
143	//
144	//
145	sumx += data[index[i]]-data[index[i-hh]];
146	sumx2 += data[index[i]]data[index[i]]-data[index[i-hh]]data[index[i-hh]];
147	}
148
149	Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
150	mean = bestmean;
151	sigma = bstd;
152	delete [] index;
153	}
154
155
156	//_____________________________________________________________________________
157	Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist
158	, Int_t *outlist, Bool_t down)
159	{
160	//
161	// Sort eleements according occurancy
162	// The size of output array has is 2*n
163	//
164
165	Int_t * sindexS = new Int_t[n]; // temp array for sorting
166	Int_t * sindexF = new Int_t[2*n];
167	for (Int_t i=0;i<n;i++) sindexF[i]=0;
168	//
169	TMath::Sort(n,inlist, sindexS, down);
170	Int_t last = inlist[sindexS[0]];
171	Int_t val = last;
172	sindexF[0] = 1;
173	sindexF[0+n] = last;
174	Int_t countPos = 0;
175	//
176	// find frequency
177	for(Int_t i=1;i<n; i++){
178	val = inlist[sindexS[i]];
179	if (last == val) sindexF[countPos]++;
180	else{
181	countPos++;
182	sindexF[countPos+n] = val;
183	sindexF[countPos]++;
184	last =val;
185	}
186	}
187	if (last==val) countPos++;
188	// sort according frequency
189	TMath::Sort(countPos, sindexF, sindexS, kTRUE);
190	for (Int_t i=0;i<countPos;i++){
191	outlist[2*i ] = sindexF[sindexS[i]+n];
192	outlist[2*i+1] = sindexF[sindexS[i]];
193	}
194	delete [] sindexS;
195	delete [] sindexF;
196
197	return countPos;
198
199	}
f6659a9d	200
	201	//___AliMathBase__________________________________________________________________________
	202	void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
	203	//
	204	//
	205	//
	206	Int_t nbins = his->GetNbinsX();
	207	Float_t nentries = his->GetEntries();
	208	Float_t sum =0;
	209	Float_t mean = 0;
	210	Float_t sigma2 = 0;
	211	Float_t ncumul=0;
	212	for (Int_t ibin=1;ibin<nbins; ibin++){
	213	ncumul+= his->GetBinContent(ibin);
	214	Float_t fraction = Float_t(ncumul)/Float_t(nentries);
	215	if (fraction>down && fraction<up){
	216	sum+=his->GetBinContent(ibin);
	217	mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
	218	sigma2+=his->GetBinCenter(ibin)his->GetBinCenter(ibin)his->GetBinContent(ibin);
	219	}
	220	}
	221	mean/=sum;
	222	sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
	223	if (param){
	224	(*param)[0] = his->GetMaximum();
	225	(*param)[1] = mean;
	226	(*param)[2] = sigma2;
	227
	228	}
	229	if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
	230	}
	231
	232	void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){
	233	//
	234	// LTM
	235	//
	236	Int_t nbins = his->GetNbinsX();
	237	Int_t nentries = (Int_t)his->GetEntries();
	238	Double_t *data = new Double_t[nentries];
	239	Int_t npoints=0;
	240	for (Int_t ibin=1;ibin<nbins; ibin++){
	241	Float_t entriesI = his->GetBinContent(ibin);
	242	Float_t xcenter= his->GetBinCenter(ibin);
	243	for (Int_t ic=0; ic<entriesI; ic++){
	244	if (npoints<nentries){
	245	data[npoints]= xcenter;
	246	npoints++;
	247	}
	248	}
	249	}
	250	Double_t mean, sigma;
	251	Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
	252	npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
	253	AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2);
	254	delete [] data;
	255	if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
	256	(*param)[0] = his->GetMaximum();
	257	(*param)[1] = mean;
	258	(*param)[2] = sigma;
	259	}
	260	}
	261
	262	Double_t AliMathBase::FitGaus(TH1F* his, TVectorD param, TMatrixD matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
	263	//
264	// Fit histogram with gaussian function
265	//
266	// Prameters:
267	// return value- chi2 - if negative ( not enough points)
268	// his - input histogram
269	// param - vector with parameters
270	// xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used
271	// Fitting:
272	// 1. Step - make logarithm
273	// 2. Linear fit (parabola) - more robust - always converge
274	// 3. In case of small statistic bins are averaged
275	//
276	static TLinearFitter fitter(3,"pol2");
277	TVectorD par(3);
278	TVectorD sigma(3);
279	TMatrixD mat(3,3);
280	if (his->GetMaximum()<4) return -1;
281	if (his->GetEntries()<12) return -1;
282	if (his->GetRMS()<mat.GetTol()) return -1;
283	Float_t maxEstimate = his->GetEntries()his->GetBinWidth(1)/TMath::Sqrt((2.TMath::Pi()*his->GetRMS()));
284	Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
285
286	if (maxEstimate<1) return -1;
287	Int_t nbins = his->GetNbinsX();
288	Int_t npoints=0;
289	//
290
291
292	if (xmin>=xmax){
293	xmin = his->GetXaxis()->GetXmin();
294	xmax = his->GetXaxis()->GetXmax();
295	}
296	for (Int_t iter=0; iter<2; iter++){
297	fitter.ClearPoints();
298	npoints=0;
299	for (Int_t ibin=1;ibin<nbins-1; ibin++){
300	Int_t countB=1;
301	Float_t entriesI = his->GetBinContent(ibin);
302	for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
303	if (ibin+delta>1 &&ibin+delta<nbins-1){
304	entriesI += his->GetBinContent(ibin+delta);
305	countB++;
306	}
307	}
308	entriesI/=countB;
309	Double_t xcenter= his->GetBinCenter(ibin);
310	if (xcenter<xmin \|\| xcenter>xmax) continue;
311	Double_t error=1./TMath::Sqrt(countB);
312	Float_t cont=2;
313	if (iter>0){
314	if (par[0]+par[1]xcenter+par[2]xcenter*xcenter>20) return 0;
315	cont = TMath::Exp(par[0]+par[1]xcenter+par[2]xcenter*xcenter);
316	if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
317	}
318	if (entriesI>1&&cont>1){
319	fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
320	npoints++;
321	}
322	}
323	if (npoints>3){
324	fitter.Eval();
325	fitter.GetParameters(par);
326	}else{
327	break;
328	}
329	}
330	if (npoints<=3){
331	return -1;
332	}
333	fitter.GetParameters(par);
334	fitter.GetCovarianceMatrix(mat);
335	if (TMath::Abs(par[1])<mat.GetTol()) return -1;
336	if (TMath::Abs(par[2])<mat.GetTol()) return -1;
337	Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
338	//fitter.GetParameters();
339	if (!param) param = new TVectorD(3);
340	if (!matrix) matrix = new TMatrixD(3,3);
341	(param)[1] = par[1]/(-2.par[2]);
342	(param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.par[2]));
343	(param)[0] = TMath::Exp(par[0]+ par[1] (param)[1] + par[2](param)[1](*param)[1]);
344	if (verbose){
345	par.Print();
346	mat.Print();
347	param->Print();
348	printf("Chi2=%f\n",chi2);
349	TF1 * f1= new TF1("f1","[0]exp(-(x-[1])^2/(2[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
350	f1->SetParameter(0, (*param)[0]);
351	f1->SetParameter(1, (*param)[1]);
352	f1->SetParameter(2, (*param)[2]);
353	f1->Draw("same");
354	}
355	return chi2;
356	}
357
358