WebSVN - f9daq - Blame - Rev 291 - /cvi/apps/l2d2

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	#include "TCanvas.h"
		7	#include "TStyle.h"
		8	#include "TPad.h"
		9	#include "TF1.h"
		10	#include "TGraph.h"
		11	#include "TSpectrum.h"
		12	#include "TGraphErrors.h"
		13	#include "TSystem.h"
		14	#include "TMath.h"
		15
		16	#define COROUT
		17	#define SHOW_SLICES
		18	#define NCH 20
		19	#define USED_CH 2
		20	#define TDC_BIN (25./1000.) //1 TDC bin in ns
		21	//~ #define FIRST_CH 0
		22
		23	double g_corpars[NCH][3];
		24
		25	int mapii[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19 };
		26	int mapij[] = { 1, 0,16,16,16,16,16,16,16,16, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 };
		27
		28	// root -l "cor.cpp(\"run045\",64,1,20,\"s\")"
		29	void cor(char fname="cscan3", int range=1, const int nslicex=64, int slmin=1, const int projmin=5, char fite="s")
		30	{
		31	double fitw=1.0;
		32	char skip='n';
		33	char fullname[256];
		34
		35	//get ROOT file with histograms
		36	char fnameroot[1024];
		37	TFile * rootfile;
		38	TDirectory *dir;
		39
		40	sprintf(fnameroot, "root/%s.root", fname);
		41	rootfile = (TFile *) gROOT->FindObject(fname);
		42	if(rootfile==NULL) rootfile = new TFile(fnameroot);
		43	if(rootfile==NULL) {
		44	printf("Cannot open root file %s!!!\n",fnameroot);
		45	return;
		46	}
		47	dir = (TDirectory*) rootfile;
		48
		49
		50	// set draw style
		51	gStyle->SetOptStat("ne");
		52	gStyle->SetPalette(1, 0);
		53
		54	gStyle->SetPaperSize(TStyle::kA4);
		55	gStyle->SetStatBorderSize(1);
		56	gStyle->SetFrameBorderMode(0);
		57	gStyle->SetFrameFillColor(0);
		58	gStyle->SetCanvasBorderMode(0);
		59	gStyle->SetPadBorderMode(0);
		60	gStyle->SetPadColor(0);
		61	gStyle->SetCanvasColor(0);
		62	gStyle->SetStatColor(0);
		63	gStyle->SetOptFit(11);
		64
		65	gStyle->SetPadRightMargin(0.1);
		66
		67	char hname[256];
		68	TH1F hp1d; TH2F hp2d;
		69
		70	int ch=0;
		71	int nbinsx, deltax, bin0, bin1;
		72	double binc;
		73	TH1F *hprojx[nslicex];
		74
		75	TF1 *f_fit = 0;
		76	TF1 *fg = new TF1("fg", "gaus(0)");
		77	fg->SetParNames("Constant","Mean","Sigma");
		78	TF1 *fgg = new TF1("fgg", "gaus(0)+gaus(3)");
		79	fgg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2");
		80	TF1 *fggg = new TF1("fgg", "gaus(0)+gaus(3)+gaus(6)");
		81	fggg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2","Constant3","Mean3","Sigma3");
		82	TF1 *fitsqrt = new TF1("fitsqrt","[0]+[1]/TMath::Sqrt(x-[2])");
		83
		84	TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)");
		85	fcg->SetParNames("C", "Constant","Mean","Sigma");
		86
		87	TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)");
		88	fcgg->SetParName(0,"Const");
		89	fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
		90	fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
		91
		92	TGraphErrors *gcor[NCH];
		93	double fitcenter, fitampl;
		94
		95	double fitsigmas[2];
		96	int fitminsigmai;
		97	double fiterror, fitmean;
		98
		99	#ifdef COROUT
		100	FILE *fp;
		101	//open file to save parameters
		102	sprintf(fullname, "data/%s_cor.txt", fname);
		103	if( (fp=fopen(fullname, "wt")) == NULL )
		104	printf("Cannot create parameter file %s !!!\n", fullname);
		105	else
		106	printf("Created parameter file %s\n", fullname);
		107	#endif
		108
		109
		110	TCanvas *c1 = new TCanvas("Cor.Fits", "Cor.Fits", 0, 0, 1300, 450);
		111	c1->Divide(2,1);
		112
		113
		114	#ifdef SHOW_SLICES
		115	TCanvas *c0 = new TCanvas("Slices", "Slices", 0, 520, 1300, 470);
		116	c0->Divide(5,5);
		117	#endif
		118	const int c0pads = 25;
		119
		120
		121	if(nslicex==0) {
		122	#ifdef COROUT
		123	for(int j=0;j<NCH; j++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);
		124	gSystem->Exit(1);
		125	return;
		126	#endif
		127	} else {
		128	#ifdef FIRST_CH
		129	ch=FIRST_CH;
		130	#else
		131	for(ch=0;ch<USED_CH; ch++)
		132	#endif
		133	{
		134
		135	gcor[ch] = new TGraphErrors(nslicex);
		136	sprintf(hname,"TDC Corection ch%d;QDC;TDC [ps]",ch);
		137	gcor[ch]->SetTitle(hname);
		138	if(range == 0) {
		139	sprintf(hname, "hdiffcor_low_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);
		140	} else if(range == 2) {
		141	sprintf(hname, "hdiffcor_hi_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);
		142	} else {
		143	sprintf(hname, "hdiffcor_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);
		144	}
		145	//~ sprintf(hname, "hcor%d",ch); hp2d = (TH2F *) dir->Get(hname);
		146	nbinsx=hp2d->GetNbinsX();
		147	deltax=nbinsx/nslicex;
		148
		149	c0->Clear();c0->Divide(8,4);
		150	for(int i=slmin; i<nslicex; i++) {
		151	bin0=i*deltax+1;
		152	bin1=(i+1)*deltax+1;
		153	binc=(hp2d->GetXaxis()->GetBinCenter(bin0)+hp2d->GetXaxis()->GetBinCenter(bin1))/2.0;
		154
		155	sprintf(hname, "hprojx%d",i);
		156	hprojx[i]=(TH1F*)hp2d->ProjectionY(hname, bin0, bin1);
		157
		158	//~ (hprojx[i]->GetXaxis())->SetRangeUser(-1, 2);
		159	if( hprojx[i]->GetMaximum() > projmin ) {
		160
		161	fitcenter = hprojx[i]->GetBinCenter(hprojx[i]->GetMaximumBin());
		162	fitampl = hprojx[i]->GetMaximum();
		163
		164	//------/
		165	if(i<1024) {
		166
		167	//~ if( (ch<8) && (i<3) && (fitcenter < 0.4) ) fitcenter = 0.7;
		168
		169	f_fit = fcg;
		170	f_fit->SetParameters(fitampl0.1, fitampl0.9, fitcenter, 0.1);
		171	//~
		172	//~ f_fit = fcgg;
		173	//~ f_fit->SetParameters(fitampl*0.01,
		174	//~ fitampl*0.4, 0, 0.4,
		175	//~ fitampl*0.4, 0, 0.2);
		176
		177	//~ f_fit = fgg;
		178	//~ f_fit->SetParameters(fitampl0.9, fitcenter, 0.1, fitampl0.1, fitcenter, 0.5);
		179	//~
		180	//~ f_fit = fg;
		181	//~ f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1);
		182
		183	//~ f_fit->SetParLimits(0, 0, fitampl*3);
		184
		185	#ifdef SHOW_SLICES
		186	c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
		187	//~ hprojx[i]->Fit(f_fit,"QL", "", fitcenter-1.0fitw, fitcenter+0.5fitw);
		188	hprojx[i]->Fit(f_fit,"QL", "");
		189	//c0->Update();
		190	//if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
		191	//if( (skip=='q') \|\| (skip=='Q') ) return;
		192	//hprojx[i]->DrawClone();
		193	//fgg->SetRange(-1,2);
		194	//fgg->DrawClone("LSAME");
		195	#else
		196	//~ hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0fitw, fitcenter+2.0fitw);
		197	hprojx[i]->Fit(f_fit,"0QL", "");
		198	#endif
		199
		200	//~ fitsigmas[0] = TMath::Abs(f_fit->GetParameter(2));
		201	//~ fitsigmas[1] = TMath::Abs(f_fit->GetParameter(5));
		202	//~ fitminsigmai = TMath::LocMin(2,fitsigmas);
		203	//~ fitmean = f_fit->GetParameter(1 + fitminsigmai*3);
		204	//~ fiterror = TMath::Abs(f_fit->GetParameter(2 + fitminsigmai*3));
		205	//~ printf("fitminsigmai = %d {%lf, %lf}\n", fitminsigmai, fitsigmas[0], fitsigmas[1]);
		206	//~
		207	//~ fitmean = f_fit->GetParameter(1);
		208	//~ fiterror = TMath::Abs(f_fit->GetParameter(2));
		209	//~
		210
		211	// fcg
		212	fitmean = f_fit->GetMaximumX(fitcenter-1.0fitw, fitcenter+0.5fitw);
		213	fiterror = TMath::Abs(f_fit->GetParameter(3));
		214
		215	//~ // fcgg
		216	//~ fitmean = f_fit->GetMaximumX(fitcenter-1.0fitw, fitcenter+1.0fitw);
		217	//~ double HalfMax = (f_fit->GetMaximum(fitcenter-1.0fitw, fitcenter+1.0fitw) + f_fit->GetParameter(0))/2.0;
		218	//~ double MinHalfMax = f_fit->GetX(HalfMax, fitcenter-1.0*fitw, fitcenter);
		219	//~ double MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitcenter+1.0*fitw);
		220	//~ double FWHMc = MaxHalfMax - MinHalfMax;
		221	//~ fiterror = TMath::Abs(FWHMc);
		222	//~ //------/
		223	} else {
		224	fgg->SetParameters(fitampl, fitcenter, 0.05,
		225	fitampl/15, fitcenter, 0.5);
		226
		227	fgg->SetParLimits(0, fitampl0.3, fitampl1.1);
		228	fgg->SetParLimits(3, fitampl0.01, fitampl0.7);
		229
		230	fgg->SetParLimits(1, -1, 1);
		231	fgg->SetParLimits(4, -1, 1);
		232
		233	fgg->SetParLimits(2, 0, 0.300);
		234	fgg->SetParLimits(5, 0.3, 4.0);
		235
		236	#ifdef SHOW_SLICES
		237	c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
		238	hprojx[i]->Fit(fgg,"QL", "", fitcenter-1.0fitw, fitcenter+2.0fitw);
		239	//c0->Update();
		240	//if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
		241	//if( (skip=='q') \|\| (skip=='Q') ) return;
		242	//hprojx[i]->DrawClone();
		243	//fgg->SetRange(-1,2);
		244	//fgg->DrawClone("LSAME");
		245	#else
		246	hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0fitw, fitcenter+2.0fitw);
		247	#endif
		248
		249	fitsigmas[0] = TMath::Abs(fgg->GetParameter(2));
		250	fitsigmas[1] = TMath::Abs(fgg->GetParameter(5));
		251	fitminsigmai = TMath::LocMin(2,fitsigmas);
		252	//printf("fitminsigmai = %d {%.1lf, %.1lf}\n", fitminsigmai, fitsigmas[0]1000, fitsigmas[1]1000);
		253	fitmean = fgg->GetParameter(1 + fitminsigmai*3);
		254	fiterror = TMath::Abs(fgg->GetParameter(2 + fitminsigmai*3));
		255	//------/
		256	} //if(i<4)
		257
		258	printf(">>> ADC %3d - %3d -> %5.1lf ps\n", bin0, bin1, fiterror*1000.);
		259
		260	gcor[ch]->SetPoint(i, binc, fitmean);
		261	// weight = slice fit sigma (default)
		262	if( strcmp(fite, "s")==0 )
		263	gcor[ch]->SetPointError(i, 0.0, fiterror);
		264	// weight = slice maximum bin content ~ n of events
		265	else if( strcmp(fite, "h")==0 )
		266	gcor[ch]->SetPointError(i, 0.0, 2./sqrt(hprojx[i]->GetMaximum()));
		267	// weight = both
		268	else if( strcmp(fite, "b")==0 )
		269	gcor[ch]->SetPointError(i, 0, 100.*fiterror/hprojx[i]->GetMaximum());
		270	else if( strcmp(fite, "k")==0 )
		271	gcor[ch]->SetPointError(i, 0, 1./sqrt(hprojx[i]->GetMaximum()));
		272
		273	} //if( hprojx[i]->GetMaximum() > projmin )
		274	}// for(int i=slmin; i<nslicex; i++)
		275	#ifdef SHOW_SLICES
		276	c0->Update();
		277	sprintf(fullname, "gif/%s_ch_%d.gif", fname, ch); c0->SaveAs(fullname);
		278	#endif
		279
		280
		281	c1->cd(ch+1); gPad->SetLogz();
		282
		283	fitsqrt->SetLineWidth(2);
		284	gStyle->SetOptStat("e");
		285	gStyle->SetOptFit(11);
		286	(gcor[ch]->GetYaxis())->SetRangeUser((hp2d->GetYaxis())->GetXmin(), (hp2d->GetYaxis())->GetXmax());
		287	//~ (gcor[ch]->GetYaxis())->SetRangeUser(-2,3);
		288	//gcor[ch]->Draw("AP");
		289
		290	fitsqrt->SetParameters(0, 10, 50);
		291
		292	fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),(gcor[ch]->GetXaxis())->GetXmax());
		293	//~ fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),3*256.);
		294	//fitsqrt->SetParameters(-0.25, 3, 55);
		295	//fitsqrt->SetParLimits(0,-4,8);
		296	//~ fitsqrt->SetParLimits(1,0,20);
		297	//~ fitsqrt->SetParLimits(2,-100,2000);
		298	fitsqrt->SetParLimits(2,0,2000);
		299	//~ #ifdef FIRST_CH
		300	//~ if(ch==FIRST_CH) gcor[ch]->Fit(fitsqrt, "0Q","");
		301	//~ #else
		302	//~ if(ch==0) gcor[ch]->Fit(fitsqrt, "0Q","");
		303	//~ #endif
		304	gcor[ch]->Fit(fitsqrt, "0QR","");
		305	//(gcor[ch]->GetYaxis())->UnZoom();
		306	sprintf(fullname, "SiPM_%d; Charge [200 fC/bin]; TDC [ns]", ch+1);
		307	hp2d->SetTitle(fullname);
		308	//~ (hp2d->GetXaxis())->SetRangeUser(0.,4*256.);
		309	//~ (hp2d->GetXaxis())->SetRangeUser(0.,3*256.);
		310	//(hp2d->GetYaxis())->SetRangeUser(-1.,3.);
		311	//~ (hp2d->GetYaxis())->SetRangeUser(-2.5,2.5);
		312	hp2d->DrawCopy("COLZ");
		313	fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),8*256.);
		314	//~ fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),(gcor[ch]->GetXaxis())->GetXmax());
		315	fitsqrt->DrawCopy("LSAME");
		316	gcor[ch]->Draw("PSAME");
		317
		318	sprintf(fullname, "gif/%s_cor.gif", fname); c1->SaveAs(fullname);
		319
		320	printf("[ch.%d - ch.%d] %lf %lf\n", mapii[ch], mapij[ch], fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
		321	#ifdef COROUT
		322	fprintf(fp,"%lf %lf %lf\n", fitsqrt->GetParameter(0), fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
		323	#endif
		324	g_corpars[ch][0]=fitsqrt->GetParameter(0);
		325	g_corpars[ch][1]=fitsqrt->GetParameter(1);
		326	g_corpars[ch][2]=fitsqrt->GetParameter(2);
		327
		328	}
		329	}
		330
		331	#ifdef COROUT
		332	for(int i=USED_CH;i<NCH; i++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);
		333	if(fp) fclose(fp);
		334	#endif
		335
		336	//rootfile->Close();
		337	}

Subversion Repositories f9daq

(root)/cvi/apps/l2d2_easyroc/cor.cpp @ 352 - Rev 291