#include "T49DedxFunction.h"
#include "TAxis.h"
#include "TMatrix.h"
#include "TH1.h"
#include "TMinuit.h"
#include "TVector.h"
/*
$Log: T49DedxFunction.C,v $
Revision 1.5  2005/02/08 10:22:20  flierl
TVector.h included

Revision 1.4  2000/07/07 14:12:11  cblume
Made ROOT v2.25 compatibel by Marco

 * Revision 1.3  2000/04/13  15:21:09  cblume
 * Bug fixes from Marco
 *
 * Revision 1.2  2000/01/30  14:44:12  cblume
 * New versions from Marco with updated comments
 *
 * Revision 1.1  2000/01/17  12:38:08  cblume
 * Add T49DedxFunction.C
 *
*/

//____________________________________________________________________
//  dE/dx fit function class
//  
//  This is an interface to a TF1, used to fit a dE/dx spectra.
//  The function consists of four gaussians, with identical relative
//  widths. Different methods are provide to fix or release the
//  positions of the different peaks. Internally, the pion peak position
//  is absolute, whereas the other peak positions are relative
//  to it.

ClassImp(T49DedxFunction)

 T49DedxFunction::T49DedxFunction(Char_t* name): 
     TF1(name,
     	"[6]*exp(-0.5*(x-[2])*(x-[2])/([2]*[2]*[0]))+[5]*exp(-0.5*(x-[1]*[2])*(x-[1]*[2])/([1]*[2]*[1]*[2]*[0]))+[7]*exp(-0.5*(x-[3]*[2])*(x-[3]*[2])/([3]*[2]*[3]*[2]*[0]))+[8]*exp(-0.5*(x-[4]*[2])*(x-[4]*[2])/([4]*[2]*[4]*[2]*[0]))")
{
  SetParNames("SqRelRes","ElRelPos","PiPos","KaRelPos","PRelPos","ElAmp","PiAmp","KaAmp","PAmp");
  SetParLimits(1,1,-1);
  SetParLimits(3,1,-1);
  SetParLimits(4,1,-1);

  fAmpFix[0]=0;
  fAmpFix[1]=0;
  fAmpFix[2]=0;
  fAmpFix[3]=0;

  fPtot = -1;

  // Parameters for Bethe Bloch-parametrization     
  fParaC =  1.597546;// 1.5624;
  fParaD =  9.8;
  fParaE =  2.38;
  fParaF =  0.2;
  
  fStored=0;

}

 void T49DedxFunction::SetFree()
{
  //  Release bounds on relative positions (N.B. the pion position 
  //  is not affected)

  if (!fAmpFix[0] || GetAmplitude(0) > 0) SetParLimits(1,0,0);
  if (!fAmpFix[2] || GetAmplitude(2) > 0) SetParLimits(3,0,0);
  if (!fAmpFix[3] || GetAmplitude(3) > 0) SetParLimits(4,0,0);
}

 void T49DedxFunction::SetFreePosition(Int_t i,Float_t x)
{
  //  Release bound on peak position i, and set position to x.

  SetPosition(i,x);
  SetParLimits(1+i,0,0);
}

 void T49DedxFunction::SetFixedPosition(Int_t i,Float_t x)
{
//  Set fixed value for peak position i

  SetPosition(i,x);
  SetParLimits(1+i,1,-1);
}

 void T49DedxFunction::SetFreePosition(Int_t i)
{
  //  Release bounds on position of peak i

  if (!fAmpFix[i] || GetAmplitude(i)>0)
    SetParLimits(1+i,0,0);
}

 void T49DedxFunction::SetFixedPosition(Int_t i)
{
  //  Fix position of peak i

  SetParLimits(1+i,1,-1);
}

 void T49DedxFunction::ResetFreePosition(Int_t i)
{
  //  Set position to Bethe Bloch value, and 
 
  if (i==1)
    SetParameter(2,GetRelRise(1));
  else
    SetParameter(1+i,GetRelRise(i)/GetRelRise(1));
  if (!fAmpFix[i] || GetAmplitude(i)>0)    
    SetParLimits(1+i,0,0);
}


 void T49DedxFunction::SetPosition(Int_t i, Float_t x, Float_t lowlim, Float_t uplim)
{
  //  Set position and bounds. (N.B.: the bounds of non-pion peaks 
  //  change with the pion peak position, since they are relative positions)

  if (i==1) {
    for (Int_t p=0;p<4;p++) {
      if (p!=1) SetParameter(1+p,GetPosition(p)/x);
    }
    SetParameter(2,x);
    if (!fAmpFix[i] || GetAmplitude(i)>0)
      SetParLimits(2,lowlim,uplim);
  }
  else if (GetParameter(2)!=0) {
    SetParameter(1+i,x/GetParameter(2));
    if (!fAmpFix[i] || GetAmplitude(i)>0)
      SetParLimits(1+i,lowlim/GetParameter(2),uplim/GetParameter(2));
  }
  else printf("First set Pion Position!! n");
}

 void T49DedxFunction::InitPositions(Float_t pipos)
{
  //  Sets pion position to pipos, and relative positions of other peaks 
  //  acoording to Bethe-Bloch parametrization

  SetPosition(1,pipos);
  for (Int_t i=0;i<4;i++) if (i!=1) SetParameter(1+i,GetRelRise(i)/GetRelRise(1));
  SetParLimits(1,1,-1);
  SetParLimits(3,1,-1);
  SetParLimits(4,1,-1);
}

 void T49DedxFunction::SetAmplitude(Int_t i, Float_t x,Float_t lowlim, Float_t uplim)
{
  SetParameter(5+i,x);
  SetParLimits(5+i,lowlim,uplim);
}

 void T49DedxFunction::SetAmplitude(Int_t i, Float_t x)
{
  SetParameter(5+i,x);
}

 void T49DedxFunction::SetFixedAmplitude(Int_t i, Float_t x)
{
  //  Facility to fix the amplitude. Mainly used to fix amplitudes at zero.

  SetParameter(5+i,x);
  SetParLimits(5+i,1,-1);
  fAmpFix[i]=1;
  if (x==0) SetFixedPosition(i);
}

 void T49DedxFunction::SetFixedAmplitude(Int_t i)
{
  SetParLimits(5+i,1,-1);
  fAmpFix[i]=1;
}

 void T49DedxFunction::SetResolution(Float_t x)
{
  SetParameter(0,x*x);
}

 Float_t T49DedxFunction::GetPosition(Int_t i)
{
  if (i==1) 
    return GetParameter(2);
  else
    return GetParameter(1+i)*GetParameter(2);
}

 Float_t T49DedxFunction::GetRelPosition(Int_t i)
{
  if (i!=1) return GetParameter(1+i);
  printf("This is probably a mistake, Relative Pion Position == 1n");
  return 1;
}

 Float_t T49DedxFunction::GetAmplitude(Int_t i)
{
  return GetParameter(5+i);
}

 Float_t T49DedxFunction::GetResolution()
{
  return sqrt(GetParameter(0));
}

 void T49DedxFunction::GetRatio(Int_t part, Float_t &val, Float_t &err)
{
  Float_t sum=0, err_sq=0, p_cnt=0, err_cnt_sq=0;
  Float_t Der[9];
  Der[0]=0;
  for (Int_t p=0; p<4; p++) {
    sum+=GetResolution()*GetPosition(p)*GetAmplitude(p);
    Der[0]+=GetPosition(p)*GetAmplitude(p);
    Der[1+2*p]=GetResolution()*GetPosition(p);
    Der[2+2*p]=GetResolution()*GetAmplitude(p);
  }
  for (Int_t i=0; i<9; i++)
    for (Int_t j=0; j<9; j++)
      err_sq+=Der[i]*fCovMat[i][j]*Der[j];

  p_cnt=GetResolution()*GetPosition(part)*GetAmplitude(part);
  for (Int_t i=0; i<9; i++) 
    Der[i]=0;

  Der[0]=GetPosition(part)*GetAmplitude(part);
  Der[1+2*part]=GetResolution()*GetPosition(part);
  Der[2+2*part]=GetResolution()*GetAmplitude(part);
  for (Int_t i=0; i<9; i++)
    for (Int_t j=0; j<9; j++)
      err_cnt_sq+=Der[i]*fCovMat[i][j]*Der[j];

  val=p_cnt/sum;
  err=sqrt(err_sq/sum/sum+err_cnt_sq/p_cnt/p_cnt);

}

 void T49DedxFunction::GetRatio2(Int_t part, Float_t &val, Float_t &err)
{
  Float_t denom=0;
  Float_t num=0;
  Float_t Der[9];
  num=GetPosition(part)*GetAmplitude(part);  
  for (Int_t p=0; p<4; p++) 
    denom+=GetPosition(p)*GetAmplitude(p);
  val=num/denom;

  Float_t denomsq=denom*denom;
  Der[0]=0;
  for (Int_t p=0; p<4; p++) 
    if (p==part) {
      Der[1+2*p]=GetPosition(p)/denom-num/denomsq*GetPosition(p);
      Der[2+2*p]=GetAmplitude(p)/denom-num/denomsq*GetAmplitude(p);
    }
    else {
      Der[1+2*p]=-num/denomsq*GetPosition(p);
      Der[2+2*p]=-num/denomsq*GetAmplitude(p);
    }

  err=0;
  for (Int_t i=1; i<9; i++)
    for (Int_t j=1; j<9; j++)
      err+=Der[i]*fCovMat[i][j]*Der[j];
  err=sqrt(err);
}

 void T49DedxFunction::StoreCovMat()
{
  Double_t *pCovMat=(Double_t *) fCovMatMinuit;
  gMinuit->mnemat(pCovMat,9);
  if (gMinuit->GetNumFixedPars()!=0) {
    Double_t val,err;
    // Remap matrix, taking care of fixed parameters
    for (Int_t i=8;i>=0;i--) {
      Int_t k=gMinuit->GetParameter(i,val,err);
      for (Int_t j=8;j>=0;j--) {
	Int_t l=gMinuit->GetParameter(j,val,err);
	if (k==0||l==0) 
	  fCovMatMinuit[i][j]=0;
	else
	  fCovMatMinuit[i][j]=fCovMatMinuit[k-1][l-1];
      }
    }
  }
  // Project onto new parameters: res, elamp,elpos,piamp,pipos etc...
  Float_t DerMat[9][9];
  for (Int_t i=0;i<9;i++)
    for (Int_t j=0;j<9;j++)
      DerMat[i][j]=0;
  //res
  DerMat[0][0]=0.5/sqrt(GetParameter(0));
  //elamp and pos
  DerMat[1][5]=1.;
  DerMat[2][1]=GetParameter(2);
  DerMat[2][2]=GetParameter(1);
  //piamp and pos
  DerMat[3][6]=1.;
  DerMat[4][2]=1.;
  //kaamp and pos
  DerMat[5][7]=1.;
  DerMat[6][3]=GetParameter(2);
  DerMat[6][2]=GetParameter(3);
  //pramp and pos
  DerMat[7][8]=1.;
  DerMat[8][4]=GetParameter(2);
  DerMat[8][2]=GetParameter(4);
  for (Int_t i=0;i<9;i++)
    for (Int_t j=0;j<9;j++) {
      fCovMat[i][j]=0;
      for (Int_t k=0;k<9;k++)
	for (Int_t l=0;l<9;l++) 
	  fCovMat[i][j]+=DerMat[i][k]*fCovMatMinuit[k][l]*DerMat[j][l];
    }
}

 void T49DedxFunction::StoreParameters()
{
  //  The class may hold one alternative set of parameters. This method 
  //  stores the parameters internally.

  GetParameters(fStoredPars);
  for (Int_t i=0; i<9; i++) 
    GetParLimits(i,fStoredLowLim[i],fStoredUpLim[i]);
  fStored=1;
}

 void T49DedxFunction::RestoreParameters()
{
  //  Sets the parameters to the values stored using StoreParameters().

  if (fStored) {
    SetParameters(fStoredPars);
    for (Int_t i=0; i<9; i++) { 
      SetParLimits(i,fStoredLowLim[i],fStoredUpLim[i]);
      if (i>5) {
	if (fStoredLowLim[i]>fStoredUpLim[i]) 
	fAmpFix[i-5]=1;
      else
	fAmpFix[i-5]=0;
      }
    }
  }
  else printf("No parameters stored for this function.n");
}

 Float_t T49DedxFunction::GetRelRise(Int_t part)
{
  //  Returns the value of dE/dx from the Bethe-Bloch-parametrization. First 
  //  call SetPtot()

  if (fPtot>0)
    return RelRise(fPtot,part);
  else {
    printf("First call SetPtot()!!! n");
    return 0;
  }
}

 void T49DedxFunction::Init(TH1F* Hist,Char_t* Opt)
{
//  This method uses the histogram to provide an estimate for the peak
//  positions, amplitudes and widths. Relies heavily on the Bethe-Bloch
//  parametrization, so first call SetPtot(). Option "V" triggers
//  verbose behaviour.

  Float_t dedx[4];
  Float_t pos[4];
  Bool_t Verbose=0;
  if (strspn("V",Opt)) Verbose=1;
  
  dedx[0] = GetRelRise(0);
  dedx[1] = GetRelRise(1);
  dedx[2] = GetRelRise(2);
  dedx[3] = GetRelRise(3);

  pos[0] = dedx[0]/dedx[1];
  pos[1] = dedx[1]/dedx[1];  
  pos[2] = dedx[2]/dedx[1];
  pos[3] = dedx[3]/dedx[1];

// Set up peak search window.

  TAxis *Xax=Hist->GetXaxis();
  Int_t Bin = (Xax->FindBin(dedx[1])+Xax->FindBin(dedx[3]))/2;
  const Int_t Window=35;
  Float_t Smooth[2*Window+1];
  Int_t LeftBin = Bin-Window;
  Int_t RightBin = Bin+Window;
  if (LeftBin<1) LeftBin=1;
  if (RightBin>Xax->GetNbins()) RightBin = Xax->GetNbins();

// Smooth data

  if (Verbose) printf("Smoothing histogram, bin %d to %d, X= %f to %f n",
                      Bin-Window,Bin+Window,Hist->GetBinCenter(Bin-Window),Hist->GetBinCenter(Bin+Window));
  for (Int_t i = Bin-Window;i <= Bin+Window;i++)
     if (i<1 || i>Xax->GetNbins()) 
        Smooth[i-Bin+Window]=0;
     else
        Smooth[i-Bin+Window]=(Hist->GetBinContent(i-1)+Hist->GetBinContent(i)+Hist->GetBinContent(i+1))/3;

// Find maxima
  
  // Maximum number of peaks is quite high, since some histos are ragged
  const Int_t MaxPeakCount = 30;
  Int_t PeakIndex[MaxPeakCount];
  Float_t Peaks[MaxPeakCount];
  Int_t PeakQual[MaxPeakCount];
  Int_t PeakCount=0;
  if (Verbose) printf("Searching for maximum, trying: n");
  for (Int_t i=1; i < 2*Window; i++) {
    if ((Smooth[i]>Smooth[i-1] && Smooth[i]>=Smooth[i+1]) ||
	(Smooth[i]>=Smooth[i-1] && Smooth[i]>Smooth[i+1])) {
      PeakCount++;
      if (PeakCount-1<MaxPeakCount) {
	PeakIndex[PeakCount-1] = i;
	Peaks[PeakCount-1] = Smooth[i];
	if ((i>1 && i<2*Window-1) &&
	    (Smooth[i]>=Smooth[i-1] && Smooth[i]>Smooth[i-2]) &&
	    (Smooth[i]>=Smooth[i+1] && Smooth[i]>Smooth[i+2])) {
	  PeakQual[PeakCount-1] = 2; 
	} else PeakQual[PeakCount-1] = 1;

	if (PeakCount > 1 && i-PeakIndex[PeakCount-2] < 5) { 
	  //if peak very close to former peak
	  Int_t l = PeakIndex[PeakCount-2];
	  if ((l>1 && i<2*Window-1) && 
	      (Smooth[l]>Smooth[l-1] && Smooth[l]>Smooth[l-2]) &&
	      (Smooth[i]>Smooth[i+1] && Smooth[i]>Smooth[i+2])) {
	    //and peaks seem to be part of one peak, try to merge peaks
	    PeakQual[PeakCount-2]=2;
	    PeakQual[PeakCount-1]=2;
	  } 
	}

	if (Verbose) printf("Maximum #%d found, at bin %d, smoothed value: %f, Quality: %dn",PeakCount,i+Bin-Window,Smooth[i],PeakQual[PeakCount-1]);
      }
      else printf("WARNING: too many maxima found in T49DedxFunction::Init. n");
    }
  }

  // Initialize xvals, yvals, etc for PolyFit:

  Float_t xvals[2*Window+1];
  Float_t yvals[2*Window+1];
  Float_t errs[2*Window+1];
  Int_t Idx;
  for (Int_t i=LeftBin;i<=RightBin;i++) {
    Idx = i-Bin+Window;
    xvals[Idx] = Xax->GetBinCenter(i);
    yvals[Idx] = Hist->GetBinContent(i);
    errs[Idx] = sqrt(yvals[Idx]); 
  }
  Float_t pars[3];

  Int_t Quality=2;
  const Int_t MaxFits=10;
  Int_t NFits=0;
  Int_t LeftBound[MaxFits], RightBound[MaxFits];
  Float_t Pos[MaxFits], SigSq[MaxFits], Amp[MaxFits];  

  while (Quality>0 && NFits==0) {
    for (Int_t i=0; i<PeakCount; i++) {
      if (PeakQual[i] == Quality && Smooth[PeakIndex[i]]>7) {
	Bool_t InRange=0;
	for (Int_t j=0; j<NFits; j++) 
	  if (PeakIndex[i]>LeftBound[j] && PeakIndex[i]<RightBound[j]) InRange=1;
	if (!InRange && NFits < MaxFits) {
	  LeftBound[NFits]=0;
	  RightBound[NFits]=0;
	  // Find region > 0.7*Max, or stop where the smoothed function starts to
	  // rise again...
	  if (Verbose) printf("Probing left boundary... Trying bin ");
	  for (Int_t k = PeakIndex[i];k>1;k--) {
	    if (Verbose) printf("%d, ",k+Bin-Window);
	    if ((Smooth[k]<0.7*Peaks[i])||
		((Smooth[k-1]>Smooth[k]) && (Smooth[k-2]>Smooth[k-1]))) {
	      LeftBound[NFits]=k;
	      if (Verbose) printf("found!! n");
	      break;
	    }
	  }
	  if (Verbose) printf("Probing right boundary... Trying bin ");
	  for (Int_t k = PeakIndex[i];k<2*Window-1;k++) {
	    if (Verbose) printf("%d, ",k+Bin-Window);
	    if ((Smooth[k]<0.7*Peaks[i])||
		((Smooth[k+1]>Smooth[k]) && (Smooth[k+2]>Smooth[k+1]))) {
	      RightBound[NFits]=k;
	      if (Verbose) printf("found!! n");
	      break;
	    }
	  }
	  if (LeftBound[NFits]==0) {
	    printf("WARNING: Peak search goes out of window, T49DedxFunction::Init(), letfb n");
	  }
	  if (RightBound[NFits]==0) {
	    printf("WARNING: Peak search goes out of window, T49DedxFunction::Init(), rigthb n");
	    RightBound[NFits]=2*Window;
	  }
	  if (Verbose) printf("Boundaries set at: left %d, %f; right %d, %f n",LeftBound[NFits]+Bin-Window,xvals[LeftBound[NFits]],RightBound[NFits]+Bin-Window,xvals[RightBound[NFits]]);
	  if (RightBound[NFits]-LeftBound[NFits]>3) {
	    PolyFit(RightBound[NFits]-LeftBound[NFits]+1,&xvals[LeftBound[NFits]],&yvals[LeftBound[NFits]],&errs[LeftBound[NFits]],2,pars);
	    //  PolyFit(rightb-leftb+1,&xvals[Idx],&yvals[Idx],2,pars);
	    Pos[NFits] = -pars[1]/2/pars[2];
	    SigSq[NFits] = -0.5/pars[2]*(pars[0]-0.25*pars[1]*pars[1]/pars[2]);
	    Amp[NFits] = -2*pars[2]*SigSq[NFits];
	    if (Verbose) printf("Fit # %d, Pos: %f, Width: %f, Amp: %f, ",NFits+1,Pos[NFits],SigSq[NFits],Amp[NFits]);
	    if (SigSq[NFits]>0 &&
		Pos[NFits]>xvals[LeftBound[NFits]] && 
		Pos[NFits]<xvals[RightBound[NFits]]) {
	      NFits++;
	      if (Verbose) printf("kept. n");
	    }
	    else if (Verbose) printf("discarded. n");
	    if (Verbose && NFits==MaxFits) printf("Maximum number of fits (%d) reached, no more peak will be fittedn",MaxFits);
	  }
	}
      }
    }
    Quality--;
    if (NFits==0) printf("WARNING: Poor peak quality in T49DedxFunction::Init. n");
  }
  Float_t PiPos,PiAmp,Res,PrAmp;
  Float_t KaAmp = 0;
  if (NFits==0) {
    printf("WARNING: No peaks found in T49DedxFunction::Init. n");
      PiPos=dedx[1];
      PiAmp=100;
      PrAmp=10;
      Res=0.001;
  }
  else {
    Int_t IdxPi=0;
    Float_t DistPi=fabs(Pos[0]-dedx[1]);
    Int_t IdxPr=-1;
    Float_t DistPr=fabs(Pos[0]-dedx[3])+0.1;
    Int_t IdxKa=-1;
    Float_t DistKa=fabs(Pos[0]-dedx[2])+0.1;
    for (Int_t i=0; i<NFits; i++) {
      if (fabs(Pos[i]-dedx[1])<DistPi) {
	DistPi=fabs(Pos[i]-dedx[1]);
	IdxPi=i;
      } 
      if (fabs(Pos[i]-dedx[3]) < DistPr && 
	  fabs(Pos[i]-dedx[1]) > fabs(Pos[i]-dedx[3])) {
	DistPr=fabs(Pos[i]-dedx[3]);
	IdxPr=i;
      }
      if (fabs(Pos[i]-dedx[2])<DistKa &&
	  fabs(Pos[i]-dedx[1]) > fabs(Pos[i]-dedx[2])) {
	DistKa=fabs(Pos[i]-dedx[2]);
	IdxKa=i;
      }
    }
    if (Verbose) printf("Proton peak: fit# %d, Kaon peak: fit# %d, Pion peak: fit# %d n",IdxPr+1,IdxKa+1,IdxPi+1);
    if (Verbose) printf("Amp[IdxPi]= %f, Amp[IdxPr]=%f n",
			IdxPi>=0?Amp[IdxPi]:0,IdxPr>=0?Amp[IdxPr]:0);
    // Initialise Pion peak
    PiPos=Pos[IdxPi];
    PiAmp=Amp[IdxPi];
    Res=SigSq[IdxPi];

    // Possibly override, due to large proton peak 
    // (also valid if largest peak is proton peak)

    if (IdxPr >= 0 && Amp[IdxPi]<=Amp[IdxPr]) {
      if (DistKa < DistPr) {
	PiPos=Pos[IdxKa]/pos[2];
	KaAmp=Amp[IdxKa];
	Res=SigSq[IdxKa]/pos[2]/pos[2];
	PrAmp=0.05*KaAmp;
      }
      else {
	PiPos=Pos[IdxPr]/pos[3];
	PrAmp=Amp[IdxPr];
	Res=SigSq[IdxPr]/pos[3]/pos[3];
	KaAmp=0.05*PrAmp;
      }
      if (IdxPr==IdxPi) 
	PiAmp=0.05*PrAmp;
      else { 
	PiAmp=Amp[IdxPi]-PrAmp*exp(-0.5*(pos[3]-1)*PiPos*(pos[3]-1)*PiPos/Res)-
	  KaAmp*exp(-0.5*(pos[2]-1)*PiPos*(pos[2]-1)*PiPos/Res);
	if (PiAmp<0)
	  PiAmp=0.05*PrAmp;
      }
    }
    else {
      //if (IdxKa==IdxPi || ( DistKa > fabs(PiPos-dedx[2]) && DistPr > fabs(PiPos-dedx[3]))) {
      //printf("WARNING: Second peak not reliable in T49DedxFunction::Init.n");
      // If no second peak, ready
      if (IdxKa < 0 && IdxPr < 0) {
	PrAmp=0.05*PiAmp;
	KaAmp=0.05*PiAmp;
      }
      else {
	// Else, recalculate proton/kaon peak
	// Start addition 991215
	//if (Verbose) printf("Subtracting pion peak, new points: n");
	for (Int_t i=0; i<2*Window+1; i++) {
	  yvals[i]=yvals[i]-PiAmp*exp(-0.5*(xvals[i]-PiPos)*(xvals[i]-PiPos)/Res);
	  //if (Verbose) printf("%f, %f n",xvals[i],yvals[i]);
	}
	// refit kaon
	if (IdxKa >= 0) {
	  PolyFit(RightBound[IdxKa]-LeftBound[IdxKa]+1,&xvals[LeftBound[IdxKa]],&yvals[LeftBound[IdxKa]],&errs[LeftBound[IdxKa]],2,pars);
	  Pos[IdxKa] = -pars[1]/2/pars[2];
	  SigSq[IdxKa] = -0.5/pars[2]*(pars[0]-0.25*pars[1]*pars[1]/pars[2]);
	  Amp[IdxKa] = -2*pars[2]*SigSq[IdxKa];
	  if (Verbose) printf("Refitted kaon (#%d), Pos: %f, Width: %f, Amp: %fn",IdxKa+1,Pos[IdxKa],SigSq[IdxKa],Amp[IdxKa]);
	  if (SigSq[IdxKa]<0 ||
	      Pos[IdxKa]<xvals[LeftBound[IdxKa]] || 
	      Pos[IdxKa]>xvals[RightBound[IdxKa]]) 
	    printf("WARNING in T49DedxFunction::Init: error for refitted kaon peak. n");
	}
	// refit proton
	if (IdxPr >=0 && IdxPr != IdxKa) {
	  PolyFit(RightBound[IdxPr]-LeftBound[IdxPr]+1,&xvals[LeftBound[IdxPr]],&yvals[LeftBound[IdxPr]],&errs[LeftBound[IdxPr]],2,pars);
	  Pos[IdxPr] = -pars[1]/2/pars[2];
	  SigSq[IdxPr] = -0.5/pars[2]*(pars[0]-0.25*pars[1]*pars[1]/pars[2]);
	  Amp[IdxPr] = -2*pars[2]*SigSq[IdxPr];
	  if (Verbose) printf("Refitted proton (#%d), Pos: %f, Width: %f, Amp: %fn",IdxPr+1,Pos[IdxPr],SigSq[IdxPr],Amp[IdxPr]);
	  if (SigSq[IdxPr]<0 ||
	      Pos[IdxPr]<xvals[LeftBound[IdxPr]] || 
	      Pos[IdxPr]>xvals[RightBound[IdxPr]]) 
	    printf("WARNING in T49DedxFunction::Init: error for refitted proton peak. n");
	}
	if (IdxPr == 0 ||
	    fabs(Pos[IdxKa]-pos[2]*PiPos) < fabs(Pos[IdxKa]-pos[3]*PiPos)) {
	  KaAmp=Amp[IdxKa];
	  //-PiAmp*exp(-0.5*(pos[2]-1)*PiPos*(pos[2]-1)*PiPos/Res);
	  PrAmp=0.01*PiAmp;
	}
	else {
	  PrAmp=Amp[IdxPr];
	  //-PiAmp*exp(-0.5*(pos[3]-1)*PiPos*(pos[3]-1)*PiPos/Res);
	  KaAmp=0.01*PiAmp;
	}
      }
    }
  }
  if (sqrt(Res)>(1-pos[3]) && (PrAmp>0.5*PiAmp && PrAmp<1.5*PiAmp)) {
    // For severely overlapping peaks:
    Res*=0.5;
    KaAmp/=2;
    PrAmp/=2;
    PiAmp/=2;
  }
  Float_t MaxAmp=PiAmp;
  if (KaAmp>MaxAmp) MaxAmp=KaAmp;
  if (PrAmp>MaxAmp) MaxAmp=PrAmp;
  InitPositions(PiPos);
  SetResolution(sqrt(Res)/PiPos);
  SetAmplitude(0,0.05*MaxAmp,0,5*MaxAmp);
  SetAmplitude(1,PiAmp,0,5*MaxAmp);
  SetAmplitude(2,KaAmp,0,5*MaxAmp);
  SetAmplitude(3,PrAmp,0,5*MaxAmp);
  SetRange(0.5,2.0);
}

 Float_t T49DedxFunction::RelRise(Float_t Ptot, Int_t Part)
{
//  Returns expected dE/dx values, according to Bethe-Bloch parametrization
//  (Called by GetRelRise())
//  ** Code provided by Christof Roland **
 
  Float_t bg, dedx;
  Float_t c, d, e, f, x1,x2, p0, dfq, d1, d2, d3;
  Double_t mass[4] = {MASS_E_2, MASS_PI_2, MASS_K_2, MASS_P_2};
  // char *names [4] = {"elec","Pions", "Kaons", "Protons"};
  
  bg = Ptot / sqrt (mass[Part]);
  
  c =  fParaC;
  d =  fParaD;
  e =  fParaE;
  f =  fParaF;

  x1 = pow(10,(e - 1.0/3.0 * sqrt(2.0*log(10.0)/(3.0 * f))));
  x2 = pow(10,(e + 2.0/3.0 * sqrt(2.0*log(10.0)/(3.0 * f))));
  
  p0 = c/(d + 2.0*log(10.0)*e - 1.0);
  
  if (bg<x1)
    dfq = 0;
  else
    {
      dfq = -2.0 * log(bg) + 2.0 * log(10.0) * e;
      if(bg<x2)
	{
	  d1 = 2.0/3.0*sqrt(2.0*log(10.0)/(3.0 * f));
	  d2 = log(bg)/log(10.0);
	  d3 = pow(( e + d1 - d2),3);
	  dfq -= f*d3;
	}
    }
  
  dedx = p0*( (1+ bg*bg)/(bg*bg) * ( d + log(1+(bg*bg)) + dfq)-1.0 );
  
  return dedx;
}

 void T49DedxFunction::PolyFit(const Int_t NPoints, Float_t* x, Float_t* y, Float_t* err, Int_t Order, Float_t* Result)
{
//  Fits polynomial of form par[0] + par[1]*x +par[2]*x^2, etc
//  Uses the provided errors as weights.
//  (Called internally by Init())

   Float_t u[NPoints];
   Float_t v[NPoints];

   Float_t MinX=x[0];
   Float_t MaxX=x[0];
   Float_t MinY=y[0];
   Float_t MaxY=y[0];
   for (Int_t i=0; i<NPoints; i++) {
      if (x[i]>MaxX) MaxX = x[i];
      if (x[i]<MinX) MinX = x[i];
      if (y[i]>MaxY) MaxY = y[i];
      if (y[i]<MinY) MinY = y[i];
   }
   Float_t XScale=1/(MaxX-MinX);
   Float_t YScale=1/(MaxY-MinY);
   for (Int_t i=0; i<NPoints; i++) {
      u[i]=XScale*(x[i]-MinX);
      v[i]=YScale*(y[i]-MinY);
   }

   TMatrix upow(NPoints,2*Order+1);
   for (Int_t i=0; i<=2*Order; i++)
     for (Int_t j=0; j<NPoints; j++){
       if (i==0) upow(j,i)=1/(err[j]*err[j]*YScale*YScale);
          else upow(j,i)=upow(j,i-1)*u[j]; }

   TVector B(Order+1);
   for (Int_t i=0; i<=Order; i++) {
     B(i)=0;
     for (Int_t j=0; j<NPoints; j++){
       B(i)+=upow(j,i)*v[j];
     }
   }
   TMatrix A(Order+1,Order+1);
   for (Int_t i=0; i<=Order; i++)
     for (Int_t j=0; j<=Order; j++) {
       if (j>=i) {
          A(i,j)=0;
          for (Int_t k=0; k<NPoints; k++) {
             A(i,j) += upow(k,j+i);
          }; 
       } else A(i,j)=A(j,i);
     }
   TMatrix Ainv(TMatrix::kInverted,A);
   B *= Ainv;
   Ainv *= A;
   if (Ainv(Order,1)>1e-4) printf("ERROR: numerical accuracy too low in T49DedxFunction::PolyFitn");
   Float_t XScalePow=1.0;
   Float_t MinXPow=1.0;
   for (Int_t k=0;k<=Order;k++) {
      Result[k]=0;
      for (Int_t j=0;j<=Order;j++) {
        if (j==0)
           XScalePow=1.0;
        else XScalePow*=XScale;
        if (j>=k) {
          if (j-k==0) 
             MinXPow=1.0;
          else MinXPow*=-MinX;
          Result[k]+=B(j)*XScalePow*MinXPow*Binom(j,k);
        }
      }
      Result[k]=1/YScale*Result[k];
      if (k==0) Result[k]+=MinY;
   }
}

 Int_t T49DedxFunction::Binom(Int_t n,Int_t k)
{
//  Calculates binomial coefficient n over k.
//  Called by PolyFit()

  if (n<k) {
    printf("ERROR in Binom: n<k, n=%d k=%dn",n,k);
    return 0;
  }
  else {
    Int_t Num=1;
    for (Int_t i=n;i>k;i--) Num*=i;
    Int_t Denom=1;
    for (Int_t i=n-k;i>1;i--) Denom*=i;
    return Num/Denom;
  }
}



















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.