// working at 27-06-17
// Script to plot PHD Vs ADC channels of CAEN Digitizer N6724a
// Optimized to fit Energy Loss Distribution of beta (from 90Sr decay) impinging over a diamond detector.
// © CC BY-NC Pietro Oliva & Rino Conte @ UniRoma3 university, Dept. of Science  
// Thanx to bellenot from https://root-forum.cern.ch for the reading section


#include <iostream>   // std ::cout
#include <fstream>    // std::ifstream
#include <vector>
#include <string>
#include "TFile.h"
#include <TCanvas.h>
#include "TLegend.h"
#include "TGraph.h"
#include "TSystemDirectory.h"
#include "TList.h"
#include "TF1.h"
#include "TAxis.h"
#include "TGAxis.h"
#include "TGButton.h"
#include "TTimer.h"

using namespace std;

TF1 *fSpectrum = 0;

void updateFitting()
{
   static TF1 *l1 = 0, *l2 = 0, *g = 0;
   Double_t param, Gamp, Gmean, Gsigma; 
   Double_t L1amp, L1mean, L1sigma;
   Double_t L2amp, L2mean, L2sigma;
   Double_t XRangeMax = 18000.0, XRangeMin = 0.0;
   if (fSpectrum == 0) return;
   L1amp=fSpectrum->GetParameter(0);
   L1mean=fSpectrum->GetParameter(1);
   L1sigma=fSpectrum->GetParameter(2);
   L2amp=fSpectrum->GetParameter(3);
   L2mean=fSpectrum->GetParameter(4);
   L2sigma=fSpectrum->GetParameter(5);
   Gamp=fSpectrum->GetParameter(6);
   Gmean=fSpectrum->GetParameter(7);
   Gsigma=fSpectrum->GetParameter(8);

   if (l1 == 0) {
      //Draw the pdfs
      l1 = new TF1("l1","[0]*(TMath::Landau(x,[1],[2],0))", XRangeMin, XRangeMax);
      l1->SetLineStyle(2);
      l1->SetLineColor(6);
   }
   l1->SetParameters(L1amp,L1mean,L1sigma);
   l1->Draw("same");

   if (l2 == 0) {
      l2 = new TF1("l2","[3]*(TMath::Landau(x,[4],[5],0))", XRangeMin, XRangeMax);
      l2->SetLineStyle(5);
      l2->SetLineColor(8);
   }
   l2->SetParameters(L2amp,L2mean,L2sigma);
   l2->Draw("same");

   if (g == 0) {
      g = new TF1("g","gaus",XRangeMin,XRangeMax);
      g->SetLineStyle(9);
      g->SetLineColor(4);
   }
   g->SetParameters(Gamp,Gmean,Gsigma);
   g->Draw("same");

   gPad->Update();
}

void startUpdateTimer()
{
   TGTextButton *btn = (TGTextButton *) gTQSender;
   if (btn) {
      TString text = btn->GetString();
      if (text == "Fit") {
         // one can adjust the time (in ms), depending on the fitting time...
         TTimer::SingleShot(1000, 0, 0, "updateFitting()");
      }
   }
}

   ////====================================
   //
   //
   //  READING SECTION useful to read multiple
   //  two column files data like .dat or .txt
   //  and plot & save in .root
   //
   //
   ////====================================

// Search for the .dat files in the declared folder and generate a vector for the loop
void list_files(vector<string> &filenames, const char *dirname=" ",
                const char *ext=".dat")
{
   TSystemDirectory dir(dirname, dirname);
   TList *files = dir.GetListOfFiles();
   if (files) {
      TSystemFile *file;
      TString fname;
      TIter next(files);
      while ((file=(TSystemFile*)next())) {
         fname = file->GetName();
         if (!file->IsDirectory() && fname.EndsWith(ext)) {
         string filename(fname.ReplaceAll(ext, "").Data());
         cout << filename << endl;
            filenames.push_back(filename);
         }
      }
   }
}

//User inputs
void N6724_eh_s()
{
   vector<string> filenames;
   string filename, infile, fileout, extin, extout, input_dir;

   extin = ".dat";
   extout = ".root";
   cout << "Please enter the directory containing the .dat files you wish to plot\n"
           "NOTE: if you are already in the right directory, just type \"./\"\n"
           "Now, Enter Directory: " << endl;
   cin >> input_dir;
   cout << "Reading files from " << input_dir << endl;    //Just a check
   list_files(filenames, input_dir.c_str(), extin.c_str());

   for (auto name : filenames) {
      infile = input_dir + name + extin;
      fileout = name + extout;
      // Graphic Section

      new TCanvas(infile.c_str(), "Digitizer N6724a");     
      TGraph *gr = new TGraph(infile.c_str(), "%lg %lg");  // reads two columns
      gr->SetLineStyle(1);                            // Aesthetics
      gr->SetLineColor(kOrange-2);                    // Orange color 
      gr->SetLineWidth(1);                            // line width
      gr->GetXaxis()->SetTitle("ADC Channels");       // X Axis title
      gr->GetYaxis()->SetTitle("Counts");             // Y Axis title
      gr->SetTitle(infile.c_str());                   // Naming the graph like the .dat file
      gr->Draw("AL");                                 // Line with no markers
      gr->GetXaxis()->CenterTitle();                
      gr->GetYaxis()->CenterTitle();
      gr->GetXaxis()->SetNoExponent(kTRUE);
      gr->GetYaxis()->SetNoExponent(kTRUE);
      TGaxis *myY = (TGaxis*)gr->GetYaxis();
      myY->SetMaxDigits(2);

      TQObject::Connect("TGTextButton", "Clicked()", 0, 0, "startUpdateTimer()");

      ////====================================
      //
      //
      //  FIT SECTION 90Sr decay beta- in matter
      //  if fitting is not needed comment out
      //  till END line
      //
      //
      ////====================================

      // Input variables declaration
      Double_t L1_Counts, L1_Location, L1_Scale;
      Double_t L2_Counts, L2_Location, L2_Scale;
      Double_t G_Counts, G_Mean, G_Sigma;

#if 0
      // FIRST GUESS: request user input 3 guess for Landau1,
      // 3 for Landau2, 3 for gauss noise.
      //=======================================
      // LANDAU Sr->Y
      //
      cout << "Max Counts expected (Landau from Sr->Y decay)? ";
      scanf("%lf", &L1_Counts);
      cout << "Location parameter (Landau from Sr->Y decay)? ";
      scanf("%lf", &L1_Location);
      cout << "Scale parameter (Landau from Sr->Y decay)? ";
      scanf("%lf", &L1_Scale);
  
      cout << "First Guess (user input): p0= " << L1_Counts
           << " p1 = " << L1_Location
           << " p2 = " << L1_Scale<<endl;
  
      // LANDAU Y->Zr
      //
      cout << "Max Counts expected (Landau from Y->Zr decay)? ";
      scanf("%lf", &L2_Counts);
      cout << "Location parameter (Landau from Y->Zr decay)? ";
      scanf("%lf", &L2_Location);
      cout << "Scale parameter (Landau from Y->Zr decay)? ";
      scanf("%lf", &L2_Scale);
  
      cout << "First Guess (user input): p3= " << L2_Counts
           << " p4 = " << L2_Location
           << " p5 = " << L2_Scale<<endl;

      // GAUSS
      //
      cout << "Max Counts expected (Gaussian noise)? ";
      scanf("%lf", &G_Counts);
      cout << "Mean (Gaussian noise)? ";
      scanf("%lf", &G_Mean);
      cout << "Sigma  (Gaussian noise)? ";
      scanf("%lf", &G_Sigma);
  
      cout << "First Guess (user input): p6= " << G_Counts
           << " p7 = " << G_Mean
           << " p8 = " << G_Sigma<<endl;
#endif

      L1_Counts = 2000.0, L1_Location = 400.0, L1_Scale = 1.0;
      L2_Counts = 1400.0, L2_Location = 400.0, L2_Scale = 1.0;
      G_Counts = 1200.0, G_Mean = 400.0, G_Sigma = 100.0;
   
      // Fit Section
      TF1 *fSignalY = new TF1("fSignalY","landau(0)", 0., 1800.); // Landau1 ruling 90Sr->90Y
      TF1 *fSignalZr = new TF1("fSignalZr","landau(3)", 0., 1800.); // Landau2 ruling 90Sr->90Zr
      TF1 *fBackground = new TF1("fBackground","gaus(6)", 0., 1800.); // Gauss noise
      fSpectrum = new TF1("fSpectrum","gaus(6)+landau(3)+landau(0)", 0., 1800.); //Sum over pdfs
      fSpectrum->SetParameters(L1_Counts, L1_Location,  L1_Scale, L2_Counts, L2_Location,  L2_Scale,  G_Counts, G_Mean, G_Sigma);
      gr->Fit("fSpectrum","R");

      // Drawing the pdfs, data, fit altogether

      updateFitting();

      ////JUST CHECKS uncomment if VERBOSE is wanted 
      //param = fSpectrum->GetNpar();
      //cout<<"N. di parametri: "<<param<<endl;
      //cout<<" landau1: p0=amplitude p1=location parameter (~the most probable value) p2=scale parameter (~sigma of truncated approx) "<<endl;
      //cout<<" landau2: p3=amplitude p4=location parameter (~the most probable value p5=scale parameter (~sigma of truncated approx) "<<endl;
      //cout<<" la gaussiana è definita come g(x) = p6*exp(-0.5*((x-p7)/p8)^2)) "<<endl;
 
      //cout<<"Landau1 Amplitude: "<<L1amp<<" Landau1 mean: "<<L1mean<<" Landau1 sigma: "<<L1sigma<<endl;
      //cout<<"Landau2 Amplitude: "<<L2amp<<" Landau2 mean: "<<L2mean<<" Landau2 sigma: "<<L2sigma<<endl;
      //cout<<"Gauss Amplitude: "<<Gamp<<" Gauss mean: "<<Gmean<<" Gauss sigma: "<<Gsigma<<endl;

      //// Uncomment if legenda is needed
      //TLegend *leg=new TLegend(0.7,0.8,0.9,0.9,"");    
      //leg->SetHeader("Contributions to fit line:","C");
      //leg->SetBorderSize(0); // get rid of the box
      //leg->AddEntry(l1,"2,28 MeV #beta^{-}","l");
      //leg->AddEntry(l2,"0,546 MeV #beta^{-}","l");
      //leg->AddEntry(g,"Gaussian background","l");
      //leg->Draw("Same");

      //===========================
      //
      // END OF FIT SECTION
      //
      //==========================

      // Open, Write, Save .root
      //TFile *f0=new TFile(fileout.c_str(),"RECREATE");  
      //gr->Write();   
      //f0->Close();   
      //delete f0;
   }
}