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

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