#include <fstream>
#include <iostream>
#include "TROOT.h"
#include "TString.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TPad.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TLine.h"
#include "TROOT.h"
#include "TString.h"
#include <TStyle.h>

#include "PlotScurvesAll-diffFitConst.h"
#include "PlotScurvesAll-diffFit.h"

//#if defined(__CINT__) && !defined(__MAKECINT__)
//#include "PlotScurveAll-diffFit.h"
//#else
//#include "PlotScurvesAll-diffFit.h"
//#endif



// for format of EC_ASIC Scurves_all_ASIC data
// ASIC A,B,C,.... : i=0,1,2,... ch 0,1,...,63 : j=0,1,...,63
// % cat scurves_PC_inj_run5.lvm |awk '{i=2;j=37;print($1,$(i*64+j+1+1))}'
//   > scurves_PC_inj_run5-C-ch37.data
// => % cd $datadir
// => % sh $srcdir/makeDataFile.sh with assigining proper fname
// PlotScurvesAll-diffFit-v6.6.C: dac_sfit[ch]=-1; in the case of fitting failure.

using namespace std;
Int_t          kXmin       = 0;

//--------------------------------------------------------------------------------------------------------
void Help()
{
    
    cout << "This is the Help of the PlotScurveAll-diffFit" <<endl;
    cout << " Version 1.0 frol September 24 2013 " <<endl;
    cout << " Main author Hiroko Miyamoto " << endl;
    cout << " ==> Please in ROOT interpreter, issue the command .L PlotScurvesAll-diffFit.cc " <<endl;
    cout << " ==> Then play with PlotScurve() function " << endl;
    cout << " ==> If you need more information, please type Menu()" << endl;
    
}
//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
void Menu()
{
    cout << "This is the Detailed menu information for running the PlotScurveAll-diffFit" <<endl;
    
    cout << " void PlotScurve( " << endl;
    cout << " \t " <<  " Int_t    kDrawCh=31, " << endl;
    cout << " \t " <<  " TString  fname=kFname, " << endl;
    cout << " \t " <<  " TString  fname_ped=kFnamePed, //ASIC_D ped_file " << endl;
    cout << " \t " <<  " Double_t dac_half_spe_fix = -1," << endl;
    cout << " \t " <<  " TString  fname_gain=\"gain-org.txt\"," << endl;
    cout << " \t " <<  " Int_t    fit_min=kFitMin," << endl;
    cout << " \t " <<  " Int_t    fit_max=kFitMax," << endl;
    cout << " \t " <<  " Int_t    fit_min_scurve=kFitMinS," << endl;
    cout << " \t " <<  " Int_t    fit_max_scurve=kFitMaxS," << endl;
    cout << " \t " <<  " Bool_t   kBoard=1, //0:ASIC_C, 1:ASIC_D " << endl;
    cout << " \t " <<  " Int_t    kChRef=kRefCh,//ASIC_C " << endl;
    cout << " \t " <<  " Bool_t   checkAve=0 //0:use ref pix, 1:check average " << endl;
    cout << " ) " << endl;
    cout << "Set the parameters on the top of the program depending on your analysis:  " << endl;
    cout << "General policy: parameters starting with a \"k\" at the beginning is fixed gloval parameter and usually set in the header file, never change through the program. " << endl;
    cout << " : " << endl;
    cout << " - kDrawCh: channel just you want to check on display. doens't effect the analysis. " << endl;
    
    cout << " - fname, fname_ped (kFname, kFnamePed): \"YOUR DATA FILE\" without .data"  << endl;
    cout << " See examples in the header file. " << endl;
    cout << " - dac_half_spe_fix: set DAC value only when you want to adjust the gain to particular pulse height (DAC value) instead of either average of the pulse height of reference pixel. Usually just leave it -1. " << endl;
    cout << " - fit_min, fit_max (kFitMin, kFitMax): leave them as default, it is not currently used for the analysis."  << endl;
    cout << " - fit_min_scurve, fit_max_scurve (kFitMinS, kFitMaxS): define the range of s-curve fitting. search proper values by yourself with running root program a few times."  << endl;
    cout << " - kBoard: don't forget to change 0(C) or 1(D) depending on which ASIC you're analysing "  << endl;
    cout << " - kChRef (kRefCh): a pixel whose gain remains 64 thorough the gain adjustment. " << endl;
    cout << " - checkAve: for the first run, set this parameter 1, and find a pixel which remains 64. In the case of parameter 1, all the gains are adjusted to the average pulse height. " << endl;
    cout <<  " Once you choose the reference pixel, set the parameter 0, as well as set kRefCh the pixel number. (better to keep in the header file for the record of each PMTs.) " << endl;
    
    cout << " Actual configuration " << endl;
    cout << " \t " <<  "TString kFname = " << kFname << endl;
    cout << " \t " <<  "TString  kFnamePed = " << kFnamePed  << endl;
    cout << " \t " <<  "Int_t kFitMin = " << kFitMin << endl;
    cout << " \t " <<  "Int_t kFitMax = " << kFitMax << endl;
    
}
//--------------------------------------------------------------------------------------------------------

//-------------------------------------------------------------------
void InitRoot()
{
    
    gStyle->SetTitleFillColor(kWhite);
    gStyle->SetTitleColor(kBlue,"xyz");
    gStyle->SetTitleSize(0.03,"xyz");
    gStyle->SetLabelSize(0.03,"xyz");
    gStyle->SetLabelColor(kBlue,"xyz");
    gStyle->SetAxisColor(kMagenta);
    gStyle->SetPalette(1,0);
    gStyle->SetOptFit(101);
    gStyle->SetOptStat(0);
    
    //  gStyle->SetOptStat("emr");
    gStyle->SetOptStat("");
}
//-------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
/*
 
 The declaration of this function is now put inside the h file with the default values to be used
 
 void PlotScurve(
 Int_t    kDrawCh=31,
 TString  fname=kFname,
 TString  fname_ped=kFnamePed, //ASIC_D ped_file
 Double_t dac_half_spe_fix = -1,
 TString  fname_gain="gain-org.txt",
 Int_t    fit_min=kFitMin,
 Int_t    fit_max=kFitMax,
 Int_t    fit_min_scurve=kFitMinS,
 Int_t    fit_max_scurve=kFitMaxS,
 Bool_t   kBoard=1, //0:ASIC_C, 1:ASIC_D
 Int_t    kChRef=kRefCh,//ASIC_C
 Bool_t   checkAve=0 //0:use ref pix, 1:check average
 )
 */

// Definition of the function  PlotScurve
void PlotScurve(
                Int_t    kDrawCh,
                TString  fname,
                TString  fname_ped, //ASIC_D ped_file
                Double_t dac_half_spe_fix,
                TString  fname_gain,
                Int_t    fit_min,
                Int_t    fit_max,
                Int_t    fit_min_scurve,
                Int_t    fit_max_scurve,
                Bool_t   kBoard, //0:ASIC_C, 1:ASIC_D
                Int_t    kChRef,//ASIC_C
                Bool_t   checkAve //0:use ref pix, 1:check average
                )


{
    //NOTE
    //Double_t dac_half_spe_fix = 177; //for EC_ASICb N101(asic N158).
    //check DAC value corresponds to voltage to 160fC input
    //(depends on DAC linearity)
    //Double_t dac_half_spe_fix = -1 //set gain to DAC average (default)
    
    /******************* Readout data file *********************/
    
    /*********************** SIGNAL RUN ********************/
    TString *input             = new TString[kNch];
    TString *input_ped         = new TString[kNch];
    Int_t    ch;
    Float_t  dac_half_sum      = 0;
    Float_t  dac_ped_sum       = 0;
    Float_t  dac_spe_sum       = 0;
    Float_t  dac_sfit_sum      = 0;
    Float_t  diff_sum          = 0;
    Float_t  mean_sum          = 0;
    Float_t  maxBinSum         = 0;
    Double_t dac_half[kNch];
    Double_t dac_ped[kNch];
    Double_t dac_spe[kNch];
    Double_t dac_sfit[kNch];
    Double_t mean_diff[kNch];
    Double_t gain_org[kNch];
    Double_t sfit_chiS[kNch];
    Double_t flag_fit[kNch];
    
    
    InitRoot();
    
    TCanvas *c1                = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
    c1->Divide(11,6);
    TCanvas *c2                = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
    c2->Divide(11,6);
    TCanvas *c5                = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
    c5->Divide(11,6);
    TCanvas *c3                = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
    c3->Divide(4,1);
    /*TCanvas *c6                = new TCanvas("c6","c6",450,150,1000,500);
     c6->Divide(4,2);
     TCanvas *c7                = new TCanvas("c7","c7",470,200,1000,500);
     c7->Divide(4,2);*/
    TCanvas *c4                = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
    
    TString  name_hist[kNch];
    TString  name_hist_ped[kNch];
    TString  name_hist_diff[kNch];
    TString  name_hist_diff_fit[kNch];
    TString  name_fit_diff[kNch];
    TH1D    *hist_scurve       = new TH1D[kNch]; // smoothed S-curves
    TH1D    *hist_scurve_ped   = new TH1D[kNch]; // S-curves for pedestal runs
    TH1D    *hist_diff         = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
    TH1D    *hist_diff_fit     = new TH1D[kNch]; // differentiated histograms of fitted s-curves
    TF1     *fit_sdiff         = new TF1[kNch];  // function fitted on differentiated histograms of fitted s-curves
    
    
    // Open the ASIC Gain file
    //-------------------------
    
    TString input_gain(fname_gain);
    ifstream in_gain(input_gain.Data());
    if(!in_gain.good()){
        cerr << "File " << input_gain << " open failed!" << endl;
        return;
    }
    
    
    
    // first loop on channels to fill histonames
    //----------------------
    for(ch=0;ch<kNch;ch++){
        dac_half[ch]=0;
        dac_ped[ch]=0;
        dac_spe[ch]=0;
        dac_sfit[ch]=0;
        gain_org[ch]=0;
        mean_diff[ch]=0;
        sfit_chiS[ch]=0;
        flag_fit[ch]=0;
        name_hist[ch]           ="hist";
        name_hist[ch]          +=ch;
        name_hist_ped[ch]       ="hist_ped";
        name_hist_ped[ch]      +=ch;
        name_hist_diff[ch]      ="hist_diff";
        name_hist_diff[ch]     +=ch;
        name_hist_diff_fit[ch]  ="hist_diff_fit";
        name_hist_diff_fit[ch] +=ch;
        name_fit_diff[ch]       ="fit_diff";
        name_fit_diff[ch]      +=ch;
    }
    
    
    // read gain
    Double_t da;
    Int_t count=0;
    ch=0;
    while(in_gain >> da){
        if(!kBoard && (count>=128 && count<192)){
            gain_org[ch]=da;
            cerr << "!kBoard, count, gain_org[" << ch << "] "
            << !kBoard << " " << count << " " << gain_org[ch] << endl;
            ch++;
        }
        if(kBoard && (count>=192 && count<256)){
            gain_org[ch]=da;
            cerr << "kBoard, count, gain_org[" << ch << "] "
            << kBoard << " " << count << " " << gain_org[ch] << endl;
            ch++;
        }
        count++;
    }// end if first loop on all channels
    in_gain.close();
    in_gain.clear();
    
    //----
    
    // Open the pedestal file and the S-Curve files for each channel and test if the file can be open
    
    Int_t trigger=0, trigger_fit=0, trigger_mean=0, trigger_sfit=0;
    // Big loop on all the channels : second loop
    //--------------------------------------------
    // build the input data file name for each channel fname-chXX.dat
    for(ch=0;ch<kNch;ch++){
        input[ch] += fname;
        input[ch] += "-ch";
        input[ch] += ch;
        input[ch] += ".dat";
        ifstream in(input[ch].Data());
        if(!in.good()){
            cerr << "File " << input[ch] << " open failed!" << endl;
            return;
        }
        
        // build the input data file name for each channel fname_ped-chXX.dat
        input_ped[ch] += fname_ped;
        input_ped[ch] += "-ch";
        input_ped[ch] += ch;
        input_ped[ch] += ".dat";
        ifstream in_ped(input_ped[ch].Data());
        if(!in_ped.good()){
            cerr << "File " << input_ped[ch] << " open failed!" << endl;
            return;
        }
        
        
        
        //10DAC~10mV
        // determine the DAC range in dac_min, dac_max, step_dac from file input_ped
        // and count the number of events
        Double_t da, dat, dac_min, dac_max, step_dac;
        count = 0;
        while (in >> da >> dat ){
            if(count==0) dac_min  = da;
            if(count==1) step_dac = da - dac_min;
            dac_max=da;
            count++;
        }
        in.close();
        in.clear();
        
        /*
         if(ch==43 || ch==47){//for T1 ASIC_C
         cerr << "CHECK " << ch << " " << fit_min_scurve << endl;
         Int_t fit_min_scurve=140;
         } else {
         Int_t fit_min_scurve=kFitMinS;
         }
         */
        
        // correct the scale for DAC range
        if(kXmin<dac_min)kXmin=dac_min;
        const Int_t kNdata    = count;
        count = (fit_max_scurve-fit_min_scurve)/step_dac;
        const Int_t kNdataFit = count;
        if(ch==0){
            cerr << "kNdata, dac_min, dac_max, step_dac: "
            << kNdata << " " << dac_min << " " << dac_max << " "
            << step_dac << endl;
        }
        
        Double_t dac_min_ped, dac_max_ped, step_dac_ped;
        count = 0;
        while (in_ped >> da >> dat ){
            if(count==0) dac_min_ped  = da;
            if(count==1) step_dac_ped = da - dac_min_ped;
            dac_max_ped=da;
            count++;
        }
        
        in_ped.close();
        in_ped.clear();
        
        
        const Int_t kNdataPed = count;
        if(ch==0){
            cerr << "kNdataPed, dac_min_ped, dac_max_ped, step_dac_ped: "
            << kNdataPed << " " << dac_min_ped << " " << dac_max_ped
            << " " << step_dac_ped << endl;
        }
        
        // book histograms
        // SDC: 27/09/2013 hist_scurve is a pointer on an array of THD, thus the syntax is correct (no new operator)
        
        
        hist_scurve[ch]      = TH1D(name_hist[ch],name_hist[ch],kNdata-1,
                                    dac_min,dac_max);
        hist_scurve_ped[ch]  = TH1D(name_hist_ped[ch],name_hist_ped[ch],
                                    kNdataPed-1,dac_min_ped,dac_max_ped);
        hist_diff[ch]        = TH1D(name_hist_diff[ch],name_hist_diff[ch],
                                    kNdata-2,dac_min,dac_max);
        hist_diff_fit[ch]    = TH1D(name_hist_diff_fit[ch],name_hist_diff_fit[ch],
                                    kNdataFit,fit_min_scurve,fit_max_scurve);
        
        
        
        
        
        // read the data for that channel for normal data
        in.open(input[ch].Data());
        Double_t       data[kNdata], dac[kNdata], dac_fit[kNdataFit], data_fit[kNdataFit];
        while (in >> da >> dat ){
            hist_scurve[ch].Fill(da,dat); // fill the histograms of the s-curve
        }
        in.close();
        in.clear();
        
        // read the data for pedestal data
        in_ped.open(input_ped[ch].Data());
        while (in_ped >> da >> dat ){
            hist_scurve_ped[ch].Fill(da,dat); // fill the histogram for pedestal data
        }
        in_ped.close();
        in_ped.clear();
        
        //dac_ped[ch]  = hist_scurve_ped[ch].GetMaximumBin()+dac_min_ped; //absolute value
        dac_ped[ch]  = hist_scurve_ped[ch].GetBinCenter(hist_scurve_ped[ch].GetMaximumBin()); //absolute value
        dac_ped_sum += dac_ped[ch]; //absolute value
        //cerr << "dac_ped[" << ch << "] dac_ped_sum " << dac_ped[ch]
        //<< " " << dac_ped_sum << endl;
        
        
        hist_scurve[ch].Smooth(3); // Smooth the S curve
        Int_t binx=0;
        //hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
        //                                  dac_max-dac_min+1,kHalfMax/2.);
        
        // ??????
        hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
                                          dac_max-dac_min+1,2000);
        dac_half[ch] = hist_scurve[ch].GetBinCenter(binx); //absolute value
        if(dac_half[ch]<kXmin)dac_half[ch]=-1;
        
        // put the histograms of the smoothed S-Curves into two arrays (dac, data)
        //-------------------------------------------------------------------------
        Int_t i=0;
        for(i=0;i<kNdata;i++){
            dac[i]  = hist_scurve[ch].GetBinCenter(i+1);
            data[i] = hist_scurve[ch].GetBinContent(i+1);
        }
        
        /*** compute derive diff on S-Curve ***/
        Double_t diff[kNdata-1], dac_diff[kNdata-1];
        Int_t trigger_diff=0;
        for(i=0;i<kNdata-1;i++){
            diff[i]     = (data[i]-data[i+1])/step_dac;
            dac_diff[i] = dac[i]+step_dac/2;
            hist_diff[ch].Fill(dac_diff[i],diff[i]); // fills differentiated histograms of smoothed s-curves which is the single photoelectron spectrum
            if(dac_diff[i]>=fit_min){
                diff_sum += diff[i]*dac_diff[i];
                //cerr << "diff[" << i << "], dac_diff[" << i << "] "
                //<< diff[i] << " " << dac_diff[i] << endl;
                trigger_diff+=diff[i];
            }
        }
        //mean_diff[ch] = diff_sum/trigger_diff;
        hist_diff[ch].GetXaxis()->SetRangeUser(fit_min,fit_max);
        mean_diff[ch] = hist_diff[ch].GetMean();
        //cerr << "trigger_diff, mean_diff[" << ch << "] "
        //<< trigger_diff << " " << mean_diff[ch] << endl;
        
        
        /******************* GRAPHICS *****************/
        
        
        // Now start the analysis which consist in finding the peak in the differentiated spectrum
        c1->cd(ch+1);
        /*** fit scurves  ***/
        cerr << "ch " << ch;
        
        //    SDC 24/09/2013 :: Need to declare outside this variable in a compiled program
        TF1 *fit_pol ;
        
        
        // Fit a polynomial curve in a given range on the s-curve directly. This is a kind of DEBUG procedure.
        
        //if(ch==43 || ch==47){ // SDC: 25/09/2013 Why a special channel : the fit pol4 is performed for another DAC range
        /*   if(ch==170){
         fit_pol      = new TF1("fit_pol","pol4",150,fit_max_scurve);
         hist_scurve[ch].Fit(fit_pol,"","",150,fit_max_scurve);
         }else{ */
        fit_pol      = new TF1("fit_pol","pol4",fit_min_scurve,fit_max_scurve);
        hist_scurve[ch].Fit(fit_pol,"","",fit_min_scurve,fit_max_scurve);
        // }
        
        
        
        Double_t fit_par_pol4[5]; // container for fit parameter results
        fit_pol->GetParameters(fit_par_pol4);
        sfit_chiS[ch] =  fit_pol->GetChisquare(); // chi squared
        TF1 *func               = hist_scurve[ch].GetFunction("fit_pol");
        Double_t data_diff_fit;
        for(i=0;i<kNdataFit;i++){
            //if(ch==43 || ch==47){
            if(ch==170){
                dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+150-dac_min;
                data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+150-dac_min); // get the value of the fitted function
            }else{
                dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min;
                data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min); // get the value of the fitted function
                //cerr << "dac_fit[" << i << "], data_fit[" << i << "] "
                //<< dac_fit[i] << " " << data_fit[i] << endl;
            }
        }
        
        // derive the fitted function on the S-curve to obtain the SPE
        Double_t dac_diff_fit;
        for(i=0;i<kNdataFit-1;i++){
            data_diff_fit = (data_fit[i]-data_fit[i+1])/step_dac;
            hist_diff_fit[ch].Fill(dac_fit[i],data_diff_fit); // Fill differentiated fitted S-curve
        }
        
        TLine *lhalf_max   = new TLine(dac[0],kHalfMax,dac[kNdata-1],kHalfMax);
        TLine *ldac_half   = new TLine(dac_half[ch],0,dac_half[ch],
                                       (Double_t)kDrawMax);
        hist_scurve[ch].SetMaximum(kDrawMax);
        hist_scurve[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
        hist_scurve[ch].Draw();
        lhalf_max->Draw("same");
        ldac_half->Draw("same");
        if(dac_half[ch]>0){
            dac_half_sum += dac_half[ch]-dac_ped[ch]; //relative absolute value
            trigger++;
        }
        
        /*** Differential plot from scurve fitting function  ***/
        c5->cd(ch+1);
        hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve,fit_max_scurve-2*step_dac);
        Double_t dac_sfit_max=hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
        //if(dac_sfit_max<=fit_min_scurve+4*step_dac){
        if(dac_sfit_max<=fit_min_scurve+2*step_dac){
            hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve+kFitSadj,fit_max_scurve-2*step_dac);
            dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
            //if(dac_sfit[ch]<=fit_min_scurve+kFitSadj+4*step_dac)dac_sfit[ch]=-1;
            if(dac_sfit[ch]<fit_min_scurve+kFitSadj+step_dac){
                dac_sfit[ch]=-1;
                //dac_sfit[ch]=dac_half[ch];
                flag_fit[ch] = -1;
                cerr << "flag_fit[" << ch << "] " << flag_fit[ch]  << endl;
            }
        }else {
            dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
        }
        
        if(dac_sfit[ch]>fit_min_scurve+5){
            dac_sfit_sum += dac_sfit[ch]-dac_ped[ch];
            trigger_sfit++;
        }else{
            //dac_sfit[ch] = fit_min_scurve+10;
            dac_sfit[ch] = -1;
        }
        //cerr << "dac_sfit[" << ch << "] " << dac_sfit[ch] << endl;
        TLine *ldac_sfit      = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMaxFit);
        TLine *ldac_sfit_copy = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMax);
        ldac_sfit->SetLineColor(kMagenta);
        ldac_sfit_copy->SetLineColor(kMagenta);
        hist_diff_fit[ch].SetMaximum(kDrawMaxFit);
        hist_diff_fit[ch].SetMinimum(0);
        hist_diff_fit[ch].Draw();
        ldac_sfit->Draw("same");
        
        /**** differential plot ****/
        c2->cd(ch+1);
        //hist_diff[ch].Rebin(8);
        hist_diff[ch].Rebin(3);
        fit_sdiff     = new TF1("fit_sdiff","gaus",dac_min,dac_max);
        //fit_sdiff->SetLineWidth(1);
        //if(ch==8 || ch==18 || ch==44 || ch==56)
        mean_sum += mean_diff[ch]-dac_ped[ch];
        trigger_mean++;
        
        // fit the derivated polynom fit
        Double_t fit_par[3];
        //hist_diff[ch].Fit(fit_sdiff,"RQ0","",fit_min,fit_max);
        hist_diff[ch].Fit(fit_sdiff,"","",fit_min,fit_max);
        hist_diff[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
        hist_diff[ch].SetMinimum(0);
        hist_diff[ch].SetMaximum(kDrawMaxDiff);
        //hist_diff[ch].Rebin();
        hist_diff[ch].Draw();
        fit_sdiff->GetParameters(fit_par); // get the parameters for this fit
        dac_spe[ch] = fit_par[1];
        if(dac_spe[ch]>fit_min){
            dac_spe_sum += dac_spe[ch]-dac_ped[ch];
            //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
            trigger_fit++;
        }else{
            dac_spe[ch] = fit_min+10;
        }
        //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
        TLine *ldac_spe    = new TLine(dac_spe[ch],0,dac_spe[ch],
                                       (Double_t)kDrawMaxDiff);
        ldac_spe->SetLineColor(kRed);
        TLine *lmean_diff  = new TLine(mean_diff[ch],0,mean_diff[ch],
                                       (Double_t)kDrawMaxDiff);
        lmean_diff->SetLineColor(kCyan);
        TLine *lfit_min    = new TLine(fit_min,0,fit_min,kDrawMaxDiff);
        lfit_min->SetLineColor(kGray);
        
        ldac_spe->Draw();
        lmean_diff->Draw("same");
        lfit_min->Draw("same");
        
        c3->cd(1);
        if(ch==0){
            hist_scurve[ch].SetMaximum(kDrawMax);
            hist_scurve[ch].Draw();
        }else{
            /*if(ch==8 || ch==18 || ch==44 || ch==56)*/
            hist_scurve[ch].Draw("same");
        }
        if(ch==kDrawCh){
            c3->cd(2);
            hist_scurve[ch].SetMaximum(kDrawMax);
            hist_scurve[kDrawCh].Draw();
            lhalf_max->Draw("same");
            ldac_half->Draw("same");
            ldac_sfit_copy->Draw("same");
            c3->cd(3);
            hist_diff[kDrawCh].Draw();
            ldac_spe->Draw();
            lmean_diff->Draw("same");
            lfit_min->Draw("same");
            c3->cd(4);
            hist_diff_fit[kDrawCh].Draw();
            ldac_sfit->Draw("same");
        }
        /*
         if(ch<32){
         c6->cd(ch/4+1);
         hist_scurve[ch].SetMaximum(kDrawMax);
         if(ch%4==0){
         hist_scurve[ch].Draw();
         }else{
         hist_scurve[ch].Draw("same");
         }
         ldac_sfit_copy->Draw("same");
         }else{
         c7->cd(ch/4-7);
         hist_scurve[ch].SetMaximum(kDrawMax);
         if(ch%4==0){
         hist_scurve[ch].Draw();
         }else{
         hist_scurve[ch].Draw("same");
         }
         ldac_sfit_copy->Draw("same");
         }
         */
        
        c4->cd();
        if(ch==0){
            hist_scurve_ped[ch].SetMaximum(kDrawMaxPed);
            hist_scurve_ped[ch].Draw("l");
        }else{
            hist_scurve_ped[ch].Draw("lsame");
        }
        
    }// end of second loop on all the channels
    
    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
    
    Float_t dac_ped_ave  ; //absolute value
    Float_t dac_half_ave  ;
    Float_t dac_spe_ave   ;
    Float_t dac_sfit_ave;
    Float_t mean_diff_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
    
    if(!checkAve){
        dac_ped_ave   = dac_ped[kChRef]; //absolute value
        dac_half_ave  = dac_half[kChRef];
        dac_spe_ave   = dac_spe[kChRef];
        dac_sfit_ave  = dac_sfit[kChRef];
        //Float_t dac_sfit_ave  = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
        mean_diff_ave = mean_diff[kChRef];
    }else{
        dac_ped_ave  = dac_ped_sum/(Float_t)kNch; //absolute value
        dac_half_ave = dac_half_sum/(Float_t)trigger+dac_ped_ave;
        dac_spe_ave  = dac_spe_sum/(Float_t)trigger_fit+dac_ped_ave;
        mean_diff_ave = mean_sum/(Float_t)trigger_mean+dac_ped_ave;
        dac_sfit_ave = dac_sfit_sum/(Float_t)trigger_sfit+dac_ped_ave;
    }
    
    
    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
    
    Double_t dac_half_spe;
    
    //trigger_mean=4;
    if(dac_half_spe_fix!=-1) {
        dac_half_spe = dac_half_spe_fix;
    } else {
        //Double_t dac_half_spe = dac_half_ave;
        dac_half_spe = dac_spe_ave;
    }
    
    TLine *ldac_ped_ave   = new TLine(dac_ped_ave,0,dac_ped_ave,
                                      (Double_t)kDrawMax);
    TLine *ldac_ped_ave_copy = new TLine(dac_ped_ave,0,dac_ped_ave,
                                         (Double_t)kDrawMaxDiff);
    TLine *ldac_half_ave  = new TLine(dac_half_ave,0,dac_half_ave,
                                      (Double_t)kDrawMax);
    TLine *ldac_half_spe  = new TLine(dac_half_spe,0,dac_half_spe,
                                      (Double_t)kDrawMax);
    TLine *lmean_diff_ave = new TLine(mean_diff_ave,0,mean_diff_ave,
                                      (Double_t)kDrawMax);
    TLine *ldac_sfit_ave  = new TLine(dac_sfit_ave,0,dac_sfit_ave,
                                      (Double_t)kDrawMax);
    ldac_ped_ave->SetLineColor(kOrange);
    ldac_half_ave->SetLineColor(kRed);
    lmean_diff_ave->SetLineColor(kCyan);
    //ldac_half_spe->SetLineColor(kGreen);
    ldac_sfit_ave->SetLineColor(kMagenta);
    
    c3->cd(1);
    ldac_ped_ave->Draw("same");
    ldac_half_ave->Draw("same");
    // lhalf_max->Draw("same");  SDC (24/09/2013) :: Not declared
    ldac_half_spe->Draw("same");
    lmean_diff_ave->Draw("same");
    ldac_sfit_ave->Draw("same");
    
    c3->cd(3);
    ldac_ped_ave_copy->SetLineColor(kOrange);
    ldac_ped_ave_copy->Draw("same");
    
    ofstream out(kFnameOut.Data());
    
    Int_t    gain_asic_diff, gain_asic_dacHalf, gain_asic_diff_mean, gain_asic_sfit;
    out << "title  ch  g_asic_dacHalf  g_asic_diff  g_asic_diff_mean  g_asic_sfit  dac_half[ch]  dac_ped[ch]  dac_spe[ch]  dac_half_ave dac_half_spe dac_ped_ave  dac_sfit_ave"
    << endl;
    
    
    // loop 3 on all channels to make final computations on the gain
    for(ch=0;ch<kNch;ch++){
        /*
         if(ch==8 || ch==44 || ch==56 || ch==18){
         gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
         (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
         }else{
         gain_asic_diff_mean=0;
         }
         */
        if(dac_half[ch]<0){
            gain_asic_dacHalf   = 0;
            gain_asic_diff      = 0;
            gain_asic_diff_mean = 0;
            gain_asic_sfit      = 0;
        }else{
            gain_asic_dacHalf  =(Int_t)((dac_half_ave-dac_ped_ave)/
                                        (dac_half[ch]-dac_ped[ch])*gain_org[ch]);
            gain_asic_diff     =(Int_t)((dac_half_spe-dac_ped_ave)/
                                        (dac_spe[ch]-dac_ped[ch])*gain_org[ch]);
            gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
                                        (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
            gain_asic_sfit     =(Int_t)((dac_sfit_ave-dac_ped_ave)/
                                        (dac_sfit[ch]-dac_ped[ch])*gain_org[ch]);
            //cerr << ch << " " << gain_org[ch] << endl;
        }
        
        
        
        out << "OUT_TABLE    "
        << ch << "    "
        << gain_asic_dacHalf   << "    "
        << gain_asic_diff      << "    "
        << gain_asic_diff_mean << "    "
        << gain_asic_sfit      << "    "
        << dac_half[ch]        << "    "
        << dac_ped[ch]         << "    "
        << (Int_t)dac_spe[ch]  << "    "
        << dac_half_ave        << "    "
        << dac_half_spe        << "    "
        << dac_ped_ave         << "    "
        << dac_sfit_ave        << "    "
        << sfit_chiS[ch]       << "    "
        << ch                  << "    "
        << flag_fit[ch]
        << endl;
        std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
    }
    
    
    
    /*
     for(ch=kNch-1;ch>=0;ch--){
     cout << "OUT_TABLE_INV    "
     << (Int_t)(dac_half_ave/(dac_half[ch]-dac_ped)*64)
     << endl;
     }
     */
    
    
}

//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
// Analysis Constructor
//--------------------------------------------------------------------------------------------------------
Analysis::Analysis():ch(0),dac_half_sum(0),dac_spe_sum(0),dac_sfit_sum(0),diff_sum(0),mean_sum(0),maxBinSum(0)
{
    if(DEBUG)
        cout << " >>>>>>> START Analysis Constructor " << endl;
    
    /*********************** SIGNAL RUN ********************/
    input             = new TString[kNch];
    input_ped         = new TString[kNch];
    
    // Init sums
    dac_half_sum      = 0;
    dac_ped_sum       = 0;
    dac_spe_sum       = 0;
    dac_sfit_sum      = 0;
    diff_sum          = 0;
    mean_sum          = 0;
    maxBinSum         = 0;
    
    if(DEBUG)
        cout << " >>>>>>> END Analysis Constructor " << endl;
}
//--------------------------------------------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
// Analysis destructor
//--------------------------------------------------------------------------------------------------------
Analysis::~Analysis()
{
    
    if(DEBUG)
        cout << " >>>>>>> START Analysis Destructor " << endl;
    
    if(input!=0) delete input;
    if(input_ped!=0) delete input_ped;
    
    if(DEBUG)
        cout << " >>>>>>> END Analysis Destructor " << endl;
    
}
//--------------------------------------------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
// Main public function of Analysis
//--------------------------------------------------------------------------------------------------------
/*
 void Analysis::PlotScurve(
 Int_t    kDrawCh=31,
 TString  fname=kFname,
 TString  fname_ped=kFnamePed, //ASIC_D ped_file
 Double_t dac_half_spe_fix = -1,
 TString  fname_gain="gain-org.txt",
 Int_t    fit_min=kFitMin,
 Int_t    fit_max=kFitMax,
 Int_t    fit_min_scurve=kFitMinS,
 Int_t    fit_max_scurve=kFitMaxS,
 Bool_t   kBoard=0, //0:ASIC_C, 1:ASIC_D
 Int_t    kChRef=kRefCh,//ASIC_C
 Bool_t   checkAve=0 //0:use ref pix, 1:check average
 )
 */
void Analysis:: PlotScurve(
                           Int_t    kDrawCh,
                           TString  fname,
                           TString  fname_ped, //ASIC_D ped_file
                           Double_t dac_half_spe_fix,
                           TString  fname_gain,
                           Int_t    fit_min,
                           Int_t    fit_max,
                           Int_t    fit_min_scurve,
                           Int_t    fit_max_scurve,
                           Bool_t   kBoard, //0:ASIC_C, 1:ASIC_D
                           Int_t    kChRef,//ASIC_C
                           Bool_t   checkAve //0:use ref pix, 1:check average
                           )



{
    
    if(DEBUG)
    {
        cout << " >>>>>>> START Analysis::PlotScurve " << endl;
        cout << " \t - kDrawCh = " << kDrawCh << endl;
        cout << " \t - fname = " << fname << endl;
        cout << " \t - fname_ped = " << fname_ped <<endl;
        cout << " \t - dac_half_spe_fix = " << dac_half_spe_fix << endl;
        
        cout << " \t - fname_gain = " << fname_gain << endl;
        cout << " \t - fit_min = " << fit_min << endl;
        cout << " \t - fit_max = " << fit_max << endl;
        cout << " \t - fit_min_scurve = " << fit_min_scurve << endl;
        cout << " \t - fit_max_scurve = " << fit_max_scurve << endl;
        cout << " \t - kBoard = " << kBoard << " 0:ASIC_C, 1:ASIC_D " << endl;
        cout << " \t - kChRef = " << kChRef << " ASIC_C " << endl;
        cout << " \t - checkAve = " << checkAve << " 0:use ref pix, 1:check average " << endl;
    }
    // reset all variables
    Init();
    
    //OpenGainFile(fname_gain);
    
    FillHistoName();
    
    ReadGainFile(fname_gain,kBoard);
    
    // Open the pedestal file and the S-Curve files for each channel and test if the file can be open
    
    trigger=0;
    trigger_fit=0;
    trigger_mean=0;
    trigger_sfit=0;
    
    // Big loop on all the channels : second loop
    //--------------------------------------------
    // build the input data file name for each channel fname-chXX.dat
    for(ch=0;ch<kNch;ch++)
    {
        ProcessSingleChannel(kDrawCh,fname,fname_ped,dac_half_spe_fix,fname_gain,
                             fit_min,fit_max,fit_min_scurve,fit_max_scurve,kBoard,kChRef);
    }
    
    
    
    ComputeDACAverage(dac_half_spe_fix ,kChRef,checkAve);
    OutputGain(kFnameOut);
    
    if(DEBUG)
        cout << " >>>>>>> END Analysis::PlotSCurve " << endl;
    
}
//--------------------------------------------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
// Initialize the program by reseting counters or variables and create objects
//--------------------------------------------------------------------------------------------------------
void Analysis::Init()
{
    
    if(DEBUG)
        cout << " >>>>>>> START Analysis::Init " << endl;
    
    //if(input!=0) delete input;
    //if(input_ped!=0) delete input_ped;
    
    // recreate the name
    input             = new TString[kNch];
    input_ped         = new TString[kNch];
    
    // Init sums
    dac_half_sum      = 0;
    dac_ped_sum       = 0;
    dac_spe_sum       = 0;
    dac_sfit_sum      = 0;
    diff_sum          = 0;
    mean_sum          = 0;
    maxBinSum         = 0;
    
    
    // create canvas
    
    c1                = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
    c1->Divide(11,6);
    c2                = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
    c2->Divide(11,6);
    c5                = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
    c5->Divide(11,6);
    c3                = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
    c3->Divide(4,1);
    c4                = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
    
    
    // create histos and function
    hist_scurve       = new TH1D[kNch]; // smoothed S-curves
    hist_scurve_ped   = new TH1D[kNch]; // S-curves for pedestal runs
    hist_diff         = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
    hist_diff_fit     = new TH1D[kNch]; // differentiated histograms of fitted s-curves
    fit_sdiff         = new TF1[kNch];  // function fitted on differentiated histograms of fitted s-curves
    
    
    for(ch=0;ch<kNch;ch++){
        dac_half[ch]=0;
        dac_ped[ch]=0;
        dac_spe[ch]=0;
        dac_sfit[ch]=0;
        gain_org[ch]=0;
        mean_diff[ch]=0;
        sfit_chiS[ch]=0;
        flag_fit[ch]=0;
    }
    
    if(DEBUG)
        cout << " >>>>>>> END Analysis::Init " << endl;
    
    
}
//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
// Initialize the program by reseting counters or variables and create objects
//--------------------------------------------------------------------------------------------------------
void Analysis::OpenGainFile(const TString fname_gain)
{
    if(DEBUG)
        cout << " >>>>>>> START Analysis::OpenGainFile " << endl;
    
    
    
    TString input_gain(fname_gain);
    ifstream in_gain(input_gain.Data());
    if(!in_gain.good()){
        cerr << "File " << input_gain << " open failed!" << endl;
        exit(-1);
    }
    in_gain.close();
    in_gain.clear();
    
    if(DEBUG)
        cout << " >>>>>>> END Analysis::OpenGainFile " << endl;
    
    
}
//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
// Fill Histogram name
//--------------------------------------------------------------------------------------------------------
void Analysis::FillHistoName()
{
    
    if(DEBUG)
        cout << " >>>>>>> START Analysis::FillHistoName " << endl;
    
    
    
    // first loop on channels to fill histonames
    //----------------------
    for(ch=0;ch<kNch;ch++){
        
        name_hist[ch]           ="hist";
        name_hist[ch]          +=ch;
        name_hist_ped[ch]       ="hist_ped";
        name_hist_ped[ch]      +=ch;
        name_hist_diff[ch]      ="hist_diff";
        name_hist_diff[ch]     +=ch;
        name_hist_diff_fit[ch]  ="hist_diff_fit";
        name_hist_diff_fit[ch] +=ch;
        name_fit_diff[ch]       ="fit_diff";
        name_fit_diff[ch]      +=ch;
    }
    
    if(DEBUG)
        cout << " >>>>>>> END Analysis::FillHistoName " << endl;
    
    
    
}
//--------------------------------------------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
// Read the gain
// input : the gain filename
// output the gain_org[ch] array
//--------------------------------------------------------------------------------------------------------

void Analysis::ReadGainFile(const TString fname_gain,Bool_t kBoard)
{
    
    if(DEBUG)
    {
        cout << " >>>>>>> START Analysis::ReadGainFile " << endl;
        cout << " \t - fname_gain = " << fname_gain << endl;
        cout << " \t - kBoard = " << kBoard << " 0:ASIC_C, 1:ASIC_D " << endl;
        
    }
    
    
    
    TString input_gain(fname_gain);
    ifstream in_gain(input_gain.Data());
    if(!in_gain.good()){
        cerr << "File " << input_gain << " open failed!" << endl;
        exit(-1);
    }
    // read gain
    Double_t da;
    Int_t count=0;
    ch=0;
    while(in_gain >> da){
        if(!kBoard && (count>=128 && count<192)){
            gain_org[ch]=da;
            cerr << "!kBoard, count, gain_org[" << ch << "] "
            << !kBoard << " " << count << " " << gain_org[ch] << endl;
            ch++;
        }
        if(kBoard && (count>=192 && count<256)){
            gain_org[ch]=da;
            cerr << "kBoard, count, gain_org[" << ch << "] "
            << kBoard << " " << count << " " << gain_org[ch] << endl;
            ch++;
        }
        count++;
    }// end if first loop on all channels
    in_gain.close();
    in_gain.clear();
    
    if(DEBUG)
    {
        cout << " >>>>>>> END Analysis::ReadGainFile " << endl;
        
    }
    
}
//--------------------------------------------------------------------------------------------------------



//--------------------------------------------------------------------------------------------------------
// ProcessSingleChannel
//--------------------------------------------------------------------------------------------------------
void Analysis::ProcessSingleChannel(
                                    Int_t    kDrawCh,
                                    TString  fname,
                                    TString  fname_ped, //ASIC_D ped_file
                                    Double_t dac_half_spe_fix,
                                    TString  fname_gain,
                                    Int_t    fit_min,
                                    Int_t    fit_max,
                                    Int_t    fit_min_scurve,
                                    Int_t    fit_max_scurve,
                                    Bool_t   kBoard, //0:ASIC_C, 1:ASIC_D
                                    Int_t    kChRef//ASIC_C
                                    //	Bool_t   checkAve //0:use ref pix, 1:check average)
                                    )
{
    
    
    if(DEBUG)
    {
        cout << " >>>>>>> START Analysis::ProcessSingleChannel ch=" << ch << endl;
        cout << " \t - kDrawCh = " << kDrawCh << endl;
        cout << " \t - fname = " << fname << endl;
        cout << " \t - fname_ped = " << fname_ped <<endl;
        cout << " \t - dac_half_spe_fix = " << dac_half_spe_fix << endl;
        
        cout << " \t - fname_gain = " << fname_gain << endl;
        cout << " \t - fit_min = " << fit_min << endl;
        cout << " \t - fit_max = " << fit_max << endl;
        cout << " \t - fit_min_scurve = " << fit_min_scurve << endl;
        cout << " \t - fit_max_scurve = " << fit_max_scurve << endl;
        cout << " \t - kBoard = " << kBoard << " 0:ASIC_C, 1:ASIC_D " << endl;
        cout << " \t - kChRef = " << kChRef << " ASIC_C " << endl;
        //    cout << " \t - checkAve = " << checkAve << " 0:use ref pix, 1:check average " << endl;L
    }
    
    
    
    input[ch] += fname;
    input[ch] += "-ch";
    input[ch] += ch;
    input[ch] += ".dat";
    ifstream in(input[ch].Data());
    if(!in.good()){
        cerr << "File " << input[ch] << " open failed!" << endl;
        return;
    }
    
    // build the input data file name for each channel fname_ped-chXX.dat
    input_ped[ch] += fname_ped;
    input_ped[ch] += "-ch";
    input_ped[ch] += ch;
    input_ped[ch] += ".dat";
    ifstream in_ped(input_ped[ch].Data());
    if(!in_ped.good()){
        cerr << "File " << input_ped[ch] << " open failed!" << endl;
        return;
    }
    
    
    
    //10DAC~10mV
    // determine the DAC range in dac_min, dac_max, step_dac from file input_ped
    // and count the number of events
    Double_t da, dat, dac_min, dac_max, step_dac;
    Int_t count = 0;
    while (in >> da >> dat ){
        if(count==0) dac_min  = da;
        if(count==1) step_dac = da - dac_min;
        dac_max=da;
        count++;
    }
    in.close();
    in.clear();
    
    /*
     if(ch==43 || ch==47){//for T1 ASIC_C
     cerr << "CHECK " << ch << " " << fit_min_scurve << endl;
     Int_t fit_min_scurve=140;
     } else {
     Int_t fit_min_scurve=kFitMinS;
     }
     */
    
    // correct the scale for DAC range
    if(kXmin<dac_min)kXmin=dac_min;
    const Int_t kNdata    = count;
    count = (fit_max_scurve-fit_min_scurve)/step_dac;
    const Int_t kNdataFit = count;
    if(ch==0){
        cerr << "kNdata, dac_min, dac_max, step_dac: "
        << kNdata << " " << dac_min << " " << dac_max << " "
        << step_dac << endl;
    }
    
    Double_t dac_min_ped, dac_max_ped, step_dac_ped;
    count = 0;
    while (in_ped >> da >> dat ){
        if(count==0) dac_min_ped  = da;
        if(count==1) step_dac_ped = da - dac_min_ped;
        dac_max_ped=da;
        count++;
    }
    
    in_ped.close();
    in_ped.clear();
    
    
    const Int_t kNdataPed = count;
    if(ch==0){
        cerr << "kNdataPed, dac_min_ped, dac_max_ped, step_dac_ped: "
        << kNdataPed << " " << dac_min_ped << " " << dac_max_ped
        << " " << step_dac_ped << endl;
    }
    
    // book histograms
    // SDC: 27/09/2013 hist_scurve is a pointer on an array of THD, thus the syntax is correct (no new operator)
    
    
    hist_scurve[ch]      = TH1D(name_hist[ch],name_hist[ch],kNdata-1,
                                dac_min,dac_max);
    hist_scurve_ped[ch]  = TH1D(name_hist_ped[ch],name_hist_ped[ch],
                                kNdataPed-1,dac_min_ped,dac_max_ped);
    hist_diff[ch]        = TH1D(name_hist_diff[ch],name_hist_diff[ch],
                                kNdata-2,dac_min,dac_max);
    hist_diff_fit[ch]    = TH1D(name_hist_diff_fit[ch],name_hist_diff_fit[ch],
                                kNdataFit,fit_min_scurve,fit_max_scurve);
    
    
    
    
    
    // read the data for that channel for normal data
    in.open(input[ch].Data());
    Double_t       data[kNdata], dac[kNdata], dac_fit[kNdataFit], data_fit[kNdataFit];
    while (in >> da >> dat ){
        hist_scurve[ch].Fill(da,dat); // fill the histograms of the s-curve
    }
    in.close();
    in.clear();
    
    // read the data for pedestal data
    in_ped.open(input_ped[ch].Data());
    while (in_ped >> da >> dat ){
        hist_scurve_ped[ch].Fill(da,dat); // fill the histogram for pedestal data
    }
    in_ped.close();
    in_ped.clear();
    
    //dac_ped[ch]  = hist_scurve_ped[ch].GetMaximumBin()+dac_min_ped; //absolute value
    dac_ped[ch]  = hist_scurve_ped[ch].GetBinCenter(hist_scurve_ped[ch].GetMaximumBin()); //absolute value
    dac_ped_sum += dac_ped[ch]; //absolute value
    //cerr << "dac_ped[" << ch << "] dac_ped_sum " << dac_ped[ch]
    //<< " " << dac_ped_sum << endl;
    
    
    hist_scurve[ch].Smooth(3); // Smooth the S curve
    Int_t binx=0;
    //hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
    //                                  dac_max-dac_min+1,kHalfMax/2.);
    
    // ??????
    hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
                                      dac_max-dac_min+1,2000);
    dac_half[ch] = hist_scurve[ch].GetBinCenter(binx); //absolute value
    if(dac_half[ch]<kXmin)dac_half[ch]=-1;
    
    // put the histograms of the smoothed S-Curves into two arrays (dac, data)
    //-------------------------------------------------------------------------
    Int_t i=0;
    for(i=0;i<kNdata;i++){
        dac[i]  = hist_scurve[ch].GetBinCenter(i+1);
        data[i] = hist_scurve[ch].GetBinContent(i+1);
    }
    
    /*** compute derive diff on S-Curve ***/
    Double_t diff[kNdata-1], dac_diff[kNdata-1];
    Int_t trigger_diff=0;
    for(i=0;i<kNdata-1;i++){
        diff[i]     = (data[i]-data[i+1])/step_dac;
        dac_diff[i] = dac[i]+step_dac/2;
        hist_diff[ch].Fill(dac_diff[i],diff[i]); // fills differentiated histograms of smoothed s-curves which is the single photoelectron spectrum
        if(dac_diff[i]>=fit_min){
            diff_sum += diff[i]*dac_diff[i];
            //cerr << "diff[" << i << "], dac_diff[" << i << "] "
            //<< diff[i] << " " << dac_diff[i] << endl;
            trigger_diff+=diff[i];
        }
    }
    //mean_diff[ch] = diff_sum/trigger_diff;
    hist_diff[ch].GetXaxis()->SetRangeUser(fit_min,fit_max);
    mean_diff[ch] = hist_diff[ch].GetMean();
    //cerr << "trigger_diff, mean_diff[" << ch << "] "
    //<< trigger_diff << " " << mean_diff[ch] << endl;
    
    
    /******************* GRAPHICS *****************/
    
    
    // Now start the analysis which consist in finding the peak in the differentiated spectrum
    c1->cd(ch+1);
    /*** fit scurves  ***/
    cerr << "ch " << ch;
    
    //    SDC 24/09/2013 :: Need to declare outside this variable in a compiled program
    TF1 *fit_pol ;
    
    
    // Fit a polynomial curve in a given range on the s-curve directly. This is a kind of DEBUG procedure.
    
    //if(ch==43 || ch==47){ // SDC: 25/09/2013 Why a special channel : the fit pol4 is performed for another DAC range
    /*   if(ch==170){
     fit_pol      = new TF1("fit_pol","pol4",150,fit_max_scurve);
     hist_scurve[ch].Fit(fit_pol,"","",150,fit_max_scurve);
     }else{ */
    fit_pol      = new TF1("fit_pol","pol4",fit_min_scurve,fit_max_scurve);
    hist_scurve[ch].Fit(fit_pol,"","",fit_min_scurve,fit_max_scurve);
    // }
    
    
    
    Double_t fit_par_pol4[5]; // container for fit parameter results
    fit_pol->GetParameters(fit_par_pol4);
    sfit_chiS[ch] =  fit_pol->GetChisquare(); // chi squared
    TF1 *func               = hist_scurve[ch].GetFunction("fit_pol");
    Double_t data_diff_fit;
    for(i=0;i<kNdataFit;i++){
        //if(ch==43 || ch==47){
        if(ch==170){
            dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+150-dac_min;
            data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+150-dac_min); // get the value of the fitted function
        }else{
            dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min;
            data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min); // get the value of the fitted function
            //cerr << "dac_fit[" << i << "], data_fit[" << i << "] "
            //<< dac_fit[i] << " " << data_fit[i] << endl;
        }
    }
    
    // derive the fitted function on the S-curve to obtain the SPE
    Double_t dac_diff_fit;
    for(i=0;i<kNdataFit-1;i++){
        data_diff_fit = (data_fit[i]-data_fit[i+1])/step_dac;
        hist_diff_fit[ch].Fill(dac_fit[i],data_diff_fit); // Fill differentiated fitted S-curve
    }
    
    TLine *lhalf_max   = new TLine(dac[0],kHalfMax,dac[kNdata-1],kHalfMax);
    TLine *ldac_half   = new TLine(dac_half[ch],0,dac_half[ch],
                                   (Double_t)kDrawMax);
    hist_scurve[ch].SetMaximum(kDrawMax);
    hist_scurve[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
    hist_scurve[ch].Draw();
    lhalf_max->Draw("same");
    ldac_half->Draw("same");
    if(dac_half[ch]>0){
        dac_half_sum += dac_half[ch]-dac_ped[ch]; //relative absolute value
        trigger++;
    }
    
    /*** Differential plot from scurve fitting function  ***/
    c5->cd(ch+1);
    hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve,fit_max_scurve-2*step_dac);
    Double_t dac_sfit_max=hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
    //if(dac_sfit_max<=fit_min_scurve+4*step_dac){
    if(dac_sfit_max<=fit_min_scurve+2*step_dac){
        hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve+kFitSadj,fit_max_scurve-2*step_dac);
        dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
        //if(dac_sfit[ch]<=fit_min_scurve+kFitSadj+4*step_dac)dac_sfit[ch]=-1;
        if(dac_sfit[ch]<fit_min_scurve+kFitSadj+step_dac){
            dac_sfit[ch]=-1;
            //dac_sfit[ch]=dac_half[ch];
            flag_fit[ch] = -1;
            cerr << "flag_fit[" << ch << "] " << flag_fit[ch]  << endl;
        }
    }else {
        dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
    }
    
    if(dac_sfit[ch]>fit_min_scurve+5){
        dac_sfit_sum += dac_sfit[ch]-dac_ped[ch];
        trigger_sfit++;
    }else{
        //dac_sfit[ch] = fit_min_scurve+10;
        dac_sfit[ch] = -1;
    }
    //cerr << "dac_sfit[" << ch << "] " << dac_sfit[ch] << endl;
    TLine *ldac_sfit      = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMaxFit);
    TLine *ldac_sfit_copy = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMax);
    ldac_sfit->SetLineColor(kMagenta);
    ldac_sfit_copy->SetLineColor(kMagenta);
    hist_diff_fit[ch].SetMaximum(kDrawMaxFit);
    hist_diff_fit[ch].SetMinimum(0);
    hist_diff_fit[ch].Draw();
    ldac_sfit->Draw("same");
    
    /**** differential plot ****/
    c2->cd(ch+1);
    //hist_diff[ch].Rebin(8);
    hist_diff[ch].Rebin(3);
    fit_sdiff     = new TF1("fit_sdiff","gaus",dac_min,dac_max);
    //fit_sdiff->SetLineWidth(1);
    //if(ch==8 || ch==18 || ch==44 || ch==56)
    mean_sum += mean_diff[ch]-dac_ped[ch];
    trigger_mean++;
    
    // fit the derivated polynom fit
    Double_t fit_par[3];
    //hist_diff[ch].Fit(fit_sdiff,"RQ0","",fit_min,fit_max);
    hist_diff[ch].Fit(fit_sdiff,"","",fit_min,fit_max);
    hist_diff[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
    hist_diff[ch].SetMinimum(0);
    hist_diff[ch].SetMaximum(kDrawMaxDiff);
    //hist_diff[ch].Rebin();
    hist_diff[ch].Draw();
    fit_sdiff->GetParameters(fit_par); // get the parameters for this fit
    dac_spe[ch] = fit_par[1];
    if(dac_spe[ch]>fit_min){
        dac_spe_sum += dac_spe[ch]-dac_ped[ch];
        //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
        trigger_fit++;
    }else{
        dac_spe[ch] = fit_min+10;
    }
    //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
    TLine *ldac_spe    = new TLine(dac_spe[ch],0,dac_spe[ch],
                                   (Double_t)kDrawMaxDiff);
    ldac_spe->SetLineColor(kRed);
    TLine *lmean_diff  = new TLine(mean_diff[ch],0,mean_diff[ch],
                                   (Double_t)kDrawMaxDiff);
    lmean_diff->SetLineColor(kCyan);
    TLine *lfit_min    = new TLine(fit_min,0,fit_min,kDrawMaxDiff);
    lfit_min->SetLineColor(kGray);
    
    ldac_spe->Draw();
    lmean_diff->Draw("same");
    lfit_min->Draw("same");
    
    c3->cd(1);
    if(ch==0){
        hist_scurve[ch].SetMaximum(kDrawMax);
        hist_scurve[ch].Draw();
    }else{
        /*if(ch==8 || ch==18 || ch==44 || ch==56)*/
        hist_scurve[ch].Draw("same");
    }
    if(ch==kDrawCh){
        c3->cd(2);
        hist_scurve[ch].SetMaximum(kDrawMax);
        hist_scurve[kDrawCh].Draw();
        lhalf_max->Draw("same");
        ldac_half->Draw("same");
        ldac_sfit_copy->Draw("same");
        c3->cd(3);
        hist_diff[kDrawCh].Draw();
        ldac_spe->Draw();
        lmean_diff->Draw("same");
        lfit_min->Draw("same");
        c3->cd(4);
        hist_diff_fit[kDrawCh].Draw();
        ldac_sfit->Draw("same");
    }
    /*
     if(ch<32){
     c6->cd(ch/4+1);
     hist_scurve[ch].SetMaximum(kDrawMax);
     if(ch%4==0){
     hist_scurve[ch].Draw();
     }else{
     hist_scurve[ch].Draw("same");
     }
     ldac_sfit_copy->Draw("same");
     }else{
     c7->cd(ch/4-7);
     hist_scurve[ch].SetMaximum(kDrawMax);
     if(ch%4==0){
     hist_scurve[ch].Draw();
     }else{
     hist_scurve[ch].Draw("same");
     }
     ldac_sfit_copy->Draw("same");
     }
     */
    
    c4->cd();
    if(ch==0){
        hist_scurve_ped[ch].SetMaximum(kDrawMaxPed);
        hist_scurve_ped[ch].Draw("l");
    }else{
        hist_scurve_ped[ch].Draw("lsame");
    }
    
    
    if(DEBUG)
    {
        cout << " >>>>>>> END Analysis::ProcessSingleChannel " << endl;
    }
    
}// end of process single  channels
//---------------------------------------------------------------------------------------------------


//---------------------------------------------------------------------------------------------------
// Compute some averages
//---------------------------------------------------------------------------------------------------
void Analysis::ComputeDACAverage(Double_t dac_half_spe_fix ,Int_t    kChRef, Bool_t   checkAve )//0:use ref pix, 1:check average
{
    
    
    if(DEBUG)
    {
        cout << " >>>>>>> START Analysis::ComputeDACAverage " << endl;
        
        cout << " \t - dac_half_spe_fix = " << dac_half_spe_fix << endl;
        cout << " \t - kChRef = " << kChRef << " ASIC_C " << endl;
        cout << " \t - checkAve = " << checkAve << " 0:use ref pix, 1:check average " << endl;
    }
    
    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
    
    
    mean_diff_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
    
    if(!checkAve){
        dac_ped_ave   = dac_ped[kChRef]; //absolute value
        dac_half_ave  = dac_half[kChRef];
        dac_spe_ave   = dac_spe[kChRef];
        dac_sfit_ave  = dac_sfit[kChRef];
        //Float_t dac_sfit_ave  = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
        mean_diff_ave = mean_diff[kChRef];
    }else{
        dac_ped_ave  = dac_ped_sum/(Float_t)kNch; //absolute value
        dac_half_ave = dac_half_sum/(Float_t)trigger+dac_ped_ave;
        dac_spe_ave  = dac_spe_sum/(Float_t)trigger_fit+dac_ped_ave;
        mean_diff_ave = mean_sum/(Float_t)trigger_mean+dac_ped_ave;
        dac_sfit_ave = dac_sfit_sum/(Float_t)trigger_sfit+dac_ped_ave;
    }
    
    
    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
    
    
    //trigger_mean=4;
    if(dac_half_spe_fix!=-1) {
        dac_half_spe = dac_half_spe_fix;
    } else {
        //Double_t dac_half_spe = dac_half_ave;
        dac_half_spe = dac_spe_ave;
    }
    
    TLine *ldac_ped_ave   = new TLine(dac_ped_ave,0,dac_ped_ave,
                                      (Double_t)kDrawMax);
    TLine *ldac_ped_ave_copy = new TLine(dac_ped_ave,0,dac_ped_ave,
                                         (Double_t)kDrawMaxDiff);
    TLine *ldac_half_ave  = new TLine(dac_half_ave,0,dac_half_ave,
                                      (Double_t)kDrawMax);
    TLine *ldac_half_spe  = new TLine(dac_half_spe,0,dac_half_spe,
                                      (Double_t)kDrawMax);
    TLine *lmean_diff_ave = new TLine(mean_diff_ave,0,mean_diff_ave,
                                      (Double_t)kDrawMax);
    TLine *ldac_sfit_ave  = new TLine(dac_sfit_ave,0,dac_sfit_ave,
                                      (Double_t)kDrawMax);
    ldac_ped_ave->SetLineColor(kOrange);
    ldac_half_ave->SetLineColor(kRed);
    lmean_diff_ave->SetLineColor(kCyan);
    //ldac_half_spe->SetLineColor(kGreen);
    ldac_sfit_ave->SetLineColor(kMagenta);
    
    c3->cd(1);
    ldac_ped_ave->Draw("same");
    ldac_half_ave->Draw("same");
    // lhalf_max->Draw("same");  SDC (24/09/2013) :: Not declared
    ldac_half_spe->Draw("same");
    lmean_diff_ave->Draw("same");
    ldac_sfit_ave->Draw("same");
    
    c3->cd(3);
    ldac_ped_ave_copy->SetLineColor(kOrange);
    ldac_ped_ave_copy->Draw("same");
    
    if(DEBUG)
    {
        cout << " >>>>>>> END Analysis::ComputeDACAverage " << endl;
    }
    
    
}

//---------------------------------------------------------------------------------------------------

//---------------------------------------------------------------------------------------------------
void Analysis::OutputGain(const TString fname_out)
{
    
    if(DEBUG)
    {
        cout << " >>>>>>> START Analysis::OutputGain " << endl;
        cout << " \t - fname_out = " << fname_out << endl;
    }
    
    ofstream out(fname_out.Data());
    
    Int_t    gain_asic_diff, gain_asic_dacHalf, gain_asic_diff_mean, gain_asic_sfit;
    out << "title  ch  g_asic_dacHalf  g_asic_diff  g_asic_diff_mean  g_asic_sfit  dac_half[ch]  dac_ped[ch]  dac_spe[ch]  dac_half_ave dac_half_spe dac_ped_ave  dac_sfit_ave"
    << endl;
    
    
    // loop 3 on all channels to make final computations on the gain 
    for(ch=0;ch<kNch;ch++){
        /*
         if(ch==8 || ch==44 || ch==56 || ch==18){
         gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
         (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
         }else{
         gain_asic_diff_mean=0;
         }
         */
        if(dac_half[ch]<0){
            gain_asic_dacHalf   = 0;
            gain_asic_diff      = 0;
            gain_asic_diff_mean = 0;
            gain_asic_sfit      = 0;
        }else{
            gain_asic_dacHalf  =(Int_t)((dac_half_ave-dac_ped_ave)/
                                        (dac_half[ch]-dac_ped[ch])*gain_org[ch]);
            gain_asic_diff     =(Int_t)((dac_half_spe-dac_ped_ave)/
                                        (dac_spe[ch]-dac_ped[ch])*gain_org[ch]);
            gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
                                        (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
            gain_asic_sfit     =(Int_t)((dac_sfit_ave-dac_ped_ave)/
                                        (dac_sfit[ch]-dac_ped[ch])*gain_org[ch]);
            //cerr << ch << " " << gain_org[ch] << endl;
        }
        
        
        
        out << "OUT_TABLE    "
        << ch << "    "
        << gain_asic_dacHalf   << "    "
        << gain_asic_diff      << "    "
        << gain_asic_diff_mean << "    "
        << gain_asic_sfit      << "    "
        << dac_half[ch]        << "    "
        << dac_ped[ch]         << "    "
        << (Int_t)dac_spe[ch]  << "    "
        << dac_half_ave        << "    "
        << dac_half_spe        << "    "
        << dac_ped_ave         << "    "
        << dac_sfit_ave        << "    "
        << sfit_chiS[ch]       << "    "
        << ch                  << "    "
        << flag_fit[ch]       
        << endl;
        
        if(DEBUG>2)
        {
            std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
        }
    }
    out.close();
    
    if(DEBUG)
    {
        cout << " >>>>>>> END Analysis::OutputGain " << endl;
        
    }
    
    
}
//--------------------------------------------------------------------------------------------------------







//--------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------
void Analysis::Help()
{
    
    cout << "This is the Help of the PlotScurveAll-diffFit" <<endl;
    cout << " Version 2.0 from September 27 2013 " <<endl;
    cout << " Main author Hiroko Miyamoto " << endl;
    cout << " Adapted for modularity and compilable program by Sylvie Dagoret " << endl;
    cout << " ==> Please from shell execute directly the compiled program :  " <<endl;
    cout << " ==> ../../src/PlotScurvesAllMain  " << endl;
    
}
//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------
void Analysis::Menu()
{
    cout << "This is the Detailed menu information for running the PlotScurveAll-diffFit" <<endl;
    
    cout << " void PlotScurve( " << endl;
    cout << " \t " <<  " Int_t    kDrawCh=31, " << endl;
    cout << " \t " <<  " TString  fname=kFname, " << endl;
    cout << " \t " <<  " TString  fname_ped=kFnamePed, //ASIC_D ped_file " << endl;
    cout << " \t " <<  " Double_t dac_half_spe_fix = -1," << endl;
    cout << " \t " <<  " TString  fname_gain=\"gain-org.txt\"," << endl;
    cout << " \t " <<  " Int_t    fit_min=kFitMin," << endl;
    cout << " \t " <<  " Int_t    fit_max=kFitMax," << endl;
    cout << " \t " <<  " Int_t    fit_min_scurve=kFitMinS," << endl;
    cout << " \t " <<  " Int_t    fit_max_scurve=kFitMaxS," << endl;
    cout << " \t " <<  " Bool_t   kBoard=1, //0:ASIC_C, 1:ASIC_D " << endl;
    cout << " \t " <<  " Int_t    kChRef=kRefCh,//ASIC_C " << endl;
    cout << " \t " <<  " Bool_t   checkAve=0 //0:use ref pix, 1:check average " << endl;
    cout << " ) " << endl;		
    cout << "Set the parameters on the top of the program depending on your analysis:  " << endl;
    cout << "General policy: parameters starting with a \"k\" at the beginning is fixed gloval parameter and usually set in the header file, never change through the program. " << endl;
    cout << " : " << endl;
    cout << " - kDrawCh: channel just you want to check on display. doens't effect the analysis. " << endl;
    
    cout << " - fname, fname_ped (kFname, kFnamePed): \"YOUR DATA FILE\" without .data"  << endl;
    cout << " See examples in the header file. " << endl;
    cout << " - dac_half_spe_fix: set DAC value only when you want to adjust the gain to particular pulse height (DAC value) instead of either average of the pulse height of reference pixel. Usually just leave it -1. " << endl;
    cout << " - fit_min, fit_max (kFitMin, kFitMax): leave them as default, it is not currently used for the analysis."  << endl;
    cout << " - fit_min_scurve, fit_max_scurve (kFitMinS, kFitMaxS): define the range of s-curve fitting. search proper values by yourself with running root program a few times."  << endl;
    cout << " - kBoard: don't forget to change 0(C) or 1(D) depending on which ASIC you're analysing "  << endl;
    cout << " - kChRef (kRefCh): a pixel whose gain remains 64 thorough the gain adjustment. " << endl;
    cout << " - checkAve: for the first run, set this parameter 1, and find a pixel which remains 64. In the case of parameter 1, all the gains are adjusted to the average pulse height. " << endl;
    cout <<  " Once you choose the reference pixel, set the parameter 0, as well as set kRefCh the pixel number. (better to keep in the header file for the record of each PMTs.) " << endl;
    
    cout << " Actual configuration " << endl;
    cout << " \t " <<  "TString kFname = " << kFname << endl;
    cout << " \t " <<  "TString  kFnamePed = " << kFnamePed  << endl;
    cout << " \t " <<  "Int_t kFitMin = " << kFitMin << endl; 
    cout << " \t " <<  "Int_t kFitMax = " << kFitMax << endl;
    
}

//-------------------------------------------------------------------------------------------------
// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "(make;) "
// End: