Rev 40 | 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; // mmconst 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 styleRTSetStyle(gStyle);// open the file with histogramschar 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");TCanvas* canvas02 = new TCanvas("canvas02","",800, 800);//canvas02->Divide(2,2);canvas02->cd();TH3F* h_biasThreshold = (TH3F*) rootfile->Get("h_biasThreshold");//h_biasThreshold->Draw();//canvas02->cd(2);//h_biasThreshold->GetXaxis()->SetRange(1,5);h_biasThreshold->Project3D("zy");h_biasThreshold_zy->DrawCopy("colz");TCanvas* canvas03 = new TCanvas("canvas03","",1200, 1200);canvas03->Divide(4,4);int x1=0;int x2=3;int y1=0;int y2=3;for(int j=x1; j<=x2; j++){for(int k=y1; k<=y2; k++){int channel = map[j][k];int canvasPosition = 4*(4-k-1) + 4 - j;canvas03->cd(canvasPosition);h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);char histoName[128];sprintf(histoName, "Ch %d", channel);h_biasThreshold->Project3D("zy");//sprintf(option, "h_biasThresholsd_zy%f", channel);h_biasThreshold_zy->SetTitle(histoName);h_biasThreshold_zy->DrawCopy("colz");}}TCanvas* canvas04 = new TCanvas("canvas04","",1200, 1200);canvas04->Divide(4,4);int x1=0;int x2=3;int y1=4;int y2=7;for(int j=x1; j<=x2; j++){for(int k=y1; k<=y2; k++){int channel = map[j][k];int canvasPosition = 4*(8-k-1) + 4 - j;canvas04->cd(canvasPosition);h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);char histoName[128];sprintf(histoName, "Ch %d", channel);h_biasThreshold->Project3D("zy");//sprintf(option, "h_biasThresholsd_zy%f", channel);h_biasThreshold_zy->SetTitle(histoName);h_biasThreshold_zy->DrawCopy("colz");}}TCanvas* canvas05 = new TCanvas("canvas05","",1200, 1200);canvas05->Divide(4,4);int x1=4;int x2=7;int y1=0;int y2=3;for(int j=x1; j<=x2; j++){for(int k=y1; k<=y2; k++){int channel = map[j][k];int canvasPosition = 4*(4-k-1) + 8 - j;canvas05->cd(canvasPosition);h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);char histoName[128];sprintf(histoName, "Ch %d", channel);h_biasThreshold->Project3D("zy");//sprintf(option, "h_biasThresholsd_zy%f", channel);h_biasThreshold_zy->SetTitle(histoName);h_biasThreshold_zy->DrawCopy("colz");}}TCanvas* canvas06 = new TCanvas("canvas06","",1200, 1200);canvas06->Divide(4,4);int x1=4;int x2=7;int y1=4;int y2=7;for(int j=x1; j<=x2; j++){for(int k=y1; k<=y2; k++){int channel = map[j][k];int canvasPosition = 4*(8-k-1) + 8 - j;canvas06->cd(canvasPosition);h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);char histoName[128];sprintf(histoName, "Ch %d", channel);h_biasThreshold->Project3D("zy");//sprintf(option, "h_biasThresholsd_zy%f", channel);h_biasThreshold_zy->SetTitle(histoName);h_biasThreshold_zy->DrawCopy("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 correctionTCanvas* 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 statisticsif (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 efficiencyint 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 parametersTH2F* 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 scalingfor(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));}