Go to most recent revision |
Blame |
Compare with Previous |
Last modification |
View Log
| RSS feed
#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
gROOT->Reset();
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "DrootHelper.cpp"
#define NCH 2
//-------------------------------------------------------------------------------------
void SetGS() {
const UInt_t Number = 2;
Double_t Red[Number] = {0.00, 1.00};
Double_t Green[Number] = {0.00, 1.00};
Double_t Blue[Number] = {0.00, 1.00};
Double_t Stops[Number] = {0.00, 1.00};
Int_t nb=50;
TColor::CreateGradientColorTable(Number, Stops, Red, Green, Blue, nb);
}
//-------------------------------------------------------------------------------------
Double_t ConstantStep(Double_t x)
{
if(x<=0) return 0.;
else return 1.;
}
// =======================================================================================================
void plotCorTDC(TH1F *hp1d, char *hTitle, double fitmin, double fitmax, double drawmin, double drawmax)
{
char sbuff[256];
TLegend *leg1 = new TLegend(0.57,0.76,0.88,0.88);
leg1->SetFillColor(0);
leg1->SetBorderSize(1);
leg1->SetTextSize(0.05);
leg1->SetTextAlign(22);
//gStyle->SetOptFit(111);
sprintf(sbuff, "%s;Time[ns];", hTitle);
hp1d->SetTitle(sbuff);
//~ hp1d->GetXaxis()->SetRangeUser(2*fitmin, 2*fitmax);
hp1d->GetXaxis()->SetRangeUser(drawmin, drawmax);
//~ hcTDC->DrawClone();
//~ TF1 *fs = new TF1("fs","pol0(0) + gaus(1)");
TF1 *fs = new TF1("fs","gaus(0) + ([3]*exp((x)/[4]) + [5]) * (0.5+0.5*TMath::Erf((x-[6])/[7]))");
fs->SetParName(0, "G_Const"); fs->SetParName(1, "G_Mean"); fs->SetParName(2, "G_Sigma");
fs->SetParName(3, "Exp_p0"); fs->SetParName(4, "Exp_p1"); fs->SetParName(5, "Exp_p2");
fs->SetParName(6, "Erf_p0"); fs->SetParName(7, "Erf_p1");
//~
TF1 *fg1 = new TF1("fg1", "gaus(0)"); fg1->SetNpx(400);
TF1 *fe1 = new TF1("fe1", "[0]*exp((x)/[1]) + [2]");
TF1 *ferf1 = new TF1("ferf1", "[0]*(0.5+0.5*TMath::Erf((x-[1])/[2]))");
TF1 *fe1erf = new TF1("fe1erf", "([0]*exp((x)/[1]) + [2]) * (0.5+0.5*TMath::Erf((x-[3])/[4]))"); fe1erf->SetNpx(400);
double fitampl = hp1d->GetMaximum();
double fitcenter = hp1d->GetBinCenter(hp1d->GetMaximumBin());
fs->SetNpx(400);
fs->SetParameters(fitampl, fitcenter, 0.01, 1e2, -0.1, 1e1, 0, 0.05);
// G + Er + Ex
//fs->SetParLimits( 1, -0.1, 0.1);
//fs->SetParLimits( 2, 0.018, 0.5);
fs->SetParLimits( 4, -0.3, 0);
//fs->SetParLimits( 6, -1e-1, 1e-1);
fs->SetParLimits( 7, 0.04, 0.6);
//~ // G
//~ fs->SetParameters(fitampl/100., fitampl, fitcenter, 0.01);
hp1d->Fit(fs,"Q","",fitmin, fitmax);
//~ printf("Chi / NDF = %lf/%lf\n", fs->GetChisquare(), fs->GetNDF());
fg1->SetParameters(fs->GetParameter(0), fs->GetParameter(1), fs->GetParameter(2));
//~ fe1->SetParameters(fs->GetParameter(3), fs->GetParameter(4), fs->GetParameter(5));
//~ ferf1->SetParameters(1e4, fs->GetParameter(6), fs->GetParameter(7));
fe1erf->SetParameters(fs->GetParameter(3), fs->GetParameter(4), fs->GetParameter(5), fs->GetParameter(6), fs->GetParameter(7));
fg1->SetLineColor(kRed); fg1->SetLineWidth(1); fg1->SetLineStyle(1); fg1->SetRange(fitmin, fitmax); fg1->DrawCopy("SAME");
//~ fe1->SetLineColor(kBlue); fe1->SetLineWidth(1); fe1->SetLineStyle(1); fe1->SetRange(fitmin, fitmax); fe1->DrawCopy("SAME");
//~ ferf1->SetLineColor(kYellow); ferf1->SetLineWidth(1);ferf1->SetLineStyle(1);ferf1->SetRange(fitmin, fitmax);ferf1->DrawCopy("SAME");
fe1erf->SetLineColor(kBlue); fe1erf->SetLineWidth(1);fe1erf->SetLineStyle(1);fe1erf->SetRange(fitmin, fitmax);fe1erf->DrawCopy("SAME");
double HalfMax = fitampl/2.0;
double MinHalfMax = fs->GetX(HalfMax, -0.5, fitcenter);
double MaxHalfMax = fs->GetX(HalfMax, fitcenter, 1.5);
double FWHM = MaxHalfMax - MinHalfMax;
double sigmaFWHM = fs->GetParameter(2);
printf("Corrected %s resolution FWHM = %.1lf [ps]\n", hTitle, FWHM*1000);
//~ fprintf(fp,"%lf ", FWHM*1000);
sprintf(sbuff, "FWHM = %.0lf ps", FWHM*1000);
leg1->AddEntry((TF1*)fs->Clone(),sbuff, "L");
leg1->Draw();
}
// =======================================================================================================
void plots2(char *fname="480", char *plopt="d", double fitcenter = 1190., double fitw = 70., double comhi_range_hi = 2100., int batch_q=0)
{
char sbuff[256];
char *fitopt="crystal1";
double inmin=-10;
double inmax=10;
double Ecut_lo = 450.;
double Ecut_hi = 625.;
int printeps = 0;
int printgif=-1, double fitlo=1.0, double fithi=1.0;
int ch=0;
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;
DrSetDrawStyle();
TCanvas *c[64];
int cc=-1;
char hname[256], htitle[256];
TH1F *hp1d, *hp1dcut; TH2F *hp2d;
TLegend *leg[64];
int legi = -1;
TH1F *legentry[64];
int legentryi = -1;
TGraph *graph[16];
int igraph=-1;
int colarr[20] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, kOrange, 10,11,12,13,14,15,16,17,18,19};
int legsty[] = {1,7,2,3,9,6,9,5,6};
TF1 *fg = new TF1("fg", "gaus"); fg->SetNpx(300);
TF1 *fgg = new TF1("fgg", "gaus(0)+gaus(3)"); fg->SetNpx(300);
fgg->SetParName(0,"G1_Const"); fgg->SetParName(3,"G2_Const");
fgg->SetParName(1,"G1_Mean" ); fgg->SetParName(4,"G2_Mean" );
fgg->SetParName(2,"G1_Sigma"); fgg->SetParName(5,"G2_Sigma");
TF1 *fg3 = new TF1("fg3", "gaus(0)+gaus(3)+gaus(6)"); fg3->SetNpx(300);
fg3->SetParName(0,"G1_Const"); fg3->SetParName(3,"G2_Const"); fg3->SetParName(6,"G3_Const");
fg3->SetParName(1,"G1_Mean" ); fg3->SetParName(4,"G2_Mean" ); fg3->SetParName(7,"G3_Mean" );
fg3->SetParName(2,"G1_Sigma"); fg3->SetParName(5,"G2_Sigma"); fg3->SetParName(8,"G3_Sigma");
TF1 *fg4 = new TF1("fg4", "gaus(0)+gaus(3)+gaus(6)+gaus(9)"); fg4->SetNpx(300);
fg4->SetParName(0,"G1_Const"); fg4->SetParName(3,"G2_Const"); fg4->SetParName(6,"G3_Const"); fg4->SetParName(9,"G4_Const");
fg4->SetParName(1,"G1_Mean" ); fg4->SetParName(4,"G2_Mean" ); fg4->SetParName(7,"G3_Mean" ); fg4->SetParName(10,"G4_Mean" );
fg4->SetParName(2,"G1_Sigma"); fg4->SetParName(5,"G2_Sigma"); fg4->SetParName(8,"G3_Sigma"); fg4->SetParName(11,"G4_Sigma");
TF1 *fg5 = new TF1("fg5", "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)"); fg5->SetNpx(300);
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");
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" );
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");
// -------------------------------------------------------------------
if( strchr(plopt, 'a') != NULL ) {
//~ gStyle->SetOptStat(1111);
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);
c[++cc] = new TCanvas("a", "a", 0, 0, 1200, 700);
c[cc]->Divide(2,2);
int draw_cuts = 1;
int rebin_tdcs = 0;
int ccccdi = 0;
int ich;
// -------------------------- qdc 0 ----------------------------------
ich = 0;
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("MPPC 1;Charge [200 fC/bin]");
hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "hadc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
//~ sprintf(hname, "hadc_cut_2%d", ich);
//~ hp1dcut = DrTH1F(dir, hname, "");
//~ hp1dcut->SetLineColor(kGreen);
//~ hp1dcut->DrawClone("SAME");
}
// -------------------------- TDC 0 ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("MPPC 1;Time [ns]");
if(rebin_tdcs) hp1d->Rebin(rebin_tdcs);
//~ hp1d->GetXaxis()->SetRangeUser(0,22);
hp1d->GetXaxis()->SetRangeUser(20,34);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "htdc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
//~ sprintf(hname, "htdc_cut_2%d", ich);
//~ hp1dcut = DrTH1F(dir, hname, "");
//~ if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
//~ hp1dcut->SetLineColor(kGreen);
//~ hp1dcut->DrawClone("SAME");
}
// -------------------------- qdc 1 ----------------------------------
ich = 1;
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("MPPC 2;Charge [200 fC/bin]");
hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "hadc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
//~ sprintf(hname, "hadc_cut_2%d", ich);
//~ hp1dcut = DrTH1F(dir, hname, "");
//~ hp1dcut->SetLineColor(kGreen);
//~ hp1dcut->DrawClone("SAME");
}
// -------------------------- TDC 1 ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("MPPC 2;Time [ns]");
if(rebin_tdcs) hp1d->Rebin(rebin_tdcs);
//~ hp1d->GetXaxis()->SetRangeUser(0,22);
hp1d->GetXaxis()->SetRangeUser(20,34);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "htdc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
//~ sprintf(hname, "htdc_cut_2%d", ich);
//~ hp1dcut = DrTH1F(dir, hname, "");
//~ if(rebin_tdcs) hp1dcut->Rebin(rebin_tdcs);
//~ hp1dcut->SetLineColor(kGreen);
//~ hp1dcut->DrawClone("SAME");
}
sprintf(fullname, "gif/a_%s.gif", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'q') != NULL ) {
gStyle->SetOptStat(1111);
gStyle->SetOptFit(0);
c[++cc] = new TCanvas("q", "q", 0, 0, 1200, 700);
c[cc]->Divide(3,2);
int draw_cuts = 1;
int rebin_tdcs = 1;
int ccccdi = 0;
int ich;
/*
ich = 2;
// -------------------------- qdc ref ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", ich);
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->GetXaxis()->SetRangeUser(0,100);
hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
sprintf(hname, "hadc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "hadc_cut_2%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
// -------------------------- TDC ref ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->Rebin(8);
hp1d->GetXaxis()->SetRangeUser(5,105);
hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
sprintf(hname, "htdc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(8);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdc_cut_2%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(8);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
*/
// -------------------------- TDC DIFF ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(0);
sprintf(hname, "htdcdiff");
hp1d = DrTH1F(dir, hname, "");
hp1d->Rebin(rebin_tdcs);
hp1d->GetXaxis()->SetRangeUser(-15,15);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*2);
hp1d->GetYaxis()->SetRangeUser(0.,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "htdcdiff_cut");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdcdiff_cut_2");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
// -------------------------- qdc MPPC ----------------------------------
ich = 0;
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->GetXaxis()->SetRangeUser(0,1500);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "hadc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "hadc_cut_2%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
// -------------------------- TDC MPPC ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->Rebin(rebin_tdcs);
hp1d->GetXaxis()->SetRangeUser(0,20);
hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "htdc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdc_cut_2%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
// -------------------------- qdc MPPC ----------------------------------
ich = 1;
ccccdi++;
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->GetXaxis()->SetRangeUser(0,1500);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "hadc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "hadc_cut_2%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
// -------------------------- TDC MPPC ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->Rebin(rebin_tdcs);
hp1d->GetXaxis()->SetRangeUser(0,20);
hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->DrawClone();
if(draw_cuts) {
sprintf(hname, "htdc_cut%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdc_cut_2%d", ich);
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
sprintf(fullname, "gif/q_%s.gif", fname); c[cc]->SaveAs(fullname);
//~ sprintf(fullname, "%s_j2.eps", fname); c[cc]->SaveAs(fullname);
if(batch_q) gSystem->Exit(1);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'n') != NULL ) {
gStyle->SetOptStat("e");
//~ gStyle->SetOptFit(1);
c[++cc] = new TCanvas("d", "d", 0, 0, 900, 600);
//~ c[cc]->Divide(2,2);
int draw_cuts = 0;
int rebin_tdcs = 1;
double drawmin = -17;
double drawmax = +17;
double fitmin = -10;
double fitmax = 10;
double fitcenter = 0.0;
double HalfMax, MinHalfMax, MaxHalfMax, FWHM, FWHMc;
int ccccdi = 0;
int ich = 0;
TF1 *f_fit;
TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)"); fcg->SetNpx(300);
fcg->SetParName(0,"Const"); fcg->SetParName(1,"G_Const");
fcg->SetParName(2,"G_Mean" ); fcg->SetParName(3,"G_Sigma" );
TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)"); fcgg->SetNpx(300);
fcgg->SetParName(0,"Const");
fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
(c[cc]->cd(++ccccdi))->SetLogy(0);
sprintf(hname, "htdcdiff");
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->Rebin(rebin_tdcs);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*1.1);
hp1d->SetTitle(";Coincidence Time [ns];Counts");
f_fit = fcgg;
f_fit->SetParameters(hp1d->GetMaximum()*0.25,
hp1d->GetMaximum()*0.7, 0, 0.4,
hp1d->GetMaximum()*0.05, 0, 0.05);
f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
f_fit->SetParLimits(2, fitmin,fitmax);
f_fit->SetParLimits(3, 0.2, 1);
f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
f_fit->SetParLimits(5, fitmin,fitmax);
f_fit->SetParLimits(6, 0.02, 0.3);
//~ f_fit = fcgg;
//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
//~ hp1d->GetMaximum()*0.5, 0, 0.05,
//~ hp1d->GetMaximum()*0.4, 0, 0.5);
//~ f_fit = fcg;
//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
//~ hp1d->GetMaximum()*0.9, 0, 0.05);
f_fit->SetRange(fitmin,fitmax);
f_fit->SetLineWidth(1);
hp1d->Fit(f_fit, "Q", "", fitmin,fitmax);
HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
FWHM = MaxHalfMax - MinHalfMax;
HalfMax = (f_fit->GetMaximum(fitmin, fitmax) + f_fit->GetParameter(0))/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
FWHMc = MaxHalfMax - MinHalfMax;
//~ printf("FWHM = %lf ns | FWHMc = %lf ns\n", FWHM, FWHMc);
printf("FWHM (peak only) = %.3lf ns\n", FWHMc);
leg[++legi] = new TLegend(0.122768, 0.740418, 0.420759, 0.878049);
leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.03); leg[legi]->SetTextAlign(22);
for(int i=0; i<6; i++) fgg->SetParameter(i, f_fit->GetParameter(i+1));
//~ fgg->SetRange(drawmin, drawmax); fgg->SetNpx(300); fgg->SetLineColor(kBlue); fgg->SetLineWidth(1); fgg->DrawClone("LSAME");
double Icherenkovs = fgg->Integral(drawmin, drawmax);
double Iall = f_fit->Integral(fitmin, fitmax);
double Rsn = Icherenkovs/(Iall-Icherenkovs);
double Ientr = hp1d->GetEntries()*hp1d->GetBinWidth(1);
double Rsne = Icherenkovs/(Ientr-Icherenkovs);
double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1)*hp1d->GetBinWidth(1);
double Rsnh = Icherenkovs/(Ihist-Icherenkovs);
//~ printf("ICherenkov = %lf | IAll = %lf | Ientr = %lf | Ihist = %lf\n", Icherenkovs, Iall, Ientr, Ihist);
//~ printf("TF1 -> S/N = %lf || TH1F(e) -> S/N = %lf || TH1F(i) -> S/N = %lf\n", Rsn, Rsne, Rsnh);
printf("S/N (vs. All) = %lf\n", Rsnh);
if(draw_cuts) {
sprintf(hname, "htdcdiff_cut");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdcdiff_cut_2");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
//~ #define USE_NOISE_FILE
#ifdef USE_NOISE_FILE
//get ROOT file with histograms
char fnameroot[1024];
TFile * rootfile2;
TDirectory *dir2;
sprintf(fname, "run054");
sprintf(fnameroot, "root/%s.root", fname);
rootfile2 = (TFile *) gROOT->FindObject(fname);
if(rootfile2==NULL) rootfile2 = new TFile(fnameroot);
if(rootfile2==NULL) {
printf("Cannot open root file 2 %s!!!\n",fnameroot);
return;
}
dir2 = (TDirectory*) rootfile2;
sprintf(hname, "htdcdiff");
hp1dcut = DrTH1F(dir2, hname, "");
hp1dcut->SetLineColor(kMagenta);
hp1dcut->DrawClone("SAME");
leg[legi]->AddEntry((TH1F*)hp1d->Clone(), "With source (30 min)", "L");
leg[legi]->AddEntry((TH1F*)hp1dcut->Clone(), "W/o source (30 min)", "L");
sprintf(fullname, "Fit FWHM = %.3lf ns", FWHMc);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "L");
leg[legi]->Draw();
#endif
sprintf(fullname, "gif/n_%s.gif", fname); c[cc]->SaveAs(fullname);
if(batch_q) gSystem->Exit(1);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'f') != NULL ) {
DrSetDrawStyle(0.05);
gStyle->SetOptStat("e");
//~ gStyle->SetOptFit(1);
gStyle->SetPadRightMargin(0.03);
gStyle->SetPadTopMargin(0.05);
gStyle->SetPadBottomMargin(0.12);
gStyle->SetStatFontSize(0.05);
gStyle->SetStatX(0.97);
gStyle->SetStatY(0.95);
c[++cc] = new TCanvas("f", "f", 0, 460, 720, 400);
//~ c[cc]->Divide(2,2);
int draw_cuts = 0;
int rebin_tdcs = 1;
double drawmin = -10;
double drawmax = +10;
//~ double fitmin = inmin;
//~ double fitmax = inmax;
double fitmin = -10;
double fitmax = 10;
double fitcenter = 0.0;
double HalfMax, MinHalfMax, MaxHalfMax, FWHM, FWHMc;
int ccccdi = 0;
int ich = 0;
TF1 *f_fit;
TF1 *f_integrate;
TF1 *fc = new TF1("fc", "pol0(0)");
TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)"); fcg->SetNpx(300);
fcg->SetParName(0,"Const"); fcg->SetParName(1,"G_Const");
fcg->SetParName(2,"G_Mean" ); fcg->SetParName(3,"G_Sigma" );
TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)"); fcgg->SetNpx(300);
fcgg->SetParName(0,"Const");
fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
TF1 *fcggg = new TF1("fcggg", "pol0(0)+gaus(1)+gaus(4)+gaus(7)"); fcggg->SetNpx(300);
fcggg->SetParName(0,"Const");
fcggg->SetParName(1,"G1_Const"); fcggg->SetParName(2,"G1_Mean" ); fcggg->SetParName(3,"G1_Sigma" );
fcggg->SetParName(4,"G2_Const"); fcggg->SetParName(5,"G2_Mean" ); fcggg->SetParName(6,"G2_Sigma" );
fcggg->SetParName(7,"G3_Const"); fcggg->SetParName(8,"G3_Mean" ); fcggg->SetParName(9,"G3_Sigma" );
//~ 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]");
TF1 *f_nocrystal = new TF1("f_nocrystal", "[0]*([1]**2*TMath::Gaus(x,[2],sqrt(2)*[3],1)+(1-[1])**2*TMath::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]");
f_nocrystal->SetParName(0, "const");
f_nocrystal->SetParName(1, "ratio");
f_nocrystal->SetParName(2, "WW_mean");
f_nocrystal->SetParName(3, "Window_sigma");
f_nocrystal->SetParName(4, "MM_mean");
f_nocrystal->SetParName(5, "MCP_sigma");
f_nocrystal->SetParName(6, "WM1_mean");
f_nocrystal->SetParName(7, "WM2_mean");
f_nocrystal->SetParName(8, "Constant");
//~ 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)))");
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)))");
f_crystal->SetParName(0, "dark");
f_crystal->SetParName(1, "const");
f_crystal->SetParName(2, "Prob_Crystal");
f_crystal->SetParName(3, "CC_mean");
f_crystal->SetParName(4, "Crystal_sigma");
f_crystal->SetParName(5, "Prob_window");
f_crystal->SetParName(6, "WW_mean");
f_crystal->SetParName(7, "WW_sigma");
f_crystal->SetParName(8, "MM_mean");
f_crystal->SetParName(9, "MM_sigma");
f_crystal->SetParName(10, "CW1_mean");
f_crystal->SetParName(11, "CW2_mean");
f_crystal->SetParName(12, "CM1_mean");
f_crystal->SetParName(13, "CM2_mean");
f_crystal->SetParName(14, "WM1_mean");
f_crystal->SetParName(15, "WM2_mean");
//~ TF1 *f_crystal1 = new TF1("f_crystal1", "[0]+[1]*([2]**2*TMath::Gaus(x,[3],sqrt(2)*[4],1)+(1.-5./3*[2])**2*TMath::Gaus(x,[5],sqrt(2)*[6],1)+(4.*[2]/6.)**2*TMath::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))");
TF1 *f_crystal1 = new TF1("f_crystal1", "[0]+[1]*([2]**2*TMath::Gaus(x,[3],sqrt(2)*[4],1)+(1.-(1.+[15])*[2])**2*TMath::Gaus(x,[5],sqrt(2)*[6],1)+([15]*[2])**2*TMath::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))");
f_crystal1->SetParName(0, "dark");
f_crystal1->SetParName(1, "const");
f_crystal1->SetParName(2, "Prob_Window");
f_crystal1->SetParName(3, "WW_mean");
f_crystal1->SetParName(4, "Window_sigma");
f_crystal1->SetParName(5, "CC_mean");
f_crystal1->SetParName(6, "Crystal_sigma");
f_crystal1->SetParName(7, "MM_mean");
f_crystal1->SetParName(8, "MCP_sigma");
f_crystal1->SetParName(9, "CW1_mean");
f_crystal1->SetParName(10, "CW2_mean");
f_crystal1->SetParName(11, "WM1_mean");
f_crystal1->SetParName(12, "WM2_mean");
f_crystal1->SetParName(13, "CM1_mean");
f_crystal1->SetParName(14, "CM2_mean");
f_crystal1->SetParName(15, "ratio_mcp/window");
f_crystal1->SetParName(16, "Sup._CW");
double ratio_r=0;
(c[cc]->cd(++ccccdi))->SetLogy(0);
sprintf(hname, "htdcdiff_cut");
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->Rebin(rebin_tdcs);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*1.1);
sprintf(fullname, ";Coincidence Time [ns];Counts");
hp1d->SetTitle(fullname);
if( strcmp(fitopt, "cgg")==0 ) {
f_fit = fcgg;
f_fit->SetParameters(hp1d->GetMaximum()*0.1,
hp1d->GetMaximum()*0.4, 0, 0.8,
hp1d->GetMaximum()*0.5, 0, 0.05);
f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
f_fit->SetParLimits(2, fitmin,fitmax);
f_fit->SetParLimits(3, 0.05, 1);
f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
f_fit->SetParLimits(5, fitmin,fitmax);
f_fit->SetParLimits(6, 0.01, 1);
} else if(strcmp(fitopt, "cggg")==0) {
f_fit = fcggg;
f_fit->SetParameters(hp1d->GetMaximum()*0.1,
hp1d->GetMaximum()*0.4, 0, 0.4,
hp1d->GetMaximum()*0.4, 0, 0.2,
hp1d->GetMaximum()*0.1, 0, 0.05);
f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
f_fit->SetParLimits(2, fitmin,fitmax);
f_fit->SetParLimits(3, 0.01, 1);
f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
f_fit->SetParLimits(5, fitmin,fitmax);
f_fit->SetParLimits(6, 0.05, 1);
f_fit->SetParLimits(7, 0, hp1d->GetMaximum());
f_fit->SetParLimits(8, fitmin,fitmax);
f_fit->SetParLimits(9, 0.2, 1);
} else if(strcmp(fitopt, "nocrystal")==0)
{ f_fit = f_nocrystal;
//~ f_fit->SetParameters(4000,0.87,0.03,0.03,0.28,0.1,200,200,0.01);
f_fit->SetParameters(200,0.8,0.0,0.03,0.0,0.1,-0.13,0.13,10.16);
f_fit->SetParLimits(0,0,10000);
f_fit->SetParLimits(1,0,1);
f_fit->SetParLimits(2,-0.02,0.02);
f_fit->SetParLimits(3,0.01,0.3);
f_fit->SetParLimits(4,-0.02,0.02);
f_fit->SetParLimits(5,0,0.8);
f_fit->SetParLimits(6,-0.5,0.0);
f_fit->SetParLimits(7,0.0,0.5);
f_fit->SetParLimits(8,0.00,9000.5);
//~ f_fit->FixParameter(6, 100);
//~ f_fit->FixParameter(7, 100);
//~
} else if(strcmp(fitopt, "crystal")==0)
{
f_fit = f_crystal;
//~ 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);
f_fit->SetParameter(0, 0.1);
f_fit->SetParameter(1, 300);
f_fit->SetParameter(2, 0.8);
f_fit->SetParameter(3, 0.);
f_fit->SetParameter(4, 0.07);
f_fit->SetParameter(5, 0.01);
f_fit->SetParameter(6, 0.0);
f_fit->SetParameter(7, 0.028);
f_fit->SetParameter(8, 0.0);
f_fit->SetParameter(9, 0.15);
f_fit->SetParameter(10, 0.0);
f_fit->SetParameter(11, 0.0);
f_fit->SetParameter(12, -0.2);
f_fit->SetParameter(13, 0.2);
f_fit->SetParameter(14, -0.2);
f_fit->SetParameter(15, 0.2);
f_fit->SetParameter(15, 0.2);
f_fit->SetParLimits(0, 0, 10);
f_fit->SetParLimits(1, 0, 1000);
f_fit->SetParLimits(2, 0, 1);
f_fit->SetParLimits(3, -0.02, 0.02); f_fit->SetParLimits(4, 0, 0.3);
f_fit->SetParLimits(5, 0, 1); f_fit->SetParLimits(6, -0.02, 0.02);
f_fit->SetParLimits(7, 0, 0.035);
f_fit->SetParLimits(8, -0.02, 0.02);
f_fit->SetParLimits(9, 0, 0.2);
f_fit->SetParLimits(10, -0.02, 0.02);
f_fit->SetParLimits(11, -0.02, 0.02);
f_fit->SetParLimits(12, -0.04, 0);
f_fit->SetParLimits(13, 0, 0.04);
f_fit->SetParLimits(14, -0.04, 0);
f_fit->SetParLimits(15, 0, 0.04);
} else if(strcmp(fitopt, "crystal1")==0)
{
f_fit = f_crystal1;
//~ 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);
f_fit->SetParameter(0, 100.1);
f_fit->SetParameter(1, 1);
f_fit->SetParameter(2, 0.2); //prob. window
f_fit->SetParameter(3, 0.0);
f_fit->SetParameter(4, 0.02); // window sigma
f_fit->SetParameter(5, 0.0);
f_fit->SetParameter(6, 120/3.3/1000); //crystal sigma
f_fit->SetParameter(7, 0.0);
f_fit->SetParameter(8, 0.15); // mcp sigma
f_fit->SetParameter(9, -0.1825); //cw
f_fit->SetParameter(10, 0.1825);
f_fit->SetParameter(11, -0.2); //wm
f_fit->SetParameter(12, 0.2);
f_fit->SetParameter(13, -0.3); //cm
f_fit->SetParameter(14, 0.3);
f_fit->SetParameter(15, 1/0.6795 -1.);
f_fit->SetParameter(16, 0.1);
f_fit->SetParameter(1, hp1d->GetMaximum()/(f_fit->Eval(0)-f_fit->GetParameter(0)));
f_fit->SetParLimits(0, 0, 400);
f_fit->SetParLimits(1, 0, 100000);
f_fit->SetParLimits(2, 0, 1.);
f_fit->SetParLimits(3, -0.025, 0.025);
f_fit->SetParLimits(4, 0.02, 0.08);
f_fit->SetParLimits(5, -0.025, 0.025);
f_fit->SetParLimits(6, 0.01, 0.2);
f_fit->SetParLimits(7, -0.025, 0.025);
f_fit->SetParLimits(8, 0, 0.5);
f_fit->SetParLimits(9, -0.6, 0.02);
f_fit->SetParLimits(10, -0.02, 0.6);
f_fit->SetParLimits(11, -0.3, -0.1);
f_fit->SetParLimits(12, 0.1, 0.3);
f_fit->SetParLimits(13, -0.7, -0.1);
f_fit->SetParLimits(14, 0.1, 0.7);
f_fit->SetParLimits(15, 0, 0.999);
f_fit->SetParLimits(16, 0, 1);
f_fit->FixParameter(4,0.02668);
//~ f_fit->FixParameter(6,0.048);
f_fit->FixParameter(8,0.0877);
f_fit->FixParameter(11,-0.19);
f_fit->FixParameter(12,0.19);
//~
ratio_r = 1/0.6795 -1.;
//~ ratio_r = 1/0.8089 -1.;
f_fit->FixParameter(15,ratio_r);
//~ f_fit->FixParameter(16,1);
//~ f_fit->FixParameter(2,0.52);
//~ f_fit->FixParameter(6,170/3.3/1000);
//~ f_fit->FixParameter(3,0);
//~ f_fit->FixParameter(5,0);
//~ f_fit->FixParameter(7,0);
//~ cout<<"hogeeeeeeeeeeee "<<f_fit->GetParameter(1)*((1.-5./3.*(f_fit->GetParameter(2)))*(4./6.*(f_fit->GetParameter(2))))<<endl;
}
//~ f_fit->SetRange(fitmin,fitmax);
f_fit->SetRange(-10,10);
f_fit->SetLineWidth(1); f_fit->SetLineStyle(2);
f_fit->SetNpx(300);
//~ hp1d->Fit(f_fit, "QWW", "", fitmin,fitmax);
//~ hp1d->Fit(f_fit, "ILM+", "", -0.4,0.4);
//~
hp1d->Fit(f_fit, "LM+", "", -.5,.5);
//~ hp1d->Draw();
//~ f_fit->Draw("same"); //f_fit->SetLineColor(6);
f_fit->SetLineColor(1);
hp1d->GetXaxis()->SetRangeUser(-1.5,2);
gStyle->SetOptFit(1111);
double wmsigma = sqrt(f_fit->GetParameter(3)**2 + f_fit->GetParameter(5)**2);
double cm1sigma = sqrt(f_fit->GetParameter(6)**2 + f_fit->GetParameter(8)**2);
double wm1sigma = sqrt(f_fit->GetParameter(4)**2 + f_fit->GetParameter(8)**2);
double cw1sigma = sqrt(f_fit->GetParameter(6)**2 + f_fit->GetParameter(4)**2);
//kobayashi
TPaveStats* st1 = (TPaveStats*) hp1d->FindObject("stats");
st1->SetX1NDC(0.63); st1->SetX2NDC(0.99);
st1->SetY1NDC(0.24); st1->SetY2NDC(0.99);
double pw = f_fit->GetParameter(2);
double pm = f_fit->GetParameter(15)*f_fit->GetParameter(2);
double pc = 1-pw - pm;
double rate = pw**2 + pm**2 + 2*pw*pm;
double sumrate = pw**2+pm**2+pc**2+2*pw*pm+2*pm*pc+2*pc*pw*f_fit->GetParameter(16)**2;
cout<<"pw = "<<pw<<" pm = "<<pm<<" pc = "<<pc<<endl;
cout<<"rate = "<<rate<<" sumrate = "<<sumrate<<" rate/sumrate = "<<rate/sumrate<<endl;
cout<<" FWHM_sim = "<<f_fit->GetParameter(6)*sqrt(2)*2.35<<endl;
if(strcmp(fitopt, "nocrystal")==0)
{
//~ TF1 *fg_WW = new TF1("fg_WW","[0]*exp(-(x-[1])**2/2/([3]**2))/sqrt(2*3.141592*[2]*[2])");
//~ TF1 *fg_MM = new TF1("fg_MM","[0]*exp(-(x-[1])**2/2/([3]**2))/sqrt(2*3.141592*[2]*[2])");
//~ TF1 *fg_WM1 = new TF1("fg_WM1","[0]*exp(-(x-[1])**2/2/([3]**2))/sqrt(2*3.141592*[2]*[2])");
//~ TF1 *fg_WM2 = new TF1("fg_WM2","[0]*exp(-(x-[1])**2/2/([3]**2))/sqrt(2*3.141592*[2]*[2])");
TF1 *fg_WW = new TF1("fg_WW","[0]*TMath::Gaus(x,[1],[2],1)");fg_WW->SetNpx(300);
TF1 *fg_MM = new TF1("fg_MM","[0]*TMath::Gaus(x,[1],[2],1)");fg_MM->SetNpx(300);
TF1 *fg_WM1 = new TF1("fg_WM1","[0]*TMath::Gaus(x,[1],[2],1)"); fg_WM1->SetNpx(300);
TF1 *fg_WM2 = new TF1("fg_WM2","[0]*TMath::Gaus(x,[1],[2],1)"); fg_WM2->SetNpx(300);
fg_WW -> SetParameters(f_fit->GetParameter(1)**2*f_fit->GetParameter(0),f_fit->GetParameter(2),sqrt(2)*f_fit->GetParameter(3));
fg_MM -> SetParameters((1-f_fit->GetParameter(1))**2*f_fit->GetParameter(0),f_fit->GetParameter(4),sqrt(2)*f_fit->GetParameter(5));
fg_WM1 -> SetParameters(f_fit->GetParameter(0)*(1-f_fit->GetParameter(1))*f_fit->GetParameter(1),f_fit->GetParameter(6),wmsigma);
fg_WM2 -> SetParameters(f_fit->GetParameter(0)*(1-f_fit->GetParameter(1))*f_fit->GetParameter(1),f_fit->GetParameter(7),wmsigma);
fg_WW->Draw("sames"); fg_WW->SetLineColor(kYellow+1); fg_WW->SetRange(-10,10);
fg_MM->Draw("sames"); fg_MM->SetLineColor(kMagenta-7); fg_MM->SetRange(-10,10);
fg_WM1->Draw("sames"); fg_WM1->SetLineColor(kRed-4); fg_WM1->SetRange(-10,10);
fg_WM2->Draw("sames"); fg_WM2->SetLineColor(kRed-4); fg_WM2->SetRange(-10,10);
double integ_WW, integ_MM, integ_WM1, integ_WM2;
integ_WW = fg_WW->Integral(-4,4);
}
if(strcmp(fitopt, "crystal1")==0)
{
TF1 *fg_CC = new TF1("fg_CC","[0]*TMath::Gaus(x,[1],[2],1)"); fg_CC->SetNpx(300);
TF1 *fg_WW = new TF1("fg_WW","[0]*TMath::Gaus(x,[1],[2],1)"); fg_WW->SetNpx(300);
TF1 *fg_MM = new TF1("fg_MM","[0]*TMath::Gaus(x,[1],[2],1)"); fg_MM->SetNpx(300);
TF1 *fg_WM1 = new TF1("fg_WM1","[0]*TMath::Gaus(x,[1],[2],1)"); fg_WM1->SetNpx(300);
TF1 *fg_WM2 = new TF1("fg_WM2","[0]*TMath::Gaus(x,[1],[2],1)"); fg_WM2->SetNpx(300);
TF1 *fg_CW1 = new TF1("fg_CW1","[0]*TMath::Gaus(x,[1],[2],1)"); fg_CW1->SetNpx(300);
TF1 *fg_CW2 = new TF1("fg_CW2","[0]*TMath::Gaus(x,[1],[2],1)"); fg_CW2->SetNpx(300);
TF1 *fg_CM1 = new TF1("fg_CM1","[0]*TMath::Gaus(x,[1],[2],1)"); fg_CM1->SetNpx(300);
TF1 *fg_CM2 = new TF1("fg_CM2","[0]*TMath::Gaus(x,[1],[2],1)"); fg_CM2->SetNpx(300);
// normalize
fg_CC -> SetParameters((1.-(1+f_fit->GetParameter(15))*f_fit->GetParameter(2))**2*f_fit->GetParameter(1),f_fit->GetParameter(5),sqrt(2)*f_fit->GetParameter(6));
fg_WW -> SetParameters(f_fit->GetParameter(2)**2*f_fit->GetParameter(1),f_fit->GetParameter(3),sqrt(2)*f_fit->GetParameter(4));
fg_MM -> SetParameters((f_fit->GetParameter(15)*f_fit->GetParameter(2))**2*f_fit->GetParameter(1),f_fit->GetParameter(7),sqrt(2)*f_fit->GetParameter(8));
fg_WM1 -> SetParameters(f_fit->GetParameter(1)*(f_fit->GetParameter(15)*f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(11),wm1sigma);
fg_WM2 -> SetParameters(f_fit->GetParameter(1)*(f_fit->GetParameter(15)*f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(12),wm1sigma);
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);
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);
fg_CM1 -> SetParameters(f_fit->GetParameter(1)*(ratio_r*f_fit->GetParameter(2))*(1.-(1+f_fit->GetParameter(15))*f_fit->GetParameter(2)),f_fit->GetParameter(13),cm1sigma);
fg_CM2 -> SetParameters(f_fit->GetParameter(1)*(ratio_r*f_fit->GetParameter(2))*(1.-(1+f_fit->GetParameter(15))*f_fit->GetParameter(2)),f_fit->GetParameter(14),cm1sigma);
//non normalize
//~ fg_CC -> SetParameters((1.-5./3.*f_fit->GetParameter(2))**2*f_fit->GetParameter(1),f_fit->GetParameter(5),sqrt(2)*f_fit->GetParameter(6));
//~ fg_WW -> SetParameters(f_fit->GetParameter(2)**2*f_fit->GetParameter(1),f_fit->GetParameter(3),sqrt(2)*f_fit->GetParameter(4));
//~ fg_MM -> SetParameters((4./6.*f_fit->GetParameter(2))**2*f_fit->GetParameter(1),f_fit->GetParameter(7),sqrt(2)*f_fit->GetParameter(8));
//~ fg_WM1 -> SetParameters(f_fit->GetParameter(1)*(4./6.*f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(11),wm1sigma);
//~ fg_WM2 -> SetParameters(f_fit->GetParameter(1)*(4./6.*f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(12),wm1sigma);
//~ fg_CW1 -> SetParameters(f_fit->GetParameter(1)*(1.-5./3*f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(9),cw1sigma);
//~ fg_CW2 -> SetParameters(f_fit->GetParameter(1)*(1.-5./3*f_fit->GetParameter(2))*f_fit->GetParameter(2),f_fit->GetParameter(10),cw1sigma);
//~ fg_CM1 -> SetParameters(f_fit->GetParameter(1)*(4./6.*f_fit->GetParameter(2))*(1.-5./3*f_fit->GetParameter(2)),f_fit->GetParameter(13),cm1sigma);
//~ fg_CM2 -> SetParameters(f_fit->GetParameter(1)*(4./6.*f_fit->GetParameter(2))*(1.-5./3*f_fit->GetParameter(2)),f_fit->GetParameter(14),cm1sigma);
fg_CC->Draw("sames"); fg_CC->SetLineColor(kCyan-5); fg_CC->SetRange(-10,10);
fg_WW->Draw("sames"); fg_WW->SetLineColor(kYellow+1); fg_WW->SetRange(-10,10);
fg_MM->Draw("sames"); fg_MM->SetLineColor(kMagenta-7); fg_MM->SetRange(-10,10);
fg_CW1->Draw("sames"); fg_CW1->SetLineColor(kGreen+2); fg_CW1->SetRange(-10,10);
fg_CW2->Draw("sames"); fg_CW2->SetLineColor(kGreen+2); fg_CW2->SetRange(-10,10);
fg_WM1->Draw("sames"); fg_WM1->SetLineColor(kRed-4); fg_WM1->SetRange(-10,10);
fg_WM2->Draw("sames"); fg_WM2->SetLineColor(kRed-4); fg_WM2->SetRange(-10,10);
fg_CM1->Draw("sames"); fg_CM1->SetLineColor(kBlue+2); fg_CM1->SetRange(-10,10);
fg_CM2->Draw("sames"); fg_CM2->SetLineColor(kBlue+2); fg_CM2->SetRange(-10,10);
}
HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
FWHM = MaxHalfMax - MinHalfMax;
HalfMax = (f_fit->GetMaximum(fitmin, fitmax) + f_fit->GetParameter(0))/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, f_fit->GetMaximumX(fitmin, fitmax));
MaxHalfMax = f_fit->GetX(HalfMax, f_fit->GetMaximumX(fitmin, fitmax), fitmax);
FWHMc = MaxHalfMax - MinHalfMax;
printf("(f) HalfMax = %.3lf ns | MinHalfMax = %.3lf ns | MaxHalfMax = %.3lf ns\n", HalfMax, MinHalfMax, MaxHalfMax);
printf("(f) FWHM (peak only) = %.3lf ns\n", FWHMc);
f_integrate = fgg;
for(int i=0; i<6; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
//~ for(int i=0; i<9; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
//~ fgg->SetRange(drawmin, drawmax); fgg->SetNpx(300); fgg->SetLineColor(kBlue); fgg->SetLineWidth(1); fgg->DrawClone("LSAME");
double Icherenkovs = f_integrate->Integral(drawmin, drawmax)/hp1d->GetBinWidth(1);
double Iall = f_fit->Integral(fitmin, fitmax)/hp1d->GetBinWidth(1);
double Rsn = Icherenkovs/(Iall-Icherenkovs);
//~ double Ientr = hp1d->GetEntries()*hp1d->GetBinWidth(1);
double Ientr = hp1d->GetEntries();
double Rsne = Icherenkovs/(Ientr-Icherenkovs);
//~ double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1)*hp1d->GetBinWidth(1);
double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1);
double Rsnh = Icherenkovs/(Ihist-Icherenkovs);
double Icherenkovs_10ns = f_integrate->Integral(-10, 10)/hp1d->GetBinWidth(1);
fc->SetParameter(0, f_fit->GetParameter(0));
double Inoise_10ns = fc->Integral(-10, 10)/hp1d->GetBinWidth(1);
double R10ns = Icherenkovs_10ns/Inoise_10ns;
double Icherenkovs_4ns = f_integrate->Integral(-4, 4)/hp1d->GetBinWidth(1);
fc->SetParameter(0, f_fit->GetParameter(0));
double Inoise_4ns = fc->Integral(-4, 4)/hp1d->GetBinWidth(1);
double R4ns = Icherenkovs_4ns/Inoise_4ns;
double Icherenkovs_2ns = f_integrate->Integral(-2, 2)/hp1d->GetBinWidth(1);
fc->SetParameter(0, f_fit->GetParameter(0));
double Inoise_2ns = fc->Integral(-2, 2)/hp1d->GetBinWidth(1);
double R2ns = Icherenkovs_2ns/Inoise_2ns;
//~ printf("ICherenkov = %lf | IAll = %lf | Ientr = %lf | Ihist = %lf\n", Icherenkovs, Iall, Ientr, Ihist);
printf("Icherenkovs_4ns = %.0lf\n", Icherenkovs_4ns);
//~ printf("TF1 -> S/N = %lf || TH1F(e) -> S/N = %lf || TH1F(i) -> S/N = %lf\n", Rsn, Rsne, Rsnh);
//~ printf("S/N (All) = %lf\n", Rsnh);
printf("S/N (10ns) = %.2lf\n", R10ns);
printf("S/N (4 ns) = %.2lf\n", R4ns);
printf("S/N (2 ns) = %.2lf\n", R2ns);
leg[++legi] = new TLegend(0.12,0.4,0.43,0.925);
leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(12);
sprintf(fullname, "FWHM_fit = %.0lf ps", 1000*FWHMc);
//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "L");
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "FWHM_c = %.0lf ps", 1000*f_fit->GetParameter(6)*2.35*sqrt(2));
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
//~ sprintf(fullname, "S/N (#pm10ns) = %.1lf", R10ns);
//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
//~ sprintf(fullname, "S/N (#pm 4ns) = %.1lf", R4ns);
//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
//~ sprintf(fullname, "S/N (#pm 2ns) = %.1lf", R2ns);
//~ leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
//kobayashi
sprintf(fullname, "R_w = %.2f" , pw);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "R_m = %.2f" , pm);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "R_c = %.2f" , pc);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "R_w&m = %.2f" , rate);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "R_sum = %.2f" , sumrate);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "R_w&m/sum = %.2f" , rate/sumrate);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
leg[legi]->Draw();
if(draw_cuts) {
sprintf(hname, "htdcdiff_cut");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdcdiff_cut_2");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
sprintf(fullname, "gif/f_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/f_%s.eps", fname); c[cc]->SaveAs(fullname);
cout<< "maximum x value = " << hp1d->GetXaxis()->GetBinCenter( hp1d->GetMaximumBin())<<endl;
if(batch_q) gSystem->Exit(1);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'd') != NULL ) {
DrSetDrawStyle(0.05);
gStyle->SetOptStat("e");
//~ gStyle->SetOptStat(0);
gStyle->SetOptFit(1);
//~ gStyle->SetOptFit(0);
gStyle->SetPadRightMargin(0.03);
gStyle->SetPadTopMargin(0.05);
gStyle->SetPadBottomMargin(0.12);
gStyle->SetStatFontSize(0.05);
gStyle->SetStatX(0.97);
gStyle->SetStatY(0.95);
c[++cc] = new TCanvas("d", "d", 0, 0, 720, 400);
//~ c[cc]->Divide(2,2);
int draw_cuts = 0;
int rebin_tdcs = 1;
double drawmin = -10;
double drawmax = +10;
//~ double fitmin = inmin;
//~ double fitmax = inmax;
double fitmin = -10;
double fitmax = 10;
double fitcenter = 0.0;
double HalfMax, MinHalfMax, MaxHalfMax, FWHM, FWHMc;
int ccccdi = 0;
int ich = 0;
TF1 *f_fit;
TF1 *f_integrate;
TF1 *fc = new TF1("fc", "pol0(0)");
TF1 *fcg = new TF1("fcg", "pol0(0)+gaus(1)"); fcg->SetNpx(300);
fcg->SetParName(0,"Const"); fcg->SetParName(1,"G_Const");
fcg->SetParName(2,"G_Mean" ); fcg->SetParName(3,"G_Sigma" );
TF1 *fcgg = new TF1("fcgg", "pol0(0)+gaus(1)+gaus(4)"); fcgg->SetNpx(300);
fcgg->SetParName(0,"Const");
fcgg->SetParName(1,"G1_Const"); fcgg->SetParName(2,"G1_Mean" ); fcgg->SetParName(3,"G1_Sigma" );
fcgg->SetParName(4,"G2_Const"); fcgg->SetParName(5,"G2_Mean" ); fcgg->SetParName(6,"G2_Sigma" );
TF1 *fcggg = new TF1("fcggg", "pol0(0)+gaus(1)+gaus(4)+gaus(7)"); fcggg->SetNpx(300);
fcggg->SetParName(0,"Const");
fcggg->SetParName(1,"G1_Const"); fcggg->SetParName(2,"G1_Mean" ); fcggg->SetParName(3,"G1_Sigma" );
fcggg->SetParName(4,"G2_Const"); fcggg->SetParName(5,"G2_Mean" ); fcggg->SetParName(6,"G2_Sigma" );
fcggg->SetParName(7,"G3_Const"); fcggg->SetParName(8,"G3_Mean" ); fcggg->SetParName(9,"G3_Sigma" );
(c[cc]->cd(++ccccdi))->SetLogy(0);
//~ sprintf(hname, "htdcdiff");
sprintf(hname, "htdcdiff_cut");
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->Rebin(rebin_tdcs);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*1.1);
//~ hp1d->GetYaxis()->SetRangeUser(0.0,hp1d->GetMaximum()*2);
sprintf(fullname, ";Coincidence Time [ns];Counts");
hp1d->SetTitle(fullname);
//kobayashi
//~ fitmin = hp1d->GetXaxis()->GetBinCenter(hp1d->GetMaximumBin()) - hp1d->GetRMS()/15;
//~ fitmax = hp1d->GetXaxis()->GetBinCenter(hp1d->GetMaximumBin()) + hp1d->GetRMS()/15;
//~ if( strcmp(fname, "run000")==0 ) {
//~ f_fit = fcg;
//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
//~ hp1d->GetMaximum()*0.9, 0, 0.05);
//~ } else if( strcmp(fname, "run061")==0 ||
//~ strcmp(fname, "run062")==0 ||
//~ strcmp(fname, "run063")==0 ||
//~ strcmp(fname, "run064")==0 ||
//~ strcmp(fname, "run065")==0) {
//~ f_fit = fcgg;
//~ f_fit->SetParameters(hp1d->GetMaximum()*0.1,
//~ hp1d->GetMaximum()*0.5, 0, 0.05,
//~ hp1d->GetMaximum()*0.4, 0, 0.5);
//~ } else {
//~ f_fit = fcgg;
//~ f_fit->SetParameters(hp1d->GetMaximum()*0.25,
//~ hp1d->GetMaximum()*0.7, 0, 0.4,
//~ hp1d->GetMaximum()*0.05, 0, 0.05);
//~ f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
//~ f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
//~ f_fit->SetParLimits(2, fitmin,fitmax);
//~ f_fit->SetParLimits(3, 0.2, 1);
//~ f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
//~ f_fit->SetParLimits(5, fitmin,fitmax);
//~ if( strcmp(fname, "run059")==0 ||
//~ strcmp(fname, "run060")==0)
//~ f_fit->SetParLimits(6, 0.02, 0.3);
//~ else
//~ f_fit->SetParLimits(6, 0.02, 0.2);
//~
//~ }
if( strcmp(fitopt, "cgg")==0 ) {
f_fit = fcgg;
f_fit->SetParameters(hp1d->GetMaximum()*0.1,
hp1d->GetMaximum()*0.4, 0, 0.9,
hp1d->GetMaximum()*0.5, 0, 0.05);
f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
f_fit->SetParLimits(1, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
f_fit->SetParLimits(2, fitmin,fitmax);
f_fit->SetParLimits(3, 0.05, 1);
f_fit->SetParLimits(4, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
f_fit->SetParLimits(5, fitmin,fitmax);
f_fit->SetParLimits(6, 0.01, 1);
} else if(strcmp(fitopt, "cggg")==0) {
f_fit = fcggg;
//original
//~ f_fit->SetParameters(hp1d->GetMaximum()*0.01,
//~ hp1d->GetMaximum()*0.4, 0, 0.4,
//~ hp1d->GetMaximum()*0.4, 0, 0.2,
//~ hp1d->GetMaximum()*0.1, 0, 0.05);
//kobayashi
double init_mean = hp1d->GetXaxis()->GetBinCenter(hp1d->GetMaximumBin());
f_fit->SetParameters(hp1d->GetMaximum()*0.0001,
hp1d->GetMaximum()*0.01, init_mean , 0.04,
hp1d->GetMaximum()*0.01, init_mean + 0.2, 0.03,
hp1d->GetMaximum()*0.01, init_mean - 0.2, 0.03);
f_fit->SetParLimits(0, 0, hp1d->GetMaximum());
f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
f_fit->SetParLimits(2, fitmin,fitmax);
f_fit->SetParLimits(3, 0.001, 1);
f_fit->SetParLimits(4, hp1d->GetMaximum()*0.1, hp1d->GetMaximum());
f_fit->SetParLimits(5, fitmin,fitmax);
f_fit->SetParLimits(6, 0.005, 3);
f_fit->SetParLimits(7, hp1d->GetMaximum()*0.1, hp1d->GetMaximum());
f_fit->SetParLimits(8, fitmin,fitmax);
f_fit->SetParLimits(9, 0.02, 5);
} else if(strcmp(fitopt, "gg")==0) {
f_fit = fgg;
f_fit->SetParameters(hp1d->GetMaximum()*0.4, 0, 0.9,
hp1d->GetMaximum()*0.5, 0, 0.05);
f_fit->SetParLimits(0, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
f_fit->SetParLimits(1, fitmin,fitmax);
f_fit->SetParLimits(2, 0.05, 1);
f_fit->SetParLimits(3, hp1d->GetMaximum()*0.01, hp1d->GetMaximum());
f_fit->SetParLimits(4, fitmin,fitmax);
f_fit->SetParLimits(5, 0.01, 1);
}
f_fit->SetRange(fitmin,fitmax);
f_fit->SetLineWidth(1); f_fit->SetLineStyle(9);
hp1d->Fit(f_fit, "QWW", "", fitmin,fitmax);
//~ hp1d->Fit(f_fit, "QWW", "", -0.5,+0.5);
//~ fg->SetRange(fitmin,fitmax); fg->SetLineWidth(1.0);
//~ fg->SetParameters(f_fit->GetParameter(0+1), f_fit->GetParameter(1+1), f_fit->GetParameter(2+1)); fg->SetLineColor(kMagenta); fg->DrawCopy("LSAME");
//~ fg->SetParameters(f_fit->GetParameter(0+4), f_fit->GetParameter(1+4), f_fit->GetParameter(2+4)); fg->SetLineColor(kGreen); fg->DrawCopy("LSAME");
//~ fg->SetParameters(f_fit->GetParameter(0+7), f_fit->GetParameter(1+7), f_fit->GetParameter(2+7)); fg->SetLineColor(kYellow); fg->DrawCopy("LSAME");
HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
FWHM = MaxHalfMax - MinHalfMax;
HalfMax = (f_fit->GetMaximum(fitmin, fitmax) + f_fit->GetParameter(0))/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, f_fit->GetMaximumX(fitmin, fitmax));
MaxHalfMax = f_fit->GetX(HalfMax, f_fit->GetMaximumX(fitmin, fitmax), fitmax);
FWHMc = MaxHalfMax - MinHalfMax;
//~ printf("FWHM = %lf ns | FWHMc = %lf ns\n", FWHM, FWHMc);
printf("(d) FWHM (peak only) = %.3lf ns\n", FWHMc);
f_integrate = fgg;
for(int i=0; i<6; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
//~ for(int i=0; i<9; i++) f_integrate->SetParameter(i, f_fit->GetParameter(i+1));
//~ fgg->SetRange(drawmin, drawmax); fgg->SetNpx(300); fgg->SetLineColor(kBlue); fgg->SetLineWidth(1); fgg->DrawClone("LSAME");
double Icherenkovs = f_integrate->Integral(drawmin, drawmax)/hp1d->GetBinWidth(1);
double Iall = f_fit->Integral(fitmin, fitmax)/hp1d->GetBinWidth(1);
double Rsn = Icherenkovs/(Iall-Icherenkovs);
//~ double Ientr = hp1d->GetEntries()*hp1d->GetBinWidth(1);
double Ientr = hp1d->GetEntries();
double Rsne = Icherenkovs/(Ientr-Icherenkovs);
//~ double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1)*hp1d->GetBinWidth(1);
double Ihist = hp1d->Integral(1, hp1d->GetNbinsX()-1);
double Rsnh = Icherenkovs/(Ihist-Icherenkovs);
double Icherenkovs_10ns = f_integrate->Integral(-10, 10)/hp1d->GetBinWidth(1);
fc->SetParameter(0, f_fit->GetParameter(0));
double Inoise_10ns = fc->Integral(-10, 10)/hp1d->GetBinWidth(1);
double R10ns = Icherenkovs_10ns/Inoise_10ns;
double Icherenkovs_4ns = f_integrate->Integral(-4, 4)/hp1d->GetBinWidth(1);
fc->SetParameter(0, f_fit->GetParameter(0));
double Inoise_4ns = fc->Integral(-4, 4)/hp1d->GetBinWidth(1);
double R4ns = Icherenkovs_4ns/Inoise_4ns;
double Icherenkovs_2ns = f_integrate->Integral(-2, 2)/hp1d->GetBinWidth(1);
fc->SetParameter(0, f_fit->GetParameter(0));
double Inoise_2ns = fc->Integral(-2, 2)/hp1d->GetBinWidth(1);
double R2ns = Icherenkovs_2ns/Inoise_2ns;
//~ printf("ICherenkov = %lf | IAll = %lf | Ientr = %lf | Ihist = %lf\n", Icherenkovs, Iall, Ientr, Ihist);
printf("Icherenkovs_4ns = %.0lf\n", Icherenkovs_4ns);
//~ printf("TF1 -> S/N = %lf || TH1F(e) -> S/N = %lf || TH1F(i) -> S/N = %lf\n", Rsn, Rsne, Rsnh);
//~ printf("S/N (All) = %lf\n", Rsnh);
printf("S/N (10ns) = %.2lf\n", R10ns);
printf("S/N (4 ns) = %.2lf\n", R4ns);
printf("S/N (2 ns) = %.2lf\n", R2ns);
//~ #define USE_NOISE_FILE
#ifdef USE_NOISE_FILE
//get ROOT file with histograms
char fnameroot[1024];
TFile * rootfile2;
TDirectory *dir2;
//~ sprintf(fname, "run054");
sprintf(fname, "run_309");
sprintf(fnameroot, "root/%s.root", fname);
rootfile2 = (TFile *) gROOT->FindObject(fname);
if(rootfile2==NULL) rootfile2 = new TFile(fnameroot);
if(rootfile2==NULL) {
printf("Cannot open root file 2 %s!!!\n",fnameroot);
return;
}
dir2 = (TDirectory*) rootfile2;
sprintf(hname, "htdcdiff");
hp1dcut = DrTH1F(dir2, hname, "");
hp1dcut->SetLineColor(40);
hp1dcut->SetLineStyle(2);
hp1dcut->DrawClone("SAME");
leg[++legi] = new TLegend(0.125,0.5,0.425,0.925);
leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(12);
leg[legi]->AddEntry((TH1F*)hp1d->Clone(), "With source", "L");
leg[legi]->AddEntry((TH1F*)hp1dcut->Clone(), "W/o source", "L");
#else
leg[++legi] = new TLegend(0.125,0.625,0.425,0.925);
leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(12);
#endif
sprintf(fullname, "FWHM = %.0lf ps", 1000*FWHMc);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "L");
sprintf(fullname, "S/N (#pm10ns) = %.1lf", R10ns);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "S/N (#pm 4ns) = %.1lf", R4ns);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
sprintf(fullname, "S/N (#pm 2ns) = %.1lf", R2ns);
leg[legi]->AddEntry((TF1*)f_fit->Clone(), fullname, "");
leg[legi]->Draw();
if(draw_cuts) {
sprintf(hname, "htdcdiff_cut");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kRed);
hp1dcut->DrawClone("SAME");
sprintf(hname, "htdcdiff_cut_2");
hp1dcut = DrTH1F(dir, hname, "");
hp1dcut->Rebin(rebin_tdcs);
hp1dcut->SetLineColor(kGreen);
hp1dcut->DrawClone("SAME");
}
sprintf(fullname, "gif/d_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/d_%s.eps", fname); c[cc]->SaveAs(fullname);
if(batch_q) gSystem->Exit(1);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'c') != NULL ) {
//~ gStyle->SetOptStat(1111);
//~ gStyle->SetOptFit(0);
c[++cc] = new TCanvas("c", "c", 0, 0, 900, 700);
//~ c[cc]->Divide(2,2);
int ccccdi = 0;
(c[cc]->cd(++ccccdi))->SetLogz(1);
sprintf(hname, "htdccor");
hp2d = DrTH2F(dir, hname, "");
//~ hp2d->GetYaxis()->SetRangeUser(0,1500);
//~ sprintf(fullname, "%s;QDC Reference 1275 keV; QDC Reference Coincidence", fname);
//~ hp2d->SetTitle(fullname);
hp2d->DrawClone("COLZ");
//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
//~ sprintf(hname, "hcor%d",0);
//~ hp2d = DrTH2F(dir, hname, "");
//~ hp2d->DrawClone("COLZ");
//~
//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
//~ sprintf(hname, "hcor%d",1);
//~ hp2d = DrTH2F(dir, hname, "");
//~ hp2d->DrawClone("COLZ");
//~
//~
//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
//~ sprintf(hname, "hdiffcor_%d_%d",0,1);
//~ hp2d = DrTH2F(dir, hname, "");
//~ hp2d->DrawClone("COLZ");
//~
//~ (c[cc]->cd(++ccccdi))->SetLogz(1);
//~ sprintf(hname, "hdiffcor_%d_%d",1,0);
//~ hp2d = DrTH2F(dir, hname, "");
//~ hp2d->DrawClone("COLZ");
sprintf(fullname, "gif/c_%s.gif", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------
if( strchr(plopt, 't') != NULL ) {
gStyle->SetOptStat(1111);
//~ gStyle->SetOptFit(0);
c[++cc] = new TCanvas("c", "c", 0, 0, 700, 400);
//~ c[cc]->Divide(2,2);
double drawmin = -1;
double drawmax = 2;
double fitmin = -1;
double fitmax = 0.75;
double fitcenter = 0.0;
double HalfMax, MinHalfMax, MaxHalfMax, FWHM;
int ccccdi = 0;
//~ (c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hctdc%d", 0);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("Vov = 1.5V; Time [ns];");
hp1d->GetXaxis()->SetRangeUser(drawmin,drawmax);
f_fit = fgg;
f_fit->SetParameters(hp1d->GetMaximum()*0.9, 0, 0.1,
hp1d->GetMaximum()*0.1, 0, 0.5);
//~ f_fit->SetParLimits(1, 0, hp1d->GetMaximum());
//~ f_fit->SetParLimits(2, fitmin,fitmax);
//~ f_fit->SetParLimits(3, 0.2, 1);
//~ f_fit->SetParLimits(4, 0, hp1d->GetMaximum());
//~ f_fit->SetParLimits(5, fitmin,fitmax);
f_fit->SetRange(fitmin,fitmax);
f_fit->SetLineWidth(1);
hp1d->Fit(f_fit, "Q", "", fitmin,fitmax);
HalfMax = f_fit->GetMaximum(fitmin, fitmax)/2.0;
MinHalfMax = f_fit->GetX(HalfMax, fitmin, fitcenter);
MaxHalfMax = f_fit->GetX(HalfMax, fitcenter, fitmax);
FWHM = MaxHalfMax - MinHalfMax;
printf("FWHM = %.3lf ns\n", FWHM);
fg->SetRange(-2,5); fg->SetLineWidth(1.0);
fg->SetParameters(f_fit->GetParameter(0), f_fit->GetParameter(1), f_fit->GetParameter(2));
fg->SetLineColor(kGreen); fg->DrawClone("LSAME");
fg->SetParameters(f_fit->GetParameter(0+3), f_fit->GetParameter(1+3), f_fit->GetParameter(2+3));
fg->SetLineColor(kMagenta); fg->DrawClone("LSAME");
sprintf(fullname, "gif/t_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/t_%s.eps", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------
if( strchr(plopt, 'e') != NULL ) {
DrSetDrawStyle(0.05);
gStyle->SetOptStat("e");
gStyle->SetOptFit(0);
gStyle->SetPadRightMargin(0.03);
gStyle->SetPadTopMargin(0.05);
gStyle->SetPadBottomMargin(0.12);
gStyle->SetStatFontSize(0.05);
gStyle->SetStatX(0.97);
gStyle->SetStatY(0.95);
c[++cc] = new TCanvas("e_ref", "e_ref", 0, 0, 720, 400);
int spectra_rebin = 0;
int yaxis_max = 0;
//~ int yaxis_max = 1200;
//~ int xaxis_max = 2500;
int xaxis_max = 3000;
//~ double fitcenter = 1190.;
//~ double fitw = 70.;
//~ double Cher_min = 10;
//~ double Cher_max = 16;
double Cher_min = -120;
double Cher_max = -110;
//~ double const_min = 20;
//~ double const_max = 80;
double const_min = -250;
double const_max = -200;
int ccccdi = 0;
if( strcmp(fname, "eff003")==0) {fitcenter=760.; fitw=100.;}
if( strcmp(fname, "eff004")==0) {fitcenter=1220.; fitw=150.; spectra_rebin=8; yaxis_max=3300;}
if( strcmp(fname, "eff005")==0) {fitcenter=1140.; fitw=60.;}
if( strcmp(fname, "eff006")==0) {fitcenter=2050.; fitw=60.; xaxis_max=4000;}
if( strcmp(fname, "eff012")==0) {Cher_min=13; Cher_max=15;}
if( strcmp(fname, "eff013")==0) {Cher_min=12; Cher_max=14;}
if( strcmp(fname, "eff014")==0) {Cher_min=12; Cher_max=14;}
if( strcmp(fname, "eff015")==0) {Cher_min=12; Cher_max=14;}
if( strcmp(fname, "eff016")==0) {Cher_min=12; Cher_max=14;}
if( strcmp(fname, "eff017")==0) {Cher_min=12; Cher_max=14;}
if( strcmp(fname, "eff018")==0) {Cher_min=12; Cher_max=14;}
sprintf(hname, "hadc%d", 0);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("; Charge [200 fC/bin];Counts");
if(spectra_rebin) hp1d->Rebin(spectra_rebin);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0,yaxis_max);
hp1d->Draw();
//~ printf("Photofraction ~ %lf\n", hp1d->Integral(1300,1700)/hp1d->Integral(0,1300));
f_fit = fg;
f_fit->SetParameters(hp1d->GetMaximum()*0.1, fitcenter, fitw);
f_fit->SetRange(fitcenter-fitw, fitcenter+fitw);
f_fit->SetLineWidth(3);
hp1d->Fit(f_fit, "0Q", "", fitcenter-fitw, fitcenter+fitw);
f_fit->DrawClone("LSAME");
double PP_1sigma_min = f_fit->GetParameter(1)-f_fit->GetParameter(2);
double PP_1sigma_max = f_fit->GetParameter(1)+f_fit->GetParameter(2);
printf("PPcenter, PPwidth = %.0lf, %.0lf\n", f_fit->GetParameter(1), f_fit->GetParameter(2));
printf("PP+/-1sigma = [%.0lf, %.0lf]\n", PP_1sigma_min, PP_1sigma_max);
printf("PP resolution = %.1lf\%% (FWHM = %.1lf\%%)\n", 100*f_fit->GetParameter(2)/f_fit->GetParameter(1), 235*f_fit->GetParameter(2)/f_fit->GetParameter(1));
// -------- 1275keV Compton estimation --------
// COMPTON EDGE
//~ TF1 *fcompt = new TF1("fcompt", "[2]*(1 - ConstantStep(((x-[0])/[1]) - (1275.*(2*1275./511./(1+2*1275./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");
//~ TF1 *fcomhi = new TF1("fcomhi", "gaus(0)+fcompt");
//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
//~ double comhi_range_low = 1670.;
//~ //double comhi_range_hi = 2130.;
//~ double comhi_range_hi = 2070.;
//~ fcomhi->SetParameters(2000., 1710., 50., 100, 2, 200);
//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
//~ hp1d->Fit(fcomhi, "0QR");
//~ //fcomhi->DrawClone("LSAME");
//~
//~ fcompt->SetNpx(300); fcompt->SetLineWidth(1); fcompt->SetLineColor(kBlue); fcompt->SetLineStyle(9);
//~ fcompt->SetParameter(0, fcomhi->GetParameter(3));
//~ fcompt->SetParameter(1, fcomhi->GetParameter(4));
//~ fcompt->SetParameter(2, fcomhi->GetParameter(5));
//~ fcompt->SetRange(0, comhi_range_hi);
//~ fcompt->DrawClone("LSAME");
//~
//~ double IComHi = fcompt->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
//~ printf("IComHi = %.0lf\n", IComHi);
// G LINEAR G
//~ TF1 *fcomhi = new TF1("fcomhi", "gaus(0)+pol0(3)+gaus(4)");
//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(kBlue);
//~ double comhi_range_low = 1700.;
//~ double comhi_range_hi = 2150.;
//~ fcomhi->SetParameters(2000., 1710., 50.,
//~ 150.,
//~ 200., 2150., 100.);
//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
//~ hp1d->Fit(fcomhi, "0QR");
//~ fcomhi->DrawClone("LSAME");
//~
//~ TF1 *fcom_lin = new TF1("fcom_lin", "pol0(0)");
//~ fcom_lin->SetNpx(300); fcom_lin->SetLineWidth(1); fcom_lin->SetLineColor(6);
//~ fcom_lin->SetParameter(0, fcomhi->GetParameter(3));
//~ fcom_lin->SetRange(0, 4000);
//~ fcom_lin->DrawClone("LSAME");
//~
//~ double IComHi = fcom_lin->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
// G LINEAR EXP
TF1 *fcomhi = new TF1("fcomhi", "gaus(0) + [3] + [4]*TMath::Exp((x-[5])/[6])");
fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
//double comhi_range_low = 1250.;
double comhi_range_low = fitcenter - 1.5*fitw;
//double comhi_range_hi = 2100.;
//~ fcomhi->SetParameters(2000., fitcenter, fitw,
//~ 170.,
//~ 85., 2000., 100.);
fcomhi->SetParameters(800., fitcenter, fitw,
40.,
5., comhi_range_hi, 100.);
//~ fcomhi->SetParLimits(2, 50., 2000.);
fcomhi->SetRange(comhi_range_low, comhi_range_hi);
hp1d->Fit(fcomhi, "0QR");
//~ fcomhi->SetRange(0, 4000);
fcomhi->DrawClone("LSAME");
TF1 *fcompt = new TF1("fcompt", "[0] + [1]*TMath::Exp((x-[2])/[3])");
fcompt->SetNpx(300); fcompt->SetLineWidth(1); fcompt->SetLineColor(kBlue); fcompt->SetLineStyle(9);
fcompt->SetParameter(0, fcomhi->GetParameter(3));
fcompt->SetParameter(1, fcomhi->GetParameter(4));
fcompt->SetParameter(2, fcomhi->GetParameter(5));
fcompt->SetParameter(3, fcomhi->GetParameter(6));
//~ fcompt->SetRange(1000, comhi_range_hi);
fcompt->SetRange(400, comhi_range_hi);
fcompt->DrawClone("LSAME");
//~
double IComHi = fcompt->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
printf("IComHi = %.0lf\n", IComHi);
//~
//~ // EXP
//~ TF1 *fcomhi = new TF1("fcomhi", "[0] + [1]*TMath::Exp((x-[2])/[3])");
//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
//~ double comhi_range_low = 1800.;
//~ double comhi_range_hi = 2000.;
//~ fcomhi->SetParameters(150., 1., 2000.,1.);
//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
//~ hp1d->Fit(fcomhi, "R");
//~ fcomhi->SetRange(0, 4000);
//~ fcomhi->DrawClone("LSAME");
//~
//~ double IComHi = fcomhi->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
//~ printf("IComHi = %.0lf\n", IComHi);
// LIN
//~ TF1 *fcomhi = new TF1("fcomhi", "pol0");
//~ fcomhi->SetNpx(500); fcomhi->SetLineWidth(1); fcomhi->SetLineColor(6);
//~ double comhi_range_low = 1800.;
//~ double comhi_range_hi = 1900.;
//~ fcomhi->SetParameter(0,150);
//~ fcomhi->SetRange(comhi_range_low, comhi_range_hi);
//~ hp1d->Fit(fcomhi, "QR");
//~ fcomhi->SetRange(0, 4000);
//~ fcomhi->DrawClone("LSAME");
//~
//~ double IComHi = fcomhi->Integral(PP_1sigma_min, PP_1sigma_max)/hp1d->GetBinWidth(1);
//~ printf("IComHi = %.0lf\n", IComHi);
// -------- 1275keV Compton estimation --------
sprintf(hname, "hadc_cut%d", 0);
hp1dcut = DrTH1F(dir, hname, "");
if(spectra_rebin) hp1dcut->Rebin(spectra_rebin);
hp1dcut->SetLineColor(kRed);
hp1dcut->SetLineColor(40);
hp1dcut->SetLineStyle(2);
hp1dcut->DrawClone("SAME");
double N511 = hp1dcut->GetEntries();
//~ double N511 = hp1dcut->Integral(PP_1sigma_min-10, PP_1sigma_max+10);
leg[++legi] = new TLegend(0.124302, 0.65, 0.6, 0.925134);
leg[legi]->SetFillColor(0); leg[legi]->SetBorderSize(1); leg[legi]->SetTextSize(0.05); leg[legi]->SetTextAlign(22);
leg[legi]->AddEntry((TH1F*)hp1d->Clone(), "Measurement", "L");
leg[legi]->AddEntry((TH1F*)hp1dcut->Clone(), "#pm1 #sigma Photopeak Cut", "L");
leg[legi]->AddEntry((TF1*)fcompt->Clone(), "1275 keV Compton Estimate", "L");
//~ leg[legi]->Draw();
//~ sprintf(hname, "hadc_cut_2%d", 0);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kGreen);
//~ hp1d->DrawClone("SAME");
sprintf(fullname, "gif/e_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/e_%s.eps", fname); c[cc]->SaveAs(fullname);
c[++cc] = new TCanvas("e_sipm", "e_sipm", 0, 470, 720, 300);
int ch_Sprttype = 8;
(c[cc]->cd(1))->SetLogy(1);
sprintf(hname, "htdc%d", ch_Sprttype);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("; Time [ns];Counts");
//~ hp1d->Rebin(4);
hp1d->SetLineColor(kBlack);
//~ hp1d->GetXaxis()->SetRangeUser(0,80);//original
hp1d->GetXaxis()->SetRangeUser(-400,50);
hp1d->GetYaxis()->SetRangeUser(0.1,hp1d->GetMaximum()*2);
hp1d->DrawClone();
sprintf(hname, "htdc_cut%d", ch_Sprttype);
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->Rebin(4);
hp1d->SetLineColor(kRed);
TF1 *fc = new TF1("fc", "pol0"); fc->SetNpx(300);
fc->SetParameter(0,1);
//~ fc->SetParameter(1,-0.1);
fc->SetRange(const_min,const_max);
hp1d->Fit(fc, "0QWWR");
hp1d->DrawClone("SAME");
fc->SetLineColor(kBlue);
fc->DrawClone("LSAME");
double Ntdc = hp1d->Integral(hp1d->FindBin(Cher_min), hp1d->FindBin(Cher_max));
double Nbckg = fc->Integral(Cher_min,Cher_max)/hp1d->GetBinWidth(1);
printf("fc par(0) = %lf | Nbckg = %.0lf\n", fc->GetParameter(0), Nbckg);
//~ sprintf(hname, "htdc_cut_2%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kGreen);
//~ hp1d->DrawClone("SAME");
printf("N511 = %.0lf | Ntdc = %.0lf | R = %.2lf \%% | R_1 = %.2lf \%% | R_2 = %.2lf \%%\n", N511, Ntdc, 100*Ntdc/N511, 100*(Ntdc-Nbckg)/N511, 100*(Ntdc-Nbckg)/(N511-IComHi));
sprintf(fullname, "gif/e2_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/e2_%s.eps", fname); c[cc]->SaveAs(fullname);
int kobayashi =1;
if(kobayashi){
double Nnohit = hp1d->Integral(hp1d->FindBin(-715), hp1d->FindBin(-685));
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));
}
//~ (c[cc]->cd(++ccccdi))->SetLogy(1);
//~ sprintf(hname, "hadc%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kBlack);
//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->DrawClone();
//~
//~ sprintf(hname, "hadc_cut%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kRed);
//~ hp1d->DrawClone("SAME");
//~
//~ sprintf(hname, "hadc_cut_2%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kGreen);
//~ hp1d->DrawClone("SAME");
//~ sprintf(fullname, "gif/e_%s.gif", fname); c[cc]->SaveAs(fullname);
//~ sprintf(fullname, "eps/e_%s.eps", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'r') != NULL ) {
//~ gStyle->SetOptStat(1111);
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);
c[++cc] = new TCanvas("e", "e", 0, 0, 1200, 450);
c[cc]->Divide(3,1);
int spectra_rebin = 0;
int yaxis_max = 0;
int xaxis_max = 3000;
double fitcenter = 700.;
double fitw = 60.;
double Cher_min = 4;
double Cher_max = 16;
double const_min = Cher_max;
double const_max = 80;
int ccccdi = 0;
if( strcmp(fname, "eff003")==0) {fitcenter=760.; fitw=100.;}
if( strcmp(fname, "eff004")==0) {fitcenter=1220.; fitw=150.; spectra_rebin=8; yaxis_max=3300;}
if( strcmp(fname, "eff005")==0) {fitcenter=1140.; fitw=60.;}
if( strcmp(fname, "eff006")==0) {fitcenter=2050.; fitw=60.; xaxis_max=4000;}
if( strcmp(fname, "eff007")==0) {fitcenter=710.; fitw=50.;}
if( strcmp(fname, "eff009")==0) {fitcenter=900.; fitw=60.; xaxis_max=2000;}
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", 0);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("Reference Coincidence; Charge [200 fC/bin];");
if(spectra_rebin) hp1d->Rebin(spectra_rebin);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0.5,yaxis_max);
f_fit = fg;
f_fit->SetParameters(hp1d->GetMaximum()*0.1, fitcenter, fitw);
f_fit->SetRange(fitcenter-fitw, fitcenter+fitw);
f_fit->SetLineWidth(1);
hp1d->Fit(f_fit, "Q", "", fitcenter-fitw, fitcenter+fitw);
printf("PP+/-1sigma = [%.0lf, %.0lf]\n", f_fit->GetParameter(1)-f_fit->GetParameter(2), f_fit->GetParameter(1)+f_fit->GetParameter(2));
printf("PP resolution = %.1lf\%% (FWHM = %.1lf\%%)\n", 100*f_fit->GetParameter(2)/f_fit->GetParameter(1), 235*f_fit->GetParameter(2)/f_fit->GetParameter(1));
sprintf(hname, "hadc_cut%d", 0);
hp1d = DrTH1F(dir, hname, "");
if(spectra_rebin) hp1d->Rebin(spectra_rebin);
hp1d->SetLineColor(kRed);
hp1d->DrawClone("SAME");
double N511 = hp1d->GetEntries();
sprintf(hname, "hadc_cut_2%d", 0);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetLineColor(kGreen);
hp1d->DrawClone("SAME");
//-----------------------------------------------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "hadc%d", 2);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("Reference 1275 keV; Charge [200 fC/bin];");
//~ if(spectra_rebin) hp1d->Rebin(spectra_rebin);
hp1d->SetLineColor(kBlack);
//~ hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
//~ if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0,yaxis_max);
hp1d->DrawClone();
sprintf(hname, "hadc_cut%d", 2);
hp1d = DrTH1F(dir, hname, "");
if(spectra_rebin) hp1d->Rebin(spectra_rebin);
hp1d->SetLineColor(kRed);
hp1d->DrawClone("SAME");
sprintf(hname, "hadc_cut_2%d", 2);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetLineColor(kGreen);
hp1d->DrawClone("SAME");
//-----------------------------------------------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", 1);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("MPPC Cherenkov; Time [ns];");
//~ hp1d->Rebin(4);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(0,80);
hp1d->GetYaxis()->SetRangeUser(0.1,hp1d->GetMaximum()*2);
hp1d->DrawClone();
sprintf(hname, "htdc_cut%d", 1);
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->Rebin(4);
hp1d->SetLineColor(kRed);
TF1 *fc = new TF1("fc", "pol0"); fc->SetNpx(300);
fc->SetParameter(0,1);
fc->SetRange(const_min,const_max);
hp1d->Fit(fc, "0QWWR");
hp1d->DrawClone("SAME");
fc->SetLineColor(kBlue);
fc->DrawClone("LSAME");
double Ntdc = hp1d->Integral(hp1d->FindBin(Cher_min), hp1d->FindBin(Cher_max));
double Nbckg = fc->Integral(Cher_min,Cher_max)/hp1d->GetBinWidth(1);
printf("fc par(0) = %lf | Nbckg = %.0lf\n", fc->GetParameter(0), Nbckg);
sprintf(hname, "htdc_cut_2%d", 1);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetLineColor(kGreen);
hp1d->DrawClone("SAME");
printf("N511 = %.0lf | Ntdc = %.0lf | R = %lf\n", N511, Ntdc, Ntdc/N511);
printf("R_1 = %lf\n", (Ntdc-Nbckg)/N511);
//~ (c[cc]->cd(++ccccdi))->SetLogy(1);
//~ sprintf(hname, "hadc%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kBlack);
//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->DrawClone();
//~
//~ sprintf(hname, "hadc_cut%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kRed);
//~ hp1d->DrawClone("SAME");
//~
//~ sprintf(hname, "hadc_cut_2%d", 1);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kGreen);
//~ hp1d->DrawClone("SAME");
sprintf(fullname, "gif/e_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/e_%s.eps", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------
if( strchr(plopt, 'p') != NULL ) {
DrSetDrawStyle(0.05);
gStyle->SetOptStat("e");
gStyle->SetOptFit(0);
gStyle->SetPadRightMargin(0.03);
gStyle->SetPadTopMargin(0.05);
gStyle->SetPadBottomMargin(0.12);
gStyle->SetStatFontSize(0.05);
gStyle->SetStatX(0.97);
gStyle->SetStatY(0.95);
c[++cc] = new TCanvas("Print", "Print", 0, 0, 720, 400);
int spectra_rebin = 0;
int yaxis_max = 1000;
int xaxis_max = 1200;
double Cher_min = 10;
double Cher_max = 16;
double const_min = 20;
double const_max = 80;
int ccccdi = 0;
sprintf(hname, "hadc%d", 0);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("; Charge [200 fC/bin];Counts");
if(spectra_rebin) hp1d->Rebin(spectra_rebin);
hp1d->SetLineColor(kBlack);
hp1d->GetXaxis()->SetRangeUser(0,xaxis_max);
if(yaxis_max) hp1d->GetYaxis()->SetRangeUser(0,yaxis_max);
hp1d->Draw();
sprintf(hname, "hadc_cut%d", 0);
hp1dcut = DrTH1F(dir, hname, "");
if(spectra_rebin) hp1dcut->Rebin(spectra_rebin);
hp1dcut->SetLineColor(kRed);
hp1dcut->Draw("SAME");
sprintf(fullname, "gif/p_%s.gif", fname); c[cc]->SaveAs(fullname);
sprintf(fullname, "eps/p_%s.eps", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'b') != NULL ) {
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);
c[++cc] = new TCanvas("b", "b", 0, 0, 1200, 700);
c[cc]->Divide(2,2);
int ccccdi = 0;
int ich;
/*
// ------------------------------------------------------------
ich = 0;
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle(hname);
//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
// ------------------------------------------------------------
ich = 1;
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle(hname);
//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
*/
// -------------------------- qdc 0 ----------------------------------
ich = 0;
(c[cc]->cd(++ccccdi))->SetLogy(1);
/*
sprintf(hname, "hadc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("SiPM;Charge [200 fC/bin]");
hp1d->GetXaxis()->SetRangeUser(0,1000);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
*/
TLegend *leg0 = new TLegend(0.57,0.66,0.88,0.88);
leg0->SetFillColor(0);
leg0->SetBorderSize(1);
leg0->SetTextSize(0.05);
leg0->SetTextAlign(22);
sprintf(hname, "hadc%d", ich);
TH1F* hADC = DrTH1F(dir, hname, "");
hADC->SetTitle("SiPM;Charge [200 fC/bin]");
hADC->GetXaxis()->SetRangeUser(0,1000);
hADC->SetLineColor(kBlack);
hADC->DrawClone();
int dbgprint=1;
const int nPeaks=3;
double pSigma=10;
double pThr=0.005;
TSpectrum *sADC = new TSpectrum(nPeaks,0.3);
sADC->Search(hADC, pSigma,"",pThr);
int foundNPeaks = sADC->GetNPeaks();
float *foundPeakX = sADC->GetPositionX();
float *foundPeakY = sADC->GetPositionY();
if(1 < foundNPeaks) {
//void Sort(Int_t n1, const Float_t *a, Int_t *index, Bool_t down)
int sortedPeakIndex[nPeaks];
TMath::Sort(foundNPeaks, foundPeakX, sortedPeakIndex, 0);
if(dbgprint) {
printf("foundNPeaks = %d\n",foundNPeaks);
for(int i=0; i<foundNPeaks; i++) printf("Peak[%d] = %f\n",i, foundPeakX[sortedPeakIndex[i]]);
}
float peakMidpointX[nPeaks+1];
float peakMidpointY[nPeaks+1];
peakMidpointX[0] = 0; peakMidpointY[0] = 1;
for(int i=1; i<foundNPeaks; i++) {
peakMidpointX[i] = (foundPeakX[sortedPeakIndex[i]] + foundPeakX[sortedPeakIndex[i-1]])/2.0;
peakMidpointY[i] = 1;
if(dbgprint) printf("peakMidpointX[%d] = %f\n",i, peakMidpointX[i]);
}
peakMidpointX[foundNPeaks] = 1120; peakMidpointY[foundNPeaks] = 1;
if(dbgprint) printf("peakMidpointX[%d] = %f\n",foundNPeaks, peakMidpointX[foundNPeaks]);
TPolyMarker *pmMidpoints = new TPolyMarker(nPeaks+1, peakMidpointX, peakMidpointY);
pmMidpoints->SetMarkerStyle(33); pmMidpoints->SetMarkerColor(9); pmMidpoints->SetMarkerSize(2);
pmMidpoints->Draw("SAME");
//~ fprintf(fp,"peakEvents:\n");
double peakEvents[nPeaks];
for(int i=0; i<foundNPeaks-1; i++) {
peakEvents[i] = hADC->Integral(hADC->FindBin(peakMidpointX[i]), hADC->FindBin(peakMidpointX[i+1]));
if(dbgprint) printf("I(%.1f, %.1f) = %.1lf\n", peakMidpointX[i], peakMidpointX[i+1], peakEvents[i]);
//~ fprintf(fp,"%lf\n", peakEvents[i]);
}
peakEvents[foundNPeaks-1] = hADC->Integral(hADC->FindBin(peakMidpointX[foundNPeaks-1]), hADC->FindBin(1130.));
if(dbgprint) printf("I(%.1f, %.1f) = %.1lf\n", peakMidpointX[foundNPeaks-1], 1130., peakEvents[foundNPeaks-1]);
//~ fprintf(fp,"%lf\n", peakEvents[foundNPeaks-1]);
double N0all = -TMath::Log(peakEvents[0]/hADC->GetEntries());
double N01 = peakEvents[1]/peakEvents[0];
double N02 = 2*peakEvents[2]/peakEvents[1];
if(dbgprint) printf("Poisson <N>: N0all -> %.3lf | N01 -> %.3lf | N12 -> %.3lf\n", N0all, N01, N02);
sprintf(sbuff, "<N> = %.3lf", N0all);
leg0->AddEntry((TH1F*)hADC->Clone(),sbuff, "");
sprintf(sbuff, "1 CPH = %.0lf", foundPeakX[sortedPeakIndex[1]]);
leg0->AddEntry((TH1F*)hADC->Clone(),sbuff, "");
leg0->Draw();
}
double N_adc_integral = hADC->Integral(hADC->FindBin(peakMidpointX[1]), hADC->FindBin(1000.));
double N_poisson = (1 - TMath::Exp(-N0all))*hADC->GetEntries();
printf("hADC->GetEntries() = %.0lf\n", hADC->GetEntries());
printf("N_adc_integral = %.0lf | N_poisson = %.0lf\n", N_adc_integral, N_poisson);
sprintf(hname, "hadc_cut%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetLineColor(kRed);
hp1d->DrawClone("SAME");
// -------------------------- TDC 0 ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(1);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
hp1d->SetTitle("SiPM;Time [ns]");
hp1d->GetXaxis()->SetRangeUser(-5,5);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
hp1d->DrawClone();
printf("hTDC->GetEntries() = %.0lf\n", hp1d->GetEntries());
double N_tdc_integral = hp1d->Integral(hp1d->FindBin(-10), hp1d->FindBin(10));
printf("N_tdc_integral = %.0lf\n", N_tdc_integral);
// -----------------------------------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(0);
sprintf(hname, "hcor%d", ich);
hp2d = DrTH2F(dir, hname, "");
hp2d->SetTitle(";Charge [200 fC/bin];cTime [ns]");
//~ hp1d->GetXaxis()->SetRangeUser(0,2000);
//~ hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
//~ hp1d->SetLineColor(kBlack);
hp2d->DrawClone("COLZ");
// -------------------------- cTDC 0 ----------------------------------
(c[cc]->cd(++ccccdi))->SetLogy(0);
sprintf(hname, "hctdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetTitle("SiPM;cTime [ns]");
//~ hp1d->GetXaxis()->SetRangeUser(-5,5);
//~ //hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
//~ hp1d->DrawClone();
plotCorTDC(hp1d, "SiPM;cTime [ns]", -0.5, 0.45, -1., 2);
printf("hcTDC->GetEntries() = %.0lf\n", hp1d->GetEntries());
double N_ctdc_integral = hp1d->Integral(hp1d->FindBin(-2), hp1d->FindBin(2));
printf("N_ctdc_integral = %.0lf\n", N_ctdc_integral);
printf("N_ctdc_integral/N_poisson = %.3lf\n", N_ctdc_integral/N_poisson);
printf("N_adc_integral/N_poisson = %.3lf\n", N_adc_integral/N_poisson);
printf("N_ctdc_integral/N_adc_integral = %.3lf\n", N_ctdc_integral/N_adc_integral);
//~ sprintf(hname, "htdc_cut%d", ich);
//~ hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetLineColor(kRed);
//~ hp1d->DrawClone("SAME");
sprintf(fullname, "gif/b_%s.gif", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------
if( strchr(plopt, 'g') != NULL ) {
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);
c[++cc] = new TCanvas("g", "g", 640, 0, 640, 360);
int ccccdi = 0;
int ich = 15;
(c[cc]->cd(++ccccdi))->SetLogy(0);
sprintf(hname, "htdc%d", ich);
hp1d = DrTH1F(dir, hname, "");
//~ hp1d->SetTitle("SiPM;cTime [ns]");
//~ hp1d->GetXaxis()->SetRangeUser(-5,5);
//~ //hp1d->GetYaxis()->SetRangeUser(0.9,hp1d->GetMaximum()*1.1);
hp1d->SetLineColor(kBlack);
//~ hp1d->DrawClone();
sprintf(hname, "Trigger(ch.%d)-Trigger(ch.31);Time [ns]", ich);
plotCorTDC(hp1d, hname, -0.25, 0.25, -0.5, 0.5);
sprintf(fullname, "gif/g_%s.gif", fname); c[cc]->SaveAs(fullname);
}
// -------------------------------------------------------------------------------
}