#define macro01_cxx

#include "macro01.h"
#include <TH2.h>
#include <TStyle.h>

void macro01::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).
   TString option = GetOption();

// ----------------------------------------------------------------------------
// input parameters   
   nChannels = 16384;                 // the number of channels the detector has

   // energy calibration
   cut = false;                       // for better peak detection the spectrum can be split in two parts; do you want to perform this cut?
   maxNumberOfPeaksPerPart = 20;      // each part can contain up to this number of peaks
   cutPosition = 100;                // the position of the cut given in ADC channels
   startCalibrationRegion = 2000;
   endCalibrationRegion = 14000;     
   endOfInterestingPart = 15000;     // end of the insteresting part of the spectrum (only relevant for display)
   peakFinderThreshold = 0.009;         // peaks with amplitude less than threshold*highest_peak are discarded; 0.05 by standard
    // peaks used for calibration https://www.radiochemistry.org/periodictable/gamma_spectra/pdf/th228.pdf
   ThalliumPeaks = {2614, 2614-511, 2614-2*511, 860, 583, 510};  // ?? am I supposed to use escape peaks for calibration ??
   BismuthPeaks = {1620, 727};
   //OtherPeaks = {238}; 

   // set output directory; the name will be chosen according to the input ROOT file
   std::string outputDirectory = "/mnt/e15/comellato/DriftTimeStudies/analysis/b19035/Th/NonCollimated/Side/Tier2_tom/";
   // end of input parameters
// ----------------------------------------------------------------------------
   // set up histograms
   hSpectrumUncalibrated = new TH1F("hSpectrumUncalibrated", "Th-228 -- uncalibrated", nChannels, 0, nChannels-1); // the first spectrum
   hSpectrumUncalibrated->SetXTitle("ADC channel");
   hSpectrumUncalibrated->SetYTitle("counts");
   // join all calibration peak energies in one array
   calibrationEnergies.insert(calibrationEnergies.end(), ThalliumPeaks.begin(), ThalliumPeaks.end());
   calibrationEnergies.insert(calibrationEnergies.end(), BismuthPeaks.begin(), BismuthPeaks.end());
   calibrationEnergies.insert(calibrationEnergies.end(), OtherPeaks.begin(), OtherPeaks.end());
   // and sort it
   sort( calibrationEnergies.begin(), calibrationEnergies.end() );
}

void macro01::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

}

Bool_t macro01::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // When processing keyed objects with PROOF, the object is already loaded
   // and is available via the fObject pointer.
   //
   // This function should contain the \"body\" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

   fReader.SetLocalEntry(entry); // the "Local" is relevant if you process a TChain and not only a TTree
   
   // here we assume, that the following array stores the relevant information at its 0th index 
   Double_t energy = GEMDEnergyGauss_energy[0];
   hSpectrumUncalibrated->Fill(energy); // entry is the number of the channel
   
   return kTRUE;
}

void macro01::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void macro01::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
  
  
  printf("-------------------------------------------------------\n");
  //--------------------- CALIBRATION
  // search the spectrum for peaks to perform a calibration on 
  std::vector<Double_t> calibrationPeaks = getCalibrationPeaks(hSpectrumUncalibrated, startCalibrationRegion, endCalibrationRegion, maxNumberOfPeaksPerPart, peakFinderThreshold);
  sort( calibrationPeaks.begin(), calibrationPeaks.end() );
  printf("expected number of peaks: %zu; found number of peaks: %zu\n", calibrationEnergies.size(), calibrationPeaks.size());
  if (calibrationEnergies.size() == calibrationPeaks.size()) { printf("--> Success!\n"); }
  else { printf("---> Oh no!\n"); }
  printf("The coarse peak positions (found by TSpectrum) are:\n");
  for(Int_t i=0; i<calibrationPeaks.size(); i++) { printf("Peak #%d, channel = %f\n", i, calibrationPeaks[i]); }
  printf("-------------------------------------------------------\n");
  hSpectrumUncalibrated->GetXaxis()->SetRange(0, endOfInterestingPart);
  // get their exact position by fitting them
  TH1F * hSpectrumCalibration = (TH1F*) hSpectrumUncalibrated->Clone("hSpectrumCalibration");
  hSpectrumCalibration->SetTitle("Th-228 -- calibration");
  std::vector<Double_t> tempCalibrationPeaks = fitPeaks(hSpectrumCalibration, calibrationPeaks, 50);
  calibrationPeaks = tempCalibrationPeaks;
  // create a calibration graph and a calibration (fit) function
  printf("And here come the fine peak positions (found by fitting) and their corresponding energies used for calibration:\n");
  TGraph * gCalibration = getCalibrationGraph(calibrationPeaks, calibrationEnergies);
  printf("-------------------------------------------------------\n");
  gCalibration->SetName("gCalibration");
  gCalibration->SetTitle("calibration graph");
  gCalibration->Fit("pol1");
  TF1 * fCalibration = gCalibration->GetFunction("pol1");
  // create a residual plot for the calibration
  std::vector<Double_t> resCalibrationEnergies( calibrationPeaks.size() );
  for(Int_t i=0; i<calibrationPeaks.size(); i++) {
    resCalibrationEnergies[i] = channel2Energy(fCalibration, calibrationPeaks[i]) - calibrationEnergies[i]; // don't know if that's correct
  }
  TGraph * gResCalibration = getGraph(calibrationPeaks, resCalibrationEnergies);
  gResCalibration->SetName("gResCalibration");
  gResCalibration->SetTitle("residuals of the calibration graph fit");
  gResCalibration->GetXaxis()->SetTitle("channel");
  gResCalibration->GetYaxis()->SetTitle("residuals: E_(calibrated) - E_(actual)");
  // use the calibration fit function to create a new, calibrated spectrum
  hSpectrumCalibrated = new TH1F("hSpectrumCalibrated", "Th-228 -- calibrated", nChannels, 0, channel2Energy(fCalibration, nChannels-1));
  for(Int_t i=0; i<nChannels; i++) {
    hSpectrumCalibrated->SetBinContent(i, hSpectrumUncalibrated->GetBinContent(i));
  }
  hSpectrumCalibrated->SetXTitle("energy (keV)");
  hSpectrumCalibrated->SetYTitle("counts");
  //--------------------

  //-------------------- FITTING
  
    
  //--------------------
  
  
  // store the results in a ROOT file
  /*std::string fileName = _file0->GetName(); doenst work for TChains
  fileName = fileName.substr(39, 27);   // truncate the path so that only the actual filename remains */
  std::string fileName = "analysis.root";
  std::string outputDirectory = "/mnt/e15/comellato/DriftTimeStudies/analysis/";
  std::string fullOutputPath = std::string(outputDirectory) + std::string(fileName);
  TFile * analysis = new TFile( fullOutputPath.c_str(), "RECREATE" );

  hSpectrumUncalibrated->Write();
  hSpectrumCalibration->Write();
  hSpectrumCalibrated->Write();
//  gCalibration->Write();
//  gResCalibration->Write();

  analysis->Close();

  printf("the analysis should have been saved to %s\n", fullOutputPath.c_str());
  
  //printf("this is the file name: %s\n", fileName.c_str());  
  /*
*/
  
  // the idea for the following code is something like
  /* 
   * create a resolution plot (i.e. find the FWHM of various peaks and plot them wtr peakEnergy)
   * a) do it for the fitted functions already created at calibration --> find the list of associated functions
   * b) fit functions at the calibrated historam as well
   *
   * fix this stupid bug which crashes root when Draw() is acted upon hSpectrumCalibrated
   *
   */ 
}