void Fit(TH1F* hmass);
TFile* ofile;
TGraph* gsigma;

//////////////////////////////////////////////////////
void GetTotalYield(){

  TFile* ifile = new TFile("histo_mass_all.root","read");

  ofile = new TFile("mass_all_fitted.root","recreate");

  gsigma = new TGraph();
  TString histo_name = "hmass_all";
  TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
  Fit(h1);
  ofile->Close();

  TCanvas* c2 = new TCanvas("c2","c2",1200,1200);
  c2->cd();
  gsigma->SetMarkerStyle(8);
  gsigma->Draw("ap");


}

//////////////////////////////////////////////////////
void Fit(TH1F* hmass){
  int rebin = 1;
  //hmass->Rebin(rebin);
  TCanvas* c1 = new TCanvas("c1","c1",1200,1200);
  c1->cd();
  hmass->Draw("E1");

  int Amin = 82;
  int Amax = 153;

  int A = Amin;

  int NumberOfA = Amax - Amin;
  Double_t para[3*NumberOfA];

  TF1* g[NumberOfA];
  double Integral[NumberOfA];
  double Integral_err[NumberOfA];
  double total_integral = 0;
  double Yield[NumberOfA];

  for(int i=0; i<NumberOfA; i++){
    g[i] = new TF1(Form("g%i",i),"gaus",A-0.4,A+0.4);   
    if(i==0){
      g[i]->SetParameter(0,200);
      g[i]->SetParameter(1,Amin);
      g[i]->SetParameter(2,0.3);
    }
    else if(i>0){
      g[i]->SetParameter(0,g[i-1]->GetParameter(0));
      g[i]->SetParameter(1,A);
      g[i]->SetParameter(2,g[i-1]->GetParameter(2));

      g[i]->SetParLimits(1,A-0.15,A+0.15);
      g[i]->SetParLimits(2,0.25,0.40);

    }

    hmass->Fit(g[i],"qBR+");

    g[i]->Draw("lsame");
    g[i]->GetParameters(&para[3*i]);

    A++;
  }


  TString total_func = "gaus(0)";
  for(int i=1; i<NumberOfA; i++){
    TString gaus_func = Form("gaus(%i)",3*i);
    total_func += "+" + gaus_func;
  }

  TF1* total = new TF1("total",total_func,Amin-0.5, Amin+NumberOfA);
  total->SetParameters(para);
  A = Amin;
  for(int i=0; i<NumberOfA; i++){
    total->SetParameter(3*i,g[i]->GetParameter(0));
    total->SetParameter(3*i+1,g[i]->GetParameter(1));
    total->SetParameter(3*i+2,g[i]->GetParameter(2));

    //cout << A << " " << g[i]->GetParameter(1) << endl;
    total->SetParLimits(3*i+1,A-0.2,A+0.2);
    total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.85,g[i]->GetParameter(2)*1.15);
    A++;
  }
  total->SetLineColor(4);
  total->SetNpx(1000);

  TFitResultPtr fit_result_total = hmass->Fit(total,"qBRS");

  TMatrixD cov = fit_result_total->GetCovarianceMatrix();
  cov.Print();


  total->Draw("same");
    
  A = Amin;
  for(int i=0; i<NumberOfA; i++){

    double Amplitude = total->GetParameter(3*i);
    double mean = total->GetParameter(3*i+1);
    double sigma = total->GetParameter(3*i+2);
    double Amplitude_err = total->GetParError(3*i);
    double sigma_err = total->GetParError(3*i+2);

    /*cout << "A=  " << A << endl;
    cout << "Amplitude= " << Amplitude << " / " << g[i]->GetParameter(0) << endl;
    cout << "Mean= " << mean << " / " << g[i]->GetParameter(1) << endl;
    cout << "Sigma= " << sigma << " / " << g[i]->GetParameter(2) << endl;*/


    Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
    //Integral_err[i] = sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
    Integral_err[i] = sqrt(Integral[i]);

    gsigma->SetPoint(i,A,sigma);

    total_integral += Integral[i];
    A++;
  }

  double histo_integral = hmass->Integral();
  hmass->Write();
  total->Write();
  gsigma->Write();

  string output_filename = "Total_Yield.dat";
  ofstream ofile;
  ofile.open(output_filename.c_str());

  // Yield calculation //
  A = Amin;
  for(int i=0; i<NumberOfA; i++){
    //Yield[i] = Integral[i]/total_integral*200;
    Yield[i] = Integral[i]/histo_integral*200;
    
    //double yield_err = Integral_err[i]/total_integral*200;
    double yield_err = Integral_err[i]/histo_integral*200;
    ofile << A << " " << Yield[i] << " " << yield_err << endl;
    A++;
  }
  ofile.close();

}