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

Last change on this file since 208 was 208, checked in by moretto, 11 years ago

addition of default gain tables for ASIC C & D needed to create the default gain-org.txt gain table

File size: 65.8 KB
Line 
1#include <fstream>
2#include <iostream>
3#include "TROOT.h"
4#include "TString.h"
5#include "TH1D.h"
6#include "TH2D.h"
7#include "TCanvas.h"
8#include "TF1.h"
9#include "TPad.h"
10#include "TGraphErrors.h"
11#include "TGraph.h"
12#include "TLine.h"
13#include "TROOT.h"
14#include "TString.h"
15#include <TStyle.h>
16
17#include "PlotScurvesAll-diffFitConst.h"
18#include "PlotScurvesAll-diffFit.h"
19
20//#if defined(__CINT__) && !defined(__MAKECINT__)
21//#include "PlotScurveAll-diffFit.h"
22//#else
23//#include "PlotScurvesAll-diffFit.h"
24//#endif
25
26
27
28// for format of EC_ASIC Scurves_all_ASIC data
29// ASIC A,B,C,.... : i=0,1,2,... ch 0,1,...,63 : j=0,1,...,63
30// % cat scurves_PC_inj_run5.lvm |awk '{i=2;j=37;print($1,$(i*64+j+1+1))}'
31//   > scurves_PC_inj_run5-C-ch37.data
32// => % cd $datadir
33// => % sh $srcdir/makeDataFile.sh with assigining proper fname
34// PlotScurvesAll-diffFit-v6.6.C: dac_sfit[ch]=-1; in the case of fitting failure.
35
36using namespace std;
37Int_t          kXmin       = 0;
38
39//--------------------------------------------------------------------------------------------------------
40void Help()
41{
42   
43    cout << "This is the Help of the PlotScurveAll-diffFit" <<endl;
44    cout << " Version 1.0 frol September 24 2013 " <<endl;
45    cout << " Main author Hiroko Miyamoto " << endl;
46    cout << " ==> Please in ROOT interpreter, issue the command .L PlotScurvesAll-diffFit.cc " <<endl;
47    cout << " ==> Then play with PlotScurve() function " << endl;
48    cout << " ==> If you need more information, please type Menu()" << endl;
49   
50}
51//--------------------------------------------------------------------------------------------------------
52
53//--------------------------------------------------------------------------------------------------------
54void Menu()
55{
56    cout << "This is the Detailed menu information for running the PlotScurveAll-diffFit" <<endl;
57   
58    cout << " void PlotScurve( " << endl;
59    cout << " \t " <<  " Int_t    kDrawCh=31, " << endl;
60    cout << " \t " <<  " TString  fname=kFname, " << endl;
61    cout << " \t " <<  " TString  fname_ped=kFnamePed, //ASIC_D ped_file " << endl;
62    cout << " \t " <<  " Double_t dac_half_spe_fix = -1," << endl;
63    cout << " \t " <<  " TString  fname_gain=\"gain-org.txt\"," << endl;
64    cout << " \t " <<  " Int_t    fit_min=kFitMin," << endl;
65    cout << " \t " <<  " Int_t    fit_max=kFitMax," << endl;
66    cout << " \t " <<  " Int_t    fit_min_scurve=kFitMinS," << endl;
67    cout << " \t " <<  " Int_t    fit_max_scurve=kFitMaxS," << endl;
68    cout << " \t " <<  " Bool_t   kBoard=1, //0:ASIC_C, 1:ASIC_D " << endl;
69    cout << " \t " <<  " Int_t    kChRef=kRefCh,//ASIC_C " << endl;
70    cout << " \t " <<  " Bool_t   checkAve=0 //0:use ref pix, 1:check average " << endl;
71    cout << " ) " << endl;
72    cout << "Set the parameters on the top of the program depending on your analysis:  " << endl;
73    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;
74    cout << " : " << endl;
75    cout << " - kDrawCh: channel just you want to check on display. doens't effect the analysis. " << endl;
76   
77    cout << " - fname, fname_ped (kFname, kFnamePed): \"YOUR DATA FILE\" without .data"  << endl;
78    cout << " See examples in the header file. " << endl;
79    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;
80    cout << " - fit_min, fit_max (kFitMin, kFitMax): leave them as default, it is not currently used for the analysis."  << endl;
81    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;
82    cout << " - kBoard: don't forget to change 0(C) or 1(D) depending on which ASIC you're analysing "  << endl;
83    cout << " - kChRef (kRefCh): a pixel whose gain remains 64 thorough the gain adjustment. " << endl;
84    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;
85    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;
86   
87    cout << " Actual configuration " << endl;
88    cout << " \t " <<  "TString kFname = " << kFname << endl;
89    cout << " \t " <<  "TString  kFnamePed = " << kFnamePed  << endl;
90    cout << " \t " <<  "Int_t kFitMin = " << kFitMin << endl;
91    cout << " \t " <<  "Int_t kFitMax = " << kFitMax << endl;
92   
93}
94//--------------------------------------------------------------------------------------------------------
95
96//-------------------------------------------------------------------
97void InitRoot()
98{
99   
100    gStyle->SetTitleFillColor(kWhite);
101    gStyle->SetTitleColor(kBlue,"xyz");
102    gStyle->SetTitleSize(0.03,"xyz");
103    gStyle->SetLabelSize(0.03,"xyz");
104    gStyle->SetLabelColor(kBlue,"xyz");
105    gStyle->SetAxisColor(kMagenta);
106    gStyle->SetPalette(1,0);
107    gStyle->SetOptFit(101);
108    gStyle->SetOptStat(0);
109   
110    //  gStyle->SetOptStat("emr");
111    gStyle->SetOptStat("");
112}
113//-------------------------------------------------------------------
114
115
116//--------------------------------------------------------------------------------------------------------
117/*
118 
119 The declaration of this function is now put inside the h file with the default values to be used
120 
121 void PlotScurve(
122 Int_t    kDrawCh=31,
123 TString  fname=kFname,
124 TString  fname_ped=kFnamePed, //ASIC_D ped_file
125 Double_t dac_half_spe_fix = -1,
126 TString  fname_gain="gain-org.txt",
127 Int_t    fit_min=kFitMin,
128 Int_t    fit_max=kFitMax,
129 Int_t    fit_min_scurve=kFitMinS,
130 Int_t    fit_max_scurve=kFitMaxS,
131 Bool_t   kBoard=1, //0:ASIC_C, 1:ASIC_D
132 Int_t    kChRef=kRefCh,//ASIC_C
133 Bool_t   checkAve=0 //0:use ref pix, 1:check average
134 )
135 */
136
137// Definition of the function  PlotScurve
138void PlotScurve(
139                Int_t    kDrawCh,
140                TString  fname,
141                TString  fname_ped, //ASIC_D ped_file
142                Double_t dac_half_spe_fix,
143                TString  fname_gain,
144                Int_t    fit_min,
145                Int_t    fit_max,
146                Int_t    fit_min_scurve,
147                Int_t    fit_max_scurve,
148                Bool_t   kBoard, //0:ASIC_C, 1:ASIC_D
149                Int_t    kChRef,//ASIC_C
150                Bool_t   checkAve //0:use ref pix, 1:check average
151                )
152
153
154{
155    //NOTE
156    //Double_t dac_half_spe_fix = 177; //for EC_ASICb N101(asic N158).
157    //check DAC value corresponds to voltage to 160fC input
158    //(depends on DAC linearity)
159    //Double_t dac_half_spe_fix = -1 //set gain to DAC average (default)
160   
161    /******************* Readout data file *********************/
162   
163    /*********************** SIGNAL RUN ********************/
164    TString *input             = new TString[kNch];
165    TString *input_ped         = new TString[kNch];
166    Int_t    ch;
167    Float_t  dac_half_sum      = 0;
168    Float_t  dac_ped_sum       = 0;
169    Float_t  dac_spe_sum       = 0;
170    Float_t  dac_sfit_sum      = 0;
171    Float_t  diff_sum          = 0;
172    Float_t  mean_sum          = 0;
173    Float_t  maxBinSum         = 0;
174    Double_t dac_half[kNch];
175    Double_t dac_ped[kNch];
176    Double_t dac_spe[kNch];
177    Double_t dac_sfit[kNch];
178    Double_t mean_diff[kNch];
179    Double_t gain_org[kNch];
180    Double_t sfit_chiS[kNch];
181    Double_t flag_fit[kNch];
182   
183   
184    InitRoot();
185   
186    TCanvas *c1                = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
187    c1->Divide(11,6);
188    TCanvas *c2                = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
189    c2->Divide(11,6);
190    TCanvas *c5                = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
191    c5->Divide(11,6);
192    TCanvas *c3                = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
193    c3->Divide(4,1);
194    /*TCanvas *c6                = new TCanvas("c6","c6",450,150,1000,500);
195     c6->Divide(4,2);
196     TCanvas *c7                = new TCanvas("c7","c7",470,200,1000,500);
197     c7->Divide(4,2);*/
198    TCanvas *c4                = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
199   
200    TString  name_hist[kNch];
201    TString  name_hist_ped[kNch];
202    TString  name_hist_diff[kNch];
203    TString  name_hist_diff_fit[kNch];
204    TString  name_fit_diff[kNch];
205    TH1D    *hist_scurve       = new TH1D[kNch]; // smoothed S-curves
206    TH1D    *hist_scurve_ped   = new TH1D[kNch]; // S-curves for pedestal runs
207    TH1D    *hist_diff         = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
208    TH1D    *hist_diff_fit     = new TH1D[kNch]; // differentiated histograms of fitted s-curves
209    TF1     *fit_sdiff         = new TF1[kNch];  // function fitted on differentiated histograms of fitted s-curves
210   
211   
212    // Open the ASIC Gain file
213    //-------------------------
214   
215    TString input_gain(fname_gain);
216    ifstream in_gain(input_gain.Data());
217    if(!in_gain.good()){
218        cerr << "File " << input_gain << " open failed!" << endl;
219        return;
220    }
221   
222   
223   
224    // first loop on channels to fill histonames
225    //----------------------
226    for(ch=0;ch<kNch;ch++){
227        dac_half[ch]=0;
228        dac_ped[ch]=0;
229        dac_spe[ch]=0;
230        dac_sfit[ch]=0;
231        gain_org[ch]=0;
232        mean_diff[ch]=0;
233        sfit_chiS[ch]=0;
234        flag_fit[ch]=0;
235        name_hist[ch]           ="hist";
236        name_hist[ch]          +=ch;
237        name_hist_ped[ch]       ="hist_ped";
238        name_hist_ped[ch]      +=ch;
239        name_hist_diff[ch]      ="hist_diff";
240        name_hist_diff[ch]     +=ch;
241        name_hist_diff_fit[ch]  ="hist_diff_fit";
242        name_hist_diff_fit[ch] +=ch;
243        name_fit_diff[ch]       ="fit_diff";
244        name_fit_diff[ch]      +=ch;
245    }
246   
247   
248    // read gain
249    Double_t da;
250    Int_t count=0;
251    ch=0;
252    while(in_gain >> da){
253        if(!kBoard && (count>=128 && count<192)){
254            gain_org[ch]=da;
255            cerr << "!kBoard, count, gain_org[" << ch << "] "
256            << !kBoard << " " << count << " " << gain_org[ch] << endl;
257            ch++;
258        }
259        if(kBoard && (count>=192 && count<256)){
260            gain_org[ch]=da;
261            cerr << "kBoard, count, gain_org[" << ch << "] "
262            << kBoard << " " << count << " " << gain_org[ch] << endl;
263            ch++;
264        }
265        count++;
266    }// end if first loop on all channels
267    in_gain.close();
268    in_gain.clear();
269   
270    //----
271   
272    // Open the pedestal file and the S-Curve files for each channel and test if the file can be open
273   
274    Int_t trigger=0, trigger_fit=0, trigger_mean=0, trigger_sfit=0;
275    // Big loop on all the channels : second loop
276    //--------------------------------------------
277    // build the input data file name for each channel fname-chXX.dat
278    for(ch=0;ch<kNch;ch++){
279        input[ch] += fname;
280        input[ch] += "-ch";
281        input[ch] += ch;
282        input[ch] += ".dat";
283        ifstream in(input[ch].Data());
284        if(!in.good()){
285            cerr << "File " << input[ch] << " open failed!" << endl;
286            return;
287        }
288       
289        // build the input data file name for each channel fname_ped-chXX.dat
290        input_ped[ch] += fname_ped;
291        input_ped[ch] += "-ch";
292        input_ped[ch] += ch;
293        input_ped[ch] += ".dat";
294        ifstream in_ped(input_ped[ch].Data());
295        if(!in_ped.good()){
296            cerr << "File " << input_ped[ch] << " open failed!" << endl;
297            return;
298        }
299       
300       
301       
302        //10DAC~10mV
303        // determine the DAC range in dac_min, dac_max, step_dac from file input_ped
304        // and count the number of events
305        Double_t da, dat, dac_min, dac_max, step_dac;
306        count = 0;
307        while (in >> da >> dat ){
308            if(count==0) dac_min  = da;
309            if(count==1) step_dac = da - dac_min;
310            dac_max=da;
311            count++;
312        }
313        in.close();
314        in.clear();
315       
316        /*
317         if(ch==43 || ch==47){//for T1 ASIC_C
318         cerr << "CHECK " << ch << " " << fit_min_scurve << endl;
319         Int_t fit_min_scurve=140;
320         } else {
321         Int_t fit_min_scurve=kFitMinS;
322         }
323         */
324       
325        // correct the scale for DAC range
326        if(kXmin<dac_min)kXmin=dac_min;
327        const Int_t kNdata    = count;
328        count = (fit_max_scurve-fit_min_scurve)/step_dac;
329        const Int_t kNdataFit = count;
330        if(ch==0){
331            cerr << "kNdata, dac_min, dac_max, step_dac: "
332            << kNdata << " " << dac_min << " " << dac_max << " "
333            << step_dac << endl;
334        }
335       
336        Double_t dac_min_ped, dac_max_ped, step_dac_ped;
337        count = 0;
338        while (in_ped >> da >> dat ){
339            if(count==0) dac_min_ped  = da;
340            if(count==1) step_dac_ped = da - dac_min_ped;
341            dac_max_ped=da;
342            count++;
343        }
344       
345        in_ped.close();
346        in_ped.clear();
347       
348       
349        const Int_t kNdataPed = count;
350        if(ch==0){
351            cerr << "kNdataPed, dac_min_ped, dac_max_ped, step_dac_ped: "
352            << kNdataPed << " " << dac_min_ped << " " << dac_max_ped
353            << " " << step_dac_ped << endl;
354        }
355       
356        // book histograms
357        // SDC: 27/09/2013 hist_scurve is a pointer on an array of THD, thus the syntax is correct (no new operator)
358       
359       
360        hist_scurve[ch]      = TH1D(name_hist[ch],name_hist[ch],kNdata-1,
361                                    dac_min,dac_max);
362        hist_scurve_ped[ch]  = TH1D(name_hist_ped[ch],name_hist_ped[ch],
363                                    kNdataPed-1,dac_min_ped,dac_max_ped);
364        hist_diff[ch]        = TH1D(name_hist_diff[ch],name_hist_diff[ch],
365                                    kNdata-2,dac_min,dac_max);
366        hist_diff_fit[ch]    = TH1D(name_hist_diff_fit[ch],name_hist_diff_fit[ch],
367                                    kNdataFit,fit_min_scurve,fit_max_scurve);
368       
369       
370       
371       
372       
373        // read the data for that channel for normal data
374        in.open(input[ch].Data());
375        Double_t       data[kNdata], dac[kNdata], dac_fit[kNdataFit], data_fit[kNdataFit];
376        while (in >> da >> dat ){
377            hist_scurve[ch].Fill(da,dat); // fill the histograms of the s-curve
378        }
379        in.close();
380        in.clear();
381       
382        // read the data for pedestal data
383        in_ped.open(input_ped[ch].Data());
384        while (in_ped >> da >> dat ){
385            hist_scurve_ped[ch].Fill(da,dat); // fill the histogram for pedestal data
386        }
387        in_ped.close();
388        in_ped.clear();
389       
390        //dac_ped[ch]  = hist_scurve_ped[ch].GetMaximumBin()+dac_min_ped; //absolute value
391        dac_ped[ch]  = hist_scurve_ped[ch].GetBinCenter(hist_scurve_ped[ch].GetMaximumBin()); //absolute value
392        dac_ped_sum += dac_ped[ch]; //absolute value
393        //cerr << "dac_ped[" << ch << "] dac_ped_sum " << dac_ped[ch]
394        //<< " " << dac_ped_sum << endl;
395       
396       
397        hist_scurve[ch].Smooth(3); // Smooth the S curve
398        Int_t binx=0;
399        //hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
400        //                                  dac_max-dac_min+1,kHalfMax/2.);
401       
402        // ??????
403        hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
404                                          dac_max-dac_min+1,2000);
405        dac_half[ch] = hist_scurve[ch].GetBinCenter(binx); //absolute value
406        if(dac_half[ch]<kXmin)dac_half[ch]=-1;
407       
408        // put the histograms of the smoothed S-Curves into two arrays (dac, data)
409        //-------------------------------------------------------------------------
410        Int_t i=0;
411        for(i=0;i<kNdata;i++){
412            dac[i]  = hist_scurve[ch].GetBinCenter(i+1);
413            data[i] = hist_scurve[ch].GetBinContent(i+1);
414        }
415       
416        /*** compute derive diff on S-Curve ***/
417        Double_t diff[kNdata-1], dac_diff[kNdata-1];
418        Int_t trigger_diff=0;
419        for(i=0;i<kNdata-1;i++){
420            diff[i]     = (data[i]-data[i+1])/step_dac;
421            dac_diff[i] = dac[i]+step_dac/2;
422            hist_diff[ch].Fill(dac_diff[i],diff[i]); // fills differentiated histograms of smoothed s-curves which is the single photoelectron spectrum
423            if(dac_diff[i]>=fit_min){
424                diff_sum += diff[i]*dac_diff[i];
425                //cerr << "diff[" << i << "], dac_diff[" << i << "] "
426                //<< diff[i] << " " << dac_diff[i] << endl;
427                trigger_diff+=diff[i];
428            }
429        }
430        //mean_diff[ch] = diff_sum/trigger_diff;
431        hist_diff[ch].GetXaxis()->SetRangeUser(fit_min,fit_max);
432        mean_diff[ch] = hist_diff[ch].GetMean();
433        //cerr << "trigger_diff, mean_diff[" << ch << "] "
434        //<< trigger_diff << " " << mean_diff[ch] << endl;
435       
436       
437        /******************* GRAPHICS *****************/
438       
439       
440        // Now start the analysis which consist in finding the peak in the differentiated spectrum
441        c1->cd(ch+1);
442        /*** fit scurves  ***/
443        cerr << "ch " << ch;
444       
445        //    SDC 24/09/2013 :: Need to declare outside this variable in a compiled program
446        TF1 *fit_pol ;
447       
448       
449        // Fit a polynomial curve in a given range on the s-curve directly. This is a kind of DEBUG procedure.
450       
451        //if(ch==43 || ch==47){ // SDC: 25/09/2013 Why a special channel : the fit pol4 is performed for another DAC range
452        /*   if(ch==170){
453         fit_pol      = new TF1("fit_pol","pol4",150,fit_max_scurve);
454         hist_scurve[ch].Fit(fit_pol,"","",150,fit_max_scurve);
455         }else{ */
456        fit_pol      = new TF1("fit_pol","pol4",fit_min_scurve,fit_max_scurve);
457        hist_scurve[ch].Fit(fit_pol,"","",fit_min_scurve,fit_max_scurve);
458        // }
459       
460       
461       
462        Double_t fit_par_pol4[5]; // container for fit parameter results
463        fit_pol->GetParameters(fit_par_pol4);
464        sfit_chiS[ch] =  fit_pol->GetChisquare(); // chi squared
465        TF1 *func               = hist_scurve[ch].GetFunction("fit_pol");
466        Double_t data_diff_fit;
467        for(i=0;i<kNdataFit;i++){
468            //if(ch==43 || ch==47){
469            if(ch==170){
470                dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+150-dac_min;
471                data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+150-dac_min); // get the value of the fitted function
472            }else{
473                dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min;
474                data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min); // get the value of the fitted function
475                //cerr << "dac_fit[" << i << "], data_fit[" << i << "] "
476                //<< dac_fit[i] << " " << data_fit[i] << endl;
477            }
478        }
479       
480        // derive the fitted function on the S-curve to obtain the SPE
481        Double_t dac_diff_fit;
482        for(i=0;i<kNdataFit-1;i++){
483            data_diff_fit = (data_fit[i]-data_fit[i+1])/step_dac;
484            hist_diff_fit[ch].Fill(dac_fit[i],data_diff_fit); // Fill differentiated fitted S-curve
485        }
486       
487        TLine *lhalf_max   = new TLine(dac[0],kHalfMax,dac[kNdata-1],kHalfMax);
488        TLine *ldac_half   = new TLine(dac_half[ch],0,dac_half[ch],
489                                       (Double_t)kDrawMax);
490        hist_scurve[ch].SetMaximum(kDrawMax);
491        hist_scurve[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
492        hist_scurve[ch].Draw();
493        lhalf_max->Draw("same");
494        ldac_half->Draw("same");
495        if(dac_half[ch]>0){
496            dac_half_sum += dac_half[ch]-dac_ped[ch]; //relative absolute value
497            trigger++;
498        }
499       
500        /*** Differential plot from scurve fitting function  ***/
501        c5->cd(ch+1);
502        hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve,fit_max_scurve-2*step_dac);
503        Double_t dac_sfit_max=hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
504        //if(dac_sfit_max<=fit_min_scurve+4*step_dac){
505        if(dac_sfit_max<=fit_min_scurve+2*step_dac){
506            hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve+kFitSadj,fit_max_scurve-2*step_dac);
507            dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
508            //if(dac_sfit[ch]<=fit_min_scurve+kFitSadj+4*step_dac)dac_sfit[ch]=-1;
509            if(dac_sfit[ch]<fit_min_scurve+kFitSadj+step_dac){
510                dac_sfit[ch]=-1;
511                //dac_sfit[ch]=dac_half[ch];
512                flag_fit[ch] = -1;
513                cerr << "flag_fit[" << ch << "] " << flag_fit[ch]  << endl;
514            }
515        }else {
516            dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
517        }
518       
519        if(dac_sfit[ch]>fit_min_scurve+5){
520            dac_sfit_sum += dac_sfit[ch]-dac_ped[ch];
521            trigger_sfit++;
522        }else{
523            //dac_sfit[ch] = fit_min_scurve+10;
524            dac_sfit[ch] = -1;
525        }
526        //cerr << "dac_sfit[" << ch << "] " << dac_sfit[ch] << endl;
527        TLine *ldac_sfit      = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMaxFit);
528        TLine *ldac_sfit_copy = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMax);
529        ldac_sfit->SetLineColor(kMagenta);
530        ldac_sfit_copy->SetLineColor(kMagenta);
531        hist_diff_fit[ch].SetMaximum(kDrawMaxFit);
532        hist_diff_fit[ch].SetMinimum(0);
533        hist_diff_fit[ch].Draw();
534        ldac_sfit->Draw("same");
535       
536        /**** differential plot ****/
537        c2->cd(ch+1);
538        //hist_diff[ch].Rebin(8);
539        hist_diff[ch].Rebin(3);
540        fit_sdiff     = new TF1("fit_sdiff","gaus",dac_min,dac_max);
541        //fit_sdiff->SetLineWidth(1);
542        //if(ch==8 || ch==18 || ch==44 || ch==56)
543        mean_sum += mean_diff[ch]-dac_ped[ch];
544        trigger_mean++;
545       
546        // fit the derivated polynom fit
547        Double_t fit_par[3];
548        //hist_diff[ch].Fit(fit_sdiff,"RQ0","",fit_min,fit_max);
549        hist_diff[ch].Fit(fit_sdiff,"","",fit_min,fit_max);
550        hist_diff[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
551        hist_diff[ch].SetMinimum(0);
552        hist_diff[ch].SetMaximum(kDrawMaxDiff);
553        //hist_diff[ch].Rebin();
554        hist_diff[ch].Draw();
555        fit_sdiff->GetParameters(fit_par); // get the parameters for this fit
556        dac_spe[ch] = fit_par[1];
557        if(dac_spe[ch]>fit_min){
558            dac_spe_sum += dac_spe[ch]-dac_ped[ch];
559            //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
560            trigger_fit++;
561        }else{
562            dac_spe[ch] = fit_min+10;
563        }
564        //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
565        TLine *ldac_spe    = new TLine(dac_spe[ch],0,dac_spe[ch],
566                                       (Double_t)kDrawMaxDiff);
567        ldac_spe->SetLineColor(kRed);
568        TLine *lmean_diff  = new TLine(mean_diff[ch],0,mean_diff[ch],
569                                       (Double_t)kDrawMaxDiff);
570        lmean_diff->SetLineColor(kCyan);
571        TLine *lfit_min    = new TLine(fit_min,0,fit_min,kDrawMaxDiff);
572        lfit_min->SetLineColor(kGray);
573       
574        ldac_spe->Draw();
575        lmean_diff->Draw("same");
576        lfit_min->Draw("same");
577       
578        c3->cd(1);
579        if(ch==0){
580            hist_scurve[ch].SetMaximum(kDrawMax);
581            hist_scurve[ch].Draw();
582        }else{
583            /*if(ch==8 || ch==18 || ch==44 || ch==56)*/
584            hist_scurve[ch].Draw("same");
585        }
586        if(ch==kDrawCh){
587            c3->cd(2);
588            hist_scurve[ch].SetMaximum(kDrawMax);
589            hist_scurve[kDrawCh].Draw();
590            lhalf_max->Draw("same");
591            ldac_half->Draw("same");
592            ldac_sfit_copy->Draw("same");
593            c3->cd(3);
594            hist_diff[kDrawCh].Draw();
595            ldac_spe->Draw();
596            lmean_diff->Draw("same");
597            lfit_min->Draw("same");
598            c3->cd(4);
599            hist_diff_fit[kDrawCh].Draw();
600            ldac_sfit->Draw("same");
601        }
602        /*
603         if(ch<32){
604         c6->cd(ch/4+1);
605         hist_scurve[ch].SetMaximum(kDrawMax);
606         if(ch%4==0){
607         hist_scurve[ch].Draw();
608         }else{
609         hist_scurve[ch].Draw("same");
610         }
611         ldac_sfit_copy->Draw("same");
612         }else{
613         c7->cd(ch/4-7);
614         hist_scurve[ch].SetMaximum(kDrawMax);
615         if(ch%4==0){
616         hist_scurve[ch].Draw();
617         }else{
618         hist_scurve[ch].Draw("same");
619         }
620         ldac_sfit_copy->Draw("same");
621         }
622         */
623       
624        c4->cd();
625        if(ch==0){
626            hist_scurve_ped[ch].SetMaximum(kDrawMaxPed);
627            hist_scurve_ped[ch].Draw("l");
628        }else{
629            hist_scurve_ped[ch].Draw("lsame");
630        }
631       
632    }// end of second loop on all the channels
633   
634    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
635   
636    Float_t dac_ped_ave  ; //absolute value
637    Float_t dac_half_ave  ;
638    Float_t dac_spe_ave   ;
639    Float_t dac_sfit_ave;
640    Float_t mean_diff_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
641   
642    if(!checkAve){
643        dac_ped_ave   = dac_ped[kChRef]; //absolute value
644        dac_half_ave  = dac_half[kChRef];
645        dac_spe_ave   = dac_spe[kChRef];
646        dac_sfit_ave  = dac_sfit[kChRef];
647        //Float_t dac_sfit_ave  = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
648        mean_diff_ave = mean_diff[kChRef];
649    }else{
650        dac_ped_ave  = dac_ped_sum/(Float_t)kNch; //absolute value
651        dac_half_ave = dac_half_sum/(Float_t)trigger+dac_ped_ave;
652        dac_spe_ave  = dac_spe_sum/(Float_t)trigger_fit+dac_ped_ave;
653        mean_diff_ave = mean_sum/(Float_t)trigger_mean+dac_ped_ave;
654        dac_sfit_ave = dac_sfit_sum/(Float_t)trigger_sfit+dac_ped_ave;
655    }
656   
657   
658    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
659   
660    Double_t dac_half_spe;
661   
662    //trigger_mean=4;
663    if(dac_half_spe_fix!=-1) {
664        dac_half_spe = dac_half_spe_fix;
665    } else {
666        //Double_t dac_half_spe = dac_half_ave;
667        dac_half_spe = dac_spe_ave;
668    }
669   
670    TLine *ldac_ped_ave   = new TLine(dac_ped_ave,0,dac_ped_ave,
671                                      (Double_t)kDrawMax);
672    TLine *ldac_ped_ave_copy = new TLine(dac_ped_ave,0,dac_ped_ave,
673                                         (Double_t)kDrawMaxDiff);
674    TLine *ldac_half_ave  = new TLine(dac_half_ave,0,dac_half_ave,
675                                      (Double_t)kDrawMax);
676    TLine *ldac_half_spe  = new TLine(dac_half_spe,0,dac_half_spe,
677                                      (Double_t)kDrawMax);
678    TLine *lmean_diff_ave = new TLine(mean_diff_ave,0,mean_diff_ave,
679                                      (Double_t)kDrawMax);
680    TLine *ldac_sfit_ave  = new TLine(dac_sfit_ave,0,dac_sfit_ave,
681                                      (Double_t)kDrawMax);
682    ldac_ped_ave->SetLineColor(kOrange);
683    ldac_half_ave->SetLineColor(kRed);
684    lmean_diff_ave->SetLineColor(kCyan);
685    //ldac_half_spe->SetLineColor(kGreen);
686    ldac_sfit_ave->SetLineColor(kMagenta);
687   
688    c3->cd(1);
689    ldac_ped_ave->Draw("same");
690    ldac_half_ave->Draw("same");
691    // lhalf_max->Draw("same");  SDC (24/09/2013) :: Not declared
692    ldac_half_spe->Draw("same");
693    lmean_diff_ave->Draw("same");
694    ldac_sfit_ave->Draw("same");
695   
696    c3->cd(3);
697    ldac_ped_ave_copy->SetLineColor(kOrange);
698    ldac_ped_ave_copy->Draw("same");
699   
700    ofstream out(kFnameOut.Data());
701   
702    Int_t    gain_asic_diff, gain_asic_dacHalf, gain_asic_diff_mean, gain_asic_sfit;
703    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"
704    << endl;
705   
706   
707    // loop 3 on all channels to make final computations on the gain
708    for(ch=0;ch<kNch;ch++){
709        /*
710         if(ch==8 || ch==44 || ch==56 || ch==18){
711         gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
712         (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
713         }else{
714         gain_asic_diff_mean=0;
715         }
716         */
717        if(dac_half[ch]<0){
718            gain_asic_dacHalf   = 0;
719            gain_asic_diff      = 0;
720            gain_asic_diff_mean = 0;
721            gain_asic_sfit      = 0;
722        }else{
723            gain_asic_dacHalf  =(Int_t)((dac_half_ave-dac_ped_ave)/
724                                        (dac_half[ch]-dac_ped[ch])*gain_org[ch]);
725            gain_asic_diff     =(Int_t)((dac_half_spe-dac_ped_ave)/
726                                        (dac_spe[ch]-dac_ped[ch])*gain_org[ch]);
727            gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
728                                        (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
729            gain_asic_sfit     =(Int_t)((dac_sfit_ave-dac_ped_ave)/
730                                        (dac_sfit[ch]-dac_ped[ch])*gain_org[ch]);
731            //cerr << ch << " " << gain_org[ch] << endl;
732        }
733       
734       
735       
736        out << "OUT_TABLE    "
737        << ch << "    "
738        << gain_asic_dacHalf   << "    "
739        << gain_asic_diff      << "    "
740        << gain_asic_diff_mean << "    "
741        << gain_asic_sfit      << "    "
742        << dac_half[ch]        << "    "
743        << dac_ped[ch]         << "    "
744        << (Int_t)dac_spe[ch]  << "    "
745        << dac_half_ave        << "    "
746        << dac_half_spe        << "    "
747        << dac_ped_ave         << "    "
748        << dac_sfit_ave        << "    "
749        << sfit_chiS[ch]       << "    "
750        << ch                  << "    "
751        << flag_fit[ch]
752        << endl;
753        std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
754    }
755   
756   
757   
758    /*
759     for(ch=kNch-1;ch>=0;ch--){
760     cout << "OUT_TABLE_INV    "
761     << (Int_t)(dac_half_ave/(dac_half[ch]-dac_ped)*64)
762     << endl;
763     }
764     */
765   
766   
767}
768
769//--------------------------------------------------------------------------------------------------------
770
771//--------------------------------------------------------------------------------------------------------
772// Analysis Constructor
773//--------------------------------------------------------------------------------------------------------
774Analysis::Analysis():ch(0),dac_half_sum(0),dac_spe_sum(0),dac_sfit_sum(0),diff_sum(0),mean_sum(0),maxBinSum(0)
775{
776    if(DEBUG)
777        cout << " >>>>>>> START Analysis Constructor " << endl;
778   
779    /*********************** SIGNAL RUN ********************/
780    input             = new TString[kNch];
781    input_ped         = new TString[kNch];
782   
783    // Init sums
784    dac_half_sum      = 0;
785    dac_ped_sum       = 0;
786    dac_spe_sum       = 0;
787    dac_sfit_sum      = 0;
788    diff_sum          = 0;
789    mean_sum          = 0;
790    maxBinSum         = 0;
791   
792    if(DEBUG)
793        cout << " >>>>>>> END Analysis Constructor " << endl;
794}
795//--------------------------------------------------------------------------------------------------------
796
797
798//--------------------------------------------------------------------------------------------------------
799// Analysis destructor
800//--------------------------------------------------------------------------------------------------------
801Analysis::~Analysis()
802{
803   
804    if(DEBUG)
805        cout << " >>>>>>> START Analysis Destructor " << endl;
806   
807    if(input!=0) delete input;
808    if(input_ped!=0) delete input_ped;
809   
810    if(DEBUG)
811        cout << " >>>>>>> END Analysis Destructor " << endl;
812   
813}
814//--------------------------------------------------------------------------------------------------------
815
816
817//--------------------------------------------------------------------------------------------------------
818// Main public function of Analysis
819//--------------------------------------------------------------------------------------------------------
820/*
821 void Analysis::PlotScurve(
822 Int_t    kDrawCh=31,
823 TString  fname=kFname,
824 TString  fname_ped=kFnamePed, //ASIC_D ped_file
825 Double_t dac_half_spe_fix = -1,
826 TString  fname_gain="gain-org.txt",
827 Int_t    fit_min=kFitMin,
828 Int_t    fit_max=kFitMax,
829 Int_t    fit_min_scurve=kFitMinS,
830 Int_t    fit_max_scurve=kFitMaxS,
831 Bool_t   kBoard=0, //0:ASIC_C, 1:ASIC_D
832 Int_t    kChRef=kRefCh,//ASIC_C
833 Bool_t   checkAve=0 //0:use ref pix, 1:check average
834 )
835 */
836void Analysis:: PlotScurve(
837                           Int_t    kDrawCh,
838                           TString  fname,
839                           TString  fname_ped, //ASIC_D ped_file
840                           Double_t dac_half_spe_fix,
841                           TString  fname_gain,
842                           Int_t    fit_min,
843                           Int_t    fit_max,
844                           Int_t    fit_min_scurve,
845                           Int_t    fit_max_scurve,
846                           Bool_t   kBoard, //0:ASIC_C, 1:ASIC_D
847                           Int_t    kChRef,//ASIC_C
848                           Bool_t   checkAve //0:use ref pix, 1:check average
849                           )
850
851
852
853{
854   
855    if(DEBUG)
856    {
857        cout << " >>>>>>> START Analysis::PlotScurve " << endl;
858        cout << " \t - kDrawCh = " << kDrawCh << endl;
859        cout << " \t - fname = " << fname << endl;
860        cout << " \t - fname_ped = " << fname_ped <<endl;
861        cout << " \t - dac_half_spe_fix = " << dac_half_spe_fix << endl;
862       
863        cout << " \t - fname_gain = " << fname_gain << endl;
864        cout << " \t - fit_min = " << fit_min << endl;
865        cout << " \t - fit_max = " << fit_max << endl;
866        cout << " \t - fit_min_scurve = " << fit_min_scurve << endl;
867        cout << " \t - fit_max_scurve = " << fit_max_scurve << endl;
868        cout << " \t - kBoard = " << kBoard << " 0:ASIC_C, 1:ASIC_D " << endl;
869        cout << " \t - kChRef = " << kChRef << " ASIC_C " << endl;
870        cout << " \t - checkAve = " << checkAve << " 0:use ref pix, 1:check average " << endl;
871    }
872    // reset all variables
873    Init();
874   
875    //OpenGainFile(fname_gain);
876   
877    FillHistoName();
878   
879    ReadGainFile(fname_gain,kBoard);
880   
881    // Open the pedestal file and the S-Curve files for each channel and test if the file can be open
882   
883    trigger=0;
884    trigger_fit=0;
885    trigger_mean=0;
886    trigger_sfit=0;
887   
888    // Big loop on all the channels : second loop
889    //--------------------------------------------
890    // build the input data file name for each channel fname-chXX.dat
891    for(ch=0;ch<kNch;ch++)
892    {
893        ProcessSingleChannel(kDrawCh,fname,fname_ped,dac_half_spe_fix,fname_gain,
894                             fit_min,fit_max,fit_min_scurve,fit_max_scurve,kBoard,kChRef);
895    }
896   
897   
898   
899    ComputeDACAverage(dac_half_spe_fix ,kChRef,checkAve);
900    OutputGain(kFnameOut);
901   
902    if(DEBUG)
903        cout << " >>>>>>> END Analysis::PlotSCurve " << endl;
904   
905}
906//--------------------------------------------------------------------------------------------------------
907
908
909//--------------------------------------------------------------------------------------------------------
910// Initialize the program by reseting counters or variables and create objects
911//--------------------------------------------------------------------------------------------------------
912void Analysis::Init()
913{
914   
915    if(DEBUG)
916        cout << " >>>>>>> START Analysis::Init " << endl;
917   
918    //if(input!=0) delete input;
919    //if(input_ped!=0) delete input_ped;
920   
921    // recreate the name
922    input             = new TString[kNch];
923    input_ped         = new TString[kNch];
924   
925    // Init sums
926    dac_half_sum      = 0;
927    dac_ped_sum       = 0;
928    dac_spe_sum       = 0;
929    dac_sfit_sum      = 0;
930    diff_sum          = 0;
931    mean_sum          = 0;
932    maxBinSum         = 0;
933   
934   
935    // create canvas
936   
937    c1                = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
938    c1->Divide(11,6);
939    c2                = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
940    c2->Divide(11,6);
941    c5                = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
942    c5->Divide(11,6);
943    c3                = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
944    c3->Divide(4,1);
945    c4                = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
946   
947   
948    // create histos and function
949    hist_scurve       = new TH1D[kNch]; // smoothed S-curves
950    hist_scurve_ped   = new TH1D[kNch]; // S-curves for pedestal runs
951    hist_diff         = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
952    hist_diff_fit     = new TH1D[kNch]; // differentiated histograms of fitted s-curves
953    fit_sdiff         = new TF1[kNch];  // function fitted on differentiated histograms of fitted s-curves
954   
955   
956    for(ch=0;ch<kNch;ch++){
957        dac_half[ch]=0;
958        dac_ped[ch]=0;
959        dac_spe[ch]=0;
960        dac_sfit[ch]=0;
961        gain_org[ch]=0;
962        mean_diff[ch]=0;
963        sfit_chiS[ch]=0;
964        flag_fit[ch]=0;
965    }
966   
967    if(DEBUG)
968        cout << " >>>>>>> END Analysis::Init " << endl;
969   
970   
971}
972//--------------------------------------------------------------------------------------------------------
973
974//--------------------------------------------------------------------------------------------------------
975// Initialize the program by reseting counters or variables and create objects
976//--------------------------------------------------------------------------------------------------------
977void Analysis::OpenGainFile(const TString fname_gain)
978{
979    if(DEBUG)
980        cout << " >>>>>>> START Analysis::OpenGainFile " << endl;
981   
982   
983   
984    TString input_gain(fname_gain);
985    ifstream in_gain(input_gain.Data());
986    if(!in_gain.good()){
987        cerr << "File " << input_gain << " open failed!" << endl;
988        exit(-1);
989    }
990    in_gain.close();
991    in_gain.clear();
992   
993    if(DEBUG)
994        cout << " >>>>>>> END Analysis::OpenGainFile " << endl;
995   
996   
997}
998//--------------------------------------------------------------------------------------------------------
999
1000//--------------------------------------------------------------------------------------------------------
1001// Fill Histogram name
1002//--------------------------------------------------------------------------------------------------------
1003void Analysis::FillHistoName()
1004{
1005   
1006    if(DEBUG)
1007        cout << " >>>>>>> START Analysis::FillHistoName " << endl;
1008   
1009   
1010   
1011    // first loop on channels to fill histonames
1012    //----------------------
1013    for(ch=0;ch<kNch;ch++){
1014       
1015        name_hist[ch]           ="hist";
1016        name_hist[ch]          +=ch;
1017        name_hist_ped[ch]       ="hist_ped";
1018        name_hist_ped[ch]      +=ch;
1019        name_hist_diff[ch]      ="hist_diff";
1020        name_hist_diff[ch]     +=ch;
1021        name_hist_diff_fit[ch]  ="hist_diff_fit";
1022        name_hist_diff_fit[ch] +=ch;
1023        name_fit_diff[ch]       ="fit_diff";
1024        name_fit_diff[ch]      +=ch;
1025    }
1026   
1027    if(DEBUG)
1028        cout << " >>>>>>> END Analysis::FillHistoName " << endl;
1029   
1030   
1031   
1032}
1033//--------------------------------------------------------------------------------------------------------
1034
1035
1036//--------------------------------------------------------------------------------------------------------
1037// Read the gain
1038// input : the gain filename
1039// output the gain_org[ch] array
1040//--------------------------------------------------------------------------------------------------------
1041
1042void Analysis::ReadGainFile(const TString fname_gain,Bool_t kBoard)
1043{
1044   
1045    if(DEBUG)
1046    {
1047        cout << " >>>>>>> START Analysis::ReadGainFile " << endl;
1048        cout << " \t - fname_gain = " << fname_gain << endl;
1049        cout << " \t - kBoard = " << kBoard << " 0:ASIC_C, 1:ASIC_D " << endl;
1050       
1051    }
1052   
1053   
1054   
1055    TString input_gain(fname_gain);
1056    ifstream in_gain(input_gain.Data());
1057    if(!in_gain.good()){
1058        cerr << "File " << input_gain << " open failed!" << endl;
1059        exit(-1);
1060    }
1061    // read gain
1062    Double_t da;
1063    Int_t count=0;
1064    ch=0;
1065    while(in_gain >> da){
1066        if(!kBoard && (count>=128 && count<192)){
1067            gain_org[ch]=da;
1068            cerr << "!kBoard, count, gain_org[" << ch << "] "
1069            << !kBoard << " " << count << " " << gain_org[ch] << endl;
1070            ch++;
1071        }
1072        if(kBoard && (count>=192 && count<256)){
1073            gain_org[ch]=da;
1074            cerr << "kBoard, count, gain_org[" << ch << "] "
1075            << kBoard << " " << count << " " << gain_org[ch] << endl;
1076            ch++;
1077        }
1078        count++;
1079    }// end if first loop on all channels
1080    in_gain.close();
1081    in_gain.clear();
1082   
1083    if(DEBUG)
1084    {
1085        cout << " >>>>>>> END Analysis::ReadGainFile " << endl;
1086       
1087    }
1088   
1089}
1090//--------------------------------------------------------------------------------------------------------
1091
1092
1093
1094//--------------------------------------------------------------------------------------------------------
1095// ProcessSingleChannel
1096//--------------------------------------------------------------------------------------------------------
1097void Analysis::ProcessSingleChannel(
1098                                    Int_t    kDrawCh,
1099                                    TString  fname,
1100                                    TString  fname_ped, //ASIC_D ped_file
1101                                    Double_t dac_half_spe_fix,
1102                                    TString  fname_gain,
1103                                    Int_t    fit_min,
1104                                    Int_t    fit_max,
1105                                    Int_t    fit_min_scurve,
1106                                    Int_t    fit_max_scurve,
1107                                    Bool_t   kBoard, //0:ASIC_C, 1:ASIC_D
1108                                    Int_t    kChRef//ASIC_C
1109                                    //  Bool_t   checkAve //0:use ref pix, 1:check average)
1110                                    )
1111{
1112   
1113   
1114    if(DEBUG)
1115    {
1116        cout << " >>>>>>> START Analysis::ProcessSingleChannel ch=" << ch << endl;
1117        cout << " \t - kDrawCh = " << kDrawCh << endl;
1118        cout << " \t - fname = " << fname << endl;
1119        cout << " \t - fname_ped = " << fname_ped <<endl;
1120        cout << " \t - dac_half_spe_fix = " << dac_half_spe_fix << endl;
1121       
1122        cout << " \t - fname_gain = " << fname_gain << endl;
1123        cout << " \t - fit_min = " << fit_min << endl;
1124        cout << " \t - fit_max = " << fit_max << endl;
1125        cout << " \t - fit_min_scurve = " << fit_min_scurve << endl;
1126        cout << " \t - fit_max_scurve = " << fit_max_scurve << endl;
1127        cout << " \t - kBoard = " << kBoard << " 0:ASIC_C, 1:ASIC_D " << endl;
1128        cout << " \t - kChRef = " << kChRef << " ASIC_C " << endl;
1129        //    cout << " \t - checkAve = " << checkAve << " 0:use ref pix, 1:check average " << endl;L
1130    }
1131   
1132   
1133   
1134    input[ch] += fname;
1135    input[ch] += "-ch";
1136    input[ch] += ch;
1137    input[ch] += ".dat";
1138    ifstream in(input[ch].Data());
1139    if(!in.good()){
1140        cerr << "File " << input[ch] << " open failed!" << endl;
1141        return;
1142    }
1143   
1144    // build the input data file name for each channel fname_ped-chXX.dat
1145    input_ped[ch] += fname_ped;
1146    input_ped[ch] += "-ch";
1147    input_ped[ch] += ch;
1148    input_ped[ch] += ".dat";
1149    ifstream in_ped(input_ped[ch].Data());
1150    if(!in_ped.good()){
1151        cerr << "File " << input_ped[ch] << " open failed!" << endl;
1152        return;
1153    }
1154   
1155   
1156   
1157    //10DAC~10mV
1158    // determine the DAC range in dac_min, dac_max, step_dac from file input_ped
1159    // and count the number of events
1160    Double_t da, dat, dac_min, dac_max, step_dac;
1161    Int_t count = 0;
1162    while (in >> da >> dat ){
1163        if(count==0) dac_min  = da;
1164        if(count==1) step_dac = da - dac_min;
1165        dac_max=da;
1166        count++;
1167    }
1168    in.close();
1169    in.clear();
1170   
1171    /*
1172     if(ch==43 || ch==47){//for T1 ASIC_C
1173     cerr << "CHECK " << ch << " " << fit_min_scurve << endl;
1174     Int_t fit_min_scurve=140;
1175     } else {
1176     Int_t fit_min_scurve=kFitMinS;
1177     }
1178     */
1179   
1180    // correct the scale for DAC range
1181    if(kXmin<dac_min)kXmin=dac_min;
1182    const Int_t kNdata    = count;
1183    count = (fit_max_scurve-fit_min_scurve)/step_dac;
1184    const Int_t kNdataFit = count;
1185    if(ch==0){
1186        cerr << "kNdata, dac_min, dac_max, step_dac: "
1187        << kNdata << " " << dac_min << " " << dac_max << " "
1188        << step_dac << endl;
1189    }
1190   
1191    Double_t dac_min_ped, dac_max_ped, step_dac_ped;
1192    count = 0;
1193    while (in_ped >> da >> dat ){
1194        if(count==0) dac_min_ped  = da;
1195        if(count==1) step_dac_ped = da - dac_min_ped;
1196        dac_max_ped=da;
1197        count++;
1198    }
1199   
1200    in_ped.close();
1201    in_ped.clear();
1202   
1203   
1204    const Int_t kNdataPed = count;
1205    if(ch==0){
1206        cerr << "kNdataPed, dac_min_ped, dac_max_ped, step_dac_ped: "
1207        << kNdataPed << " " << dac_min_ped << " " << dac_max_ped
1208        << " " << step_dac_ped << endl;
1209    }
1210   
1211    // book histograms
1212    // SDC: 27/09/2013 hist_scurve is a pointer on an array of THD, thus the syntax is correct (no new operator)
1213   
1214   
1215    hist_scurve[ch]      = TH1D(name_hist[ch],name_hist[ch],kNdata-1,
1216                                dac_min,dac_max);
1217    hist_scurve_ped[ch]  = TH1D(name_hist_ped[ch],name_hist_ped[ch],
1218                                kNdataPed-1,dac_min_ped,dac_max_ped);
1219    hist_diff[ch]        = TH1D(name_hist_diff[ch],name_hist_diff[ch],
1220                                kNdata-2,dac_min,dac_max);
1221    hist_diff_fit[ch]    = TH1D(name_hist_diff_fit[ch],name_hist_diff_fit[ch],
1222                                kNdataFit,fit_min_scurve,fit_max_scurve);
1223   
1224   
1225   
1226   
1227   
1228    // read the data for that channel for normal data
1229    in.open(input[ch].Data());
1230    Double_t       data[kNdata], dac[kNdata], dac_fit[kNdataFit], data_fit[kNdataFit];
1231    while (in >> da >> dat ){
1232        hist_scurve[ch].Fill(da,dat); // fill the histograms of the s-curve
1233    }
1234    in.close();
1235    in.clear();
1236   
1237    // read the data for pedestal data
1238    in_ped.open(input_ped[ch].Data());
1239    while (in_ped >> da >> dat ){
1240        hist_scurve_ped[ch].Fill(da,dat); // fill the histogram for pedestal data
1241    }
1242    in_ped.close();
1243    in_ped.clear();
1244   
1245    //dac_ped[ch]  = hist_scurve_ped[ch].GetMaximumBin()+dac_min_ped; //absolute value
1246    dac_ped[ch]  = hist_scurve_ped[ch].GetBinCenter(hist_scurve_ped[ch].GetMaximumBin()); //absolute value
1247    dac_ped_sum += dac_ped[ch]; //absolute value
1248    //cerr << "dac_ped[" << ch << "] dac_ped_sum " << dac_ped[ch]
1249    //<< " " << dac_ped_sum << endl;
1250   
1251   
1252    hist_scurve[ch].Smooth(3); // Smooth the S curve
1253    Int_t binx=0;
1254    //hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
1255    //                                  dac_max-dac_min+1,kHalfMax/2.);
1256   
1257    // ??????
1258    hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
1259                                      dac_max-dac_min+1,2000);
1260    dac_half[ch] = hist_scurve[ch].GetBinCenter(binx); //absolute value
1261    if(dac_half[ch]<kXmin)dac_half[ch]=-1;
1262   
1263    // put the histograms of the smoothed S-Curves into two arrays (dac, data)
1264    //-------------------------------------------------------------------------
1265    Int_t i=0;
1266    for(i=0;i<kNdata;i++){
1267        dac[i]  = hist_scurve[ch].GetBinCenter(i+1);
1268        data[i] = hist_scurve[ch].GetBinContent(i+1);
1269    }
1270   
1271    /*** compute derive diff on S-Curve ***/
1272    Double_t diff[kNdata-1], dac_diff[kNdata-1];
1273    Int_t trigger_diff=0;
1274    for(i=0;i<kNdata-1;i++){
1275        diff[i]     = (data[i]-data[i+1])/step_dac;
1276        dac_diff[i] = dac[i]+step_dac/2;
1277        hist_diff[ch].Fill(dac_diff[i],diff[i]); // fills differentiated histograms of smoothed s-curves which is the single photoelectron spectrum
1278        if(dac_diff[i]>=fit_min){
1279            diff_sum += diff[i]*dac_diff[i];
1280            //cerr << "diff[" << i << "], dac_diff[" << i << "] "
1281            //<< diff[i] << " " << dac_diff[i] << endl;
1282            trigger_diff+=diff[i];
1283        }
1284    }
1285    //mean_diff[ch] = diff_sum/trigger_diff;
1286    hist_diff[ch].GetXaxis()->SetRangeUser(fit_min,fit_max);
1287    mean_diff[ch] = hist_diff[ch].GetMean();
1288    //cerr << "trigger_diff, mean_diff[" << ch << "] "
1289    //<< trigger_diff << " " << mean_diff[ch] << endl;
1290   
1291   
1292    /******************* GRAPHICS *****************/
1293   
1294   
1295    // Now start the analysis which consist in finding the peak in the differentiated spectrum
1296    c1->cd(ch+1);
1297    /*** fit scurves  ***/
1298    cerr << "ch " << ch;
1299   
1300    //    SDC 24/09/2013 :: Need to declare outside this variable in a compiled program
1301    TF1 *fit_pol ;
1302   
1303   
1304    // Fit a polynomial curve in a given range on the s-curve directly. This is a kind of DEBUG procedure.
1305   
1306    //if(ch==43 || ch==47){ // SDC: 25/09/2013 Why a special channel : the fit pol4 is performed for another DAC range
1307    /*   if(ch==170){
1308     fit_pol      = new TF1("fit_pol","pol4",150,fit_max_scurve);
1309     hist_scurve[ch].Fit(fit_pol,"","",150,fit_max_scurve);
1310     }else{ */
1311    fit_pol      = new TF1("fit_pol","pol4",fit_min_scurve,fit_max_scurve);
1312    hist_scurve[ch].Fit(fit_pol,"","",fit_min_scurve,fit_max_scurve);
1313    // }
1314   
1315   
1316   
1317    Double_t fit_par_pol4[5]; // container for fit parameter results
1318    fit_pol->GetParameters(fit_par_pol4);
1319    sfit_chiS[ch] =  fit_pol->GetChisquare(); // chi squared
1320    TF1 *func               = hist_scurve[ch].GetFunction("fit_pol");
1321    Double_t data_diff_fit;
1322    for(i=0;i<kNdataFit;i++){
1323        //if(ch==43 || ch==47){
1324        if(ch==170){
1325            dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+150-dac_min;
1326            data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+150-dac_min); // get the value of the fitted function
1327        }else{
1328            dac_fit[i]  = hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min;
1329            data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min); // get the value of the fitted function
1330            //cerr << "dac_fit[" << i << "], data_fit[" << i << "] "
1331            //<< dac_fit[i] << " " << data_fit[i] << endl;
1332        }
1333    }
1334   
1335    // derive the fitted function on the S-curve to obtain the SPE
1336    Double_t dac_diff_fit;
1337    for(i=0;i<kNdataFit-1;i++){
1338        data_diff_fit = (data_fit[i]-data_fit[i+1])/step_dac;
1339        hist_diff_fit[ch].Fill(dac_fit[i],data_diff_fit); // Fill differentiated fitted S-curve
1340    }
1341   
1342    TLine *lhalf_max   = new TLine(dac[0],kHalfMax,dac[kNdata-1],kHalfMax);
1343    TLine *ldac_half   = new TLine(dac_half[ch],0,dac_half[ch],
1344                                   (Double_t)kDrawMax);
1345    hist_scurve[ch].SetMaximum(kDrawMax);
1346    hist_scurve[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
1347    hist_scurve[ch].Draw();
1348    lhalf_max->Draw("same");
1349    ldac_half->Draw("same");
1350    if(dac_half[ch]>0){
1351        dac_half_sum += dac_half[ch]-dac_ped[ch]; //relative absolute value
1352        trigger++;
1353    }
1354   
1355    /*** Differential plot from scurve fitting function  ***/
1356    c5->cd(ch+1);
1357    hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve,fit_max_scurve-2*step_dac);
1358    Double_t dac_sfit_max=hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
1359    //if(dac_sfit_max<=fit_min_scurve+4*step_dac){
1360    if(dac_sfit_max<=fit_min_scurve+2*step_dac){
1361        hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve+kFitSadj,fit_max_scurve-2*step_dac);
1362        dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
1363        //if(dac_sfit[ch]<=fit_min_scurve+kFitSadj+4*step_dac)dac_sfit[ch]=-1;
1364        if(dac_sfit[ch]<fit_min_scurve+kFitSadj+step_dac){
1365            dac_sfit[ch]=-1;
1366            //dac_sfit[ch]=dac_half[ch];
1367            flag_fit[ch] = -1;
1368            cerr << "flag_fit[" << ch << "] " << flag_fit[ch]  << endl;
1369        }
1370    }else {
1371        dac_sfit[ch]  =  hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
1372    }
1373   
1374    if(dac_sfit[ch]>fit_min_scurve+5){
1375        dac_sfit_sum += dac_sfit[ch]-dac_ped[ch];
1376        trigger_sfit++;
1377    }else{
1378        //dac_sfit[ch] = fit_min_scurve+10;
1379        dac_sfit[ch] = -1;
1380    }
1381    //cerr << "dac_sfit[" << ch << "] " << dac_sfit[ch] << endl;
1382    TLine *ldac_sfit      = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMaxFit);
1383    TLine *ldac_sfit_copy = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMax);
1384    ldac_sfit->SetLineColor(kMagenta);
1385    ldac_sfit_copy->SetLineColor(kMagenta);
1386    hist_diff_fit[ch].SetMaximum(kDrawMaxFit);
1387    hist_diff_fit[ch].SetMinimum(0);
1388    hist_diff_fit[ch].Draw();
1389    ldac_sfit->Draw("same");
1390   
1391    /**** differential plot ****/
1392    c2->cd(ch+1);
1393    //hist_diff[ch].Rebin(8);
1394    hist_diff[ch].Rebin(3);
1395    fit_sdiff     = new TF1("fit_sdiff","gaus",dac_min,dac_max);
1396    //fit_sdiff->SetLineWidth(1);
1397    //if(ch==8 || ch==18 || ch==44 || ch==56)
1398    mean_sum += mean_diff[ch]-dac_ped[ch];
1399    trigger_mean++;
1400   
1401    // fit the derivated polynom fit
1402    Double_t fit_par[3];
1403    //hist_diff[ch].Fit(fit_sdiff,"RQ0","",fit_min,fit_max);
1404    hist_diff[ch].Fit(fit_sdiff,"","",fit_min,fit_max);
1405    hist_diff[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
1406    hist_diff[ch].SetMinimum(0);
1407    hist_diff[ch].SetMaximum(kDrawMaxDiff);
1408    //hist_diff[ch].Rebin();
1409    hist_diff[ch].Draw();
1410    fit_sdiff->GetParameters(fit_par); // get the parameters for this fit
1411    dac_spe[ch] = fit_par[1];
1412    if(dac_spe[ch]>fit_min){
1413        dac_spe_sum += dac_spe[ch]-dac_ped[ch];
1414        //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
1415        trigger_fit++;
1416    }else{
1417        dac_spe[ch] = fit_min+10;
1418    }
1419    //cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
1420    TLine *ldac_spe    = new TLine(dac_spe[ch],0,dac_spe[ch],
1421                                   (Double_t)kDrawMaxDiff);
1422    ldac_spe->SetLineColor(kRed);
1423    TLine *lmean_diff  = new TLine(mean_diff[ch],0,mean_diff[ch],
1424                                   (Double_t)kDrawMaxDiff);
1425    lmean_diff->SetLineColor(kCyan);
1426    TLine *lfit_min    = new TLine(fit_min,0,fit_min,kDrawMaxDiff);
1427    lfit_min->SetLineColor(kGray);
1428   
1429    ldac_spe->Draw();
1430    lmean_diff->Draw("same");
1431    lfit_min->Draw("same");
1432   
1433    c3->cd(1);
1434    if(ch==0){
1435        hist_scurve[ch].SetMaximum(kDrawMax);
1436        hist_scurve[ch].Draw();
1437    }else{
1438        /*if(ch==8 || ch==18 || ch==44 || ch==56)*/
1439        hist_scurve[ch].Draw("same");
1440    }
1441    if(ch==kDrawCh){
1442        c3->cd(2);
1443        hist_scurve[ch].SetMaximum(kDrawMax);
1444        hist_scurve[kDrawCh].Draw();
1445        lhalf_max->Draw("same");
1446        ldac_half->Draw("same");
1447        ldac_sfit_copy->Draw("same");
1448        c3->cd(3);
1449        hist_diff[kDrawCh].Draw();
1450        ldac_spe->Draw();
1451        lmean_diff->Draw("same");
1452        lfit_min->Draw("same");
1453        c3->cd(4);
1454        hist_diff_fit[kDrawCh].Draw();
1455        ldac_sfit->Draw("same");
1456    }
1457    /*
1458     if(ch<32){
1459     c6->cd(ch/4+1);
1460     hist_scurve[ch].SetMaximum(kDrawMax);
1461     if(ch%4==0){
1462     hist_scurve[ch].Draw();
1463     }else{
1464     hist_scurve[ch].Draw("same");
1465     }
1466     ldac_sfit_copy->Draw("same");
1467     }else{
1468     c7->cd(ch/4-7);
1469     hist_scurve[ch].SetMaximum(kDrawMax);
1470     if(ch%4==0){
1471     hist_scurve[ch].Draw();
1472     }else{
1473     hist_scurve[ch].Draw("same");
1474     }
1475     ldac_sfit_copy->Draw("same");
1476     }
1477     */
1478   
1479    c4->cd();
1480    if(ch==0){
1481        hist_scurve_ped[ch].SetMaximum(kDrawMaxPed);
1482        hist_scurve_ped[ch].Draw("l");
1483    }else{
1484        hist_scurve_ped[ch].Draw("lsame");
1485    }
1486   
1487   
1488    if(DEBUG)
1489    {
1490        cout << " >>>>>>> END Analysis::ProcessSingleChannel " << endl;
1491    }
1492   
1493}// end of process single  channels
1494//---------------------------------------------------------------------------------------------------
1495
1496
1497//---------------------------------------------------------------------------------------------------
1498// Compute some averages
1499//---------------------------------------------------------------------------------------------------
1500void Analysis::ComputeDACAverage(Double_t dac_half_spe_fix ,Int_t    kChRef, Bool_t   checkAve )//0:use ref pix, 1:check average
1501{
1502   
1503   
1504    if(DEBUG)
1505    {
1506        cout << " >>>>>>> START Analysis::ComputeDACAverage " << endl;
1507       
1508        cout << " \t - dac_half_spe_fix = " << dac_half_spe_fix << endl;
1509        cout << " \t - kChRef = " << kChRef << " ASIC_C " << endl;
1510        cout << " \t - checkAve = " << checkAve << " 0:use ref pix, 1:check average " << endl;
1511    }
1512   
1513    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
1514   
1515   
1516    mean_diff_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
1517   
1518    if(!checkAve){
1519        dac_ped_ave   = dac_ped[kChRef]; //absolute value
1520        dac_half_ave  = dac_half[kChRef];
1521        dac_spe_ave   = dac_spe[kChRef];
1522        dac_sfit_ave  = dac_sfit[kChRef];
1523        //Float_t dac_sfit_ave  = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
1524        mean_diff_ave = mean_diff[kChRef];
1525    }else{
1526        dac_ped_ave  = dac_ped_sum/(Float_t)kNch; //absolute value
1527        dac_half_ave = dac_half_sum/(Float_t)trigger+dac_ped_ave;
1528        dac_spe_ave  = dac_spe_sum/(Float_t)trigger_fit+dac_ped_ave;
1529        mean_diff_ave = mean_sum/(Float_t)trigger_mean+dac_ped_ave;
1530        dac_sfit_ave = dac_sfit_sum/(Float_t)trigger_sfit+dac_ped_ave;
1531    }
1532   
1533   
1534    // SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
1535   
1536   
1537    //trigger_mean=4;
1538    if(dac_half_spe_fix!=-1) {
1539        dac_half_spe = dac_half_spe_fix;
1540    } else {
1541        //Double_t dac_half_spe = dac_half_ave;
1542        dac_half_spe = dac_spe_ave;
1543    }
1544   
1545    TLine *ldac_ped_ave   = new TLine(dac_ped_ave,0,dac_ped_ave,
1546                                      (Double_t)kDrawMax);
1547    TLine *ldac_ped_ave_copy = new TLine(dac_ped_ave,0,dac_ped_ave,
1548                                         (Double_t)kDrawMaxDiff);
1549    TLine *ldac_half_ave  = new TLine(dac_half_ave,0,dac_half_ave,
1550                                      (Double_t)kDrawMax);
1551    TLine *ldac_half_spe  = new TLine(dac_half_spe,0,dac_half_spe,
1552                                      (Double_t)kDrawMax);
1553    TLine *lmean_diff_ave = new TLine(mean_diff_ave,0,mean_diff_ave,
1554                                      (Double_t)kDrawMax);
1555    TLine *ldac_sfit_ave  = new TLine(dac_sfit_ave,0,dac_sfit_ave,
1556                                      (Double_t)kDrawMax);
1557    ldac_ped_ave->SetLineColor(kOrange);
1558    ldac_half_ave->SetLineColor(kRed);
1559    lmean_diff_ave->SetLineColor(kCyan);
1560    //ldac_half_spe->SetLineColor(kGreen);
1561    ldac_sfit_ave->SetLineColor(kMagenta);
1562   
1563    c3->cd(1);
1564    ldac_ped_ave->Draw("same");
1565    ldac_half_ave->Draw("same");
1566    // lhalf_max->Draw("same");  SDC (24/09/2013) :: Not declared
1567    ldac_half_spe->Draw("same");
1568    lmean_diff_ave->Draw("same");
1569    ldac_sfit_ave->Draw("same");
1570   
1571    c3->cd(3);
1572    ldac_ped_ave_copy->SetLineColor(kOrange);
1573    ldac_ped_ave_copy->Draw("same");
1574   
1575    if(DEBUG)
1576    {
1577        cout << " >>>>>>> END Analysis::ComputeDACAverage " << endl;
1578    }
1579   
1580   
1581}
1582
1583//---------------------------------------------------------------------------------------------------
1584
1585//---------------------------------------------------------------------------------------------------
1586void Analysis::OutputGain(const TString fname_out)
1587{
1588   
1589    if(DEBUG)
1590    {
1591        cout << " >>>>>>> START Analysis::OutputGain " << endl;
1592        cout << " \t - fname_out = " << fname_out << endl;
1593    }
1594   
1595    ofstream out(fname_out.Data());
1596   
1597    Int_t    gain_asic_diff, gain_asic_dacHalf, gain_asic_diff_mean, gain_asic_sfit;
1598    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"
1599    << endl;
1600   
1601   
1602    // loop 3 on all channels to make final computations on the gain
1603    for(ch=0;ch<kNch;ch++){
1604        /*
1605         if(ch==8 || ch==44 || ch==56 || ch==18){
1606         gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
1607         (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
1608         }else{
1609         gain_asic_diff_mean=0;
1610         }
1611         */
1612        if(dac_half[ch]<0){
1613            gain_asic_dacHalf   = 0;
1614            gain_asic_diff      = 0;
1615            gain_asic_diff_mean = 0;
1616            gain_asic_sfit      = 0;
1617        }else{
1618            gain_asic_dacHalf  =(Int_t)((dac_half_ave-dac_ped_ave)/
1619                                        (dac_half[ch]-dac_ped[ch])*gain_org[ch]);
1620            gain_asic_diff     =(Int_t)((dac_half_spe-dac_ped_ave)/
1621                                        (dac_spe[ch]-dac_ped[ch])*gain_org[ch]);
1622            gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
1623                                        (mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
1624            gain_asic_sfit     =(Int_t)((dac_sfit_ave-dac_ped_ave)/
1625                                        (dac_sfit[ch]-dac_ped[ch])*gain_org[ch]);
1626            //cerr << ch << " " << gain_org[ch] << endl;
1627        }
1628       
1629       
1630       
1631        out << "OUT_TABLE    "
1632        << ch << "    "
1633        << gain_asic_dacHalf   << "    "
1634        << gain_asic_diff      << "    "
1635        << gain_asic_diff_mean << "    "
1636        << gain_asic_sfit      << "    "
1637        << dac_half[ch]        << "    "
1638        << dac_ped[ch]         << "    "
1639        << (Int_t)dac_spe[ch]  << "    "
1640        << dac_half_ave        << "    "
1641        << dac_half_spe        << "    "
1642        << dac_ped_ave         << "    "
1643        << dac_sfit_ave        << "    "
1644        << sfit_chiS[ch]       << "    "
1645        << ch                  << "    "
1646        << flag_fit[ch]       
1647        << endl;
1648       
1649        if(DEBUG>2)
1650        {
1651            std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
1652        }
1653    }
1654    out.close();
1655   
1656    if(DEBUG)
1657    {
1658        cout << " >>>>>>> END Analysis::OutputGain " << endl;
1659       
1660    }
1661   
1662   
1663}
1664//--------------------------------------------------------------------------------------------------------
1665
1666
1667
1668
1669
1670
1671
1672//--------------------------------------------------------------------------------------------------------
1673//--------------------------------------------------------------------------------------------------------
1674void Analysis::Help()
1675{
1676   
1677    cout << "This is the Help of the PlotScurveAll-diffFit" <<endl;
1678    cout << " Version 2.0 from September 27 2013 " <<endl;
1679    cout << " Main author Hiroko Miyamoto " << endl;
1680    cout << " Adapted for modularity and compilable program by Sylvie Dagoret " << endl;
1681    cout << " ==> Please from shell execute directly the compiled program :  " <<endl;
1682    cout << " ==> ../../src/PlotScurvesAllMain  " << endl;
1683   
1684}
1685//--------------------------------------------------------------------------------------------------------
1686
1687//--------------------------------------------------------------------------------------------------------
1688//--------------------------------------------------------------------------------------------------------
1689void Analysis::Menu()
1690{
1691    cout << "This is the Detailed menu information for running the PlotScurveAll-diffFit" <<endl;
1692   
1693    cout << " void PlotScurve( " << endl;
1694    cout << " \t " <<  " Int_t    kDrawCh=31, " << endl;
1695    cout << " \t " <<  " TString  fname=kFname, " << endl;
1696    cout << " \t " <<  " TString  fname_ped=kFnamePed, //ASIC_D ped_file " << endl;
1697    cout << " \t " <<  " Double_t dac_half_spe_fix = -1," << endl;
1698    cout << " \t " <<  " TString  fname_gain=\"gain-org.txt\"," << endl;
1699    cout << " \t " <<  " Int_t    fit_min=kFitMin," << endl;
1700    cout << " \t " <<  " Int_t    fit_max=kFitMax," << endl;
1701    cout << " \t " <<  " Int_t    fit_min_scurve=kFitMinS," << endl;
1702    cout << " \t " <<  " Int_t    fit_max_scurve=kFitMaxS," << endl;
1703    cout << " \t " <<  " Bool_t   kBoard=1, //0:ASIC_C, 1:ASIC_D " << endl;
1704    cout << " \t " <<  " Int_t    kChRef=kRefCh,//ASIC_C " << endl;
1705    cout << " \t " <<  " Bool_t   checkAve=0 //0:use ref pix, 1:check average " << endl;
1706    cout << " ) " << endl;             
1707    cout << "Set the parameters on the top of the program depending on your analysis:  " << endl;
1708    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;
1709    cout << " : " << endl;
1710    cout << " - kDrawCh: channel just you want to check on display. doens't effect the analysis. " << endl;
1711   
1712    cout << " - fname, fname_ped (kFname, kFnamePed): \"YOUR DATA FILE\" without .data"  << endl;
1713    cout << " See examples in the header file. " << endl;
1714    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;
1715    cout << " - fit_min, fit_max (kFitMin, kFitMax): leave them as default, it is not currently used for the analysis."  << endl;
1716    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;
1717    cout << " - kBoard: don't forget to change 0(C) or 1(D) depending on which ASIC you're analysing "  << endl;
1718    cout << " - kChRef (kRefCh): a pixel whose gain remains 64 thorough the gain adjustment. " << endl;
1719    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;
1720    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;
1721   
1722    cout << " Actual configuration " << endl;
1723    cout << " \t " <<  "TString kFname = " << kFname << endl;
1724    cout << " \t " <<  "TString  kFnamePed = " << kFnamePed  << endl;
1725    cout << " \t " <<  "Int_t kFitMin = " << kFitMin << endl; 
1726    cout << " \t " <<  "Int_t kFitMax = " << kFitMax << endl;
1727   
1728}
1729
1730//-------------------------------------------------------------------------------------------------
1731// Configure (x)emacs for this file ...
1732// Local Variables:
1733// mode: c++
1734// compile-command: "(make;) "
1735// End:
Note: See TracBrowser for help on using the repository browser.