source: branches/removing_cpp11_features/macros/pedestalstudies.C@ 19693

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