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

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