// Author: Gero Flucke <mailto:gero.flucke@desy.de>
//________________________________________
//
//   last update: $Date: 2005/11/26 20:04:56 $
//   by:          $Author: flucke $
//
// fitting class for D*-Analysis:
//                  fits delta-mass-plots (m_{Kpipis} - m_{Kpi})
//  (for use in ROOT 3.00_06 or higher only!) 
//
// performs the 'unsmart' fits for DstarDmFitter,
// so please use that class! 


#include "DstarDmFitterUnsmart.h"
#include "DstarDmFitterSpecial_I.h"

#include <TROOT.h>
#include <TError.h>
#include <TH1.h>
#include <TF1.h>
#include <TMath.h>
// #include <TVirtualFitter.h>

#include <iostream>
#include <iomanip>

using namespace std;

ClassImp(DstarDmFitterUnsmart)

//********************************************************
//******* constructor, destructor, set and get ***********
//********************************************************

 DstarDmFitterUnsmart::DstarDmFitterUnsmart(TH1* hist, TH1* hWc)
{
  // contructor: hist as fitted histogram (may not be a TH2/3*!)

  // initialise some data members:
  fBackGround = fFixMean = fFixSigma = fFixBackgrExp = fFixBackgrSqr = kFALSE;
  fFixWcScale = kFALSE;
  fLimitBackgrSqr = kFALSE;
  fLimitBackgrExp = kFALSE;
  fWcScale = 1.;

  this->SetUpFitLimit();

  this->SetUpperBackgrSqrt();
  this->SetLowerBackgrSqrt();
  this->SetUpperBackgrExp ();
  this->SetLowerBackgrExp ();

  fDstarFitFormu = new TF1("dstarFitFormuUnsm", FitFct, 0.135, fUpFitLimit,6);
  fDstarFitFormu->SetBit(TFormula::kNotGlobal);
  gROOT->GetListOfFunctions()->Remove(fDstarFitFormu);

  this->SetHist(hist, hWc, kTRUE); // default fit values are set!
}

//********************************************************

 DstarDmFitterUnsmart::DstarDmFitterUnsmart(const DstarDmFitterUnsmart& source)
{
  fHist = source.GetHistPtr();
  fHistWc = source.GetHistWcPtr();
  this->SetUpFitLimit(source.GetUpFitLimit());
  fIsMC = source.HistIsMC();
  this->SetDetectIsMC(source.fDetectIsMC);

  fFixMean = source.fFixMean;
  fFixSigma = source.fFixSigma;
  fFixBackgrExp = source.fFixBackgrExp;
  fFixBackgrSqr = source.fFixBackgrExp;
  fFixWcScale = source.fFixWcScale;
  fLimitBackgrSqr = source.fLimitBackgrSqr;
  fUpperBackgrSqr = source.fUpperBackgrSqr;
  fLowerBackgrSqr = source.fLowerBackgrSqr;
  fLimitBackgrExp = source.fLimitBackgrExp;
  fUpperBackgrExp = source.fUpperBackgrExp;
  fLowerBackgrExp = source.fLowerBackgrExp;

  fWcScale = source.fWcScale;

  // potential source for bug: 
  //     -which parameters are fixed?
  //     -name clash???????
  // :
  fDstarFitFormu = new TF1(*(source.fDstarFitFormu));
  gROOT->GetListOfFunctions()->Remove(fDstarFitFormu);
  fDstarFitFormu->SetBit(TFormula::kNotGlobal);
}

//********************************************************

 DstarDmFitterUnsmart::DstarDmFitterUnsmart(const DstarDmFitterSpecial_I& source)
{
  fHist = source.GetHistPtr();
  fHistWc = source.GetHistWcPtr();
  this->SetUpFitLimit(source.GetUpFitLimit());
  fIsMC = source.HistIsMC();
  this->SetDetectIsMC(source.IsDetectIsMC());

  fFitOption = source.GetFitOption();

  fFixMean = fFixSigma = fFixBackgrExp = fFixBackgrSqr = kFALSE;
  fFixWcScale = kFALSE;
  fLimitBackgrSqr = kFALSE;
  fLimitBackgrExp = kFALSE;
  fWcScale = 1.;

  this->SetUpperBackgrSqrt();
  this->SetLowerBackgrSqrt();
  this->SetUpperBackgrExp ();
  this->SetLowerBackgrExp ();
  
  fDstarFitFormu = new TF1("dstarFitFormuUnsm", FitFct, 0.135, fUpFitLimit,6);
  fDstarFitFormu->SetParameters(source.GetBackgrNorm(),	source.GetBackgrExp(),
				source.GetBackgrSqr(), source.GetMean(),
				source.GetSigma(), source.GetNDstar());
  fDstarFitFormu->SetParError(0, source.GetBackgrNormErr());
  fDstarFitFormu->SetParError(5, source.GetNDstarErr());
  gROOT->GetListOfFunctions()->Remove(fDstarFitFormu);
  fDstarFitFormu->SetBit(TFormula::kNotGlobal);
}

DstarDmFitterUnsmart& DstarDmFitterUnsmart::operator=
(const DstarDmFitterUnsmart& rhs)
{
  // get the same/aequivalent data member values from rhs
  if (this != &rhs) {
    fHist = rhs.GetHistPtr();
    fHistWc = rhs.GetHistWcPtr();
    this->SetUpFitLimit(rhs.GetUpFitLimit());
    fIsMC = rhs.HistIsMC();
    this->SetDetectIsMC(rhs.IsDetectIsMC());
    fFitOption = rhs.fFitOption;

    fFixMean = rhs.fFixMean;
    fFixSigma = rhs.fFixSigma;
    fFixBackgrExp = rhs.fFixBackgrExp;
    fFixBackgrSqr = rhs.fFixBackgrExp;
    fLimitBackgrExp = rhs.fLimitBackgrExp;
    fLimitBackgrSqr = rhs.fLimitBackgrSqr;
    fFixWcScale = rhs.fFixWcScale;

    fDstarFitFormu->SetParameter(0, rhs.GetBackgrNorm());
    fDstarFitFormu->SetParError(0, rhs.GetBackgrNormErr());

    this->SetBackgrExp(rhs.GetBackgrExp());
    this->SetBackgrSqr(rhs.GetBackgrSqr());
    this->SetMean(rhs.GetMean());
    this->SetSigma(rhs.GetSigma());

    this->SetUpperBackgrSqrt(rhs.GetUpperBackgrSqr());
    this->SetLowerBackgrSqrt(rhs.GetLowerBackgrSqr());
    this->SetUpperBackgrExp (rhs.GetUpperBackgrExp());
    this->SetLowerBackgrExp (rhs.GetLowerBackgrExp());

    fDstarFitFormu->SetParameter(5, rhs.GetNDstar());
    fDstarFitFormu->SetParError(5, rhs.GetNDstarErr());
  }
  return (*this);
}

DstarDmFitterUnsmart& DstarDmFitterUnsmart::operator=
(const DstarDmFitterSpecial_I& rhs)
{
  // get the same/aequivalent data member values from rhs
  // fit mode booleans are not set to false!
  if (this != &rhs) {
    fHist = rhs.GetHistPtr();
    fHistWc = rhs.GetHistWcPtr();
    this->SetUpFitLimit(rhs.GetUpFitLimit());
    fIsMC = rhs.HistIsMC();
    this->SetDetectIsMC(rhs.IsDetectIsMC());
    
    fFixMean = fFixSigma = fFixBackgrExp = fFixBackgrSqr = fFixWcScale = kFALSE;
    fLimitBackgrSqr = kFALSE;
    fLimitBackgrExp = kFALSE;

    fDstarFitFormu->SetParameter(0, rhs.GetBackgrNorm());
    fDstarFitFormu->SetParError(0, rhs.GetBackgrNormErr());

    this->SetBackgrExp(rhs.GetBackgrExp());
    this->SetBackgrSqr(rhs.GetBackgrSqr());
    this->SetMean(rhs.GetMean());
    this->SetSigma(rhs.GetSigma());

    this->SetUpperBackgrSqrt();
    this->SetLowerBackgrSqrt();
    this->SetUpperBackgrExp ();
    this->SetLowerBackgrExp ();

    fDstarFitFormu->SetParameter(5, rhs.GetNDstar());
    fDstarFitFormu->SetParError(5, rhs.GetNDstarErr());
  }
  return (*this);
}


/********************************************************/

 DstarDmFitterUnsmart::~DstarDmFitterUnsmart()
{
  //deletes all possibly created objects that members point at
  if(fDstarFitFormu)  {delete fDstarFitFormu; fDstarFitFormu = 0;}
}

/********************************************************/

 void DstarDmFitterUnsmart::SetHistReal(TH1* histo, TH1* hWc,
				Bool_t defaultSigmaMeanExpSqr)
{
  // by default defaultSigmaMeanExpSqr is kFALSE, if kTRUE, 
  // Sigma, Mean, BackgrExp and BackgrSqr are set to the default.

  fHist = histo;
  fHistWc = hWc;
  if(fHist){
    cout << "Fitted histogram is " << fHist->GetName();
    if(fHistWc) cout << ", WC is " << fHistWc->GetName(); 
    cout << endl;
  }

  if(defaultSigmaMeanExpSqr){ // sets defaults, 
    fDstarFitFormu->SetParameter(1, fgBackgrExpDefault);
    fDstarFitFormu->SetParameter(2, fgBackgrSqrDefault);
    fDstarFitFormu->SetParameter(3, fgMeanDefault);
    fDstarFitFormu->SetParameter(4, fgSigmaDefault);
  } // else last values are used for next fit
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetLimitBackground(Bool_t l)
{
  // switch limits for background parameters
  SetLimitBackgrExp (l);
  SetLimitBackgrSqrt(l);
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetLimitBackgrExp(Bool_t l)
{
  // switch limits for exponent of background
  fLimitBackgrExp = l;
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetLimitBackgrSqrt(Bool_t l)
{
  // switch limits for sqr-correction of background
  fLimitBackgrSqr = l;
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetUpperBackgrExp(Double_t e)
{
  // set upper limit for exponent of background
  fUpperBackgrExp = e;
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetLowerBackgrExp(Double_t e)
{
  // set lower limit for exponent of background
  fLowerBackgrExp = e;
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetUpperBackgrSqrt(Double_t s)
{
  // set upper limit for sqr-correction of background 
  fUpperBackgrSqr = s;
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetLowerBackgrSqrt(Double_t s)
{
  // set lower limit for sqr-correction of background 
  fLowerBackgrSqr = s;
}

//******************************************************
  
 void DstarDmFitterUnsmart::SetBackgrExp(Double_t e)
{
  // set exponent of background
  fDstarFitFormu->SetParameter(1, e);
}

//******************************************************

 void DstarDmFitterUnsmart::SetBackgrSqr(Double_t s)
{
  // set sqr-correction of background
  fDstarFitFormu->SetParameter(2, s);
}

//******************************************************

 void DstarDmFitterUnsmart::SetMean(Double_t m)
{
  // set mean of D*-Peak
  if(fBackGround) {
    fBackGround = kFALSE;
    this->MakeFitFunction(this->GetBackgrNorm(), 0.);
  }
  fDstarFitFormu->SetParameter(3, m);
}

//******************************************************

 void DstarDmFitterUnsmart::SetSigma(Double_t s)
{
  // set sigma of D*-peak
  if(fBackGround) {
    fBackGround = kFALSE;
    this->MakeFitFunction(this->GetBackgrNorm(), 0.);
  }
  fDstarFitFormu->SetParameter(4, s);
}

//******************************************************

 void DstarDmFitterUnsmart::SetWcScale(Double_t s)
{
  // set scale to be applied for Wc background
  fWcScale = s;
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetBackgrNorm() const
{
  // get normalisation factor of background
  return fDstarFitFormu->GetParameter(0);
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetBackgrNormErr() const
{
  // get sigma of normalisation factor of background
  return fDstarFitFormu->GetParError(0);
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetBackgrExp() const
{
  // get exponent of background
  return fDstarFitFormu->GetParameter(1);
}

//**********************************************

 Bool_t DstarDmFitterUnsmart::GetLimitBackgrExp() const
{
  // limit for exponent of backgr. switched on?
  return fLimitBackgrExp;
}

//**********************************************

 Double_t DstarDmFitterUnsmart::GetUpperBackgrExp() const
{
  // get upper limit for exponent of backgr.
  return fUpperBackgrExp;
}

//**********************************************

 Double_t DstarDmFitterUnsmart::GetLowerBackgrExp() const
{
  // get lower limit for exponent of backgr.
  return fLowerBackgrExp;
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetBackgrSqr() const
{
  // get sqr correction of background
  return fDstarFitFormu->GetParameter(2);
}

//**********************************************

 Bool_t DstarDmFitterUnsmart::GetLimitBackgrSqr() const
{
  // limit for sqrt correction of backgr. switched on?
  return fLimitBackgrSqr;
}

//**********************************************

 Double_t DstarDmFitterUnsmart::GetUpperBackgrSqr() const
{
  // get upper limit for sqrt correction of backgr.
  return fUpperBackgrSqr;
}

//**********************************************

 Double_t DstarDmFitterUnsmart::GetLowerBackgrSqr() const
{
  // get lower limit for sqrt correction of backgr.
  return fLowerBackgrSqr;
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetMean() const
{
  // get mean of D*-Peak
  return fBackGround ? fgMeanDefault : fDstarFitFormu->GetParameter(3);
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetSigma() const
{
  // get sigma of D*-peak
  return fBackGround ? fgSigmaDefault : fDstarFitFormu->GetParameter(4);
}

//**********************************************

 Double_t DstarDmFitterUnsmart::GetNDstar() const
{
  // get N(D*)
  return fBackGround ? 0. : fDstarFitFormu->GetParameter(5);
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetNDstarErr() const
{
  // get sigma(N(D*))
  return fBackGround ? 0. : fDstarFitFormu->GetParError(5);
}

//**********************************************

 Double_t DstarDmFitterUnsmart::GetWcScale() const
{
  // get U_n / U_n,wc
  return fWcScale;
}

//******************************************************

 Double_t DstarDmFitterUnsmart::GetSignalBack() const
{
  // Get ratio of signal over background of last fit
  // in (mean +- 2*sigma)-region.
  // Gives 0. if MC or background or any other strange result,
  // but doesn't check whether you have not fitted yet.
  if(fIsMC || !fDstarFitFormu || fBackGround) return 0.;

  Double_t signalRegionLow = this->GetMean() - 2.*this->GetSigma();
  Double_t signalRegionHigh = this->GetMean() + 2.*this->GetSigma();

  TF1 backGr("backGrUnsm", BackgrFct, 0.135, fUpFitLimit, 3);
  backGr.SetParameters(this->GetBackgrNorm(), this->GetBackgrExp(), 
		       this->GetBackgrSqr());

  fgBinSize = GetHistBinSize(); //needed to ensure 
  fgUpFitLimit = fUpFitLimit;   // dito

  Double_t backgrIntegr = backGr.Integral(signalRegionLow, signalRegionHigh);
  Double_t signalIntegr = fDstarFitFormu->Integral(signalRegionLow,
						   signalRegionHigh);
  return 
    backgrIntegr ? (signalIntegr - backgrIntegr) / backgrIntegr : 0.0;
}

//******************************************************

 TF1* DstarDmFitterUnsmart::GetFittedFunc() const 
{
  // pointer to the the fitted function associated with the histogram
  TH1* h = this->GetHistPtr();
  return (h ? h->GetFunction("dstarFitFormuUnsm") : NULL);
}

//******************************************************
 TF1* DstarDmFitterUnsmart::GetFittedBackgrFunc() const
{
  // pointer to the background part of the fitted function associated with the histogram
  TH1* h = this->GetHistPtr();
  return (h ? h->GetFunction("backGrUnsm") : NULL);
}

/*********************************************************/
/****************** the fit: *****************************/
/*********************************************************/

 Bool_t DstarDmFitterUnsmart::Fit(Int_t mode)
{
  // fits the 6(7) parameter fitting fct to the histogram
  // 
  // the last four digits of 'mode' determine which parameters
  // of the fitting function should be fixed:
  // If the digit == 0 ==> parameter is fitted
  //          else     ==> parameter is fixed to actual value
  // 1st digit:    Mean
  // 2nd :         Sigma
  // 3rd :         BackgrExp
  // 4th :         BackgrSqr
  // 5th :         WcScale (ignored if no WC hist)
  //
  // Digit are counted from the right, so 
  //    mode = 1011 means: fix Mean, Sigma and BackgrSqr, 
  //                       but fit BackgrExp (and WcScale)
  //                      (and fit BackgrNorm, NDstar)
  //
  // never use leading 0's: that will interpret mode as octal!
  // so mode = 10   means: fix only sigma, the rest is fitted
  // default mode is 0: fit all 6 parameters
  //
  // if mode < 0: fit background only:
  //       mode = -100 means: N(D*) == 0, fit background with fixed BackgrExp


  if(!fHist) {
    this->Error("Fit","No histogram! Fit canceled!");
    return kFALSE;
  }
  const Double_t backgrNormStart = this->DetermBackgrStartParam(); // also fIsMC set!
  const Double_t nDstarStart = this->DetermNDstarStart(backgrNormStart);
  this->DetermFitMode(mode); // which parameter should be fixed?

  const Int_t fitResult = this->DoTheFit(backgrNormStart, nDstarStart);
  if(fitResult) return kFALSE;
  // add background function to hist 
  if(!fBackGround) this->AddBackgrToHist();
  return kTRUE;
}

//********************************************************

 void DstarDmFitterUnsmart::AddBackgrToHist()
{
  TF1* backGr = new TF1("backGrUnsm", BackgrFct, 
			fgMassPion, fUpFitLimit, 3);
  gROOT->GetListOfFunctions()->Remove(backGr);
  backGr->SetBit(TFormula::kNotGlobal);

  backGr->SetParameters(this->GetBackgrNorm(), this->GetBackgrExp(), 
			this->GetBackgrSqr());
  backGr->SetParent(fHist);
  backGr->SetLineStyle(2);
  backGr->SetLineWidth(2);
  backGr->SetFillStyle(0); // no filling: Why do we need to set it here...? Since 3.04_02!
  TList* hFuncList = fHist->GetListOfFunctions();
  hFuncList->Add(backGr); // now fHist is responsible for deleting backGr
}


/********************************************************/
/******** private member functions: *********************/
/********************************************************/

 void DstarDmFitterUnsmart::DetermFitMode(Int_t mode)
{
  //interprets mode and saves the result in the foreseen Bool_t's
  if(mode < 0){ //background only!
    fBackGround = kTRUE;
    mode /= 10;
    fFixMean = fFixSigma = kFALSE; // for security...
  } else {
    fBackGround = kFALSE;
    fFixMean      =  mode        % 10;
    fFixSigma     = (mode /= 10) % 10;
    fFixWcScale   = (mode /1000) % 10;
  }
  fFixBackgrExp = (mode /= 10) % 10;
  fFixBackgrSqr = (mode /= 10) % 10;
}

//********************************************************

 Int_t DstarDmFitterUnsmart::DoTheFit(Double_t backgrNormStart, 
				    Double_t nDstarStart)
{
  // does the fit having found start values

  TH1* histToFit = this->MakeFitFunction(backgrNormStart, nDstarStart); // before set parlimits!

  if(fIsMC) fDstarFitFormu->SetParLimits(0,1,1);
  if(fFixBackgrExp || fIsMC) fDstarFitFormu->SetParLimits(1,1,1);
  if(fFixBackgrSqr || fIsMC) fDstarFitFormu->SetParLimits(2,1,1);
  if(fFixMean) fDstarFitFormu->SetParLimits(3,1,1);
  if(fFixSigma) fDstarFitFormu->SetParLimits(4,1,1);
  if(fFixWcScale) fDstarFitFormu->SetParLimits(6,1,1);

  if(fLimitBackgrExp && !fFixBackgrExp){
    if(fLowerBackgrExp <= fUpperBackgrExp){

      fDstarFitFormu->SetParLimits(1, fLowerBackgrExp, fUpperBackgrExp);
    } else{
      this->Error("DoTheFit", "Lower Bound %f exceeds upper bound %f, "
		  "limits for background exponent not set!", 
		  fLowerBackgrExp, fUpperBackgrExp);
    }
  }
  if(fLimitBackgrSqr && !fFixBackgrSqr){
    if(fLowerBackgrSqr <= fUpperBackgrSqr){
      fDstarFitFormu->SetParLimits(2, fLowerBackgrSqr, fUpperBackgrSqr);
    } else{
      this->Error("DoTheFit", "Lower Bound %f exceeds upper bound %f, " 
		  "limits for background sqrt not set!", 
		  fLowerBackgrSqr, fUpperBackgrSqr);
    }
  }

  fgBinSize = this->GetHistBinSize(); // important!!!!!!!!!!
  fgUpFitLimit = fUpFitLimit;   // dito
  if(fgBinSize <= 0.) {
    if(histToFit != fHist) delete histToFit;
    return -1;
  }
 
  Double_t oldParams[20]; // we have max 7 parameters...
  fDstarFitFormu->GetParameters(oldParams);
  const Int_t fitResult = this->RealFit(histToFit, fDstarFitFormu);
  if(fitResult != 0)  fDstarFitFormu->SetParameters(oldParams);

  if(histToFit != fHist){// e.g. WC hist given in addition
    if(!fFitOption.Contains("+")) {// remove old functions as TH1::Fit does, too:
      TIter next(fHist->GetListOfFunctions(), kIterBackward);
      while (TObject *obj = next()) {
	if (obj->InheritsFrom(TF1::Class())) delete fHist->GetListOfFunctions()->Remove(obj);
      }
    }
    if(!fitResult){    // add fitted func to histogram:
      TF1* fit = new TF1(*histToFit->GetFunction("dstarFitFormuUnsm"));
      gROOT->GetListOfFunctions()->Remove(fit);
      fit->SetBit(TFormula::kNotGlobal);

      fHist->GetListOfFunctions()->Add(fit);
      fit->SetParent(fHist);
      fit->SetRange(fgMassPion, fUpFitLimit); // to restrict drawn function 
      fWcScale = fit->GetParameter(6);
    }
    delete histToFit;
  } else {
    fWcScale = 1.;
  }
  return fitResult;
}

//********************************************************

 TH1* DstarDmFitterUnsmart::MakeFitFunction(Double_t backgrNormStart, 
				   Double_t nDstarStart)
{
  // makes fDstarFitFormu fit to environment...
  // and if background hist is given creates the corresponding hist that merges
  // the two hists - this hist has to be deleted afterwards!
  // (else the normal fHist is returned - do not delete that!)

  TF1* f = NULL;
  TH1* hist = NULL;
  if(fHistWc == NULL){
    if(fBackGround){
      f = new TF1("dstarFitFormuUnsm", BackgrFct, fgMassPion, fUpFitLimit, 3);
      f->SetParameters(backgrNormStart, this->GetBackgrExp(), this->GetBackgrSqr());
      f->SetParNames("U_{n}","U_{exp}","U_{sqr}");
    } else {
      f = new TF1("dstarFitFormuUnsm", FitFct, fgMassPion, fUpFitLimit, 6);
      f->SetParameters(backgrNormStart, this->GetBackgrExp(), this->GetBackgrSqr(),
		       this->GetMean(), this->GetSigma(), nDstarStart);
      f->SetParNames("U_{n}","U_{exp}","U_{sqr}","#mu","#sigma","N(D*)");
    }
    hist = fHist;
  } else {
    hist = this->CreateMergedHist(fHist, fHistWc, fgShiftWcHist); // fgShiftWcHist is set here!
    fgEndRightCharge = fHist->GetXaxis()->GetXmax();
    f = new TF1("dstarFitFormuUnsm", FitFctWithWc, fgMassPion, fgShiftWcHist+fUpFitLimit, 7);
    f->SetParameters(backgrNormStart, this->GetBackgrExp(), this->GetBackgrSqr(),
		     this->GetMean(), this->GetSigma(), nDstarStart, fWcScale);
    f->SetParNames("U_{n}","U_{exp}","U_{sqr}","#mu","#sigma","N(D*)", "U_{n}/U_{n,wc}");
  }


  delete fDstarFitFormu;
  fDstarFitFormu = f;
  gROOT->GetListOfFunctions()->Remove(fDstarFitFormu);
  fDstarFitFormu->SetBit(TFormula::kNotGlobal);

  return hist;
}

//********************************************************

 Double_t DstarDmFitterUnsmart::DetermNDstarStart(Double_t backgrNormStart)
{
  // determines start value of N(D*)
  Int_t binDsLow = fHist->FindBin(fgMeanDefault - 2.*fgSigmaDefault);
  Int_t binDsUp = fHist->FindBin(fgMeanDefault + 2.*fgSigmaDefault);
  Double_t integral = fHist->Integral(binDsLow, binDsUp);

  TF1 backGr("backGrUnsm", BackgrFct, fgMassPion, fUpFitLimit, 3);
  backGr.SetParNames("U_{n}", "U_{exp}", "U_{sqr}");
  backGr.SetParameters(backgrNormStart,
		       this->GetBackgrExp(),
		       this->GetBackgrSqr());
  // 		       fDstarFitFormu->GetParameter(1),
  // 		       fDstarFitFormu->GetParameter(2));
  // 		       fgBackgrExpDefault,
  // 		       fgBackgrSqrDefault);

  fgBinSize = this->GetHistBinSize(); // fgBinSize is used in BackgrFct!
  fgUpFitLimit = fUpFitLimit;         // dito

  Double_t backgrIntegr = 
    backGr.Integral(fgMeanDefault - 2.*fgSigmaDefault, 
		    fgMeanDefault + 2.*fgSigmaDefault) / fgBinSize;
  
  Double_t nDstarStart = integral-backgrIntegr > 0. ? 
    integral-backgrIntegr : 30.;
//   int oldPrecision = cout.precision(2);
  cout << "Start parameter for N(D*): " << nDstarStart
       << endl;
//   cout.precision(oldPrecision);
  return nDstarStart;
}


//********************************************************
//************ functions used for fits: ******************
//********************************************************

 Double_t DstarDmFitterUnsmart::BackgrFct(Double_t *x, Double_t *par)
{
  // background function:
  // par[0]: normalisation
  // par[1]: exponent
  // par[2]: square factor
  // x[0]:   deltaM

  if(x[0] <= fgMassPion) { // lower than mPion not possible for deltaM!
    TF1::RejectPoint(); // in case of 'WC-included' fit we are at point between both hists!
    return 0.;
  }

  // calculate normalisation factor (cf. log book 2003-04-05)
  // (integral from fgMassPion to upfitlimit == par[0])
  const Double_t relLimit = fgUpFitLimit - fgMassPion;
  const Double_t relPowP1Plus1 = TMath::Power(relLimit, par[1]+1.); 
  Bool_t print = kFALSE;
  Double_t norm = 0.;
  if(par[1] != -1. && par[1] != -2. && par[1] != -3.) {
    norm = (1.-par[2]*fgMassPion*fgMassPion)    *relPowP1Plus1 / (par[1]+1.)
           - par[2] * (relPowP1Plus1        *relLimit*relLimit / (par[1]+3.)
		       + 2.*fgMassPion* relPowP1Plus1*relLimit / (par[1]+2.)
		  );
    if(norm == 0.) {
      print = kTRUE;
      norm = DBL_MIN;
      ::Warning("DstarDmFitterUnsmart::BackgrFct", "Set norm = %g", norm); 
    }
  } else {
    // we might improve here by taking care of these special exponents...
    ::Warning("DstarDmFitterUnsmart::BackgrFct", "exponent = %.2f, reject point", par[1]); 
    TF1::RejectPoint();
    return 0.;
  }
  
  Double_t power = TMath::Power((x[0] - fgMassPion), par[1]);
  if(power == 0.) {
    print = kTRUE;
    power = DBL_MIN;
    static Long_t count = 0;
    if(count%100000 == 0)
      ::Warning("DstarDmFitterUnsmart::BackgrFct", "Set power = %g (%d)", power, count++); 
  }
  const Double_t result = (par[0] * fgBinSize / norm) * power * (1. - par[2] * x[0]*x[0]);
  if(fgDebug && print) {
    ::Info("DstarDmFitterUnsmart::BackgrFct", "result %.2g, norm %.2g, exp %.2f,"
	   "sqr %.2f, delta m %.2f, power %.2g, norm %.2g, binsize %.2g",
	   result, par[0], par[1], par[2], x[0], power, norm, fgBinSize);
  }
  
  return result;
}

//********************************************************

 Double_t DstarDmFitterUnsmart::FitFct(Double_t* x, Double_t* par)
{
  // fit function:
  //par[0]: normalisation of background
  //par[1]: exponent of background
  //par[2]: square factor of background
  //par[3]: gauss' mean
  //par[4]: gauss' sigma
  //par[5]: N(D*)
  //x[0]:   deltaM
  if(x[0] <= fgMassPion) { // lower than mPion not possible for deltaM!
    TF1::RejectPoint();
    return 0.;
  }
  if(par[4]==0.){ // no width ==> senseless!
    ::Warning("DstarDmFitterUnsmart::FitFct", "width = 0, reject point");
    TF1::RejectPoint();
    return 0.;
  }

  // fgBinSize is a static member of DstarDmFitterUnsmart 
  // and should be determined before every fit (&Integral)
  // to ensure that it has the value of the histogram's
  // binSize (and not that of another histogram that is 
  // fitted by another instance of DstarDmFitterUnsmart)!
  // [It has to be static because it is used 
  // in this static function DstarDmFitterUnsmart::FitFct(...)!]

//   if(x[0] >= 0.1475 && x[0] < 0.1495) {
//     TF1::RejectPoint();
//     return 0.;
//   }
  // DstarDmFitterSpecial_I::kSqrtOf2Pi = sqrt(2*pi) in Double_t-precision
  const Double_t dx = (x[0]-par[3]) / par[4];
  const Double_t gaus = 1./(DstarDmFitterSpecial_I::kSqrtOf2Pi*par[4]) * fgBinSize * par[5]
    * TMath::Exp(-dx*dx / 2.);
  return gaus + BackgrFct(x, par);
}

//********************************************************

 Double_t DstarDmFitterUnsmart::FitFctWithWc(Double_t* x, Double_t* par)
{
  // fit function:
  //par[0]: normalisation of background
  //par[1]: exponent of background
  //par[2]: square factor of background
  //par[3]: gauss' mean
  //par[4]: gauss' sigma
  //par[5]: N(D*)
  //par[6]: relative factor between par[0] and normalisation of WC background
  //x[0]:   deltaM

//   if(x[0] >= 0.1475 && x[0] < 0.1495) {
//     TF1::RejectPoint();
//     return 0.;
//   }
  if(x[0] > fgUpFitLimit){
    if(x[0] <= fgEndRightCharge) {
      TF1::RejectPoint();
      return 0.;
    } else {
      if(par[6] == 0.) {
	::Warning("DstarDmFitterUnsmart::FitFctWithWc", "U_{n}/U_{n,wc} = 0, reject point");
	TF1::RejectPoint();
	return 0.;
      }
      Double_t parWc[3] = {par[0]/par[6], par[1], par[2]};
      Double_t xWc = x[0] - fgShiftWcHist;
      return BackgrFct(&xWc, parWc);
    }
  } else {
    return FitFct(x, par);
  }
}