WebSVN – f9daq – Blame – /cvi/apps/l2d2_easyroc/plots2.cpp

Rev	Author	Line No.	Line
291	f9daq	1	#include "TROOT.h"
		2	#include "TFile.h"
		3	#include "TBenchmark.h"
		4	#include "TH1F.h"
		5	#include "TH2F.h"
		6	gROOT->Reset();
		7
		8	#include "TCanvas.h"
		9	#include "TStyle.h"
		10	#include "TPad.h"
		11	#include "TF1.h"
		12	#include "TGraph.h"
		13	#include "TSpectrum.h"
		14
		15	#include "DrootHelper.cpp"
		16
		17
		18	#define NCH 2
		19
		20	//-------------------------------------------------------------------------------------
		21
		22	void SetGS() {
		23	const UInt_t Number = 2;
		24	Double_t Red[Number] = {0.00, 1.00};
		25	Double_t Green[Number] = {0.00, 1.00};
		26	Double_t Blue[Number] = {0.00, 1.00};
		27	Double_t Stops[Number] = {0.00, 1.00};
		28
		29	Int_t nb=50;
		30	TColor::CreateGradientColorTable(Number, Stops, Red, Green, Blue, nb);
		31	}
		32
		33
		34	//-------------------------------------------------------------------------------------
		35	Double_t ConstantStep(Double_t x)
		36	{
		37	if(x<=0) return 0.;
		38	else return 1.;
		39	}
		40
		41	// =======================================================================================================
		42
		43	void plotCorTDC(TH1F hp1d, char hTitle, double fitmin, double fitmax, double drawmin, double drawmax)
		44	{
		45	char sbuff[256];
		46
		47	TLegend *leg1 = new TLegend(0.57,0.76,0.88,0.88);
		48	leg1->SetFillColor(0);
		49	leg1->SetBorderSize(1);
		50	leg1->SetTextSize(0.05);
		51	leg1->SetTextAlign(22);
		52
		53	//gStyle->SetOptFit(111);
		54	sprintf(sbuff, "%s;Time[ns];", hTitle);
		55	hp1d->SetTitle(sbuff);
		56	//~ hp1d->GetXaxis()->SetRangeUser(2fitmin, 2fitmax);
		57	hp1d->GetXaxis()->SetRangeUser(drawmin, drawmax);
		58	//~ hcTDC->DrawClone();
		59
		60	//~ TF1 *fs = new TF1("fs","pol0(0) + gaus(1)");
		61
		62	TF1 fs = new TF1("fs","gaus(0) + ([3]exp((x)/[4]) + [5]) * (0.5+0.5*TMath::Erf((x-[6])/[7]))");
		63	fs->SetParName(0, "G_Const"); fs->SetParName(1, "G_Mean"); fs->SetParName(2, "G_Sigma");
		64	fs->SetParName(3, "Exp_p0"); fs->SetParName(4, "Exp_p1"); fs->SetParName(5, "Exp_p2");
		65	fs->SetParName(6, "Erf_p0"); fs->SetParName(7, "Erf_p1");
		66	//~
		67	TF1 *fg1 = new TF1("fg1", "gaus(0)"); fg1->SetNpx(400);
		68	TF1 fe1 = new TF1("fe1", "[0]exp((x)/[1]) + [2]");
		69	TF1 ferf1 = new TF1("ferf1", "[0](0.5+0.5*TMath::Erf((x-[1])/[2]))");
		70	TF1 fe1erf = new TF1("fe1erf", "([0]exp((x)/[1]) + [2]) * (0.5+0.5*TMath::Erf((x-[3])/[4]))"); fe1erf->SetNpx(400);
		71
		72	double fitampl = hp1d->GetMaximum();
		73	double fitcenter = hp1d->GetBinCenter(hp1d->GetMaximumBin());
		74
		75	fs->SetNpx(400);
		76	fs->SetParameters(fitampl, fitcenter, 0.01, 1e2, -0.1, 1e1, 0, 0.05);
		77
		78	// G + Er + Ex
		79	//fs->SetParLimits( 1, -0.1, 0.1);
		80	//fs->SetParLimits( 2, 0.018, 0.5);
		81	fs->SetParLimits( 4, -0.3, 0);
		82	//fs->SetParLimits( 6, -1e-1, 1e-1);
		83	fs->SetParLimits( 7, 0.04, 0.6);
		84
		85	//~ // G
		86	//~ fs->SetParameters(fitampl/100., fitampl, fitcenter, 0.01);
		87
		88	hp1d->Fit(fs,"Q","",fitmin, fitmax);
		89
		90	//~ printf("Chi / NDF = %lf/%lf\n", fs->GetChisquare(), fs->GetNDF());
		91
		92	fg1->SetParameters(fs->GetParameter(0), fs->GetParameter(1), fs->GetParameter(2));
		93	//~ fe1->SetParameters(fs->GetParameter(3), fs->GetParameter(4), fs->GetParameter(5));
		94	//~ ferf1->SetParameters(1e4, fs->GetParameter(6), fs->GetParameter(7));
		95	fe1erf->SetParameters(fs->GetParameter(3), fs->GetParameter(4), fs->GetParameter(5), fs->GetParameter(6), fs->GetParameter(7));
		96
		97	fg1->SetLineColor(kRed); fg1->SetLineWidth(1); fg1->SetLineStyle(1); fg1->SetRange(fitmin, fitmax); fg1->DrawCopy("SAME");
		98	//~ fe1->SetLineColor(kBlue); fe1->SetLineWidth(1); fe1->SetLineStyle(1); fe1->SetRange(fitmin, fitmax); fe1->DrawCopy("SAME");
		99	//~ ferf1->SetLineColor(kYellow); ferf1->SetLineWidth(1);ferf1->SetLineStyle(1);ferf1->SetRange(fitmin, fitmax);ferf1->DrawCopy("SAME");
		100	fe1erf->SetLineColor(kBlue); fe1erf->SetLineWidth(1);fe1erf->SetLineStyle(1);fe1erf->SetRange(fitmin, fitmax);fe1erf->DrawCopy("SAME");
		101
		102	double HalfMax = fitampl/2.0;
		103	double MinHalfMax = fs->GetX(HalfMax, -0.5, fitcenter);
		104	double MaxHalfMax = fs->GetX(HalfMax, fitcenter, 1.5);
		105	double FWHM = MaxHalfMax - MinHalfMax;
		106	double sigmaFWHM = fs->GetParameter(2);
		107	printf("Corrected %s resolution FWHM = %.1lf [ps]\n", hTitle, FWHM*1000);
		108	//~ fprintf(fp,"%lf ", FWHM*1000);
		109
		110	sprintf(sbuff, "FWHM = %.0lf ps", FWHM*1000);
		111	leg1->AddEntry((TF1*)fs->Clone(),sbuff, "L");
		112	leg1->Draw();
		113	}
		114
		115	// =======================================================================================================
		116
		117	void plots2(char fname="480", char plopt="d", double fitcenter = 1190., double fitw = 70., double comhi_range_hi = 2100., int batch_q=0)
		118	{
		119	char sbuff[256];
		120
		121	char *fitopt="crystal1";
		122	double inmin=-10;
		123	double inmax=10;
		124	double Ecut_lo = 450.;
		125	double Ecut_hi = 625.;
		126
		127
		128	int printeps = 0;
		129	int printgif=-1, double fitlo=1.0, double fithi=1.0;
		130	int ch=0;
		131	char fullname[256];
		132
		133	//get ROOT file with histograms
		134	char fnameroot[1024];
		135	TFile * rootfile;
		136	TDirectory *dir;
		137
		138	sprintf(fnameroot, "root/%s.root", fname);
		139	rootfile = (TFile *) gROOT->FindObject(fname);
		140	if(rootfile==NULL) rootfile = new TFile(fnameroot);
		141	if(rootfile==NULL) {
		142	printf("Cannot open root file %s!!!\n",fnameroot);
		143	return;
		144	}
		145	dir = (TDirectory*) rootfile;
		146
		147	DrSetDrawStyle();
		148
		149	TCanvas *c[64];
		150	int cc=-1;
		151	char hname[256], htitle[256];
		152	TH1F hp1d, hp1dcut; TH2F *hp2d;
		153
		154	TLegend *leg[64];
		155	int legi = -1;
		156	TH1F *legentry[64];
		157	int legentryi = -1;
		158
		159	TGraph *graph[16];
		160	int igraph=-1;
		161
		162	int colarr[20] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, kOrange, 10,11,12,13,14,15,16,17,18,19};
		163	int legsty[] = {1,7,2,3,9,6,9,5,6};
		164
		165	TF1 *fg = new TF1("fg", "gaus"); fg->SetNpx(300);
		166
		167	TF1 *fgg = new TF1("fgg", "gaus(0)+gaus(3)"); fg->SetNpx(300);
		168	fgg->SetParName(0,"G1_Const"); fgg->SetParName(3,"G2_Const");
		169	fgg->SetParName(1,"G1_Mean" ); fgg->SetParName(4,"G2_Mean" );
		170	fgg->SetParName(2,"G1_Sigma"); fgg->SetParName(5,"G2_Sigma");
		171
		172	TF1 *fg3 = new TF1("fg3", "gaus(0)+gaus(3)+gaus(6)"); fg3->SetNpx(300);
		173	fg3->SetParName(0,"G1_Const"); fg3->SetParName(3,"G2_Const"); fg3->SetParName(6,"G3_Const");
		174	fg3->SetParName(1,"G1_Mean" ); fg3->SetParName(4,"G2_Mean" ); fg3->SetParName(7,"G3_Mean" );
		175	fg3->SetParName(2,"G1_Sigma"); fg3->SetParName(5,"G2_Sigma"); fg3->SetParName(8,"G3_Sigma");
		176
		177	TF1 *fg4 = new TF1("fg4", "gaus(0)+gaus(3)+gaus(6)+gaus(9)"); fg4->SetNpx(300);
		178	fg4->SetParName(0,"G1_Const"); fg4->SetParName(3,"G2_Const"); fg4->SetParName(6,"G3_Const"); fg4->SetParName(9,"G4_Const");
		179	fg4->SetParName(1,"G1_Mean" ); fg4->SetParName(4,"G2_Mean" ); fg4->SetParName(7,"G3_Mean" ); fg4->SetParName(10,"G4_Mean" );
		180	fg4->SetParName(2,"G1_Sigma"); fg4->SetParName(5,"G2_Sigma"); fg4->SetParName(8,"G3_Sigma"); fg4->SetParName(11,"G4_Sigma");
		181
		182	TF1 *fg5 = new TF1("fg5", "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)"); fg5->SetNpx(300);
		183	fg5->SetParName(0,"G1_Const"); fg5->SetParName(3,"G2_Const"); fg5->SetParName(6,"G3_Const"); fg5->SetParName(9,"G4_Const"); fg5->SetParName(12,"G5_Const");
		184	fg5->SetParName(1,"G1_Mean" ); fg5->SetParName(4,"G2_Mean" ); fg5->SetParName(7,"G3_Mean" ); fg5->SetParName(10,"G4_Mean" ); fg5->SetParName(13,"G5_Mean" );
		185	fg5->SetParName(2,"G1_Sigma"); fg5->SetParName(5,"G2_Sigma"); fg5->SetParName(8,"G3_Sigma"); fg5->SetParName(11,"G4_Sigma"); fg5->SetParName(14,"G5_Sigma");
		186
		187	// -------------------------------------------------------------------
		188
		189	if( strchr(plopt, 'a') != NULL ) {
		190
		191	//~ gStyle->SetOptStat(1111);
		192	gStyle->SetOptStat(0);
		193	gStyle->SetOptFit(0);
		194
		195	c[++cc] = new TCanvas("a", "a", 0, 0, 1200, 700);
		196	c[cc]->Divide(2,2);
		197
		198	int draw_cuts = 1;
		199	int rebin_tdcs = 0;
		200
		201	int ccccdi = 0;
		202	int ich;
		203
		204	// -------------------------- qdc 0 ----------------------------------
		205	ich = 0;
		206	(c[cc]->cd(++ccccdi))->SetLogy(1);
		207	sprintf(hname, "hadc%d", ich);
		208	hp1d = DrTH1F(dir, hname, "");
		209	hp1d->SetTitle("MPPC 1;Charge [200 fC/bin]");
		210	hp1d->GetXaxis()->SetRangeUser(0,2000);
		211	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		212	hp1d->SetLineColor(kBlack);
		213	hp1d->DrawClone();
		214
		215	if(draw_cuts) {
		216	sprintf(hname, "hadc_cut%d", ich);
		217	hp1dcut = DrTH1F(dir, hname, "");
		218	hp1dcut->SetLineColor(kRed);
		219	hp1dcut->DrawClone("SAME");
		220
		221	//~ sprintf(hname, "hadc_cut_2%d", ich);
		222	//~ hp1dcut = DrTH1F(dir, hname, "");
		223	//~ hp1dcut->SetLineColor(kGreen);
		224	//~ hp1dcut->DrawClone("SAME");
		225	}
		226
		227	// -------------------------- TDC 0 ----------------------------------
		228	(c[cc]->cd(++ccccdi))->SetLogy(1);
		229	sprintf(hname, "htdc%d", ich);
		230	hp1d = DrTH1F(dir, hname, "");
		231	hp1d->SetTitle("MPPC 1;Time [ns]");
		232	if(rebin_tdcs) hp1d->Rebin(rebin_tdcs);
		233	//~ hp1d->GetXaxis()->SetRangeUser(0,22);
		234	hp1d->GetXaxis()->SetRangeUser(20,34);
		235	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		236	hp1d->SetLineColor(kBlack);
		237	hp1d->DrawClone();
		238
		239	if(draw_cuts) {
		240	sprintf(hname, "htdc_cut%d", ich);
		241	hp1dcut = DrTH1F(dir, hname, "");
		242	if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
		243	hp1dcut->SetLineColor(kRed);
		244	hp1dcut->DrawClone("SAME");
		245
		246	//~ sprintf(hname, "htdc_cut_2%d", ich);
		247	//~ hp1dcut = DrTH1F(dir, hname, "");
		248	//~ if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
		249	//~ hp1dcut->SetLineColor(kGreen);
		250	//~ hp1dcut->DrawClone("SAME");
		251	}
		252
		253	// -------------------------- qdc 1 ----------------------------------
		254	ich = 1;
		255	(c[cc]->cd(++ccccdi))->SetLogy(1);
		256	sprintf(hname, "hadc%d", ich);
		257	hp1d = DrTH1F(dir, hname, "");
		258	hp1d->SetTitle("MPPC 2;Charge [200 fC/bin]");
		259	hp1d->GetXaxis()->SetRangeUser(0,2000);
		260	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		261	hp1d->SetLineColor(kBlack);
		262	hp1d->DrawClone();
		263
		264	if(draw_cuts) {
		265	sprintf(hname, "hadc_cut%d", ich);
		266	hp1dcut = DrTH1F(dir, hname, "");
		267	hp1dcut->SetLineColor(kRed);
		268	hp1dcut->DrawClone("SAME");
		269
		270	//~ sprintf(hname, "hadc_cut_2%d", ich);
		271	//~ hp1dcut = DrTH1F(dir, hname, "");
		272	//~ hp1dcut->SetLineColor(kGreen);
		273	//~ hp1dcut->DrawClone("SAME");
		274	}
		275
		276	// -------------------------- TDC 1 ----------------------------------
		277	(c[cc]->cd(++ccccdi))->SetLogy(1);
		278	sprintf(hname, "htdc%d", ich);
		279	hp1d = DrTH1F(dir, hname, "");
		280	hp1d->SetTitle("MPPC 2;Time [ns]");
		281	if(rebin_tdcs) hp1d->Rebin(rebin_tdcs);
		282	//~ hp1d->GetXaxis()->SetRangeUser(0,22);
		283	hp1d->GetXaxis()->SetRangeUser(20,34);
		284	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		285	hp1d->SetLineColor(kBlack);
		286	hp1d->DrawClone();
		287
		288	if(draw_cuts) {
		289	sprintf(hname, "htdc_cut%d", ich);
		290	hp1dcut = DrTH1F(dir, hname, "");
		291	if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
		292	hp1dcut->SetLineColor(kRed);
		293	hp1dcut->DrawClone("SAME");
		294
		295	//~ sprintf(hname, "htdc_cut_2%d", ich);
		296	//~ hp1dcut = DrTH1F(dir, hname, "");
		297	//~ if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
		298	//~ hp1dcut->SetLineColor(kGreen);
		299	//~ hp1dcut->DrawClone("SAME");
		300	}
		301
		302	sprintf(fullname, "gif/a_%s.gif", fname); c[cc]->SaveAs(fullname);
		303	}
		304
		305	// -------------------------------------------------------------------
		306
		307	if( strchr(plopt, 'q') != NULL ) {
		308
		309	gStyle->SetOptStat(1111);
		310	gStyle->SetOptFit(0);
		311
		312	c[++cc] = new TCanvas("q", "q", 0, 0, 1200, 700);
		313	c[cc]->Divide(3,2);
		314
		315	int draw_cuts = 1;
		316	int rebin_tdcs = 1;
		317
		318	int ccccdi = 0;
		319	int ich;
		320	/*
		321	ich = 2;
		322	// -------------------------- qdc ref ----------------------------------
		323	(c[cc]->cd(++ccccdi))->SetLogy(1);
		324	sprintf(hname, "hadc%d", ich);
		325	hp1d = DrTH1F(dir, hname, "");
		326	//~ hp1d->GetXaxis()->SetRangeUser(0,100);
		327	hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		328	hp1d->DrawClone();
		329
		330	sprintf(hname, "hadc_cut%d", ich);
		331	hp1dcut = DrTH1F(dir, hname, "");
		332	hp1dcut->SetLineColor(kRed);
		333	hp1dcut->DrawClone("SAME");
		334
		335	sprintf(hname, "hadc_cut_2%d", ich);
		336	hp1dcut = DrTH1F(dir, hname, "");
		337	hp1dcut->SetLineColor(kGreen);
		338	hp1dcut->DrawClone("SAME");
		339
		340	// -------------------------- TDC ref ----------------------------------
		341	(c[cc]->cd(++ccccdi))->SetLogy(1);
		342	sprintf(hname, "htdc%d", ich);
		343	hp1d = DrTH1F(dir, hname, "");
		344	hp1d->Rebin(8);
		345	hp1d->GetXaxis()->SetRangeUser(5,105);
		346	hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		347	hp1d->DrawClone();
		348
		349	sprintf(hname, "htdc_cut%d", ich);
		350	hp1dcut = DrTH1F(dir, hname, "");
		351	hp1dcut->Rebin(8);
		352	hp1dcut->SetLineColor(kRed);
		353	hp1dcut->DrawClone("SAME");
		354
		355	sprintf(hname, "htdc_cut_2%d", ich);
		356	hp1dcut = DrTH1F(dir, hname, "");
		357	hp1dcut->Rebin(8);
		358	hp1dcut->SetLineColor(kGreen);
		359	hp1dcut->DrawClone("SAME");
		360	*/
		361
		362	// -------------------------- TDC DIFF ----------------------------------
		363	(c[cc]->cd(++ccccdi))->SetLogy(0);
		364	sprintf(hname, "htdcdiff");
		365	hp1d = DrTH1F(dir, hname, "");
		366	hp1d->Rebin(rebin_tdcs);
		367	hp1d->GetXaxis()->SetRangeUser(-15,15);
		368	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*2);
		369	hp1d->GetYaxis()->SetRangeUser(0.,hp1d->GetMaximum()*1.1);
		370	hp1d->DrawClone();
		371
		372	if(draw_cuts) {
		373	sprintf(hname, "htdcdiff_cut");
		374	hp1dcut = DrTH1F(dir, hname, "");
		375	hp1dcut->Rebin(rebin_tdcs);
		376	hp1dcut->SetLineColor(kRed);
		377	hp1dcut->DrawClone("SAME");
		378
		379	sprintf(hname, "htdcdiff_cut_2");
		380	hp1dcut = DrTH1F(dir, hname, "");
		381	hp1dcut->Rebin(rebin_tdcs);
		382	hp1dcut->SetLineColor(kGreen);
		383	hp1dcut->DrawClone("SAME");
		384	}
		385
		386	// -------------------------- qdc MPPC ----------------------------------
		387	ich = 0;
		388	(c[cc]->cd(++ccccdi))->SetLogy(1);
		389	sprintf(hname, "hadc%d", ich);
		390	hp1d = DrTH1F(dir, hname, "");
		391	hp1d->GetXaxis()->SetRangeUser(0,1500);
		392	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		393	hp1d->DrawClone();
		394
		395	if(draw_cuts) {
		396	sprintf(hname, "hadc_cut%d", ich);
		397	hp1dcut = DrTH1F(dir, hname, "");
		398	hp1dcut->SetLineColor(kRed);
		399	hp1dcut->DrawClone("SAME");
		400
		401	sprintf(hname, "hadc_cut_2%d", ich);
		402	hp1dcut = DrTH1F(dir, hname, "");
		403	hp1dcut->SetLineColor(kGreen);
		404	hp1dcut->DrawClone("SAME");
		405	}
		406
		407	// -------------------------- TDC MPPC ----------------------------------
		408	(c[cc]->cd(++ccccdi))->SetLogy(1);
		409	sprintf(hname, "htdc%d", ich);
		410	hp1d = DrTH1F(dir, hname, "");
		411	hp1d->Rebin(rebin_tdcs);
		412	hp1d->GetXaxis()->SetRangeUser(0,20);
		413	hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		414	hp1d->DrawClone();
		415
		416	if(draw_cuts) {
		417	sprintf(hname, "htdc_cut%d", ich);
		418	hp1dcut = DrTH1F(dir, hname, "");
		419	hp1dcut->Rebin(rebin_tdcs);
		420	hp1dcut->SetLineColor(kRed);
		421	hp1dcut->DrawClone("SAME");
		422
		423	sprintf(hname, "htdc_cut_2%d", ich);
		424	hp1dcut = DrTH1F(dir, hname, "");
		425	hp1dcut->Rebin(rebin_tdcs);
		426	hp1dcut->SetLineColor(kGreen);
		427	hp1dcut->DrawClone("SAME");
		428	}
		429
		430	// -------------------------- qdc MPPC ----------------------------------
		431	ich = 1;
		432	ccccdi++;
		433	(c[cc]->cd(++ccccdi))->SetLogy(1);
		434	sprintf(hname, "hadc%d", ich);
		435	hp1d = DrTH1F(dir, hname, "");
		436	hp1d->GetXaxis()->SetRangeUser(0,1500);
		437	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		438	hp1d->DrawClone();
		439
		440	if(draw_cuts) {
		441	sprintf(hname, "hadc_cut%d", ich);
		442	hp1dcut = DrTH1F(dir, hname, "");
		443	hp1dcut->SetLineColor(kRed);
		444	hp1dcut->DrawClone("SAME");
		445
		446	sprintf(hname, "hadc_cut_2%d", ich);
		447	hp1dcut = DrTH1F(dir, hname, "");
		448	hp1dcut->SetLineColor(kGreen);
		449	hp1dcut->DrawClone("SAME");
		450	}
		451
		452	// -------------------------- TDC MPPC ----------------------------------
		453	(c[cc]->cd(++ccccdi))->SetLogy(1);
		454	sprintf(hname, "htdc%d", ich);
		455	hp1d = DrTH1F(dir, hname, "");
		456	hp1d->Rebin(rebin_tdcs);
		457	hp1d->GetXaxis()->SetRangeUser(0,20);
		458	hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		459	hp1d->DrawClone();
		460
		461	if(draw_cuts) {
		462	sprintf(hname, "htdc_cut%d", ich);
		463	hp1dcut = DrTH1F(dir, hname, "");
		464	hp1dcut->Rebin(rebin_tdcs);
		465	hp1dcut->SetLineColor(kRed);
		466	hp1dcut->DrawClone("SAME");
		467
		468	sprintf(hname, "htdc_cut_2%d", ich);
		469	hp1dcut = DrTH1F(dir, hname, "");
		470	hp1dcut->Rebin(rebin_tdcs);
		471	hp1dcut->SetLineColor(kGreen);
		472	hp1dcut->DrawClone("SAME");
		473	}
		474
		475	sprintf(fullname, "gif/q_%s.gif", fname); c[cc]->SaveAs(fullname);
		476	//~ sprintf(fullname, "%s_j2.eps", fname); c[cc]->SaveAs(fullname);
		477
		478	if(batch_q) gSystem->Exit(1);
		479	}
		480	// -------------------------------------------------------------------
		481
		482	if( strchr(plopt, 'n') != NULL ) {
		483
		484	gStyle->SetOptStat("e");
		485	//~ gStyle->SetOptFit(1);
		486
		487	c[++cc] = new TCanvas("d", "d", 0, 0, 900, 600);
		488	//~ c[cc]->Divide(2,2);
		489
		490	int draw_cuts = 0;
		491	int rebin_tdcs = 1;
		492
		493	double drawmin = -17;
		494	double drawmax = +17;
		495	double fitmin = -10;
		496	double fitmax = 10;
		497	double fitcenter = 0.0;
		498	double HalfMax, MinHalfMax, MaxHalfMax, FWHM, FWHMc;
		499
		500	int ccccdi = 0;
		501	int ich = 0;
		502
		503	TF1 *f_fit;
		504
		505	TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)"); fcg->SetNpx(300);
		506	fcg->SetParName(0,"Const"); fcg->SetParName(1,"G_Const");
		507	fcg->SetParName(2,"G_Mean" ); fcg->SetParName(3,"G_Sigma" );
		508
		509	TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)"); fcgg->SetNpx(300);
		510	fcgg->SetParName(0,"Const");
		511	fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
		512	fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
		513
		514	(c[cc]->cd(++ccccdi))->SetLogy(0);
		515	sprintf(hname, "htdcdiff");
		516	hp1d = DrTH1F(dir, hname, "");
		517	//~ hp1d->Rebin(rebin_tdcs);
		518	hp1d->SetLineColor(kBlack);
		519	hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
		520	hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*1.1);
		521	hp1d->SetTitle(";Coincidence Time [ns];Counts");
		522
		523
		524	f_fit = fcgg;
		525	f_fit->SetParameters(hp1d->GetMaximum()*0.25,
		526	hp1d->GetMaximum()*0.7, 0, 0.4,
		527	hp1d->GetMaximum()*0.05, 0, 0.05);
		528	f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
		529	f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
		530	f_fit->SetParLimits(2, fitmin,fitmax);
		531	f_fit->SetParLimits(3, 0.2, 1);
		532	f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
		533	f_fit->SetParLimits(5, fitmin,fitmax);
		534	f_fit->SetParLimits(6, 0.02, 0.3);
		535
		536	//~ f_fit = fcgg;
		537	//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		538	//~ hp1d->GetMaximum()*0.5, 0, 0.05,
		539	//~ hp1d->GetMaximum()*0.4, 0, 0.5);
		540
		541	//~ f_fit = fcg;
		542	//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		543	//~ hp1d->GetMaximum()*0.9, 0, 0.05);
		544
		545	f_fit->SetRange(fitmin,fitmax);
		546	f_fit->SetLineWidth(1);
		547	hp1d->Fit(f_fit, "Q", "", fitmin,fitmax);
		548
		549
		550	HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
		551	MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
		552	MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
		553	FWHM = MaxHalfMax - MinHalfMax;
		554
		555	HalfMax = (f_fit->GetMaximum(fitmin, fitmax) + f_fit->GetParameter(0))/2.0;
		556	MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
		557	MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
		558	FWHMc = MaxHalfMax - MinHalfMax;
		559
		560	//~ printf("FWHM = %lf ns \| FWHMc = %lf ns\n", FWHM, FWHMc);
		561	printf("FWHM (peak only) = %.3lf ns\n", FWHMc);
		562
		563
		564	leg[++legi] = new TLegend(0.122768, 0.740418, 0.420759, 0.878049);
		565	leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.03); leg[legi]->SetTextAlign(22);
		566
		567
		568	for(int i=0; i<6; i++) fgg->SetParameter(i, f_fit->GetParameter(i+1));
		569	//~ fgg->SetRange(drawmin, drawmax); fgg->SetNpx(300); fgg->SetLineColor(kBlue); fgg->SetLineWidth(1); fgg->DrawClone("LSAME");
		570	double Icherenkovs = fgg->Integral(drawmin, drawmax);
		571
		572	double Iall = f_fit->Integral(fitmin, fitmax);
		573	double Rsn = Icherenkovs/(Iall-Icherenkovs);
		574
		575	double Ientr = hp1d->GetEntries()*hp1d->GetBinWidth(1);
		576	double Rsne = Icherenkovs/(Ientr-Icherenkovs);
		577
		578	double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1)*hp1d->GetBinWidth(1);
		579	double Rsnh = Icherenkovs/(Ihist-Icherenkovs);
		580
		581	//~ printf("ICherenkov = %lf \| IAll = %lf \| Ientr = %lf \| Ihist = %lf\n", Icherenkovs, Iall, Ientr, Ihist);
		582	//~ printf("TF1 -> S/N = %lf \|\| TH1F(e) -> S/N = %lf \|\| TH1F(i) -> S/N = %lf\n", Rsn, Rsne, Rsnh);
		583	printf("S/N (vs. All) = %lf\n", Rsnh);
		584
		585	if(draw_cuts) {
		586	sprintf(hname, "htdcdiff_cut");
		587	hp1dcut = DrTH1F(dir, hname, "");
		588	hp1dcut->Rebin(rebin_tdcs);
		589	hp1dcut->SetLineColor(kRed);
		590	hp1dcut->DrawClone("SAME");
		591
		592	sprintf(hname, "htdcdiff_cut_2");
		593	hp1dcut = DrTH1F(dir, hname, "");
		594	hp1dcut->Rebin(rebin_tdcs);
		595	hp1dcut->SetLineColor(kGreen);
		596	hp1dcut->DrawClone("SAME");
		597	}
		598
		599	//~ #define USE_NOISE_FILE
		600
		601	#ifdef USE_NOISE_FILE
		602	//get ROOT file with histograms
		603	char fnameroot[1024];
		604	TFile * rootfile2;
		605	TDirectory *dir2;
		606
		607	sprintf(fname, "run054");
		608	sprintf(fnameroot, "root/%s.root", fname);
		609	rootfile2 = (TFile *) gROOT->FindObject(fname);
		610	if(rootfile2==NULL) rootfile2 = new TFile(fnameroot);
		611	if(rootfile2==NULL) {
		612	printf("Cannot open root file 2 %s!!!\n",fnameroot);
		613	return;
		614	}
		615	dir2 = (TDirectory*) rootfile2;
		616
		617	sprintf(hname, "htdcdiff");
		618	hp1dcut = DrTH1F(dir2, hname, "");
		619	hp1dcut->SetLineColor(kMagenta);
		620	hp1dcut->DrawClone("SAME");
		621
		622	leg[legi]->AddEntry((TH1F*)hp1d->Clone(), "With source (30 min)", "L");
		623	leg[legi]->AddEntry((TH1F*)hp1dcut->Clone(), "W/o source (30 min)", "L");
		624	sprintf(fullname, "Fit FWHM = %.3lf ns", FWHMc);
		625	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "L");
		626	leg[legi]->Draw();
		627	#endif
		628
		629	sprintf(fullname, "gif/n_%s.gif", fname); c[cc]->SaveAs(fullname);
		630
		631	if(batch_q) gSystem->Exit(1);
		632	}
		633	// -------------------------------------------------------------------
		634
		635	if( strchr(plopt, 'f') != NULL ) {
		636
		637	DrSetDrawStyle(0.05);
		638
		639	gStyle->SetOptStat("e");
		640	//~ gStyle->SetOptFit(1);
		641
		642	gStyle->SetPadRightMargin(0.03);
		643	gStyle->SetPadTopMargin(0.05);
		644	gStyle->SetPadBottomMargin(0.12);
		645	gStyle->SetStatFontSize(0.05);
		646	gStyle->SetStatX(0.97);
		647	gStyle->SetStatY(0.95);
		648
		649	c[++cc] = new TCanvas("f", "f", 0, 460, 720, 400);
		650	//~ c[cc]->Divide(2,2);
		651
		652	int draw_cuts = 0;
		653	int rebin_tdcs = 1;
		654
		655	double drawmin = -10;
		656	double drawmax = +10;
		657	//~ double fitmin = inmin;
		658	//~ double fitmax = inmax;
		659	double fitmin = -10;
		660	double fitmax = 10;
		661	double fitcenter = 0.0;
		662	double HalfMax, MinHalfMax, MaxHalfMax, FWHM, FWHMc;
		663
		664	int ccccdi = 0;
		665	int ich = 0;
		666
		667	TF1 *f_fit;
		668	TF1 *f_integrate;
		669
		670	TF1 *fc = new TF1("fc", "pol0(0)");
		671
		672	TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)"); fcg->SetNpx(300);
		673	fcg->SetParName(0,"Const"); fcg->SetParName(1,"G_Const");
		674	fcg->SetParName(2,"G_Mean" ); fcg->SetParName(3,"G_Sigma" );
		675
		676	TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)"); fcgg->SetNpx(300);
		677	fcgg->SetParName(0,"Const");
		678	fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
		679	fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
		680
		681	TF1 *fcggg = new TF1("fcggg", "pol0(0)+gaus(1)+gaus(4)+gaus(7)"); fcggg->SetNpx(300);
		682	fcggg->SetParName(0,"Const");
		683	fcggg->SetParName(1,"G1_Const"); fcggg->SetParName(2,"G1_Mean" ); fcggg->SetParName(3,"G1_Sigma" );
		684	fcggg->SetParName(4,"G2_Const"); fcggg->SetParName(5,"G2_Mean" ); fcggg->SetParName(6,"G2_Sigma" );
		685	fcggg->SetParName(7,"G3_Const"); fcggg->SetParName(8,"G3_Mean" ); fcggg->SetParName(9,"G3_Sigma" );
		686
		687	//~ TF1 f_nocrystal = new TF1("f_nocrystal", "[0]exp(-0.5((x-[1])/(sqrt(2)[2]))*2)+[3]exp(-0.5((x-[4])/(sqrt(2)[5]))*2)+[6]exp(-1((x-[7])2)/(2(([2])2+([5])2)))+[8]exp(-1(((x-[9])*2)/(2(([2])2+([5])2))))+[10]");
		688	TF1 f_nocrystal = new TF1("f_nocrystal", "[0]([1]*2TMath::Gaus(x,[2],sqrt(2)[3],1)+(1-[1])2TMath::Gaus(x,[4],sqrt(2)[5],1)+[1](1-[1])TMath::Gaus(x,[6],sqrt([3]2+[5]2),1)+[1](1-[1])TMath::Gaus(x,[7],sqrt([3]2+[5]*2),1))+[8]");
		689	f_nocrystal->SetParName(0, "const");
		690	f_nocrystal->SetParName(1, "ratio");
		691	f_nocrystal->SetParName(2, "WW_mean");
		692	f_nocrystal->SetParName(3, "Window_sigma");
		693	f_nocrystal->SetParName(4, "MM_mean");
		694	f_nocrystal->SetParName(5, "MCP_sigma");
		695	f_nocrystal->SetParName(6, "WM1_mean");
		696	f_nocrystal->SetParName(7, "WM2_mean");
		697	f_nocrystal->SetParName(8, "Constant");
		698
		699	//~ TF1 f_crystal = new TF1("f_crystal", "[0]+[1]([2]*2(TMath::Gaus(x,[3],sqrt(2)[4],1))+([5])2(TMath::Gaus(x,[6],sqrt(2)[7],1))+(4[5]/6)*2(TMath::Gaus(x,[8],sqrt(2)[9],1))+[2][5](TMath::Gaus(x,[10],sqrt([2]2+[5]2),1))+[5][2](TMath::Gaus(x,[11],sqrt([2]2+[5]2),1))+[2](4/6[5])(TMath::Gaus(x,[12],sqrt([2]*2+(4/6[5])*2),1))+[2](4/6[5])(TMath::Gaus(x,[13],sqrt([2]*2+(4/6[5])*2),1))+[5](4/6[5])(TMath::Gaus(x,[14],sqrt([5]*2+(4/6[5])*2),1))+[5](4/6[5])(TMath::Gaus(x,[15],sqrt([5]*2+(4/6[5])**2),1)))");
		700	TF1 f_crystal = new TF1("f_crystal", "[0]+[1]([2]*2(TMath::Gaus(x,[3],sqrt(2)[4],1))+([5])2(TMath::Gaus(x,[6],sqrt(2)[7],1))+(4[5]/6)*2(TMath::Gaus(x,[8],sqrt(2)[9],1))+[2][5](TMath::Gaus(x,[10],sqrt([4]2+[7]2),1))+[5][2](TMath::Gaus(x,[11],sqrt([4]2+[7]2),1))+[2](4/6[5])(TMath::Gaus(x,[12],sqrt([4]2+[9]2),1))+[2](4/6[5])(TMath::Gaus(x,[13],sqrt([4]2+[9]2),1))+[5](4/6[5])(TMath::Gaus(x,[14],sqrt([7]2+[9]2),1))+[5](4/6[5])(TMath::Gaus(x,[15],sqrt([7]2+[9]*2),1)))");
		701	f_crystal->SetParName(0, "dark");
		702	f_crystal->SetParName(1, "const");
		703	f_crystal->SetParName(2, "Prob_Crystal");
		704	f_crystal->SetParName(3, "CC_mean");
		705	f_crystal->SetParName(4, "Crystal_sigma");
		706	f_crystal->SetParName(5, "Prob_window");
		707	f_crystal->SetParName(6, "WW_mean");
		708	f_crystal->SetParName(7, "WW_sigma");
		709	f_crystal->SetParName(8, "MM_mean");
		710	f_crystal->SetParName(9, "MM_sigma");
		711	f_crystal->SetParName(10, "CW1_mean");
		712	f_crystal->SetParName(11, "CW2_mean");
		713	f_crystal->SetParName(12, "CM1_mean");
		714	f_crystal->SetParName(13, "CM2_mean");
		715	f_crystal->SetParName(14, "WM1_mean");
		716	f_crystal->SetParName(15, "WM2_mean");
		717
		718	//~ TF1 f_crystal1 = new TF1("f_crystal1", "[0]+[1]([2]*2TMath::Gaus(x,[3],sqrt(2)[4],1)+(1.-5./3[2])*2TMath::Gaus(x,[5],sqrt(2)[6],1)+(4.[2]/6.)*2TMath::Gaus(x,[7],sqrt(2)[8],1)+[2](1.-(5./3)[2])TMath::Gaus(x,[9],sqrt([4]2+[6]2),1)+(1.-(5./3)[2])[2]TMath::Gaus(x,[10],sqrt([4]2+[6]2),1)+[2]4./6.[2]TMath::Gaus(x,[11],sqrt([4]2+[8]2),1)+[2]4/6.[2]TMath::Gaus(x,[12],sqrt([4]2+[8]2),1)+4./6[2](1.-5./3[2])TMath::Gaus(x,[13],sqrt([6]2+[8]2),1)+(1.-5./3[2])4./6[2]TMath::Gaus(x,[14],sqrt([8]2+[6]*2),1))");
		719	TF1 f_crystal1 = new TF1("f_crystal1", "[0]+[1]([2]*2TMath::Gaus(x,[3],sqrt(2)[4],1)+(1.-(1.+[15])[2])*2TMath::Gaus(x,[5],sqrt(2)[6],1)+([15][2])*2TMath::Gaus(x,[7],sqrt(2)[8],1)+[16]2[2](1.-(1+[15])[2])TMath::Gaus(x,[9],sqrt([4]2+[6]2),1)+[16]2(1.-(1+[15])[2])[2]TMath::Gaus(x,[10],sqrt([4]2+[6]2),1)+[2][15][2]TMath::Gaus(x,[11],sqrt([4]2+[8]2),1)+[2][15][2]TMath::Gaus(x,[12],sqrt([4]2+[8]2),1)+[15][2](1.-(1+[15])[2])TMath::Gaus(x,[13],sqrt([6]2+[8]2),1)+(1.-(1.+[15])[2])[15][2]TMath::Gaus(x,[14],sqrt([8]2+[6]*2),1))");
		720	f_crystal1->SetParName(0, "dark");
		721	f_crystal1->SetParName(1, "const");
		722	f_crystal1->SetParName(2, "Prob_Window");
		723	f_crystal1->SetParName(3, "WW_mean");
		724	f_crystal1->SetParName(4, "Window_sigma");
		725	f_crystal1->SetParName(5, "CC_mean");
		726	f_crystal1->SetParName(6, "Crystal_sigma");
		727	f_crystal1->SetParName(7, "MM_mean");
		728	f_crystal1->SetParName(8, "MCP_sigma");
		729	f_crystal1->SetParName(9, "CW1_mean");
		730	f_crystal1->SetParName(10, "CW2_mean");
		731	f_crystal1->SetParName(11, "WM1_mean");
		732	f_crystal1->SetParName(12, "WM2_mean");
		733	f_crystal1->SetParName(13, "CM1_mean");
		734	f_crystal1->SetParName(14, "CM2_mean");
		735	f_crystal1->SetParName(15, "ratio_mcp/window");
		736	f_crystal1->SetParName(16, "Sup._CW");
		737
		738	double ratio_r=0;
		739
		740	(c[cc]->cd(++ccccdi))->SetLogy(0);
		741	sprintf(hname, "htdcdiff_cut");
		742	hp1d = DrTH1F(dir, hname, "");
		743	//~ hp1d->Rebin(rebin_tdcs);
		744	hp1d->SetLineColor(kBlack);
		745	hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
		746	hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*1.1);
		747	sprintf(fullname, ";Coincidence Time [ns];Counts");
		748	hp1d->SetTitle(fullname);
		749
		750	if( strcmp(fitopt, "cgg")==0 ) {
		751	f_fit = fcgg;
		752	f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		753	hp1d->GetMaximum()*0.4, 0, 0.8,
		754	hp1d->GetMaximum()*0.5, 0, 0.05);
		755	f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
		756	f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
		757	f_fit->SetParLimits(2, fitmin,fitmax);
		758	f_fit->SetParLimits(3, 0.05, 1);
		759	f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
		760	f_fit->SetParLimits(5, fitmin,fitmax);
		761	f_fit->SetParLimits(6, 0.01, 1);
		762	} else if(strcmp(fitopt, "cggg")==0) {
		763	f_fit = fcggg;
		764	f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		765	hp1d->GetMaximum()*0.4, 0, 0.4,
		766	hp1d->GetMaximum()*0.4, 0, 0.2,
		767	hp1d->GetMaximum()*0.1, 0, 0.05);
		768	f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
		769	f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
		770	f_fit->SetParLimits(2, fitmin,fitmax);
		771	f_fit->SetParLimits(3, 0.01, 1);
		772	f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
		773	f_fit->SetParLimits(5, fitmin,fitmax);
		774	f_fit->SetParLimits(6, 0.05, 1);
		775	f_fit->SetParLimits(7, 0, hp1d->GetMaximum());
		776	f_fit->SetParLimits(8, fitmin,fitmax);
		777	f_fit->SetParLimits(9, 0.2, 1);
		778	} else if(strcmp(fitopt, "nocrystal")==0)
		779	{ f_fit = f_nocrystal;
		780
		781	//~ f_fit->SetParameters(4000,0.87,0.03,0.03,0.28,0.1,200,200,0.01);
		782	f_fit->SetParameters(200,0.8,0.0,0.03,0.0,0.1,-0.13,0.13,10.16);
		783	f_fit->SetParLimits(0,0,10000);
		784	f_fit->SetParLimits(1,0,1);
		785	f_fit->SetParLimits(2,-0.02,0.02);
		786	f_fit->SetParLimits(3,0.01,0.3);
		787	f_fit->SetParLimits(4,-0.02,0.02);
		788	f_fit->SetParLimits(5,0,0.8);
		789	f_fit->SetParLimits(6,-0.5,0.0);
		790	f_fit->SetParLimits(7,0.0,0.5);
		791	f_fit->SetParLimits(8,0.00,9000.5);
		792	//~ f_fit->FixParameter(6, 100);
		793	//~ f_fit->FixParameter(7, 100);
		794	//~
		795
		796	} else if(strcmp(fitopt, "crystal")==0)
		797	{
		798	f_fit = f_crystal;
		799	//~ f_fit->SetParameters(0.1,300,0.8,0,0.07,0.01,0,0.028,0,0.15,0,0,-0.2,0.2,-0.2,0.2);
		800	f_fit->SetParameter(0, 0.1);
		801	f_fit->SetParameter(1, 300);
		802	f_fit->SetParameter(2, 0.8);
		803	f_fit->SetParameter(3, 0.);
		804	f_fit->SetParameter(4, 0.07);
		805	f_fit->SetParameter(5, 0.01);
		806	f_fit->SetParameter(6, 0.0);
		807	f_fit->SetParameter(7, 0.028);
		808	f_fit->SetParameter(8, 0.0);
		809	f_fit->SetParameter(9, 0.15);
		810	f_fit->SetParameter(10, 0.0);
		811	f_fit->SetParameter(11, 0.0);
		812	f_fit->SetParameter(12, -0.2);
		813	f_fit->SetParameter(13, 0.2);
		814	f_fit->SetParameter(14, -0.2);
		815	f_fit->SetParameter(15, 0.2);
		816	f_fit->SetParameter(15, 0.2);
		817
		818	f_fit->SetParLimits(0, 0, 10);
		819	f_fit->SetParLimits(1, 0, 1000);
		820	f_fit->SetParLimits(2, 0, 1);
		821	f_fit->SetParLimits(3, -0.02, 0.02); f_fit->SetParLimits(4, 0, 0.3);
		822	f_fit->SetParLimits(5, 0, 1); f_fit->SetParLimits(6, -0.02, 0.02);
		823	f_fit->SetParLimits(7, 0, 0.035);
		824	f_fit->SetParLimits(8, -0.02, 0.02);
		825	f_fit->SetParLimits(9, 0, 0.2);
		826	f_fit->SetParLimits(10, -0.02, 0.02);
		827	f_fit->SetParLimits(11, -0.02, 0.02);
		828	f_fit->SetParLimits(12, -0.04, 0);
		829	f_fit->SetParLimits(13, 0, 0.04);
		830	f_fit->SetParLimits(14, -0.04, 0);
		831	f_fit->SetParLimits(15, 0, 0.04);
		832
		833	} else if(strcmp(fitopt, "crystal1")==0)
		834	{
		835	f_fit = f_crystal1;
		836	//~ f_fit->SetParameters(0.1,300,0.8,0,0.07,0.01,0,0.028,0,0.15,0,0,-0.2,0.2,-0.2,0.2);
		837	f_fit->SetParameter(0, 100.1);
		838	f_fit->SetParameter(1, 1);
		839	f_fit->SetParameter(2, 0.2); //prob. window
		840	f_fit->SetParameter(3, 0.0);
		841	f_fit->SetParameter(4, 0.02); // window sigma
		842	f_fit->SetParameter(5, 0.0);
		843	f_fit->SetParameter(6, 120/3.3/1000); //crystal sigma
		844	f_fit->SetParameter(7, 0.0);
		845	f_fit->SetParameter(8, 0.15); // mcp sigma
		846	f_fit->SetParameter(9, -0.1825); //cw
		847	f_fit->SetParameter(10, 0.1825);
		848	f_fit->SetParameter(11, -0.2); //wm
		849	f_fit->SetParameter(12, 0.2);
		850	f_fit->SetParameter(13, -0.3); //cm
		851	f_fit->SetParameter(14, 0.3);
		852	f_fit->SetParameter(15, 1/0.6795 -1.);
		853	f_fit->SetParameter(16, 0.1);
		854	f_fit->SetParameter(1, hp1d->GetMaximum()/(f_fit->Eval(0)-f_fit->GetParameter(0)));
		855
		856
		857	f_fit->SetParLimits(0, 0, 400);
		858	f_fit->SetParLimits(1, 0, 100000);
		859	f_fit->SetParLimits(2, 0, 1.);
		860	f_fit->SetParLimits(3, -0.025, 0.025);
		861	f_fit->SetParLimits(4, 0.02, 0.08);
		862	f_fit->SetParLimits(5, -0.025, 0.025);
		863	f_fit->SetParLimits(6, 0.01, 0.2);
		864	f_fit->SetParLimits(7, -0.025, 0.025);
		865	f_fit->SetParLimits(8, 0, 0.5);
		866	f_fit->SetParLimits(9, -0.6, 0.02);
		867	f_fit->SetParLimits(10, -0.02, 0.6);
		868	f_fit->SetParLimits(11, -0.3, -0.1);
		869	f_fit->SetParLimits(12, 0.1, 0.3);
		870	f_fit->SetParLimits(13, -0.7, -0.1);
		871	f_fit->SetParLimits(14, 0.1, 0.7);
		872	f_fit->SetParLimits(15, 0, 0.999);
		873	f_fit->SetParLimits(16, 0, 1);
		874	f_fit->FixParameter(4,0.02668);
		875	//~ f_fit->FixParameter(6,0.048);
		876	f_fit->FixParameter(8,0.0877);
		877	f_fit->FixParameter(11,-0.19);
		878	f_fit->FixParameter(12,0.19);
		879	//~
		880	ratio_r = 1/0.6795 -1.;
		881	//~ ratio_r = 1/0.8089 -1.;
		882	f_fit->FixParameter(15,ratio_r);
		883	//~ f_fit->FixParameter(16,1);
		884	//~ f_fit->FixParameter(2,0.52);
		885	//~ f_fit->FixParameter(6,170/3.3/1000);
		886	//~ f_fit->FixParameter(3,0);
		887	//~ f_fit->FixParameter(5,0);
		888	//~ f_fit->FixParameter(7,0);
		889
		890	//~ cout<<"hogeeeeeeeeeeee "<<f_fit->GetParameter(1)((1.-5./3.(f_fit->GetParameter(2)))(4./6.(f_fit->GetParameter(2))))<<endl;
		891	}
		892	//~ f_fit->SetRange(fitmin,fitmax);
		893	f_fit->SetRange(-10,10);
		894	f_fit->SetLineWidth(1); f_fit->SetLineStyle(2);
		895	f_fit->SetNpx(300);
		896	//~ hp1d->Fit(f_fit, "QWW", "", fitmin,fitmax);
		897	//~ hp1d->Fit(f_fit, "ILM+", "", -0.4,0.4);
		898	//~
		899	hp1d->Fit(f_fit, "LM+", "", -.5,.5);
		900	//~ hp1d->Draw();
		901	//~ f_fit->Draw("same"); //f_fit->SetLineColor(6);
		902
		903	f_fit->SetLineColor(1);
		904	hp1d->GetXaxis()->SetRangeUser(-1.5,2);
		905
		906	gStyle->SetOptFit(1111);
		907	double wmsigma = sqrt(f_fit->GetParameter(3)2 + f_fit->GetParameter(5)2);
		908
		909	double cm1sigma = sqrt(f_fit->GetParameter(6)2 + f_fit->GetParameter(8)2);
		910	double wm1sigma = sqrt(f_fit->GetParameter(4)2 + f_fit->GetParameter(8)2);
		911	double cw1sigma = sqrt(f_fit->GetParameter(6)2 + f_fit->GetParameter(4)2);
		912
		913
		914	//kobayashi
		915	TPaveStats* st1 = (TPaveStats*) hp1d->FindObject("stats");
		916	st1->SetX1NDC(0.63); st1->SetX2NDC(0.99);
		917	st1->SetY1NDC(0.24); st1->SetY2NDC(0.99);
		918
		919	double pw = f_fit->GetParameter(2);
		920	double pm = f_fit->GetParameter(15)*f_fit->GetParameter(2);
		921	double pc = 1-pw - pm;
		922
		923	double rate = pw2 + pm2 + 2pwpm;
		924	double sumrate = pw2+pm2+pc*2+2pwpm+2pmpc+2pcpwf_fit->GetParameter(16)**2;
		925
		926	cout<<"pw = "<<pw<<" pm = "<<pm<<" pc = "<<pc<<endl;
		927	cout<<"rate = "<<rate<<" sumrate = "<<sumrate<<" rate/sumrate = "<<rate/sumrate<<endl;
		928	cout<<" FWHM_sim = "<<f_fit->GetParameter(6)sqrt(2)2.35<<endl;
		929
		930	if(strcmp(fitopt, "nocrystal")==0)
		931	{
		932	//~ TF1 fg_WW = new TF1("fg_WW","[0]exp(-(x-[1])2/2/([3]2))/sqrt(23.141592[2]*[2])");
		933	//~ TF1 fg_MM = new TF1("fg_MM","[0]exp(-(x-[1])2/2/([3]2))/sqrt(23.141592[2]*[2])");
		934	//~ TF1 fg_WM1 = new TF1("fg_WM1","[0]exp(-(x-[1])2/2/([3]2))/sqrt(23.141592[2]*[2])");
		935	//~ TF1 fg_WM2 = new TF1("fg_WM2","[0]exp(-(x-[1])2/2/([3]2))/sqrt(23.141592[2]*[2])");
		936	TF1 fg_WW = new TF1("fg_WW","[0]TMath::Gaus(x,[1],[2],1)");fg_WW->SetNpx(300);
		937	TF1 fg_MM = new TF1("fg_MM","[0]TMath::Gaus(x,[1],[2],1)");fg_MM->SetNpx(300);
		938	TF1 fg_WM1 = new TF1("fg_WM1","[0]TMath::Gaus(x,[1],[2],1)"); fg_WM1->SetNpx(300);
		939	TF1 fg_WM2 = new TF1("fg_WM2","[0]TMath::Gaus(x,[1],[2],1)"); fg_WM2->SetNpx(300);
		940	fg_WW -> SetParameters(f_fit->GetParameter(1)*2f_fit->GetParameter(0),f_fit->GetParameter(2),sqrt(2)*f_fit->GetParameter(3));
		941	fg_MM -> SetParameters((1-f_fit->GetParameter(1))*2f_fit->GetParameter(0),f_fit->GetParameter(4),sqrt(2)*f_fit->GetParameter(5));
		942	fg_WM1 -> SetParameters(f_fit->GetParameter(0)(1-f_fit->GetParameter(1))f_fit->GetParameter(1),f_fit->GetParameter(6),wmsigma);
		943	fg_WM2 -> SetParameters(f_fit->GetParameter(0)(1-f_fit->GetParameter(1))f_fit->GetParameter(1),f_fit->GetParameter(7),wmsigma);
		944
		945	fg_WW->Draw("sames"); fg_WW->SetLineColor(kYellow+1); fg_WW->SetRange(-10,10);
		946	fg_MM->Draw("sames"); fg_MM->SetLineColor(kMagenta-7); fg_MM->SetRange(-10,10);
		947	fg_WM1->Draw("sames"); fg_WM1->SetLineColor(kRed-4); fg_WM1->SetRange(-10,10);
		948	fg_WM2->Draw("sames"); fg_WM2->SetLineColor(kRed-4); fg_WM2->SetRange(-10,10);
		949
		950	double integ_WW, integ_MM, integ_WM1, integ_WM2;
		951	integ_WW = fg_WW->Integral(-4,4);
		952	}
		953	if(strcmp(fitopt, "crystal1")==0)
		954	{
		955	TF1 fg_CC = new TF1("fg_CC","[0]TMath::Gaus(x,[1],[2],1)"); fg_CC->SetNpx(300);
		956	TF1 fg_WW = new TF1("fg_WW","[0]TMath::Gaus(x,[1],[2],1)"); fg_WW->SetNpx(300);
		957	TF1 fg_MM = new TF1("fg_MM","[0]TMath::Gaus(x,[1],[2],1)"); fg_MM->SetNpx(300);
		958	TF1 fg_WM1 = new TF1("fg_WM1","[0]TMath::Gaus(x,[1],[2],1)"); fg_WM1->SetNpx(300);
		959	TF1 fg_WM2 = new TF1("fg_WM2","[0]TMath::Gaus(x,[1],[2],1)"); fg_WM2->SetNpx(300);
		960	TF1 fg_CW1 = new TF1("fg_CW1","[0]TMath::Gaus(x,[1],[2],1)"); fg_CW1->SetNpx(300);
		961	TF1 fg_CW2 = new TF1("fg_CW2","[0]TMath::Gaus(x,[1],[2],1)"); fg_CW2->SetNpx(300);
		962	TF1 fg_CM1 = new TF1("fg_CM1","[0]TMath::Gaus(x,[1],[2],1)"); fg_CM1->SetNpx(300);
		963	TF1 fg_CM2 = new TF1("fg_CM2","[0]TMath::Gaus(x,[1],[2],1)"); fg_CM2->SetNpx(300);
		964
		965
		966	// normalize
		967	fg_CC -> SetParameters((1.-(1+f_fit->GetParameter(15))f_fit->GetParameter(2))2f_fit->GetParameter(1),f_fit->GetParameter(5),sqrt(2)*f_fit->GetParameter(6));
		968	fg_WW -> SetParameters(f_fit->GetParameter(2)*2f_fit->GetParameter(1),f_fit->GetParameter(3),sqrt(2)*f_fit->GetParameter(4));
		969	fg_MM -> SetParameters((f_fit->GetParameter(15)f_fit->GetParameter(2))2f_fit->GetParameter(1),f_fit->GetParameter(7),sqrt(2)*f_fit->GetParameter(8));
		970	fg_WM1 -> SetParameters(f_fit->GetParameter(1)(f_fit->GetParameter(15)f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(11),wm1sigma);
		971	fg_WM2 -> SetParameters(f_fit->GetParameter(1)(f_fit->GetParameter(15)f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(12),wm1sigma);
		972	fg_CW1 -> SetParameters(f_fit->GetParameter(1)f_fit->GetParameter(16)(1.-(1+f_fit->GetParameter(15))f_fit->GetParameter(2))f_fit->GetParameter(2),f_fit->GetParameter(9),cw1sigma);
		973	fg_CW2 -> SetParameters(f_fit->GetParameter(1)f_fit->GetParameter(16)(1.-(1+f_fit->GetParameter(15))f_fit->GetParameter(2))f_fit->GetParameter(2),f_fit->GetParameter(10),cw1sigma);
		974	fg_CM1 -> SetParameters(f_fit->GetParameter(1)(ratio_rf_fit->GetParameter(2))(1.-(1+f_fit->GetParameter(15))f_fit->GetParameter(2)),f_fit->GetParameter(13),cm1sigma);
		975	fg_CM2 -> SetParameters(f_fit->GetParameter(1)(ratio_rf_fit->GetParameter(2))(1.-(1+f_fit->GetParameter(15))f_fit->GetParameter(2)),f_fit->GetParameter(14),cm1sigma);
		976
		977	//non normalize
		978	//~ fg_CC -> SetParameters((1.-5./3.f_fit->GetParameter(2))2f_fit->GetParameter(1),f_fit->GetParameter(5),sqrt(2)*f_fit->GetParameter(6));
		979	//~ fg_WW -> SetParameters(f_fit->GetParameter(2)*2f_fit->GetParameter(1),f_fit->GetParameter(3),sqrt(2)*f_fit->GetParameter(4));
		980	//~ fg_MM -> SetParameters((4./6.f_fit->GetParameter(2))2f_fit->GetParameter(1),f_fit->GetParameter(7),sqrt(2)*f_fit->GetParameter(8));
		981	//~ fg_WM1 -> SetParameters(f_fit->GetParameter(1)(4./6.f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(11),wm1sigma);
		982	//~ fg_WM2 -> SetParameters(f_fit->GetParameter(1)(4./6.f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(12),wm1sigma);
		983	//~ fg_CW1 -> SetParameters(f_fit->GetParameter(1)(1.-5./3f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(9),cw1sigma);
		984	//~ fg_CW2 -> SetParameters(f_fit->GetParameter(1)(1.-5./3f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(10),cw1sigma);
		985	//~ fg_CM1 -> SetParameters(f_fit->GetParameter(1)(4./6.f_fit->GetParameter(2))(1.-5./3f_fit->GetParameter(2)),f_fit->GetParameter(13),cm1sigma);
		986	//~ fg_CM2 -> SetParameters(f_fit->GetParameter(1)(4./6.f_fit->GetParameter(2))(1.-5./3f_fit->GetParameter(2)),f_fit->GetParameter(14),cm1sigma);
		987
		988
		989	fg_CC->Draw("sames"); fg_CC->SetLineColor(kCyan-5); fg_CC->SetRange(-10,10);
		990	fg_WW->Draw("sames"); fg_WW->SetLineColor(kYellow+1); fg_WW->SetRange(-10,10);
		991	fg_MM->Draw("sames"); fg_MM->SetLineColor(kMagenta-7); fg_MM->SetRange(-10,10);
		992	fg_CW1->Draw("sames"); fg_CW1->SetLineColor(kGreen+2); fg_CW1->SetRange(-10,10);
		993	fg_CW2->Draw("sames"); fg_CW2->SetLineColor(kGreen+2); fg_CW2->SetRange(-10,10);
		994	fg_WM1->Draw("sames"); fg_WM1->SetLineColor(kRed-4); fg_WM1->SetRange(-10,10);
		995	fg_WM2->Draw("sames"); fg_WM2->SetLineColor(kRed-4); fg_WM2->SetRange(-10,10);
		996	fg_CM1->Draw("sames"); fg_CM1->SetLineColor(kBlue+2); fg_CM1->SetRange(-10,10);
		997	fg_CM2->Draw("sames"); fg_CM2->SetLineColor(kBlue+2); fg_CM2->SetRange(-10,10);
		998	}
		999
		1000	HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
		1001	MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
		1002	MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
		1003	FWHM = MaxHalfMax - MinHalfMax;
		1004
		1005	HalfMax = (f_fit->GetMaximum(fitmin, fitmax) + f_fit->GetParameter(0))/2.0;
		1006	MinHalfMax = f_fit->GetX(HalfMax, fitmin, f_fit->GetMaximumX(fitmin, fitmax));
		1007	MaxHalfMax = f_fit->GetX(HalfMax, f_fit->GetMaximumX(fitmin, fitmax), fitmax);
		1008	FWHMc = MaxHalfMax - MinHalfMax;
		1009
		1010	printf("(f) HalfMax = %.3lf ns \| MinHalfMax = %.3lf ns \| MaxHalfMax = %.3lf ns\n", HalfMax, MinHalfMax, MaxHalfMax);
		1011	printf("(f) FWHM (peak only) = %.3lf ns\n", FWHMc);
		1012
		1013	f_integrate = fgg;
		1014	for(int i=0; i<6; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
		1015	//~ for(int i=0; i<9; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
		1016	//~ fgg->SetRange(drawmin, drawmax); fgg->SetNpx(300); fgg->SetLineColor(kBlue); fgg->SetLineWidth(1); fgg->DrawClone("LSAME");
		1017	double Icherenkovs = f_integrate->Integral(drawmin, drawmax)/hp1d->GetBinWidth(1);
		1018
		1019	double Iall = f_fit->Integral(fitmin, fitmax)/hp1d->GetBinWidth(1);
		1020	double Rsn = Icherenkovs/(Iall-Icherenkovs);
		1021
		1022	//~ double Ientr = hp1d->GetEntries()*hp1d->GetBinWidth(1);
		1023	double Ientr = hp1d->GetEntries();
		1024	double Rsne = Icherenkovs/(Ientr-Icherenkovs);
		1025
		1026	//~ double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1)*hp1d->GetBinWidth(1);
		1027	double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1);
		1028	double Rsnh = Icherenkovs/(Ihist-Icherenkovs);
		1029
		1030	double Icherenkovs_10ns = f_integrate->Integral(-10, 10)/hp1d->GetBinWidth(1);
		1031	fc->SetParameter(0, f_fit->GetParameter(0));
		1032	double Inoise_10ns = fc->Integral(-10, 10)/hp1d->GetBinWidth(1);
		1033	double R10ns = Icherenkovs_10ns/Inoise_10ns;
		1034
		1035	double Icherenkovs_4ns = f_integrate->Integral(-4, 4)/hp1d->GetBinWidth(1);
		1036	fc->SetParameter(0, f_fit->GetParameter(0));
		1037	double Inoise_4ns = fc->Integral(-4, 4)/hp1d->GetBinWidth(1);
		1038	double R4ns = Icherenkovs_4ns/Inoise_4ns;
		1039
		1040	double Icherenkovs_2ns = f_integrate->Integral(-2, 2)/hp1d->GetBinWidth(1);
		1041	fc->SetParameter(0, f_fit->GetParameter(0));
		1042	double Inoise_2ns = fc->Integral(-2, 2)/hp1d->GetBinWidth(1);
		1043	double R2ns = Icherenkovs_2ns/Inoise_2ns;
		1044
		1045	//~ printf("ICherenkov = %lf \| IAll = %lf \| Ientr = %lf \| Ihist = %lf\n", Icherenkovs, Iall, Ientr, Ihist);
		1046	printf("Icherenkovs_4ns = %.0lf\n", Icherenkovs_4ns);
		1047	//~ printf("TF1 -> S/N = %lf \|\| TH1F(e) -> S/N = %lf \|\| TH1F(i) -> S/N = %lf\n", Rsn, Rsne, Rsnh);
		1048	//~ printf("S/N (All) = %lf\n", Rsnh);
		1049	printf("S/N (10ns) = %.2lf\n", R10ns);
		1050	printf("S/N (4 ns) = %.2lf\n", R4ns);
		1051	printf("S/N (2 ns) = %.2lf\n", R2ns);
		1052
		1053	leg[++legi] = new TLegend(0.12,0.4,0.43,0.925);
		1054	leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(12);
		1055
		1056	sprintf(fullname, "FWHM_fit = %.0lf ps", 1000*FWHMc);
		1057	//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "L");
		1058	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1059
		1060	sprintf(fullname, "FWHM_c = %.0lf ps", 1000f_fit->GetParameter(6)2.35*sqrt(2));
		1061	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1062
		1063	//~ sprintf(fullname, "S/N (#pm10ns) = %.1lf", R10ns);
		1064	//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1065	//~ sprintf(fullname, "S/N (#pm 4ns) = %.1lf", R4ns);
		1066	//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1067	//~ sprintf(fullname, "S/N (#pm 2ns) = %.1lf", R2ns);
		1068	//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1069
		1070	//kobayashi
		1071	sprintf(fullname, "R_w = %.2f" , pw);
		1072	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1073	sprintf(fullname, "R_m = %.2f" , pm);
		1074	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1075	sprintf(fullname, "R_c = %.2f" , pc);
		1076	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1077	sprintf(fullname, "R_w&m = %.2f" , rate);
		1078	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1079	sprintf(fullname, "R_sum = %.2f" , sumrate);
		1080	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1081	sprintf(fullname, "R_w&m/sum = %.2f" , rate/sumrate);
		1082	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1083
		1084
		1085	leg[legi]->Draw();
		1086
		1087	if(draw_cuts) {
		1088	sprintf(hname, "htdcdiff_cut");
		1089	hp1dcut = DrTH1F(dir, hname, "");
		1090	hp1dcut->Rebin(rebin_tdcs);
		1091	hp1dcut->SetLineColor(kRed);
		1092	hp1dcut->DrawClone("SAME");
		1093
		1094	sprintf(hname, "htdcdiff_cut_2");
		1095	hp1dcut = DrTH1F(dir, hname, "");
		1096	hp1dcut->Rebin(rebin_tdcs);
		1097	hp1dcut->SetLineColor(kGreen);
		1098	hp1dcut->DrawClone("SAME");
		1099	}
		1100
		1101	sprintf(fullname, "gif/f_%s.gif", fname); c[cc]->SaveAs(fullname);
		1102	sprintf(fullname, "eps/f_%s.eps", fname); c[cc]->SaveAs(fullname);
		1103	cout<< "maximum x value = " << hp1d->GetXaxis()->GetBinCenter( hp1d->GetMaximumBin())<<endl;
		1104	if(batch_q) gSystem->Exit(1);
		1105	}
		1106
		1107	// -------------------------------------------------------------------
		1108
		1109	if( strchr(plopt, 'd') != NULL ) {
		1110
		1111	DrSetDrawStyle(0.05);
		1112
		1113	gStyle->SetOptStat("e");
		1114	//~ gStyle->SetOptStat(0);
		1115	gStyle->SetOptFit(1);
		1116	//~ gStyle->SetOptFit(0);
		1117
		1118	gStyle->SetPadRightMargin(0.03);
		1119	gStyle->SetPadTopMargin(0.05);
		1120	gStyle->SetPadBottomMargin(0.12);
		1121	gStyle->SetStatFontSize(0.05);
		1122	gStyle->SetStatX(0.97);
		1123	gStyle->SetStatY(0.95);
		1124
		1125	c[++cc] = new TCanvas("d", "d", 0, 0, 720, 400);
		1126	//~ c[cc]->Divide(2,2);
		1127
		1128	int draw_cuts = 0;
		1129	int rebin_tdcs = 1;
		1130
		1131	double drawmin = -10;
		1132	double drawmax = +10;
		1133	//~ double fitmin = inmin;
		1134	//~ double fitmax = inmax;
		1135	double fitmin = -10;
		1136	double fitmax = 10;
		1137	double fitcenter = 0.0;
		1138	double HalfMax, MinHalfMax, MaxHalfMax, FWHM, FWHMc;
		1139
		1140	int ccccdi = 0;
		1141	int ich = 0;
		1142
		1143	TF1 *f_fit;
		1144	TF1 *f_integrate;
		1145
		1146	TF1 *fc = new TF1("fc", "pol0(0)");
		1147
		1148	TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)"); fcg->SetNpx(300);
		1149	fcg->SetParName(0,"Const"); fcg->SetParName(1,"G_Const");
		1150	fcg->SetParName(2,"G_Mean" ); fcg->SetParName(3,"G_Sigma" );
		1151
		1152	TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)"); fcgg->SetNpx(300);
		1153	fcgg->SetParName(0,"Const");
		1154	fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
		1155	fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
		1156
		1157	TF1 *fcggg = new TF1("fcggg", "pol0(0)+gaus(1)+gaus(4)+gaus(7)"); fcggg->SetNpx(300);
		1158	fcggg->SetParName(0,"Const");
		1159	fcggg->SetParName(1,"G1_Const"); fcggg->SetParName(2,"G1_Mean" ); fcggg->SetParName(3,"G1_Sigma" );
		1160	fcggg->SetParName(4,"G2_Const"); fcggg->SetParName(5,"G2_Mean" ); fcggg->SetParName(6,"G2_Sigma" );
		1161	fcggg->SetParName(7,"G3_Const"); fcggg->SetParName(8,"G3_Mean" ); fcggg->SetParName(9,"G3_Sigma" );
		1162
		1163	(c[cc]->cd(++ccccdi))->SetLogy(0);
		1164	//~ sprintf(hname, "htdcdiff");
		1165	sprintf(hname, "htdcdiff_cut");
		1166	hp1d = DrTH1F(dir, hname, "");
		1167	//~ hp1d->Rebin(rebin_tdcs);
		1168	hp1d->SetLineColor(kBlack);
		1169	hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
		1170	hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*1.1);
		1171	//~ hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*2);
		1172	sprintf(fullname, ";Coincidence Time [ns];Counts");
		1173	hp1d->SetTitle(fullname);
		1174
		1175	//kobayashi
		1176	//~ fitmin = hp1d->GetXaxis()->GetBinCenter(hp1d->GetMaximumBin()) - hp1d->GetRMS()/15;
		1177	//~ fitmax = hp1d->GetXaxis()->GetBinCenter(hp1d->GetMaximumBin()) + hp1d->GetRMS()/15;
		1178
		1179	//~ if( strcmp(fname, "run000")==0 ) {
		1180	//~ f_fit = fcg;
		1181	//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		1182	//~ hp1d->GetMaximum()*0.9, 0, 0.05);
		1183	//~ } else if( strcmp(fname, "run061")==0 \|\|
		1184	//~ strcmp(fname, "run062")==0 \|\|
		1185	//~ strcmp(fname, "run063")==0 \|\|
		1186	//~ strcmp(fname, "run064")==0 \|\|
		1187	//~ strcmp(fname, "run065")==0) {
		1188	//~ f_fit = fcgg;
		1189	//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		1190	//~ hp1d->GetMaximum()*0.5, 0, 0.05,
		1191	//~ hp1d->GetMaximum()*0.4, 0, 0.5);
		1192	//~ } else {
		1193	//~ f_fit = fcgg;
		1194	//~ f_fit->SetParameters(hp1d->GetMaximum()*0.25,
		1195	//~ hp1d->GetMaximum()*0.7, 0, 0.4,
		1196	//~ hp1d->GetMaximum()*0.05, 0, 0.05);
		1197	//~ f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
		1198	//~ f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
		1199	//~ f_fit->SetParLimits(2, fitmin,fitmax);
		1200	//~ f_fit->SetParLimits(3, 0.2, 1);
		1201	//~ f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
		1202	//~ f_fit->SetParLimits(5, fitmin,fitmax);
		1203	//~ if( strcmp(fname, "run059")==0 \|\|
		1204	//~ strcmp(fname, "run060")==0)
		1205	//~ f_fit->SetParLimits(6, 0.02, 0.3);
		1206	//~ else
		1207	//~ f_fit->SetParLimits(6, 0.02, 0.2);
		1208	//~
		1209	//~ }
		1210
		1211	if( strcmp(fitopt, "cgg")==0 ) {
		1212	f_fit = fcgg;
		1213	f_fit->SetParameters(hp1d->GetMaximum()*0.1,
		1214	hp1d->GetMaximum()*0.4, 0, 0.9,
		1215	hp1d->GetMaximum()*0.5, 0, 0.05);
		1216	f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
		1217	f_fit->SetParLimits(1, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
		1218	f_fit->SetParLimits(2, fitmin,fitmax);
		1219	f_fit->SetParLimits(3, 0.05, 1);
		1220	f_fit->SetParLimits(4, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
		1221	f_fit->SetParLimits(5, fitmin,fitmax);
		1222	f_fit->SetParLimits(6, 0.01, 1);
		1223	} else if(strcmp(fitopt, "cggg")==0) {
		1224	f_fit = fcggg;
		1225	//original
		1226	//~ f_fit->SetParameters(hp1d->GetMaximum()*0.01,
		1227	//~ hp1d->GetMaximum()*0.4, 0, 0.4,
		1228	//~ hp1d->GetMaximum()*0.4, 0, 0.2,
		1229	//~ hp1d->GetMaximum()*0.1, 0, 0.05);
		1230	//kobayashi
		1231	double init_mean = hp1d->GetXaxis()->GetBinCenter(hp1d->GetMaximumBin());
		1232	f_fit->SetParameters(hp1d->GetMaximum()*0.0001,
		1233	hp1d->GetMaximum()*0.01, init_mean , 0.04,
		1234	hp1d->GetMaximum()*0.01, init_mean + 0.2, 0.03,
		1235	hp1d->GetMaximum()*0.01, init_mean - 0.2, 0.03);
		1236	f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
		1237	f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
		1238	f_fit->SetParLimits(2, fitmin,fitmax);
		1239	f_fit->SetParLimits(3, 0.001, 1);
		1240	f_fit->SetParLimits(4, hp1d->GetMaximum()*0.1, hp1d->GetMaximum());
		1241	f_fit->SetParLimits(5, fitmin,fitmax);
		1242	f_fit->SetParLimits(6, 0.005, 3);
		1243	f_fit->SetParLimits(7, hp1d->GetMaximum()*0.1, hp1d->GetMaximum());
		1244	f_fit->SetParLimits(8, fitmin,fitmax);
		1245	f_fit->SetParLimits(9, 0.02, 5);
		1246	} else if(strcmp(fitopt, "gg")==0) {
		1247	f_fit = fgg;
		1248	f_fit->SetParameters(hp1d->GetMaximum()*0.4, 0, 0.9,
		1249	hp1d->GetMaximum()*0.5, 0, 0.05);
		1250	f_fit->SetParLimits(0, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
		1251	f_fit->SetParLimits(1, fitmin,fitmax);
		1252	f_fit->SetParLimits(2, 0.05, 1);
		1253	f_fit->SetParLimits(3, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
		1254	f_fit->SetParLimits(4, fitmin,fitmax);
		1255	f_fit->SetParLimits(5, 0.01, 1);
		1256	}
		1257
		1258	f_fit->SetRange(fitmin,fitmax);
		1259	f_fit->SetLineWidth(1); f_fit->SetLineStyle(9);
		1260	hp1d->Fit(f_fit, "QWW", "", fitmin,fitmax);
		1261	//~ hp1d->Fit(f_fit, "QWW", "", -0.5,+0.5);
		1262
		1263	//~ fg->SetRange(fitmin,fitmax); fg->SetLineWidth(1.0);
		1264	//~ fg->SetParameters(f_fit->GetParameter(0+1), f_fit->GetParameter(1+1), f_fit->GetParameter(2+1)); fg->SetLineColor(kMagenta); fg->DrawCopy("LSAME");
		1265	//~ fg->SetParameters(f_fit->GetParameter(0+4), f_fit->GetParameter(1+4), f_fit->GetParameter(2+4)); fg->SetLineColor(kGreen); fg->DrawCopy("LSAME");
		1266	//~ fg->SetParameters(f_fit->GetParameter(0+7), f_fit->GetParameter(1+7), f_fit->GetParameter(2+7)); fg->SetLineColor(kYellow); fg->DrawCopy("LSAME");
		1267
		1268
		1269	HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
		1270	MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
		1271	MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
		1272	FWHM = MaxHalfMax - MinHalfMax;
		1273
		1274	HalfMax = (f_fit->GetMaximum(fitmin, fitmax) + f_fit->GetParameter(0))/2.0;
		1275	MinHalfMax = f_fit->GetX(HalfMax, fitmin, f_fit->GetMaximumX(fitmin, fitmax));
		1276	MaxHalfMax = f_fit->GetX(HalfMax, f_fit->GetMaximumX(fitmin, fitmax), fitmax);
		1277	FWHMc = MaxHalfMax - MinHalfMax;
		1278
		1279	//~ printf("FWHM = %lf ns \| FWHMc = %lf ns\n", FWHM, FWHMc);
		1280	printf("(d) FWHM (peak only) = %.3lf ns\n", FWHMc);
		1281
		1282	f_integrate = fgg;
		1283	for(int i=0; i<6; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
		1284	//~ for(int i=0; i<9; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
		1285	//~ fgg->SetRange(drawmin, drawmax); fgg->SetNpx(300); fgg->SetLineColor(kBlue); fgg->SetLineWidth(1); fgg->DrawClone("LSAME");
		1286	double Icherenkovs = f_integrate->Integral(drawmin, drawmax)/hp1d->GetBinWidth(1);
		1287
		1288	double Iall = f_fit->Integral(fitmin, fitmax)/hp1d->GetBinWidth(1);
		1289	double Rsn = Icherenkovs/(Iall-Icherenkovs);
		1290
		1291	//~ double Ientr = hp1d->GetEntries()*hp1d->GetBinWidth(1);
		1292	double Ientr = hp1d->GetEntries();
		1293	double Rsne = Icherenkovs/(Ientr-Icherenkovs);
		1294
		1295	//~ double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1)*hp1d->GetBinWidth(1);
		1296	double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1);
		1297	double Rsnh = Icherenkovs/(Ihist-Icherenkovs);
		1298
		1299	double Icherenkovs_10ns = f_integrate->Integral(-10, 10)/hp1d->GetBinWidth(1);
		1300	fc->SetParameter(0, f_fit->GetParameter(0));
		1301	double Inoise_10ns = fc->Integral(-10, 10)/hp1d->GetBinWidth(1);
		1302	double R10ns = Icherenkovs_10ns/Inoise_10ns;
		1303
		1304	double Icherenkovs_4ns = f_integrate->Integral(-4, 4)/hp1d->GetBinWidth(1);
		1305	fc->SetParameter(0, f_fit->GetParameter(0));
		1306	double Inoise_4ns = fc->Integral(-4, 4)/hp1d->GetBinWidth(1);
		1307	double R4ns = Icherenkovs_4ns/Inoise_4ns;
		1308
		1309	double Icherenkovs_2ns = f_integrate->Integral(-2, 2)/hp1d->GetBinWidth(1);
		1310	fc->SetParameter(0, f_fit->GetParameter(0));
		1311	double Inoise_2ns = fc->Integral(-2, 2)/hp1d->GetBinWidth(1);
		1312	double R2ns = Icherenkovs_2ns/Inoise_2ns;
		1313
		1314	//~ printf("ICherenkov = %lf \| IAll = %lf \| Ientr = %lf \| Ihist = %lf\n", Icherenkovs, Iall, Ientr, Ihist);
		1315	printf("Icherenkovs_4ns = %.0lf\n", Icherenkovs_4ns);
		1316	//~ printf("TF1 -> S/N = %lf \|\| TH1F(e) -> S/N = %lf \|\| TH1F(i) -> S/N = %lf\n", Rsn, Rsne, Rsnh);
		1317	//~ printf("S/N (All) = %lf\n", Rsnh);
		1318	printf("S/N (10ns) = %.2lf\n", R10ns);
		1319	printf("S/N (4 ns) = %.2lf\n", R4ns);
		1320	printf("S/N (2 ns) = %.2lf\n", R2ns);
		1321
		1322	//~ #define USE_NOISE_FILE
		1323
		1324	#ifdef USE_NOISE_FILE
		1325	//get ROOT file with histograms
		1326	char fnameroot[1024];
		1327	TFile * rootfile2;
		1328	TDirectory *dir2;
		1329
		1330	//~ sprintf(fname, "run054");
		1331	sprintf(fname, "run_309");
		1332	sprintf(fnameroot, "root/%s.root", fname);
		1333	rootfile2 = (TFile *) gROOT->FindObject(fname);
		1334	if(rootfile2==NULL) rootfile2 = new TFile(fnameroot);
		1335	if(rootfile2==NULL) {
		1336	printf("Cannot open root file 2 %s!!!\n",fnameroot);
		1337	return;
		1338	}
		1339	dir2 = (TDirectory*) rootfile2;
		1340
		1341	sprintf(hname, "htdcdiff");
		1342	hp1dcut = DrTH1F(dir2, hname, "");
		1343	hp1dcut->SetLineColor(40);
		1344	hp1dcut->SetLineStyle(2);
		1345	hp1dcut->DrawClone("SAME");
		1346
		1347	leg[++legi] = new TLegend(0.125,0.5,0.425,0.925);
		1348	leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(12);
		1349
		1350	leg[legi]->AddEntry((TH1F*)hp1d->Clone(), "With source", "L");
		1351	leg[legi]->AddEntry((TH1F*)hp1dcut->Clone(), "W/o source", "L");
		1352	#else
		1353	leg[++legi] = new TLegend(0.125,0.625,0.425,0.925);
		1354	leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(12);
		1355	#endif
		1356
		1357	sprintf(fullname, "FWHM = %.0lf ps", 1000*FWHMc);
		1358	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "L");
		1359	sprintf(fullname, "S/N (#pm10ns) = %.1lf", R10ns);
		1360	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1361	sprintf(fullname, "S/N (#pm 4ns) = %.1lf", R4ns);
		1362	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1363	sprintf(fullname, "S/N (#pm 2ns) = %.1lf", R2ns);
		1364	leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
		1365	leg[legi]->Draw();
		1366
		1367	if(draw_cuts) {
		1368	sprintf(hname, "htdcdiff_cut");
		1369	hp1dcut = DrTH1F(dir, hname, "");
		1370	hp1dcut->Rebin(rebin_tdcs);
		1371	hp1dcut->SetLineColor(kRed);
		1372	hp1dcut->DrawClone("SAME");
		1373
		1374	sprintf(hname, "htdcdiff_cut_2");
		1375	hp1dcut = DrTH1F(dir, hname, "");
		1376	hp1dcut->Rebin(rebin_tdcs);
		1377	hp1dcut->SetLineColor(kGreen);
		1378	hp1dcut->DrawClone("SAME");
		1379	}
		1380
		1381	sprintf(fullname, "gif/d_%s.gif", fname); c[cc]->SaveAs(fullname);
		1382	sprintf(fullname, "eps/d_%s.eps", fname); c[cc]->SaveAs(fullname);
		1383
		1384	if(batch_q) gSystem->Exit(1);
		1385	}
		1386
		1387	// -------------------------------------------------------------------
		1388
		1389	if( strchr(plopt, 'c') != NULL ) {
		1390
		1391	//~ gStyle->SetOptStat(1111);
		1392	//~ gStyle->SetOptFit(0);
		1393
		1394	c[++cc] = new TCanvas("c", "c", 0, 0, 900, 700);
		1395	//~ c[cc]->Divide(2,2);
		1396
		1397	int ccccdi = 0;
		1398
		1399	(c[cc]->cd(++ccccdi))->SetLogz(1);
		1400	sprintf(hname, "htdccor");
		1401	hp2d = DrTH2F(dir, hname, "");
		1402	//~ hp2d->GetYaxis()->SetRangeUser(0,1500);
		1403	//~ sprintf(fullname, "%s;QDC Reference 1275 keV; QDC Reference Coincidence", fname);
		1404	//~ hp2d->SetTitle(fullname);
		1405	hp2d->DrawClone("COLZ");
		1406
		1407	//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
		1408	//~ sprintf(hname, "hcor%d",0);
		1409	//~ hp2d = DrTH2F(dir, hname, "");
		1410	//~ hp2d->DrawClone("COLZ");
		1411	//~
		1412	//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
		1413	//~ sprintf(hname, "hcor%d",1);
		1414	//~ hp2d = DrTH2F(dir, hname, "");
		1415	//~ hp2d->DrawClone("COLZ");
		1416	//~
		1417	//~
		1418	//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
		1419	//~ sprintf(hname, "hdiffcor_%d_%d",0,1);
		1420	//~ hp2d = DrTH2F(dir, hname, "");
		1421	//~ hp2d->DrawClone("COLZ");
		1422	//~
		1423	//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
		1424	//~ sprintf(hname, "hdiffcor_%d_%d",1,0);
		1425	//~ hp2d = DrTH2F(dir, hname, "");
		1426	//~ hp2d->DrawClone("COLZ");
		1427
		1428	sprintf(fullname, "gif/c_%s.gif", fname); c[cc]->SaveAs(fullname);
		1429	}
		1430
		1431	// -------------------------------------------------------------------------------
		1432	// -------------------------------------------------------------------
		1433
		1434	if( strchr(plopt, 't') != NULL ) {
		1435
		1436	gStyle->SetOptStat(1111);
		1437	//~ gStyle->SetOptFit(0);
		1438
		1439	c[++cc] = new TCanvas("c", "c", 0, 0, 700, 400);
		1440	//~ c[cc]->Divide(2,2);
		1441
		1442	double drawmin = -1;
		1443	double drawmax = 2;
		1444	double fitmin = -1;
		1445	double fitmax = 0.75;
		1446	double fitcenter = 0.0;
		1447	double HalfMax, MinHalfMax, MaxHalfMax, FWHM;
		1448
		1449	int ccccdi = 0;
		1450
		1451	//~ (c[cc]->cd(++ccccdi))->SetLogy(1);
		1452	sprintf(hname, "hctdc%d", 0);
		1453	hp1d = DrTH1F(dir, hname, "");
		1454	hp1d->SetTitle("Vov = 1.5V; Time [ns];");
		1455	hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
		1456
		1457	f_fit = fgg;
		1458	f_fit->SetParameters(hp1d->GetMaximum()*0.9, 0, 0.1,
		1459	hp1d->GetMaximum()*0.1, 0, 0.5);
		1460	//~ f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
		1461	//~ f_fit->SetParLimits(2, fitmin,fitmax);
		1462	//~ f_fit->SetParLimits(3, 0.2, 1);
		1463	//~ f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
		1464	//~ f_fit->SetParLimits(5, fitmin,fitmax);
		1465
		1466	f_fit->SetRange(fitmin,fitmax);
		1467	f_fit->SetLineWidth(1);
		1468	hp1d->Fit(f_fit, "Q", "", fitmin,fitmax);
		1469
		1470
		1471	HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
		1472	MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
		1473	MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
		1474	FWHM = MaxHalfMax - MinHalfMax;
		1475	printf("FWHM = %.3lf ns\n", FWHM);
		1476
		1477	fg->SetRange(-2,5); fg->SetLineWidth(1.0);
		1478	fg->SetParameters(f_fit->GetParameter(0), f_fit->GetParameter(1), f_fit->GetParameter(2));
		1479	fg->SetLineColor(kGreen); fg->DrawClone("LSAME");
		1480	fg->SetParameters(f_fit->GetParameter(0+3), f_fit->GetParameter(1+3), f_fit->GetParameter(2+3));
		1481	fg->SetLineColor(kMagenta); fg->DrawClone("LSAME");
		1482
		1483	sprintf(fullname, "gif/t_%s.gif", fname); c[cc]->SaveAs(fullname);
		1484	sprintf(fullname, "eps/t_%s.eps", fname); c[cc]->SaveAs(fullname);
		1485	}
		1486
		1487	// -------------------------------------------------------------------------------
		1488	// -------------------------------------------------------------------
		1489
		1490	if( strchr(plopt, 'e') != NULL ) {
		1491
		1492
		1493	DrSetDrawStyle(0.05);
		1494
		1495	gStyle->SetOptStat("e");
		1496	gStyle->SetOptFit(0);
		1497
		1498	gStyle->SetPadRightMargin(0.03);
		1499	gStyle->SetPadTopMargin(0.05);
		1500	gStyle->SetPadBottomMargin(0.12);
		1501	gStyle->SetStatFontSize(0.05);
		1502	gStyle->SetStatX(0.97);
		1503	gStyle->SetStatY(0.95);
		1504
		1505	c[++cc] = new TCanvas("e_ref", "e_ref", 0, 0, 720, 400);
		1506
		1507
		1508	int spectra_rebin = 0;
		1509	int yaxis_max = 0;
		1510	//~ int yaxis_max = 1200;
		1511	//~ int xaxis_max = 2500;
		1512	int xaxis_max = 3000;
		1513
		1514	//~ double fitcenter = 1190.;
		1515	//~ double fitw = 70.;
		1516	//~ double Cher_min = 10;
		1517	//~ double Cher_max = 16;
		1518	double Cher_min = -120;
		1519	double Cher_max = -110;
		1520	//~ double const_min = 20;
		1521	//~ double const_max = 80;
		1522	double const_min = -250;
		1523	double const_max = -200;
		1524	int ccccdi = 0;
		1525
		1526	if( strcmp(fname, "eff003")==0) {fitcenter=760.; fitw=100.;}
		1527	if( strcmp(fname, "eff004")==0) {fitcenter=1220.; fitw=150.; spectra_rebin=8; yaxis_max=3300;}
		1528	if( strcmp(fname, "eff005")==0) {fitcenter=1140.; fitw=60.;}
		1529	if( strcmp(fname, "eff006")==0) {fitcenter=2050.; fitw=60.; xaxis_max=4000;}
		1530
		1531	if( strcmp(fname, "eff012")==0) {Cher_min=13; Cher_max=15;}
		1532	if( strcmp(fname, "eff013")==0) {Cher_min=12; Cher_max=14;}
		1533	if( strcmp(fname, "eff014")==0) {Cher_min=12; Cher_max=14;}
		1534	if( strcmp(fname, "eff015")==0) {Cher_min=12; Cher_max=14;}
		1535	if( strcmp(fname, "eff016")==0) {Cher_min=12; Cher_max=14;}
		1536	if( strcmp(fname, "eff017")==0) {Cher_min=12; Cher_max=14;}
		1537	if( strcmp(fname, "eff018")==0) {Cher_min=12; Cher_max=14;}
		1538
		1539	sprintf(hname, "hadc%d", 0);
		1540	hp1d = DrTH1F(dir, hname, "");
		1541	hp1d->SetTitle("; Charge [200 fC/bin];Counts");
		1542	if(spectra_rebin) hp1d->Rebin(spectra_rebin);
		1543	hp1d->SetLineColor(kBlack);
		1544	hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
		1545	if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0,yaxis_max);
		1546	hp1d->Draw();
		1547
		1548	//~ printf("Photofraction ~ %lf\n", hp1d->Integral(1300,1700)/hp1d->Integral(0,1300));
		1549
		1550	f_fit = fg;
		1551	f_fit->SetParameters(hp1d->GetMaximum()*0.1, fitcenter, fitw);
		1552
		1553	f_fit->SetRange(fitcenter-fitw, fitcenter+fitw);
		1554	f_fit->SetLineWidth(3);
		1555	hp1d->Fit(f_fit, "0Q", "", fitcenter-fitw, fitcenter+fitw);
		1556	f_fit->DrawClone("LSAME");
		1557
		1558	double PP_1sigma_min = f_fit->GetParameter(1)-f_fit->GetParameter(2);
		1559	double PP_1sigma_max = f_fit->GetParameter(1)+f_fit->GetParameter(2);
		1560	printf("PPcenter, PPwidth = %.0lf, %.0lf\n", f_fit->GetParameter(1), f_fit->GetParameter(2));
		1561	printf("PP+/-1sigma = [%.0lf, %.0lf]\n", PP_1sigma_min, PP_1sigma_max);
		1562	printf("PP resolution = %.1lf\%% (FWHM = %.1lf\%%)\n", 100f_fit->GetParameter(2)/f_fit->GetParameter(1), 235f_fit->GetParameter(2)/f_fit->GetParameter(1));
		1563
		1564	// -------- 1275keV Compton estimation --------
		1565
		1566	// COMPTON EDGE
		1567	//~ TF1 fcompt = new TF1("fcompt", "[2](1 - ConstantStep(((x-[0])/[1]) - (1275.(21275./511./(1+21275./511.)))))(2 + (((x-[0])/[1])/1275.)^2/((1275./511.)^2(1 - (((x-[0])/[1])/1275.))^2) + (((x-[0])/[1])/1275.)((((x-[0])/[1])/1275.) - 2/(1275./511.))/(1 - (((x-[0])/[1])/1275.)))/(1275./511.)^2");
		1568	//~ TF1 *fcomhi = new TF1("fcomhi", "gaus(0)+fcompt");
		1569	//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
		1570	//~ double comhi_range_low = 1670.;
		1571	//~ //double comhi_range_hi = 2130.;
		1572	//~ double comhi_range_hi = 2070.;
		1573	//~ fcomhi->SetParameters(2000., 1710., 50., 100, 2, 200);
		1574	//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
		1575	//~ hp1d->Fit(fcomhi, "0QR");
		1576	//~ //fcomhi->DrawClone("LSAME");
		1577	//~
		1578	//~ fcompt->SetNpx(300); fcompt->SetLineWidth(1); fcompt->SetLineColor(kBlue); fcompt->SetLineStyle(9);
		1579	//~ fcompt->SetParameter(0, fcomhi->GetParameter(3));
		1580	//~ fcompt->SetParameter(1, fcomhi->GetParameter(4));
		1581	//~ fcompt->SetParameter(2, fcomhi->GetParameter(5));
		1582	//~ fcompt->SetRange(0, comhi_range_hi);
		1583	//~ fcompt->DrawClone("LSAME");
		1584	//~
		1585	//~ double IComHi = fcompt->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
		1586	//~ printf("IComHi = %.0lf\n", IComHi);
		1587
		1588	// G LINEAR G
		1589	//~ TF1 *fcomhi = new TF1("fcomhi", "gaus(0)+pol0(3)+gaus(4)");
		1590	//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(kBlue);
		1591	//~ double comhi_range_low = 1700.;
		1592	//~ double comhi_range_hi = 2150.;
		1593	//~ fcomhi->SetParameters(2000., 1710., 50.,
		1594	//~ 150.,
		1595	//~ 200., 2150., 100.);
		1596	//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
		1597	//~ hp1d->Fit(fcomhi, "0QR");
		1598	//~ fcomhi->DrawClone("LSAME");
		1599	//~
		1600	//~ TF1 *fcom_lin = new TF1("fcom_lin", "pol0(0)");
		1601	//~ fcom_lin->SetNpx(300); fcom_lin->SetLineWidth(1); fcom_lin->SetLineColor(6);
		1602	//~ fcom_lin->SetParameter(0, fcomhi->GetParameter(3));
		1603	//~ fcom_lin->SetRange(0, 4000);
		1604	//~ fcom_lin->DrawClone("LSAME");
		1605	//~
		1606	//~ double IComHi = fcom_lin->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
		1607	// G LINEAR EXP
		1608	TF1 fcomhi = new TF1("fcomhi", "gaus(0) + [3] + [4]TMath::Exp((x-[5])/[6])");
		1609	fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
		1610	//double comhi_range_low = 1250.;
		1611	double comhi_range_low = fitcenter - 1.5*fitw;
		1612	//double comhi_range_hi = 2100.;
		1613	//~ fcomhi->SetParameters(2000., fitcenter, fitw,
		1614	//~ 170.,
		1615	//~ 85., 2000., 100.);
		1616	fcomhi->SetParameters(800., fitcenter, fitw,
		1617	40.,
		1618	5., comhi_range_hi, 100.);
		1619
		1620	//~ fcomhi->SetParLimits(2, 50., 2000.);
		1621	fcomhi->SetRange(comhi_range_low, comhi_range_hi);
		1622	hp1d->Fit(fcomhi, "0QR");
		1623	//~ fcomhi->SetRange(0, 4000);
		1624	fcomhi->DrawClone("LSAME");
		1625
		1626	TF1 fcompt = new TF1("fcompt", "[0] + [1]TMath::Exp((x-[2])/[3])");
		1627	fcompt->SetNpx(300); fcompt->SetLineWidth(1); fcompt->SetLineColor(kBlue); fcompt->SetLineStyle(9);
		1628	fcompt->SetParameter(0, fcomhi->GetParameter(3));
		1629	fcompt->SetParameter(1, fcomhi->GetParameter(4));
		1630	fcompt->SetParameter(2, fcomhi->GetParameter(5));
		1631	fcompt->SetParameter(3, fcomhi->GetParameter(6));
		1632	//~ fcompt->SetRange(1000, comhi_range_hi);
		1633	fcompt->SetRange(400, comhi_range_hi);
		1634	fcompt->DrawClone("LSAME");
		1635	//~
		1636	double IComHi = fcompt->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
		1637	printf("IComHi = %.0lf\n", IComHi);
		1638	//~
		1639	//~ // EXP
		1640	//~ TF1 fcomhi = new TF1("fcomhi", "[0] + [1]TMath::Exp((x-[2])/[3])");
		1641	//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
		1642	//~ double comhi_range_low = 1800.;
		1643	//~ double comhi_range_hi = 2000.;
		1644	//~ fcomhi->SetParameters(150., 1., 2000.,1.);
		1645	//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
		1646	//~ hp1d->Fit(fcomhi, "R");
		1647	//~ fcomhi->SetRange(0, 4000);
		1648	//~ fcomhi->DrawClone("LSAME");
		1649	//~
		1650	//~ double IComHi = fcomhi->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
		1651	//~ printf("IComHi = %.0lf\n", IComHi);
		1652	// LIN
		1653	//~ TF1 *fcomhi = new TF1("fcomhi", "pol0");
		1654	//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
		1655	//~ double comhi_range_low = 1800.;
		1656	//~ double comhi_range_hi = 1900.;
		1657	//~ fcomhi->SetParameter(0,150);
		1658	//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
		1659	//~ hp1d->Fit(fcomhi, "QR");
		1660	//~ fcomhi->SetRange(0, 4000);
		1661	//~ fcomhi->DrawClone("LSAME");
		1662	//~
		1663	//~ double IComHi = fcomhi->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
		1664	//~ printf("IComHi = %.0lf\n", IComHi);
		1665
		1666	// -------- 1275keV Compton estimation --------
		1667
		1668	sprintf(hname, "hadc_cut%d", 0);
		1669	hp1dcut = DrTH1F(dir, hname, "");
		1670	if(spectra_rebin) hp1dcut->Rebin(spectra_rebin);
		1671	hp1dcut->SetLineColor(kRed);
		1672
		1673	hp1dcut->SetLineColor(40);
		1674	hp1dcut->SetLineStyle(2);
		1675	hp1dcut->DrawClone("SAME");
		1676	double N511 = hp1dcut->GetEntries();
		1677	//~ double N511 = hp1dcut->Integral(PP_1sigma_min-10, PP_1sigma_max+10);
		1678
		1679	leg[++legi] = new TLegend(0.124302, 0.65, 0.6, 0.925134);
		1680	leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(22);
		1681	leg[legi]->AddEntry((TH1F*)hp1d->Clone(), "Measurement", "L");
		1682	leg[legi]->AddEntry((TH1F*)hp1dcut->Clone(), "#pm1 #sigma Photopeak Cut", "L");
		1683	leg[legi]->AddEntry((TF1*)fcompt->Clone(), "1275 keV Compton Estimate", "L");
		1684
		1685	//~ leg[legi]->Draw();
		1686
		1687	//~ sprintf(hname, "hadc_cut_2%d", 0);
		1688	//~ hp1d = DrTH1F(dir, hname, "");
		1689	//~ hp1d->SetLineColor(kGreen);
		1690	//~ hp1d->DrawClone("SAME");
		1691
		1692
		1693	sprintf(fullname, "gif/e_%s.gif", fname); c[cc]->SaveAs(fullname);
		1694	sprintf(fullname, "eps/e_%s.eps", fname); c[cc]->SaveAs(fullname);
		1695
		1696
		1697
		1698	c[++cc] = new TCanvas("e_sipm", "e_sipm", 0, 470, 720, 300);
		1699	int ch_Sprttype = 8;
		1700
		1701	(c[cc]->cd(1))->SetLogy(1);
		1702	sprintf(hname, "htdc%d", ch_Sprttype);
		1703	hp1d = DrTH1F(dir, hname, "");
		1704	hp1d->SetTitle("; Time [ns];Counts");
		1705	//~ hp1d->Rebin(4);
		1706	hp1d->SetLineColor(kBlack);
		1707	//~ hp1d->GetXaxis()->SetRangeUser(0,80);//original
		1708	hp1d->GetXaxis()->SetRangeUser(-400,50);
		1709	hp1d->GetYaxis()->SetRangeUser(0.1,hp1d->GetMaximum()*2);
		1710	hp1d->DrawClone();
		1711
		1712	sprintf(hname, "htdc_cut%d", ch_Sprttype);
		1713	hp1d = DrTH1F(dir, hname, "");
		1714	//~ hp1d->Rebin(4);
		1715	hp1d->SetLineColor(kRed);
		1716
		1717	TF1 *fc = new TF1("fc", "pol0"); fc->SetNpx(300);
		1718	fc->SetParameter(0,1);
		1719	//~ fc->SetParameter(1,-0.1);
		1720	fc->SetRange(const_min,const_max);
		1721	hp1d->Fit(fc, "0QWWR");
		1722
		1723	hp1d->DrawClone("SAME");
		1724	fc->SetLineColor(kBlue);
		1725	fc->DrawClone("LSAME");
		1726
		1727	double Ntdc = hp1d->Integral(hp1d->FindBin(Cher_min), hp1d->FindBin(Cher_max));
		1728	double Nbckg = fc->Integral(Cher_min,Cher_max)/hp1d->GetBinWidth(1);
		1729	printf("fc par(0) = %lf \| Nbckg = %.0lf\n", fc->GetParameter(0), Nbckg);
		1730
		1731
		1732	//~ sprintf(hname, "htdc_cut_2%d", 1);
		1733	//~ hp1d = DrTH1F(dir, hname, "");
		1734	//~ hp1d->SetLineColor(kGreen);
		1735	//~ hp1d->DrawClone("SAME");
		1736
		1737	printf("N511 = %.0lf \| Ntdc = %.0lf \| R = %.2lf \%% \| R_1 = %.2lf \%% \| R_2 = %.2lf \%%\n", N511, Ntdc, 100Ntdc/N511, 100(Ntdc-Nbckg)/N511, 100*(Ntdc-Nbckg)/(N511-IComHi));
		1738
		1739	sprintf(fullname, "gif/e2_%s.gif", fname); c[cc]->SaveAs(fullname);
		1740	sprintf(fullname, "eps/e2_%s.eps", fname); c[cc]->SaveAs(fullname);
		1741
		1742	int kobayashi =1;
		1743	if(kobayashi){
		1744	double Nnohit = hp1d->Integral(hp1d->FindBin(-715), hp1d->FindBin(-685));
		1745
		1746	printf("N511 = %.0lf \| Nnohit = %.0lf \| R = %.2lf \%% \| R_1 = %.2lf \%% \| R_2 = %.2lf \%%\n", N511, Nnohit, 100(N511-Nnohit)/N511, 100((N511-Nnohit)-Nbckg)/N511, 100*((N511-Nnohit)-Nbckg)/(N511-IComHi));
		1747
		1748	}
		1749
		1750	//~ (c[cc]->cd(++ccccdi))->SetLogy(1);
		1751	//~ sprintf(hname, "hadc%d", 1);
		1752	//~ hp1d = DrTH1F(dir, hname, "");
		1753	//~ hp1d->SetLineColor(kBlack);
		1754	//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
		1755	//~ hp1d->DrawClone();
		1756	//~
		1757	//~ sprintf(hname, "hadc_cut%d", 1);
		1758	//~ hp1d = DrTH1F(dir, hname, "");
		1759	//~ hp1d->SetLineColor(kRed);
		1760	//~ hp1d->DrawClone("SAME");
		1761	//~
		1762	//~ sprintf(hname, "hadc_cut_2%d", 1);
		1763	//~ hp1d = DrTH1F(dir, hname, "");
		1764	//~ hp1d->SetLineColor(kGreen);
		1765	//~ hp1d->DrawClone("SAME");
		1766
		1767	//~ sprintf(fullname, "gif/e_%s.gif", fname); c[cc]->SaveAs(fullname);
		1768	//~ sprintf(fullname, "eps/e_%s.eps", fname); c[cc]->SaveAs(fullname);
		1769	}
		1770	// -------------------------------------------------------------------
		1771
		1772	if( strchr(plopt, 'r') != NULL ) {
		1773
		1774	//~ gStyle->SetOptStat(1111);
		1775	gStyle->SetOptStat(0);
		1776	gStyle->SetOptFit(0);
		1777
		1778	c[++cc] = new TCanvas("e", "e", 0, 0, 1200, 450);
		1779	c[cc]->Divide(3,1);
		1780
		1781	int spectra_rebin = 0;
		1782	int yaxis_max = 0;
		1783	int xaxis_max = 3000;
		1784
		1785	double fitcenter = 700.;
		1786	double fitw = 60.;
		1787	double Cher_min = 4;
		1788	double Cher_max = 16;
		1789	double const_min = Cher_max;
		1790	double const_max = 80;
		1791	int ccccdi = 0;
		1792
		1793	if( strcmp(fname, "eff003")==0) {fitcenter=760.; fitw=100.;}
		1794	if( strcmp(fname, "eff004")==0) {fitcenter=1220.; fitw=150.; spectra_rebin=8; yaxis_max=3300;}
		1795	if( strcmp(fname, "eff005")==0) {fitcenter=1140.; fitw=60.;}
		1796	if( strcmp(fname, "eff006")==0) {fitcenter=2050.; fitw=60.; xaxis_max=4000;}
		1797	if( strcmp(fname, "eff007")==0) {fitcenter=710.; fitw=50.;}
		1798	if( strcmp(fname, "eff009")==0) {fitcenter=900.; fitw=60.; xaxis_max=2000;}
		1799
		1800	(c[cc]->cd(++ccccdi))->SetLogy(1);
		1801	sprintf(hname, "hadc%d", 0);
		1802	hp1d = DrTH1F(dir, hname, "");
		1803	hp1d->SetTitle("Reference Coincidence; Charge [200 fC/bin];");
		1804	if(spectra_rebin) hp1d->Rebin(spectra_rebin);
		1805	hp1d->SetLineColor(kBlack);
		1806	hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
		1807	if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0.5,yaxis_max);
		1808
		1809	f_fit = fg;
		1810	f_fit->SetParameters(hp1d->GetMaximum()*0.1, fitcenter, fitw);
		1811
		1812	f_fit->SetRange(fitcenter-fitw, fitcenter+fitw);
		1813	f_fit->SetLineWidth(1);
		1814	hp1d->Fit(f_fit, "Q", "", fitcenter-fitw, fitcenter+fitw);
		1815
		1816	printf("PP+/-1sigma = [%.0lf, %.0lf]\n", f_fit->GetParameter(1)-f_fit->GetParameter(2), f_fit->GetParameter(1)+f_fit->GetParameter(2));
		1817	printf("PP resolution = %.1lf\%% (FWHM = %.1lf\%%)\n", 100f_fit->GetParameter(2)/f_fit->GetParameter(1), 235f_fit->GetParameter(2)/f_fit->GetParameter(1));
		1818
		1819	sprintf(hname, "hadc_cut%d", 0);
		1820	hp1d = DrTH1F(dir, hname, "");
		1821	if(spectra_rebin) hp1d->Rebin(spectra_rebin);
		1822	hp1d->SetLineColor(kRed);
		1823	hp1d->DrawClone("SAME");
		1824	double N511 = hp1d->GetEntries();
		1825
		1826	sprintf(hname, "hadc_cut_2%d", 0);
		1827	hp1d = DrTH1F(dir, hname, "");
		1828	hp1d->SetLineColor(kGreen);
		1829	hp1d->DrawClone("SAME");
		1830
		1831	//-----------------------------------------------------------------------
		1832
		1833	(c[cc]->cd(++ccccdi))->SetLogy(1);
		1834	sprintf(hname, "hadc%d", 2);
		1835	hp1d = DrTH1F(dir, hname, "");
		1836	hp1d->SetTitle("Reference 1275 keV; Charge [200 fC/bin];");
		1837	//~ if(spectra_rebin) hp1d->Rebin(spectra_rebin);
		1838	hp1d->SetLineColor(kBlack);
		1839	//~ hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
		1840	//~ if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0,yaxis_max);
		1841	hp1d->DrawClone();
		1842
		1843	sprintf(hname, "hadc_cut%d", 2);
		1844	hp1d = DrTH1F(dir, hname, "");
		1845	if(spectra_rebin) hp1d->Rebin(spectra_rebin);
		1846	hp1d->SetLineColor(kRed);
		1847	hp1d->DrawClone("SAME");
		1848
		1849	sprintf(hname, "hadc_cut_2%d", 2);
		1850	hp1d = DrTH1F(dir, hname, "");
		1851	hp1d->SetLineColor(kGreen);
		1852	hp1d->DrawClone("SAME");
		1853
		1854	//-----------------------------------------------------------------------
		1855
		1856	(c[cc]->cd(++ccccdi))->SetLogy(1);
		1857	sprintf(hname, "htdc%d", 1);
		1858	hp1d = DrTH1F(dir, hname, "");
		1859	hp1d->SetTitle("MPPC Cherenkov; Time [ns];");
		1860	//~ hp1d->Rebin(4);
		1861	hp1d->SetLineColor(kBlack);
		1862	hp1d->GetXaxis()->SetRangeUser(0,80);
		1863	hp1d->GetYaxis()->SetRangeUser(0.1,hp1d->GetMaximum()*2);
		1864	hp1d->DrawClone();
		1865
		1866	sprintf(hname, "htdc_cut%d", 1);
		1867	hp1d = DrTH1F(dir, hname, "");
		1868	//~ hp1d->Rebin(4);
		1869	hp1d->SetLineColor(kRed);
		1870
		1871	TF1 *fc = new TF1("fc", "pol0"); fc->SetNpx(300);
		1872	fc->SetParameter(0,1);
		1873	fc->SetRange(const_min,const_max);
		1874	hp1d->Fit(fc, "0QWWR");
		1875
		1876	hp1d->DrawClone("SAME");
		1877	fc->SetLineColor(kBlue);
		1878	fc->DrawClone("LSAME");
		1879
		1880	double Ntdc = hp1d->Integral(hp1d->FindBin(Cher_min), hp1d->FindBin(Cher_max));
		1881	double Nbckg = fc->Integral(Cher_min,Cher_max)/hp1d->GetBinWidth(1);
		1882	printf("fc par(0) = %lf \| Nbckg = %.0lf\n", fc->GetParameter(0), Nbckg);
		1883
		1884
		1885	sprintf(hname, "htdc_cut_2%d", 1);
		1886	hp1d = DrTH1F(dir, hname, "");
		1887	hp1d->SetLineColor(kGreen);
		1888	hp1d->DrawClone("SAME");
		1889
		1890	printf("N511 = %.0lf \| Ntdc = %.0lf \| R = %lf\n", N511, Ntdc, Ntdc/N511);
		1891
		1892	printf("R_1 = %lf\n", (Ntdc-Nbckg)/N511);
		1893
		1894	//~ (c[cc]->cd(++ccccdi))->SetLogy(1);
		1895	//~ sprintf(hname, "hadc%d", 1);
		1896	//~ hp1d = DrTH1F(dir, hname, "");
		1897	//~ hp1d->SetLineColor(kBlack);
		1898	//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
		1899	//~ hp1d->DrawClone();
		1900	//~
		1901	//~ sprintf(hname, "hadc_cut%d", 1);
		1902	//~ hp1d = DrTH1F(dir, hname, "");
		1903	//~ hp1d->SetLineColor(kRed);
		1904	//~ hp1d->DrawClone("SAME");
		1905	//~
		1906	//~ sprintf(hname, "hadc_cut_2%d", 1);
		1907	//~ hp1d = DrTH1F(dir, hname, "");
		1908	//~ hp1d->SetLineColor(kGreen);
		1909	//~ hp1d->DrawClone("SAME");
		1910
		1911	sprintf(fullname, "gif/e_%s.gif", fname); c[cc]->SaveAs(fullname);
		1912	sprintf(fullname, "eps/e_%s.eps", fname); c[cc]->SaveAs(fullname);
		1913	}
		1914
		1915	// -------------------------------------------------------------------------------
		1916	// -------------------------------------------------------------------
		1917
		1918	if( strchr(plopt, 'p') != NULL ) {
		1919
		1920
		1921	DrSetDrawStyle(0.05);
		1922
		1923	gStyle->SetOptStat("e");
		1924	gStyle->SetOptFit(0);
		1925
		1926	gStyle->SetPadRightMargin(0.03);
		1927	gStyle->SetPadTopMargin(0.05);
		1928	gStyle->SetPadBottomMargin(0.12);
		1929	gStyle->SetStatFontSize(0.05);
		1930	gStyle->SetStatX(0.97);
		1931	gStyle->SetStatY(0.95);
		1932
		1933	c[++cc] = new TCanvas("Print", "Print", 0, 0, 720, 400);
		1934
		1935
		1936	int spectra_rebin = 0;
		1937	int yaxis_max = 1000;
		1938	int xaxis_max = 1200;
		1939
		1940
		1941	double Cher_min = 10;
		1942	double Cher_max = 16;
		1943	double const_min = 20;
		1944	double const_max = 80;
		1945	int ccccdi = 0;
		1946
		1947	sprintf(hname, "hadc%d", 0);
		1948	hp1d = DrTH1F(dir, hname, "");
		1949	hp1d->SetTitle("; Charge [200 fC/bin];Counts");
		1950	if(spectra_rebin) hp1d->Rebin(spectra_rebin);
		1951	hp1d->SetLineColor(kBlack);
		1952	hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
		1953	if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0,yaxis_max);
		1954	hp1d->Draw();
		1955
		1956	sprintf(hname, "hadc_cut%d", 0);
		1957	hp1dcut = DrTH1F(dir, hname, "");
		1958	if(spectra_rebin) hp1dcut->Rebin(spectra_rebin);
		1959	hp1dcut->SetLineColor(kRed);
		1960	hp1dcut->Draw("SAME");
		1961
		1962	sprintf(fullname, "gif/p_%s.gif", fname); c[cc]->SaveAs(fullname);
		1963	sprintf(fullname, "eps/p_%s.eps", fname); c[cc]->SaveAs(fullname);
		1964	}
		1965
		1966	// -------------------------------------------------------------------
		1967
		1968	if( strchr(plopt, 'b') != NULL ) {
		1969	gStyle->SetOptStat(0);
		1970	gStyle->SetOptFit(0);
		1971
		1972	c[++cc] = new TCanvas("b", "b", 0, 0, 1200, 700);
		1973	c[cc]->Divide(2,2);
		1974
		1975	int ccccdi = 0;
		1976	int ich;
		1977	/*
		1978	// ------------------------------------------------------------
		1979	ich = 0;
		1980	(c[cc]->cd(++ccccdi))->SetLogy(1);
		1981	sprintf(hname, "htdc%d", ich);
		1982	hp1d = DrTH1F(dir, hname, "");
		1983	hp1d->SetTitle(hname);
		1984	//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
		1985	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		1986	hp1d->SetLineColor(kBlack);
		1987	hp1d->DrawClone();
		1988
		1989	// ------------------------------------------------------------
		1990	ich = 1;
		1991	(c[cc]->cd(++ccccdi))->SetLogy(1);
		1992	sprintf(hname, "htdc%d", ich);
		1993	hp1d = DrTH1F(dir, hname, "");
		1994	hp1d->SetTitle(hname);
		1995	//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
		1996	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		1997	hp1d->SetLineColor(kBlack);
		1998	hp1d->DrawClone();
		1999	*/
		2000	// -------------------------- qdc 0 ----------------------------------
		2001	ich = 0;
		2002	(c[cc]->cd(++ccccdi))->SetLogy(1);
		2003	/*
		2004	sprintf(hname, "hadc%d", ich);
		2005	hp1d = DrTH1F(dir, hname, "");
		2006	hp1d->SetTitle("SiPM;Charge [200 fC/bin]");
		2007	hp1d->GetXaxis()->SetRangeUser(0,1000);
		2008	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		2009	hp1d->SetLineColor(kBlack);
		2010	hp1d->DrawClone();
		2011	*/
		2012
		2013	TLegend *leg0 = new TLegend(0.57,0.66,0.88,0.88);
		2014	leg0->SetFillColor(0);
		2015	leg0->SetBorderSize(1);
		2016	leg0->SetTextSize(0.05);
		2017	leg0->SetTextAlign(22);
		2018
		2019	sprintf(hname, "hadc%d", ich);
		2020	TH1F* hADC = DrTH1F(dir, hname, "");
		2021
		2022	hADC->SetTitle("SiPM;Charge [200 fC/bin]");
		2023	hADC->GetXaxis()->SetRangeUser(0,1000);
		2024	hADC->SetLineColor(kBlack);
		2025	hADC->DrawClone();
		2026
		2027	int dbgprint=1;
		2028	const int nPeaks=3;
		2029	double pSigma=10;
		2030	double pThr=0.005;
		2031
		2032	TSpectrum *sADC = new TSpectrum(nPeaks,0.3);
		2033	sADC->Search(hADC, pSigma,"",pThr);
		2034	int foundNPeaks = sADC->GetNPeaks();
		2035	float *foundPeakX = sADC->GetPositionX();
		2036	float *foundPeakY = sADC->GetPositionY();
		2037
		2038	if(1 < foundNPeaks) {
		2039	//void Sort(Int_t n1, const Float_t a, Int_t index, Bool_t down)
		2040	int sortedPeakIndex[nPeaks];
		2041	TMath::Sort(foundNPeaks, foundPeakX, sortedPeakIndex, 0);
		2042	if(dbgprint) {
		2043	printf("foundNPeaks = %d\n",foundNPeaks);
		2044	for(int i=0; i<foundNPeaks; i++) printf("Peak[%d] = %f\n",i, foundPeakX[sortedPeakIndex[i]]);
		2045	}
		2046	float peakMidpointX[nPeaks+1];
		2047	float peakMidpointY[nPeaks+1];
		2048
		2049	peakMidpointX[0] = 0; peakMidpointY[0] = 1;
		2050	for(int i=1; i<foundNPeaks; i++) {
		2051	peakMidpointX[i] = (foundPeakX[sortedPeakIndex[i]] + foundPeakX[sortedPeakIndex[i-1]])/2.0;
		2052	peakMidpointY[i] = 1;
		2053
		2054	if(dbgprint) printf("peakMidpointX[%d] = %f\n",i, peakMidpointX[i]);
		2055
		2056	}
		2057	peakMidpointX[foundNPeaks] = 1120; peakMidpointY[foundNPeaks] = 1;
		2058
		2059	if(dbgprint) printf("peakMidpointX[%d] = %f\n",foundNPeaks, peakMidpointX[foundNPeaks]);
		2060
		2061	TPolyMarker *pmMidpoints = new TPolyMarker(nPeaks+1, peakMidpointX, peakMidpointY);
		2062	pmMidpoints->SetMarkerStyle(33); pmMidpoints->SetMarkerColor(9); pmMidpoints->SetMarkerSize(2);
		2063	pmMidpoints->Draw("SAME");
		2064
		2065	//~ fprintf(fp,"peakEvents:\n");
		2066
		2067	double peakEvents[nPeaks];
		2068	for(int i=0; i<foundNPeaks-1; i++) {
		2069	peakEvents[i] = hADC->Integral(hADC->FindBin(peakMidpointX[i]), hADC->FindBin(peakMidpointX[i+1]));
		2070
		2071	if(dbgprint) printf("I(%.1f, %.1f) = %.1lf\n", peakMidpointX[i], peakMidpointX[i+1], peakEvents[i]);
		2072	//~ fprintf(fp,"%lf\n", peakEvents[i]);
		2073
		2074	}
		2075	peakEvents[foundNPeaks-1] = hADC->Integral(hADC->FindBin(peakMidpointX[foundNPeaks-1]), hADC->FindBin(1130.));
		2076
		2077	if(dbgprint) printf("I(%.1f, %.1f) = %.1lf\n", peakMidpointX[foundNPeaks-1], 1130., peakEvents[foundNPeaks-1]);
		2078	//~ fprintf(fp,"%lf\n", peakEvents[foundNPeaks-1]);
		2079
		2080	double N0all = -TMath::Log(peakEvents[0]/hADC->GetEntries());
		2081	double N01 = peakEvents[1]/peakEvents[0];
		2082	double N02 = 2*peakEvents[2]/peakEvents[1];
		2083
		2084	if(dbgprint) printf("Poisson <N>: N0all -> %.3lf \| N01 -> %.3lf \| N12 -> %.3lf\n", N0all, N01, N02);
		2085
		2086
		2087
		2088	sprintf(sbuff, "<N> = %.3lf", N0all);
		2089	leg0->AddEntry((TH1F*)hADC->Clone(),sbuff, "");
		2090	sprintf(sbuff, "1 CPH = %.0lf", foundPeakX[sortedPeakIndex[1]]);
		2091	leg0->AddEntry((TH1F*)hADC->Clone(),sbuff, "");
		2092	leg0->Draw();
		2093	}
		2094
		2095	double N_adc_integral = hADC->Integral(hADC->FindBin(peakMidpointX[1]), hADC->FindBin(1000.));
		2096	double N_poisson = (1 - TMath::Exp(-N0all))*hADC->GetEntries();
		2097	printf("hADC->GetEntries() = %.0lf\n", hADC->GetEntries());
		2098	printf("N_adc_integral = %.0lf \| N_poisson = %.0lf\n", N_adc_integral, N_poisson);
		2099
		2100	sprintf(hname, "hadc_cut%d", ich);
		2101	hp1d = DrTH1F(dir, hname, "");
		2102	hp1d->SetLineColor(kRed);
		2103	hp1d->DrawClone("SAME");
		2104
		2105	// -------------------------- TDC 0 ----------------------------------
		2106	(c[cc]->cd(++ccccdi))->SetLogy(1);
		2107	sprintf(hname, "htdc%d", ich);
		2108	hp1d = DrTH1F(dir, hname, "");
		2109	hp1d->SetTitle("SiPM;Time [ns]");
		2110	hp1d->GetXaxis()->SetRangeUser(-5,5);
		2111	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		2112	hp1d->SetLineColor(kBlack);
		2113	hp1d->DrawClone();
		2114
		2115	printf("hTDC->GetEntries() = %.0lf\n", hp1d->GetEntries());
		2116	double N_tdc_integral = hp1d->Integral(hp1d->FindBin(-10), hp1d->FindBin(10));
		2117	printf("N_tdc_integral = %.0lf\n", N_tdc_integral);
		2118
		2119	// -----------------------------------------------------------
		2120
		2121	(c[cc]->cd(++ccccdi))->SetLogy(0);
		2122	sprintf(hname, "hcor%d", ich);
		2123	hp2d = DrTH2F(dir, hname, "");
		2124	hp2d->SetTitle(";Charge [200 fC/bin];cTime [ns]");
		2125	//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
		2126	//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		2127	//~ hp1d->SetLineColor(kBlack);
		2128	hp2d->DrawClone("COLZ");
		2129
		2130	// -------------------------- cTDC 0 ----------------------------------
		2131	(c[cc]->cd(++ccccdi))->SetLogy(0);
		2132	sprintf(hname, "hctdc%d", ich);
		2133	hp1d = DrTH1F(dir, hname, "");
		2134	//~ hp1d->SetTitle("SiPM;cTime [ns]");
		2135	//~ hp1d->GetXaxis()->SetRangeUser(-5,5);
		2136	//~ //hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		2137	hp1d->SetLineColor(kBlack);
		2138	//~ hp1d->DrawClone();
		2139
		2140	plotCorTDC(hp1d, "SiPM;cTime [ns]", -0.5, 0.45, -1., 2);
		2141
		2142	printf("hcTDC->GetEntries() = %.0lf\n", hp1d->GetEntries());
		2143	double N_ctdc_integral = hp1d->Integral(hp1d->FindBin(-2), hp1d->FindBin(2));
		2144	printf("N_ctdc_integral = %.0lf\n", N_ctdc_integral);
		2145
		2146	printf("N_ctdc_integral/N_poisson = %.3lf\n", N_ctdc_integral/N_poisson);
		2147	printf("N_adc_integral/N_poisson = %.3lf\n", N_adc_integral/N_poisson);
		2148	printf("N_ctdc_integral/N_adc_integral = %.3lf\n", N_ctdc_integral/N_adc_integral);
		2149
		2150	//~ sprintf(hname, "htdc_cut%d", ich);
		2151	//~ hp1d = DrTH1F(dir, hname, "");
		2152	//~ hp1d->SetLineColor(kRed);
		2153	//~ hp1d->DrawClone("SAME");
		2154
		2155	sprintf(fullname, "gif/b_%s.gif", fname); c[cc]->SaveAs(fullname);
		2156	}
		2157
		2158
		2159	// -------------------------------------------------------------------
		2160
		2161	if( strchr(plopt, 'g') != NULL ) {
		2162	gStyle->SetOptStat(0);
		2163	gStyle->SetOptFit(0);
		2164
		2165	c[++cc] = new TCanvas("g", "g", 640, 0, 640, 360);
		2166
		2167	int ccccdi = 0;
		2168	int ich = 15;
		2169
		2170	(c[cc]->cd(++ccccdi))->SetLogy(0);
		2171	sprintf(hname, "htdc%d", ich);
		2172	hp1d = DrTH1F(dir, hname, "");
		2173	//~ hp1d->SetTitle("SiPM;cTime [ns]");
		2174	//~ hp1d->GetXaxis()->SetRangeUser(-5,5);
		2175	//~ //hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
		2176	hp1d->SetLineColor(kBlack);
		2177
		2178	//~ hp1d->DrawClone();
		2179
		2180	sprintf(hname, "Trigger(ch.%d)-Trigger(ch.31);Time [ns]", ich);
		2181	plotCorTDC(hp1d, hname, -0.25, 0.25, -0.5, 0.5);
		2182
		2183	sprintf(fullname, "gif/g_%s.gif", fname); c[cc]->SaveAs(fullname);
		2184	}
		2185	// -------------------------------------------------------------------------------
		2186	}

Subversion Repositories f9daq

(root)/cvi/apps/l2d2_easyroc/plots2.cpp – Rev 291