source: trunk/MagicSoft/Mars/macros/calibrate_data.C@ 3534

Last change on this file since 3534 was 3489, checked in by gaug, 21 years ago
*** empty log message ***
File size: 23.8 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): Hendrik Bartko, 03/2004 <mailto:hbartko@mppmu.mpg.de>
19! Markus Gaug, 03/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25#include "MAGIC.h"
26
27const TString defpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/";
28const TString defrout = "output_test.root";
29
30const Int_t defpedr [] = {12461};
31const Int_t defcalr [] = {12526,12527,12528,12529};
32const Int_t defdatar[] = {12517,12518};
33
34void calibrate_data(const TString inpath=defpath,
35 const Int_t psize=1, const Int_t pedruns[]=defpedr,
36 const Int_t csize=4, const Int_t calruns[]=defcalr,
37 const Int_t dsize=2, const Int_t dataruns[]=defdatar,
38 const TString resname=defrout)
39
40{
41
42 MRunIter pruns;
43 MRunIter cruns;
44 MRunIter druns;
45
46 for (Int_t i=0;i<psize;i++) {
47 cout << "Adding pedestal run: " << pedruns[i] << endl;
48 pruns.AddRun(pedruns[i],inpath);
49 }
50 for (Int_t i=0;i<csize;i++) {
51 cout << "Adding calibration run: " << calruns[i] << endl;
52 cruns.AddRun(calruns[i],inpath);
53 }
54 for (Int_t i=0;i<dsize;i++) {
55 cout << "Adding data run: " << dataruns[i] << endl;
56 druns.AddRun(dataruns[i],inpath);
57 }
58
59 MStatusDisplay *display = new MStatusDisplay;
60 display->SetUpdateTime(3000);
61 display->Resize(850,700);
62
63 gStyle->SetOptStat(1111);
64 gStyle->SetOptFit();
65
66 /************************************/
67 /* FIRST LOOP: PEDESTAL COMPUTATION */
68 /************************************/
69
70 MParList plist1;
71 MTaskList tlist1;
72 plist1.AddToList(&tlist1);
73
74 // containers
75 MPedestalCam pedcam;
76 MBadPixelsCam badcam;
77
78 //
79 // for excluding pixels from the beginning:
80 //
81 // badcam1.AsciiRead("badpixels.dat");
82
83 plist1.AddToList(&pedcam);
84 plist1.AddToList(&badcam);
85
86 //tasks
87 MReadMarsFile read("Events");
88 read.DisableAutoScheme();
89 static_cast<MRead&>(read).AddFiles(pruns);
90
91 MGeomApply geomapl;
92 MPedCalcPedRun pedcalc;
93 MGeomCamMagic geomcam;
94
95 tlist1.AddToList(&read);
96 tlist1.AddToList(&geomapl);
97 tlist1.AddToList(&pedcalc);
98
99 // Create and execute the event looper
100 MEvtLoop pedloop;
101 pedloop.SetParList(&plist1);
102 pedloop.SetDisplay(display);
103
104 cout << "*************************" << endl;
105 cout << "** COMPUTING PEDESTALS **" << endl;
106 cout << "*************************" << endl;
107
108 if (!pedloop.Eventloop())
109 return;
110
111 tlist1.PrintStatistics();
112
113 //
114 // Now the short version:
115 //
116/*
117 MJCalibration calloop;
118 calloop.SetInput(&cruns);
119 calloop.SetDisplay(display);
120 if (!calloop.Process(pedloop.GetPedestalCam()))
121 return;
122#if 0
123*/
124 //
125 // The longer version:
126 //
127
128 //
129 // Create a empty Parameter List and an empty Task List
130 //
131 MParList plist2;
132 MTaskList tlist2;
133 plist2.AddToList(&tlist2);
134 plist2.AddToList(&pedcam);
135 plist2.AddToList(&badcam);
136
137 gLog << endl;;
138 gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
139 gLog << endl;
140
141 MReadMarsFile read2("Events");
142 read2.DisableAutoScheme();
143 static_cast<MRead&>(read2).AddFiles(cruns);
144
145 MGeomCamMagic geomcam;
146 MExtractedSignalCam sigcam;
147 MArrivalTimeCam timecam;
148 MCalibrationChargeCam calcam;
149 // MCalibrationChargePINDiode pindiode;
150 // MCalibrationChargeBlindPix blindpix;
151
152 MHCalibrationRelTimeCam histtime;
153 MHCalibrationChargeCam histcharge;
154 // MHCalibrationChargePINDiode histpin;
155 // MHCalibrationChargeBlindPix histblind;
156 //
157 //
158 // Get the previously created MPedestalCam into the new Parameter List
159 //
160 plist2.AddToList(&geomcam);
161 plist2.AddToList(&sigcam);
162 plist2.AddToList(&timecam);
163 plist2.AddToList(&calcam);
164 plist2.AddToList(&histtime);
165 plist2.AddToList(&histcharge);
166// plist2.AddToList(&histpin);
167// plist2.AddToList(&histblind);
168
169 //
170 // We saw that the signal jumps between slices,
171 // thus take the sliding window
172 //
173 MExtractSignal2 sigcalc;
174 MExtractPINDiode pincalc;
175 MExtractBlindPixel blindcalc;
176
177 MArrivalTimeCalc2 timecalc;
178 MCalibrationChargeCalc calcalc;
179 MGeomApply geomapl;
180
181 MFillH filltime( "MHCalibrationRelTimeCam" , "MArrivalTimeCam");
182// MFillH fillpin ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
183// MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
184 MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam");
185
186 //
187 // Skip the HiGain vs. LoGain calibration
188 //
189 calcalc.SkipHiLoGainCalibration();
190
191 //
192 // Apply a filter against cosmics
193 // (was directly in MCalibrationCalc in earlier versions)
194 //
195 MFCosmics cosmics;
196 MContinue cont(&cosmics);
197
198 tlist2.AddToList(&read2);
199 tlist2.AddToList(&geomapl);
200 tlist2.AddToList(&sigcalc);
201 // tlist2.AddToList(&blindcalc);
202 // tlist2.AddToList(&pincalc);
203 //
204 // In case, you want to skip the cosmics rejection,
205 // uncomment the next line
206 //
207 tlist2.AddToList(&cont);
208 //
209 // In case, you want to skip the somewhat lengthy calculation
210 // of the arrival times using a spline, uncomment the next two lines
211 //
212 tlist2.AddToList(&timecalc);
213 tlist2.AddToList(&filltime);
214// tlist2.AddToList(&fillpin);
215// tlist2.AddToList(&fillblind);
216 tlist2.AddToList(&fillcam);
217 //
218 tlist2.AddToList(&calcalc);
219 //
220 // Create and setup the eventloop
221 //
222 MEvtLoop evtloop;
223 evtloop.SetParList(&plist2);
224 evtloop.SetDisplay(display);
225
226 cout << "***************************" << endl;
227 cout << "** COMPUTING CALIBRATION **" << endl;
228 cout << "***************************" << endl;
229
230 //
231 // Execute second analysis
232 //
233 if (!evtloop.Eventloop())
234 return;
235
236 tlist2.PrintStatistics();
237
238 //
239 // print the most important results of all pixels to a file
240 //
241 // MLog gauglog;
242 // gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
243 // calcam.SetLogStream(&gauglog);
244 // calcam.Print();
245 // calcam.SetLogStream(&gLog);
246 //
247 // just one example how to get the plots of individual pixels
248 //
249 // histblind.DrawClone("all");
250 // histcharge[5].DrawClone("all");
251 // histcharge(5).DrawClone("all");
252 // histtime[5].DrawClone("fourierevents");
253
254 // Create histograms to display
255 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
256 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
257 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
258 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
259 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
260 MHCamera disp6 (geomcam, "Cal;FFactorPh", "Nr. of Photo-electrons (F-Factor Method)");
261 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (F-Factor Method)");
262 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
263 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)");
264 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)");
265 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
266 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)");
267 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor tp photons (PIN Diode Method)");
268 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
269 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
270 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
271 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
272 MHCamera disp18 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels HI Gain");
273 MHCamera disp19 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO Gain");
274 MHCamera disp20 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain");
275 MHCamera disp21 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain");
276 MHCamera disp22 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
277 MHCamera disp23 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration");
278 MHCamera disp24 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration");
279 MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration");
280
281 MHCamera disp26 (geomcam, "Cal;Ped", "Pedestals");
282 MHCamera disp27 (geomcam, "Cal;PedRms", "Pedestal RMS");
283
284 MHCamera disp28 (geomcam, "time;Time", "Rel. Arrival Times");
285 MHCamera disp29 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times");
286 MHCamera disp30 (geomcam, "time;TimeProb", "Probability of Time Fit");
287 MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
288 MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels");
289
290 MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
291 MHCamera disp34 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times");
292
293 // Fitted charge means and sigmas
294 disp1.SetCamContent(calcam, 0);
295 disp1.SetCamError( calcam, 1);
296 disp2.SetCamContent(calcam, 2);
297 disp2.SetCamError( calcam, 3);
298
299 // Fit probabilities
300 disp3.SetCamContent(calcam, 4);
301
302 // Reduced Sigmas and reduced sigmas per charge
303 disp4.SetCamContent(calcam, 5);
304 disp4.SetCamError( calcam, 6);
305 disp5.SetCamContent(calcam, 7);
306 disp5.SetCamError( calcam, 8);
307
308 // F-Factor Method
309 disp6.SetCamContent(calcam, 9);
310 disp6.SetCamError( calcam, 10);
311 disp7.SetCamContent(calcam, 11);
312 disp7.SetCamError( calcam, 12);
313 disp8.SetCamContent(calcam, 13);
314 disp8.SetCamError( calcam, 14);
315
316 /// Blind Pixel Method
317 disp9.SetCamContent(calcam, 15);
318 disp9.SetCamError( calcam, 16);
319 disp10.SetCamContent(calcam,17);
320 disp10.SetCamError( calcam,18);
321 disp11.SetCamContent(calcam,19);
322 disp11.SetCamError( calcam,20);
323
324 // PIN Diode Method
325 disp12.SetCamContent(calcam,21);
326 disp12.SetCamError( calcam,22);
327 disp13.SetCamContent(calcam,23);
328 disp13.SetCamError( calcam,24);
329 disp14.SetCamContent(calcam,25);
330 disp14.SetCamError( calcam,26);
331
332 // Pixels with defects
333 disp15.SetCamContent(calcam,27);
334 disp16.SetCamContent(calcam,28);
335 disp17.SetCamContent(badcam,9);
336 disp18.SetCamContent(badcam,16);
337 disp19.SetCamContent(badcam,15);
338 disp20.SetCamContent(calcam,29);
339 disp21.SetCamContent(calcam,30);
340
341 // Lo Gain calibration
342 disp22.SetCamContent(calcam,31);
343
344 // Valid flags
345 disp23.SetCamContent(calcam,32);
346 disp24.SetCamContent(calcam,33);
347 disp25.SetCamContent(calcam,34);
348
349 // Pedestals
350 disp26.SetCamContent(calcam,35);
351 disp26.SetCamError( calcam,36);
352 disp27.SetCamContent(calcam,37);
353 disp27.SetCamError( calcam,38);
354
355
356 // Relative Times
357 disp28.SetCamContent(histtime,0);
358 disp28.SetCamError( histtime,1);
359 disp29.SetCamContent(histtime,2);
360 disp29.SetCamError( histtime,3);
361 disp30.SetCamContent(histtime,4);
362 disp31.SetCamContent(histtime,5);
363 disp32.SetCamContent(histtime,6);
364
365 // Absolute Times
366 disp33.SetCamContent(calcam,39);
367 disp33.SetCamError( calcam,40);
368 disp34.SetCamContent(calcam,41);
369
370 disp1.SetYTitle("Charge [FADC units]");
371 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
372 disp3.SetYTitle("P_{Charge} [1]");
373
374 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
375 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
376
377 disp6.SetYTitle("Nr. Photo-electrons [1]");
378 disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
379 disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");
380
381 disp9.SetYTitle("Photon flux [ph/mm^2]");
382 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
383 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
384
385 disp12.SetYTitle("Photon flux [ph/mm^2]");
386 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
387 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
388
389 disp15.SetYTitle("[1]");
390 disp16.SetYTitle("[1]");
391 disp17.SetYTitle("[1]");
392 disp18.SetYTitle("[1]");
393 disp19.SetYTitle("[1]");
394 disp20.SetYTitle("[1]");
395 disp21.SetYTitle("[1]");
396 disp22.SetYTitle("[1]");
397 disp23.SetYTitle("[1]");
398 disp24.SetYTitle("[1]");
399 disp25.SetYTitle("[1]");
400
401 disp26.SetYTitle("Ped [FADC Counts ]");
402 disp27.SetYTitle("RMS_{Ped} [FADC Counts ]");
403
404 disp28.SetYTitle("Time Offset [ns]");
405 disp29.SetYTitle("Timing resolution [ns]");
406 disp30.SetYTitle("P_{Time} [1]");
407
408 disp31.SetYTitle("[1]");
409 disp32.SetYTitle("[1]");
410
411 disp33.SetYTitle("Mean Abs. Time [FADC slice]");
412 disp34.SetYTitle("RMS Abs. Time [FADC slices]");
413
414 gStyle->SetOptStat(1111);
415 gStyle->SetOptFit();
416
417 // Charges
418 TCanvas &c1 = display->AddTab("Fit.Charge");
419 c1.Divide(2, 3);
420
421 CamDraw(c1, disp1,calcam,1, 2 , 2);
422 CamDraw(c1, disp2,calcam,2, 2 , 2);
423
424 // Fit Probability
425 TCanvas &c2 = display->AddTab("Fit.Prob");
426 c2.Divide(1,3);
427
428 CamDraw(c2, disp3,calcam,1, 1 , 4);
429
430 // Reduced Sigmas
431 TCanvas &c3 = display->AddTab("Red.Sigma");
432 c3.Divide(2,3);
433
434 CamDraw(c3, disp4,calcam,1, 2 , 2);
435 CamDraw(c3, disp5,calcam,2, 2 , 2);
436
437 // F-Factor Method
438 TCanvas &c4 = display->AddTab("F-Factor");
439 c4.Divide(3,3);
440
441 CamDraw(c4, disp6,calcam,1, 3 , 2);
442 CamDraw(c4, disp7,calcam,2, 3 , 2);
443 CamDraw(c4, disp8,calcam,3, 3 , 2);
444
445 // Blind Pixel Method
446 TCanvas &c5 = display->AddTab("BlindPix");
447 c5.Divide(3, 3);
448
449 CamDraw(c5, disp9,calcam,1, 3 , 9);
450 CamDraw(c5, disp10,calcam,2, 3 , 2);
451 CamDraw(c5, disp11,calcam,3, 3 , 2);
452
453 // PIN Diode Method
454 TCanvas &c6 = display->AddTab("PINDiode");
455 c6.Divide(3,3);
456
457 CamDraw(c6, disp12,calcam,1, 3 , 9);
458 CamDraw(c6, disp13,calcam,2, 3 , 2);
459 CamDraw(c6, disp14,calcam,3, 3 , 2);
460
461 // Defects
462 TCanvas &c7 = display->AddTab("Defects");
463 c7.Divide(4,2);
464
465 CamDraw(c7, disp15,calcam,1,4, 0);
466 CamDraw(c7, disp16,calcam,2,4, 0);
467 CamDraw(c7, disp20,calcam,3,4, 0);
468 CamDraw(c7, disp21,calcam,4,4, 0);
469
470 // BadCam
471 TCanvas &c8 = display->AddTab("Defects");
472 c8.Divide(3,2);
473
474 CamDraw(c8, disp17,badcam,1,3, 0);
475 CamDraw(c8, disp18,badcam,2,3, 0);
476 CamDraw(c8, disp19,badcam,3,3, 0);
477
478
479 // Valid flags
480 TCanvas &c9 = display->AddTab("Validity");
481 c9.Divide(4,2);
482
483 CamDraw(c9, disp22,calcam,1,4,0);
484 CamDraw(c9, disp23,calcam,2,4,0);
485 CamDraw(c9, disp24,calcam,3,4,0);
486 CamDraw(c9, disp25,calcam,4,4,0);
487
488
489 // Pedestals
490 TCanvas &c10 = display->AddTab("Pedestals");
491 c10.Divide(2,3);
492
493 CamDraw(c10,disp26,calcam,1,2,1);
494 CamDraw(c10,disp27,calcam,2,2,2);
495
496 // Rel. Times
497 TCanvas &c11 = display->AddTab("Fitted Rel. Times");
498 c11.Divide(3,3);
499
500 CamDraw(c11,disp28,calcam,1,3,2);
501 CamDraw(c11,disp29,calcam,2,3,2);
502 CamDraw(c11,disp30,calcam,3,3,4);
503
504 // Time Defects
505 TCanvas &c12 = display->AddTab("Time Def.");
506 c12.Divide(2,2);
507
508 CamDraw(c12, disp31,calcam,1,2, 0);
509 CamDraw(c12, disp32,calcam,2,2, 0);
510
511 // Abs. Times
512 TCanvas &c13 = display->AddTab("Abs. Times");
513 c13.Divide(2,3);
514
515 CamDraw(c13,disp33,calcam,1,2,2);
516 CamDraw(c13,disp34,calcam,2,2,2);
517
518 /************************************************************************/
519 /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */
520 /************************************************************************/
521
522 // Create an empty Parameter List and an empty Task List
523 MParList plist3;
524 MTaskList tlist3;
525 plist3.AddToList(&tlist3);
526
527 // containers
528 MCerPhotEvt photevt;
529 MPedPhotCam pedphotcam;
530 MSrcPosCam srccam;
531 MRawRunHeader runhead;
532
533 plist3.AddToList(&geomcam );
534 plist3.AddToList(&pedcam );
535 plist3.AddToList(&calcam );
536 plist3.AddToList(&badcam );
537 plist3.AddToList(&timecam );
538 plist3.AddToList(&sigcam );
539 plist3.AddToList(&histtime);
540 plist3.AddToList(&photevt);
541 plist3.AddToList(&pedphotcam);
542 plist3.AddToList(&srccam);
543 plist3.AddToList(&runhead);
544
545 //tasks
546 MReadMarsFile read3("Events");
547 read3.DisableAutoScheme();
548 static_cast<MRead&>(read3).AddFiles(druns);
549
550 MCalibrateData photcalc;
551 photcalc.SetCalibrationMode(MCalibrateData::kFfactor); // !!! was only MCalibrate
552 // MPedPhotCalc pedphotcalc; // already done by MCalibrate Data
553 // MCerPhotCalc cerphotcalc; // already done by MCalibrate Data
554
555 tlist3.AddToList(&read3);
556 tlist3.AddToList(&geomapl);
557 tlist3.AddToList(&sigcalc);
558 tlist3.AddToList(&timecalc);
559 // tlist3.AddToList(&cerphotcalc); // already done by MCalibrate Data
560 tlist3.AddToList(&photcalc);
561 // tlist3.AddToList(&pedphotcalc); // already done by MCalibrate Data
562
563 MWriteRootFile write(resname);
564
565 write.AddContainer("MGeomCam" , "RunHeaders");
566 write.AddContainer("MRawRunHeader" , "RunHeaders");
567 write.AddContainer("MSrcPosCam" , "RunHeaders");
568 write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
569 // write.AddContainer("MPedPhotCam","RunHeaders"); // Attention, was in Events - Tree!!
570 write.AddContainer("MPedestalCam" , "RunHeaders");
571 write.AddContainer("MHCalibrationRelTimeCam","RunHeaders");
572
573 write.AddContainer("MCerPhotEvt" , "Events");
574 write.AddContainer("MRawEvtHeader" , "Events");
575 write.AddContainer("MBadPixelsCam" , "Events");
576 write.AddContainer("MPedPhotCam" , "Events");
577
578 tlist3.AddToList(&write);
579
580 // Create and execute eventloop
581 MEvtLoop evtloop3;
582 evtloop3.SetParList(&plist3);
583
584 cout << "*************************************************************" << endl;
585 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
586 cout << "*************************************************************" << endl;
587
588 if (!evtloop3.Eventloop())
589 return;
590 tlist3.PrintStatistics();
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(obj1.GetName());
612
613// obj2->Sumw2();
614 obj2->Draw();
615 obj2->SetBit(kCanDelete);
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 TArrayI s0(3);
701 s0[0] = 6;
702 s0[1] = 1;
703 s0[2] = 2;
704
705 TArrayI s1(3);
706 s1[0] = 3;
707 s1[1] = 4;
708 s1[2] = 5;
709
710 TArrayI inner(1);
711 inner[0] = 0;
712
713 TArrayI outer(1);
714 outer[0] = 1;
715
716 // Just to get the right (maximum) binning
717 TH1D *half[4];
718 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
719 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
720 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
721 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
722
723 for (int i=0; i<4; i++)
724 {
725 half[i]->SetLineColor(kRed+i);
726 half[i]->SetDirectory(0);
727 half[i]->SetBit(kCanDelete);
728 half[i]->Draw("same");
729 }
730
731 gPad->Modified();
732 gPad->Update();
733
734 }
735}
736
737
Note: See TracBrowser for help on using the repository browser.