#include "fitfunction.hh"


Double_t MyQuadroPeaks(Double_t* x, Double_t* par)
{
  //parameters are in this order:
  //general: sigma, tail_percentage, lambda, eta 
  //gauss1: area1, position1,
  //gauss2: area2, position2,
  //gauss3: area3, position3,
  //gauss3: area4, position4, sigma4
  //background: factor, offset
  Double_t peak_pos1 = x[0]-par[5];
  Double_t gauss1;
  if(par[0]!=0)
    gauss1 = par[4]/(par[0]*sqrt(2*(TMath::Pi()))) * TMath::Exp(-0.5*(peak_pos1*peak_pos1)/(par[0]*par[0]));
  else
    gauss1 = par[4] * TMath::Exp(-0.5*peak_pos1*peak_pos1);
  Double_t tail1;
  tail1 = par[1]*par[4] * TMath::Exp(par[2]*(x[0]-(par[5]-par[3]*2.35482005*par[0])));
  if(par[0] != 0)
    tail1 *= TMath::Erfc((x[0]-(par[5]-par[3]*2.35482005*par[0]))/par[0]);
  
  Double_t peak_pos2 = x[0]-par[7];
  Double_t gauss2;
  if(par[0]!=0)
    gauss2 = par[6]/(par[0]*sqrt(2*(TMath::Pi()))) * TMath::Exp(-0.5*(peak_pos2*peak_pos2)/(par[0]*par[0]));
  else
    gauss2 = par[6] * TMath::Exp(-0.5*peak_pos2*peak_pos2);
  Double_t tail2;
  tail2 = par[1]*par[6] * TMath::Exp(par[2]*(x[0]-(par[7]-par[3]*2.35482005*par[0])));
  if(par[0] != 0)
    tail2 *= TMath::Erfc((x[0]-(par[7]-par[3]*2.35482005*par[0]))/par[0]);

  Double_t peak_pos3 = x[0]-par[9];
  Double_t gauss3;
  if(par[0]!=0)
    gauss3 = par[8]/(par[0]*sqrt(2*(TMath::Pi()))) * TMath::Exp(-0.5*(peak_pos3*peak_pos3)/(par[0]*par[0]));
  else
    gauss3 = par[8] * TMath::Exp(-0.5*peak_pos3*peak_pos3);
  
  Double_t peak_pos4 = x[0]-par[11];
  Double_t gauss4;
  if(par[12]!=0)
    gauss4 = par[10]/(par[12]*sqrt(2*(TMath::Pi()))) * TMath::Exp(-0.5*(peak_pos4*peak_pos4)/(par[12]*par[12]));
  else
    gauss4 = par[10] * TMath::Exp(-0.5*peak_pos4*peak_pos4);
  
  Double_t linear = par[13]*x[0]+par[14];

  return (gauss1+tail1+gauss2+tail2+gauss3+gauss4+linear);
}

Double_t MyQuadroTails(Double_t* x, Double_t* par)
{
  //parameters are in this order:
  //general: sigma, tail_percentage, lambda, eta 
  //gauss1: area1, position1,
  //gauss2: area2, position2,
  //gauss3: area3, position3,
  //gauss3: area4, position4, sigma4
  //background: factor, offset
  Double_t tail1;
  tail1 = par[1]*par[4] * TMath::Exp(par[2]*(x[0]-(par[5]-par[3]*2.35482005*par[0])));
  if(par[0] != 0)
    tail1 *= TMath::Erfc((x[0]-(par[5]-par[3]*2.35482005*par[0]))/par[0]);

  Double_t tail2;
  tail2 = par[1]*par[6] * TMath::Exp(par[2]*(x[0]-(par[7]-par[3]*2.35482005*par[0])));
  if(par[0] != 0)
    tail2 *= TMath::Erfc((x[0]-(par[7]-par[3]*2.35482005*par[0]))/par[0]);

  return (tail1+tail2);
}

Double_t FitMyQuadroPeaks(TH1F* spec, Double_t* par, Double_t* par_err, 
			  Double_t range1, Double_t range2, Double_t range3, Double_t range4,
			  Double_t min_area1, Double_t min_area2, Double_t min_area3, Double_t min_area4,
			  Double_t fit_low, Double_t fit_high, Int_t noBG,
			  Int_t color, Int_t integrate, TCanvas* canvas)
{
  //parameters are in this order:
  //general: sigma, tail_percentage, lambda, eta 
  //gauss1: area1, position1,
  //gauss2: area2, position2,
  //gauss3: area3, position3,
  //gauss3: area4, position4, sigma4
  //background: factor, offset

  Int_t peak1_fit=false;
  Int_t peak2_fit=false;
  Int_t peak3_fit=false;
  Int_t peak4_fit=false;

  if(fit_low == -1)
    fit_low = spec->GetXaxis()->GetXmin();
  if(fit_high == -1)
    fit_high = spec->GetXaxis()->GetXmax();

  Double_t range_low, range_high;

  Int_t npar = 15;
  Int_t i;

  char fit_option[10];
  strncpy(fit_option, "RBN", 9);
  if(color <= 1)
    strncat(fit_option, "O", 2);
  if(integrate == 1)
    strncat(fit_option, "I", 2);

  //needs at least the positions of the three peaks
  if(par[5] == 0 || par[7] == 0 || par[9] == 0 || par[11] == 0)
    {
      cerr<<"Parameters are not set!"<<endl;
      return -1;
    }

  //set other parameters to default values if they are 0
  if(par[0] == 0)
    par[0] = 1;
  if(par[1] == 0)
    par[1] = 0.1;
  if(par[2] == 0)
    par[2] = 1;
  if(par[3] == 0)
    par[3] = 1;
  if(par[4] == 0)
    par[4] = spec->GetBinContent(spec->FindBin(par[5]));
  if(par[6] == 0)
    par[6] = spec->GetBinContent(spec->FindBin(par[7]));
  if(par[8] == 0)
    par[8] = spec->GetBinContent(spec->FindBin(par[9]));
  if(par[10] == 0)
    par[10] = spec->GetBinContent(spec->FindBin(par[11]));
  if(par[12] == 0)
    par[12] = 1;
  if(par[13] == 0)
    par[13] = 0;
  if(par[14] == 0)
    par[14] = 0;

  if(noBG == 1)
    {
      par[13] = 0;
      par[14] = 0;
    }

  TF1* fitfun;
  fitfun = new TF1("fitfun", MyQuadroPeaks, fit_low, fit_high, npar);
  fitfun->SetParNames("sigma", "tail_p", "lambda", "eta", "area1", "pos1", "area2", "pos2", "area3", "pos3");
  fitfun->SetLineColor(color);
  fitfun->SetLineWidth(2);

  // fit highest peak
  if(par[10]>par[8] && par[10]>par[6] && par[10]>par[4])
    {
      for(i=0;i<npar;i++)
	{
	  if(10<=i && i<=12)
	    {
	      fitfun->ReleaseParameter(i);
	      fitfun->SetParameter(i, par[i]);
	    }
	  else
	    fitfun->FixParameter(i, par[i]);
	}
      if(range4 == -1)
	{
	  range_low = par[11]-(par[11]-par[9])/2.;
	  range_high = par[11]+(par[11]-par[9])/2.;
	}
      else
	{
	  range_low = (par[11]-(par[11]-par[9])/2 > par[11]-range4/2.)?
	              (par[11]-(par[11]-par[9])/2):(par[11]-range4/2.);
	  range_high = par[11]+range4/2.;
	}
      cout << "subfitrange: " << range_low << " to " << range_high << endl;
      fitfun->SetRange(range_low, range_high);

      fitfun->SetParLimits(10, min_area4, MAX_HEIGHT);
      fitfun->SetParLimits(11, par[11]-par[9], par[11]+par[9]);
      fitfun->SetParLimits(12, 0., 2.);
      
      spec->Fit("fitfun", fit_option);
      
      cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
      
      fitfun->FixParameter(10, fitfun->GetParameter(10));
      fitfun->FixParameter(11, fitfun->GetParameter(11));
      fitfun->FixParameter(12, fitfun->GetParameter(12));

      if(par != NULL)
	{
	  fitfun->GetParameters(par);
	}
      
      peak4_fit = true;
    }
  else if(par[8]>par[6] && par[8]>par[4]) 
    {
      for(i=0;i<npar;i++)
	{
	  if(i == 0 || (8<=i && i<=9))
	    {
	      fitfun->ReleaseParameter(i);
	      fitfun->SetParameter(i, par[i]);
	    }
	  else
	    fitfun->FixParameter(i, par[i]);
	}
      if(range3 == -1)
	{
	  range_low = par[9]-(par[9]-par[7])/2.;
	  range_high = par[9]+(par[11]-par[9])/2.;
	}
      else
	{
	  range_low = (par[9]-(par[9]-par[7])/2 > par[9]-range3/2.)?
	              (par[9]-(par[9]-par[7])/2):(par[9]-range3/2.);
	  range_high = (par[9]+(par[11]-par[9])/2 > par[9]+range3/2.)?
	              (par[9]+(par[11]-par[9])/2):(par[9]+range3/2.);
	}
      cout << "subfitrange: " << range_low << " to " << range_high << endl;
      fitfun->SetRange(range_low, range_high);
      
      fitfun->SetParLimits(0, 0.3, 1.5);
      fitfun->SetParLimits(8, min_area3, MAX_HEIGHT);
      fitfun->SetParLimits(9, range_low, range_high);
      
      spec->Fit("fitfun", fit_option);
      
      cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
      
      fitfun->FixParameter(0, fitfun->GetParameter(0));
      fitfun->FixParameter(8, fitfun->GetParameter(8));
      fitfun->FixParameter(9, fitfun->GetParameter(9));
      
      if(par != NULL)
	{
	  fitfun->GetParameters(par);
	}
      
      peak3_fit = true;
    } 
  else if(par[6]>par[4]) 
    {
      for(i=0;i<npar;i++)
	{
	  if(i == 0 || (6<=i && i<=7))
	    {
	      fitfun->ReleaseParameter(i);
	      fitfun->SetParameter(i, par[i]);
	    }
	  else
	    fitfun->FixParameter(i, par[i]);
	}
      if(range2 == -1)
	{
	  range_low = par[7]-(par[7]-par[5])/2.;
	  range_high = par[7]+(par[9]-par[7])/2.;
	}
      else
	{
	  range_low = (par[7]-(par[7]-par[5])/2 > par[7]-range2/2.)?
	              (par[7]-(par[7]-par[5])/2):(par[7]-range2/2.);
	  range_high = (par[7]+(par[9]-par[7])/2 > par[7]+range2/2.)?
	              (par[7]+(par[9]-par[7])/2):(par[7]+range2/2.);
	}
      cout << "subfitrange: " << range_low << " to " << range_high << endl;
      fitfun->SetRange(range_low, range_high);
      
      fitfun->SetParLimits(0, 0.3, 1.5);
      fitfun->SetParLimits(1, 0, 0.3);
      fitfun->SetParLimits(2, 0.1, 2);
      fitfun->SetParLimits(3, 0., 1.);
      fitfun->SetParLimits(6, min_area2, MAX_HEIGHT);
      fitfun->SetParLimits(7, range_low, range_high);
      
      spec->Fit("fitfun", fit_option);
      
      cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
      
      fitfun->FixParameter(0, fitfun->GetParameter(0));
      fitfun->FixParameter(1, fitfun->GetParameter(1));
      fitfun->FixParameter(2, fitfun->GetParameter(2));
      fitfun->FixParameter(3, fitfun->GetParameter(3));
      fitfun->FixParameter(6, fitfun->GetParameter(6));
      fitfun->FixParameter(7, fitfun->GetParameter(7));
      
      if(par != NULL)
	{
	  fitfun->GetParameters(par);
	}
      
      peak2_fit = true;
    } 
  else 
    {
      for(i=0;i<npar;i++)
	{
	  if(i<=5)
	    {
	      fitfun->ReleaseParameter(i);
	      fitfun->SetParameter(i, par[i]);
	    }
	  else
	    fitfun->FixParameter(i, par[i]);
	}
      if(range1 == -1)
	{
	  range_low = par[5]-(par[7]-par[5])/2.;
	  range_high = par[5]+(par[7]-par[5])/2.;
	}
      else
	{
	  range_low = (par[5]-(par[7]-par[5])/2. > par[5]-range1/2.)?
	              (par[5]-(par[7]-par[5])/2.):(par[5]-range1/2.);
	  range_high = (par[5]+(par[7]-par[5])/2. > par[5]+range1/2.)?
	              (par[5]+(par[7]-par[5])/2.):(par[5]+range1/2.);
	}
      cout << "subfitrange: " << range_low << " to " << range_high << endl;
      fitfun->SetRange(range_low, range_high);

      fitfun->SetParLimits(0, 0.3, 1.5);
      fitfun->SetParLimits(1, 0, 0.3);
      fitfun->SetParLimits(2, 0.1, 2);
      fitfun->SetParLimits(3, 0., 1.);
      fitfun->SetParLimits(4, min_area1, MAX_HEIGHT);
      fitfun->SetParLimits(5, range_low, range_high);
      
      spec->Fit("fitfun", fit_option);
      
      cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
      
      for(i=0;i<=5;i++)
	fitfun->FixParameter(i, fitfun->GetParameter(i));
      
      if(par != NULL)
	{
	  fitfun->GetParameters(par);
	}
      
      peak1_fit = true;
    }

  // fit second and third peak
  while( !peak1_fit || !peak2_fit || !peak3_fit || !peak4_fit) 
    {
      if( (peak4_fit && par[8]>=par[4] && par[8]>=par[4] && !peak3_fit) ||
	  (peak2_fit && par[8]>=par[10] && par[8]>=par[4] && !peak3_fit) ||
	  (peak1_fit && par[8]>=par[10] && par[8]>=par[6] && !peak3_fit) ||
	  (peak4_fit && peak2_fit && par[8]>=par[4] && !peak3_fit) ||
	  (peak4_fit && peak1_fit && par[8]>=par[6] && !peak3_fit) ||
	  (peak2_fit && peak1_fit && par[8]>=par[10] && !peak3_fit) ||
	  peak4_fit && peak2_fit && peak1_fit)
	{
	  if(!peak1_fit && !peak2_fit)
	    {
	      fitfun->ReleaseParameter(0);
	      fitfun->SetParameter(0, par[0]);
	    }
	  else
	    {
	      fitfun->FixParameter(0, par[0]);
	    }
	  for(i=1;i<npar;i++)
	    {
	      if(8<=i && i<=9)
		{
		  fitfun->ReleaseParameter(i);
		  fitfun->SetParameter(i, par[i]);
		}
	      else
		fitfun->FixParameter(i, par[i]);
	    }

	  if(range3 == -1)
	    {
	      range_low = par[9]-(par[9]-par[7])/2.;
	      range_high = par[9]+(par[11]-par[9])/2.;
	    }
	  else
	    {
	      range_low = (par[9]-(par[9]-par[7])/2 > par[9]-range3/2.)?
		          (par[9]-(par[9]-par[7])/2):(par[9]-range3/2.);
	      range_high = (par[9]+(par[11]-par[9])/2 > par[9]+range3/2.)?
		           (par[9]+(par[11]-par[9])/2):(par[9]+range3/2.);
	    }
	  cout << "subfitrange: " << range_low << " to " << range_high << endl;
	  fitfun->SetRange(range_low, range_high);
	  
	  if(!peak1_fit && !peak2_fit)
	    fitfun->SetParLimits(0, 0.3, 1.5);
	  fitfun->SetParLimits(8, min_area3, MAX_HEIGHT);
	  fitfun->SetParLimits(9, range_low, range_high);
	  
	  spec->Fit("fitfun", fit_option);
	  
	  cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
	  
	  if(!peak1_fit && !peak2_fit)
	    {
	      fitfun->FixParameter(0, fitfun->GetParameter(0));
	    }
	  fitfun->FixParameter(8, fitfun->GetParameter(8));
	  fitfun->FixParameter(9, fitfun->GetParameter(9));
	  
	  if(par != NULL)
	    {
	      fitfun->GetParameters(par);
	    }
	  
	  peak3_fit = true;
	} 
      else if( (peak4_fit && par[6]>=par[8] && par[6]>=par[4] && !peak2_fit) ||
	       (peak3_fit && par[6]>=par[10] && par[6]>=par[4] && !peak2_fit) ||
	       (peak1_fit && par[6]>=par[10] && par[6]>=par[8] && !peak2_fit) ||
	       (peak4_fit && peak3_fit && par[6]>=par[4] && !peak2_fit) ||
	       (peak4_fit && peak1_fit && par[6]>=par[8] && !peak2_fit) ||
	       (peak3_fit && peak1_fit && par[6]>=par[10] && !peak2_fit) ||
	       peak4_fit && peak3_fit && peak1_fit)
	{
	  if(!peak1_fit && !peak3_fit)
	    {
	      fitfun->ReleaseParameter(0);
	      fitfun->SetParameter(0, par[0]);
	    }
	  else
	    {
	      fitfun->FixParameter(0, par[0]);
	    }
	  if(!peak1_fit)
	    {
	      fitfun->ReleaseParameter(1);
	      fitfun->ReleaseParameter(2);
	      fitfun->ReleaseParameter(3);
	      fitfun->SetParameter(1, par[1]);
	      fitfun->SetParameter(2, par[2]);
	      fitfun->SetParameter(3, par[3]);
	    }
	  else
	    {
	      fitfun->FixParameter(1, par[1]);
	      fitfun->FixParameter(2, par[2]);
	      fitfun->FixParameter(3, par[3]);
	    }
	  for(i=4;i<npar;i++)
	    {
	      if(6<=i && i<=7)
		{
		  fitfun->ReleaseParameter(i);
		  fitfun->SetParameter(i, par[i]);
		}
	      else
		fitfun->FixParameter(i, par[i]);
	    }

	  if(range2 == -1)
	    {
	      range_low = par[7]-(par[7]-par[5])/2.;
	      range_high = par[7]+(par[9]-par[7])/2.;
	    }
	  else
	    {
	      range_low = (par[7]-(par[7]-par[5])/2 > par[7]-range2/2.)?
		          (par[7]-(par[7]-par[5])/2):(par[7]-range2/2.);
	      range_high = (par[7]+(par[9]-par[7])/2 > par[7]+range2/2.)?
		           (par[7]+(par[9]-par[7])/2):(par[7]+range2/2.);
	    }
	  cout << "subfitrange: " << range_low << " to " << range_high << endl;
	  fitfun->SetRange(range_low, range_high);
	  
	  if(!peak1_fit && !peak3_fit)
	    fitfun->SetParLimits(0, 0.3, 1.5);
	  if(!peak1_fit)
	    {
	      fitfun->SetParLimits(1, 0, 0.3);
	      fitfun->SetParLimits(2, 0.1, 2);
	      fitfun->SetParLimits(3, 0., 1.);
	    }
	  
	  fitfun->SetParLimits(6, min_area2, MAX_HEIGHT);
	  fitfun->SetParLimits(7, range_low, range_high);
	  
	  spec->Fit("fitfun", fit_option);
	  
	  cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
	  
	  if(!peak1_fit && !peak3_fit)
	    {
	      fitfun->FixParameter(0, fitfun->GetParameter(0));
	    }
	  if(!peak1_fit)
	    {
	      fitfun->FixParameter(1, fitfun->GetParameter(1));
	      fitfun->FixParameter(2, fitfun->GetParameter(2));
	      fitfun->FixParameter(3, fitfun->GetParameter(3));
	    }
	  fitfun->FixParameter(6, fitfun->GetParameter(6));
	  fitfun->FixParameter(7, fitfun->GetParameter(7));
	  
	  if(par != NULL)
	    {
	      fitfun->GetParameters(par);
	    }
	  
	  peak2_fit = true;
	} 
      else if( (peak4_fit && par[4]>=par[8] && par[4]>=par[6] && !peak1_fit) ||
	       (peak3_fit && par[4]>=par[10] && par[4]>=par[6] && !peak1_fit) ||
	       (peak2_fit && par[4]>=par[10] && par[4]>=par[8] && !peak1_fit) ||
	       (peak4_fit && peak3_fit && par[4]>=par[6] && !peak1_fit) ||
	       (peak4_fit && peak2_fit && par[4]>=par[8] && !peak1_fit) ||
	       (peak3_fit && peak2_fit && par[4]>=par[10] && !peak1_fit) ||
	       peak4_fit && peak3_fit && peak2_fit)
	{
	  if(!peak2_fit && !peak3_fit)
	    {
	      fitfun->ReleaseParameter(0);
	      fitfun->SetParameter(0, par[0]);
	    }
	  else
	    {
	      fitfun->FixParameter(0, par[0]);
	    }
	  if(!peak2_fit)
	    {
	      fitfun->ReleaseParameter(1);
	      fitfun->ReleaseParameter(2);
	      fitfun->ReleaseParameter(3);
	      fitfun->SetParameter(1, par[1]);
	      fitfun->SetParameter(2, par[2]);
	      fitfun->SetParameter(3, par[3]);
	    }
	  else
	    {
	      fitfun->FixParameter(1, par[1]);
	      fitfun->FixParameter(2, par[2]);
	      fitfun->FixParameter(3, par[3]);
	    }
	  for(i=4;i<npar;i++)
	    {
	      if(4<=i && i<=5)
		{
		  fitfun->ReleaseParameter(i);
		  fitfun->SetParameter(i, par[i]);
		}
	      else
		fitfun->FixParameter(i, par[i]);
	    }

	  if(range1 == -1)
	    {
	      range_low = par[5]-(par[4]-par[5])/2.;
	      range_high = par[5]+(par[4]-par[5])/2.;
	    }
	  else
	    {
	      range_low = (par[5]-(par[7]-par[5])/2 > par[5]-range1/2.)?
		          (par[5]-(par[7]-par[5])/2):(par[5]-range1/2.);
	      range_high = (par[5]+(par[7]-par[5])/2 > par[5]+range1/2.)?
		           (par[5]+(par[7]-par[5])/2):(par[5]+range1/2.);
	    }
	  cout << "subfitrange: " << range_low << " to " << range_high << endl;
	  fitfun->SetRange(range_low, range_high);
	  
	  if(!peak2_fit && !peak3_fit)
	    fitfun->SetParLimits(0, 0.3, 1.5);
	  if(!peak2_fit)
	    {
	      fitfun->SetParLimits(1, 0, 0.3);
	      fitfun->SetParLimits(2, 0.1, 2);
	      fitfun->SetParLimits(3, 0., 1.);
	    }

	  fitfun->SetParLimits(4, min_area1, MAX_HEIGHT);
	  fitfun->SetParLimits(5, range_low, range_high);
	  
	  spec->Fit("fitfun", fit_option);
	  
	  cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
	  
	  if(!peak2_fit && !peak3_fit)
	    {
	      fitfun->FixParameter(0, fitfun->GetParameter(0));
	    }
	  if(!peak2_fit)
	    {
	      fitfun->FixParameter(1, fitfun->GetParameter(1));
	      fitfun->FixParameter(2, fitfun->GetParameter(2));
	      fitfun->FixParameter(3, fitfun->GetParameter(3));
	    }
	  fitfun->FixParameter(4, fitfun->GetParameter(4));
	  fitfun->FixParameter(5, fitfun->GetParameter(5));
	  
	  if(par != NULL)
	    {
	      fitfun->GetParameters(par);
	    }
	  
	  peak1_fit = true;
	} 
      else 
	{
	  for(i=0;i<npar;i++)
	    {
	      if(10<=i && i<=12)
		{
		  fitfun->ReleaseParameter(i);
		  fitfun->SetParameter(i, par[i]);
		}
	      else
		fitfun->FixParameter(i, par[i]);
	    }

	  if(range4 == -1)
	    {
	      range_low = par[11]-(par[11]-par[9])/2.;
	      range_high = par[11]+(par[11]-par[9])/2.;
	    }
	  else
	    {
	      range_low = (par[11]-(par[11]-par[9])/2 > par[11]-range4/2.)?
		          (par[11]-(par[11]-par[9])/2):(par[11]-range4/2.);
	      range_high = par[11]+range4/2.;
	    }
	  cout << "subfitrange: " << range_low << " to " << range_high << endl;
	  fitfun->SetRange(range_low, range_high);
	  
	  fitfun->SetParLimits(10, min_area4, MAX_HEIGHT);
	  fitfun->SetParLimits(11, range_low, range_high);
	  fitfun->SetParLimits(12, 0.3, 2.);
	  
	  spec->Fit("fitfun", fit_option);
	  
	  cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
	  
	  fitfun->FixParameter(10, fitfun->GetParameter(10));
	  fitfun->FixParameter(11, fitfun->GetParameter(11));
	  fitfun->FixParameter(12, fitfun->GetParameter(12));
	  
	  if(par != NULL)
	    {
	      fitfun->GetParameters(par);
	    }
	  
	  peak4_fit = true;
	}
    }
  
  //fit background
  fitfun->SetRange(fit_low, fit_high);
  if(noBG == 0)
    {
      fitfun->ReleaseParameter(13);
      fitfun->ReleaseParameter(14);
      fitfun->SetParLimits(13, -1, 1);
      fitfun->SetParLimits(14, -10, 10);
      
      cout << "fitrange: " << fit_low << " to " << fit_high << endl;

      spec->Fit("fitfun", fit_option);
      
      cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;
    }

  //fit all
  for(i = 0; i < npar; i++)
    fitfun->ReleaseParameter(i);
  fitfun->SetParLimits(0, 0.3, 1.);
  fitfun->SetParLimits(1, 0, 0.3);
  fitfun->SetParLimits(2, 0.1, 2);
  fitfun->SetParLimits(3, 0., 1.);
  fitfun->SetParLimits(4, min_area1, MAX_HEIGHT);
  fitfun->SetParLimits(5, par[5]-par[0], par[5]-par[0]);
  fitfun->SetParLimits(6, min_area2, MAX_HEIGHT);
  fitfun->SetParLimits(7, par[7]-par[0], par[7]-par[0]);
  fitfun->SetParLimits(8, min_area3, MAX_HEIGHT);
  fitfun->SetParLimits(9, par[9]-par[0], par[9]-par[0]);
  fitfun->SetParLimits(10, min_area4, MAX_HEIGHT);
  fitfun->SetParLimits(11, par[11]-par[12], par[11]+par[12]);
  fitfun->SetParLimits(12, 0.3, 2.);
  if(noBG == 0)
    {
      fitfun->SetParLimits(13, -1, 1);
      fitfun->SetParLimits(14, -10, 10);
    }
  else
    {
      fitfun->FixParameter(13, 0);
      fitfun->FixParameter(14, 0);
    }

  cout << "fitrange: " << fit_low << " to " << fit_high << endl;

  strncat(fit_option, "EM+", 3);
  spec->Fit("fitfun", fit_option);

  cout << "Reduced chisquare: " << fitfun->GetChisquare()/fitfun->GetNDF() << endl;

  if(par != NULL)
    {
      fitfun->GetParameters(par);
    }
  
  if(par_err != NULL)
    {
      for(Int_t i=0;i<npar;i++)
	{
	  par_err[i] = fitfun->GetParError(i);
	}
    }
  
  if(color > 1)
    {
      TCanvas* MyCanvas;
      MyCanvas = (TCanvas*) new TCanvas("MyCanvas", "MyCanvas", 1);
      spec->Draw();
      fit_low = spec->GetXaxis()->GetXmin();
      fit_high = spec->GetXaxis()->GetXmax();
      TF1* function;
      function = new TF1("function", MyQuadroPeaks, fit_low, fit_high, npar);
      function->SetParameters(par);
      function->SetLineColor(2);
      TF1* gauss1;
      gauss1 = new TF1("gauss1", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss1->SetParameters(par);
      gauss1->SetParameter(1,0);
      gauss1->SetParameter(6,0);
      gauss1->SetParameter(8,0);
      gauss1->SetParameter(10,0);
      gauss1->SetLineColor(color+1);
      gauss1->Draw("CSAME");
      TF1* gauss2;
      gauss2 = new TF1("gauss2", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss2->SetParameters(par);
      gauss2->SetParameter(1,0);
      gauss2->SetParameter(4,0);
      gauss2->SetParameter(8,0);
      gauss2->SetParameter(10,0);
      gauss2->SetLineColor(color+1);
      gauss2->Draw("CSAME");
      TF1* gauss3;
      gauss3 = new TF1("gauss3", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss3->SetParameters(par);
      gauss3->SetParameter(1,0);
      gauss3->SetParameter(4,0);
      gauss3->SetParameter(6,0);
      gauss3->SetParameter(10,0);
      gauss3->SetLineColor(color+1);
      gauss3->Draw("CSAME");
      TF1* gauss4;
      gauss4 = new TF1("gauss4", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss4->SetParameters(par);
      gauss4->SetParameter(1,0);
      gauss4->SetParameter(4,0);
      gauss4->SetParameter(6,0);
      gauss4->SetParameter(8,0);
      gauss4->SetLineColor(color+1);
      gauss4->Draw("CSAME");
      TF1* tail1;
      tail1 = new TF1("tail1", MyQuadroTails, fit_low, fit_high, npar);
      tail1->SetParameters(par);
      tail1->SetParameter(6,0);
      tail1->SetParameter(8,0);
      tail1->SetParameter(10,0);
      tail1->SetLineColor(color+2);
      tail1->Draw("CSAME");
      TF1* tail2;
      tail2 = new TF1("tail2", MyQuadroTails, fit_low, fit_high, npar);
      tail2->SetParameters(par);
      tail2->SetParameter(4,0);
      tail2->SetParameter(8,0);
      tail2->SetParameter(10,0);
      tail2->SetLineColor(color+2);
      tail2->Draw("CSAME");
      function->Draw("CSAME");
    }
  if(color == -1 && canvas != NULL)
    {
      canvas->cd();

      spec->SetTitle("");
      spec->Draw();
      fit_low = spec->GetXaxis()->GetXmin();
      fit_high = spec->GetXaxis()->GetXmax();
      TF1* function;
      function = new TF1("function", MyQuadroPeaks, fit_low, fit_high, npar);
      function->SetParameters(par);
      function->SetLineColor(2);
      TF1* gauss1;
      gauss1 = new TF1("gauss1", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss1->SetParameters(par);
      gauss1->SetParameter(1,0);
      gauss1->SetParameter(6,0);
      gauss1->SetParameter(8,0);
      gauss1->SetParameter(10,0);
      gauss1->SetLineColor(4);
      gauss1->Draw("CSAME");
      TF1* gauss2;
      gauss2 = new TF1("gauss2", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss2->SetParameters(par);
      gauss2->SetParameter(1,0);
      gauss2->SetParameter(4,0);
      gauss2->SetParameter(8,0);
      gauss2->SetParameter(10,0);
      gauss2->SetLineColor(4);
      gauss2->Draw("CSAME");
      TF1* gauss3;
      gauss3 = new TF1("gauss3", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss3->SetParameters(par);
      gauss3->SetParameter(1,0);
      gauss3->SetParameter(4,0);
      gauss3->SetParameter(6,0);
      gauss3->SetParameter(10,0);
      gauss3->SetLineColor(4);
      gauss3->Draw("CSAME");
      TF1* gauss4;
      gauss4 = new TF1("gauss4", MyQuadroPeaks, fit_low, fit_high, npar);
      gauss4->SetParameters(par);
      gauss4->SetParameter(1,0);
      gauss4->SetParameter(4,0);
      gauss4->SetParameter(6,0);
      gauss4->SetParameter(8,0);
      gauss4->SetLineColor(4);
      gauss4->Draw("CSAME");
      TF1* tail1;
      tail1 = new TF1("tail1", MyQuadroTails, fit_low, fit_high, npar);
      tail1->SetParameters(par);
      tail1->SetParameter(6,0);
      tail1->SetParameter(8,0);
      tail1->SetParameter(10,0);
      tail1->SetLineColor(3);
      tail1->Draw("CSAME");
      TF1* tail2;
      tail2 = new TF1("tail2", MyQuadroTails, fit_low, fit_high, npar);
      tail2->SetParameters(par);
      tail2->SetParameter(4,0);
      tail2->SetParameter(8,0);
      tail2->SetParameter(10,0);
      tail2->SetLineColor(3);
      tail2->Draw("CSAME");
      function->Draw("CSAME");

      canvas->Print(canvas->GetTitle(),"ps");
      delete function;
      delete gauss1;
      delete gauss2;
      delete gauss3;
      delete gauss4;
    }
  
  Double_t red_chisqr = (fitfun->GetChisquare()/fitfun->GetNDF());

  delete fitfun;

  return red_chisqr;  
}