double g2(Double_t *x, Double_t *par) {
   double r1 = (x[0]-par[1])/par[2];
   return par[0]*TMath::Exp(-0.5*(r1*r1));
}
double totalfunc(Double_t *x, Double_t *par) {
   double par_gaus1[3];
   par_gaus1[0] = par[0]; // height
   par_gaus1[1] = par[1]; // mean
   par_gaus1[2] = par[2]; // sigma

   double par_gaus2[3];
   par_gaus2[0] = par[3]; // height
   par_gaus2[1] = par[4]; // mean
   par_gaus2[2] = par[5]; // sigma

   double par_gaus3[3];
   par_gaus3[0] = par[6]; // height
   par_gaus3[1] = par[7]; // mean
   par_gaus3[2] = par[8]; // sigma

   double par_gaus4[3];
   par_gaus4[0] = par[3]; // height
   par_gaus4[1] = par[9]; // mean
   par_gaus4[2] = par[5]; // sigma

   double par_gaus5[3];
   par_gaus5[0] = par[0]; // height
   par_gaus5[1] = par[10]; // mean
   par_gaus5[2] = par[2]; // sigma

   return g2(x,par_gaus1)+g2(x,par_gaus2)+g2(x,par_gaus3)+g2(x,par_gaus4)+g2(x,par_gaus5);
}

void fit_5Gaussians()
{

   TCanvas *c1 = new TCanvas();
   Double_t x[29] = { 33.83562, 35.75342, 37.94521, 39.58904, 41.78082, 43.42466, 45.34247, 47.26027, 49.17808,
   51.09589, 53.28767, 55.20548, 56.84932, 59.0411, 60.9589, 62.60274, 64.79452,   66.71233, 68.63014, 70.54795, 72.46575,
   74.38356, 76.30137, 78.21918, 80.13699, 82.05479, 83.9726, 85.89041, 87.80822 };
   Double_t y[29] = { 0.4141745, 2.099813, 3.643964, 11.80912, 8.987073, 15.74378, 8.837867, 10.10097, 4.885202,
    5.584925, 4.312174, 6.561194, 5.571419, 7.39726, 6.265998, 7.3889, 5.552769, 6.534182, 4.276159, 5.539263, 4.971381,
    10.0373, 8.906039, 15.6621, 8.897035, 11.70943, 3.676764, 2.122966, 0.4283234 };

   TGraph *graph = new TGraph(29,x,y);
   graph->SetMarkerStyle(8);
   graph->Draw("APL");

   // Get starting values by fitting a Gaussin inside a range
   double g2par[3];
   g2par[0] = 10;
   g2par[2] = 10;

   TF1 *fit1 = new TF1("fit1",g2, 34, 38,3);
   fit1->SetLineColor(kMagenta);
   fit1->SetLineStyle(2);
   g2par[1] = 36.0;
   fit1->SetParLimits(1,34,38);
   fit1->SetParameters(g2par);
   graph->Fit(fit1, "R0+");

   TF1 *fit2 = new TF1("fit2",g2, 38, 50,3);
   fit2->SetLineColor(kRed);
   fit2->SetLineStyle(2);
   g2par[1] = 44.0;
   fit2->SetParameters(g2par);
   graph->Fit(fit2, "R0+");

   TF1 *fit3 = new TF1("fit3",g2, 50, 72,3);
   fit3->SetLineColor(kGreen);
   fit3->SetLineStyle(2);
   g2par[1] = 61.0;
   fit3->SetParameters(g2par);
   graph->Fit(fit3, "R0+");

   TF1 *fit4 = new TF1("fit4",g2, 72, 84,3);
   fit4->SetLineColor(kRed);
   fit4->SetLineStyle(2);
   g2par[1] = 80.0;
   fit4->SetParameters(g2par);
   graph->Fit(fit4, "R0+");

   TF1 *fit5 = new TF1("fit5",g2, 84, 88,3);
   fit5->SetLineColor(kMagenta);
   fit5->SetLineStyle(2);
   g2par[1] = 86.0;
   fit5->SetParLimits(1,84,88);
   fit5->SetParameters(g2par);
   graph->Fit(fit5, "R0+");

   double par[15];
   fit1->GetParameters(&par[0]);
   fit2->GetParameters(&par[3]);
   fit3->GetParameters(&par[6]);
   fit4->GetParameters(&par[9]);
   fit5->GetParameters(&par[12]);

   double partot[11];

   partot[0] = (par[0]+par[12])/2; // height gaus1,5   0-2,3-5,6-8,9-11,12-14
   partot[1] = par[1];             // mean   gaus1
   partot[2] = (par[2]+par[14])/2; // sigma  gaus1,5

   partot[3] = (par[3]+par[9])/2;  // height gaus2,4
   partot[4] = par[4];             // mean   gaus2
   partot[5] = (par[5]+par[11])/2; // sigma  gaus2,4

   partot[6] = par[6];             // height gaus3
   partot[7] = par[7];             // mean   gaus3
   partot[8] = par[8];             // sigma  gaus3

   partot[9] = par[10];            // mean   gaus4

   partot[10] = par[13];           // mean   gaus5

   TF1 *totalFit = new TF1("total",totalfunc,0,90,11);
   totalFit->SetLineColor(kBlue);
   totalFit->SetParameters(partot);
    
   // height limits gaus 1,2,4,5
   totalFit->SetParLimits(0,0.8*partot[0],1.2*partot[0]);
   totalFit->SetParLimits(3,0.8*partot[3],1.2*partot[3]);

   // mean limits gaus 1,2,4,5
   totalFit->SetParLimits(1,partot[1]-2,partot[1]+2);
   totalFit->SetParLimits(4,partot[4]-5,partot[4]+5);
   totalFit->SetParLimits(9,partot[9]-5,partot[9]+5);
   totalFit->SetParLimits(10,partot[10]-2,partot[10]+2);

   // sigma limits gaus 1,5
   totalFit->SetParLimits(2,0.8*partot[2],1.2*partot[2]);

   graph->Fit(totalFit, "IR+");

   //Draw the individual peaks
   fit1->SetParameter(0, totalFit->GetParameter(0)); //0
   fit1->SetParameter(1, totalFit->GetParameter(1)); //1
   fit1->SetParameter(2, totalFit->GetParameter(2)); //2

   fit2->SetParameter(0, totalFit->GetParameter(3)); //3
   fit2->SetParameter(1, totalFit->GetParameter(4)); //4
   fit2->SetParameter(2, totalFit->GetParameter(5)); //5

   fit3->SetParameter(0, totalFit->GetParameter(6)); //6
   fit3->SetParameter(1, totalFit->GetParameter(7)); //7
   fit3->SetParameter(2, totalFit->GetParameter(8)); //8

   fit4->SetParameter(0, totalFit->GetParameter(3)); //9
   fit4->SetParameter(1, totalFit->GetParameter(9)); //10
   fit4->SetParameter(2, totalFit->GetParameter(5)); //11

   fit5->SetParameter(0, totalFit->GetParameter(0)); //12
   fit5->SetParameter(1, totalFit->GetParameter(10)); //13
   fit5->SetParameter(2, totalFit->GetParameter(2)); //14

   fit1->SetRange(0,90);
   fit1->Draw("same");
   fit2->SetRange(0,90);
   fit2->Draw("same");
   fit3->SetRange(0,90);
   fit3->Draw("same");
   fit4->SetRange(0,90);
   fit4->Draw("same");
   fit5->SetRange(0,90);
   fit5->Draw("same");
}