Context Navigation

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

Last change on this file since 198 was 198, checked in by dagoret, 11 years ago
does not succed to remove the error in the interpreter with defines but keep the statements in a comment for later
File size: 33.3 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
36	using namespace std;
37	Int_t kXmin = 0;
38
39	//--------------------------------------------------------------------------------------------------------
40	void 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	//--------------------------------------------------------------------------------------------------------
54	void 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	//-------------------------------------------------------------------
97	void 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
138	void 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	ofstream out(kFnameOut.Data());
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 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
358	hist_scurve[ch] = TH1D(name_hist[ch],name_hist[ch],kNdata-1,
359	dac_min,dac_max);
360	hist_scurve_ped[ch] = TH1D(name_hist_ped[ch],name_hist_ped[ch],
361	kNdataPed-1,dac_min_ped,dac_max_ped);
362	hist_diff[ch] = TH1D(name_hist_diff[ch],name_hist_diff[ch],
363	kNdata-2,dac_min,dac_max);
364	hist_diff_fit[ch] = TH1D(name_hist_diff_fit[ch],name_hist_diff_fit[ch],
365	kNdataFit,fit_min_scurve,fit_max_scurve);
366
367	// read the data for that channel for normal data
368	in.open(input[ch].Data());
369	Double_t data[kNdata], dac[kNdata], dac_fit[kNdataFit], data_fit[kNdataFit];
370	while (in >> da >> dat ){
371	hist_scurve[ch].Fill(da,dat); // fill the histograms of the s-curve
372	}
373	in.close();
374	in.clear();
375
376	// read the data for pedestal data
377	in_ped.open(input_ped[ch].Data());
378	while (in_ped >> da >> dat ){
379	hist_scurve_ped[ch].Fill(da,dat); // fill the histogram for pedestal data
380	}
381	in_ped.close();
382	in_ped.clear();
383
384	//dac_ped[ch] = hist_scurve_ped[ch].GetMaximumBin()+dac_min_ped; //absolute value
385	dac_ped[ch] = hist_scurve_ped[ch].GetBinCenter(hist_scurve_ped[ch].GetMaximumBin()); //absolute value
386	dac_ped_sum += dac_ped[ch]; //absolute value
387	//cerr << "dac_ped[" << ch << "] dac_ped_sum " << dac_ped[ch]
388	//<< " " << dac_ped_sum << endl;
389
390
391	hist_scurve[ch].Smooth(3); // Smooth the S curve
392	Int_t binx=0;
393	//hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
394	// dac_max-dac_min+1,kHalfMax/2.);
395
396	// ??????
397	hist_scurve[ch].GetBinWithContent(kHalfMax,binx,kXmin-dac_min+1,
398	dac_max-dac_min+1,2000);
399	dac_half[ch] = hist_scurve[ch].GetBinCenter(binx); //absolute value
400	if(dac_half[ch]<kXmin)dac_half[ch]=-1;
401
402	// put the histograms of the smoothed S-Curves into two arrays (dac, data)
403	//-------------------------------------------------------------------------
404	Int_t i=0;
405	for(i=0;i<kNdata;i++){
406	dac[i] = hist_scurve[ch].GetBinCenter(i+1);
407	data[i] = hist_scurve[ch].GetBinContent(i+1);
408	}
409
410	/* compute derive diff on S-Curve */
411	Double_t diff[kNdata-1], dac_diff[kNdata-1];
412	Int_t trigger_diff=0;
413	for(i=0;i<kNdata-1;i++){
414	diff[i] = (data[i]-data[i+1])/step_dac;
415	dac_diff[i] = dac[i]+step_dac/2;
416	hist_diff[ch].Fill(dac_diff[i],diff[i]); // fills differentiated histograms of smoothed s-curves which is the single photoelectron spectrum
417	if(dac_diff[i]>=fit_min){
418	diff_sum += diff[i]*dac_diff[i];
419	//cerr << "diff[" << i << "], dac_diff[" << i << "] "
420	//<< diff[i] << " " << dac_diff[i] << endl;
421	trigger_diff+=diff[i];
422	}
423	}
424	//mean_diff[ch] = diff_sum/trigger_diff;
425	hist_diff[ch].GetXaxis()->SetRangeUser(fit_min,fit_max);
426	mean_diff[ch] = hist_diff[ch].GetMean();
427	//cerr << "trigger_diff, mean_diff[" << ch << "] "
428	//<< trigger_diff << " " << mean_diff[ch] << endl;
429
430
431	/***************** GRAPHICS ***************/
432
433
434	// Now start the analysis which consist in finding the peak in the differentiated spectrum
435	c1->cd(ch+1);
436	/* fit scurves */
437	cerr << "ch " << ch;
438
439	// SDC 24/09/2013 :: Need to declare outside this variable in a compiled program
440	TF1 *fit_pol ;
441
442
443	// Fit a polynomial curve in a given range on the s-curve directly. This is a kind of DEBUG procedure.
444
445	//if(ch==43 \|\| ch==47){ // SDC: 25/09/2013 Why a special channel : the fit pol4 is performed for another DAC range
446	/* if(ch==170){
447	fit_pol = new TF1("fit_pol","pol4",150,fit_max_scurve);
448	hist_scurve[ch].Fit(fit_pol,"","",150,fit_max_scurve);
449	}else{ */
450	fit_pol = new TF1("fit_pol","pol4",fit_min_scurve,fit_max_scurve);
451	hist_scurve[ch].Fit(fit_pol,"","",fit_min_scurve,fit_max_scurve);
452	// }
453
454
455
456	Double_t fit_par_pol4[5]; // container for fit parameter results
457	fit_pol->GetParameters(fit_par_pol4);
458	sfit_chiS[ch] = fit_pol->GetChisquare(); // chi squared
459	TF1 *func = hist_scurve[ch].GetFunction("fit_pol");
460	Double_t data_diff_fit;
461	for(i=0;i<kNdataFit;i++){
462	//if(ch==43 \|\| ch==47){
463	if(ch==170){
464	dac_fit[i] = hist_scurve[ch].GetBinCenter(i)+150-dac_min;
465	data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+150-dac_min); // get the value of the fitted function
466	}else{
467	dac_fit[i] = hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min;
468	data_fit[i] = func->Eval(hist_scurve[ch].GetBinCenter(i)+fit_min_scurve-dac_min); // get the value of the fitted function
469	//cerr << "dac_fit[" << i << "], data_fit[" << i << "] "
470	//<< dac_fit[i] << " " << data_fit[i] << endl;
471	}
472	}
473
474	// derive the fitted function on the S-curve to obtain the SPE
475	Double_t dac_diff_fit;
476	for(i=0;i<kNdataFit-1;i++){
477	data_diff_fit = (data_fit[i]-data_fit[i+1])/step_dac;
478	hist_diff_fit[ch].Fill(dac_fit[i],data_diff_fit); // Fill differentiated fitted S-curve
479	}
480
481	TLine *lhalf_max = new TLine(dac[0],kHalfMax,dac[kNdata-1],kHalfMax);
482	TLine *ldac_half = new TLine(dac_half[ch],0,dac_half[ch],
483	(Double_t)kDrawMax);
484	hist_scurve[ch].SetMaximum(kDrawMax);
485	hist_scurve[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
486	hist_scurve[ch].Draw();
487	lhalf_max->Draw("same");
488	ldac_half->Draw("same");
489	if(dac_half[ch]>0){
490	dac_half_sum += dac_half[ch]-dac_ped[ch]; //relative absolute value
491	trigger++;
492	}
493
494	/* Differential plot from scurve fitting function */
495	c5->cd(ch+1);
496	hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve,fit_max_scurve-2*step_dac);
497	Double_t dac_sfit_max=hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
498	//if(dac_sfit_max<=fit_min_scurve+4*step_dac){
499	if(dac_sfit_max<=fit_min_scurve+2*step_dac){
500	hist_diff_fit[ch].GetXaxis()->SetRangeUser(fit_min_scurve+kFitSadj,fit_max_scurve-2*step_dac);
501	dac_sfit[ch] = hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
502	//if(dac_sfit[ch]<=fit_min_scurve+kFitSadj+4*step_dac)dac_sfit[ch]=-1;
503	if(dac_sfit[ch]<fit_min_scurve+kFitSadj+step_dac){
504	dac_sfit[ch]=-1;
505	//dac_sfit[ch]=dac_half[ch];
506	flag_fit[ch] = -1;
507	cerr << "flag_fit[" << ch << "] " << flag_fit[ch] << endl;
508	}
509	}else {
510	dac_sfit[ch] = hist_diff_fit[ch].GetBinCenter(hist_diff_fit[ch].GetMaximumBin());
511	}
512
513	if(dac_sfit[ch]>fit_min_scurve+5){
514	dac_sfit_sum += dac_sfit[ch]-dac_ped[ch];
515	trigger_sfit++;
516	}else{
517	//dac_sfit[ch] = fit_min_scurve+10;
518	dac_sfit[ch] = -1;
519	}
520	//cerr << "dac_sfit[" << ch << "] " << dac_sfit[ch] << endl;
521	TLine *ldac_sfit = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMaxFit);
522	TLine *ldac_sfit_copy = new TLine(dac_sfit[ch],0,dac_sfit[ch],kDrawMax);
523	ldac_sfit->SetLineColor(kMagenta);
524	ldac_sfit_copy->SetLineColor(kMagenta);
525	hist_diff_fit[ch].SetMaximum(kDrawMaxFit);
526	hist_diff_fit[ch].SetMinimum(0);
527	hist_diff_fit[ch].Draw();
528	ldac_sfit->Draw("same");
529
530	/** differential plot **/
531	c2->cd(ch+1);
532	//hist_diff[ch].Rebin(8);
533	hist_diff[ch].Rebin(3);
534	fit_sdiff = new TF1("fit_sdiff","gaus",dac_min,dac_max);
535	//fit_sdiff->SetLineWidth(1);
536	//if(ch==8 \|\| ch==18 \|\| ch==44 \|\| ch==56)
537	mean_sum += mean_diff[ch]-dac_ped[ch];
538	trigger_mean++;
539
540	// fit the derivated polynom fit
541	Double_t fit_par[3];
542	//hist_diff[ch].Fit(fit_sdiff,"RQ0","",fit_min,fit_max);
543	hist_diff[ch].Fit(fit_sdiff,"","",fit_min,fit_max);
544	hist_diff[ch].GetXaxis()->SetRangeUser(kDrawDacMin,kDrawDacMax);
545	hist_diff[ch].SetMinimum(0);
546	hist_diff[ch].SetMaximum(kDrawMaxDiff);
547	//hist_diff[ch].Rebin();
548	hist_diff[ch].Draw();
549	fit_sdiff->GetParameters(fit_par); // get the parameters for this fit
550	dac_spe[ch] = fit_par[1];
551	if(dac_spe[ch]>fit_min){
552	dac_spe_sum += dac_spe[ch]-dac_ped[ch];
553	//cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
554	trigger_fit++;
555	}else{
556	dac_spe[ch] = fit_min+10;
557	}
558	//cerr << "dac_spe[" << ch << "] " << dac_spe[ch] << endl;
559	TLine *ldac_spe = new TLine(dac_spe[ch],0,dac_spe[ch],
560	(Double_t)kDrawMaxDiff);
561	ldac_spe->SetLineColor(kRed);
562	TLine *lmean_diff = new TLine(mean_diff[ch],0,mean_diff[ch],
563	(Double_t)kDrawMaxDiff);
564	lmean_diff->SetLineColor(kCyan);
565	TLine *lfit_min = new TLine(fit_min,0,fit_min,kDrawMaxDiff);
566	lfit_min->SetLineColor(kGray);
567
568	ldac_spe->Draw();
569	lmean_diff->Draw("same");
570	lfit_min->Draw("same");
571
572	c3->cd(1);
573	if(ch==0){
574	hist_scurve[ch].SetMaximum(kDrawMax);
575	hist_scurve[ch].Draw();
576	}else{
577	/if(ch==8 \|\| ch==18 \|\| ch==44 \|\| ch==56)/
578	hist_scurve[ch].Draw("same");
579	}
580	if(ch==kDrawCh){
581	c3->cd(2);
582	hist_scurve[ch].SetMaximum(kDrawMax);
583	hist_scurve[kDrawCh].Draw();
584	lhalf_max->Draw("same");
585	ldac_half->Draw("same");
586	ldac_sfit_copy->Draw("same");
587	c3->cd(3);
588	hist_diff[kDrawCh].Draw();
589	ldac_spe->Draw();
590	lmean_diff->Draw("same");
591	lfit_min->Draw("same");
592	c3->cd(4);
593	hist_diff_fit[kDrawCh].Draw();
594	ldac_sfit->Draw("same");
595	}
596	/*
597	if(ch<32){
598	c6->cd(ch/4+1);
599	hist_scurve[ch].SetMaximum(kDrawMax);
600	if(ch%4==0){
601	hist_scurve[ch].Draw();
602	}else{
603	hist_scurve[ch].Draw("same");
604	}
605	ldac_sfit_copy->Draw("same");
606	}else{
607	c7->cd(ch/4-7);
608	hist_scurve[ch].SetMaximum(kDrawMax);
609	if(ch%4==0){
610	hist_scurve[ch].Draw();
611	}else{
612	hist_scurve[ch].Draw("same");
613	}
614	ldac_sfit_copy->Draw("same");
615	}
616	*/
617
618	c4->cd();
619	if(ch==0){
620	hist_scurve_ped[ch].SetMaximum(kDrawMaxPed);
621	hist_scurve_ped[ch].Draw("l");
622	}else{
623	hist_scurve_ped[ch].Draw("lsame");
624	}
625
626	}// end of second loop on all the channels
627
628	// SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
629
630	Float_t dac_ped_ave ; //absolute value
631	Float_t dac_half_ave ;
632	Float_t dac_spe_ave ;
633	Float_t dac_sfit_ave;
634	Float_t mean_diff_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
635
636	if(!checkAve){
637	dac_ped_ave = dac_ped[kChRef]; //absolute value
638	dac_half_ave = dac_half[kChRef];
639	dac_spe_ave = dac_spe[kChRef];
640	dac_sfit_ave = dac_sfit[kChRef];
641	//Float_t dac_sfit_ave = 252.5; //dac_sfit of ch31 of ASIC_C with the 180 for lower gain pixels
642	mean_diff_ave = mean_diff[kChRef];
643	}else{
644	dac_ped_ave = dac_ped_sum/(Float_t)kNch; //absolute value
645	dac_half_ave = dac_half_sum/(Float_t)trigger+dac_ped_ave;
646	dac_spe_ave = dac_spe_sum/(Float_t)trigger_fit+dac_ped_ave;
647	mean_diff_ave = mean_sum/(Float_t)trigger_mean+dac_ped_ave;
648	dac_sfit_ave = dac_sfit_sum/(Float_t)trigger_sfit+dac_ped_ave;
649	}
650
651
652	// SDC 24/09/2013 :: In a compiled program, I have to declare the following variables outside the parenthesis
653
654	Double_t dac_half_spe;
655
656	//trigger_mean=4;
657	if(dac_half_spe_fix!=-1) {
658	dac_half_spe = dac_half_spe_fix;
659	} else {
660	//Double_t dac_half_spe = dac_half_ave;
661	dac_half_spe = dac_spe_ave;
662	}
663
664	TLine *ldac_ped_ave = new TLine(dac_ped_ave,0,dac_ped_ave,
665	(Double_t)kDrawMax);
666	TLine *ldac_ped_ave_copy = new TLine(dac_ped_ave,0,dac_ped_ave,
667	(Double_t)kDrawMaxDiff);
668	TLine *ldac_half_ave = new TLine(dac_half_ave,0,dac_half_ave,
669	(Double_t)kDrawMax);
670	TLine *ldac_half_spe = new TLine(dac_half_spe,0,dac_half_spe,
671	(Double_t)kDrawMax);
672	TLine *lmean_diff_ave = new TLine(mean_diff_ave,0,mean_diff_ave,
673	(Double_t)kDrawMax);
674	TLine *ldac_sfit_ave = new TLine(dac_sfit_ave,0,dac_sfit_ave,
675	(Double_t)kDrawMax);
676	ldac_ped_ave->SetLineColor(kOrange);
677	ldac_half_ave->SetLineColor(kRed);
678	lmean_diff_ave->SetLineColor(kCyan);
679	//ldac_half_spe->SetLineColor(kGreen);
680	ldac_sfit_ave->SetLineColor(kMagenta);
681
682	c3->cd(1);
683	ldac_ped_ave->Draw("same");
684	ldac_half_ave->Draw("same");
685	// lhalf_max->Draw("same"); SDC (24/09/2013) :: Not declared
686	ldac_half_spe->Draw("same");
687	lmean_diff_ave->Draw("same");
688	ldac_sfit_ave->Draw("same");
689
690	c3->cd(3);
691	ldac_ped_ave_copy->SetLineColor(kOrange);
692	ldac_ped_ave_copy->Draw("same");
693
694	Int_t gain_asic_diff, gain_asic_dacHalf, gain_asic_diff_mean, gain_asic_sfit;
695	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"
696	<< endl;
697
698
699	// loop 3 on all channels to make final computations on the gain
700	for(ch=0;ch<kNch;ch++){
701	/*
702	if(ch==8 \|\| ch==44 \|\| ch==56 \|\| ch==18){
703	gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
704	(mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
705	}else{
706	gain_asic_diff_mean=0;
707	}
708	*/
709	if(dac_half[ch]<0){
710	gain_asic_dacHalf = 0;
711	gain_asic_diff = 0;
712	gain_asic_diff_mean = 0;
713	gain_asic_sfit = 0;
714	}else{
715	gain_asic_dacHalf =(Int_t)((dac_half_ave-dac_ped_ave)/
716	(dac_half[ch]-dac_ped[ch])*gain_org[ch]);
717	gain_asic_diff =(Int_t)((dac_half_spe-dac_ped_ave)/
718	(dac_spe[ch]-dac_ped[ch])*gain_org[ch]);
719	gain_asic_diff_mean=(Int_t)((mean_diff_ave-dac_ped_ave)/
720	(mean_diff[ch]-dac_ped[ch])*gain_org[ch]);
721	gain_asic_sfit =(Int_t)((dac_sfit_ave-dac_ped_ave)/
722	(dac_sfit[ch]-dac_ped[ch])*gain_org[ch]);
723	//cerr << ch << " " << gain_org[ch] << endl;
724	}
725
726
727
728	out << "OUT_TABLE "
729	<< ch << " "
730	<< gain_asic_dacHalf << " "
731	<< gain_asic_diff << " "
732	<< gain_asic_diff_mean << " "
733	<< gain_asic_sfit << " "
734	<< dac_half[ch] << " "
735	<< dac_ped[ch] << " "
736	<< (Int_t)dac_spe[ch] << " "
737	<< dac_half_ave << " "
738	<< dac_half_spe << " "
739	<< dac_ped_ave << " "
740	<< dac_sfit_ave << " "
741	<< sfit_chiS[ch] << " "
742	<< ch << " "
743	<< flag_fit[ch]
744	<< endl;
745	std::cout << "channel " << ch << " old gain " << gain_org[ch] << " new gain " << gain_asic_sfit << std::endl;
746	}
747
748
749
750	/*
751	for(ch=kNch-1;ch>=0;ch--){
752	cout << "OUT_TABLE_INV "
753	<< (Int_t)(dac_half_ave/(dac_half[ch]-dac_ped)*64)
754	<< endl;
755	}
756	*/
757
758
759	}
760
761	//--------------------------------------------------------------------------------------------------------
762
763	//--------------------------------------------------------------------------------------------------------
764	// Analysis Constructor
765	//--------------------------------------------------------------------------------------------------------
766	Analysis::Analysis():ch(0),dac_half_sum(0),dac_spe_sum(0),dac_sfit_sum(0),diff_sum(0),mean_sum(0),maxBinSum(0)
767	{
768
769
770
771	/********************* SIGNAL RUN ******************/
772	input = new TString[kNch];
773	input_ped = new TString[kNch];
774
775	// Init sums
776	dac_half_sum = 0;
777	dac_ped_sum = 0;
778	dac_spe_sum = 0;
779	dac_sfit_sum = 0;
780	diff_sum = 0;
781	mean_sum = 0;
782	maxBinSum = 0;
783
784
785	}
786	//--------------------------------------------------------------------------------------------------------
787
788
789	//--------------------------------------------------------------------------------------------------------
790	// Analysis destructor
791	//--------------------------------------------------------------------------------------------------------
792	Analysis::~Analysis()
793	{
794
795	if(input!=0) delete input;
796	if(input_ped!=0) delete input_ped;
797
798
799	}
800	//--------------------------------------------------------------------------------------------------------
801
802
803	//--------------------------------------------------------------------------------------------------------
804	// Main public function of Analysis
805	//--------------------------------------------------------------------------------------------------------
806	void Analysis::PlotScurve(
807	Int_t kDrawCh=31,
808	TString fname=kFname,
809	TString fname_ped=kFnamePed, //ASIC_D ped_file
810	Double_t dac_half_spe_fix = -1,
811	TString fname_gain="gain-org.txt",
812	Int_t fit_min=kFitMin,
813	Int_t fit_max=kFitMax,
814	Int_t fit_min_scurve=kFitMinS,
815	Int_t fit_max_scurve=kFitMaxS,
816	Bool_t kBoard=0, //0:ASIC_C, 1:ASIC_D
817	Int_t kChRef=kRefCh,//ASIC_C
818	Bool_t checkAve=0 //0:use ref pix, 1:check average*/
819	)
820	{
821
822	// reset all variables
823	Init();
824
825	//OpenGainFile(fname_gain);
826
827	FillHistoName();
828
829	ReadGainFile(fname_gain,kBoard);
830	}
831	//--------------------------------------------------------------------------------------------------------
832
833
834	//--------------------------------------------------------------------------------------------------------
835	// Initialize the program by reseting counters or variables and create objects
836	//--------------------------------------------------------------------------------------------------------
837	void Analysis::Init()
838	{
839
840	if(input!=0) delete input;
841	if(input_ped!=0) delete input_ped;
842
843	// recreate the name
844	input = new TString[kNch];
845	input_ped = new TString[kNch];
846
847	// Init sums
848	dac_half_sum = 0;
849	dac_ped_sum = 0;
850	dac_spe_sum = 0;
851	dac_sfit_sum = 0;
852	diff_sum = 0;
853	mean_sum = 0;
854	maxBinSum = 0;
855
856
857	// create canvas
858
859	c1 = new TCanvas(kCanvasName1,kCanvasTitle1,200,20,1100,600); //smoothed S-curves
860	c1->Divide(11,6);
861	c2 = new TCanvas(kCanvasName2,kCanvasTitle2,300,50,1100,600); //differentiated histograms of smoothed s-curves
862	c2->Divide(11,6);
863	c5 = new TCanvas(kCanvasName5,kCanvasTitle5,400,70,1100,600);
864	c5->Divide(11,6);
865	c3 = new TCanvas(kCanvasName3,kCanvasTitle3,500,100,1200,300);
866	c3->Divide(4,1);
867	c4 = new TCanvas(kCanvasName4,kCanvasTitle4,900,400,300,300);
868
869
870	// create histos and function
871	hist_scurve = new TH1D[kNch]; // smoothed S-curves
872	hist_scurve_ped = new TH1D[kNch]; // S-curves for pedestal runs
873	hist_diff = new TH1D[kNch]; // differentiated histograms of smoothed s-curves
874	hist_diff_fit = new TH1D[kNch]; // differentiated histograms of fitted s-curves
875	fit_sdiff = new TF1[kNch]; // function fitted on differentiated histograms of fitted s-curves
876
877
878	for(ch=0;ch<kNch;ch++){
879	dac_half[ch]=0;
880	dac_ped[ch]=0;
881	dac_spe[ch]=0;
882	dac_sfit[ch]=0;
883	gain_org[ch]=0;
884	mean_diff[ch]=0;
885	sfit_chiS[ch]=0;
886	flag_fit[ch]=0;
887	}
888
889
890
891	}
892	//--------------------------------------------------------------------------------------------------------
893
894	//--------------------------------------------------------------------------------------------------------
895	// Initialize the program by reseting counters or variables and create objects
896	//--------------------------------------------------------------------------------------------------------
897	void Analysis::OpenGainFile(const TString fname_gain)
898	{
899	TString input_gain(fname_gain);
900	ifstream in_gain(input_gain.Data());
901	if(!in_gain.good()){
902	cerr << "File " << input_gain << " open failed!" << endl;
903	exit(-1);
904	}
905	in_gain.close();
906	in_gain.clear();
907
908	}
909	//--------------------------------------------------------------------------------------------------------
910
911	//--------------------------------------------------------------------------------------------------------
912	// Fill Histogram name
913	//--------------------------------------------------------------------------------------------------------
914	void Analysis::FillHistoName()
915	{
916
917	// first loop on channels to fill histonames
918	//----------------------
919	for(ch=0;ch<kNch;ch++){
920
921	name_hist[ch] ="hist";
922	name_hist[ch] +=ch;
923	name_hist_ped[ch] ="hist_ped";
924	name_hist_ped[ch] +=ch;
925	name_hist_diff[ch] ="hist_diff";
926	name_hist_diff[ch] +=ch;
927	name_hist_diff_fit[ch] ="hist_diff_fit";
928	name_hist_diff_fit[ch] +=ch;
929	name_fit_diff[ch] ="fit_diff";
930	name_fit_diff[ch] +=ch;
931	}
932
933	}
934	//--------------------------------------------------------------------------------------------------------
935
936
937	//--------------------------------------------------------------------------------------------------------
938	// Read the gain
939	//--------------------------------------------------------------------------------------------------------
940
941	void Analysis::ReadGainFile(const TString fname_gain,Bool_t kBoard)
942	{
943	TString input_gain(fname_gain);
944	ifstream in_gain(input_gain.Data());
945	if(!in_gain.good()){
946	cerr << "File " << input_gain << " open failed!" << endl;
947	exit(-1);
948	}
949	// read gain
950	Double_t da;
951	Int_t count=0;
952	ch=0;
953	while(in_gain >> da){
954	if(!kBoard && (count>=128 && count<192)){
955	gain_org[ch]=da;
956	cerr << "!kBoard, count, gain_org[" << ch << "] "
957	<< !kBoard << " " << count << " " << gain_org[ch] << endl;
958	ch++;
959	}
960	if(kBoard && (count>=192 && count<256)){
961	gain_org[ch]=da;
962	cerr << "kBoard, count, gain_org[" << ch << "] "
963	<< kBoard << " " << count << " " << gain_org[ch] << endl;
964	ch++;
965	}
966	count++;
967	}// end if first loop on all channels
968	in_gain.close();
969	in_gain.clear();
970
971
972	}
973	//--------------------------------------------------------------------------------------------------------

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

Download in other formats: