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-2004
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 | #include "MAGIC.h"
|
---|
25 |
|
---|
26 | void calibration()
|
---|
27 | {
|
---|
28 |
|
---|
29 | // const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/";
|
---|
30 | const TString inpath = "/mnt/Data/rootdata/Miscellaneous/2004_03_03/";
|
---|
31 | // const TString inpath = "/home/rootdata/BlindPixel/";
|
---|
32 |
|
---|
33 | MRunIter pruns;
|
---|
34 | MRunIter cruns;
|
---|
35 |
|
---|
36 | pruns.AddRun(20132,inpath);
|
---|
37 | cruns.AddRun(20134,inpath);
|
---|
38 | // cruns.AddRun(16774,inpath);
|
---|
39 |
|
---|
40 | gStyle->SetOptStat(1111);
|
---|
41 | gStyle->SetOptFit();
|
---|
42 |
|
---|
43 | MStatusDisplay *display = new MStatusDisplay;
|
---|
44 | display->SetUpdateTime(3000);
|
---|
45 | display->Resize(850,700);
|
---|
46 |
|
---|
47 | MBadPixelsCam badcam;
|
---|
48 | badcam.AsciiRead("badpixels.dat");
|
---|
49 |
|
---|
50 | /************************************/
|
---|
51 | /* FIRST LOOP: PEDESTAL COMPUTATION */
|
---|
52 | /************************************/
|
---|
53 |
|
---|
54 | MJPedestal pedloop;
|
---|
55 | pedloop.SetInput(&pruns);
|
---|
56 | pedloop.SetDisplay(display);
|
---|
57 | pedloop.SetBadPixels(badcam);
|
---|
58 |
|
---|
59 | if (!pedloop.Process())
|
---|
60 | return;
|
---|
61 |
|
---|
62 | //
|
---|
63 | // Create a empty Parameter List and an empty Task List
|
---|
64 | // The tasklist is identified in the eventloop by its name
|
---|
65 | //
|
---|
66 | MParList plist0;
|
---|
67 | MTaskList tlist0;
|
---|
68 | plist0.AddToList(&tlist0);
|
---|
69 |
|
---|
70 | //
|
---|
71 | // Now setup the tasks and tasklist for the pedestals:
|
---|
72 | // ---------------------------------------------------
|
---|
73 | //
|
---|
74 | MReadMarsFile read("Events");
|
---|
75 | read.DisableAutoScheme();
|
---|
76 | static_cast<MRead&>(read).AddFiles(pruns);
|
---|
77 |
|
---|
78 | MGeomApply geomapl;
|
---|
79 | MExtractSignal sigcalc;
|
---|
80 |
|
---|
81 | //
|
---|
82 | // Additionally to calculating the pedestals,
|
---|
83 | // you can fill histograms and look at them
|
---|
84 | //
|
---|
85 | MFillH fill("MHPedestalCam", "MExtractedSignalCam");
|
---|
86 |
|
---|
87 | tlist0.AddToList(&read);
|
---|
88 | tlist0.AddToList(&geomapl);
|
---|
89 | tlist0.AddToList(&sigcalc);
|
---|
90 | tlist0.AddToList(&fill);
|
---|
91 |
|
---|
92 | MGeomCamMagic geomcam;
|
---|
93 | MHPedestalCam hpedcam;
|
---|
94 | plist0.AddToList(&geomcam);
|
---|
95 | plist0.AddToList(&pedloop.GetPedestalCam());
|
---|
96 | plist0.AddToList(&badcam);
|
---|
97 | plist0.AddToList(&hpedcam);
|
---|
98 |
|
---|
99 | //
|
---|
100 | // Create and setup the eventloop
|
---|
101 | //
|
---|
102 | MEvtLoop evtloop0;
|
---|
103 |
|
---|
104 | evtloop0.SetParList(&plist0);
|
---|
105 | evtloop0.SetDisplay(display);
|
---|
106 |
|
---|
107 | //
|
---|
108 | // Execute first analysis
|
---|
109 | //
|
---|
110 | if (!evtloop0.Eventloop())
|
---|
111 | return;
|
---|
112 |
|
---|
113 | tlist0.PrintStatistics();
|
---|
114 |
|
---|
115 | for (Int_t aidx=0;aidx<2;aidx++)
|
---|
116 | {
|
---|
117 | hpedcam.GetAverageHiGainArea(aidx).DrawClone("fourierevents");
|
---|
118 | hpedcam.GetAverageLoGainArea(aidx).DrawClone("fourierevents");
|
---|
119 | }
|
---|
120 |
|
---|
121 | for (Int_t sector=1;sector<7;sector++)
|
---|
122 | {
|
---|
123 | hpedcam.GetAverageHiGainSector(sector).DrawClone("fourierevents");
|
---|
124 | hpedcam.GetAverageLoGainSector(sector).DrawClone("fourierevents");
|
---|
125 | }
|
---|
126 |
|
---|
127 |
|
---|
128 | //
|
---|
129 | // Now the short version:
|
---|
130 | //
|
---|
131 | /*
|
---|
132 | MJCalibration calloop;
|
---|
133 | calloop.SetInput(&cruns);
|
---|
134 | calloop.SetDisplay(display);
|
---|
135 | calloop.SetBadPixels(pedloop.GetBadPixelsCam());
|
---|
136 | if (!calloop.Process(pedloop.GetPedestalCam()))
|
---|
137 | return;
|
---|
138 | #if 0
|
---|
139 | */
|
---|
140 | //
|
---|
141 | // The longer version:
|
---|
142 | //
|
---|
143 |
|
---|
144 | //
|
---|
145 | // Create a empty Parameter List and an empty Task List
|
---|
146 | //
|
---|
147 | MParList plist;
|
---|
148 | MTaskList tlist;
|
---|
149 | plist.AddToList(&tlist);
|
---|
150 | plist.AddToList(&pedloop.GetPedestalCam());
|
---|
151 | plist.AddToList(&pedloop.GetBadPixels());
|
---|
152 |
|
---|
153 | gLog << endl;;
|
---|
154 | gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
|
---|
155 | gLog << endl;
|
---|
156 |
|
---|
157 | MReadMarsFile read("Events");
|
---|
158 | read.DisableAutoScheme();
|
---|
159 | static_cast<MRead&>(read).AddFiles(cruns);
|
---|
160 |
|
---|
161 | MGeomCamMagic geomcam;
|
---|
162 | MExtractedSignalCam sigcam;
|
---|
163 | MArrivalTimeCam timecam;
|
---|
164 | MCalibrationRelTimeCam relcam;
|
---|
165 | MCalibrationChargeCam calcam;
|
---|
166 | MCalibrationChargePINDiode pindiode;
|
---|
167 | MCalibrationChargeBlindPix blindpix;
|
---|
168 |
|
---|
169 | MHCalibrationRelTimeCam histtime;
|
---|
170 | MHCalibrationChargeCam histcharge;
|
---|
171 | MHCalibrationChargePINDiode histpin;
|
---|
172 | MHCalibrationChargeBlindPix histblind;
|
---|
173 | histcharge.SetPulserFrequency(500);
|
---|
174 | histblind.SetSinglePheCut(600);
|
---|
175 | //
|
---|
176 | // As long, as we don't have digital modules,
|
---|
177 | // we have to set the color of the pulser LED by hand
|
---|
178 | //
|
---|
179 | blindpix.SetColor(kCT1);
|
---|
180 | // pindiode.SetColor(kCT1);
|
---|
181 | //
|
---|
182 | // Get the previously created MPedestalCam into the new Parameter List
|
---|
183 | //
|
---|
184 | plist.AddToList(&geomcam);
|
---|
185 | plist.AddToList(&sigcam);
|
---|
186 | plist.AddToList(&timecam);
|
---|
187 | plist.AddToList(&relcam);
|
---|
188 | plist.AddToList(&calcam);
|
---|
189 | plist.AddToList(&histtime);
|
---|
190 | plist.AddToList(&histcharge);
|
---|
191 | // plist.AddToList(&histpin);
|
---|
192 | plist.AddToList(&histblind);
|
---|
193 |
|
---|
194 | //
|
---|
195 | // We saw that the signal jumps between slices,
|
---|
196 | // thus take the sliding window
|
---|
197 | //
|
---|
198 | MExtractSignal2 sigcalc2;
|
---|
199 | MExtractPINDiode pincalc;
|
---|
200 | MExtractBlindPixel blindcalc;
|
---|
201 | sigcalc2.SetRange(2,15,6,5,14,6);
|
---|
202 | blindcalc.SetRange(12,17);
|
---|
203 |
|
---|
204 | MArrivalTimeCalc2 timecalc;
|
---|
205 | MCalibrationChargeCalc calcalc;
|
---|
206 | MGeomApply geomapl;
|
---|
207 |
|
---|
208 | MFillH filltime( "MHCalibrationRelTimeCam" , "MArrivalTimeCam");
|
---|
209 | // MFillH fillpin ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
|
---|
210 | MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
|
---|
211 | MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam");
|
---|
212 |
|
---|
213 | //
|
---|
214 | // Skip the HiGain vs. LoGain calibration
|
---|
215 | //
|
---|
216 | calcalc.SkipHiLoGainCalibration();
|
---|
217 |
|
---|
218 | //
|
---|
219 | // Apply a filter against cosmics
|
---|
220 | // (was directly in MCalibrationCalc in earlier versions)
|
---|
221 | //
|
---|
222 | MFCosmics cosmics;
|
---|
223 | MContinue cont(&cosmics);
|
---|
224 |
|
---|
225 | tlist.AddToList(&read);
|
---|
226 | tlist.AddToList(&geomapl);
|
---|
227 | tlist.AddToList(&sigcalc2);
|
---|
228 | tlist.AddToList(&blindcalc);
|
---|
229 | // tlist.AddToList(&pincalc);
|
---|
230 | //
|
---|
231 | // In case, you want to skip the cosmics rejection,
|
---|
232 | // uncomment the next line
|
---|
233 | //
|
---|
234 | tlist.AddToList(&cont);
|
---|
235 | //
|
---|
236 | // In case, you want to skip the somewhat lengthy calculation
|
---|
237 | // of the arrival times using a spline, uncomment the next two lines
|
---|
238 | //
|
---|
239 | tlist.AddToList(&timecalc);
|
---|
240 | tlist.AddToList(&filltime);
|
---|
241 | // tlist.AddToList(&fillpin);
|
---|
242 | tlist.AddToList(&fillblind);
|
---|
243 | tlist.AddToList(&fillcam);
|
---|
244 | //
|
---|
245 | tlist.AddToList(&calcalc);
|
---|
246 | //
|
---|
247 | // Create and setup the eventloop
|
---|
248 | //
|
---|
249 | MEvtLoop evtloop;
|
---|
250 | evtloop.SetParList(&plist);
|
---|
251 | evtloop.SetDisplay(display);
|
---|
252 |
|
---|
253 | //
|
---|
254 | // Execute second analysis
|
---|
255 | //
|
---|
256 | if (!evtloop.Eventloop())
|
---|
257 | return;
|
---|
258 |
|
---|
259 | tlist.PrintStatistics();
|
---|
260 |
|
---|
261 | MBadPixelsCam *badpixels = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
|
---|
262 |
|
---|
263 | //
|
---|
264 | // print the most important results of all pixels to a file
|
---|
265 | //
|
---|
266 | /*
|
---|
267 | MLog gauglog;
|
---|
268 | gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
|
---|
269 | calcam.SetLogStream(&gauglog);
|
---|
270 | badpixels->Print();
|
---|
271 | calcam.SetLogStream(&gLog);
|
---|
272 | */
|
---|
273 | //
|
---|
274 | // just one example how to get the plots of individual pixels
|
---|
275 | //
|
---|
276 | // histblind.DrawClone("all");
|
---|
277 | // histcharge[400].DrawClone("all");
|
---|
278 | // histcharge(5).DrawClone("all");
|
---|
279 | // histtime[5].DrawClone("fourierevents");
|
---|
280 | for (Int_t aidx=0;aidx<2;aidx++)
|
---|
281 | {
|
---|
282 | histcharge.GetAverageHiGainArea(aidx).DrawClone("all");
|
---|
283 | histcharge.GetAverageLoGainArea(aidx).DrawClone("all");
|
---|
284 | }
|
---|
285 |
|
---|
286 | for (Int_t sector=1;sector<7;sector++)
|
---|
287 | {
|
---|
288 | histcharge.GetAverageHiGainSector(sector).DrawClone("all");
|
---|
289 | histcharge.GetAverageLoGainSector(sector).DrawClone("all");
|
---|
290 | }
|
---|
291 |
|
---|
292 |
|
---|
293 | // Create histograms to display
|
---|
294 | MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
|
---|
295 | MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
|
---|
296 | MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
|
---|
297 | MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
|
---|
298 | MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
|
---|
299 | MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Photo-electrons (F-Factor Method)");
|
---|
300 | MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (F-Factor Method)");
|
---|
301 | MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
|
---|
302 | MHCamera disp9 (geomcam, "Cal;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)");
|
---|
303 | MHCamera disp10 (geomcam, "Cal;BlindPixFFactor", "Total F-Factor (Blind Pixel Method)");
|
---|
304 | MHCamera disp11 (geomcam, "Cal;PINDiodeConv", "Conversion Factor tp photons (PIN Diode Method)");
|
---|
305 | MHCamera disp12 (geomcam, "Cal;PINDiodeFFactor", "Total F-Factor (PIN Diode Method)");
|
---|
306 | MHCamera disp13 (geomcam, "Cal;Excluded", "Pixels previously excluded");
|
---|
307 | MHCamera disp14 (geomcam, "Cal;Unsuited", "Unsuited Pixels ");
|
---|
308 | MHCamera disp15 (geomcam, "Cal;Unreliable", "Unreliable Pixels");
|
---|
309 | MHCamera disp16 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels High Gain");
|
---|
310 | MHCamera disp17 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels Low Gain");
|
---|
311 | MHCamera disp18 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain");
|
---|
312 | MHCamera disp19 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain");
|
---|
313 | MHCamera disp20 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
|
---|
314 | MHCamera disp21 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration");
|
---|
315 | MHCamera disp22 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration");
|
---|
316 | MHCamera disp23 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration");
|
---|
317 |
|
---|
318 | MHCamera disp24 (geomcam, "Cal;Ped", "Pedestals");
|
---|
319 | MHCamera disp25 (geomcam, "Cal;PedRms", "Pedestal RMS");
|
---|
320 |
|
---|
321 | MHCamera disp26 (geomcam, "time;Time", "Rel. Arrival Times");
|
---|
322 | MHCamera disp27 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times");
|
---|
323 | MHCamera disp28 (geomcam, "time;TimeProb", "Probability of Time Fit");
|
---|
324 | MHCamera disp29 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
|
---|
325 | MHCamera disp30 (geomcam, "time;Oscillating", "Oscillating Pixels");
|
---|
326 |
|
---|
327 | MHCamera disp31 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
|
---|
328 | MHCamera disp32 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times");
|
---|
329 |
|
---|
330 | // Fitted charge means and sigmas
|
---|
331 | disp1.SetCamContent(calcam, 0);
|
---|
332 | disp1.SetCamError( calcam, 1);
|
---|
333 | disp2.SetCamContent(calcam, 2);
|
---|
334 | disp2.SetCamError( calcam, 3);
|
---|
335 |
|
---|
336 | // Fit probabilities
|
---|
337 | disp3.SetCamContent(calcam, 4);
|
---|
338 |
|
---|
339 | // Reduced Sigmas and reduced sigmas per charge
|
---|
340 | disp4.SetCamContent(calcam, 5);
|
---|
341 | disp4.SetCamError( calcam, 6);
|
---|
342 | disp5.SetCamContent(calcam, 7);
|
---|
343 | disp5.SetCamError( calcam, 8);
|
---|
344 |
|
---|
345 | // F-Factor Method
|
---|
346 | disp6.SetCamContent(calcam, 9);
|
---|
347 | disp6.SetCamError( calcam, 10);
|
---|
348 | disp7.SetCamContent(calcam, 11);
|
---|
349 | disp7.SetCamError( calcam, 12);
|
---|
350 | disp8.SetCamContent(calcam, 13);
|
---|
351 | disp8.SetCamError( calcam, 14);
|
---|
352 |
|
---|
353 | // Blind Pixel Method
|
---|
354 | disp9.SetCamContent(calcam, 16);
|
---|
355 | disp9.SetCamError( calcam, 17);
|
---|
356 | disp10.SetCamContent(calcam,18);
|
---|
357 | disp10.SetCamError( calcam,19);
|
---|
358 |
|
---|
359 | // PIN Diode Method
|
---|
360 | disp11.SetCamContent(calcam,21);
|
---|
361 | disp11.SetCamError( calcam,22);
|
---|
362 | disp12.SetCamContent(calcam,23);
|
---|
363 | disp12.SetCamError( calcam,24);
|
---|
364 |
|
---|
365 | // Pixels with defects
|
---|
366 | disp13.SetCamContent(calcam,26);
|
---|
367 | disp14.SetCamContent(*badpixels,1);
|
---|
368 | disp15.SetCamContent(*badpixels,3);
|
---|
369 | disp16.SetCamContent(*badpixels,10);
|
---|
370 | disp17.SetCamContent(*badpixels,11);
|
---|
371 | disp18.SetCamContent(calcam,27);
|
---|
372 | disp19.SetCamContent(calcam,28);
|
---|
373 |
|
---|
374 | // Lo Gain calibration
|
---|
375 | disp20.SetCamContent(calcam,29);
|
---|
376 |
|
---|
377 | // Valid flags
|
---|
378 | disp21.SetCamContent(calcam,15);
|
---|
379 | disp22.SetCamContent(calcam,20);
|
---|
380 | disp23.SetCamContent(calcam,25);
|
---|
381 |
|
---|
382 | // Pedestals
|
---|
383 | disp24.SetCamContent(calcam,30);
|
---|
384 | disp24.SetCamError( calcam,31);
|
---|
385 | disp25.SetCamContent(calcam,32);
|
---|
386 | disp25.SetCamError( calcam,33);
|
---|
387 |
|
---|
388 | // Relative Times
|
---|
389 | disp26.SetCamContent(histtime,0);
|
---|
390 | disp26.SetCamError( histtime,1);
|
---|
391 | disp27.SetCamContent(histtime,2);
|
---|
392 | disp27.SetCamError( histtime,3);
|
---|
393 | disp28.SetCamContent(histtime,4);
|
---|
394 | disp29.SetCamContent(histtime,5);
|
---|
395 | disp30.SetCamContent(histtime,6);
|
---|
396 |
|
---|
397 | // Absolute Times
|
---|
398 | disp31.SetCamContent(calcam,34);
|
---|
399 | disp31.SetCamError( calcam,35);
|
---|
400 | disp32.SetCamContent(calcam,35);
|
---|
401 |
|
---|
402 | disp1.SetYTitle("Mean Charge [FADC Counts]");
|
---|
403 | disp2.SetYTitle("\\sigma_{Charge} [FADC Counts]");
|
---|
404 | disp3.SetYTitle("P_{Sum} [1]");
|
---|
405 |
|
---|
406 | disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
|
---|
407 | disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
|
---|
408 |
|
---|
409 | disp6.SetYTitle("Nr. Photo-electrons [1]");
|
---|
410 | disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
|
---|
411 | disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");
|
---|
412 |
|
---|
413 | disp9.SetYTitle("Conversion Factor [Phot/FADC Count]");
|
---|
414 | disp10.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
|
---|
415 |
|
---|
416 | disp11.SetYTitle("Conversion Factor [Phot/FADC Count]");
|
---|
417 | disp12.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
|
---|
418 |
|
---|
419 | disp13.SetYTitle("[1]");
|
---|
420 | disp14.SetYTitle("[1]");
|
---|
421 | disp15.SetYTitle("[1]");
|
---|
422 | disp16.SetYTitle("[1]");
|
---|
423 | disp17.SetYTitle("[1]");
|
---|
424 | disp18.SetYTitle("[1]");
|
---|
425 | disp19.SetYTitle("[1]");
|
---|
426 | disp20.SetYTitle("[1]");
|
---|
427 | disp21.SetYTitle("[1]");
|
---|
428 | disp22.SetYTitle("[1]");
|
---|
429 | disp23.SetYTitle("[1]");
|
---|
430 |
|
---|
431 | disp24.SetYTitle("Ped [FADC Counts ]");
|
---|
432 | disp25.SetYTitle("RMS_{Ped} [FADC Counts ]");
|
---|
433 |
|
---|
434 | disp26.SetYTitle("Time Offset [ns]");
|
---|
435 | disp27.SetYTitle("Timing resolution [ns]");
|
---|
436 | disp28.SetYTitle("P_{Time} [1]");
|
---|
437 |
|
---|
438 | disp29.SetYTitle("[1]");
|
---|
439 | disp30.SetYTitle("[1]");
|
---|
440 |
|
---|
441 | disp31.SetYTitle("Mean Abs. Time [FADC slice]");
|
---|
442 | disp32.SetYTitle("RMS Abs. Time [FADC slices]");
|
---|
443 |
|
---|
444 | gStyle->SetOptStat(1111);
|
---|
445 | gStyle->SetOptFit();
|
---|
446 |
|
---|
447 | // Charges
|
---|
448 | TCanvas &c1 = display->AddTab("Fit.Charge");
|
---|
449 | c1.Divide(2, 4);
|
---|
450 |
|
---|
451 | CamDraw(c1, disp1,calcam,1, 2 , 2);
|
---|
452 | CamDraw(c1, disp2,calcam,2, 2 , 2);
|
---|
453 |
|
---|
454 | // Fit Probability
|
---|
455 | TCanvas &c2 = display->AddTab("Fit.Prob");
|
---|
456 | c2.Divide(1,4);
|
---|
457 |
|
---|
458 | CamDraw(c2, disp3,calcam,1,1,4);
|
---|
459 |
|
---|
460 | // Reduced Sigmas
|
---|
461 | TCanvas &c3 = display->AddTab("Red.Sigma");
|
---|
462 | c3.Divide(2,4);
|
---|
463 |
|
---|
464 | CamDraw(c3, disp4,calcam,1, 2 , 2);
|
---|
465 | CamDraw(c3, disp5,calcam,2, 2 , 2);
|
---|
466 |
|
---|
467 |
|
---|
468 | // F-Factor Method
|
---|
469 | TCanvas &c4 = display->AddTab("F-Factor");
|
---|
470 | c4.Divide(3,4);
|
---|
471 |
|
---|
472 | CamDraw(c4, disp6,calcam,1, 3 , 2);
|
---|
473 | CamDraw(c4, disp7,calcam,2, 3 , 2);
|
---|
474 | CamDraw(c4, disp8,calcam,3, 3 , 2);
|
---|
475 |
|
---|
476 |
|
---|
477 | // Blind Pixel Method
|
---|
478 | TCanvas &c5 = display->AddTab("BlindPix");
|
---|
479 | c5.Divide(2, 4);
|
---|
480 |
|
---|
481 | CamDraw(c5, disp9 ,calcam,1,2, 2);
|
---|
482 | CamDraw(c5, disp10,calcam,2,2, 2);
|
---|
483 |
|
---|
484 | // PIN Diode Method
|
---|
485 | TCanvas &c6 = display->AddTab("PINDiode");
|
---|
486 | c6.Divide(2,4);
|
---|
487 |
|
---|
488 | CamDraw(c6, disp11,calcam,1,2, 2);
|
---|
489 | CamDraw(c6, disp12,calcam,2,2, 2);
|
---|
490 |
|
---|
491 | // Defects
|
---|
492 | TCanvas &c7 = display->AddTab("Defects");
|
---|
493 | c7.Divide(4,2);
|
---|
494 |
|
---|
495 | CamDraw(c7, disp13,calcam,1,4, 0);
|
---|
496 | CamDraw(c7, disp14,calcam,2,4, 0);
|
---|
497 | CamDraw(c7, disp18,calcam,3,4, 0);
|
---|
498 | CamDraw(c7, disp19,calcam,4,4, 0);
|
---|
499 |
|
---|
500 | // BadCam
|
---|
501 | TCanvas &c8 = display->AddTab("Defects");
|
---|
502 | c8.Divide(3,2);
|
---|
503 |
|
---|
504 | CamDraw(c8, disp15,badcam,1,3, 0);
|
---|
505 | CamDraw(c8, disp16,badcam,2,3, 0);
|
---|
506 | CamDraw(c8, disp17,badcam,3,3, 0);
|
---|
507 |
|
---|
508 | // Valid flags
|
---|
509 | TCanvas &c9 = display->AddTab("Validity");
|
---|
510 | c9.Divide(4,2);
|
---|
511 |
|
---|
512 | CamDraw(c9, disp20,calcam,1,4,0);
|
---|
513 | CamDraw(c9, disp21,calcam,2,4,0);
|
---|
514 | CamDraw(c9, disp22,calcam,3,4,0);
|
---|
515 | CamDraw(c9, disp23,calcam,4,4,0);
|
---|
516 |
|
---|
517 | // Pedestals
|
---|
518 | TCanvas &c10 = display->AddTab("Pedestals");
|
---|
519 | c10.Divide(2,4);
|
---|
520 |
|
---|
521 | CamDraw(c10,disp24,calcam,1,2,1);
|
---|
522 | CamDraw(c10,disp25,calcam,2,2,2);
|
---|
523 |
|
---|
524 | // Rel. Times
|
---|
525 | TCanvas &c11 = display->AddTab("Fitted Rel. Times");
|
---|
526 | c11.Divide(3,4);
|
---|
527 |
|
---|
528 | CamDraw(c11,disp26,calcam,1,3,2);
|
---|
529 | CamDraw(c11,disp27,calcam,2,3,2);
|
---|
530 | CamDraw(c11,disp28,calcam,3,3,4);
|
---|
531 |
|
---|
532 | // Time Defects
|
---|
533 | TCanvas &c12 = display->AddTab("Time Def.");
|
---|
534 | c12.Divide(2,2);
|
---|
535 |
|
---|
536 | CamDraw(c12, disp29,calcam,1,2, 0);
|
---|
537 | CamDraw(c12, disp30,calcam,2,2, 0);
|
---|
538 |
|
---|
539 | // Abs. Times
|
---|
540 | TCanvas &c13 = display->AddTab("Abs. Times");
|
---|
541 | c13.Divide(2,4);
|
---|
542 |
|
---|
543 | CamDraw(c13,disp31,calcam,1,2,2);
|
---|
544 | CamDraw(c13,disp32,calcam,2,2,2);
|
---|
545 | #endif
|
---|
546 |
|
---|
547 | }
|
---|
548 |
|
---|
549 |
|
---|
550 | void CamDraw(TCanvas &c, MHCamera &cam, TObject &evt, Int_t i, Int_t j, Int_t fit)
|
---|
551 | {
|
---|
552 |
|
---|
553 | TArrayI s0(6);
|
---|
554 | s0[0] = 1;
|
---|
555 | s0[1] = 2;
|
---|
556 | s0[2] = 3;
|
---|
557 | s0[3] = 4;
|
---|
558 | s0[4] = 5;
|
---|
559 | s0[5] = 6;
|
---|
560 |
|
---|
561 | TArrayI s1(3);
|
---|
562 | s1[0] = 6;
|
---|
563 | s1[1] = 1;
|
---|
564 | s1[2] = 2;
|
---|
565 |
|
---|
566 | TArrayI s2(3);
|
---|
567 | s2[0] = 3;
|
---|
568 | s2[1] = 4;
|
---|
569 | s2[2] = 5;
|
---|
570 |
|
---|
571 | TArrayI inner(1);
|
---|
572 | inner[0] = 0;
|
---|
573 |
|
---|
574 | TArrayI outer(1);
|
---|
575 | outer[0] = 1;
|
---|
576 |
|
---|
577 | c.cd(i);
|
---|
578 | gPad->SetBorderMode(0);
|
---|
579 | gPad->SetTicks();
|
---|
580 | cam.GetXaxis()->SetLabelOffset(0.005);
|
---|
581 | cam.GetXaxis()->SetLabelSize(0.06);
|
---|
582 | cam.GetYaxis()->SetLabelOffset(0.005);
|
---|
583 | cam.GetYaxis()->SetLabelSize(0.06);
|
---|
584 | cam.GetXaxis()->SetTitleOffset(0.85);
|
---|
585 | cam.GetXaxis()->SetTitleSize(0.06);
|
---|
586 | cam.GetYaxis()->SetTitleOffset(0.7);
|
---|
587 | cam.GetYaxis()->SetTitleSize(0.06);
|
---|
588 | MHCamera *obj1 = (MHCamera*)cam.DrawCopy("hist");
|
---|
589 | obj1->SetDirectory(NULL);
|
---|
590 |
|
---|
591 |
|
---|
592 | c.cd(i+j);
|
---|
593 | // obj1->AddNotify(&evt);
|
---|
594 | obj1->SetPrettyPalette();
|
---|
595 | obj1->Draw();
|
---|
596 |
|
---|
597 | if (fit != 0)
|
---|
598 | {
|
---|
599 | c.cd(i+2*j);
|
---|
600 | gPad->SetBorderMode(0);
|
---|
601 | gPad->SetTicks();
|
---|
602 | TProfile *obj2 = obj1->RadialProfile(Form("%s%s",obj1->GetName(),"_rad"));
|
---|
603 | obj2->SetDirectory(NULL);
|
---|
604 | obj2->GetXaxis()->SetLabelOffset(0.005);
|
---|
605 | obj2->GetXaxis()->SetLabelSize(0.06);
|
---|
606 | obj2->GetYaxis()->SetLabelOffset(0.005);
|
---|
607 | obj2->GetYaxis()->SetLabelSize(0.06);
|
---|
608 | obj2->GetXaxis()->SetTitleOffset(0.85);
|
---|
609 | obj2->GetXaxis()->SetTitleSize(0.06);
|
---|
610 | obj2->GetYaxis()->SetTitleOffset(0.7);
|
---|
611 | obj2->GetYaxis()->SetTitleSize(0.06);
|
---|
612 | obj2->Draw();
|
---|
613 | obj2->SetBit(kCanDelete);
|
---|
614 |
|
---|
615 | TProfile *hprof[2];
|
---|
616 | hprof[0] = obj1->RadialProfileS(s0, inner,Form("%s%s",obj1->GetName(), "Inner"));
|
---|
617 | hprof[1] = obj1->RadialProfileS(s0, outer,Form("%s%s",obj1->GetName(), "Outer"));
|
---|
618 |
|
---|
619 |
|
---|
620 | for (Int_t k=0; k<2; k++)
|
---|
621 | {
|
---|
622 | Double_t min = cam.GetGeomCam().GetMinRadius(k);
|
---|
623 | Double_t max = cam.GetGeomCam().GetMaxRadius(k);
|
---|
624 |
|
---|
625 | hprof[k]->SetLineColor(kRed+k);
|
---|
626 | hprof[k]->SetDirectory(0);
|
---|
627 | hprof[k]->SetBit(kCanDelete);
|
---|
628 | hprof[k]->Draw("same");
|
---|
629 | hprof[k]->Fit("pol1","Q","",min,max);
|
---|
630 | hprof[k]->GetFunction("pol1")->SetLineColor(kRed+k);
|
---|
631 | hprof[k]->GetFunction("pol1")->SetLineWidth(1);
|
---|
632 | }
|
---|
633 |
|
---|
634 | gPad->Modified();
|
---|
635 | gPad->Update();
|
---|
636 |
|
---|
637 | c.cd(i+3*j);
|
---|
638 | gPad->SetBorderMode(0);
|
---|
639 | gPad->SetTicks();
|
---|
640 | TH1D *obj3 = (TH1D*)obj1->Projection(Form("%s%s",obj1->GetName(),"_py"));
|
---|
641 | obj3->SetDirectory(NULL);
|
---|
642 | // obj3->Sumw2();
|
---|
643 | obj3->GetXaxis()->SetLabelOffset(0.005);
|
---|
644 | obj3->GetXaxis()->SetLabelSize(0.06);
|
---|
645 | obj3->GetYaxis()->SetLabelOffset(0.005);
|
---|
646 | obj3->GetYaxis()->SetLabelSize(0.06);
|
---|
647 | obj3->GetXaxis()->SetTitleOffset(0.85);
|
---|
648 | obj3->GetXaxis()->SetTitleSize(0.06);
|
---|
649 | obj3->GetYaxis()->SetTitleOffset(0.7);
|
---|
650 | obj3->GetYaxis()->SetTitleSize(0.06);
|
---|
651 | obj3->Draw();
|
---|
652 | obj3->SetBit(kCanDelete);
|
---|
653 |
|
---|
654 | gPad->Modified();
|
---|
655 | gPad->Update();
|
---|
656 |
|
---|
657 | const Double_t min = obj3->GetBinCenter(obj3->GetXaxis()->GetFirst());
|
---|
658 | const Double_t max = obj3->GetBinCenter(obj3->GetXaxis()->GetLast());
|
---|
659 | const Double_t integ = obj3->Integral("width")/2.5066283;
|
---|
660 | const Double_t mean = obj3->GetMean();
|
---|
661 | const Double_t rms = obj3->GetRMS();
|
---|
662 | const Double_t width = max-min;
|
---|
663 |
|
---|
664 | if (rms == 0. || width == 0. )
|
---|
665 | return;
|
---|
666 |
|
---|
667 | switch (fit)
|
---|
668 | {
|
---|
669 | case 1:
|
---|
670 | TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
|
---|
671 | sgaus->SetBit(kCanDelete);
|
---|
672 | sgaus->SetParNames("Area","#mu","#sigma");
|
---|
673 | sgaus->SetParameters(integ/rms,mean,rms);
|
---|
674 | sgaus->SetParLimits(0,0.,integ);
|
---|
675 | sgaus->SetParLimits(1,min,max);
|
---|
676 | sgaus->SetParLimits(2,0,width/1.5);
|
---|
677 | obj3->Fit("sgaus","QLR");
|
---|
678 | obj3->GetFunction("sgaus")->SetLineColor(kYellow);
|
---|
679 | break;
|
---|
680 |
|
---|
681 | case 2:
|
---|
682 | TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
|
---|
683 | dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
|
---|
684 | TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
|
---|
685 | dgaus->SetBit(kCanDelete);
|
---|
686 | dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
|
---|
687 | dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
|
---|
688 | integ/width/2.,(max+mean)/2.,width/4.);
|
---|
689 | // The left-sided Gauss
|
---|
690 | dgaus->SetParLimits(0,integ-1.5,integ+1.5);
|
---|
691 | dgaus->SetParLimits(1,min+(width/10.),mean);
|
---|
692 | dgaus->SetParLimits(2,0,width/2.);
|
---|
693 | // The right-sided Gauss
|
---|
694 | dgaus->SetParLimits(3,0,integ);
|
---|
695 | dgaus->SetParLimits(4,mean,max-(width/10.));
|
---|
696 | dgaus->SetParLimits(5,0,width/2.);
|
---|
697 | obj3->Fit("dgaus","QLRM");
|
---|
698 | obj3->GetFunction("dgaus")->SetLineColor(kYellow);
|
---|
699 | break;
|
---|
700 |
|
---|
701 | case 3:
|
---|
702 | TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
|
---|
703 | tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
|
---|
704 | tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
|
---|
705 | TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
|
---|
706 | tgaus->SetBit(kCanDelete);
|
---|
707 | tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
|
---|
708 | "A_{2}","#mu_{2}","#sigma_{2}",
|
---|
709 | "A_{3}","#mu_{3}","#sigma_{3}");
|
---|
710 | tgaus->SetParameters(integ,(min+mean)/2,width/4.,
|
---|
711 | integ/width/3.,(max+mean)/2.,width/4.,
|
---|
712 | integ/width/3.,mean,width/2.);
|
---|
713 | // The left-sided Gauss
|
---|
714 | tgaus->SetParLimits(0,integ-1.5,integ+1.5);
|
---|
715 | tgaus->SetParLimits(1,min+(width/10.),mean);
|
---|
716 | tgaus->SetParLimits(2,width/15.,width/2.);
|
---|
717 | // The right-sided Gauss
|
---|
718 | tgaus->SetParLimits(3,0.,integ);
|
---|
719 | tgaus->SetParLimits(4,mean,max-(width/10.));
|
---|
720 | tgaus->SetParLimits(5,width/15.,width/2.);
|
---|
721 | // The Gauss describing the outliers
|
---|
722 | tgaus->SetParLimits(6,0.,integ);
|
---|
723 | tgaus->SetParLimits(7,min,max);
|
---|
724 | tgaus->SetParLimits(8,width/4.,width/1.5);
|
---|
725 | obj3->Fit("tgaus","QLRM");
|
---|
726 | obj3->GetFunction("tgaus")->SetLineColor(kYellow);
|
---|
727 | break;
|
---|
728 | case 4:
|
---|
729 | obj3->Fit("pol0","Q");
|
---|
730 | obj3->GetFunction("pol0")->SetLineColor(kYellow);
|
---|
731 | break;
|
---|
732 | case 9:
|
---|
733 | break;
|
---|
734 | default:
|
---|
735 | obj3->Fit("gaus","Q");
|
---|
736 | obj3->GetFunction("gaus")->SetLineColor(kYellow);
|
---|
737 | break;
|
---|
738 | }
|
---|
739 |
|
---|
740 |
|
---|
741 |
|
---|
742 | // Just to get the right (maximum) binning
|
---|
743 | TH1D *half[4];
|
---|
744 | half[0] = (TH1D*)obj1->ProjectionS(s1, inner, "Sector 6-1-2 Inner");
|
---|
745 | half[1] = (TH1D*)obj1->ProjectionS(s2, inner, "Sector 3-4-5 Inner");
|
---|
746 | half[2] = (TH1D*)obj1->ProjectionS(s1, outer, "Sector 6-1-2 Outer");
|
---|
747 | half[3] = (TH1D*)obj1->ProjectionS(s2, outer, "Sector 3-4-5 Outer");
|
---|
748 |
|
---|
749 | for (Int_t k=0; k<4; k++)
|
---|
750 | {
|
---|
751 | half[k]->SetLineColor(kRed+k);
|
---|
752 | half[k]->SetDirectory(0);
|
---|
753 | half[k]->SetBit(kCanDelete);
|
---|
754 | half[k]->Draw("same");
|
---|
755 | }
|
---|
756 |
|
---|
757 | gPad->Modified();
|
---|
758 | gPad->Update();
|
---|
759 |
|
---|
760 | }
|
---|
761 | }
|
---|
762 |
|
---|