TF1 * funGausQuad=NULL;
TF1 * funSig=NULL;
TF1 * funQuad=NULL;
TF1 * funBkg=NULL;

TF1 * funDGausPol0=NULL;
TF1 * funPedestal = NULL;

Float_t dGausPol0Fun(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  return par[0] 
    + par[1]*(par[2]*TMath::Gaus(xx, 0, par[3], kTRUE)+(1-par[2])*TMath::Gaus(xx, 0, par[4], kTRUE));
}

Float_t gausFun(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  return par[0]*TMath::Gaus(xx, 0, par[1], kTRUE);
}

Float_t pedestalFun(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  return par[0]+par[1]*TMath::Gaus(xx, 0, par[2], kTRUE);
}


Float_t quadFun(Double_t *x, Double_t *par){
  if (x[0] > par[3] && x[0] < par[4]){
    TF1::RejectPoint();
    return 0;
  }
  return par[0] + par[1]*x[0] + par[2]*x[0]*x[0] ; //a+bx+cx**2      
}

void fitHistDEta(TH1F* h1, Double_t mean, Double_t res,
		 Double_t &yield, Double_t &yieldErr){

  Double_t fitRangeLow = mean - 8.*res;
  Double_t fitRangeHigh = mean + 8.*res;
  
  funDGausPol0 = new TF1("funDGausPol0", dGausPol0Fun, fitRangeLow, fitRangeHigh, 5);
  funDGausPol0->SetParameter(0, h1->GetBinContent(5));
  funDGausPol0->SetParameter(1, h1->Integral(h1->FindBin(mean-3.*res), h1->FindBin(mean+3.*res))*0.6);
  funDGausPol0->SetParameter(2, 0.9);
  funDGausPol0->SetParameter(3, res);
  funDGausPol0->SetParameter(4, 3.*res);
  
  h1->Fit(funDGausPol0, "ER+");

  Double_t f = funDGausPol0->GetParameter(2);
  
  funSig = new TF1("funSig", gausFun, fitRangeLow, fitRangeHigh, 2);
  funSig->SetLineStyle(2);

  funPedestal = new TF1("funPedestal", pedestalFun, fitRangeLow, fitRangeHigh, 3);
  funPedestal->SetLineStyle(3);
  funPedestal->SetParameter(0, funDGausPol0->GetParameter(0));
  funPedestal->SetParError(0, funDGausPol0->GetParError(0));

  if(f > 0.5){
    funSig->SetParameters(0, f*funDGausPol0->GetParameter(1));
    funSig->SetParameters(1, funDGausPol0->GetParameter(3));
    
    funPedestal->SetParameter(1, (1.-f)*funDGausPol0->GetParameter(1));
    funPedestal->SetParameter(2, funDGausPol0->GetParameter(4));
    funPedestal->SetParError(1, (1.-f)*funDGausPol0->GetParError(1));
    funPedestal->SetParError(2, funDGausPol0->GetParError(4));
  }
  else{
    funSig->SetParameters(0, (1.-f)*funDGausPol0->GetParameter(1));
    funSig->SetParameters(1, funDGausPol0->GetParameter(4));

    funPedestal->SetParameter(1, f*funDGausPol0->GetParameter(1));
    funPedestal->SetParameter(2, funDGausPol0->GetParameter(3));
    funPedestal->SetParError(1, f*funDGausPol0->GetParError(1));
    funPedestal->SetParError(2, funDGausPol0->GetParError(3));
  }

  //  funSig->Draw("same");

  Double_t sigma =  funSig->GetParameter(1);
  //Double_t C = funDGausPol0->GetParameter(0);

  Double_t bkg, sig, sigFit;
  Double_t bkgErr, sigErr, sigFitErr;
  Int_t binLow = h1->FindBin(-3.*sigma);
  Int_t binHigh = h1->FindBin(3.*sigma);
  Double_t binLowEdge = h1->GetXaxis()->GetBinLowEdge(binLow);
  Double_t binHighEdge = h1->GetXaxis()->GetBinUpEdge(binHigh);

  bkg = funPedestal->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1);
  bkgErr = funPedestal->IntegralError(binLowEdge, binHighEdge)/h1->GetBinWidth(1);
  sig = h1->IntegralAndError(binLow, binHigh, sigErr)-bkg;

  sigFit = funDGausPol0->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1)-bkg;

  Double_t sigFitTrue
    = funDGausPol0->GetParameter(1)/h1->GetBinWidth(1);
  

  cout<<"signal = "<<sig<<"; "
      <<"signal from fit = "<< sigFit<<"; "
      <<"signal from fit parameter ="<<sigFitTrue<<"; "
      <<"background = " << bkg <<endl;
  
  yield = sig; //sigFitTrue;
  //eff = ?
  //Double_t eff = 0.3; //should apply to get correct statistical uncertainties
  yieldErr = TMath::Sqrt(sigErr*sigErr+bkgErr*bkgErr); //TMath::Sqrt((sig+2.*bkg)*eff);

  delete funDGausPol0;
}


void fitHist(TH1F* h1, Double_t massMean, Double_t massRes, 
	     Double_t &yield, Double_t &yieldErr, 
	     Bool_t pol0 = kFALSE, Bool_t bPhi = kFALSE){
  Double_t fitRangeLow = massMean - 15.*massRes;
  Double_t fitRangeHigh = massMean + 15.*massRes;
  if(pol0){
    fitRangeLow = massMean - 8.*massRes;
    fitRangeHigh = massMean + 8.*massRes;
    if(bPhi){
      fitRangeLow = -0.5*TMath::Pi();
      fitRangeHigh = .5*TMath::Pi();
    }
  }
  funGausQuad = new TF1("gausQuad", "gaus(0)+[3]+[4]*x +[5]*x*x",
			fitRangeLow, fitRangeHigh);
  funGausQuad->SetLineWidth(1);
  funGausQuad->SetLineColor(kRed);
  
  funQuad = new TF1("quad", quadFun, fitRangeLow, fitRangeHigh, 5);
  funQuad->SetLineWidth(2);
  funQuad->SetLineColor(kBlue);
  funQuad->SetLineStyle(2);

  funBkg = new TF1("bkg", "[0]+[1]*x+[2]*x*x", fitRangeLow, fitRangeHigh);
  funBkg->SetLineWidth(2);
  funBkg->SetLineColor(kRed);
  funBkg->SetLineStyle(2);
  funBkg->SetFillStyle(3);
  funBkg->SetFillColor(kRed);

  funSig = new TF1("funSig", "gaus");

  funQuad->SetParameter(0, 0/*h1->GetBinContent(1)*/);
  if(pol0)
    funQuad->SetParameter(0, h1->GetBinContent(5));
  funQuad->SetParameter(1, 0);
  funQuad->SetParameter(2, 0);
  if(pol0){
    funQuad->FixParameter(1, 0);
    funQuad->FixParameter(2, 0);
  }
  funQuad->FixParameter(3, massMean-3.*massRes);
  funQuad->FixParameter(4, massMean+3.*massRes);

  Int_t fitStatus0 = h1->Fit(funQuad, "WLRIN");
  if(fitStatus0 != 0)
    fitStatus0 = h1->Fit(funQuad, "WLMRIN");
  //  h1->Fit(funQuad, "RIN");

  /*
  Double_t paramA = (h1->Integral(h1->FindBin(massMean-3.*massRes),
				  h1->FindBin(massMean+3.*massRes))
		     -funQuad->GetParameter(0)*6.*massRes/h1->GetBinWidth(1))
    /2.5/massRes;
  */
  funGausQuad->SetParameter(0, h1->GetEntries()/20.);

  if(pol0)
    funGausQuad->SetParameter(0, (h1->Integral(h1->FindBin(massMean-3.*massRes), h1->FindBin(massMean+3.*massRes))
			                            -funQuad->GetParameter(0)*6.*massRes/h1->GetBinWidth(1))*.3);
  funGausQuad->SetParameter(1, massMean);
  funGausQuad->SetParameter(2, massRes);
  funGausQuad->SetParameter(3, funQuad->GetParameter(0));
  funGausQuad->SetParameter(4, funQuad->GetParameter(1));
  funGausQuad->SetParameter(5, funQuad->GetParameter(2));
  if(pol0){
    funGausQuad->FixParameter(4, funQuad->GetParameter(1));
    funGausQuad->FixParameter(5, funQuad->GetParameter(2));
  }
  
  Int_t fitStatus1 =  h1->Fit(funGausQuad, "WLRI+");
  if(fitStatus1 != 0)
    fitStatus1 =  h1->Fit(funGausQuad, "WLMRI+");
  //h1->Fit(funGausQuad, "RI+");

  funBkg->SetParameter(0, funGausQuad->GetParameter(3));
  funBkg->SetParameter(1, funGausQuad->GetParameter(4));
  funBkg->SetParameter(2, funGausQuad->GetParameter(5));
  if(pol0){
    funBkg->FixParameter(1, funGausQuad->GetParameter(4));
    funBkg->FixParameter(2, funGausQuad->GetParameter(5));
  }

  funSig->SetParameter(0, funGausQuad->GetParameter(0));
  funSig->SetParameter(1, funGausQuad->GetParameter(1));
  funSig->SetParameter(2, funGausQuad->GetParameter(2));

  Double_t mean = funGausQuad->GetParameter(1);
  Double_t sigma =TMath::Abs(funGausQuad->GetParameter(2));

  cout<< "fitted mean = "<< mean <<", fitted sigma = " << sigma <<endl;

  Double_t bkg, sig, sigFit;
  Double_t sigErr;
  Int_t binLow = h1->FindBin(mean-3.*sigma);
  Int_t binHigh = h1->FindBin(mean+3.*sigma);
  Double_t binLowEdge = h1->GetXaxis()->GetBinLowEdge(binLow);
  Double_t binHighEdge = h1->GetXaxis()->GetBinUpEdge(binHigh);

  bkg = funBkg->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1);
  sig = h1->IntegralAndError(binLow, binHigh, sigErr)-bkg;
  
  sigFit = funSig->Integral(binLowEdge, binHighEdge)/h1->GetBinWidth(1);

  Double_t sigFitTrue 
    = funGausQuad->GetParameter(0)*TMath::Sqrt(2.*TMath::Pi())*sigma/h1->GetBinWidth(1);
  
  cout<<"signal = "<<sig<<"; "
      <<"signal from fit = "<< sigFit<<"; "
      <<"signal from fit parameter ="<<sigFitTrue<<"; "
      <<"background = " << bkg <<endl;


  yield = sig;
  //  if(bkg<=0) bkg = 0;  //???
  Double_t eff = (sig+bkg)/sigErr/sigErr;
  if(fitStatus0 != 0 || fitStatus1 != 0){
    binLow = h1->FindBin(massMean-3.*massRes);
    binHigh = h1->FindBin(massMean+3.*massRes);
    sig = h1->IntegralAndError(binLow, binHigh, sigErr);
    
    Double_t bkgErr1, bkgErr2;
    Double_t bkg1 = h1->IntegralAndError(h1->FindBin(massMean-8.*massRes), 
					                               h1->FindBin(massMean-5.*massRes),
					                               bkgErr1);
    Double_t bkg2 = h1->IntegralAndError(h1->FindBin(massMean+5.*massRes),
                                         h1->FindBin(massMean+8.*massRes),
                                         bkgErr2);
    yield = sig-bkg1-bkg2;
    
    yieldErr = TMath::Sqrt(sigErr*sigErr+(bkgErr1*bkgErr1+bkgErr2*bkgErr2));
  }
  else
    yieldErr = TMath::Sqrt(sigErr*sigErr+bkg/eff);//to get the correct statistics error when efficiency appplied

  delete funGausQuad;
  delete funQuad;
  delete funBkg;
  delete funSig;
}