#include <iostream>
#include <fstream>
#include <vector>
#include <string>

#include "TTree.h"
#include "TBranch.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TBranch.h"
#include "TCut.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TString.h"
#include "TMultiGraph.h"

using namespace std;

//finito il terzo ciclo di campi magnetici, divido il file in due (i primi 3 sono con FC nostra e quelli dopo con FC di MilanoBi)
//prendo la serie delle tre migliore 
//sovrappongo i due plot e devono coincidere 

Int_t npeaks = 100;

Double_t fpeaks(Double_t *x, Double_t *par) {
   Double_t result = par[0] + par[1]*x[0];
   for (Int_t p=0;p<npeaks;p++) {
      Double_t norm  = par[3*p+2]; // "height" or "area"
      Double_t mean  = par[3*p+3];
      Double_t sigma = par[3*p+4];
#if defined(__PEAKS_C_FIT_AREAS__)
      norm /= sigma * (TMath::Sqrt(TMath::TwoPi())); // "area"
#endif /* defined(__PEAKS_C_FIT_AREAS__) */
      result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
}

void graph_xy_error_fit() {

Double_t par[3000];

  TFile *out_file_1 = new TFile("2023-04-27_scan_fit.root", "recreate");

  ifstream infile;
  infile.open("2023-04-27_scan_02.txt");


 Int_t Npts = 0;
      string line;
      while(getline(infile,line))
	{
	  Npts++;
	}

   //cout << Npts << endl; 
  
     infile.close();
     infile.open("2023-04-27_scan_02.txt");

     Npts = Npts/2; 

     cout << Npts << endl; 
  
  Double_t B[Npts];
      Double_t eB[Npts];
      Double_t i_DF[Npts];
      Double_t ei_DF[Npts];
      Double_t i_MiB[Npts];
      Double_t ei_MiB[Npts];
      Double_t i_TES[Npts];
      Double_t ei_TES[Npts];
      Double_t dummy;
      
      for(Int_t i=0; i< Npts; i++)
	{
	  infile >> B[i] >> eB[i] >> i_DF[i] >> ei_DF[i] >> dummy >> dummy >> i_MiB[i] >> ei_MiB[i] >> i_TES[i] >> ei_TES[i]; 
	  //cout << B[i] << " " << eB[i] << " " << i_DF[i] << endl;  //ok, funziona 
    //cout << "ciao" << endl;
	}


  TCanvas *c1 = new TCanvas();
  c1->Divide(1,2);
  //c1->cd(1);
  
  TGraphErrors *g = new TGraphErrors(Npts, B, i_TES, eB, ei_TES);
  //g->Draw("AP*");
  TH1F *h=new TH1F("h","test",490,4780,5050);
  h->Clear();
  h->Reset();
  auto nPoints = g->GetN(); // number of points in your TGraph
  for(int i=0; i < nPoints; ++i) {
  double x = g->GetX()[i];
  double y = g->GetY()[i];
  int k = h->FindBin(x);
  h->SetBinContent(k,y);
  //h->SetBinError(k, ei_TES);
  }

   h->Draw("HIST");
  TH1F *h2 = (TH1F*)h->Clone("h2");
   // Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h2,2,"",0.001);
   printf("Found %d candidate peaks to fit\n",nfound);

   TF1 *fline = new TF1("fline","pol1",4650,4750);

   h->Fit("fline","qn");
   // Loop on all found peaks. Eliminate peaks at the background level
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);

  npeaks = 0;
  Double_t *xpeaks;
   xpeaks = s->GetPositionX();
   for (Int_t p=0;p<nfound;p++) {
      Double_t xp = xpeaks[p];
      Int_t bin = h2->GetXaxis()->FindBin(xp);
      Double_t yp = h2->GetBinContent(bin);
      par[3*npeaks+2] = yp; // "height"
      par[3*npeaks+3] = xp; // "mean"
      par[3*npeaks+4] =2; // "sigma"
#if defined(__PEAKS_C_FIT_AREAS__)
      par[3*npeaks+2] *= par[3*npeaks+4] * (TMath::Sqrt(TMath::TwoPi())); // "area"
#endif /* defined(__PEAKS_C_FIT_AREAS__) */
      npeaks++;
   }
   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Now fitting: Be patient\n");
   TF1 *fit = new TF1("fit",fpeaks,4700,5060,3*npeaks+2);
   // We may have more than the default 25 parameters
   TVirtualFitter::Fitter(h2,npeaks);
   fit->SetParameters(par);
   //fit->SetNpx(1000);
   h2->Fit("fit");

}