WebSVN – f9daq – Diff – /sipmScan/examples/plotThresholdScan.C


#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "stdio.h"
 
#include "include/RTUtil.h"
 
double getNoise(TH2F*, int, int);
 
int plotThreshold(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)
{
  const int c_nChannels = 64;
  const double c_xOffset = 1; // mm
  const double c_yOffset = 0.7;
 
  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 style
        gStyle->SetOptStat("ne");
        gStyle->SetPalette(1, 0);
       
        gStyle->SetPaperSize(TStyle::kA4);
        gStyle->SetStatBorderSize(1);
        gStyle->SetFrameBorderMode(0);
        gStyle->SetFrameFillColor(0);
        gStyle->SetCanvasBorderMode(0);
        gStyle->SetPadBorderMode(0);
        gStyle->SetPadColor(0);
        gStyle->SetCanvasColor(0);
        gStyle->SetStatColor(0);
        gStyle->SetOptFit(11);
        gStyle->SetOptStat();
        gStyle->SetPadRightMargin(0.15);
        gStyle->SetPadLeftMargin(0.12);
        //gStyle->SetTitleYOffset(1.4);
 
  if (strstr(plopt, "th") != NULL) {
   TCanvas* canvas1 = new TCanvas("canvas1","",500,500);
   canvas1->cd();
   TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
   h_threshold->Draw("colz");
  }
 
  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();
  }
 
  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();
  }
 
  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=chXstart; i<=chXend; 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 == chXstart)
      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[chXstart][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);
  }
 
  /** 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. hits
    TH2F* 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=chXstart; j<=chXend; 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 = k*(xUpUser-xLowUser)/double(binsX);
          //double y = l*(yUpUser-yLowUser)/double(binsY);
          h_line[j]->Fill(x-c_xOffset, signal);
        //}
      }
      if (j == chXstart) {
        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 scan
    TCanvas* 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=chXstart; i<=chXend; i++) {
      //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
      //int canvasPosition = nX*(chYend-i)+chXstart+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);
            h_corrected->Fill(x-c_xOffset, y-c_yOffset, 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
    */
  if (strstr(plopt, "2d") != NULL) {
   
    int nX = chXend - chXstart + 1;
    int nY = chYend - chYstart + 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=chYstart; i<=chYend; i++) {
      for(int j=chXstart; j<=chXend; j++) {
      //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
      int canvasPosition = nX*(chYend-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 = chXend - chXstart + 1;
    int nY = chYend - chYstart + 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=chYstart; i<=chYend; 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-chYstart)+(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");
  }
 
  /** 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", 36);
  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])",17.2,17.50);
  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) {
    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 ", map[j][i]);
      printf("\n");
    }
     
  }
   
  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));
}
 
 

Rev 30	Rev 34
1	#include "TROOT.h"	1	#include "TROOT.h"
2	#include "TFile.h"	2	#include "TFile.h"
3	#include "TBenchmark.h"	3	#include "TBenchmark.h"
4	#include "TH1F.h"	4	#include "TH1F.h"
5	#include "TH2F.h"	5	#include "TH2F.h"
6	#include "TCanvas.h"	6	#include "TCanvas.h"
7	#include "TStyle.h"	7	#include "TStyle.h"
8	#include "TPad.h"	8	#include "TPad.h"
9	#include "TF1.h"	9	#include "TF1.h"
10	#include "TGraph.h"	10	#include "TGraph.h"
11	#include "TSpectrum.h"	11	#include "TSpectrum.h"
12	#include "stdio.h"	12	#include "stdio.h"
13		13
14	#include "RTUtil.h"	14	#include "include/RTUtil.h"
15		15
16	double getNoise(TH2F*, int, int);	16	double getNoise(TH2F*, int, int);
17		17
18	int plotThreshold(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)	18	int plotThreshold(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)
19	{	19	{
20	const int c_nChannels = 64;	20	const int c_nChannels = 64;
21	const double c_xOffset = 1; // mm	21	const double c_xOffset = 1; // mm
22	const double c_yOffset = 0.7;	22	const double c_yOffset = 0.7;
23		23
24	int map[8][8]={{32,34,53,55,40,42,61,63},	24	int map[8][8]={{32,34,53,55,40,42,61,63},
25	{48,50,37,39,56,58,45,47},	25	{48,50,37,39,56,58,45,47},
26	{33,35,52,54,41,43,60,62},	26	{33,35,52,54,41,43,60,62},
27	{49,51,36,38,57,59,44,46},	27	{49,51,36,38,57,59,44,46},
28	{17,19,4,6,25,27,12,14},	28	{17,19,4,6,25,27,12,14},
29	{1,3,20,22,9,11,28,30},	29	{1,3,20,22,9,11,28,30},
30	{16,18,5,7,24,26,13,15},	30	{16,18,5,7,24,26,13,15},
31	{0,2,21,23,8,10,29,31}	31	{0,2,21,23,8,10,29,31}
32	};	32	};
33		33
34	char fnameroot[256];	34	char fnameroot[256];
35	TFile* rootfile;	35	TFile* rootfile;
36	sprintf(fnameroot, "root/%s.root", filename);	36	sprintf(fnameroot, "root/%s.root", filename);
37	rootfile = (TFile *) gROOT->FindObject(filename);	37	rootfile = (TFile *) gROOT->FindObject(filename);
38	if(rootfile==NULL) rootfile = new TFile(fnameroot);	38	if(rootfile==NULL) rootfile = new TFile(fnameroot);
39	if(rootfile==NULL) {	39	if(rootfile==NULL) {
40	printf("Cannot open root file %s!!!\n",fnameroot);	40	printf("Cannot open root file %s!!!\n",fnameroot);
41	return(0);	41	return(0);
42	}	42	}
43		43
44	// set draw style	44	// set draw style
45	gStyle->SetOptStat("ne");	45	gStyle->SetOptStat("ne");
46	gStyle->SetPalette(1, 0);	46	gStyle->SetPalette(1, 0);
47		47
48	gStyle->SetPaperSize(TStyle::kA4);	48	gStyle->SetPaperSize(TStyle::kA4);
49	gStyle->SetStatBorderSize(1);	49	gStyle->SetStatBorderSize(1);
50	gStyle->SetFrameBorderMode(0);	50	gStyle->SetFrameBorderMode(0);
51	gStyle->SetFrameFillColor(0);	51	gStyle->SetFrameFillColor(0);
52	gStyle->SetCanvasBorderMode(0);	52	gStyle->SetCanvasBorderMode(0);
53	gStyle->SetPadBorderMode(0);	53	gStyle->SetPadBorderMode(0);
54	gStyle->SetPadColor(0);	54	gStyle->SetPadColor(0);
55	gStyle->SetCanvasColor(0);	55	gStyle->SetCanvasColor(0);
56	gStyle->SetStatColor(0);	56	gStyle->SetStatColor(0);
57	gStyle->SetOptFit(11);	57	gStyle->SetOptFit(11);
58	gStyle->SetOptStat();	58	gStyle->SetOptStat();
59	gStyle->SetPadRightMargin(0.15);	59	gStyle->SetPadRightMargin(0.15);
60	gStyle->SetPadLeftMargin(0.12);	60	gStyle->SetPadLeftMargin(0.12);
61	//gStyle->SetTitleYOffset(1.4);	61	//gStyle->SetTitleYOffset(1.4);
62		62
63	if (strstr(plopt, "th") != NULL) {	63	if (strstr(plopt, "th") != NULL) {
64	TCanvas* canvas1 = new TCanvas("canvas1","",500,500);	64	TCanvas* canvas1 = new TCanvas("canvas1","",500,500);
65	canvas1->cd();	65	canvas1->cd();
66	TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");	66	TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
67	h_threshold->Draw("colz");	67	h_threshold->Draw("colz");
68	}	68	}
69		69
70	if( strstr(plopt, "all") != NULL ) {	70	if( strstr(plopt, "all") != NULL ) {
71	TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);	71	TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);
72	TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);	72	TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);
73	canvas2->Divide(8,8);	73	canvas2->Divide(8,8);
74	canvas3->Divide(8,8);	74	canvas3->Divide(8,8);
75	TH1F* h_hitsx;	75	TH1F* h_hitsx;
76	TH1F* h_hitsy;	76	TH1F* h_hitsy;
77	for(int i=0; i<c_nChannels; i++) {	77	for(int i=0; i<c_nChannels; i++) {
78	canvas2->cd(i+1);	78	canvas2->cd(i+1);
79	char hname[128];	79	char hname[128];
80	sprintf(hname, "hnhitsx%d", i);	80	sprintf(hname, "hnhitsx%d", i);
81	h_hitsx = (TH1F*)rootfile->Get(hname);	81	h_hitsx = (TH1F*)rootfile->Get(hname);
82	h_hitsx->Draw();	82	h_hitsx->Draw();
83	canvas3->cd(i+1);	83	canvas3->cd(i+1);
84	sprintf(hname, "hnhitsy%d", i);	84	sprintf(hname, "hnhitsy%d", i);
85	h_hitsy = (TH1F*)rootfile->Get(hname);	85	h_hitsy = (TH1F*)rootfile->Get(hname);
86	h_hitsy->Draw();	86	h_hitsy->Draw();
87	}	87	}
88	}	88	}
89		89
90	if( strstr(plopt, "x") != NULL ) {	90	if( strstr(plopt, "x") != NULL ) {
91	TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);	91	TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);
92	TH1F* h_hitsx;	92	TH1F* h_hitsx;
93	canvas10->cd();	93	canvas10->cd();
94	char hname[128];	94	char hname[128];
95	sprintf(hname, "hnhitsx%d", chXstart);	95	sprintf(hname, "hnhitsx%d", chXstart);
96	h_hitsx = (TH1F*)rootfile->Get(hname);	96	h_hitsx = (TH1F*)rootfile->Get(hname);
97	h_hitsx->Draw();	97	h_hitsx->Draw();
98	}	98	}
99		99
100	if( strstr(plopt, "y") != NULL ) {	100	if( strstr(plopt, "y") != NULL ) {
101	TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);	101	TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);
102	TH1F* h_hitsy;	102	TH1F* h_hitsy;
103	canvas11->cd();	103	canvas11->cd();
104	char hname[128];	104	char hname[128];
105	sprintf(hname, "hnhitsy%d", chXstart);	105	sprintf(hname, "hnhitsy%d", chXstart);
106	h_hitsy = (TH1F*)rootfile->Get(hname);	106	h_hitsy = (TH1F*)rootfile->Get(hname);
107	h_hitsy->Draw();	107	h_hitsy->Draw();
108	}	108	}
109		109
110	if( strstr(plopt, "share") != NULL ) {	110	if( strstr(plopt, "share") != NULL ) {
111	/TCanvas canvas4 = new TCanvas("canvas1","canvas1",1000,1000);	111	/TCanvas canvas4 = new TCanvas("canvas1","canvas1",1000,1000);
112	int nChannels = chYend-chYstart+1;	112	int nChannels = chYend-chYstart+1;
113	int ncols = nChannels/2;	113	int ncols = nChannels/2;
114	printf("nch %d nch\\2 %d\n", nChannels, ncols);	114	printf("nch %d nch\\2 %d\n", nChannels, ncols);
115	canvas4->Divide(2,ncols);	115	canvas4->Divide(2,ncols);
116	TH1F* h_hitsy;	116	TH1F* h_hitsy;
117	for(int i=chYstart; i<=chYend; i++){	117	for(int i=chYstart; i<=chYend; i++){
118	canvas4->cd(i-chYstart+1);	118	canvas4->cd(i-chYstart+1);
119	char hname[128];	119	char hname[128];
120	int chPosition = map[0][i];	120	int chPosition = map[0][i];
121	sprintf(hname, "hnhitsy%d", chPosition);	121	sprintf(hname, "hnhitsy%d", chPosition);
122	h_hitsy = (TH1F*)rootfile->Get(hname);	122	h_hitsy = (TH1F*)rootfile->Get(hname);
123	h_hitsy->Draw();	123	h_hitsy->Draw();
124	}*/	124	}*/
125		125
126	TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);	126	TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);
127	canvas4->cd();	127	canvas4->cd();
128	for(int i=chXstart; i<=chXend; i++) {	128	for(int i=chXstart; i<=chXend; i++) {
129	TH1F* h_hitsx;	129	TH1F* h_hitsx;
130	char hname[128];	130	char hname[128];
131	int chPosition = map[i][chYstart];	131	int chPosition = map[i][chYstart];
132	sprintf(hname, "hnhitsx%d", chPosition);	132	sprintf(hname, "hnhitsx%d", chPosition);
133	h_hitsx = (TH1F*)rootfile->Get(hname);	133	h_hitsx = (TH1F*)rootfile->Get(hname);
134	h_hitsx->SetTitle("Scan X;x [mm]; Entries");	134	h_hitsx->SetTitle("Scan X;x [mm]; Entries");
135	h_hitsx->GetYaxis()->SetTitleOffset(1.3);	135	h_hitsx->GetYaxis()->SetTitleOffset(1.3);
136	h_hitsx->SetStats(0);	136	h_hitsx->SetStats(0);
137	if (i == chXstart)	137	if (i == chXstart)
138	h_hitsx->Draw();	138	h_hitsx->Draw();
139	else {	139	else {
140	h_hitsx->SetLineColor(i+1);	140	h_hitsx->SetLineColor(i+1);
141	h_hitsx->Draw("same");	141	h_hitsx->Draw("same");
142	}	142	}
143	}	143	}
144	//sprintf(fullname, "ps/%s_Yshare.eps", filename);	144	//sprintf(fullname, "ps/%s_Yshare.eps", filename);
145	//canvas4->SaveAs(fullname);	145	//canvas4->SaveAs(fullname);
146		146
147	TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);	147	TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);
148	canvas5->cd();	148	canvas5->cd();
149	for(int i=chYstart; i<=chYend; i++) {	149	for(int i=chYstart; i<=chYend; i++) {
150	TH1F* h_hitsy;	150	TH1F* h_hitsy;
151	char hname[128];	151	char hname[128];
152	int chPosition = map[chXstart][i];	152	int chPosition = map[chXstart][i];
153	sprintf(hname, "hnhitsy%d", chPosition);	153	sprintf(hname, "hnhitsy%d", chPosition);
154	h_hitsy = (TH1F*)rootfile->Get(hname);	154	h_hitsy = (TH1F*)rootfile->Get(hname);
155	h_hitsy->SetTitle("Scan Y;y [mm]; Entries");	155	h_hitsy->SetTitle("Scan Y;y [mm]; Entries");
156	h_hitsy->GetYaxis()->SetTitleOffset(1.3);	156	h_hitsy->GetYaxis()->SetTitleOffset(1.3);
157	h_hitsy->SetStats(0);	157	h_hitsy->SetStats(0);
158	if (i == chYstart)	158	if (i == chYstart)
159	h_hitsy->Draw();	159	h_hitsy->Draw();
160	else {	160	else {
161	h_hitsy->SetLineColor(i+1);	161	h_hitsy->SetLineColor(i+1);
162	h_hitsy->Draw("same");	162	h_hitsy->Draw("same");
163	}	163	}
164	}	164	}
165	//sprintf(fullname, "ps/%s_Yshare.eps", filename);	165	//sprintf(fullname, "ps/%s_Yshare.eps", filename);
166	//canvas5->SaveAs(fullname);	166	//canvas5->SaveAs(fullname);
167	}	167	}
168		168
169	/** Draws the signal from 8 channels in x-row	169	/** Draws the signal from 8 channels in x-row
170	* for one specific y bin, so the background and cross-talk	170	* for one specific y bin, so the background and cross-talk
171	* can be estimated.	171	* can be estimated.
172	* Draws also a 2d scan of these channels.	172	* Draws also a 2d scan of these channels.
173	*/	173	*/
174	if (strstr(plopt, "line") != NULL) {	174	if (strstr(plopt, "line") != NULL) {
175	TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);	175	TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);
176	canvas6->cd(0);	176	canvas6->cd(0);
177	gStyle->SetOptStat(0);	177	gStyle->SetOptStat(0);
178		178
179	TH2F* h0 = (TH2F*) rootfile->Get("h2d0");	179	TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
180	Int_t binsX = h0->GetXaxis()->GetNbins();	180	Int_t binsX = h0->GetXaxis()->GetNbins();
181	Int_t minX = h0->GetXaxis()->GetFirst();	181	Int_t minX = h0->GetXaxis()->GetFirst();
182	Int_t maxX = h0->GetXaxis()->GetLast()+1;	182	Int_t maxX = h0->GetXaxis()->GetLast()+1;
183	Int_t binsY = h0->GetYaxis()->GetNbins();	183	Int_t binsY = h0->GetYaxis()->GetNbins();
184	Int_t minY = h0->GetYaxis()->GetFirst();	184	Int_t minY = h0->GetYaxis()->GetFirst();
185	Int_t maxY = h0->GetYaxis()->GetLast()+1;	185	Int_t maxY = h0->GetYaxis()->GetLast()+1;
186	Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);	186	Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
187	Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);	187	Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
188	Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);	188	Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
189	Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);	189	Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
190	if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);	190	if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
191		191
192	//! 1-dimension position in x vs. hits	192	//! 1-dimension position in x vs. hits
193	TH2F* h[8];	193	TH2F* h[8];
194	TH1F* h_line[8];	194	TH1F* h_line[8];
195	for(int j=0; j<8; j++) {	195	for(int j=0; j<8; j++) {
196	h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);	196	h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);
197	}	197	}
198		198
199	for(int j=chXstart; j<=chXend; j++) {	199	for(int j=chXstart; j<=chXend; j++) {
200	int chPosition = map[j][chYstart];	200	int chPosition = map[j][chYstart];
201	char hname[128];	201	char hname[128];
202	sprintf(hname, "h2d%d", chPosition);	202	sprintf(hname, "h2d%d", chPosition);
203	int histogram = j;	203	int histogram = j;
204	h[histogram] = (TH2F *) rootfile->Get(hname);	204	h[histogram] = (TH2F *) rootfile->Get(hname);
205	int noise = getNoise(h[histogram], 1, 160);	205	int noise = getNoise(h[histogram], 1, 160);
206	for(int k=minX; k<=maxX; k++) {	206	for(int k=minX; k<=maxX; k++) {
207	int l=chYstart*20+12;	207	int l=chYstart*20+12;
208	//for(int l=12; l<=16; l++) {	208	//for(int l=12; l<=16; l++) {
209	double signal = h[histogram]->GetBinContent(k,l);	209	double signal = h[histogram]->GetBinContent(k,l);
210	//signal -= noise;	210	//signal -= noise;
211	//signal /= 5*10000.0;	211	//signal /= 5*10000.0;
212	double eta = -log(1 - signal);	212	double eta = -log(1 - signal);
213	double x = k*(xUpUser-xLowUser)/double(binsX);	213	double x = k*(xUpUser-xLowUser)/double(binsX);
214	//double y = l*(yUpUser-yLowUser)/double(binsY);	214	//double y = l*(yUpUser-yLowUser)/double(binsY);
215	h_line[j]->Fill(x-c_xOffset, signal);	215	h_line[j]->Fill(x-c_xOffset, signal);
216	//}	216	//}
217	}	217	}
218	if (j == chXstart) {	218	if (j == chXstart) {
219	h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");	219	h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");
220	//h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);	220	//h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);
221	//h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);	221	//h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);
222	h_line[j]->Draw("");	222	h_line[j]->Draw("");
223	}	223	}
224	else {	224	else {
225	h_line[j]->SetLineColor(j+1);	225	h_line[j]->SetLineColor(j+1);
226	h_line[j]->Draw("same");	226	h_line[j]->Draw("same");
227	}	227	}
228	}	228	}
229		229
230		230
231	//! 2d scan	231	//! 2d scan
232	TCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);	232	TCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);
233	canvas61->cd();	233	canvas61->cd();
234	TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);	234	TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
235	for(int i=chXstart; i<=chXend; i++) {	235	for(int i=chXstart; i<=chXend; i++) {
236	//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;	236	//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
237	//int canvasPosition = nX*(chYend-i)+chXstart+1;	237	//int canvasPosition = nX*(chYend-i)+chXstart+1;
238	//if (debug) printf("canvas %d\n",canvasPosition);	238	//if (debug) printf("canvas %d\n",canvasPosition);
239	int chPosition = map[i][chYstart];	239	int chPosition = map[i][chYstart];
240	char hname[128];	240	char hname[128];
241	sprintf(hname, "h2d%d", chPosition);	241	sprintf(hname, "h2d%d", chPosition);
242	int histogram = i;	242	int histogram = i;
243	h[histogram] = (TH2F *) rootfile->Get(hname);	243	h[histogram] = (TH2F *) rootfile->Get(hname);
244	int noise = getNoise(h[histogram], 1, 100);	244	int noise = getNoise(h[histogram], 1, 100);
245		245
246	for(int k=minX; k<=maxX; k++) {	246	for(int k=minX; k<=maxX; k++) {
247	for(int l=minY; l<=maxY; l++) {	247	for(int l=minY; l<=maxY; l++) {
248	int signal = h[histogram]->GetBinContent(k,l); // detected	248	int signal = h[histogram]->GetBinContent(k,l); // detected
249	//p /= 10000.;	249	//p /= 10000.;
250	//double p0 = 1.0 - p; // events with zero photons	250	//double p0 = 1.0 - p; // events with zero photons
251	//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);	251	//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
252	//double eta = -log(p0);	252	//double eta = -log(p0);
253	//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);	253	//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
254	//double signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;	254	//double signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;
255	double p = signal - noise;	255	double p = signal - noise;
256	p /= 10000.0;	256	p /= 10000.0;
257	double eta = -log(1 - p);	257	double eta = -log(1 - p);
258	double x = k*(xUpUser-xLowUser)/double(binsX);	258	double x = k*(xUpUser-xLowUser)/double(binsX);
259	double y = l*(yUpUser-yLowUser)/double(binsY);	259	double y = l*(yUpUser-yLowUser)/double(binsY);
260	h_corrected->Fill(x-c_xOffset, y-c_yOffset, eta);	260	h_corrected->Fill(x-c_xOffset, y-c_yOffset, eta);
261	}	261	}
262	}	262	}
263	}	263	}
264	h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");	264	h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");
265	h_corrected->GetZaxis()->SetRangeUser(-0.05,0.30);	265	h_corrected->GetZaxis()->SetRangeUser(-0.05,0.30);
266	h_corrected->Draw("colz");	266	h_corrected->Draw("colz");
267		267
268	}	268	}
269		269
270	/** Draws the sum of the channels	270	/** Draws the sum of the channels
271	* Each channel is a 2d plot	271	* Each channel is a 2d plot
272	* Intended for the study of 1 channel	272	* Intended for the study of 1 channel
273	*/	273	*/
274	if (strstr(plopt, "2d") != NULL) {	274	if (strstr(plopt, "2d") != NULL) {
275		275
276	int nX = chXend - chXstart + 1;	276	int nX = chXend - chXstart + 1;
277	int nY = chYend - chYstart + 1;	277	int nY = chYend - chYstart + 1;
278	TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX400,nY400);	278	TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX400,nY400);
279	printf("nx %d ny %d\n",nX,nY);	279	printf("nx %d ny %d\n",nX,nY);
280	canvas7->Divide(nX,nY);	280	canvas7->Divide(nX,nY);
281	for(int i=chYstart; i<=chYend; i++) {	281	for(int i=chYstart; i<=chYend; i++) {
282	for(int j=chXstart; j<=chXend; j++) {	282	for(int j=chXstart; j<=chXend; j++) {
283	//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;	283	//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
284	int canvasPosition = nX*(chYend-i)+(j-chXstart)+1;	284	int canvasPosition = nX*(chYend-i)+(j-chXstart)+1;
285	if (debug) printf("canvas %d\n",canvasPosition);	285	if (debug) printf("canvas %d\n",canvasPosition);
286	canvas7->cd(canvasPosition);	286	canvas7->cd(canvasPosition);
287	char hname[128];	287	char hname[128];
288	int chPosition = map[j][i];	288	int chPosition = map[j][i];
289	sprintf(hname, "h2d%d", chPosition);	289	sprintf(hname, "h2d%d", chPosition);
290	TH2F* h_2d = (TH2F*)rootfile->Get(hname);	290	TH2F* h_2d = (TH2F*)rootfile->Get(hname);
291	h_2d->Draw("colz");	291	h_2d->Draw("colz");
292	} //x	292	} //x
293	}	293	}
294		294
295	// Number of photoelectrons - Poissonian correction	295	// Number of photoelectrons - Poissonian correction
296	TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);	296	TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);
297	canvas8->cd();	297	canvas8->cd();
298	gStyle->SetOptStat(0);	298	gStyle->SetOptStat(0);
299	char hname[128];	299	char hname[128];
300	int chPosition = map[1][2];	300	int chPosition = map[1][2];
301	sprintf(hname, "h2d%d", chPosition);	301	sprintf(hname, "h2d%d", chPosition);
302	TH2F* h_2d = (TH2F*)rootfile->Get(hname);	302	TH2F* h_2d = (TH2F*)rootfile->Get(hname);
303		303
304	Int_t binsX = h_2d->GetXaxis()->GetNbins();	304	Int_t binsX = h_2d->GetXaxis()->GetNbins();
305	Int_t minX = h_2d->GetXaxis()->GetFirst();	305	Int_t minX = h_2d->GetXaxis()->GetFirst();
306	Int_t maxX = h_2d->GetXaxis()->GetLast()+1;	306	Int_t maxX = h_2d->GetXaxis()->GetLast()+1;
307	Int_t binsY = h_2d->GetYaxis()->GetNbins();	307	Int_t binsY = h_2d->GetYaxis()->GetNbins();
308	Int_t minY = h_2d->GetYaxis()->GetFirst();	308	Int_t minY = h_2d->GetYaxis()->GetFirst();
309	Int_t maxY = h_2d->GetYaxis()->GetLast()+1;	309	Int_t maxY = h_2d->GetYaxis()->GetLast()+1;
310	Double_t xLowUser = h_2d->GetXaxis()->GetBinLowEdge(minX);	310	Double_t xLowUser = h_2d->GetXaxis()->GetBinLowEdge(minX);
311	Double_t xUpUser = h_2d->GetXaxis()->GetBinUpEdge(maxX);	311	Double_t xUpUser = h_2d->GetXaxis()->GetBinUpEdge(maxX);
312	Double_t yLowUser = h_2d->GetYaxis()->GetBinLowEdge(minY);	312	Double_t yLowUser = h_2d->GetYaxis()->GetBinLowEdge(minY);
313	Double_t yUpUser = h_2d->GetYaxis()->GetBinUpEdge(maxY);	313	Double_t yUpUser = h_2d->GetYaxis()->GetBinUpEdge(maxY);
314	if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);	314	if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
315	TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);	315	TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);
316		316
317	double noise = getNoise(h_2d, 1, 89);	317	double noise = getNoise(h_2d, 1, 89);
318	if(debug) printf("Noise = %f\n", noise);	318	if(debug) printf("Noise = %f\n", noise);
319	for(int k=minX; k<=maxX; k++) {	319	for(int k=minX; k<=maxX; k++) {
320	for(int j=minY; j<=maxY; j++) {	320	for(int j=minY; j<=maxY; j++) {
321	double signal = h_2d->GetBinContent(k,j); // detected	321	double signal = h_2d->GetBinContent(k,j); // detected
322	//double p = ((signal - noise) > 1) ? (signal-noise) : 1;	322	//double p = ((signal - noise) > 1) ? (signal-noise) : 1;
323	double p = signal - noise;	323	double p = signal - noise;
324	p /= 10000.;	324	p /= 10000.;
325	double p0 = 1.0 - p; // events with zero photons	325	double p0 = 1.0 - p; // events with zero photons
326	//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);	326	//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
327	double eta = -log(p0); // constant of the poissonian statistics	327	double eta = -log(p0); // constant of the poissonian statistics
328	if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);	328	if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
329	double x = k*(xUpUser-xLowUser)/double(binsX);	329	double x = k*(xUpUser-xLowUser)/double(binsX);
330	double y = j*(yUpUser-yLowUser)/double(binsY);	330	double y = j*(yUpUser-yLowUser)/double(binsY);
331	h_corrected->Fill(x+3,y+8, eta);	331	h_corrected->Fill(x+3,y+8, eta);
332	}	332	}
333	}	333	}
334	h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");	334	h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");
335	h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);	335	h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);
336	h_corrected->Draw("colz");	336	h_corrected->Draw("colz");
337		337
338	// collection efficiency	338	// collection efficiency
339	int nPoints =0;	339	int nPoints =0;
340	double efficiency=0;	340	double efficiency=0;
341	for (int i=18; i<=58; i++) {	341	for (int i=18; i<=58; i++) {
342	for (int j=19; j<=59; j++) {	342	for (int j=19; j<=59; j++) {
343	double signal = h_corrected->GetBinContent(i,j);	343	double signal = h_corrected->GetBinContent(i,j);
344	if(debug) printf("signal %f\n",signal);	344	if(debug) printf("signal %f\n",signal);
345	efficiency += signal;	345	efficiency += signal;
346	nPoints++;	346	nPoints++;
347	}	347	}
348	}	348	}
349	printf("Signal sum = %f\n # of points = %d\n",efficiency,nPoints);	349	printf("Signal sum = %f\n # of points = %d\n",efficiency,nPoints);
350	}	350	}
351		351
352	/** Draws the sum of channel signals	352	/** Draws the sum of channel signals
353	* Each channel is a 2d ('h2d') histogram	353	* Each channel is a 2d ('h2d') histogram
354	* Suitable for 8x8 chs scan	354	* Suitable for 8x8 chs scan
355	*/	355	*/
356	if( strstr(plopt, "sum") != NULL ) {	356	if( strstr(plopt, "sum") != NULL ) {
357	int nX = chXend - chXstart + 1;	357	int nX = chXend - chXstart + 1;
358	int nY = chYend - chYstart + 1;	358	int nY = chYend - chYstart + 1;
359	TCanvas* canvas12 = new TCanvas("canvas12","c2",8200, 8200);	359	TCanvas* canvas12 = new TCanvas("canvas12","c2",8200, 8200);
360	canvas12->cd();	360	canvas12->cd();
361	gStyle->SetOptStat(0);	361	gStyle->SetOptStat(0);
362		362
363	// final histogram parameters	363	// final histogram parameters
364	TH2F* h0 = (TH2F*) rootfile->Get("h2d0");	364	TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
365	Int_t binsX = h0->GetXaxis()->GetNbins();	365	Int_t binsX = h0->GetXaxis()->GetNbins();
366	Int_t minX = h0->GetXaxis()->GetFirst();	366	Int_t minX = h0->GetXaxis()->GetFirst();
367	Int_t maxX = h0->GetXaxis()->GetLast()+1;	367	Int_t maxX = h0->GetXaxis()->GetLast()+1;
368	Int_t binsY = h0->GetYaxis()->GetNbins();	368	Int_t binsY = h0->GetYaxis()->GetNbins();
369	Int_t minY = h0->GetYaxis()->GetFirst();	369	Int_t minY = h0->GetYaxis()->GetFirst();
370	Int_t maxY = h0->GetYaxis()->GetLast()+1;	370	Int_t maxY = h0->GetYaxis()->GetLast()+1;
371	Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);	371	Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
372	Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);	372	Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
373	Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);	373	Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
374	Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);	374	Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
375	if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);	375	if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
376	TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);	376	TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
377	TH2F* h[9];	377	TH2F* h[9];
378		378
379	// 2d histogram noise subtraction and poisson scaling	379	// 2d histogram noise subtraction and poisson scaling
380	for(int i=chYstart; i<=chYend; i++) {	380	for(int i=chYstart; i<=chYend; i++) {
381	for(int j=chXstart; j<=chXend; j++) {	381	for(int j=chXstart; j<=chXend; j++) {
382	int chPosition = map[j][i];	382	int chPosition = map[j][i];
383	char hname[128];	383	char hname[128];
384	sprintf(hname, "h2d%d", chPosition);	384	sprintf(hname, "h2d%d", chPosition);
385	int histogram = nX*(i-chYstart)+(j-chXstart);	385	int histogram = nX*(i-chYstart)+(j-chXstart);
386	h[histogram] = (TH2F *) rootfile->Get(hname);	386	h[histogram] = (TH2F *) rootfile->Get(hname);
387	int noise = getNoise(h[histogram], 1, 170);	387	int noise = getNoise(h[histogram], 1, 170);
388	if (debug) printf("noise: %d\n",noise);	388	if (debug) printf("noise: %d\n",noise);
389	for(int k=minX; k<=maxX; k++) {	389	for(int k=minX; k<=maxX; k++) {
390	for(int l=minY; l<=maxY; l++) {	390	for(int l=minY; l<=maxY; l++) {
391	int signal = h[histogram]->GetBinContent(k,l); // detected	391	int signal = h[histogram]->GetBinContent(k,l); // detected
392	//double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;	392	//double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;
393	double p = signal - noise;	393	double p = signal - noise;
394	p /= 10000.;	394	p /= 10000.;
395	double p0 = 1.0 - p; // events with zero photons	395	double p0 = 1.0 - p; // events with zero photons
396	//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);	396	//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
397	double eta = -log(p0);	397	double eta = -log(p0);
398	//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);	398	//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
399	double x = k*(xUpUser-xLowUser)/double(binsX);	399	double x = k*(xUpUser-xLowUser)/double(binsX);
400	double y = l*(yUpUser-yLowUser)/double(binsY);	400	double y = l*(yUpUser-yLowUser)/double(binsY);
401	h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);	401	h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);
402	}	402	}
403	}	403	}
404		404
405	}	405	}
406	}	406	}
407		407
408	h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");	408	h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");
409	//h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);	409	//h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);
410	h_corrected->Draw("colz");	410	h_corrected->Draw("colz");
411		411
412	TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);	412	TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);
413	canvas13->Divide(2);	413	canvas13->Divide(2);
414	canvas13->cd(1);	414	canvas13->cd(1);
415	h[16]->Draw("colz");	415	h[16]->Draw("colz");
416	canvas13->cd(2);	416	canvas13->cd(2);
417	h[8]->Draw("colz");	417	h[8]->Draw("colz");
418	}	418	}
419		419
420	/** Draws the beam profile and fits it with error function	420	/** Draws the beam profile and fits it with error function
421	* on some background function	421	* on some background function
422	*/	422	*/
423	if (strstr(plopt, "beam") != NULL) {	423	if (strstr(plopt, "beam") != NULL) {
424		424
425	TCanvas* canvas9 = new TCanvas("canvas9","canvas9", 500,500);	425	TCanvas* canvas9 = new TCanvas("canvas9","canvas9", 500,500);
426	canvas9->cd();	426	canvas9->cd();
427	char hname[128];	427	char hname[128];
428	sprintf(hname, "hnhitsx%d", 36);	428	sprintf(hname, "hnhitsx%d", 36);
429	TH1F* h_laser = (TH1F*)rootfile->Get(hname);	429	TH1F* h_laser = (TH1F*)rootfile->Get(hname);
430	h_laser->Draw();	430	h_laser->Draw();
431	h_laser->SetStats(1);	431	h_laser->SetStats(1);
432		432
433	TF1* err = new TF1("err","[0]+[1]*TMath::Erf((x-[2])/[3])",17.2,17.50);	433	TF1* err = new TF1("err","[0]+[1]*TMath::Erf((x-[2])/[3])",17.2,17.50);
434	err->SetParameter(0,2500);	434	err->SetParameter(0,2500);
435	err->SetParameter(1, h_laser->GetMaximum());	435	err->SetParameter(1, h_laser->GetMaximum());
436	err->SetParameter(2, h_laser->GetBinCenter(h_laser->GetMaximumBin()));	436	err->SetParameter(2, h_laser->GetBinCenter(h_laser->GetMaximumBin()));
437	err->SetParameter(3, 0.001);	437	err->SetParameter(3, 0.001);
438	h_laser->Fit(err,"qr");	438	h_laser->Fit(err,"qr");
439	h_laser->Fit(err,"lr");	439	h_laser->Fit(err,"lr");
440	double sigma = err->GetParameter(3);	440	double sigma = err->GetParameter(3);
441	printf("sigma = %2.0f um, FWHM = %2.0f um\n", sigma1000, 2.35sigma*1000);	441	printf("sigma = %2.0f um, FWHM = %2.0f um\n", sigma1000, 2.35sigma*1000);
442	}	442	}
443		443
444	if (strstr(plopt, "map") != NULL) {	444	if (strstr(plopt, "map") != NULL) {
445	for (int i=7; i>=0; i--) {	445	for (int i=7; i>=0; i--) {
446	//for (int j=7; j>=0; j--) printf("(%d, %d) ", j,i);	446	//for (int j=7; j>=0; j--) printf("(%d, %d) ", j,i);
447	for (int j=7; j>=0; j--) printf("%2d ", map[j][i]);	447	for (int j=7; j>=0; j--) printf("%2d ", map[j][i]);
448	printf("\n");	448	printf("\n");
449	}	449	}
450		450
451	}	451	}
452		452
453	return(0);	453	return(0);
454	}	454	}
455		455
456	/** Function calculates the noise from one channel	456	/** Function calculates the noise from one channel
457	* it runs through the bins along x and returns the average value	457	* it runs through the bins along x and returns the average value
458	*/	458	*/
459	double getNoise(TH2F* histogram, int yStart, int yEnd)	459	double getNoise(TH2F* histogram, int yStart, int yEnd)
460	{	460	{
461	double noise=0;	461	double noise=0;
462	int count=0;	462	int count=0;
463	for(int j=yStart; j<yEnd; j++) {	463	for(int j=yStart; j<yEnd; j++) {
464	double value = histogram->GetBinContent(j,2);	464	double value = histogram->GetBinContent(j,2);
465	//if (noise < value) noise = value;	465	//if (noise < value) noise = value;
466	noise += value;	466	noise += value;
467	count++;	467	count++;
468	}	468	}
469	return (noise/double(count));	469	return (noise/double(count));
470	}	470	}
471		471
472		472

Subversion Repositories f9daq

(root)/sipmScan/examples/plotThresholdScan.C @ 40 – Rev 30 → 34