Rev 50 | Blame | Compare with Previous | Last modification | View Log | RSS feed
// *****************************************************************************// Jozef Stefan Institute// Dino Tahirovic// 2012-2014//// Options:// 2d - 1 channel analysis// s - save an image (png) after drawinv the canvas// p - save every canvas as page in pdf//// *****************************************************************************#include "TROOT.h"#include "TFile.h"#include "TBenchmark.h"#include "TH1F.h"#include "TH2F.h"#include "TCanvas.h"#include "TStyle.h"#include "TPad.h"#include "TPaveText.h"//#include "TLabel.h"#include "TF1.h"#include "TGraph.h"#include "TSpectrum.h"#include "stdio.h"#include "include/RTUtil.h"double getNoise(TH2F*, int, int);int sipm(char filename[256] = "test", char plopt[256]="all", int parameter1=0, int parameter2=7, int chYstart=0, int chYend=7, bool debug = false){const int c_nChannels = 64;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}};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);}// set draw styleRTSetStyle(gStyle);// Print results to .pdf (option "p")if(strstr(plopt, "p") != NULL) {TCanvas* canvasFirst = new TCanvas("canvasFirst","",500,500);TPaveText* text = new TPaveText(.05,.1,.95,.8);text->AddText("SiPM Scan results");canvasFirst->cd();text->Draw();canvasFirst->Print("result.pdf(");}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();int chPosition = map[parameter1][parameter2];char hname[128];sprintf(hname, "hnhitsx%d", chPosition);h_hitsx = (TH1F*)rootfile->Get(hname);h_hitsx->Draw();if(strstr(plopt, "p") != NULL) {canvas10->Print("result.pdf");}}if( strstr(plopt, "y") != NULL ) {TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);TH1F* h_hitsy;canvas11->cd();int chPosition = map[parameter1][parameter2];char hname[128];sprintf(hname, "hnhitsy%d", chPosition);h_hitsy = (TH1F*)rootfile->Get(hname);h_hitsy->Draw();if(strstr(plopt, "p") != NULL) {canvas11->Print("result.pdf");}}if( strstr(plopt, "share") != NULL ) {/*TCanvas *canvas4 = new TCanvas("canvas1","canvas1",1000,1000);int nChannels = chYend-chYstart+1;int ncols = nChannels/2;printf("nch %d nch\\2 %d\n", nChannels, ncols);canvas4->Divide(2,ncols);TH1F* h_hitsy;for(int i=chYstart; i<=chYend; i++){canvas4->cd(i-chYstart+1);char hname[128];int chPosition = map[0][i];sprintf(hname, "hnhitsy%d", chPosition);h_hitsy = (TH1F*)rootfile->Get(hname);h_hitsy->Draw();}*/TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);canvas4->cd();for(int i=parameter1; i<=parameter2; i++) {TH1F* h_hitsx;char hname[128];int chPosition = map[i][chYstart];sprintf(hname, "hnhitsx%d", chPosition);h_hitsx = (TH1F*)rootfile->Get(hname);h_hitsx->SetTitle("Scan X;x [mm]; Entries");h_hitsx->GetYaxis()->SetTitleOffset(1.3);h_hitsx->SetStats(0);if (i == parameter1)h_hitsx->Draw();else {h_hitsx->SetLineColor(i+1);h_hitsx->Draw("same");}}//sprintf(fullname, "ps/%s_Yshare.eps", filename);//canvas4->SaveAs(fullname);TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);canvas5->cd();for(int i=chYstart; i<=chYend; i++) {TH1F* h_hitsy;char hname[128];int chPosition = map[parameter1][i];sprintf(hname, "hnhitsy%d", chPosition);h_hitsy = (TH1F*)rootfile->Get(hname);h_hitsy->SetTitle("Scan Y;y [mm]; Entries");h_hitsy->GetYaxis()->SetTitleOffset(1.3);h_hitsy->SetStats(0);if (i == chYstart)h_hitsy->Draw();else {h_hitsy->SetLineColor(i+1);h_hitsy->Draw("same");}}//sprintf(fullname, "ps/%s_Yshare.eps", filename);//canvas5->SaveAs(fullname);}/** Draw one channel in x and y* Without and with the noise subtraction*/if (strstr(plopt, "noise") != NULL) {TCanvas *canvas51 = new TCanvas("canvas51","canvas51",500,500);canvas51->cd();// Get the filled histogram from scanTH1F* h_x;char name[32];sprintf(name, "hnhitsx%d", parameter1);h_x = (TH1F*)rootfile->Get(name);h_x->DrawCopy();TH1F* h_y;sprintf(name, "hnhitsy%d", parameter1);h_y = (TH1F*)rootfile->Get(name);// Create and fill corrected histogramInt_t binsX = h_x->GetXaxis()->GetNbins();if(debug) printf("nBins: %d\n", binsX);//Double_t xLowUser = h_x->GetXaxis()->GetBinLowEdge(1);//Double_t xUpUser = h_x->GetXaxis()->GetBinUpEdge(binsX);Double_t xLowUser = h_x->GetXaxis()->GetXmin();Double_t xUpUser = h_x->GetXaxis()->GetXmax();TH1F* h_xCorrected = new TH1F("h_xCorrected",";;",binsX, xLowUser, xUpUser);if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);double noise = h_x->GetBinContent(2);for (int jk = 0; jk < binsX; jk++) {double signal = h_x->GetBinContent(jk+1);signal -= noise;double x = h_x->GetXaxis()->GetBinCenter(jk+1);h_xCorrected->Fill(x, signal);}h_xCorrected->SetLineColor(2);h_xCorrected->DrawCopy("same");//Create and fill corrected histogramInt_t binsY = h_y->GetXaxis()->GetNbins();Double_t yLowUser = h_y->GetXaxis()->GetXmin();Double_t yUpUser = h_y->GetXaxis()->GetXmax();TH1F* h_yCorrected = new TH1F("h_yCorrected",";;",binsY, yLowUser, yUpUser);noise = h_y->GetBinContent(2);for (int jk = 0; jk < binsY; jk++) {double signal = h_y->GetBinContent(jk+1);signal -= noise;double x = h_y->GetXaxis()->GetBinCenter(jk+1);h_yCorrected->Fill(x, signal);}TCanvas *canvas52 = new TCanvas("canvas52","canvas52",500,500);canvas52->cd();h_y->DrawCopy();h_yCorrected->SetLineColor(2);h_yCorrected->DrawCopy("same");}/** Draws the signal from 8 channels in x-row* for one specific y bin, so the background and cross-talk* can be estimated.* Draws also a 2d scan of these channels.*/if (strstr(plopt, "line") != NULL) {TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);canvas6->cd(0);gStyle->SetOptStat(0);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);//! 1-dimension position in x vs. hitsTH2F* h[8];TH1F* h_line[8];for(int j=0; j<8; j++) {h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);}for(int j=parameter1; j<=parameter2; j++) {int chPosition = map[j][chYstart];char hname[128];sprintf(hname, "h2d%d", chPosition);int histogram = j;h[histogram] = (TH2F *) rootfile->Get(hname);int noise = getNoise(h[histogram], 1, 160);for(int k=minX; k<=maxX; k++) {int l=chYstart*20+12;//for(int l=12; l<=16; l++) {double signal = h[histogram]->GetBinContent(k,l);//signal -= noise;//signal /= 5*10000.0;double eta = -log(1 - signal);double x = xLowUser + k*(xUpUser-xLowUser)/double(binsX);//double y = l*(yUpUser-yLowUser)/double(binsY);h_line[j]->Fill(x, signal);//}}if (j == parameter1) {h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");//h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);//h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);h_line[j]->Draw("");}else {h_line[j]->SetLineColor(j+1);h_line[j]->Draw("same");}}//! 2d scanTCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);canvas61->cd();TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);for(int i=parameter1; i<=parameter2; i++) {//int canvasPosition = nX*(i-chYstart)+(j-parameter1)+1;//int canvasPosition = nX*(chYend-i)+parameter1+1;//if (debug) printf("canvas %d\n",canvasPosition);int chPosition = map[i][chYstart];char hname[128];sprintf(hname, "h2d%d", chPosition);int histogram = i;h[histogram] = (TH2F *) rootfile->Get(hname);int noise = getNoise(h[histogram], 1, 100);for(int k=minX; k<=maxX; k++) {for(int l=minY; l<=maxY; l++) {int signal = h[histogram]->GetBinContent(k,l); // detected//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 signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;double p = signal - noise;p /= 10000.0;double eta = -log(1 - p);//double x = k*(xUpUser-xLowUser)/double(binsX);//double y = l*(yUpUser-yLowUser)/double(binsY);double x = h[histogram]->GetXaxis()->GetBinCenter(k);double y = h[histogram]->GetYaxis()->GetBinCenter(l);h_corrected->Fill(x, y, eta);}}}h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");h_corrected->GetZaxis()->SetRangeUser(-0.05,0.30);h_corrected->Draw("colz");}/** Draws the sum of the channels* Each channel is a 2d plot* Intended for the study of 1 channel** Input: scanned channel coordinates (x,y)*/if (strstr(plopt, "2d") != NULL) {int nX = 3;int nY = 3;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=parameter2-1; i<=parameter2+1; i++) {for(int j=parameter1-1; j<=parameter1+1; j++) {//int canvasPosition = nX*(i-chYstart)+(j-parameter1)+1;int canvasPosition = nX*(parameter2+1-i) + (j-parameter1+1) +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");if ( !((i == parameter2) and (j == parameter1)) ) h_2d->GetZaxis()->SetRangeUser(0,200);} //x}// Number of photoelectrons - Poissonian correctionTCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1600,1600);canvas8->cd();gStyle->SetOptStat(0);char hname[128];int chPosition = map[parameter1][parameter2];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()->GetXmin();Double_t xUpUser = h_2d->GetXaxis()->GetXmax();Double_t yLowUser = h_2d->GetYaxis()->GetXmin();Double_t yUpUser = h_2d->GetYaxis()->GetXmax();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 diffX = xUpUser - xLowUser;double diffY = yUpUser - yLowUser;//TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, 3.00, 3.00+diffX, binsY, 8.00, 8.00+diffY);//TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, minX, maxX, binsY, minY, maxY);double noise = getNoise(h_2d, 1, 70);noise = 61;printf("[2d scan] Average 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); // detecteddouble p = signal - noise;p /= 10000.;double p0 = 1.0 - p; // events with zero photonsdouble 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 = h_2d->GetXaxis()->GetBinCenter(k);double y = h_2d->GetYaxis()->GetBinCenter(j);if (debug) printf("x=%f y=%f\n", x, y);h_corrected->Fill(x, y, eta);}}//h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");h_corrected->SetTitle(";x[mm];y[mm]");//gStyle->SetPalette(52,0); // black and white for print//gStyle->SetPalette(1); //black and white 2nd option//SetGS(); // inverse//h_corrected->GetZaxis()->SetRangeUser(-0.01,0.12);h_corrected->SetContour(50);h_corrected->Draw("colz");// collection efficiencyint nPoints =0;double efficiency=0;for (int i=minX; i<=maxX; i++) {for (int j=minY; j<=maxY; j++) {double signal = h_corrected->GetBinContent(i,j);if(debug) printf("signal %f\n",signal);if (signal > 0.02) efficiency += signal;nPoints++;}}printf("Signal sum = %f\n # of points = %d\n",efficiency,nPoints);if(strstr(plopt, "s") != NULL) {char photo[128];sprintf(photo, "ps/%s.png", filename);//canvas8->Print("result.pdf");canvas8->SaveAs(photo);}}/** 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 = parameter2 - parameter1 + 1;int nY = chYend - chYstart + 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[nX*nY];// 2d histogram noise subtraction and poisson scalingfor(int i=chYstart; i<=chYend; i++) {for(int j=parameter1; j<=parameter2; j++) {int chPosition = map[j][i];char hname[128];sprintf(hname, "h2d%d", chPosition);int histogram = nX*(i-chYstart)+(j-parameter1);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, y, eta);}}}}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");}/** Draws the beam profile and fits it with error function* on some background function*/if (strstr(plopt, "beam") != NULL) {TCanvas* canvas9 = new TCanvas("canvas9","canvas9", 500,500);canvas9->cd();char hname[128];sprintf(hname, "hnhitsx%d", 49);TH1F* h_laser = (TH1F*)rootfile->Get(hname);h_laser->Draw();h_laser->SetStats(1);TF1* err = new TF1("err","[0]+[1]*TMath::Erf((x-[2])/[3])",parameter1,parameter2);err->SetParameter(0,2500);err->SetParameter(1, h_laser->GetMaximum());err->SetParameter(2, h_laser->GetBinCenter(h_laser->GetMaximumBin()));err->SetParameter(3, 0.001);h_laser->Fit(err,"qr");h_laser->Fit(err,"lr");double sigma = err->GetParameter(3);printf("sigma = %2.0f um, FWHM = %2.0f um\n", sigma*1000, 2.35*sigma*1000);}if (strstr(plopt, "map") != NULL) {TCanvas* canvas1 = new TCanvas("canvas1","canvas1",1000,1000);canvas1->cd();TH2I* h_map = new TH2I("h_map","h_map",8,-0.5,7.5,8,-0.5,7.5);for (int i=7; i>=0; i--) {//for (int j=7; j>=0; j--) printf("(%d, %d) ", j,i);for (int j=7; j>=0; j--) {printf("%2d (%d %d)\n", map[j][i], j*10080+5040, i*10080+5040);h_map->Fill(j,i, map[j][i]);}printf("\n");}gStyle->SetOptStat(0);h_map->Draw("text");}if(strstr(plopt, "p") != NULL) {TCanvas* canvasLast = new TCanvas("canvasLast","",500,500);canvasLast->Print("result.pdf)");}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,1);//if (noise < value) noise = value;noise += value;count++;}return (noise/double(count));}