Subversion Repositories f9daq

Rev

Rev 50 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
35 f9daq 1
#include "TROOT.h"
2
#include "TFile.h"
3
#include "TBenchmark.h"
4
#include "TH1F.h"
5
#include "TH2F.h"
6
#include "TH3F.h"
7
#include "TCanvas.h"
8
#include "TStyle.h"
9
#include "TPad.h"
10
#include "TF1.h"
11
#include "TGraph.h"
12
#include "TSpectrum.h"
13
#include "stdio.h"
38 f9daq 14
#include "THStack.h"
40 f9daq 15
#include "TPaveText.h"
35 f9daq 16
 
17
#include "include/RTUtil.h"
18
 
19
double getNoise(TH2F*, int, int);
20
 
40 f9daq 21
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)
35 f9daq 22
{
23
  const int c_nChannels = 64;
40 f9daq 24
  const double c_xOffset = 0; // mm
25
  const double c_yOffset = 0;
35 f9daq 26
 
27
  int map[8][8]={{32,34,53,55,40,42,61,63},
28
                 {48,50,37,39,56,58,45,47},
29
                 {33,35,52,54,41,43,60,62},
30
                 {49,51,36,38,57,59,44,46},
31
                 {17,19,4,6,25,27,12,14},
32
                 {1,3,20,22,9,11,28,30},
33
                 {16,18,5,7,24,26,13,15},
34
                 {0,2,21,23,8,10,29,31}
35
                };
36
 
38 f9daq 37
  // Set draw style
35 f9daq 38
  RTSetStyle(gStyle);
39
 
40
  // open the file with histograms
41
  char fnameroot[256];
42
  TFile* rootfile;
43
	sprintf(fnameroot, "root/%s.root", filename);
44
	rootfile = (TFile *) gROOT->FindObject(filename);
45
	if(rootfile==NULL) rootfile = new TFile(fnameroot);
46
	if(rootfile==NULL) {
47
	  printf("Cannot open root file %s!!!\n",fnameroot);
48
	  return(0);
49
	}
50
 
51
  if(strstr(plopt, "tdc") != NULL) {
52
    TCanvas *canvas21 = new TCanvas("canvas21","canvas21",1600,800);
53
    TH3F* h0 = (TH3F*) rootfile->Get("h_correctedTDC");
54
	  Int_t binsX = h0->GetXaxis()->GetNbins();
55
    Int_t minX  = h0->GetXaxis()->GetFirst();
56
    Int_t maxX  = h0->GetXaxis()->GetLast()+1;
57
    Int_t binsY = h0->GetYaxis()->GetNbins();
58
    Int_t minY  = h0->GetYaxis()->GetFirst();
59
    Int_t maxY  = h0->GetYaxis()->GetLast()+1;
60
    Int_t binsZ = h0->GetZaxis()->GetNbins();
61
    Int_t minZ  = h0->GetZaxis()->GetFirst();
62
    Int_t maxZ  = h0->GetZaxis()->GetLast()+1;
63
    Double_t xLowUser  = h0->GetXaxis()->GetBinLowEdge(minX);
64
    Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
65
    Double_t yLowUser  = h0->GetYaxis()->GetBinLowEdge(minY);
66
    Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
67
 
68
    canvas21->Divide(2);
69
    //h0->Draw();
70
    canvas21->cd(1);
71
    //h0->Project3D("xy")->Draw("colz");
72
    //canvas21->cd(2);
73
    //h0->Project3D("yz")->Draw("colz");
74
    //canvas21->cd(3);
75
    TH2D* h_correctedTDC = (TH2D*) h0->Project3D("xz");
76
    //h0->GetZaxis()->SetRangeUser(-5,5);
77
    //h0->Project3D("xzo")->Draw("colz");
40 f9daq 78
    //h_correctedTDC->SetTitle("; t [ns]; Channel");
35 f9daq 79
    h_correctedTDC->Draw("colz");
80
    /*
81
    TH1F* tdc1 = new TH1F("tdc1",";TDC [ns];Events",binsZ, minZ, maxZ);
82
    //for(int j=minY; j<maxY; j++) {
83
      for(int k=minZ; k<maxZ; k++) {
84
        double signal = h0->GetBinContent(1,1,k);
85
        tdc1->Fill(k,signal);
86
        //}
87
        }
88
    tdc1->Draw();
89
    */
90
 
91
    TH2D* h_timeWalk = (TH2D*) h0->Project3D("zy");
92
    canvas21->cd(2);
40 f9daq 93
    //h_timeWalk->SetTitle(";Threshold [V]; t [ns]");
35 f9daq 94
    h_timeWalk->Draw("colz");
95
  }
96
 
36 f9daq 97
  if (strstr(plopt, "bias") != NULL) {
98
    TCanvas* canvas01 = new TCanvas("canvas01","",800, 800);
99
    canvas01->cd();
100
    TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");
101
    h_bias->Draw("colz");
50 f9daq 102
 
103
    TCanvas* canvas02 = new TCanvas("canvas02","",800, 800);
104
    //canvas02->Divide(2,2);
105
    canvas02->cd();
106
    TH3F* h_biasThreshold = (TH3F*) rootfile->Get("h_biasThreshold");
107
    //h_biasThreshold->Draw();
108
    //canvas02->cd(2);
109
    //h_biasThreshold->GetXaxis()->SetRange(1,5);
110
    h_biasThreshold->Project3D("zy");
91 f9daq 111
    h_biasThreshold->SetTitle("");
50 f9daq 112
    h_biasThreshold_zy->DrawCopy("colz");
113
 
114
    TCanvas* canvas03 = new TCanvas("canvas03","",1200, 1200);
115
    canvas03->Divide(4,4);
116
    int x1=0;
117
    int x2=3;
118
    int y1=0;
119
    int y2=3;
120
    for(int j=x1; j<=x2; j++){
121
    	for(int k=y1; k<=y2; k++){
122
    		int channel = map[j][k];
123
    		int canvasPosition = 4*(4-k-1) + 4 - j;
124
    		canvas03->cd(canvasPosition);
125
    		h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);
126
    		char histoName[128];
127
    		sprintf(histoName, "Ch %d", channel);
128
    		h_biasThreshold->Project3D("zy");
129
    		//sprintf(option, "h_biasThresholsd_zy%f", channel);
130
    		h_biasThreshold_zy->SetTitle(histoName);
131
    		h_biasThreshold_zy->DrawCopy("colz");
132
 
133
    	}
134
    }
135
    TCanvas* canvas04 = new TCanvas("canvas04","",1200, 1200);
136
    canvas04->Divide(4,4);
137
    int x1=0;
138
        int x2=3;
139
        int y1=4;
140
        int y2=7;
141
        for(int j=x1; j<=x2; j++){
142
        	for(int k=y1; k<=y2; k++){
143
        		int channel = map[j][k];
144
        		int canvasPosition = 4*(8-k-1) + 4 - j;
145
        		canvas04->cd(canvasPosition);
146
        		h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);
147
        		char histoName[128];
148
        		sprintf(histoName, "Ch %d", channel);
149
        		h_biasThreshold->Project3D("zy");
150
        		//sprintf(option, "h_biasThresholsd_zy%f", channel);
151
        		h_biasThreshold_zy->SetTitle(histoName);
152
        		h_biasThreshold_zy->DrawCopy("colz");
153
 
154
        	}
155
        }
156
        TCanvas* canvas05 = new TCanvas("canvas05","",1200, 1200);
157
        canvas05->Divide(4,4);
158
        int x1=4;
159
        int x2=7;
160
        int y1=0;
161
        int y2=3;
162
        for(int j=x1; j<=x2; j++){
163
        	for(int k=y1; k<=y2; k++){
164
        		int channel = map[j][k];
165
        		int canvasPosition = 4*(4-k-1) + 8 - j;
166
        		canvas05->cd(canvasPosition);
167
        		h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);
168
        		char histoName[128];
169
        		sprintf(histoName, "Ch %d", channel);
170
        		h_biasThreshold->Project3D("zy");
171
        		//sprintf(option, "h_biasThresholsd_zy%f", channel);
172
        		h_biasThreshold_zy->SetTitle(histoName);
173
        		h_biasThreshold_zy->DrawCopy("colz");
174
 
175
        	}
176
        }
177
 
178
        TCanvas* canvas06 = new TCanvas("canvas06","",1200, 1200);
179
        canvas06->Divide(4,4);
180
        int x1=4;
181
        int x2=7;
182
        int y1=4;
183
        int y2=7;
184
        for(int j=x1; j<=x2; j++){
185
        	for(int k=y1; k<=y2; k++){
186
        		int channel = map[j][k];
187
        		int canvasPosition = 4*(8-k-1) + 8 - j;
188
        		canvas06->cd(canvasPosition);
189
        		h_biasThreshold->GetXaxis()->SetRange(channel+1,channel+1);
190
        		char histoName[128];
191
        		sprintf(histoName, "Ch %d", channel);
192
        		h_biasThreshold->Project3D("zy");
193
        		//sprintf(option, "h_biasThresholsd_zy%f", channel);
194
        		h_biasThreshold_zy->SetTitle(histoName);
195
        		h_biasThreshold_zy->DrawCopy("colz");
196
 
197
        	}
198
        }
199
 
36 f9daq 200
  }
201
 
202
  if (strstr(plopt, "th") != NULL) {
203
   TCanvas* canvas1 = new TCanvas("canvas1","",800,800);
204
   canvas1->cd();
205
   //gPad->SetLogz();
206
   //canvas1->SetLogz();
207
   TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
208
   h_threshold->SetTitleOffset(1.5,"y");
209
   //h_threshold->GetZaxis()->SetRangeUser(0,1000);
210
   h_threshold->SetContour(20);
211
 
212
   h_threshold->Draw("colz");
213
 
214
   TCanvas* canvas2 = new TCanvas("canvas2","",1600,800);
215
   canvas2->Divide(2);
216
   TH1D* h_projection1 = h_threshold->ProjectionY("Ch 37",38,38);
217
   canvas2->cd(1);
218
   //h_projection1->GetYaxis()->SetRangeUser(0,1000);
219
   //gPad->SetLogy();
220
   h_projection1->Draw();
221
 
222
   TH1D* h_projection2 = h_threshold->ProjectionY("Ch 38",39,39);
223
   canvas2->cd(2);
224
   //gPad->SetLogy();
91 f9daq 225
   h_projection2->Draw();
226
 
227
   TCanvas* canvas21 = new TCanvas("canvas21","",1200,1200);
228
   canvas21->Divide(3,3);
229
   TH3F* h_biasThreshold = (TH3F*) rootfile->Get("h_biasThreshold");
230
       //h_biasThreshold->Draw();
231
       //canvas02->cd(2);
232
       //h_biasThreshold->GetXaxis()->SetRange(1,5);
233
       h_biasThreshold->Project3D("zy");
234
       canvas21->cd(1);
235
       h_biasThreshold_zy->DrawCopy("colz");
236
       canvas21->cd(2);
237
       for (int jk=0; jk<9; jk++){
238
    	   h_biasThreshold_zy->ProjectionY("_py",jk*2+1,jk*2+1);
239
    	   canvas21->cd(jk+1);
240
    	   char histoName[128];
241
    	   double bias = h_biasThreshold_zy->GetXaxis()->GetBinCenter(jk*2+1);
242
    	   sprintf(histoName, "Bias = %4.2f V", bias);
243
    	   h_biasThreshold_zy_py->SetTitle(histoName);
244
    	   h_biasThreshold_zy_py->DrawCopy();
245
       }
246
 
36 f9daq 247
  }
248
 
35 f9daq 249
  if(strstr(plopt,"16") != NULL)
250
  {
36 f9daq 251
  TH2F* h_threshold = (TH2F*) rootfile->Get("h_bias");
35 f9daq 252
  TCanvas* canvas = new TCanvas("canvas","",1600,1600);
38 f9daq 253
  TCanvas* canvas11 = new TCanvas("canvas11","",800,800);
35 f9daq 254
  canvas->cd();
255
  TVirtualPad *main = new TPad("main","main",0,0,1,1,10,1);
256
  main->Draw();
38 f9daq 257
  main->cd();
35 f9daq 258
  main->Divide(4,4);
259
  for(int i=chXstart; i<chXend; i++) {
40 f9daq 260
    for(int j=(int)par1; j<(int)par2; j++) {
35 f9daq 261
      int channel = map[i][j];
262
      TH1D* h_projection = h_threshold->ProjectionY("",channel+1,channel+1);
40 f9daq 263
      int canvasPosition = i-chXend+4*((int)par2-j)+1;
35 f9daq 264
      printf(" %d ", canvasPosition);
265
      main->cd(canvasPosition);
38 f9daq 266
      gPad->SetMargin(0.08, 0.08, 0.08, 0.08);
35 f9daq 267
      char name[128];
268
      sprintf(name,"Channel %d",channel);
269
      h_projection->SetTitle(name);
38 f9daq 270
      h_projection->GetYaxis()->SetRangeUser(0,10000);
271
      h_projection->DrawCopy("pe1x0");
272
      canvas11->cd();
273
      if (canvasPosition == 1) {
274
        h_projection->SetMarkerStyle(kFullDotMedium);
275
        h_projection->GetYaxis()->SetRangeUser(0,10000);
276
        h_projection->DrawCopy("Pl");
277
        }
278
      else {
279
        h_projection->SetMarkerStyle(kFullDotMedium);
280
        h_projection->SetLineColor(canvasPosition);
281
        h_projection->SetMarkerColor(canvasPosition);
282
        h_projection->SetTitle("");
283
        h_projection->DrawCopy("plsame");
284
        }
35 f9daq 285
    }
38 f9daq 286
 
35 f9daq 287
  }
288
  printf("\n");
289
  }
290
 
36 f9daq 291
  if(strstr(plopt, "1ch") != NULL)
292
  {
293
    TH2F* h_bias = (TH2F*) rootfile->Get("h_bias");
38 f9daq 294
    TCanvas* canvas = new TCanvas("canvas","",800,800);
36 f9daq 295
    canvas->cd();
296
    int channel = map[chXstart][chXend];
297
    TH1D* h_projection = h_bias->ProjectionY("",channel+1,channel+1);
298
    char name[128];
299
    sprintf(name,"Channel %d",channel);
300
    h_projection->SetTitle(name);
38 f9daq 301
    h_projection->SetMarkerStyle(kFullDotMedium);
302
    //h_projection->SetMarkerSize(8);
303
    gPad->SetMargin(0.1,0.1,0.1,0.1);
40 f9daq 304
 
305
 
306
    TF1* f_linear = new TF1("f_linear", "[0] + [1]*x", par1, par2);
307
    f_linear->SetParameter(0, -1);
308
    f_linear->SetParameter(1, 0.1);
309
    f_linear->SetParNames("p0", "p1");
310
    h_projection->Fit(f_linear,"q");
311
    h_projection->Fit(f_linear,"r");
38 f9daq 312
    h_projection->Draw("pe1x0");
40 f9daq 313
 
314
    double p0 = f_linear->GetParameter(0);
315
    double p1 = f_linear->GetParameter(1);
316
    printf("p0 %f p1 %f V_b = %f [V]\n", p0, p1, -p0/p1);
317
    TPaveText *text = new TPaveText(0.15,0.5,0.35,0.56);
318
    text->SetFillColor(0);
319
    text->SetBorderSize(0.1);
320
    text->SetTextSize(0.030);
321
    char string[128];
322
    sprintf(string, "V_b = %2.2f", -p0/p1);
323
    text->AddText(string);
324
    text->Draw("same");
325
    text->Paint();
36 f9daq 326
  }
327
 
35 f9daq 328
  if( strstr(plopt, "all") != NULL ) {
329
    TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);
330
    TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);
331
    canvas2->Divide(8,8);
332
    canvas3->Divide(8,8);
333
    TH1F* h_hitsx;
334
    TH1F* h_hitsy;
335
    for(int i=0; i<c_nChannels; i++) {
336
      canvas2->cd(i+1);
337
      char hname[128];
338
      sprintf(hname, "hnhitsx%d", i);
339
      h_hitsx = (TH1F*)rootfile->Get(hname);
340
      h_hitsx->Draw();
341
      canvas3->cd(i+1);
342
      sprintf(hname, "hnhitsy%d", i);
343
      h_hitsy = (TH1F*)rootfile->Get(hname);
344
      h_hitsy->Draw();
345
    }
346
  }
347
 
348
  if( strstr(plopt, "x") != NULL ) {
349
    TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);
350
    TH1F* h_hitsx;
351
    canvas10->cd();
352
    char hname[128];
353
    sprintf(hname, "hnhitsx%d", chXstart);
354
    h_hitsx = (TH1F*)rootfile->Get(hname);
355
    h_hitsx->Draw("colz");
356
  }
357
 
358
  if( strstr(plopt, "y") != NULL ) {
359
    TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);
360
    TH1F* h_hitsy;
361
    canvas11->cd();
362
    char hname[128];
363
    sprintf(hname, "hnhitsy%d", chXstart);
364
    h_hitsy = (TH1F*)rootfile->Get(hname);
365
    h_hitsy->Draw();
366
  }
40 f9daq 367
 
35 f9daq 368
 
369
  /** Draws the sum of the channels
370
    * Each channel is a 2d plot
371
    * Intended for the study of 1 channel
372
    */
373
  if (strstr(plopt, "2d") != NULL) {
374
 
375
    int nX = chXend - chXstart + 1;
40 f9daq 376
    int nY = (int)par2 - (int)par1 + 1;
35 f9daq 377
    TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX*400,nY*400);
378
    printf("nx %d ny %d\n",nX,nY);
379
    canvas7->Divide(nX,nY);
40 f9daq 380
    for(int i=(int)par1; i<=(int)par2; i++) {
35 f9daq 381
      for(int j=chXstart; j<=chXend; j++) {
382
      //int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
40 f9daq 383
      int canvasPosition = nX*((int)par2-i)+(j-chXstart)+1;
35 f9daq 384
      if (debug) printf("canvas %d\n",canvasPosition);
385
        canvas7->cd(canvasPosition);
386
        char hname[128];
387
        int chPosition = map[j][i];
388
        sprintf(hname, "h2d%d", chPosition);
389
        TH2F* h_2d = (TH2F*)rootfile->Get(hname);
390
        h_2d->Draw("colz");
391
      } //x
392
    }
393
 
394
    // Number of photoelectrons - Poissonian correction
395
    TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);
396
    canvas8->cd();
397
    gStyle->SetOptStat(0);
398
      char hname[128];
399
      int chPosition = map[1][2];
400
      sprintf(hname, "h2d%d", chPosition);
401
      TH2F* h_2d = (TH2F*)rootfile->Get(hname);
402
 
403
      Int_t binsX = h_2d->GetXaxis()->GetNbins();
404
      Int_t minX  = h_2d->GetXaxis()->GetFirst();
405
      Int_t maxX  = h_2d->GetXaxis()->GetLast()+1;
406
      Int_t binsY = h_2d->GetYaxis()->GetNbins();
407
      Int_t minY  = h_2d->GetYaxis()->GetFirst();
408
      Int_t maxY  = h_2d->GetYaxis()->GetLast()+1;
409
      Double_t xLowUser  = h_2d->GetXaxis()->GetBinLowEdge(minX);
410
      Double_t xUpUser   = h_2d->GetXaxis()->GetBinUpEdge(maxX);
411
      Double_t yLowUser  = h_2d->GetYaxis()->GetBinLowEdge(minY);
412
      Double_t yUpUser =   h_2d->GetYaxis()->GetBinUpEdge(maxY);
413
      if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
414
      TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);
415
 
416
      double noise = getNoise(h_2d, 1, 89);
417
      if(debug) printf("Noise = %f\n", noise);
418
      for(int k=minX; k<=maxX; k++) {
419
        for(int j=minY; j<=maxY; j++) {
420
          double signal = h_2d->GetBinContent(k,j); // detected
421
          //double p = ((signal - noise) > 1) ? (signal-noise) : 1;
422
          double p = signal - noise;
423
          p /= 10000.;
424
          double p0 = 1.0 - p; // events with zero photons
425
          //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
426
          double eta = -log(p0); // constant of the poissonian statistics
427
          if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
428
          double x = k*(xUpUser-xLowUser)/double(binsX);
429
          double y = j*(yUpUser-yLowUser)/double(binsY);
430
          h_corrected->Fill(x+3,y+8, eta);
431
        }
432
       }
433
      h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");
434
      h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);
435
      h_corrected->Draw("colz");
436
 
437
      // collection efficiency
438
      int nPoints =0;
439
      double efficiency=0;
440
      for (int i=18; i<=58; i++) {
441
        for (int j=19; j<=59; j++) {
442
          double signal = h_corrected->GetBinContent(i,j);
443
          if(debug) printf("signal %f\n",signal);
444
          efficiency += signal;
445
          nPoints++;
446
        }
447
      }
448
    printf("Signal sum = %f\n   # of points = %d\n",efficiency,nPoints);
449
  }
450
 
451
  /** Draws the sum of channel signals
452
    * Each channel is a 2d ('h2d') histogram
453
    * Suitable for 8x8 chs scan
454
    */
455
  if( strstr(plopt, "sum") != NULL ) {
40 f9daq 456
    int nX = 7 - 0 + 1;
457
    int nY = 7 - 0 + 1;
35 f9daq 458
	  TCanvas* canvas12 = new TCanvas("canvas12","c2",8*200, 8*200);
459
	  canvas12->cd();
460
	  gStyle->SetOptStat(0);
461
 
462
    // final histogram parameters
463
    TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
464
	  Int_t binsX = h0->GetXaxis()->GetNbins();
465
    Int_t minX  = h0->GetXaxis()->GetFirst();
466
    Int_t maxX  = h0->GetXaxis()->GetLast()+1;
467
    Int_t binsY = h0->GetYaxis()->GetNbins();
468
    Int_t minY  = h0->GetYaxis()->GetFirst();
469
    Int_t maxY  = h0->GetYaxis()->GetLast()+1;
470
    Double_t xLowUser  = h0->GetXaxis()->GetBinLowEdge(minX);
471
    Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
472
    Double_t yLowUser  = h0->GetYaxis()->GetBinLowEdge(minY);
473
    Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
474
    if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
475
    TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
476
    TH2F* h[9];
477
 
478
    // 2d histogram noise subtraction and poisson scaling
40 f9daq 479
    for(int i=0; i<=7; i++) {
35 f9daq 480
      for(int j=chXstart; j<=chXend; j++) {
481
        int chPosition = map[j][i];
482
        char hname[128];
483
	      sprintf(hname, "h2d%d", chPosition);
40 f9daq 484
	      int histogram = nX*(i-0)+(j-chXstart);
35 f9daq 485
	      h[histogram] = (TH2F *) rootfile->Get(hname);
486
	      int noise = getNoise(h[histogram], 1, 170);
487
	      if (debug) printf("noise: %d\n",noise);
488
	      for(int k=minX; k<=maxX; k++) {
489
          for(int l=minY; l<=maxY; l++) {
490
            int signal = h[histogram]->GetBinContent(k,l); // detected
491
            //double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;
492
            double p = signal - noise;
493
            p /= 10000.;
494
            double p0 = 1.0 - p; // events with zero photons
495
            //double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
496
            double eta = -log(p0);
497
            //printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
498
            double x = k*(xUpUser-xLowUser)/double(binsX);
499
            double y = l*(yUpUser-yLowUser)/double(binsY);
500
            h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);
501
          }
502
        }
503
 
504
      }
505
    }
506
 
507
    h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");
508
    //h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);
509
    h_corrected->Draw("colz");
510
 
511
    TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);
512
    canvas13->Divide(2);
513
    canvas13->cd(1);
514
    //h[16]->Draw("colz");
515
    canvas13->cd(2);
516
    //h[8]->Draw("colz");
517
  }
518
 
519
  return(0);
520
}
521
 
522
/** Function calculates the noise from one channel
523
  * it runs through the bins along x and returns the average value
524
  */
525
double getNoise(TH2F* histogram, int yStart, int yEnd)
526
{
527
  double noise=0;
528
  int count=0;
529
  for(int j=yStart; j<yEnd; j++) {
530
    double value = histogram->GetBinContent(j,2);
531
    //if (noise < value) noise = value;
532
    noise += value;
533
    count++;
534
    }
535
  return (noise/double(count));
536
}
537