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

Last change on this file since 3704 was 3651, checked in by gaug, 20 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-2004
21!
22!
23\* ======================================================================== */
24#include "MAGIC.h"
25
26void 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
550void 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
Note: See TracBrowser for help on using the repository browser.