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