source: trunk/MagicSoft/Mars/macros/calibration.C@ 3188

Last change on this file since 3188 was 3185, checked in by gaug, 21 years ago
*** empty log message ***
File size: 23.3 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//const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root";
26//const TString calfile = "/mnt/users/mdoro/Mars/Data/20040201_1441*_C_OffMrk421-1_E.root";
27
28const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root";
29const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14608_C_CrabNebula_E.root";
30
31//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root";
32//const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_1041*_C_Crab-On_E.root";
33
34//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
35//const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root";
36
37void calibration(TString pedname=pedfile,
38 TString calname=calfile)
39{
40
41 gStyle->SetOptStat(1111);
42 gStyle->SetOptFit();
43
44 MStatusDisplay *display = new MStatusDisplay;
45 display->SetUpdateTime(3000);
46 display->Resize(850,700);
47
48 //
49 // Create a empty Parameter List and an empty Task List
50 // The tasklist is identified in the eventloop by its name
51 //
52 MParList plist;
53
54 MTaskList tlist;
55 plist.AddToList(&tlist);
56
57 //
58 // Now setup the tasks and tasklist for the pedestals:
59 // ---------------------------------------------------
60 //
61
62 MReadMarsFile read("Events", pedname);
63 read.DisableAutoScheme();
64
65 MGeomApply geomapl;
66 MExtractSignal sigcalc;
67
68
69 MExtractPINDiode pidcalc;
70 MExtractedSignalPINDiode pindiode;
71 // MHCalibrationChargePINDiode pinhist;
72 //
73 // Set the extraction range higher:
74 //
75 //sigcalc.SetRange(1,14,1,14);
76
77 MPedCalcPedRun pedcalc;
78
79 //
80 // Additionally to calculating the pedestals,
81 // you can fill histograms and look at them
82 //
83 MFillH fill("MHPedestalCam", "MExtractedSignalCam");
84 // MFillH fillp("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
85
86 tlist.AddToList(&read);
87 tlist.AddToList(&geomapl);
88 tlist.AddToList(&sigcalc);
89 tlist.AddToList(&pidcalc);
90 tlist.AddToList(&pedcalc);
91 tlist.AddToList(&fill);
92 // tlist.AddToList(&fillp);
93
94 MGeomCamMagic geomcam;
95 MPedestalCam pedcam;
96 MHPedestalCam hpedcam;
97 plist.AddToList(&geomcam);
98 plist.AddToList(&pedcam);
99 plist.AddToList(&hpedcam);
100 plist.AddToList(&pindiode);
101 // plist.AddToList(&pinhist);
102
103 //
104 // Create and setup the eventloop
105 //
106 MEvtLoop evtloop;
107 evtloop.SetParList(&plist);
108 evtloop.SetDisplay(display);
109
110 //
111 // Execute first analysis
112 //
113 if (!evtloop.Eventloop())
114 return;
115
116 tlist.PrintStatistics();
117
118 //
119 // Look at one specific pixel, after all the histogram manipulations:
120 //
121 hpedcam[9].DrawClone("fourier");
122
123
124 MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice");
125 MHCamera dispped1 (geomcam, "Ped;PedestalErr", "Mean Error per Slice");
126 MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice");
127 MHCamera dispped3 (geomcam, "Ped;PedestalRmsErr", "RMS Error per Slice");
128
129 MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice");
130 MHCamera dispped5 (geomcam, "Ped;MeanErr", "Fitted Error of Mean per Slice");
131 MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice");
132 MHCamera dispped7 (geomcam, "Ped;SigmaErr", "Fitted Error of Sigma per Slice");
133 MHCamera dispped8 (geomcam, "Ped;Prob", "Probability of Fit");
134 MHCamera dispped9 (geomcam, "Ped;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)");
135 MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)");
136 MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)");
137 MHCamera dispped12 (geomcam, "Ped;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)");
138 MHCamera dispped13 (geomcam, "Ped;FitOK", "Gaus Fit not OK");
139 MHCamera dispped14 (geomcam, "Ped;FourierOK", "Fourier Analysis not OK");
140
141 dispped0.SetCamContent( pedcam, 0);
142 dispped0.SetCamError( pedcam, 1);
143 dispped1.SetCamContent( pedcam, 1);
144 dispped2.SetCamContent( pedcam, 2);
145 dispped2.SetCamError( pedcam, 3);
146 dispped3.SetCamContent( pedcam, 3);
147
148 dispped4.SetCamContent( hpedcam, 0);
149 dispped4.SetCamError( hpedcam, 1);
150 dispped5.SetCamContent( hpedcam, 1);
151 dispped6.SetCamContent( hpedcam, 2);
152 dispped6.SetCamError( hpedcam, 3);
153 dispped7.SetCamContent( hpedcam, 3);
154 dispped8.SetCamContent( hpedcam, 4);
155 dispped9.SetCamContent( hpedcam, 5);
156 dispped9.SetCamError( hpedcam, 6);
157 dispped10.SetCamContent(hpedcam, 7);
158 dispped11.SetCamContent(hpedcam, 8);
159 dispped11.SetCamError( hpedcam, 9);
160 dispped12.SetCamContent(hpedcam, 10);
161 dispped13.SetCamContent(hpedcam, 11);
162 dispped14.SetCamContent(hpedcam, 12);
163
164 dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
165 dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]");
166 dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
167 dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]");
168 dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
169 dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]");
170 dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
171 dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]");
172 dispped8.SetYTitle("Fit Probability [1]");
173 dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]");
174 dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]");
175 dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]");
176 dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]");
177 dispped13.SetYTitle("[1]");
178 dispped14.SetYTitle("[1]");
179
180 // Histogram values
181 TCanvas &b1 = display->AddTab("Ped.Calc.");
182 b1.Divide(4,3);
183
184 CamDraw(b1,dispped0,pedcam,1,4,1);
185 CamDraw(b1,dispped1,pedcam,2,4,2);
186 CamDraw(b1,dispped2,pedcam,3,4,2);
187 CamDraw(b1,dispped3,pedcam,4,4,2);
188
189 // Fitted values
190 TCanvas &b2 = display->AddTab("Ped.Fit");
191 b2.Divide(4,3);
192
193 CamDraw(b2,dispped4,hpedcam,1,4,1);
194 CamDraw(b2,dispped5,hpedcam,2,4,2);
195 CamDraw(b2,dispped6,hpedcam,3,4,2);
196 CamDraw(b2,dispped7,hpedcam,4,4,2);
197
198
199 // Fits Probability
200 TCanvas &b3 = display->AddTab("Ped.Fit Prob.");
201 b3.Divide(1,3);
202
203 CamDraw(b3,dispped8,hpedcam,1,1,3);
204
205 // Differences
206 TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit");
207 c4.Divide(4,3);
208
209 CamDraw(c4,dispped9,hpedcam,1,4,1);
210 CamDraw(c4,dispped10,hpedcam,2,4,1);
211 CamDraw(c4,dispped11,hpedcam,3,4,1);
212 CamDraw(c4,dispped12,hpedcam,4,4,1);
213
214 // Defects
215 TCanvas &c5 = display->AddTab("Defects");
216 c5.Divide(2,2);
217
218 CamDraw(c5,dispped13,hpedcam,1,2,0);
219 CamDraw(c5,dispped14,hpedcam,2,2,0);
220
221 //
222 // Create a empty Parameter List and an empty Task List
223 //
224 MParList plist2;
225 MTaskList tlist2;
226 plist2.AddToList(&tlist2);
227
228 MExtractedSignalCam sigcam;
229 MCalibrationCam calcam;
230 MHCalibrationRelTimeCam timecam;
231 //
232 // Get the previously created MPedestalCam into the new Parameter List
233 //
234 plist2.AddToList(&geomcam);
235 plist2.AddToList(&pedcam);
236 plist2.AddToList(&sigcam);
237 plist2.AddToList(&calcam);
238 plist2.AddToList(&timecam);
239 //
240 // Get the MAGIC geometry
241 //
242 tlist2.AddToList(&geomapl);
243 //
244 // Now setup the new tasks and tasklist for the calibration
245 // ---------------------------------------------------
246 //
247
248 MReadMarsFile read2("Events", calname);
249 read2.DisableAutoScheme();
250
251 //
252 // We saw that the signal jumps between slices,
253 // thus take the sliding window
254 //
255 MExtractSignal2 sigcalc2;
256 MArrivalTimeCalc2 timecalc;
257 MCalibrationCalc calcalc;
258
259 MFillH filltime("MHCalibrationRelTimeCam", "MArrivalTime");
260 //
261 // Set the range (other than default)
262 // of FADC slices for the blind pixel
263 //
264 // calcalc.SetBlindPixelRange(10,25);
265
266 //
267 // Set the cut upon which a superposition of the blind pixel
268 // FADC slices will be filled into the SinglePHE histogram
269 //
270 calcalc.SetBlindPixelSinglePheCut(500);
271
272 //
273 // Skip the HiGain vs. LoGain calibration
274 //
275 calcalc.SkipHiLoGainCalibration();
276
277 //
278 // As long, as we don't have digital modules,
279 // we have to set the color of the pulser LED by hand
280 //
281 calcalc.SetPulserColor(MCalibrationCalc::kECT1);
282
283 //
284 // In case, we want to exclude a pre-defined list of bad pixels:
285 // (This is a preliminary feature)
286 //
287 // calcalc.ExcludePixelsFromAsciiFile("badpixels.dat");
288
289 //
290 // In case, you want to skip the blind pixel method:
291 // (NOT RECOMMENDED!!!)
292 //
293 // calcalc.SkipBlindPixelFit();
294
295 //
296 // In case, you want to skip the quality checks
297 // (NOT RECOMMENDED!!!)
298 //
299 // calcalc.SkipQualityChecks();
300
301 //
302 // In case, we want to apply another fit function to the
303 // blind pixel
304 //
305 MCalibrationBlindPix *bp = calcam.GetBlindPixel();
306 bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson5);
307
308 //
309 // Apply a filter against cosmics
310 // (was directly in MCalibrationCalc in earlier versions)
311 //
312 MFCosmics cosmics;
313 MContinue cont(&cosmics);
314
315 tlist2.AddToList(&read2);
316 tlist2.AddToList(&sigcalc2);
317 //
318 // In case, you want to skip the cosmics rejection,
319 // uncomment the next line
320 //
321 tlist2.AddToList(&cont);
322 //
323 // In case, you want to skip the somewhat lengthy calculation
324 // of the arrival times using a spline, uncomment the next two lines
325 //
326 tlist2.AddToList(&timecalc);
327 tlist2.AddToList(&filltime);
328 //
329 tlist2.AddToList(&calcalc);
330 //
331 // Create and setup the eventloop
332 //
333 MEvtLoop evtloop2;
334 evtloop2.SetParList(&plist2);
335 evtloop2.SetDisplay(display);
336
337 //
338 // Execute second analysis
339 //
340 if (!evtloop2.Eventloop())
341 return;
342
343 tlist2.PrintStatistics();
344
345 //
346 // print the most important results of all pixels
347 //
348 // calcam.Print();
349
350 //
351 // just one example how to get the plots of individual pixels
352 //
353 // calcam[563].DrawClone();
354 // calcam[564].DrawClone();
355
356 // Create histograms to display
357 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
358 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
359 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
360 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
361 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
362 MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)");
363 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)");
364 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
365 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)");
366 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)");
367 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
368 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)");
369 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)");
370 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
371 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
372 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
373 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
374 MHCamera disp18 (geomcam, "Cal;Oscillating", "Oscillating Pixels");
375 MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
376
377 MHCamera disp20 (geomcam, "Cal;Ped", "Pedestals");
378 MHCamera disp21 (geomcam, "Cal;PedRms", "Pedestal RMS");
379
380 MHCamera disp22 (geomcam, "time;Time", "Rel. Arrival Times");
381 MHCamera disp23 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times");
382 MHCamera disp24 (geomcam, "time;TimeProb", "Probability of Time Fit");
383 MHCamera disp25 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
384 MHCamera disp26 (geomcam, "time;Oscillating", "Oscillating Pixels");
385
386 MHCamera disp27 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
387 MHCamera disp28 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times");
388
389 // Fitted charge means and sigmas
390 disp1.SetCamContent(calcam, 0);
391 disp1.SetCamError( calcam, 1);
392 disp2.SetCamContent(calcam, 2);
393 disp2.SetCamError( calcam, 3);
394
395 // Fit probabilities
396 disp3.SetCamContent(calcam, 4);
397
398 // Reduced Sigmas and reduced sigmas per charge
399 disp4.SetCamContent(calcam, 5);
400 disp4.SetCamError( calcam, 6);
401 disp5.SetCamContent(calcam, 7);
402 disp5.SetCamError( calcam, 8);
403
404 // F-Factor Method
405 disp6.SetCamContent(calcam, 9);
406 disp6.SetCamError( calcam, 10);
407 disp7.SetCamContent(calcam, 11);
408 disp7.SetCamError( calcam, 12);
409 disp8.SetCamContent(calcam, 13);
410 disp8.SetCamError( calcam, 14);
411
412 /// Blind Pixel Method
413 disp9.SetCamContent(calcam, 15);
414 disp9.SetCamError( calcam, 16);
415 disp10.SetCamContent(calcam,17);
416 disp10.SetCamError( calcam,18);
417 disp11.SetCamContent(calcam,19);
418 disp11.SetCamError( calcam,20);
419
420 // PIN Diode Method
421 disp12.SetCamContent(calcam,21);
422 disp12.SetCamError( calcam,22);
423 disp13.SetCamContent(calcam,23);
424 disp13.SetCamError( calcam,24);
425 disp14.SetCamContent(calcam,25);
426 disp14.SetCamError( calcam,26);
427
428 // Pixels with defects
429 disp15.SetCamContent(calcam,27);
430 disp16.SetCamContent(calcam,28);
431 disp17.SetCamContent(calcam,29);
432 disp18.SetCamContent(calcam,30);
433
434 // Lo Gain calibration
435 disp19.SetCamContent(calcam,31);
436
437
438 // Pedestals
439 disp20.SetCamContent(calcam,35);
440 disp20.SetCamError( calcam,36);
441 disp21.SetCamContent(calcam,37);
442 disp21.SetCamError( calcam,38);
443
444
445 // Relative Times
446 disp22.SetCamContent(timecam,0);
447 disp22.SetCamError( timecam,1);
448 disp23.SetCamContent(timecam,2);
449 disp23.SetCamError( timecam,3);
450 disp24.SetCamContent(timecam,4);
451 disp25.SetCamContent(timecam,5);
452 disp26.SetCamContent(timecam,6);
453
454 // Absolute Times
455 disp27.SetCamContent(calcam,39);
456 disp27.SetCamError( calcam,40);
457 disp28.SetCamContent(calcam,41);
458 disp28.SetCamError( calcam,42);
459
460
461 disp1.SetYTitle("Charge [FADC units]");
462 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
463 disp3.SetYTitle("P_{Charge} [1]");
464
465 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
466 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
467
468 disp6.SetYTitle("Nr. Photo-Electrons [1]");
469 disp7.SetYTitle("Conversion Factor [PhE/FADC Count]");
470 disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
471
472 disp9.SetYTitle("Photon flux [ph/mm^2]");
473 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
474 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
475
476 disp12.SetYTitle("Photon flux [ph/mm^2]");
477 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
478 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
479
480 disp15.SetYTitle("[1]");
481 disp16.SetYTitle("[1]");
482 disp17.SetYTitle("[1]");
483 disp18.SetYTitle("[1]");
484 disp19.SetYTitle("[1]");
485
486 disp20.SetYTitle("Ped [FADC Counts ]");
487 disp21.SetYTitle("RMS_{Ped} [FADC Counts ]");
488
489 disp22.SetYTitle("Time Offset [ns]");
490 disp23.SetYTitle("Timing resolution [ns]");
491 disp24.SetYTitle("P_{Time} [1]");
492
493 disp25.SetYTitle("[1]");
494 disp26.SetYTitle("[1]");
495
496 disp27.SetYTitle("Mean Abs. Time [FADC slice]");
497 disp28.SetYTitle("RMS Abs. Time [FADC slices]");
498
499 gStyle->SetOptStat(1111);
500 gStyle->SetOptFit();
501
502 // Charges
503 TCanvas &c1 = display->AddTab("Fit.Charge");
504 c1.Divide(2, 3);
505
506 CamDraw(c1, disp1,calcam,1, 2 , 2);
507 CamDraw(c1, disp2,calcam,2, 2 , 2);
508
509 // Fit Probability
510 TCanvas &c2 = display->AddTab("Fit.Prob");
511 c2.Divide(1,3);
512
513 CamDraw(c2, disp3,calcam,1, 1 , 4);
514
515 // Reduced Sigmas
516 TCanvas &c3 = display->AddTab("Red.Sigma");
517 c3.Divide(2,3);
518
519 CamDraw(c3, disp4,calcam,1, 2 , 2);
520 CamDraw(c3, disp5,calcam,2, 2 , 2);
521
522 // F-Factor Method
523 TCanvas &c4 = display->AddTab("F-Factor");
524 c4.Divide(2,3);
525
526 CamDraw(c4, disp6,calcam,1, 2 , 2);
527 CamDraw(c4, disp7,calcam,2, 2 , 2);
528 CamDraw(c4, disp8,calcam,3, 3 , 2);
529
530 // Blind Pixel Method
531 TCanvas &c5 = display->AddTab("BlindPix");
532 c5.Divide(3, 3);
533
534 CamDraw(c5, disp9,calcam,1, 3 , 9);
535 CamDraw(c5, disp10,calcam,2, 3 , 2);
536 CamDraw(c5, disp11,calcam,3, 3 , 2);
537
538 // PIN Diode Method
539 TCanvas &c6 = display->AddTab("PINDiode");
540 c6.Divide(3,3);
541
542 CamDraw(c6, disp12,calcam,1, 3 , 9);
543 CamDraw(c6, disp13,calcam,2, 3 , 2);
544 CamDraw(c6, disp14,calcam,3, 3 , 2);
545
546 // Defects
547 TCanvas &c7 = display->AddTab("Defects");
548 c7.Divide(4,2);
549
550 CamDraw(c7, disp15,calcam,1,4, 0);
551 CamDraw(c7, disp16,calcam,2,4, 0);
552 CamDraw(c7, disp17,calcam,3,4, 0);
553 CamDraw(c7, disp18,calcam,4,4,0);
554
555 // Lo Gain Calibration
556 TCanvas &c8 = display->AddTab("LowGain");
557 c8.Divide(1,3);
558
559 CamDraw(c8, disp19,calcam,1,4,0);
560
561
562 // Pedestals
563 TCanvas &c9 = display->AddTab("Pedestals");
564 c9.Divide(2,3);
565
566 CamDraw(c9,disp20,calcam,1,2,1);
567 CamDraw(c9,disp21,calcam,2,2,2);
568
569
570 // Rel. Times
571 TCanvas &c10 = display->AddTab("Fitted Rel. Times");
572 c10.Divide(3,3);
573
574 CamDraw(c10,disp22,calcam,1,3,2);
575 CamDraw(c10,disp23,calcam,2,3,2);
576 CamDraw(c10,disp24,calcam,3,3,4);
577
578 // Time Defects
579 TCanvas &c11 = display->AddTab("Time Def.");
580 c11.Divide(2,2);
581
582 CamDraw(c11, disp25,calcam,1,2, 0);
583 CamDraw(c11, disp26,calcam,2,2, 0);
584
585 // Abs. Times
586 TCanvas &c12 = display->AddTab("Abs. Times");
587 c12.Divide(2,3);
588
589 CamDraw(c12,disp27,calcam,1,2,2);
590 CamDraw(c12,disp28,calcam,2,2,2);
591
592}
593
594void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
595{
596
597 c.cd(i);
598 gPad->SetBorderMode(0);
599 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
600 // obj1->AddNotify(evt);
601
602 c.cd(i+j);
603 gPad->SetBorderMode(0);
604 obj1->Draw();
605 ((MHCamera*)obj1)->SetPrettyPalette();
606
607 if (fit != 0)
608 {
609 c.cd(i+2*j);
610 gPad->SetBorderMode(0);
611 TH1D *obj2 = (TH1D*)obj1->Projection();
612
613 obj2->Draw();
614 obj2->SetBit(kCanDelete);
615
616
617 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
618 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
619 const Double_t integ = obj2->Integral("width")/2.5066283;
620 const Double_t mean = obj2->GetMean();
621 const Double_t rms = obj2->GetRMS();
622 const Double_t width = max-min;
623
624 if (rms == 0. || width == 0. )
625 return;
626
627 switch (fit)
628 {
629 case 1:
630 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
631 sgaus->SetBit(kCanDelete);
632 sgaus->SetParNames("Area","#mu","#sigma");
633 sgaus->SetParameters(integ/rms,mean,rms);
634 sgaus->SetParLimits(0,0.,integ);
635 sgaus->SetParLimits(1,min,max);
636 sgaus->SetParLimits(2,0,width/1.5);
637 obj2->Fit("sgaus","QLR");
638 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
639 break;
640
641 case 2:
642 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
643 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
644 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
645 dgaus->SetBit(kCanDelete);
646 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
647 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
648 integ/width/2.,(max+mean)/2.,width/4.);
649 // The left-sided Gauss
650 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
651 dgaus->SetParLimits(1,min+(width/10.),mean);
652 dgaus->SetParLimits(2,0,width/2.);
653 // The right-sided Gauss
654 dgaus->SetParLimits(3,0,integ);
655 dgaus->SetParLimits(4,mean,max-(width/10.));
656 dgaus->SetParLimits(5,0,width/2.);
657 obj2->Fit("dgaus","QLRM");
658 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
659 break;
660
661 case 3:
662 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
663 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
664 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
665 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
666 tgaus->SetBit(kCanDelete);
667 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
668 "A_{2}","#mu_{2}","#sigma_{2}",
669 "A_{3}","#mu_{3}","#sigma_{3}");
670 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
671 integ/width/3.,(max+mean)/2.,width/4.,
672 integ/width/3.,mean,width/2.);
673 // The left-sided Gauss
674 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
675 tgaus->SetParLimits(1,min+(width/10.),mean);
676 tgaus->SetParLimits(2,width/15.,width/2.);
677 // The right-sided Gauss
678 tgaus->SetParLimits(3,0.,integ);
679 tgaus->SetParLimits(4,mean,max-(width/10.));
680 tgaus->SetParLimits(5,width/15.,width/2.);
681 // The Gauss describing the outliers
682 tgaus->SetParLimits(6,0.,integ);
683 tgaus->SetParLimits(7,min,max);
684 tgaus->SetParLimits(8,width/4.,width/1.5);
685 obj2->Fit("tgaus","QLRM");
686 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
687 break;
688 case 4:
689 obj2->Fit("pol0","Q");
690 obj2->GetFunction("pol0")->SetLineColor(kYellow);
691 break;
692 case 9:
693 break;
694 default:
695 obj2->Fit("gaus","Q");
696 obj2->GetFunction("gaus")->SetLineColor(kYellow);
697 break;
698 }
699
700 gPad->Modified();
701 gPad->Update();
702
703 }
704}
Note: See TracBrowser for help on using the repository browser.