Context Navigation

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

Last change on this file since 197 was 197, checked in by dagoret, 11 years ago
bug due to a comment left corrected
File size: 33.2 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
21
22	// for format of EC_ASIC Scurves_all_ASIC data
23	// ASIC A,B,C,.... : i=0,1,2,... ch 0,1,...,63 : j=0,1,...,63
24	// % cat scurves_PC_inj_run5.lvm \|awk '{i=2;j=37;print($1,$(i*64+j+1+1))}'
25	// > scurves_PC_inj_run5-C-ch37.data
26	// => % cd $datadir
27	// => % sh $srcdir/makeDataFile.sh with assigining proper fname
28	// PlotScurvesAll-diffFit-v6.6.C: dac_sfit[ch]=-1; in the case of fitting failure.
29
30	using namespace std;
31	Int_t kXmin = 0;
32
33	//--------------------------------------------------------------------------------------------------------
34	void Help()
35	{
36
37	cout << "This is the Help of the PlotScurveAll-diffFit" <<endl;
38	cout << " Version 1.0 frol September 24 2013 " <<endl;
39	cout << " Main author Hiroko Miyamoto " << endl;
40	cout << " ==> Please in ROOT interpreter, issue the command .L PlotScurvesAll-diffFit.cc " <<endl;
41	cout << " ==> Then play with PlotScurve() function " << endl;
42	cout << " ==> If you need more information, please type Menu()" << endl;
43
44	}
45	//--------------------------------------------------------------------------------------------------------
46
47	//--------------------------------------------------------------------------------------------------------
48	void Menu()
49	{
50	cout << "This is the Detailed menu information for running the PlotScurveAll-diffFit" <<endl;
51
52	cout << " void PlotScurve( " << endl;
53	cout << " \t " << " Int_t kDrawCh=31, " << endl;
54	cout << " \t " << " TString fname=kFname, " << endl;
55	cout << " \t " << " TString fname_ped=kFnamePed, //ASIC_D ped_file " << endl;
56	cout << " \t " << " Double_t dac_half_spe_fix = -1," << endl;
57	cout << " \t " << " TString fname_gain=\"gain-org.txt\"," << endl;
58	cout << " \t " << " Int_t fit_min=kFitMin," << endl;
59	cout << " \t " << " Int_t fit_max=kFitMax," << endl;
60	cout << " \t " << " Int_t fit_min_scurve=kFitMinS," << endl;
61	cout << " \t " << " Int_t fit_max_scurve=kFitMaxS," << endl;
62	cout << " \t " << " Bool_t kBoard=1, //0:ASIC_C, 1:ASIC_D " << endl;
63	cout << " \t " << " Int_t kChRef=kRefCh,//ASIC_C " << endl;
64	cout << " \t " << " Bool_t checkAve=0 //0:use ref pix, 1:check average " << endl;
65	cout << " ) " << endl;
66	cout << "Set the parameters on the top of the program depending on your analysis: " << endl;
67	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;
68	cout << " : " << endl;
69	cout << " - kDrawCh: channel just you want to check on display. doens't effect the analysis. " << endl;
70
71	cout << " - fname, fname_ped (kFname, kFnamePed): \"YOUR DATA FILE\" without .data" << endl;
72	cout << " See examples in the header file. " << endl;
73	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;
74	cout << " - fit_min, fit_max (kFitMin, kFitMax): leave them as default, it is not currently used for the analysis." << endl;
75	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;
76	cout << " - kBoard: don't forget to change 0(C) or 1(D) depending on which ASIC you're analysing " << endl;
77	cout << " - kChRef (kRefCh): a pixel whose gain remains 64 thorough the gain adjustment. " << endl;
78	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;
79	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;
80
81	cout << " Actual configuration " << endl;
82	cout << " \t " << "TString kFname = " << kFname << endl;
83	cout << " \t " << "TString kFnamePed = " << kFnamePed << endl;
84	cout << " \t " << "Int_t kFitMin = " << kFitMin << endl;
85	cout << " \t " << "Int_t kFitMax = " << kFitMax << endl;
86
87	}
88	//--------------------------------------------------------------------------------------------------------
89
90	//-------------------------------------------------------------------
91	void InitRoot()
92	{
93
94	gStyle->SetTitleFillColor(kWhite);
95	gStyle->SetTitleColor(kBlue,"xyz");
96	gStyle->SetTitleSize(0.03,"xyz");
97	gStyle->SetLabelSize(0.03,"xyz");
98	gStyle->SetLabelColor(kBlue,"xyz");
99	gStyle->SetAxisColor(kMagenta);
100	gStyle->SetPalette(1,0);
101	gStyle->SetOptFit(101);
102	gStyle->SetOptStat(0);
103
104	// gStyle->SetOptStat("emr");
105	gStyle->SetOptStat("");
106	}
107	//-------------------------------------------------------------------
108
109
110	//--------------------------------------------------------------------------------------------------------
111	/*
112
113	The declaration of this function is now put inside the h file with the default values to be used
114
115	void PlotScurve(
116	Int_t kDrawCh=31,
117	TString fname=kFname,
118	TString fname_ped=kFnamePed, //ASIC_D ped_file
119	Double_t dac_half_spe_fix = -1,
120	TString fname_gain="gain-org.txt",
121	Int_t fit_min=kFitMin,
122	Int_t fit_max=kFitMax,
123	Int_t fit_min_scurve=kFitMinS,
124	Int_t fit_max_scurve=kFitMaxS,
125	Bool_t kBoard=1, //0:ASIC_C, 1:ASIC_D
126	Int_t kChRef=kRefCh,//ASIC_C
127	Bool_t checkAve=0 //0:use ref pix, 1:check average
128	)
129	*/
130
131	// Definition of the function PlotScurve
132	void PlotScurve(
133	Int_t kDrawCh,
134	TString fname,
135	TString fname_ped, //ASIC_D ped_file
136	Double_t dac_half_spe_fix,
137	TString fname_gain,
138	Int_t fit_min,
139	Int_t fit_max,
140	Int_t fit_min_scurve,
141	Int_t fit_max_scurve,
142	Bool_t kBoard, //0:ASIC_C, 1:ASIC_D
143	Int_t kChRef,//ASIC_C
144	Bool_t checkAve //0:use ref pix, 1:check average
145	)
146
147
148	{
149	//NOTE
150	//Double_t dac_half_spe_fix = 177; //for EC_ASICb N101(asic N158).
151	//check DAC value corresponds to voltage to 160fC input
152	//(depends on DAC linearity)
153	//Double_t dac_half_spe_fix = -1 //set gain to DAC average (default)
154
155	/***************** Readout data file *******************/
156
157	/********************* SIGNAL RUN ******************/
158	TString *input = new TString[kNch];
159	TString *input_ped = new TString[kNch];
160	Int_t ch;
161	Float_t dac_half_sum = 0;
162	Float_t dac_ped_sum = 0;
163	Float_t dac_spe_sum = 0;
164	Float_t dac_sfit_sum = 0;
165	Float_t diff_sum = 0;
166	Float_t mean_sum = 0;
167	Float_t maxBinSum = 0;
168	Double_t dac_half[kNch];
169	Double_t dac_ped[kNch];
170	Double_t dac_spe[kNch];
171	Double_t dac_sfit[kNch];
172	Double_t mean_diff[kNch];
173	Double_t gain_org[kNch];
174	Double_t sfit_chiS[kNch];
175	Double_t flag_fit[kNch];
176
177
178	InitRoot();
179
180	TCanvas *c1 = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
181	c1->Divide(11,6);
182	TCanvas *c2 = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
183	c2->Divide(11,6);
184	TCanvas *c5 = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
185	c5->Divide(11,6);
186	TCanvas *c3 = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
187	c3->Divide(4,1);
188	/TCanvas c6 = new TCanvas("c6","c6",450,150,1000,500);
189	c6->Divide(4,2);
190	TCanvas *c7 = new TCanvas("c7","c7",470,200,1000,500);
191	c7->Divide(4,2);*/
192	TCanvas *c4 = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
193
194	TString name_hist[kNch];
195	TString name_hist_ped[kNch];
196	TString name_hist_diff[kNch];
197	TString name_hist_diff_fit[kNch];
198	TString name_fit_diff[kNch];
199	TH1D *hist_scurve = new TH1D[kNch]; // smoothed S-curves
200	TH1D *hist_scurve_ped = new TH1D[kNch]; // S-curves for pedestal runs
201	TH1D *hist_diff = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
202	TH1D *hist_diff_fit = new TH1D[kNch]; // differentiated histograms of fitted s-curves
203	TF1 *fit_sdiff = new TF1[kNch]; // function fitted on differentiated histograms of fitted s-curves
204
205
206	// Open the ASIC Gain file
207	//-------------------------
208
209	TString input_gain(fname_gain);
210	ifstream in_gain(input_gain.Data());
211	if(!in_gain.good()){
212	cerr << "File " << input_gain << " open failed!" << endl;
213	return;
214	}
215
216	ofstream out(kFnameOut.Data());
217
218	// first loop on channels to fill histonames
219	//----------------------
220	for(ch=0;ch<kNch;ch++){
221	dac_half[ch]=0;
222	dac_ped[ch]=0;
223	dac_spe[ch]=0;
224	dac_sfit[ch]=0;
225	gain_org[ch]=0;
226	mean_diff[ch]=0;
227	sfit_chiS[ch]=0;
228	flag_fit[ch]=0;
229	name_hist[ch] ="hist";
230	name_hist[ch] +=ch;
231	name_hist_ped[ch] ="hist_ped";
232	name_hist_ped[ch] +=ch;
233	name_hist_diff[ch] ="hist_diff";
234	name_hist_diff[ch] +=ch;
235	name_hist_diff_fit[ch] ="hist_diff_fit";
236	name_hist_diff_fit[ch] +=ch;
237	name_fit_diff[ch] ="fit_diff";
238	name_fit_diff[ch] +=ch;
239	}
240
241
242	// read gain
243	Double_t da;
244	Int_t count=0;
245	ch=0;
246	while(in_gain >> da){
247	if(!kBoard && (count>=128 && count<192)){
248	gain_org[ch]=da;
249	cerr << "!kBoard, count, gain_org[" << ch << "] "
250	<< !kBoard << " " << count << " " << gain_org[ch] << endl;
251	ch++;
252	}
253	if(kBoard && (count>=192 && count<256)){
254	gain_org[ch]=da;
255	cerr << "kBoard, count, gain_org[" << ch << "] "
256	<< kBoard << " " << count << " " << gain_org[ch] << endl;
257	ch++;
258	}
259	count++;
260	}// end if first loop on all channels
261	in_gain.close();
262	in_gain.clear();
263
264
265
266	// Open the S-Curve files for each channel and test if the file can be open
267
268	Int_t trigger=0, trigger_fit=0, trigger_mean=0, trigger_sfit=0;
269	// Big loop on all the channels : second loop
270	//--------------------------------------------
271	// build the input data file name for each channel fname-chXX.dat
272	for(ch=0;ch<kNch;ch++){
273	input[ch] += fname;
274	input[ch] += "-ch";
275	input[ch] += ch;
276	input[ch] += ".dat";
277	ifstream in(input[ch].Data());
278	if(!in.good()){
279	cerr << "File " << input[ch] << " open failed!" << endl;
280	return;
281	}
282
283	// build the input data file name for each channel fname_ped-chXX.dat
284	input_ped[ch] += fname_ped;
285	input_ped[ch] += "-ch";
286	input_ped[ch] += ch;
287	input_ped[ch] += ".dat";
288	ifstream in_ped(input_ped[ch].Data());
289	if(!in_ped.good()){
290	cerr << "File " << input_ped[ch] << " open failed!" << endl;
291	return;
292	}
293
294
295
296	//10DAC~10mV
297	// determine the DAC range in dac_min, dac_max, step_dac from file input_ped
298	// and count the number of events
299	Double_t da, dat, dac_min, dac_max, step_dac;
300	count = 0;
301	while (in >> da >> dat ){
302	if(count==0) dac_min = da;
303	if(count==1) step_dac = da - dac_min;
304	dac_max=da;
305	count++;
306	}
307	in.close();
308	in.clear();
309
310	/*
311	if(ch==43 \|\| ch==47){//for T1 ASIC_C
312	cerr << "CHECK " << ch << " " << fit_min_scurve << endl;
313	Int_t fit_min_scurve=140;
314	} else {
315	Int_t fit_min_scurve=kFitMinS;
316	}
317	*/
318
319	// correct the scale for DAC range
320	if(kXmin<dac_min)kXmin=dac_min;
321	const Int_t kNdata = count;
322	count = (fit_max_scurve-fit_min_scurve)/step_dac;
323	const Int_t kNdataFit = count;
324	if(ch==0){
325	cerr << "kNdata, dac_min, dac_max, step_dac: "
326	<< kNdata << " " << dac_min << " " << dac_max << " "
327	<< step_dac << endl;
328	}
329
330	Double_t dac_min_ped, dac_max_ped, step_dac_ped;
331	count = 0;
332	while (in_ped >> da >> dat ){
333	if(count==0) dac_min_ped = da;
334	if(count==1) step_dac_ped = da - dac_min_ped;
335	dac_max_ped=da;
336	count++;
337	}
338
339	in_ped.close();
340	in_ped.clear();
341
342
343	const Int_t kNdataPed = count;
344	if(ch==0){
345	cerr << "kNdataPed, dac_min_ped, dac_max_ped, step_dac_ped: "
346	<< kNdataPed << " " << dac_min_ped << " " << dac_max_ped
347	<< " " << step_dac_ped << endl;
348	}
349
350	// book histograms
351
352	hist_scurve[ch] = TH1D(name_hist[ch],name_hist[ch],kNdata-1,
353	dac_min,dac_max);
354	hist_scurve_ped[ch] = TH1D(name_hist_ped[ch],name_hist_ped[ch],
355	kNdataPed-1,dac_min_ped,dac_max_ped);
356	hist_diff[ch] = TH1D(name_hist_diff[ch],name_hist_diff[ch],
357	kNdata-2,dac_min,dac_max);
358	hist_diff_fit[ch] = TH1D(name_hist_diff_fit[ch],name_hist_diff_fit[ch],
359	kNdataFit,fit_min_scurve,fit_max_scurve);
360
361	// read the data for that channel for normal data
362	in.open(input[ch].Data());
363	Double_t data[kNdata], dac[kNdata], dac_fit[kNdataFit], data_fit[kNdataFit];
364	while (in >> da >> dat ){
365	hist_scurve[ch].Fill(da,dat); // fill the histograms of the s-curve
366	}
367	in.close();
368	in.clear();
369
370	// read the data for pedestal data
371	in_ped.open(input_ped[ch].Data());
372	while (in_ped >> da >> dat ){
373	hist_scurve_ped[ch].Fill(da,dat); // fill the histogram for pedestal data
374	}
375	in_ped.close();
376	in_ped.clear();
377
378	//dac_ped[ch] = hist_scurve_ped[ch].GetMaximumBin()+dac_min_ped; //absolute value
379	dac_ped[ch] = hist_scurve_ped[ch].GetBinCenter(hist_scurve_ped[ch].GetMaximumBin()); //absolute value
380	dac_ped_sum += dac_ped[ch]; //absolute value
381	//cerr << "dac_ped[" << ch << "] dac_ped_sum " << dac_ped[ch]
382	//<< " " << dac_ped_sum << endl;
383
384
385	hist_scurve[ch].Smooth(3); // Smooth the S curve
386	Int_t binx=0;
387	//hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
388	// dac_max-dac_min+1,kHalfMax/2.);
389
390	// ??????
391	hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
392	dac_max-dac_min+1,2000);
393	dac_half[ch] = hist_scurve[ch].GetBinCenter(binx); //absolute value
394	if(dac_half[ch]<kXmin)dac_half[ch]=-1;
395
396	// put the histograms of the smoothed S-Curves into two arrays (dac, data)
397	//-------------------------------------------------------------------------
398	Int_t i=0;
399	for(i=0;i<kNdata;i++){
400	dac[i] = hist_scurve[ch].GetBinCenter(i+1);
401	data[i] = hist_scurve[ch].GetBinContent(i+1);
402	}
403
404	/* compute derive diff on S-Curve */
405	Double_t diff[kNdata-1], dac_diff[kNdata-1];
406	Int_t trigger_diff=0;
407	for(i=0;i<kNdata-1;i++){
408	diff[i] = (data[i]-data[i+1])/step_dac;
409	dac_diff[i] = dac[i]+step_dac/2;
410	hist_diff[ch].Fill(dac_diff[i],diff[i]); // fills differentiated histograms of smoothed s-curves which is the single photoelectron spectrum
411	if(dac_diff[i]>=fit_min){
412	diff_sum += diff[i]*dac_diff[i];
413	//cerr << "diff[" << i << "], dac_diff[" << i << "] "
414	//<< diff[i] << " " << dac_diff[i] << endl;
415	trigger_diff+=diff[i];
416	}
417	}
418	//mean_diff[ch] = diff_sum/trigger_diff;
419	hist_diff[ch].GetXaxis()->SetRangeUser(fit_min,fit_max);
420	mean_diff[ch] = hist_diff[ch].GetMean();
421	//cerr << "trigger_diff, mean_diff[" << ch << "] "
422	//<< trigger_diff << " " << mean_diff[ch] << endl;
423
424
425	/***************** GRAPHICS ***************/
426
427
428	// Now start the analysis which consist in finding the peak in the differentiated spectrum
429	c1->cd(ch+1);
430	/* fit scurves */
431	cerr << "ch " << ch;
432
433	// SDC 24/09/2013 :: Need to declare outside this variable in a compiled program
434	TF1 *fit_pol ;
435
436
437	// Fit a polynomial curve in a given range on the s-curve directly. This is a kind of DEBUG procedure.
438
439	//if(ch==43 \|\| ch==47){ // SDC: 25/09/2013 Why a special channel : the fit pol4 is performed for another DAC range
440	/* if(ch==170){
441	fit_pol = new TF1("fit_pol","pol4",150,fit_max_scurve);
442	hist_scurve[ch].Fit(fit_pol,"","",150,fit_max_scurve);
443	}else{ */
444	fit_pol = new TF1("fit_pol","pol4",fit_min_scurve,fit_max_scurve);
445	hist_scurve[ch].Fit(fit_pol,"","",fit_min_scurve,fit_max_scurve);
446	// }
447
448
449
450	Double_t fit_par_pol4[5]; // container for fit parameter results
451	fit_pol->GetParameters(fit_par_pol4);
452	sfit_chiS[ch] = fit_pol->GetChisquare(); // chi squared
453	TF1 *func = hist_scurve[ch].GetFunction("fit_pol");
454	Double_t data_diff_fit;
455	for(i=0;i<kNdataFit;i++){
456	//if(ch==43 \|\| ch==47){
457	if(ch==170){
458	dac_fit[i] = hist_scurve[ch].GetBinCenter(i)+150-dac_min;
459	data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+150-dac_min); // get the value of the fitted function
460	}else{
461	dac_fit[i] = hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min;
462	data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min); // get the value of the fitted function
463	//cerr << "dac_fit[" << i << "], data_fit[" << i << "] "
464	//<< dac_fit[i] << " " << data_fit[i] << endl;
465	}
466	}
467
468	// derive the fitted function on the S-curve to obtain the SPE
469	Double_t dac_diff_fit;
470	for(i=0;i<kNdataFit-1;i++){
471	data_diff_fit = (data_fit[i]-data_fit[i+1])/step_dac;
472	hist_diff_fit[ch].Fill(dac_fit[i],data_diff_fit); // Fill differentiated fitted S-curve
473	}
474
475	TLine *lhalf_max = new TLine(dac[0],kHalfMax,dac[kNdata-1],kHalfMax);
476	TLine *ldac_half = new TLine(dac_half[ch],0,dac_half[ch],
477	(Double_t)kDrawMax);
478	hist_scurve[ch].SetMaximum(kDrawMax);
479	hist_scurve[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
480	hist_scurve[ch].Draw();
481	lhalf_max->Draw("same");
482	ldac_half->Draw("same");
483	if(dac_half[ch]>0){
484	dac_half_sum += dac_half[ch]-dac_ped[ch]; //relative absolute value
485	trigger++;
486	}
487
488	/* Differential plot from scurve fitting function */
489	c5->cd(ch+1);
490	hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve,fit_max_scurve-2*step_dac);
491	Double_t dac_sfit_max=hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
492	//if(dac_sfit_max<=fit_min_scurve+4*step_dac){
493	if(dac_sfit_max<=fit_min_scurve+2*step_dac){
494	hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve+kFitSadj,fit_max_scurve-2*step_dac);
495	dac_sfit[ch] = hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
496	//if(dac_sfit[ch]<=fit_min_scurve+kFitSadj+4*step_dac)dac_sfit[ch]=-1;
497	if(dac_sfit[ch]<fit_min_scurve+kFitSadj+step_dac){
498	dac_sfit[ch]=-1;
499	//dac_sfit[ch]=dac_half[ch];
500	flag_fit[ch] = -1;
501	cerr << "flag_fit[" << ch << "] " << flag_fit[ch] << endl;
502	}
503	}else {
504	dac_sfit[ch] = hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
505	}
506
507	if(dac_sfit[ch]>fit_min_scurve+5){
508	dac_sfit_sum += dac_sfit[ch]-dac_ped[ch];
509	trigger_sfit++;
510	}else{
511	//dac_sfit[ch] = fit_min_scurve+10;
512	dac_sfit[ch] = -1;
513	}
514	//cerr << "dac_sfit[" << ch << "] " << dac_sfit[ch] << endl;
515	TLine *ldac_sfit = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMaxFit);
516	TLine *ldac_sfit_copy = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMax);
517	ldac_sfit->SetLineColor(kMagenta);
518	ldac_sfit_copy->SetLineColor(kMagenta);
519	hist_diff_fit[ch].SetMaximum(kDrawMaxFit);
520	hist_diff_fit[ch].SetMinimum(0);
521	hist_diff_fit[ch].Draw();
522	ldac_sfit->Draw("same");
523
524	/** differential plot **/
525	c2->cd(ch+1);
526	//hist_diff[ch].Rebin(8);
527	hist_diff[ch].Rebin(3);
528	fit_sdiff = new TF1("fit_sdiff","gaus",dac_min,dac_max);
529	//fit_sdiff->SetLineWidth(1);
530	//if(ch==8 \|\| ch==18 \|\| ch==44 \|\| ch==56)
531	mean_sum += mean_diff[ch]-dac_ped[ch];
532	trigger_mean++;
533
534	// fit the derivated polynom fit
535	Double_t fit_par[3];
536	//hist_diff[ch].Fit(fit_sdiff,"RQ0","",fit_min,fit_max);
537	hist_diff[ch].Fit(fit_sdiff,"","",fit_min,fit_max);
538	hist_diff[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
539	hist_diff[ch].SetMinimum(0);
540	hist_diff[ch].SetMaximum(kDrawMaxDiff);
541	//hist_diff[ch].Rebin();
542	hist_diff[ch].Draw();
543	fit_sdiff->GetParameters(fit_par); // get the parameters for this fit
544	dac_spe[ch] = fit_par[1];
545	if(dac_spe[ch]>fit_min){
546	dac_spe_sum += dac_spe[ch]-dac_ped[ch];
547	//cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
548	trigger_fit++;
549	}else{
550	dac_spe[ch] = fit_min+10;
551	}
552	//cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
553	TLine *ldac_spe = new TLine(dac_spe[ch],0,dac_spe[ch],
554	(Double_t)kDrawMaxDiff);
555	ldac_spe->SetLineColor(kRed);
556	TLine *lmean_diff = new TLine(mean_diff[ch],0,mean_diff[ch],
557	(Double_t)kDrawMaxDiff);
558	lmean_diff->SetLineColor(kCyan);
559	TLine *lfit_min = new TLine(fit_min,0,fit_min,kDrawMaxDiff);
560	lfit_min->SetLineColor(kGray);
561
562	ldac_spe->Draw();
563	lmean_diff->Draw("same");
564	lfit_min->Draw("same");
565
566	c3->cd(1);
567	if(ch==0){
568	hist_scurve[ch].SetMaximum(kDrawMax);
569	hist_scurve[ch].Draw();
570	}else{
571	/if(ch==8 \|\| ch==18 \|\| ch==44 \|\| ch==56)/
572	hist_scurve[ch].Draw("same");
573	}
574	if(ch==kDrawCh){
575	c3->cd(2);
576	hist_scurve[ch].SetMaximum(kDrawMax);
577	hist_scurve[kDrawCh].Draw();
578	lhalf_max->Draw("same");
579	ldac_half->Draw("same");
580	ldac_sfit_copy->Draw("same");
581	c3->cd(3);
582	hist_diff[kDrawCh].Draw();
583	ldac_spe->Draw();
584	lmean_diff->Draw("same");
585	lfit_min->Draw("same");
586	c3->cd(4);
587	hist_diff_fit[kDrawCh].Draw();
588	ldac_sfit->Draw("same");
589	}
590	/*
591	if(ch<32){
592	c6->cd(ch/4+1);
593	hist_scurve[ch].SetMaximum(kDrawMax);
594	if(ch%4==0){
595	hist_scurve[ch].Draw();
596	}else{
597	hist_scurve[ch].Draw("same");
598	}
599	ldac_sfit_copy->Draw("same");
600	}else{
601	c7->cd(ch/4-7);
602	hist_scurve[ch].SetMaximum(kDrawMax);
603	if(ch%4==0){
604	hist_scurve[ch].Draw();
605	}else{
606	hist_scurve[ch].Draw("same");
607	}
608	ldac_sfit_copy->Draw("same");
609	}
610	*/
611
612	c4->cd();
613	if(ch==0){
614	hist_scurve_ped[ch].SetMaximum(kDrawMaxPed);
615	hist_scurve_ped[ch].Draw("l");
616	}else{
617	hist_scurve_ped[ch].Draw("lsame");
618	}
619
620	}// end of second loop on all the channels
621
622	// SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
623
624	Float_t dac_ped_ave ; //absolute value
625	Float_t dac_half_ave ;
626	Float_t dac_spe_ave ;
627	Float_t dac_sfit_ave;
628	Float_t mean_diff_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
629
630	if(!checkAve){
631	dac_ped_ave = dac_ped[kChRef]; //absolute value
632	dac_half_ave = dac_half[kChRef];
633	dac_spe_ave = dac_spe[kChRef];
634	dac_sfit_ave = dac_sfit[kChRef];
635	//Float_t dac_sfit_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
636	mean_diff_ave = mean_diff[kChRef];
637	}else{
638	dac_ped_ave = dac_ped_sum/(Float_t)kNch; //absolute value
639	dac_half_ave = dac_half_sum/(Float_t)trigger+dac_ped_ave;
640	dac_spe_ave = dac_spe_sum/(Float_t)trigger_fit+dac_ped_ave;
641	mean_diff_ave = mean_sum/(Float_t)trigger_mean+dac_ped_ave;
642	dac_sfit_ave = dac_sfit_sum/(Float_t)trigger_sfit+dac_ped_ave;
643	}
644
645
646	// SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
647
648	Double_t dac_half_spe;
649
650	//trigger_mean=4;
651	if(dac_half_spe_fix!=-1) {
652	dac_half_spe = dac_half_spe_fix;
653	} else {
654	//Double_t dac_half_spe = dac_half_ave;
655	dac_half_spe = dac_spe_ave;
656	}
657
658	TLine *ldac_ped_ave = new TLine(dac_ped_ave,0,dac_ped_ave,
659	(Double_t)kDrawMax);
660	TLine *ldac_ped_ave_copy = new TLine(dac_ped_ave,0,dac_ped_ave,
661	(Double_t)kDrawMaxDiff);
662	TLine *ldac_half_ave = new TLine(dac_half_ave,0,dac_half_ave,
663	(Double_t)kDrawMax);
664	TLine *ldac_half_spe = new TLine(dac_half_spe,0,dac_half_spe,
665	(Double_t)kDrawMax);
666	TLine *lmean_diff_ave = new TLine(mean_diff_ave,0,mean_diff_ave,
667	(Double_t)kDrawMax);
668	TLine *ldac_sfit_ave = new TLine(dac_sfit_ave,0,dac_sfit_ave,
669	(Double_t)kDrawMax);
670	ldac_ped_ave->SetLineColor(kOrange);
671	ldac_half_ave->SetLineColor(kRed);
672	lmean_diff_ave->SetLineColor(kCyan);
673	//ldac_half_spe->SetLineColor(kGreen);
674	ldac_sfit_ave->SetLineColor(kMagenta);
675
676	c3->cd(1);
677	ldac_ped_ave->Draw("same");
678	ldac_half_ave->Draw("same");
679	// lhalf_max->Draw("same"); SDC (24/09/2013) :: Not declared
680	ldac_half_spe->Draw("same");
681	lmean_diff_ave->Draw("same");
682	ldac_sfit_ave->Draw("same");
683
684	c3->cd(3);
685	ldac_ped_ave_copy->SetLineColor(kOrange);
686	ldac_ped_ave_copy->Draw("same");
687
688	Int_t gain_asic_diff, gain_asic_dacHalf, gain_asic_diff_mean, gain_asic_sfit;
689	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"
690	<< endl;
691
692
693	// loop 3 on all channels to make final computations on the gain
694	for(ch=0;ch<kNch;ch++){
695	/*
696	if(ch==8 \|\| ch==44 \|\| ch==56 \|\| ch==18){
697	gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
698	(mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
699	}else{
700	gain_asic_diff_mean=0;
701	}
702	*/
703	if(dac_half[ch]<0){
704	gain_asic_dacHalf = 0;
705	gain_asic_diff = 0;
706	gain_asic_diff_mean = 0;
707	gain_asic_sfit = 0;
708	}else{
709	gain_asic_dacHalf =(Int_t)((dac_half_ave-dac_ped_ave)/
710	(dac_half[ch]-dac_ped[ch])*gain_org[ch]);
711	gain_asic_diff =(Int_t)((dac_half_spe-dac_ped_ave)/
712	(dac_spe[ch]-dac_ped[ch])*gain_org[ch]);
713	gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
714	(mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
715	gain_asic_sfit =(Int_t)((dac_sfit_ave-dac_ped_ave)/
716	(dac_sfit[ch]-dac_ped[ch])*gain_org[ch]);
717	//cerr << ch << " " << gain_org[ch] << endl;
718	}
719
720
721
722	out << "OUT_TABLE "
723	<< ch << " "
724	<< gain_asic_dacHalf << " "
725	<< gain_asic_diff << " "
726	<< gain_asic_diff_mean << " "
727	<< gain_asic_sfit << " "
728	<< dac_half[ch] << " "
729	<< dac_ped[ch] << " "
730	<< (Int_t)dac_spe[ch] << " "
731	<< dac_half_ave << " "
732	<< dac_half_spe << " "
733	<< dac_ped_ave << " "
734	<< dac_sfit_ave << " "
735	<< sfit_chiS[ch] << " "
736	<< ch << " "
737	<< flag_fit[ch]
738	<< endl;
739	std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
740	}
741
742
743
744	/*
745	for(ch=kNch-1;ch>=0;ch--){
746	cout << "OUT_TABLE_INV "
747	<< (Int_t)(dac_half_ave/(dac_half[ch]-dac_ped)*64)
748	<< endl;
749	}
750	*/
751
752
753	}
754
755	//--------------------------------------------------------------------------------------------------------
756
757	//--------------------------------------------------------------------------------------------------------
758	// Analysis Constructor
759	//--------------------------------------------------------------------------------------------------------
760	Analysis::Analysis():ch(0),dac_half_sum(0),dac_spe_sum(0),dac_sfit_sum(0),diff_sum(0),mean_sum(0),maxBinSum(0)
761	{
762
763
764
765	/********************* SIGNAL RUN ******************/
766	input = new TString[kNch];
767	input_ped = new TString[kNch];
768
769	// Init sums
770	dac_half_sum = 0;
771	dac_ped_sum = 0;
772	dac_spe_sum = 0;
773	dac_sfit_sum = 0;
774	diff_sum = 0;
775	mean_sum = 0;
776	maxBinSum = 0;
777
778
779	}
780	//--------------------------------------------------------------------------------------------------------
781
782
783	//--------------------------------------------------------------------------------------------------------
784	// Analysis destructor
785	//--------------------------------------------------------------------------------------------------------
786	Analysis::~Analysis()
787	{
788
789	if(input!=0) delete input;
790	if(input_ped!=0) delete input_ped;
791
792
793	}
794	//--------------------------------------------------------------------------------------------------------
795
796
797	//--------------------------------------------------------------------------------------------------------
798	// Main public function of Analysis
799	//--------------------------------------------------------------------------------------------------------
800	void Analysis::PlotScurve(
801	Int_t kDrawCh=31,
802	TString fname=kFname,
803	TString fname_ped=kFnamePed, //ASIC_D ped_file
804	Double_t dac_half_spe_fix = -1,
805	TString fname_gain="gain-org.txt",
806	Int_t fit_min=kFitMin,
807	Int_t fit_max=kFitMax,
808	Int_t fit_min_scurve=kFitMinS,
809	Int_t fit_max_scurve=kFitMaxS,
810	Bool_t kBoard=0, //0:ASIC_C, 1:ASIC_D
811	Int_t kChRef=kRefCh,//ASIC_C
812	Bool_t checkAve=0 //0:use ref pix, 1:check average*/
813	)
814	{
815
816	// reset all variables
817	Init();
818
819	//OpenGainFile(fname_gain);
820
821	FillHistoName();
822
823	ReadGainFile(fname_gain,kBoard);
824	}
825	//--------------------------------------------------------------------------------------------------------
826
827
828	//--------------------------------------------------------------------------------------------------------
829	// Initialize the program by reseting counters or variables and create objects
830	//--------------------------------------------------------------------------------------------------------
831	void Analysis::Init()
832	{
833
834	if(input!=0) delete input;
835	if(input_ped!=0) delete input_ped;
836
837	// recreate the name
838	input = new TString[kNch];
839	input_ped = new TString[kNch];
840
841	// Init sums
842	dac_half_sum = 0;
843	dac_ped_sum = 0;
844	dac_spe_sum = 0;
845	dac_sfit_sum = 0;
846	diff_sum = 0;
847	mean_sum = 0;
848	maxBinSum = 0;
849
850
851	// create canvas
852
853	c1 = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
854	c1->Divide(11,6);
855	c2 = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
856	c2->Divide(11,6);
857	c5 = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
858	c5->Divide(11,6);
859	c3 = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
860	c3->Divide(4,1);
861	c4 = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
862
863
864	// create histos and function
865	hist_scurve = new TH1D[kNch]; // smoothed S-curves
866	hist_scurve_ped = new TH1D[kNch]; // S-curves for pedestal runs
867	hist_diff = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
868	hist_diff_fit = new TH1D[kNch]; // differentiated histograms of fitted s-curves
869	fit_sdiff = new TF1[kNch]; // function fitted on differentiated histograms of fitted s-curves
870
871
872	for(ch=0;ch<kNch;ch++){
873	dac_half[ch]=0;
874	dac_ped[ch]=0;
875	dac_spe[ch]=0;
876	dac_sfit[ch]=0;
877	gain_org[ch]=0;
878	mean_diff[ch]=0;
879	sfit_chiS[ch]=0;
880	flag_fit[ch]=0;
881	}
882
883
884
885	}
886	//--------------------------------------------------------------------------------------------------------
887
888	//--------------------------------------------------------------------------------------------------------
889	// Initialize the program by reseting counters or variables and create objects
890	//--------------------------------------------------------------------------------------------------------
891	void Analysis::OpenGainFile(const TString fname_gain)
892	{
893	TString input_gain(fname_gain);
894	ifstream in_gain(input_gain.Data());
895	if(!in_gain.good()){
896	cerr << "File " << input_gain << " open failed!" << endl;
897	exit(-1);
898	}
899	in_gain.close();
900	in_gain.clear();
901
902	}
903	//--------------------------------------------------------------------------------------------------------
904
905	//--------------------------------------------------------------------------------------------------------
906	// Fill Histogram name
907	//--------------------------------------------------------------------------------------------------------
908	void Analysis::FillHistoName()
909	{
910
911	// first loop on channels to fill histonames
912	//----------------------
913	for(ch=0;ch<kNch;ch++){
914
915	name_hist[ch] ="hist";
916	name_hist[ch] +=ch;
917	name_hist_ped[ch] ="hist_ped";
918	name_hist_ped[ch] +=ch;
919	name_hist_diff[ch] ="hist_diff";
920	name_hist_diff[ch] +=ch;
921	name_hist_diff_fit[ch] ="hist_diff_fit";
922	name_hist_diff_fit[ch] +=ch;
923	name_fit_diff[ch] ="fit_diff";
924	name_fit_diff[ch] +=ch;
925	}
926
927	}
928	//--------------------------------------------------------------------------------------------------------
929
930
931	//--------------------------------------------------------------------------------------------------------
932	// Read the gain
933	//--------------------------------------------------------------------------------------------------------
934
935	void Analysis::ReadGainFile(const TString fname_gain,Bool_t kBoard)
936	{
937	TString input_gain(fname_gain);
938	ifstream in_gain(input_gain.Data());
939	if(!in_gain.good()){
940	cerr << "File " << input_gain << " open failed!" << endl;
941	exit(-1);
942	}
943	// read gain
944	Double_t da;
945	Int_t count=0;
946	ch=0;
947	while(in_gain >> da){
948	if(!kBoard && (count>=128 && count<192)){
949	gain_org[ch]=da;
950	cerr << "!kBoard, count, gain_org[" << ch << "] "
951	<< !kBoard << " " << count << " " << gain_org[ch] << endl;
952	ch++;
953	}
954	if(kBoard && (count>=192 && count<256)){
955	gain_org[ch]=da;
956	cerr << "kBoard, count, gain_org[" << ch << "] "
957	<< kBoard << " " << count << " " << gain_org[ch] << endl;
958	ch++;
959	}
960	count++;
961	}// end if first loop on all channels
962	in_gain.close();
963	in_gain.clear();
964
965
966	}
967	//--------------------------------------------------------------------------------------------------------

Note: See TracBrowser for help on using the repository browser.

Download in other formats: