//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+   Root Macro: GainFit2 for FACT                                           +
//+---------------------------------------------------------------------------+
//+   by: jens.buss@tu-dortmund.de                                            +
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//
//  - analyse gain of pixels from Dark-Count-Photon-Spektrum (pedestal on)
//  - read a root-file produced by "fpeak.C" and get Amplitudespektrum
//    of all pixel
//  - read a root-file produced by "fbsl.C" and get baseline
//  - calculate pixelgain from multi-gaussian-fit
//    over each pixels amplitudespektrum
//  - draw distribution of gains
//  - draw gain vs. pixel
//  - find Pixels with deviating gainvalue (e.g. dead , crazy Pixels)
//    and write their number to txt-file (e.g. date_run_gainfit.txt)
//  - save histograms into root-file (e.g. date_run_gainfit.root)
//  - Works also with ROI 300
//
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#include <TFile.h>
#include <TH2F.h>
#include <TH1D.h>
#include <TROOT.h>
#include <TF1.h>
#include <TTimer.h>
#include <TString.h>
#include <Getline.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TH1.h>
#include <TProfile.h>
#include <TStyle.h>

#include <iostream>
#include <fstream>
using namespace std;

//-----------------------------------------------------------------------------
// Decleration of variables
//-----------------------------------------------------------------------------
UInt_t NumberOfPixels;
#define MAX_SIGMA 4
  typedef struct{
  UInt_t PID;
  vector<float> fitParameters[9];
  Double_t amplSD;
  Double_t amplDT;
  Double_t cross;
  Double_t ampl_fspek_fit;
  Double_t ampl_fspek_fit2;
  Double_t bsl;
  bool crazy;
  } proj;
vector<proj> channel;

typedef struct{
  Double_t constant;
  Double_t mean;
  Double_t sigma;
} dist_t;
dist_t distribution;

//Verbosity for Histograms
TString verbos[4]={"0Q","Q", "", "V"};
//-----------------------------------------------------------------------------
// Decleration of Histograms
//-----------------------------------------------------------------------------

//Amplification Distribution
TH1F *hAmplDistribution = NULL;
TH1F *hAmplDistribution2 = NULL;
TH1F *hGainFitBsl = NULL;

//Lists of Histograms that have to be saved
TObjArray hList;
TObjArray hListPeakSpektra;
TObjArray hListCrazy;

//-----------------------------------------------------------------------------
// Functions
//-----------------------------------------------------------------------------
Double_t spek_fit(Double_t *, Double_t *);
Double_t spek_fit2(Double_t *, Double_t *);

void BookHistos();
void SaveHistograms(const char * );
bool SearchAnormal(UInt_t, Double_t, int verbosityLevel );

//-----------------------------------------------------------------------------
// Main function
//-----------------------------------------------------------------------------


int gainfit2(
        const char * InputRootFileName = "../analysis/fpeak/20111110_005/20111110_005_fpeak.root",
        const char * InputBaselineFileName = "../analysis/fpeak/20111110_005/20111110_005_fbsl.root",
        const char * OutTextFileName = "../analysis/fpeak/20111110_005/20111110_005_gainfit.txt",
        const char * RootOutFileName = "../analysis/fpeak/20111110_005/20111110_005_gainfit.root",
        bool showHistos = false,
        bool debug = false,
        int verbosityLevel = 4)
{

//-----------------------------------------------------------------------------
// Fit Verbosity
//-----------------------------------------------------------------------------

  UInt_t fitverb = verbosityLevel;
  if (verbosityLevel >= 3) fitverb =3;

//-----------------------------------------------------------------------------
// Histograms, Canvases and Fit Functions
//-----------------------------------------------------------------------------

//Open Rootfiles Files
    TFile * baseline = TFile::Open( InputBaselineFileName );
    TFile * file = TFile::Open( InputRootFileName );

//Baseline
    TH2F * hbaseline = (TH2F*)baseline->Get("histo_mean");
//Amplitude spectrum
    TH2F *hAmplSpek = (TH2F*)file->Get("hAmplSpek_cfd");
    NumberOfPixels = hAmplSpek->GetNbinsX();
    channel.resize(hAmplSpek->GetNbinsX());

//Amplitude Spectrum of a Pixel
    TH1D* hPixelAmplSpek[1440]; ///++++ TODO: Warning hardcoded pixelnumber, change that +++

//Book Histograms
    BookHistos();

//Define Fit Functions
    Double_t par[6];
    Double_t par2[6];
    TF1 *g1    = new TF1("g1","gaus",5,15);
    TF1 *g2    = new TF1("g2","gaus",15,25);
    //TF1 *g3    = new TF1("g3","gaus",-5,5);
    TF1 *g4    = new TF1("g4","gaus");
    //TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",7,60);
    TF1 *fspek_fit = new TF1("fspek_fit", spek_fit, 5, 60, 6);
    fspek_fit->SetParNames("1.Max", "Gain", "Sigma", "2.Max", "Baseline", "Xtalk");
    TF1 *fspek_fit2 = new TF1("fspek_fit2g3", spek_fit2, 5, 60, 4);
    fspek_fit2->SetParNames("1.Max", "Gain", "Sigma", "2.Max");
//Create Canvases
    TCanvas *cGainFit = new TCanvas("cGainFit", "Fit of Amplitude Spektra", 0,0, 1200,400);
    cGainFit->Clear();
    cGainFit->Divide(3, 2);

    TCanvas *cGraphs = new TCanvas("cGraphs", "gain vs. channel", 0,401, 1200,400);
    cGraphs->Clear();
    cGraphs->Divide(1, 2);

//Draw imported Histograms
    cGainFit->cd(1);
    gStyle->SetPalette(1,0);
    gROOT->SetStyle("Plain");
    hAmplSpek->SetTitle("Amplitude Spectrum of all Pixels");
    hAmplSpek->SetAxisRange(0, 100, "Y");
    hAmplSpek->Draw("COLZ");


    cGainFit->cd(3);
    hbaseline->SetTitle("Baselinedistribution of all Pixels");
    Double_t Xrange = hbaseline->GetMean(1);
    hbaseline->SetAxisRange(Xrange-1.5, Xrange + 1.5, "X");
    hbaseline->SetAxisRange(0, (hbaseline->GetBinContent(hbaseline->GetMaximumBin())+300), "Y");
    hbaseline->Draw();

    //-----------------------------------------------------------------------------
    // Graphs
    //-----------------------------------------------------------------------------

    Double_t simple_gain[NumberOfPixels];
    Double_t gain[NumberOfPixels];
    Double_t pixelvec[NumberOfPixels];
    for (UInt_t i=0; i<NumberOfPixels; i++){
       pixelvec[i]=i;
       simple_gain[i]=0;
       gain[i]=0;
    }

    TGraph *gr1 = new TGraph(NumberOfPixels,pixelvec,simple_gain);
    TGraph *gr2 = new TGraph(NumberOfPixels,pixelvec,gain);

    gr1->SetMarkerStyle(5);
    gr2->SetMarkerStyle(2);
    gr1->SetMarkerColor(kRed);
    gr2->SetMarkerColor(kBlue);
    gr1->SetMarkerSize(0.5);
    gr2->SetMarkerSize(0.5);
    gr1->SetTitle("G-APD Amplification vs. Channel (with Baseline Parameter)");
    gr2->SetTitle("G-APD Amplification vs. Channel (without Baseline Parameter)");
    TMultiGraph *mg = new TMultiGraph();
    mg->Add(gr1);
    mg->Add(gr2);

//-----------------------------------------------------------------------------
// Loop over all Pixels and fit Distribution of singles, doubles and tripples
//-----------------------------------------------------------------------------

    if (verbosityLevel > 1) cout << "Number of Pixels: " << NumberOfPixels << endl;
    // Begin of Loop
    for (UInt_t pixel = 0; pixel < NumberOfPixels; pixel++){

      channel[ pixel ].PID = pixel+0;

    //Title
         TString histotitle;
         histotitle += "hPixelPro_";
         histotitle += pixel+0 ;

      if (verbosityLevel > 1){
        cout << "-----------------------------------------------------------" << endl;
        cout << " PID: " << pixel+0 << endl;
        cout << "-----------------------------------------------------------" << endl;
      }
  //Amplitude Spectrum of a Pixel
        hPixelAmplSpek[ pixel ] = new TH1D;
        hPixelAmplSpek[ pixel ] = hAmplSpek->ProjectionY( histotitle , pixel+1, pixel+1);	 //Projectipon of each Pixel out of Ampl.Spectrum
        hPixelAmplSpek[ pixel ]->SetAxisRange(0, 60, "X");

        histotitle = "Amplitude Spectrum of Pixel #";
        histotitle += pixel+0 ;

        hPixelAmplSpek[ pixel ]->SetTitle(histotitle);
        hPixelAmplSpek[ pixel ]->SetYTitle("Counts [a.u.]");
        hListPeakSpektra.Add(hPixelAmplSpek[ pixel ]);


//-----------------------------------------------------------------------------
// Single Gaussian Fits and Tripple-Gaussian Fit Funciton
//-----------------------------------------------------------------------------

        //total->SetLineColor(2);
        fspek_fit->SetLineColor(kRed);
        fspek_fit2->SetLineColor(kBlue);
        fspek_fit2->SetLineStyle(2);

        cGainFit->cd(4);
        gPad->SetLogy(1);

        hPixelAmplSpek[ pixel ]->Fit(g1,verbos[fitverb]+"R");
        hPixelAmplSpek[ pixel ]->Fit(g2,verbos[fitverb]+"0R+");
        //hPixelAmplSpek[ pixel ]->Fit(g3,verbos[fitverb]+"0R+");
        g1->GetParameters(&par[0]);
        g1->GetParameters(&par2[0]);
        g2->GetParameters(&par[3]);
        g2->GetParameters(&par2[3]);

        //g3->GetParameters(&par[6]);

        //total->SetParameters(par);

        //hPixelAmplSpek[ pixel ]->Fit(total,verbos[fitverb]+"R+");
        //channel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1);
        //channel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4);
        //channel[ pixel ].amplMean = fspek_fit->GetParameter(1);

        par2[4]=par[4] = hbaseline->GetMean();
        cout << "starvalue of Baseline:" << par[4] << endl;
        fspek_fit->SetParameters(par);
        if (verbosityLevel > 2) cout << "spek_fit" << endl;
        hPixelAmplSpek[ pixel ]->Fit(fspek_fit,verbos[fitverb]+"R+");
        channel[ pixel ].ampl_fspek_fit = fspek_fit->GetParameter(1);
        channel[ pixel ].bsl = fspek_fit->GetParameter(4);
        if (verbosityLevel > 2) cout << "Baseline:" << channel[ pixel ].bsl << endl;

        fspek_fit2->SetParameters(par2);
        if (verbosityLevel > 2) cout << "fspek_fit2" << endl;
        hPixelAmplSpek[ pixel ]->Fit(fspek_fit2,verbos[fitverb]+"R+");
        channel[ pixel ].ampl_fspek_fit2 = fspek_fit2->GetParameter(1);

        if( showHistos ){
        hPixelAmplSpek[ pixel ]->Draw();
        }






//-----------------------------------------------------------------------------
// Draw Histograms
//-----------------------------------------------------------------------------

        cGainFit->cd(2);
        gPad->SetLogy(1);
        g4->SetLineColor(2);
        hAmplDistribution->Fill( channel[ pixel ].ampl_fspek_fit ) ;
        hAmplDistribution->Fit(g4, verbos[fitverb]);
        distribution.sigma = g4->GetParameter(2);
        distribution.mean = g4->GetParameter(1);
        distribution.constant = g4->GetParameter(0);

        cGainFit->cd(5);
        gPad->SetLogy(1);
        g4->SetLineColor(kBlue);
        hAmplDistribution2->Fill(channel[ pixel ].ampl_fspek_fit2) ;
        hAmplDistribution2->Fit(g4, verbos[fitverb]);

        cGainFit->cd(3);
        hGainFitBsl->Fill(channel[ pixel ].bsl);
        hGainFitBsl->Draw("same");
        if( showHistos ){
            if (pixel%40==0){
            cGainFit->cd(2);
            hAmplDistribution->Draw();
            cGainFit->cd(5);
            hAmplDistribution2->Draw();
            }

            cGainFit->cd(6);
            hGainFitBsl->Draw();
          }

//-----------------------------------------------------------------------------
// Draw Graphs of Distributions
//-----------------------------------------------------------------------------

        gr1->SetPoint(pixel, channel[pixel].PID, channel[pixel].ampl_fspek_fit);
        gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].ampl_fspek_fit2);

    if( showHistos ){
        if (pixel%60==0){
          cGraphs->cd(1);
          gr2->Draw("APL");
          cGraphs->cd(2);
          gr1->Draw("APL");
        }

        cGainFit->Modified();
        cGainFit->Update();
        cGraphs->Modified();
        cGraphs->Update();
    }

    //-----------------------------------------------------------------------------
    // Debug Modus
    //-----------------------------------------------------------------------------
        if(debug){

            TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
            timer.TurnOn();
            TString input = Getline("Type 'q' to exit, <return> to go on: ");
            timer.TurnOff();
            if (input=="q\n") {
              break;
            }
        }
    }
//-----------------------------------------------------------------------------
// End of Loop
//-----------------------------------------------------------------------------


//   ( cGainFit->cd(4);
//    hPixelAmplSpek[ NumberOfPixels ]->Draw();

    cGainFit->cd(2);
    hAmplDistribution->Draw();

    cGainFit->cd(5);
    hAmplDistribution2->Draw();

    cGainFit->cd(6);
    hGainFitBsl->Draw();
    hbaseline->Draw("same");

    cGainFit->Modified();
    cGainFit->Update();

    cGraphs->cd(1);
    gr2->Draw("APL");
    cGraphs->cd(2);
    gr1->Draw("APL");
    hList.Add(gr1);
    hList.Add(gr2);
    cGraphs->Modified();
    cGraphs->Update();

//-----------------------------------------------------------------------------
// Write into logfile
//-----------------------------------------------------------------------------
    string filename =  InputRootFileName;
    size_t lastslash = filename.find_last_of('/');
    filename = filename.substr(lastslash+1);
    size_t firstdot = filename.find_first_of('.');
    filename = filename.substr(0, firstdot);

    string::iterator it;
    for (it=filename.begin(); it < filename.end(); ){

            if (!isdigit(*it)){
                    it=filename.erase(it);
            }
            else
                    ++it;
    }
if (verbosityLevel > 0)    cout << filename << endl;

    ofstream out;
    out.open( OutTextFileName );
    out << "pixelID/I:gainSD/F:gainDT:filename/l" << endl;
    out << "# ---- the 1st line helps to use TTree::ReadFile, for reading again this txt file with root. " << endl;
    out << "#" << InputRootFileName << endl;
    out << "# PixelID : gain between Single and Double Peak : gain between Double and Tripple Peak " << endl;
    out << "# gain is derived from the parameters of a three gaussian fit function " << endl;
    for (UInt_t pixel=0; pixel<NumberOfPixels; pixel++){
            out << channel[ pixel ].PID << ";\t"    ;
            out << channel[ pixel ].ampl_fspek_fit << ";\t" ;
           // out << channel[ pixel ].amplDT << ";\t" ;
            out << channel[ pixel ].ampl_fspek_fit2 << "\t" ;

    //check if pixel gain is compareable to gain of the others and save result to logfile
            channel[ pixel ].crazy = SearchAnormal(pixel, MAX_SIGMA, verbosityLevel);
            if (channel[ pixel ].crazy) hListCrazy.Add(hPixelAmplSpek[ pixel ]);
            out << channel[ pixel ].crazy << "\t" ;

            out << filename << endl;

    }
    out.close();
SaveHistograms( RootOutFileName );
return( 0 );
}

//-----------------------------------------------------------------------------
// End of Main
//-----------------------------------------------------------------------------




//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+                                                                           +
//+                     Function definitions                                  +
//+                                                                           +
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



//-----------------------------------------------------------------------------
// Function: Book Histograms
//-----------------------------------------------------------------------------

void BookHistos(){

  //Amplification Distribution - Single2Double
  hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (with Baseline Parameter)", 500, 0 , 20 );
  hAmplDistribution->SetXTitle("Amplification [mV]");
  hAmplDistribution->SetYTitle("Counts");
  hAmplDistribution->SetAxisRange(7.5, 12.5, "X");
  hList.Add( hAmplDistribution );

  //Amplification Distribution2 - Double2Tripple
  hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(without Baseline Parameter)", 500, 0 , 20 );
  hAmplDistribution2->SetXTitle("Amplification [mV]");
  hAmplDistribution2->SetYTitle("Counts");
  hAmplDistribution2->SetAxisRange(7.5, 12.5, "X");
  hList.Add( hAmplDistribution2 );

  //cross vs Channels
  hGainFitBsl = new TH1F("hGainFitBsl","BaselineDistribution from plot", 400,-99.5,100.5);
  hGainFitBsl->SetYTitle("Counts");
  hGainFitBsl->SetXTitle("Voltage [mV]");
  hGainFitBsl->SetLineColor(2);
  hList.Add( hGainFitBsl );

}

//-----------------------------------------------------------------------------
// Save Histograms
//-----------------------------------------------------------------------------

void SaveHistograms( const char * loc_fname){
    // create the histogram file (replace if already existing)
    TFile tf( loc_fname, "RECREATE");

    hList.Write(); // write the major histograms into the top level directory
    tf.mkdir("PeakSpektra"); tf.cd("PeakSpektra"); // go to new subdirectory
    hListPeakSpektra.Write(); // write histos into subdirectory
    tf.mkdir("CrazyPixels"); tf.cd("CrazyPixels"); // go to new subdirectory
    hListCrazy.Write(); // write histos into subdirector
    tf.Close(); // close the file
}

//-----------------------------------------------------------------------------
// Fit-Function: Sum of several Gauß-Function with falling Maximum Amplitude
//-----------------------------------------------------------------------------

Double_t spek_fit(Double_t *x, Double_t *par){
  Double_t arg = 0;
  Double_t arg2 = 0;
  Double_t arg3 = 0;
  Double_t arg4 = 0;
  Double_t cross = 0;
  if (par[0] != 0) cross = par[3]/par[0];
  if (par[2] != 0) arg = (x[0] - par[1]-par[4])/par[2];
  if (par[2] != 0) arg2 = (x[0] - 2*par[1]-par[4])/par[2];
  if (par[2] != 0) arg3 = (x[0] - 3*par[1]-par[4])/par[2];
  if (par[2] != 0) arg4 = (x[0] - 4*par[1]-par[4])/par[2];
  Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg);
  Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
  Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
  Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
  par[5]=cross;
  return peak1 + peak2 + peak3 + peak4;
}

//-----------------------------------------------------------------------------
// Fit-Function: Sum of several Gauß-Function with overlying potential decay
//-----------------------------------------------------------------------------

Double_t spek_fit2(Double_t *x, Double_t *par2){

  Double_t arg = 0;
  Double_t arg2 = 0;
  Double_t arg3 = 0;
  Double_t arg4 = 0;
  Double_t cross = 0;
  if (par2[0] != 0) cross = par2[3]/par2[0];
  if (par2[2] != 0) arg = (x[0] - par2[1])/par2[2];
  if (par2[2] != 0) arg2 = (x[0] - 2*par2[1])/par2[2];
  if (par2[2] != 0) arg3 = (x[0] - 3*par2[1])/par2[2];
  if (par2[2] != 0) arg4 = (x[0] - 4*par2[1])/par2[2];
  Double_t peak1 = par2[0]*TMath::Exp(-0.5 * arg * arg);
  Double_t peak2 = cross*par2[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
  Double_t peak3 = cross*cross*par2[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
  Double_t peak4 = cross*cross*cross*par2[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
  return peak1 + peak2 + peak3 + peak4;
}

//-----------------------------------------------------------------------------
// Function: Search for Pixel that do not have a compareable gain to the Other
//-----------------------------------------------------------------------------

bool SearchAnormal(UInt_t pixel, Double_t distance, int verbosityLevel){
  Double_t check = channel[ pixel ].ampl_fspek_fit;
  // check = TMath::Sqrt(check*check);
  if ( check > ( distribution.mean - distance * distribution.sigma )
    && check < ( distribution.mean + distance * distribution.sigma )){
    return false;
    }
  else{
    if (verbosityLevel > 2) cout << "Not normal behaviour in Pixel: " << pixel << endl;
    return true;
  }
}