source: trunk/MagicSoft/Mars/macros/pedestalstudies.C@ 3831

Last change on this file since 3831 was 3775, checked in by gaug, 20 years ago
*** empty log message ***
File size: 18.0 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Markus Gaug, 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24//////////////////////////////////////////////////////////////////////////////
25//
26// pedestalstudies.C
27//
28// macro to observe the pedestals and pedestalRMS with the number of FADC
29// slices summed up.
30//
31// In order to use this macro, you have to uncomment the following
32// line in MPedCalcPedRun (line 214):
33//
34// fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1;
35//
36//
37/////////////////////////////////////////////////////////////////////////////////
38const TString pedfile = "./20040303_20123_P_NewCalBoxTestLidOpen_E.root";
39
40void pedestalstudies(const TString pedname=pedfile)
41{
42
43 Int_t loops = 13;
44 Int_t stepsize = 2;
45
46 gStyle->SetOptStat(1111);
47 gStyle->SetOptFit();
48
49 TArrayF *hmeandiffinn = new TArrayF(loops);
50 TArrayF *hrmsdiffinn = new TArrayF(loops);
51 TArrayF *hmeandiffout = new TArrayF(loops);
52 TArrayF *hrmsdiffout = new TArrayF(loops);
53 TArrayF *hmeaninn = new TArrayF(loops);
54 TArrayF *hmeanout = new TArrayF(loops);
55 TArrayF *hrmsinn = new TArrayF(loops);
56 TArrayF *hrmsout = new TArrayF(loops);
57 TArrayF *hmuinn = new TArrayF(loops);
58 TArrayF *hmuout = new TArrayF(loops);
59 TArrayF *hsigmainn = new TArrayF(loops);
60 TArrayF *hsigmaout = new TArrayF(loops);
61
62 TArrayF *hmeandiffinnerr = new TArrayF(loops);
63 TArrayF *hrmsdiffinnerr = new TArrayF(loops);
64 TArrayF *hmeandiffouterr = new TArrayF(loops);
65 TArrayF *hrmsdiffouterr = new TArrayF(loops);
66 TArrayF *hmeaninnerr = new TArrayF(loops);
67 TArrayF *hmeanouterr = new TArrayF(loops);
68 TArrayF *hrmsinnerr = new TArrayF(loops);
69 TArrayF *hrmsouterr = new TArrayF(loops);
70 TArrayF *hmuinnerr = new TArrayF(loops);
71 TArrayF *hmuouterr = new TArrayF(loops);
72 TArrayF *hsigmainnerr = new TArrayF(loops);
73 TArrayF *hsigmaouterr = new TArrayF(loops);
74
75
76 MStatusDisplay *display = new MStatusDisplay;
77 display->SetUpdateTime(500);
78 display->Resize(850,700);
79
80 //
81 // Create a empty Parameter List and an empty Task List
82 // The tasklist is identified in the eventloop by its name
83 //
84 MParList plist;
85 MTaskList tlist;
86 plist.AddToList(&tlist);
87
88 for (Int_t samples=2;samples<stepsize*loops+1;samples=samples+stepsize)
89 {
90
91 plist.Reset();
92 tlist.Reset();
93
94 //
95 // Now setup the tasks and tasklist for the pedestals:
96 // ---------------------------------------------------
97 //
98
99 MReadMarsFile read("Events", pedname);
100 read.DisableAutoScheme();
101
102 MGeomApply geomapl;
103 //
104 // Set the extraction range higher:
105 //
106 //sigcalc.SetRange(1,14,1,14);
107 MPedCalcPedRun pedcalc;
108 pedcalc.SetNumHiGainSamples((Byte_t)samples);
109
110 MExtractSignal sigcalc;
111 sigcalc.SetRange(0,samples-1,0,samples-1);
112
113 //
114 // Additionally to calculating the pedestals,
115 // you can fill histograms and look at them
116 //
117 MFillH fill("MHPedestalCam", "MExtractedSignalCam");
118 fill.SetNameTab(Form("%s%2d","PedCam",samples));
119
120 tlist.AddToList(&read);
121 tlist.AddToList(&geomapl);
122 tlist.AddToList(&sigcalc);
123 tlist.AddToList(&pedcalc);
124 tlist.AddToList(&fill);
125
126 MGeomCamMagic geomcam;
127 MPedestalCam pedcam;
128 MBadPixelsCam badcam;
129 badcam.AsciiRead("badpixels.dat");
130
131 MHPedestalCam hpedcam;
132 MCalibrationPedCam cpedcam;
133
134 plist.AddToList(&geomcam);
135 plist.AddToList(&pedcam);
136 plist.AddToList(&hpedcam);
137 plist.AddToList(&cpedcam);
138 plist.AddToList(&badcam);
139
140 //
141 // Create and setup the eventloop
142 //
143 MEvtLoop evtloop;
144
145 evtloop.SetParList(&plist);
146 evtloop.SetDisplay(display);
147
148 //
149 // Execute first analysis
150 //
151 if (!evtloop.Eventloop())
152 return;
153
154 //
155 // Look at one specific pixel, after all the histogram manipulations:
156 //
157 /*
158 MHGausEvents &hpix = hpedcam.GetAverageHiGainArea(0);
159 hpix.DrawClone("fourierevents");
160
161 MHGausEvents &lpix = hpedcam.GetAverageHiGainArea(1);
162 lpix.DrawClone("fourierevents");
163
164 hpedcam[170].DrawClone("fourierevents");
165
166 */
167
168 MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice");
169 MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice");
170 MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice");
171 MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice");
172 MHCamera dispped9 (geomcam, "Ped;DeltaPedMean", "Rel. Diff. Mean per Slice (Fit-Calc.)");
173 MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Fit-Calc.)");
174
175 dispped0.SetCamContent( pedcam, 0);
176 dispped0.SetCamError( pedcam, 1);
177 dispped2.SetCamContent( pedcam, 2);
178 dispped2.SetCamError( pedcam, 3);
179
180 dispped4.SetCamContent( hpedcam, 0);
181 dispped4.SetCamError( hpedcam, 1);
182 dispped6.SetCamContent( hpedcam, 2);
183 dispped6.SetCamError( hpedcam, 3);
184 dispped9.SetCamContent( hpedcam, 5);
185 dispped9.SetCamError( hpedcam, 6);
186 dispped11.SetCamContent(hpedcam, 8);
187 dispped11.SetCamError( hpedcam, 9);
188
189 dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
190 dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
191 dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
192 dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
193 dispped9.SetYTitle("Rel. Diff. Pedestal per slice Fit-Calc [1]");
194 dispped11.SetYTitle("Rel. Diff. Pedestal RMS per slice Fit-Calc [1]");
195
196
197 // Histogram values
198 TCanvas &b1 = display->AddTab(Form("%s%d","MeanRMS",samples));
199 b1.Divide(4,3);
200
201 CamDraw(b1,dispped0,1,4,*hmeaninn,*hmeanout,*hmeaninnerr,*hmeanouterr,samples,stepsize);
202 CamDraw(b1,dispped2,2,4,*hrmsinn,*hrmsout,*hrmsinnerr,*hrmsouterr,samples,stepsize);
203 CamDraw(b1,dispped4,3,4,*hmuinn,*hmuout,*hmuinnerr,*hmuouterr,samples,stepsize);
204 CamDraw(b1,dispped6,4,4,*hsigmainn,*hsigmaout,*hsigmainnerr,*hsigmaouterr,samples,stepsize);
205
206 display->SaveAsGIF(3*((samples-1)/stepsize)+2,Form("%s%d","MeanRmsSamples",samples));
207
208 // Differences
209 TCanvas &c4 = display->AddTab(Form("%s%d","RelDiff",samples));
210 c4.Divide(2,3);
211
212 CamDraw(c4,dispped9,1,2,*hmeandiffinn,*hmeandiffout,*hmeandiffinnerr,*hmeandiffouterr,samples,stepsize);
213 CamDraw(c4,dispped11,2,2,*hrmsdiffinn,*hrmsdiffout,*hrmsdiffinnerr,*hrmsdiffouterr,samples,stepsize);
214
215 display->SaveAsGIF(3*((samples-1)/stepsize)+3,Form("%s%d","RelDiffSamples",samples));
216
217 }
218
219 TF1 *logg = new TF1("logg","[1]+TMath::Log(x-[0])",1.,30.,2);
220 logg->SetParameters(1.,3.5);
221 logg->SetParLimits(0,-1.,3.);
222 logg->SetParLimits(1,-1.,7.);
223 logg->SetLineColor(kRed);
224
225 TCanvas *canvas = new TCanvas("PedstudInner","Pedestal Studies Inner Pixels",600,900);
226 canvas->Divide(2,3);
227 canvas->cd(1);
228
229 TGraphErrors *gmeaninn = new TGraphErrors(hmeaninn->GetSize(),
230 CreateXaxis(hmeaninn->GetSize(),stepsize),hmeaninn->GetArray(),
231 CreateXaxisErr(hmeaninnerr->GetSize(),stepsize),hmeaninnerr->GetArray());
232 gmeaninn->Draw("A*");
233 gmeaninn->SetTitle("Calculated Mean per Slice Inner Pixels");
234 gmeaninn->GetXaxis()->SetTitle("Nr. added FADC slices");
235 gmeaninn->GetYaxis()->SetTitle("Calculated Mean per slice");
236 gmeaninn->Fit("pol0");
237 gmeaninn->GetFunction("pol0")->SetLineColor(kGreen);
238 // gmeaninn->Fit(logg);
239
240 canvas->cd(2);
241
242 TGraphErrors *gmuinn = new TGraphErrors(hmuinn->GetSize(),
243 CreateXaxis(hmuinn->GetSize(),stepsize),hmuinn->GetArray(),
244 CreateXaxisErr(hmuinnerr->GetSize(),stepsize),hmuinnerr->GetArray());
245 gmuinn->Draw("A*");
246 gmuinn->SetTitle("Fitted Mean per Slice Inner Pixels");
247 gmuinn->GetXaxis()->SetTitle("Nr. added FADC slices");
248 gmuinn->GetYaxis()->SetTitle("Fitted Mean per Slice");
249 gmuinn->Fit("pol0");
250 gmuinn->GetFunction("pol0")->SetLineColor(kGreen);
251 //gmuinn->Fit(logg);
252
253
254 canvas->cd(3);
255
256 TGraphErrors *grmsinn = new TGraphErrors(hrmsinn->GetSize(),
257 CreateXaxis(hrmsinn->GetSize(),stepsize),hrmsinn->GetArray(),
258 CreateXaxisErr(hrmsinnerr->GetSize(),stepsize),hrmsinnerr->GetArray());
259 grmsinn->Draw("A*");
260 grmsinn->SetTitle("Calculated Rms per Slice Inner Pixels");
261 grmsinn->GetXaxis()->SetTitle("Nr. added FADC slices");
262 grmsinn->GetYaxis()->SetTitle("Calculated Rms per Slice");
263 //grmsinn->Fit("pol2");
264 //grmsinn->GetFunction("pol2")->SetLineColor(kRed);
265 grmsinn->Fit(logg);
266
267 canvas->cd(4);
268
269 TGraphErrors *gsigmainn = new TGraphErrors(hsigmainn->GetSize(),
270 CreateXaxis(hsigmainn->GetSize(),stepsize),hsigmainn->GetArray(),
271 CreateXaxisErr(hsigmainnerr->GetSize(),stepsize),hsigmainnerr->GetArray());
272 gsigmainn->Draw("A*");
273 gsigmainn->SetTitle("Fitted Sigma per Slice Inner Pixels");
274 gsigmainn->GetXaxis()->SetTitle("Nr. added FADC slices");
275 gsigmainn->GetYaxis()->SetTitle("Fitted Sigma per Slice");
276 // gsigmainn->Fit("pol2");
277 // gsigmainn->GetFunction("pol2")->SetLineColor(kRed);
278 gsigmainn->Fit(logg);
279
280 canvas->cd(5);
281
282 TGraphErrors *gmeandiffinn = new TGraphErrors(hmeandiffinn->GetSize(),
283 CreateXaxis(hmeandiffinn->GetSize(),stepsize),hmeandiffinn->GetArray(),
284 CreateXaxisErr(hmeandiffinnerr->GetSize(),stepsize),hmeandiffinnerr->GetArray());
285 gmeandiffinn->Draw("A*");
286 gmeandiffinn->SetTitle("Rel. Difference Mean per Slice Inner Pixels");
287 gmeandiffinn->GetXaxis()->SetTitle("Nr. added FADC slices");
288 gmeandiffinn->GetYaxis()->SetTitle("Rel. Difference Mean per Slice");
289 //gmeandiffinn->Fit("pol2");
290 //gmeandiffinn->GetFunction("pol2")->SetLineColor(kBlue);
291 gmeandiffinn->Fit(logg);
292
293
294 canvas->cd(6);
295
296 TGraphErrors *grmsdiffinn = new TGraphErrors(hrmsdiffinn->GetSize(),
297 CreateXaxis(hrmsdiffinn->GetSize(),stepsize),hrmsdiffinn->GetArray(),
298 CreateXaxisErr(hrmsdiffinnerr->GetSize(),stepsize),hrmsdiffinnerr->GetArray());
299 grmsdiffinn->Draw("A*");
300 grmsdiffinn->SetTitle("Rel. Difference Sigma per Slice-RMS Inner Pixels");
301 grmsdiffinn->GetXaxis()->SetTitle("Nr. added FADC slices");
302 grmsdiffinn->GetYaxis()->SetTitle("Rel. Difference Sigma per Slice-RMS");
303 //grmsdiffinn->Fit("pol2");
304 //grmsdiffinn->GetFunction("pol2")->SetLineColor(kBlue);
305 grmsdiffinn->Fit(logg);
306
307
308 TCanvas *canvas2 = new TCanvas("PedstudOut","Pedestal Studies Outer Pixels",600,900);
309 canvas2->Divide(2,3);
310 canvas2->cd(1);
311
312
313 canvas2->cd(1);
314
315 TGraphErrors *gmeanout = new TGraphErrors(hmeanout->GetSize(),
316 CreateXaxis(hmeanout->GetSize(),stepsize),hmeanout->GetArray(),
317 CreateXaxisErr(hmeanouterr->GetSize(),stepsize),hmeanouterr->GetArray());
318 gmeanout->Draw("A*");
319 gmeanout->SetTitle("Calculated Mean per Slice Outer Pixels");
320 gmeanout->GetXaxis()->SetTitle("Nr. added FADC slices");
321 gmeanout->GetYaxis()->SetTitle("Calculated Mean per Slice");
322 gmeanout->Fit("pol0");
323 gmeanout->GetFunction("pol0")->SetLineColor(kGreen);
324 //gmeanout->Fit(logg);
325
326 canvas2->cd(2);
327
328 TGraphErrors *gmuout = new TGraphErrors(hmuout->GetSize(),
329 CreateXaxis(hmuout->GetSize(),stepsize),hmuout->GetArray(),
330 CreateXaxisErr(hmuouterr->GetSize(),stepsize),hmuouterr->GetArray());
331 gmuout->Draw("A*");
332 gmuout->SetTitle("Fitted Mean per Slice Outer Pixels");
333 gmuout->GetXaxis()->SetTitle("Nr. added FADC slices");
334 gmuout->GetYaxis()->SetTitle("Fitted Mean per Slice");
335 gmuout->Fit("pol0");
336 gmuout->GetFunction("pol0")->SetLineColor(kGreen);
337 //gmuout->Fit(logg);
338
339 canvas2->cd(3);
340
341 TGraphErrors *grmsout = new TGraphErrors(hrmsout->GetSize(),
342 CreateXaxis(hrmsout->GetSize(),stepsize),hrmsout->GetArray(),
343 CreateXaxisErr(hrmsouterr->GetSize(),stepsize),hrmsouterr->GetArray());
344 grmsout->Draw("A*");
345 grmsout->SetTitle("Calculated Rms per Slice Outer Pixels");
346 grmsout->GetXaxis()->SetTitle("Nr. added FADC slices");
347 grmsout->GetYaxis()->SetTitle("Calculated Rms per Slice");
348 //grmsout->Fit("pol2");
349 //grmsout->GetFunction("pol2")->SetLineColor(kRed);
350 grmsout->Fit(logg);
351
352 canvas2->cd(4);
353
354 TGraphErrors *gsigmaout = new TGraphErrors(hsigmaout->GetSize(),
355 CreateXaxis(hsigmaout->GetSize(),stepsize),hsigmaout->GetArray(),
356 CreateXaxisErr(hsigmaouterr->GetSize(),stepsize),hsigmaouterr->GetArray());
357 gsigmaout->Draw("A*");
358 gsigmaout->SetTitle("Fitted Sigma per Slice Outer Pixels");
359 gsigmaout->GetXaxis()->SetTitle("Nr. added FADC slices");
360 gsigmaout->GetYaxis()->SetTitle("Fitted Sigma per Slice");
361 //gsigmaout->Fit("pol2");
362 //gsigmaout->GetFunction("pol2")->SetLineColor(kRed);
363 gsigmaout->Fit(logg);
364
365
366 canvas2->cd(5);
367
368 TGraphErrors *gmeandiffout = new TGraphErrors(hmeandiffout->GetSize(),
369 CreateXaxis(hmeandiffout->GetSize(),stepsize),hmeandiffout->GetArray(),
370 CreateXaxisErr(hmeandiffouterr->GetSize(),stepsize),hmeandiffouterr->GetArray());
371 gmeandiffout->Draw("A*");
372 gmeandiffout->SetTitle("Rel. Difference Mean per Slice Outer Pixels");
373 gmeandiffout->GetXaxis()->SetTitle("Nr. added FADC slices");
374 gmeandiffout->GetYaxis()->SetTitle("Rel. Difference Mean per Slice");
375 //gmeandiffout->Fit("pol2");
376 //gmeandiffout->GetFunction("pol2")->SetLineColor(kBlue);
377 gmeandiffout->Fit(logg);
378
379 canvas2->cd(6);
380
381 TGraphErrors *grmsdiffout = new TGraphErrors(hrmsdiffout->GetSize(),
382 CreateXaxis(hrmsdiffout->GetSize(),stepsize),hrmsdiffout->GetArray(),
383 CreateXaxisErr(hrmsdiffouterr->GetSize(),stepsize),hrmsdiffouterr->GetArray());
384 grmsdiffout->Draw("A*");
385 grmsdiffout->SetTitle("Rel. Difference Sigma per Slice-RMS Outer Pixels");
386 grmsdiffout->GetXaxis()->SetTitle("Nr. added FADC slices");
387 grmsdiffout->GetYaxis()->SetTitle("Rel. Difference Sigma per Slice-RMS");
388 //grmsdiffout->Fit("pol2");
389 //grmsdiffout->GetFunction("pol2")->SetLineColor(kBlue);
390 grmsdiffout->Fit(logg);
391
392
393}
394
395
396void CamDraw(TCanvas &c, MHCamera &cam, Int_t i, Int_t j, TArrayF &a1, TArrayF &a2,
397 TArrayF &a1err, TArrayF &a2err, Int_t samp, Int_t stepsize)
398{
399
400 c.cd(i);
401 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
402 obj1->SetDirectory(NULL);
403
404 c.cd(i+j);
405 obj1->Draw();
406 ((MHCamera*)obj1)->SetPrettyPalette();
407
408 c.cd(i+2*j);
409 TH1D *obj2 = (TH1D*)obj1->Projection();
410 obj2->SetDirectory(NULL);
411
412 // obj2->Sumw2();
413 obj2->Draw();
414 obj2->SetBit(kCanDelete);
415
416 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
417 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
418 const Double_t integ = obj2->Integral("width")/2.5066283;
419 const Double_t mean = obj2->GetMean();
420 const Double_t rms = obj2->GetRMS();
421 const Double_t width = max-min;
422
423 if (rms == 0. || width == 0. )
424 return;
425
426 TArrayI s0(6);
427 s0[0] = 6;
428 s0[1] = 1;
429 s0[2] = 2;
430 s0[3] = 3;
431 s0[4] = 4;
432 s0[5] = 5;
433
434 TArrayI inner(1);
435 inner[0] = 0;
436
437 TArrayI outer(1);
438 outer[0] = 1;
439
440 // Just to get the right (maximum) binning
441 TH1D *half[2];
442 half[0] = obj1->ProjectionS(s0, inner, "Inner");
443 half[1] = obj1->ProjectionS(s0, outer, "Outer");
444
445 half[0]->SetDirectory(NULL);
446 half[1]->SetDirectory(NULL);
447
448 for (int i=0; i<2; i++)
449 {
450 half[i]->SetLineColor(kRed+i);
451 half[i]->SetDirectory(0);
452 half[i]->SetBit(kCanDelete);
453 half[i]->Draw("same");
454 half[i]->Fit("gaus","Q+");
455
456 if (i==0)
457 {
458 a1[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParameter(1);
459 a1err[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParError(1);
460 if (a1err[(samp-1)/stepsize] > 3.)
461 a1err[(samp-1)/stepsize] = 1.;
462 }
463 if (i==1)
464 {
465 a2[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParameter(1);
466 a2err[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParError(1);
467 if (a2err[(samp-1)/stepsize] > 3.)
468 a2err[(samp-1)/stepsize] = 1.;
469 }
470 }
471
472
473}
474
475// -----------------------------------------------------------------------------
476//
477// Create the x-axis for the event graph
478//
479Float_t *CreateXaxis(Int_t n, Int_t step)
480{
481
482 Float_t *xaxis = new Float_t[n];
483
484 for (Int_t i=0;i<n;i++)
485 xaxis[i] = 2. + step*i;
486
487 return xaxis;
488
489}
490
491// -----------------------------------------------------------------------------
492//
493// Create the x-axis for the event graph
494//
495Float_t *CreateXaxisErr(Int_t n, Int_t step)
496{
497
498 Float_t *xaxis = new Float_t[n];
499
500 for (Int_t i=0;i<n;i++)
501 xaxis[i] = step/2.;
502
503 return xaxis;
504
505}
Note: See TracBrowser for help on using the repository browser.