Changeset 204 in JEM-EUSO

Timestamp:

Sep 27, 2013, 1:23:30 PM (11 years ago)

Author:

dagoret

Message:

working class analysis

Location:

equalization_gain/branches/analclass-Sylvie/src

Files:

• PlotScurvesAll-diffFit.cc (modified) (21 diffs)
• PlotScurvesAll-diffFit.h (modified) (6 diffs)
• PlotScurvesAll-diffFitConst.h (modified) (3 diffs)
• PlotScurvesAllMain.cc (modified) (1 diff)

equalization_gain/branches/analclass-Sylvie/src/PlotScurvesAll-diffFit.cc

@@ -17 +17 @@
 #include "PlotScurvesAll-diffFitConst.h"
 #include "PlotScurvesAll-diffFit.h"
+
+//#if defined(__CINT__) && !defined(__MAKECINT__) 
+//#include "PlotScurveAll-diffFit.h"
+//#else 
+//#include "PlotScurvesAll-diffFit.h"
+//#endif
 
 
@@ -214 +220 @@
   }
   
-  ofstream out(kFnameOut.Data()); 
+ 
 
   // first loop on channels to fill histonames
@@ -262 +268 @@
   in_gain.clear();
  
-
+  //----
  
-  // Open the S-Curve files for each channel and test if the file can be open
+  // 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;
@@ -349 +355 @@
     
     // 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,
@@ -358 +366 @@
     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
@@ -686 +698 @@
   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" 
@@ -737 +751 @@
          << flag_fit[ch]       
          << endl;
+      std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
   }
 
@@ -759 +774 @@
 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 ********************/
@@ -775 +790 @@
   maxBinSum         = 0;
 
-
+ if(DEBUG)
+    cout << " >>>>>>> END Analysis Constructor " << endl;
 }
 //--------------------------------------------------------------------------------------------------------
@@ -786 +802 @@
 {
  
+ 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;
 
 }
@@ -797 +818 @@
 // Main public function of Analysis
 //--------------------------------------------------------------------------------------------------------
+/*
 void Analysis::PlotScurve(
                 Int_t    kDrawCh=31,
@@ -809 +831 @@
                 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*/
+                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();
@@ -822 +878 @@
 
   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;
+
 }
 //--------------------------------------------------------------------------------------------------------
@@ -832 +913 @@
 {
 
-  if(input!=0) delete input;
-  if(input_ped!=0) delete input_ped;
+if(DEBUG)
+    cout << " >>>>>>> START Analysis::Init " << endl;
+
+//if(input!=0) delete input;
+//if(input_ped!=0) delete input_ped;
 
   // recreate the name 
@@ -881 +965 @@
  }
 
+if(DEBUG)
+    cout << " >>>>>>> END Analysis::Init " << endl;
 
   
@@ -891 +977 @@
 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());
@@ -900 +991 @@
   in_gain.clear();
 
+if(DEBUG)
+    cout << " >>>>>>> END Analysis::OpenGainFile " << endl;
+
+
 }
 //--------------------------------------------------------------------------------------------------------
@@ -908 +1003 @@
 void Analysis::FillHistoName()
 {
+
+if(DEBUG)
+    cout << " >>>>>>> START Analysis::FillHistoName " << endl;
+
+
   
 // first loop on channels to fill histonames
@@ -925 +1025 @@
   }
 
+if(DEBUG)
+    cout << " >>>>>>> END Analysis::FillHistoName " << endl;
+
+
+
 }
 //--------------------------------------------------------------------------------------------------------
@@ -931 +1036 @@
 //--------------------------------------------------------------------------------------------------------
 // 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());
@@ -963 +1081 @@
   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;
+      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:

equalization_gain/branches/analclass-Sylvie/src/PlotScurvesAll-diffFit.h

@@ -47 +47 @@
   // Main analysis function to be called from outside
 
+  void Menu();
+  void Help();
+  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=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 PlotScurve(
                 Int_t    kDrawCh,
@@ -62 +79 @@
                 );
 
-
+*/
 
 
@@ -71 +88 @@
  void FillHistoName();
  void ReadGainFile(const TString,Bool_t);
+ void 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
+         );
+ void ComputeDACAverage(Double_t dac_half_spe_fix ,Int_t    kChRef,Bool_t checkAve); //0:use ref pix, 1:check average)
+ void OutputGain(const TString fname_out);
 
 // common variables to the analysis
      TString *input;
      TString *input_ped;
-     Int_t    ch;
+     Int_t    ch; // index on current channel under analysis
      Float_t  dac_half_sum;  
      Float_t  dac_ped_sum;   
@@ -90 +122 @@
      Double_t dac_sfit[kNch];
      Double_t mean_diff[kNch];
-     Double_t gain_org[kNch];
+     Double_t gain_org[kNch]; // array containing the initial gain values 
      Double_t sfit_chiS[kNch];
      Double_t flag_fit[kNch];
@@ -115 +147 @@
      TF1     *fit_sdiff;  // function fitted on differentiated histograms of fitted s-curves
 
-
+     Int_t trigger; 
+     Int_t trigger_fit; 
+     Int_t trigger_mean;
+     Int_t trigger_sfit;
      
+     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 ;//dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
+  
+     Double_t dac_half_spe;
 
 };
@@ -122 +164 @@
 #endif
 
+//-------------------------------------------------------------------------------------------------
+// Configure (x)emacs for this file ...
+// Local Variables:
+// mode: c++
+// compile-command: "(make;) "
+// End:

equalization_gain/branches/analclass-Sylvie/src/PlotScurvesAll-diffFitConst.h

@@ -35 +35 @@
 */
 
+
+const int DEBUG=1;
 
 /*** signal file ***/
@@ -73 +75 @@
 
 // Output file
-const TString kFnameOut="output.dat";
+const TString kFnameOut="gain-output.dat";
 
 //Int_t kFitMin = 160; //T1 ASIC_D it0, default
@@ -140 +142 @@
 
 #endif
+//-------------------------------------------------------------------------------------------------
+// Configure (x)emacs for this file ...
+// Local Variables:
+// mode: c++
+// compile-command: "(make;) "
+// End:

equalization_gain/branches/analclass-Sylvie/src/PlotScurvesAllMain.cc

@@ -20 +20 @@
    TRint *theApp=new TRint("ROOT-PlotScurvesAll",&argc,argv);
     
-   Help();
    
-   PlotScurve();
+   // SDC (27/09/2013) :: Old fashion of calling the function directly
+   //PlotScurve();
 
-   //Analysis* pa = new Analysis;
-   //pa->Print();
+   // SDC (27/09/2013) :: class analysis working now 
+   Analysis* ana = new Analysis;
+
+   ana->Help();
+   ana->Menu();
+   ana->PlotScurve();
 
     theApp->Run();
 }
+
+//-------------------------------------------------------------------------------------------------
+// Configure (x)emacs for this file ...
+// Local Variables:
+// mode: c++
+// compile-command: "(make;) "
+// End: