#include <cstring>
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <cstdlib>
#include <chrono>
#include <math.h>       /* modf */
#include <cassert>
#include <cmath>
#include <stdalign.h>
#include </usr/local/include/fftw3.h>
#include <fstream>
#include <algorithm>    // std::find
#include <vector>       // std::vector
#include <iomanip>
#include <limits>
#include <execution>
#include <numeric>
#include <utility>
#include <map>
#include <stdexcept>
#include <stdint.h>
#include <TGraph.h>
#include <TF1.h>
#include <TString.h>
#include </usr/local/include/snappy.h>
#include "TMath.h"

#define RESET   "\033[0m"
#define BLACK   "\033[30m"      /* Black */
#define RED     "\033[31m"      /* Red */
#define GREEN   "\033[32m"      /* Green */
#define YELLOW  "\033[33m"      /* Yellow */
#define BLUE    "\033[34m"      /* Blue */
#define MAGENTA "\033[35m"      /* Magenta */
#define CYAN    "\033[36m"      /* Cyan */
#define WHITE   "\033[37m"      /* White */
#define BOLDBLACK   "\033[1m\033[30m"      /* Bold Black */
#define BOLDRED     "\033[1m\033[31m"      /* Bold Red */
#define BOLDGREEN   "\033[1m\033[32m"      /* Bold Green */
#define BOLDYELLOW  "\033[1m\033[33m"      /* Bold Yellow */
#define BOLDBLUE    "\033[1m\033[34m"      /* Bold Blue */
#define BOLDMAGENTA "\033[1m\033[35m"      /* Bold Magenta */
#define BOLDCYAN    "\033[1m\033[36m"      /* Bold Cyan */
#define BOLDWHITE   "\033[1m\033[37m"      /* Bold White */

std::ostream& bold_on(std::ostream& os)
{
    return os << "\e[1m";
}

std::ostream& bold_off(std::ostream& os)
{
    return os << "\e[0m";
}

#include "TArray.h"
#include "TROOT.h"
#include "TCutG.h"
#include "TBuffer.h"
#include "TVirtualPad.h"
#include "TPaveText.h"
#include "TH2.h"
#include "TH3.h"
#include "strlcpy.h"
#include "gaussian_fit.C"
//#include "FitStiefel.C"

int Gaussian_fit_calibrated_RF()
{

  TFile *fData1 = TFile::Open("Energy_calibration_file.root");
  TTree *Calibrated = (TTree*)fData1->Get("Calibrated");

//number of records per file
//detector 97 was connected to channel B in the digitizer
//detector 99 was connected to channel C in the digitizer

  const int run = 7;
  const char channel = 'B';

  //***********                      0                             1                            2                           3                           4                      5                           6                             7                            8                         9                     10                        11                         12                         13                         14                        15                         16                       17                         18                          19                         20                         21                       22                 23 
  const char *setting[24] = {"A 100  Z 40 Q 20 E 80"  ,  "A 100  Z 40 Q 20 E 120"  ,  "A 100  Z 40 Q 20 E 110"  ,  "A 100 Z 40  Q 20 E 100" , "A 100 Z 40  Q 20 E 90", "A 100  Z 40 Q 20 E 80"  ,  "A 100  Z 40 Q 20 E 70"  ,  "A 100  Z 40 Q 20 E 60"  ,  "A 130  Z 55 Q 20 E 30" , "A 130  Z 55 Q 20 E 40" , "A 130  Z 55 Q 20 E 50" , "A 130  Z 55 Q 20 E 60" , "A 130  Z 55 Q 20 E 70" , "A 130  Z 55 Q 20 E 80" , "A 136  Z 55 Q 21 E 80" , "A 136  Z 55 Q 21  E 70" , "A 136  Z 55 Q 21  E 60" , "A 136  Z 55 Q 21  E 50" , "A 136  Z 55 Q 21  E 40" , "A 140  Z 58 Q 21  E 50" , "A 140  Z 58 Q 21  E 60" , "A 140  Z 58 Q 21  E 70" , "A 4  Z 2  E 4.75" , "A 4  Z 2  E 12"};


  /*Tuning parameters for gaussian fit*/
  
  const double sig = 5.; //sigma run 5, 6, 7 --->  sig = 5, limit = 0.008 run 2, 3, 4 ---->  sig = 10, limit = 0.01 run 8 sig = 10, limit = 0.1
  const double limit = 0.008; //height*limit

  TCanvas *PlotPH_FF = new TCanvas("PlotPH_FF", "PlotPH_FF", 1550, 700, 600, 550);
  TH1F *PH_FF=new TH1F("PH_FF", setting[run], 300, -1, -1);
  PH_FF->SetLineColor(kOrange);
  Calibrated->Draw("Energy_calibration>>PH_FF", "", "");

  PH_FF->Print();

  TH1F *PH_FF_clone = (TH1F*)PH_FF->Clone("PH_FF_clone");
  TSpectrum *s1 = new TSpectrum(2*npeaks);
  Int_t npeaks = s1->Search(PH_FF_clone,sig,"",limit);
  printf("npeaks %i\n", npeaks);

  Int_t numbpeaks[1]={0};

  Double_t meanFF[npeaks];  
  Double_t sigmaFF[npeaks];
  Double_t heightFF[npeaks];

  gaussian_fit(npeaks, sig, limit, PH_FF, meanFF, sigmaFF, heightFF, numbpeaks);

  return 0;
}