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

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