source: JEM-EUSO/equalization_gain/trunk/src/PlotScurvesAll-diffFit.cc @ 201

#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) 
{
 
  cout << " analysis initialisation " << 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;


}
//--------------------------------------------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
// Analysis destructor
//--------------------------------------------------------------------------------------------------------
Analysis::~Analysis()
{
 
  if(input!=0) delete input;
  if(input_ped!=0) delete input_ped;


}
//--------------------------------------------------------------------------------------------------------


//-------------------------------------------------------------------------------------------------------- 
// 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
                )



{

  // 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);



}
//--------------------------------------------------------------------------------------------------------


//--------------------------------------------------------------------------------------------------------
// Initialize the program by reseting counters or variables and create objects
//--------------------------------------------------------------------------------------------------------
void Analysis::Init()
{

  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;
 }


  
}
//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
// Initialize the program by reseting counters or variables and create objects
//--------------------------------------------------------------------------------------------------------
void Analysis::OpenGainFile(const TString fname_gain)
{
  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();

}
//--------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------
// Fill Histogram name
//--------------------------------------------------------------------------------------------------------
void Analysis::FillHistoName()
{
  
// 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;
  }

}
//--------------------------------------------------------------------------------------------------------


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

void Analysis::ReadGainFile(const TString fname_gain,Bool_t kBoard)
{
  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();
 

}
//--------------------------------------------------------------------------------------------------------



//--------------------------------------------------------------------------------------------------------
// 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)
                                    )
{
    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");
    }
    
  }// end of process single  channels

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


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

//---------------------------------------------------------------------------------------------------
void Analysis::ComputeDACAverage(Double_t dac_half_spe_fix ,Int_t    kChRef, Bool_t   checkAve )//0:use ref pix, 1:check average
{

 // 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");

 

}

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

//---------------------------------------------------------------------------------------------------
void Analysis::OutputGain(const TString fname_out)
{

 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;
      std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
  }
  out.close();

}