#include "TF1.h"
#include "TH1.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
#include "TStyle.h"

void proton_peaks() { 

   TH1D * h1 = new TH1D("h1","h1",1000,0,10000);


   TH1F *h= new TH1F("h", "trkmm gaussian fit", 1000, 0, 3000);
   ((TF1*)gROOT->GetFunction("gaus"))->SetParameters(1., 1835, 10);
  h->FillRandom("gaus");
  ((TF1*)gROOT->GetFunction("gaus"))->SetParameters(1, 1902, 120.);
  h->FillRandom("gaus");
   
  TF1 *g1 = new TF1 ("m1", "gaus", 1600, 2000);
  g1->SetLineColor(kRed);
  TF1 *g2 = new TF1 ("m2", "gaus", 1800 ,2100);
  g2->SetLineColor(kGreen);
  TF1 *f1 = new TF1("double_gaus", "0.1*gaus(0) + gaus(3)", 1600, 2100);
  f1->SetParNames("Constant 1", "Mean 1", "Sigma 1",
                  "Constant 2", "Mean 2", "Sigma 2");
  f1->SetLineColor(kBlue);
  
  gStyle->SetOptFit(1);
  
  h->Fit(g1, "R");
  h->Fit(g2, "R+");
  
  Double_t par[6];
  g1->GetParameters(&par[0]);
  g2->GetParameters(&par[3]);
  f1->SetParameters(par);
  
  h->Fit(f1, "R+");
  // h->Fit(f1, "+", "e1", 4, 16);
  h->Draw("e1");
  
  g1->Draw("SAME");
  g2->Draw("SAME");
  f1->Draw("SAME");
  
  std::cout << std::endl; g1->Print();
  std::cout << std::endl; g2->Print();
  std::cout << std::endl; f1->Print();



  //*****************TRY**********************************************
   /*  for (int i = 0; i < 5000; ++i) { */
   /*  h1->Fill(gRandom->Gaus(1836,50)); // Fill with 5000 gaussian number with mean 3 and width 1.5. */
   /*  h1->Fill(gRandom->Gaus(1902,60)); */
   /*  } */
   /* // fit first a single gaussian in range [0,5] */
   /* TFitResultPtr r1 = h1->Fit("gaus","S","",1750,1900);  // first fit */
   /* TFitResultPtr r2 = h1->Fit("gaus","S","",1800,2000);  // first fit */

   /* TF1 * f1 = new TF1("fitFunc","gaus(0)+gaus(3)"); */
   /* // get parameters and set in global TF1 */
  
   /* // parameters of first gaussian */
   /* f1->SetParameter(0,r1->Parameter(0)); */
   /* f1->SetParameter(1,r1->Parameter(1)); */
   /* f1->SetParameter(2,r1->Parameter(2)); */
   /* // parameters of second gaussian */
   /* f1->SetParameter(3,r2->Parameter(0)); */
   /* f1->SetParameter(4,r2->Parameter(1)); */
   /* f1->SetParameter(5,r2->Parameter(2)); */


   /* h1->Fit(f1); */
   /*******************************************************************
}