Rev 38 | Rev 50 | 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"
#include "TH3F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "stdio.h"
#include "THStack.h"
#include "TPaveText.h"
#include "include/RTUtil.h"
double getNoise(TH2F*, int, int);
int biasScan(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, double par1=0.0, double par2=0.0, bool debug = false)
{
const int c_nChannels = 64;
const double c_xOffset = 0; // mm
const double c_yOffset = 0;
int map[8][8]={{32,34,53,55,40,42,61,63},
{48,50,37,39,56,58,45,47},
{33,35,52,54,41,43,60,62},
{49,51,36,38,57,59,44,46},
{17,19,4,6,25,27,12,14},
{1,3,20,22,9,11,28,30},
{16,18,5,7,24,26,13,15},
{0,2,21,23,8,10,29,31}
};
// Set draw style
RTSetStyle(gStyle);
// open the file with histograms
char fnameroot[256];
TFile* rootfile;
sprintf(fnameroot, "root/%s.root", filename);
rootfile = (TFile *) gROOT->FindObject(filename);
if(rootfile==NULL) rootfile = new TFile(fnameroot);
if(rootfile==NULL) {
printf("Cannot open root file %s!!!\n",fnameroot);
return(0);
}
if(strstr(plopt, "tdc") != NULL) {
TCanvas *canvas21 = new TCanvas("canvas21","canvas21",1600,800);
TH3F* h0 = (TH3F*) rootfile->Get("h_correctedTDC");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Int_t binsZ = h0->GetZaxis()->GetNbins();
Int_t minZ = h0->GetZaxis()->GetFirst();
Int_t maxZ = h0->GetZaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
canvas21->Divide(2);
//h0->Draw();
canvas21->cd(1);
//h0->Project3D("xy")->Draw("colz");
//canvas21->cd(2);
//h0->Project3D("yz")->Draw("colz");
//canvas21->cd(3);
TH2D* h_correctedTDC = (TH2D*) h0->Project3D("xz");
//h0->GetZaxis()->SetRangeUser(-5,5);
//h0->Project3D("xzo")->Draw("colz");
//h_correctedTDC->SetTitle("; t [ns]; Channel");
h_correctedTDC->Draw("colz");
/*
TH1F* tdc1 = new TH1F("tdc1",";TDC [ns];Events",binsZ, minZ, maxZ);
//for(int j=minY; j<maxY; j++) {
for(int k=minZ; k<maxZ; k++) {
double signal = h0->GetBinContent(1,1,k);
tdc1->Fill(k,signal);
//}
}
tdc1->Draw();
*/
TH2D* h_timeWalk = (TH2D*) h0->Project3D("zy");
canvas21->cd(2);
//h_timeWalk->SetTitle(";Threshold [V]; t [ns]");
h_timeWalk->Draw("colz");
}
if (strstr(plopt, "bias") != NULL) {
TCanvas* canvas01 = new TCanvas("canvas01","",800, 800);
canvas01->cd();
TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");
h_bias->Draw("colz");
}
if (strstr(plopt, "th") != NULL) {
TCanvas* canvas1 = new TCanvas("canvas1","",800,800);
canvas1->cd();
//gPad->SetLogz();
//canvas1->SetLogz();
TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
h_threshold->SetTitleOffset(1.5,"y");
//h_threshold->GetZaxis()->SetRangeUser(0,1000);
h_threshold->SetContour(20);
h_threshold->Draw("colz");
TCanvas* canvas2 = new TCanvas("canvas2","",1600,800);
canvas2->Divide(2);
TH1D* h_projection1 = h_threshold->ProjectionY("Ch 37",38,38);
canvas2->cd(1);
//h_projection1->GetYaxis()->SetRangeUser(0,1000);
//gPad->SetLogy();
h_projection1->Draw();
TH1D* h_projection2 = h_threshold->ProjectionY("Ch 38",39,39);
canvas2->cd(2);
//gPad->SetLogy();
h_projection2->Draw();
}
if(strstr(plopt,"16") != NULL)
{
TH2F* h_threshold = (TH2F*) rootfile->Get("h_bias");
TCanvas* canvas = new TCanvas("canvas","",1600,1600);
TCanvas* canvas11 = new TCanvas("canvas11","",800,800);
canvas->cd();
TVirtualPad *main = new TPad("main","main",0,0,1,1,10,1);
main->Draw();
main->cd();
main->Divide(4,4);
for(int i=chXstart; i<chXend; i++) {
for(int j=(int)par1; j<(int)par2; j++) {
int channel = map[i][j];
TH1D* h_projection = h_threshold->ProjectionY("",channel+1,channel+1);
int canvasPosition = i-chXend+4*((int)par2-j)+1;
printf(" %d ", canvasPosition);
main->cd(canvasPosition);
gPad->SetMargin(0.08, 0.08, 0.08, 0.08);
char name[128];
sprintf(name,"Channel %d",channel);
h_projection->SetTitle(name);
h_projection->GetYaxis()->SetRangeUser(0,10000);
h_projection->DrawCopy("pe1x0");
canvas11->cd();
if (canvasPosition == 1) {
h_projection->SetMarkerStyle(kFullDotMedium);
h_projection->GetYaxis()->SetRangeUser(0,10000);
h_projection->DrawCopy("Pl");
}
else {
h_projection->SetMarkerStyle(kFullDotMedium);
h_projection->SetLineColor(canvasPosition);
h_projection->SetMarkerColor(canvasPosition);
h_projection->SetTitle("");
h_projection->DrawCopy("plsame");
}
}
}
printf("\n");
}
if(strstr(plopt, "1ch") != NULL)
{
TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");
TCanvas* canvas = new TCanvas("canvas","",800,800);
canvas->cd();
int channel = map[chXstart][chXend];
TH1D* h_projection = h_bias->ProjectionY("",channel+1,channel+1);
char name[128];
sprintf(name,"Channel %d",channel);
h_projection->SetTitle(name);
h_projection->SetMarkerStyle(kFullDotMedium);
//h_projection->SetMarkerSize(8);
gPad->SetMargin(0.1,0.1,0.1,0.1);
TF1* f_linear = new TF1("f_linear", "[0] + [1]*x", par1, par2);
f_linear->SetParameter(0, -1);
f_linear->SetParameter(1, 0.1);
f_linear->SetParNames("p0", "p1");
h_projection->Fit(f_linear,"q");
h_projection->Fit(f_linear,"r");
h_projection->Draw("pe1x0");
double p0 = f_linear->GetParameter(0);
double p1 = f_linear->GetParameter(1);
printf("p0 %f p1 %f V_b = %f [V]\n", p0, p1, -p0/p1);
TPaveText *text = new TPaveText(0.15,0.5,0.35,0.56);
text->SetFillColor(0);
text->SetBorderSize(0.1);
text->SetTextSize(0.030);
char string[128];
sprintf(string, "V_b = %2.2f", -p0/p1);
text->AddText(string);
text->Draw("same");
text->Paint();
}
if( strstr(plopt, "all") != NULL ) {
TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);
TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);
canvas2->Divide(8,8);
canvas3->Divide(8,8);
TH1F* h_hitsx;
TH1F* h_hitsy;
for(int i=0; i<c_nChannels; i++) {
canvas2->cd(i+1);
char hname[128];
sprintf(hname, "hnhitsx%d", i);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->Draw();
canvas3->cd(i+1);
sprintf(hname, "hnhitsy%d", i);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}
}
if( strstr(plopt, "x") != NULL ) {
TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);
TH1F* h_hitsx;
canvas10->cd();
char hname[128];
sprintf(hname, "hnhitsx%d", chXstart);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->Draw("colz");
}
if( strstr(plopt, "y") != NULL ) {
TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);
TH1F* h_hitsy;
canvas11->cd();
char hname[128];
sprintf(hname, "hnhitsy%d", chXstart);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}
/** Draws the sum of the channels
* Each channel is a 2d plot
* Intended for the study of 1 channel
*/
if (strstr(plopt, "2d") != NULL) {
int nX = chXend - chXstart + 1;
int nY = (int)par2 - (int)par1 + 1;
TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX*400,nY*400);
printf("nx %d ny %d\n",nX,nY);
canvas7->Divide(nX,nY);
for(int i=(int)par1; i<=(int)par2; i++) {
for(int j=chXstart; j<=chXend; j++) {
//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
int canvasPosition = nX*((int)par2-i)+(j-chXstart)+1;
if (debug) printf("canvas %d\n",canvasPosition);
canvas7->cd(canvasPosition);
char hname[128];
int chPosition = map[j][i];
sprintf(hname, "h2d%d", chPosition);
TH2F* h_2d = (TH2F*)rootfile->Get(hname);
h_2d->Draw("colz");
} //x
}
// Number of photoelectrons - Poissonian correction
TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);
canvas8->cd();
gStyle->SetOptStat(0);
char hname[128];
int chPosition = map[1][2];
sprintf(hname, "h2d%d", chPosition);
TH2F* h_2d = (TH2F*)rootfile->Get(hname);
Int_t binsX = h_2d->GetXaxis()->GetNbins();
Int_t minX = h_2d->GetXaxis()->GetFirst();
Int_t maxX = h_2d->GetXaxis()->GetLast()+1;
Int_t binsY = h_2d->GetYaxis()->GetNbins();
Int_t minY = h_2d->GetYaxis()->GetFirst();
Int_t maxY = h_2d->GetYaxis()->GetLast()+1;
Double_t xLowUser = h_2d->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h_2d->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h_2d->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h_2d->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);
double noise = getNoise(h_2d, 1, 89);
if(debug) printf("Noise = %f\n", noise);
for(int k=minX; k<=maxX; k++) {
for(int j=minY; j<=maxY; j++) {
double signal = h_2d->GetBinContent(k,j); // detected
//double p = ((signal - noise) > 1) ? (signal-noise) : 1;
double p = signal - noise;
p /= 10000.;
double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
double eta = -log(p0); // constant of the poissonian statistics
if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
double x = k*(xUpUser-xLowUser)/double(binsX);
double y = j*(yUpUser-yLowUser)/double(binsY);
h_corrected->Fill(x+3,y+8, eta);
}
}
h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");
h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);
h_corrected->Draw("colz");
// collection efficiency
int nPoints =0;
double efficiency=0;
for (int i=18; i<=58; i++) {
for (int j=19; j<=59; j++) {
double signal = h_corrected->GetBinContent(i,j);
if(debug) printf("signal %f\n",signal);
efficiency += signal;
nPoints++;
}
}
printf("Signal sum = %f\n # of points = %d\n",efficiency,nPoints);
}
/** Draws the sum of channel signals
* Each channel is a 2d ('h2d') histogram
* Suitable for 8x8 chs scan
*/
if( strstr(plopt, "sum") != NULL ) {
int nX = 7 - 0 + 1;
int nY = 7 - 0 + 1;
TCanvas* canvas12 = new TCanvas("canvas12","c2",8*200, 8*200);
canvas12->cd();
gStyle->SetOptStat(0);
// final histogram parameters
TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
TH2F* h[9];
// 2d histogram noise subtraction and poisson scaling
for(int i=0; i<=7; i++) {
for(int j=chXstart; j<=chXend; j++) {
int chPosition = map[j][i];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = nX*(i-0)+(j-chXstart);
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 170);
if (debug) printf("noise: %d\n",noise);
for(int k=minX; k<=maxX; k++) {
for(int l=minY; l<=maxY; l++) {
int signal = h[histogram]->GetBinContent(k,l); // detected
//double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;
double p = signal - noise;
p /= 10000.;
double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
double eta = -log(p0);
//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
double x = k*(xUpUser-xLowUser)/double(binsX);
double y = l*(yUpUser-yLowUser)/double(binsY);
h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);
}
}
}
}
h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");
//h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);
h_corrected->Draw("colz");
TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);
canvas13->Divide(2);
canvas13->cd(1);
//h[16]->Draw("colz");
canvas13->cd(2);
//h[8]->Draw("colz");
}
return(0);
}
/** Function calculates the noise from one channel
* it runs through the bins along x and returns the average value
*/
double getNoise(TH2F* histogram, int yStart, int yEnd)
{
double noise=0;
int count=0;
for(int j=yStart; j<yEnd; j++) {
double value = histogram->GetBinContent(j,2);
//if (noise < value) noise = value;
noise += value;
count++;
}
return (noise/double(count));
}