Subversion Repositories f9daq

Rev

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