#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "TGraphErrors.h"
#define COROUT
#define SHOW_SLICES
#define NCH 20
#define USED_CH 1
#define TDC_BIN (25./1000.) //1 TDC bin in ns
//~ #define FIRST_CH 1
double g_corpars[NCH][3];
int mapii[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19 };
int mapij[] = { 1, 0,16,16,16,16,16,16,16,16, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 };
// root -l "cor.cpp(\"run045\",64,1,20,\"s\")"
void cor_laser(char *fname="cscan3", const int nslicex=64, int slmin=1, const int projmin=5, char *fite="s")
{
double fitw=1.0;
char skip='n';
char fullname[256];
//get ROOT file with histograms
char fnameroot[1024];
TFile * rootfile;
TDirectory *dir;
sprintf(fnameroot, "root/%s.root", fname);
rootfile = (TFile *) gROOT->FindObject(fname);
if(rootfile==NULL) rootfile = new TFile(fnameroot);
if(rootfile==NULL) {
printf("Cannot open root file %s!!!\n",fnameroot);
return;
}
dir = (TDirectory*) rootfile;
// set draw style
gStyle->SetOptStat("ne");
gStyle->SetPalette(1, 0);
gStyle->SetPaperSize(TStyle::kA4);
gStyle->SetStatBorderSize(1);
gStyle->SetFrameBorderMode(0);
gStyle->SetFrameFillColor(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetCanvasColor(0);
gStyle->SetStatColor(0);
gStyle->SetOptFit(11);
gStyle->SetPadRightMargin(0.1);
char hname[256];
TH1F *hp1d; TH2F *hp2d;
int ch=0;
int nbinsx, deltax, bin0, bin1;
double binc;
TH1F *hprojx[nslicex];
TF1 *f_fit = 0;
TF1 *fg = new TF1("fg", "gaus(0)");
fg->SetParNames("Constant","Mean","Sigma");
TF1 *fgg = new TF1("fgg", "gaus(0)+gaus(3)");
fgg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2");
TF1 *fggg = new TF1("fgg", "gaus(0)+gaus(3)+gaus(6)");
fggg->SetParNames("Constant","Mean","Sigma","Constant2","Mean2","Sigma2","Constant3","Mean3","Sigma3");
TF1 *fitsqrt = new TF1("fitsqrt","[0]+[1]/TMath::Sqrt(x-[2])");
TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)");
fcg->SetParNames("C", "Constant","Mean","Sigma");
TGraphErrors *gcor[NCH];
double fitcenter, fitampl;
double fitsigmas[2];
int fitminsigmai;
double fiterror, fitmean;
#ifdef COROUT
FILE *fp;
//open file to save parameters
sprintf(fullname, "data/%s_cor.txt", fname);
if( (fp=fopen(fullname, "wt")) == NULL )
printf("Cannot create parameter file %s !!!\n", fullname);
else
printf("Created parameter file %s\n", fullname);
#endif
TCanvas *c1 = new TCanvas("Cor.Fits", "Cor.Fits", 0, 0, 1300, 700);
c1->Divide(2,1);
#ifdef SHOW_SLICES
TCanvas *c0 = new TCanvas("Slices", "Slices", 0, 0, 1300, 700);
c0->Divide(5,5);
#endif
const int c0pads = 25;
if(nslicex==0) {
#ifdef COROUT
for(int j=0;j<NCH; j++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);
gSystem->Exit(1);
#endif
} else {
#ifdef FIRST_CH
ch=FIRST_CH;
#else
for(ch=0;ch<USED_CH; ch++)
#endif
{
gcor[ch] = new TGraphErrors(nslicex);
sprintf(hname,"TDC Corection ch%d;QDC;TDC [ps]",ch);
gcor[ch]->SetTitle(hname);
//~ sprintf(hname, "hdiffcor_%d_%d",mapii[ch],mapij[ch]); hp2d = (TH2F *) dir->Get(hname);
sprintf(hname, "hcor%d",ch); hp2d = (TH2F *) dir->Get(hname);
nbinsx=hp2d->GetNbinsX();
deltax=nbinsx/nslicex;
c0->Clear();c0->Divide(8,4);
for(int i=slmin; i<nslicex-1; i++) {
bin0=i*deltax+1;
bin1=(i+1)*deltax+1;
binc=(hp2d->GetXaxis()->GetBinCenter(bin0)+hp2d->GetXaxis()->GetBinCenter(bin1))/2.0;
sprintf(hname, "hprojx%d",i);
hprojx[i]=(TH1F*)hp2d->ProjectionY(hname, bin0, bin1);
//~ (hprojx[i]->GetXaxis())->SetRangeUser(-1, 2);
if( hprojx[i]->GetMaximum() > projmin ) {
fitcenter = hprojx[i]->GetBinCenter(hprojx[i]->GetMaximumBin());
fitampl = hprojx[i]->GetMaximum();
//------/
if(i<1024) {
//~ if( (ch<8) && (i<3) && (fitcenter < 0.4) ) fitcenter = 0.7;
//~ f_fit = fcg;
//~ f_fit->SetParameters(fitampl*0.1, fitampl*0.9, fitcenter, 0.1);
f_fit = fgg;
f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1, fitampl*0.1, fitcenter, 0.5);
//~
//~ f_fit = fg;
//~ f_fit->SetParameters(fitampl*0.9, fitcenter, 0.1);
//~ f_fit->SetParLimits(0, 0, fitampl*3);
#ifdef SHOW_SLICES
c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
hprojx[i]->Fit(f_fit,"QL", "", fitcenter-1.0*fitw, fitcenter+0.5*fitw);
//~ hprojx[i]->Fit(f_fit,"QL", "");
//c0->Update();
//if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
//if( (skip=='q') || (skip=='Q') ) return;
//hprojx[i]->DrawClone();
//fgg->SetRange(-1,2);
//fgg->DrawClone("LSAME");
#else
//~ hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
hprojx[i]->Fit(f_fit,"0QL", "");
#endif
//~ fitsigmas[0] = TMath::Abs(f_fit->GetParameter(2));
//~ fitsigmas[1] = TMath::Abs(f_fit->GetParameter(5));
//~ fitminsigmai = TMath::LocMin(2,fitsigmas);
//~ fitmean = f_fit->GetParameter(1 + fitminsigmai*3);
//~ fiterror = TMath::Abs(f_fit->GetParameter(2 + fitminsigmai*3));
//~ printf("fitminsigmai = %d {%lf, %lf}\n", fitminsigmai, fitsigmas[0], fitsigmas[1]);
//~
//~ fitmean = f_fit->GetParameter(1);
//~ fiterror = TMath::Abs(f_fit->GetParameter(2));
//~
fitmean = f_fit->GetMaximumX(fitcenter-1.0*fitw, fitcenter+0.5*fitw);
fiterror = TMath::Abs(f_fit->GetParameter(2));
//------/
} else {
fgg->SetParameters(fitampl, fitcenter, 0.05,
fitampl/15, fitcenter, 0.5);
fgg->SetParLimits(0, fitampl*0.3, fitampl*1.1);
fgg->SetParLimits(3, fitampl*0.01, fitampl*0.7);
fgg->SetParLimits(1, -1, 1);
fgg->SetParLimits(4, -1, 1);
fgg->SetParLimits(2, 0, 0.300);
fgg->SetParLimits(5, 0.3, 4.0);
#ifdef SHOW_SLICES
c0->cd(i%c0pads + 1); //c0->cd(i%4 + 1)->SetLogy();
hprojx[i]->Fit(fgg,"QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
//c0->Update();
//if( !((i+1)%9) && (skip!='y') ) {c0->Update(); printf("Skip?: "); scanf("%c",&skip);}
//if( (skip=='q') || (skip=='Q') ) return;
//hprojx[i]->DrawClone();
//fgg->SetRange(-1,2);
//fgg->DrawClone("LSAME");
#else
hprojx[i]>Fit(fgg,"0QL", "", fitcenter-1.0*fitw, fitcenter+2.0*fitw);
#endif
fitsigmas[0] = TMath::Abs(fgg->GetParameter(2));
fitsigmas[1] = TMath::Abs(fgg->GetParameter(5));
fitminsigmai = TMath::LocMin(2,fitsigmas);
//printf("fitminsigmai = %d {%.1lf, %.1lf}\n", fitminsigmai, fitsigmas[0]*1000, fitsigmas[1]*1000);
fitmean = fgg->GetParameter(1 + fitminsigmai*3);
fiterror = TMath::Abs(fgg->GetParameter(2 + fitminsigmai*3));
//------/
} //if(i<4)
printf(">>> ADC %3d - %3d -> %5.1lf ps\n", bin0, bin1, fiterror*1000.);
gcor[ch]->SetPoint(i, binc, fitmean);
// weight = slice fit sigma (default)
if( strcmp(fite, "s")==0 )
gcor[ch]->SetPointError(i, 0.0, fiterror);
// weight = slice maximum bin content ~ n of events
else if( strcmp(fite, "h")==0 )
gcor[ch]->SetPointError(i, 0.0, 10./sqrt(hprojx[i]->GetMaximum()));
// weight = both
else if( strcmp(fite, "b")==0 )
gcor[ch]->SetPointError(i, 0, 100.*fiterror/hprojx[i]->GetMaximum());
} //if( hprojx[i]->GetMaximum() > projmin )
}// for(int i=slmin; i<nslicex; i++)
#ifdef SHOW_SLICES
c0->Update();
sprintf(fullname, "gif/%s_ch_%d.gif", fname, ch); c0->SaveAs(fullname);
#endif
c1->cd(ch+1); gPad->SetLogz();
fitsqrt->SetLineWidth(2);
gStyle->SetOptFit();
(gcor[ch]->GetYaxis())->SetRangeUser((hp2d->GetYaxis())->GetXmin(), (hp2d->GetYaxis())->GetXmax());
//gcor[ch]->Draw("AP");
fitsqrt->SetParameters(0, 10, 50);
//fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),(gcor[ch]->GetXaxis())->GetXmax());
fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),3*256.);
//fitsqrt->SetParameters(-0.25, 3, 55);
//fitsqrt->SetParLimits(0,-4,8);
//~ fitsqrt->SetParLimits(1,0,20);
//~ fitsqrt->SetParLimits(2,0,200);
//~ #ifdef FIRST_CH
//~ if(ch==FIRST_CH) gcor[ch]->Fit(fitsqrt, "0Q","");
//~ #else
//~ if(ch==0) gcor[ch]->Fit(fitsqrt, "0Q","");
//~ #endif
gcor[ch]->Fit(fitsqrt, "0QR","");
//(gcor[ch]->GetYaxis())->UnZoom();
hp2d->SetTitle("TDCvsQDC; QDC; TDC [ns]");
//~ (hp2d->GetXaxis())->SetRangeUser(0.,4*256.);
//(hp2d->GetYaxis())->SetRangeUser(-1.,3.);
//~ (hp2d->GetYaxis())->SetRangeUser(-3.,3.);
hp2d->DrawCopy("COLZ");
fitsqrt->SetRange(TMath::Ceil(fitsqrt->GetParameter(2)),4*256.);
fitsqrt->DrawCopy("LSAME");
gcor[ch]->Draw("PSAME");
sprintf(fullname, "gif/%s_cor.gif", fname); c1->SaveAs(fullname);
printf("[ch.%d - ch.%d] %lf %lf\n", mapii[ch], mapij[ch], fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
#ifdef COROUT
fprintf(fp,"%lf %lf %lf\n", fitsqrt->GetParameter(0), fitsqrt->GetParameter(1), fitsqrt->GetParameter(2));
#endif
g_corpars[ch][0]=fitsqrt->GetParameter(0);
g_corpars[ch][1]=fitsqrt->GetParameter(1);
g_corpars[ch][2]=fitsqrt->GetParameter(2);
}
}
#ifdef COROUT
for(int i=USED_CH;i<NCH; i++) fprintf(fp,"%lf %lf %lf\n", 0.0,0.005,0.0);
if(fp) fclose(fp);
#endif
//rootfile->Close();
}