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

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