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

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