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 |