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

Last change on this file since 3446 was 3436, checked in by gaug, 22 years ago
*** empty log message ***
File size: 19.9 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
26
27void calibration()
28{
29
30// const char *inpath = "/remote/home/pc2/operator/NewCalBoxTest/";
31// const char *inpath = "/remote/home/pc2/operator/Crab20040214/";
32// const char *inpath = "./";
33 //const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_02_10/";
34 //const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_01_26/";
35 const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/";
36 //const TString inpath = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/";
37
38 MRunIter pruns;
39 MRunIter cruns;
40
41 pruns.AddRun(12386,inpath);
42 cruns.AddRun(12525,inpath);
43 // pruns.AddRun(14400,inpath);
44 // cruns.AddRuns(14401,14409,inpath);
45// pruns.AddRuns(15487,15494,inpath);
46// cruns.AddRuns(15495,15508,inpath);
47// pruns.AddRun(19816,inpath); // CalboxTest
48// cruns.AddRun(19820,inpath); // 3blue
49// cruns.AddRun(19820,inpath); // 1blue
50
51 gStyle->SetOptStat(1111);
52 gStyle->SetOptFit();
53
54 MStatusDisplay *display = new MStatusDisplay;
55 display->SetUpdateTime(3000);
56 display->Resize(850,700);
57
58 MJPedestal pedloop;
59 pedloop.SetInput(&pruns);
60 pedloop.SetDisplay(display);
61
62 if (!pedloop.Process())
63 return;
64
65 //
66 // Now the short version:
67 //
68/*
69 MJCalibration calloop;
70 calloop.SetInput(&cruns);
71 calloop.SetDisplay(display);
72 if (!calloop.Process(pedloop.GetPedestalCam()))
73 return;
74#if 0
75*/
76 //
77 // The longer version:
78 //
79
80 //
81 // Create a empty Parameter List and an empty Task List
82 //
83 MParList plist;
84 MTaskList tlist;
85 plist.AddToList(&tlist);
86 plist.AddToList(&pedloop.GetPedestalCam());
87
88 gLog << endl;;
89 gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
90 gLog << endl;
91
92 MReadMarsFile read("Events");
93 read.DisableAutoScheme();
94 static_cast<MRead&>(read).AddFiles(cruns);
95
96 MGeomCamMagic geomcam;
97 MExtractedSignalCam sigcam;
98 MArrivalTimeCam timecam;
99 MCalibrationChargeCam calcam;
100 MCalibrationChargePINDiode pindiode;
101 MCalibrationChargeBlindPix blindpix;
102
103 MHCalibrationRelTimeCam histtime;
104 MHCalibrationChargeCam histcharge;
105 MHCalibrationChargePINDiode histpin;
106 MHCalibrationChargeBlindPix histblind;
107// histblind.SetSinglePheCut(80);
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 plist.AddToList(&pedloop.GetBadPixels());
126
127 //
128 // We saw that the signal jumps between slices,
129 // thus take the sliding window
130 //
131 MExtractSignal2 sigcalc2;
132 MExtractPINDiode pincalc;
133 MExtractBlindPixel blindcalc;
134// blindcalc.SetRange(11,16);
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 //
194 // print the most important results of all pixels to a file
195 //
196 MLog gauglog;
197 gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
198 calcam.SetLogStream(&gauglog);
199 calcam.Print();
200 calcam.SetLogStream(&gLog);
201 //
202 // just one example how to get the plots of individual pixels
203 //
204 // histblind.DrawClone("all");
205 // histcharge[5].DrawClone("all");
206 // histcharge(5).DrawClone("all");
207 // histtime[5].DrawClone("fourierevents");
208
209 // Create histograms to display
210 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
211 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
212 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
213 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
214 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
215 MHCamera disp6 (geomcam, "Cal;FFactorPh", "Nr. of Photo-electrons (F-Factor Method)");
216 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (F-Factor Method)");
217 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
218 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)");
219 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)");
220 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
221 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)");
222 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor tp photons (PIN Diode Method)");
223 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
224 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
225 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
226 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
227 MHCamera disp18 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels HI Gain");
228 MHCamera disp19 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO Gain");
229 MHCamera disp20 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain");
230 MHCamera disp21 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain");
231 MHCamera disp22 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
232 MHCamera disp23 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration");
233 MHCamera disp24 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration");
234 MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration");
235
236 MHCamera disp26 (geomcam, "Cal;Ped", "Pedestals");
237 MHCamera disp27 (geomcam, "Cal;PedRms", "Pedestal RMS");
238
239 MHCamera disp28 (geomcam, "time;Time", "Rel. Arrival Times");
240 MHCamera disp29 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times");
241 MHCamera disp30 (geomcam, "time;TimeProb", "Probability of Time Fit");
242 MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
243 MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels");
244
245 MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
246 MHCamera disp34 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times");
247
248 // Fitted charge means and sigmas
249 disp1.SetCamContent(calcam, 0);
250 disp1.SetCamError( calcam, 1);
251 disp2.SetCamContent(calcam, 2);
252 disp2.SetCamError( calcam, 3);
253
254 // Fit probabilities
255 disp3.SetCamContent(calcam, 4);
256
257 // Reduced Sigmas and reduced sigmas per charge
258 disp4.SetCamContent(calcam, 5);
259 disp4.SetCamError( calcam, 6);
260 disp5.SetCamContent(calcam, 7);
261 disp5.SetCamError( calcam, 8);
262
263 // F-Factor Method
264 disp6.SetCamContent(calcam, 9);
265 disp6.SetCamError( calcam, 10);
266 disp7.SetCamContent(calcam, 11);
267 disp7.SetCamError( calcam, 12);
268 disp8.SetCamContent(calcam, 13);
269 disp8.SetCamError( calcam, 14);
270
271 /// Blind Pixel Method
272 disp9.SetCamContent(calcam, 15);
273 disp9.SetCamError( calcam, 16);
274 disp10.SetCamContent(calcam,17);
275 disp10.SetCamError( calcam,18);
276 disp11.SetCamContent(calcam,19);
277 disp11.SetCamError( calcam,20);
278
279 // PIN Diode Method
280 disp12.SetCamContent(calcam,21);
281 disp12.SetCamError( calcam,22);
282 disp13.SetCamContent(calcam,23);
283 disp13.SetCamError( calcam,24);
284 disp14.SetCamContent(calcam,25);
285 disp14.SetCamError( calcam,26);
286
287 // Pixels with defects
288 disp15.SetCamContent(calcam,27);
289 disp16.SetCamContent(calcam,28);
290 disp17.SetCamContent(calcam,29);
291 disp18.SetCamContent(calcam,30);
292 disp19.SetCamContent(calcam,31);
293 disp20.SetCamContent(calcam,32);
294 disp21.SetCamContent(calcam,33);
295
296 // Lo Gain calibration
297 disp22.SetCamContent(calcam,34);
298
299 // Valid flags
300 disp23.SetCamContent(calcam,35);
301 disp24.SetCamContent(calcam,36);
302 disp25.SetCamContent(calcam,37);
303
304 // Pedestals
305 disp26.SetCamContent(calcam,38);
306 disp26.SetCamError( calcam,39);
307 disp27.SetCamContent(calcam,40);
308 disp27.SetCamError( calcam,41);
309
310 // Relative Times
311 disp28.SetCamContent(timecam,0);
312 disp28.SetCamError( timecam,1);
313 disp29.SetCamContent(timecam,2);
314 disp29.SetCamError( timecam,3);
315 disp30.SetCamContent(timecam,4);
316 disp31.SetCamContent(timecam,5);
317 disp32.SetCamContent(timecam,6);
318
319 // Absolute Times
320 disp33.SetCamContent(calcam,42);
321 disp33.SetCamError( calcam,43);
322 disp34.SetCamContent(calcam,43);
323
324 disp1.SetYTitle("Charge [FADC units]");
325 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
326 disp3.SetYTitle("P_{Charge} [1]");
327
328 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
329 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
330
331 disp6.SetYTitle("Nr. Photo-electrons [1]");
332 disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
333 disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");
334
335 disp9.SetYTitle("Photon flux [ph/mm^2]");
336 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
337 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
338
339 disp12.SetYTitle("Photon flux [ph/mm^2]");
340 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
341 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
342
343 disp15.SetYTitle("[1]");
344 disp16.SetYTitle("[1]");
345 disp17.SetYTitle("[1]");
346 disp18.SetYTitle("[1]");
347 disp19.SetYTitle("[1]");
348 disp20.SetYTitle("[1]");
349 disp21.SetYTitle("[1]");
350 disp22.SetYTitle("[1]");
351 disp23.SetYTitle("[1]");
352 disp24.SetYTitle("[1]");
353 disp25.SetYTitle("[1]");
354
355 disp26.SetYTitle("Ped [FADC Counts ]");
356 disp27.SetYTitle("RMS_{Ped} [FADC Counts ]");
357
358 disp28.SetYTitle("Time Offset [ns]");
359 disp29.SetYTitle("Timing resolution [ns]");
360 disp30.SetYTitle("P_{Time} [1]");
361
362 disp31.SetYTitle("[1]");
363 disp32.SetYTitle("[1]");
364
365 disp33.SetYTitle("Mean Abs. Time [FADC slice]");
366 disp34.SetYTitle("RMS Abs. Time [FADC slices]");
367
368 gStyle->SetOptStat(1111);
369 gStyle->SetOptFit();
370
371 // Charges
372 TCanvas &c1 = display->AddTab("Fit.Charge");
373 c1.Divide(2, 3);
374
375 CamDraw(c1, disp1,calcam,1, 2 , 2);
376 CamDraw(c1, disp2,calcam,2, 2 , 2);
377
378 // Fit Probability
379 TCanvas &c2 = display->AddTab("Fit.Prob");
380 c2.Divide(1,3);
381
382 CamDraw(c2, disp3,calcam,1, 1 , 4);
383
384 // Reduced Sigmas
385 TCanvas &c3 = display->AddTab("Red.Sigma");
386 c3.Divide(2,3);
387
388 CamDraw(c3, disp4,calcam,1, 2 , 2);
389 CamDraw(c3, disp5,calcam,2, 2 , 2);
390
391 // F-Factor Method
392 TCanvas &c4 = display->AddTab("F-Factor");
393 c4.Divide(3,3);
394
395 CamDraw(c4, disp6,calcam,1, 3 , 2);
396 CamDraw(c4, disp7,calcam,2, 3 , 2);
397 CamDraw(c4, disp8,calcam,3, 3 , 2);
398
399 // Blind Pixel Method
400 TCanvas &c5 = display->AddTab("BlindPix");
401 c5.Divide(3, 3);
402
403 CamDraw(c5, disp9,calcam,1, 3 , 9);
404 CamDraw(c5, disp10,calcam,2, 3 , 2);
405 CamDraw(c5, disp11,calcam,3, 3 , 2);
406
407 // PIN Diode Method
408 TCanvas &c6 = display->AddTab("PINDiode");
409 c6.Divide(3,3);
410
411 CamDraw(c6, disp12,calcam,1, 3 , 9);
412 CamDraw(c6, disp13,calcam,2, 3 , 2);
413 CamDraw(c6, disp14,calcam,3, 3 , 2);
414
415 // Defects
416 TCanvas &c7 = display->AddTab("Defects");
417 c7.Divide(7,2);
418
419 CamDraw(c7, disp15,calcam,1,7, 0);
420 CamDraw(c7, disp16,calcam,2,7, 0);
421 CamDraw(c7, disp17,calcam,3,7, 0);
422 CamDraw(c7, disp18,calcam,4,7, 0);
423 CamDraw(c7, disp19,calcam,5,7, 0);
424 CamDraw(c7, disp20,calcam,6,7, 0);
425 CamDraw(c7, disp21,calcam,7,7, 0);
426
427 // Valid flags
428 TCanvas &c8 = display->AddTab("Validity");
429 c8.Divide(4,2);
430
431 CamDraw(c8, disp22,calcam,1,4,0);
432 CamDraw(c8, disp23,calcam,2,4,0);
433 CamDraw(c8, disp24,calcam,3,4,0);
434 CamDraw(c8, disp25,calcam,4,4,0);
435
436
437 // Pedestals
438 TCanvas &c9 = display->AddTab("Pedestals");
439 c9.Divide(2,3);
440
441 CamDraw(c9,disp26,calcam,1,2,1);
442 CamDraw(c9,disp27,calcam,2,2,2);
443
444 // Rel. Times
445 TCanvas &c10 = display->AddTab("Fitted Rel. Times");
446 c10.Divide(3,3);
447
448 CamDraw(c10,disp28,calcam,1,3,2);
449 CamDraw(c10,disp29,calcam,2,3,2);
450 CamDraw(c10,disp30,calcam,3,3,4);
451
452 // Time Defects
453 TCanvas &c11 = display->AddTab("Time Def.");
454 c11.Divide(2,2);
455
456 CamDraw(c11, disp31,calcam,1,2, 0);
457 CamDraw(c11, disp32,calcam,2,2, 0);
458
459 // Abs. Times
460 TCanvas &c12 = display->AddTab("Abs. Times");
461 c12.Divide(2,3);
462
463 CamDraw(c12,disp33,calcam,1,2,2);
464 CamDraw(c12,disp34,calcam,2,2,2);
465
466#endif
467}
468
469void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
470{
471
472 c.cd(i);
473 gPad->SetBorderMode(0);
474 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
475 // obj1->AddNotify(evt);
476
477 c.cd(i+j);
478 gPad->SetBorderMode(0);
479 obj1->Draw();
480 ((MHCamera*)obj1)->SetPrettyPalette();
481
482 if (fit != 0)
483 {
484 c.cd(i+2*j);
485 gPad->SetBorderMode(0);
486 TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
487
488// obj2->Sumw2();
489 obj2->Draw();
490 obj2->SetBit(kCanDelete);
491
492 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
493 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
494 const Double_t integ = obj2->Integral("width")/2.5066283;
495 const Double_t mean = obj2->GetMean();
496 const Double_t rms = obj2->GetRMS();
497 const Double_t width = max-min;
498
499 if (rms == 0. || width == 0. )
500 return;
501
502 switch (fit)
503 {
504 case 1:
505 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
506 sgaus->SetBit(kCanDelete);
507 sgaus->SetParNames("Area","#mu","#sigma");
508 sgaus->SetParameters(integ/rms,mean,rms);
509 sgaus->SetParLimits(0,0.,integ);
510 sgaus->SetParLimits(1,min,max);
511 sgaus->SetParLimits(2,0,width/1.5);
512 obj2->Fit("sgaus","QLR");
513 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
514 break;
515
516 case 2:
517 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
518 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
519 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
520 dgaus->SetBit(kCanDelete);
521 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
522 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
523 integ/width/2.,(max+mean)/2.,width/4.);
524 // The left-sided Gauss
525 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
526 dgaus->SetParLimits(1,min+(width/10.),mean);
527 dgaus->SetParLimits(2,0,width/2.);
528 // The right-sided Gauss
529 dgaus->SetParLimits(3,0,integ);
530 dgaus->SetParLimits(4,mean,max-(width/10.));
531 dgaus->SetParLimits(5,0,width/2.);
532 obj2->Fit("dgaus","QLRM");
533 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
534 break;
535
536 case 3:
537 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
538 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
539 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
540 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
541 tgaus->SetBit(kCanDelete);
542 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
543 "A_{2}","#mu_{2}","#sigma_{2}",
544 "A_{3}","#mu_{3}","#sigma_{3}");
545 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
546 integ/width/3.,(max+mean)/2.,width/4.,
547 integ/width/3.,mean,width/2.);
548 // The left-sided Gauss
549 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
550 tgaus->SetParLimits(1,min+(width/10.),mean);
551 tgaus->SetParLimits(2,width/15.,width/2.);
552 // The right-sided Gauss
553 tgaus->SetParLimits(3,0.,integ);
554 tgaus->SetParLimits(4,mean,max-(width/10.));
555 tgaus->SetParLimits(5,width/15.,width/2.);
556 // The Gauss describing the outliers
557 tgaus->SetParLimits(6,0.,integ);
558 tgaus->SetParLimits(7,min,max);
559 tgaus->SetParLimits(8,width/4.,width/1.5);
560 obj2->Fit("tgaus","QLRM");
561 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
562 break;
563 case 4:
564 obj2->Fit("pol0","Q");
565 obj2->GetFunction("pol0")->SetLineColor(kYellow);
566 break;
567 case 9:
568 break;
569 default:
570 obj2->Fit("gaus","Q");
571 obj2->GetFunction("gaus")->SetLineColor(kYellow);
572 break;
573 }
574
575 TArrayI s0(3);
576 s0[0] = 6;
577 s0[1] = 1;
578 s0[2] = 2;
579
580 TArrayI s1(3);
581 s1[0] = 3;
582 s1[1] = 4;
583 s1[2] = 5;
584
585 TArrayI inner(1);
586 inner[0] = 0;
587
588 // Just to get the right (maximum) binning
589 TH1D *half[2];
590 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2");
591 half[1] = obj1->ProjectionS(s0, inner, "Sector 3-4-5");
592
593 for (int i=0; i<2; i++)
594 {
595 half[i]->SetLineColor(kRed+i);
596 half[i]->SetDirectory(0);
597 half[i]->SetBit(kCanDelete);
598 half[i]->Draw("same");
599 }
600
601 gPad->Modified();
602 gPad->Update();
603
604 }
605}
606
607
Note: See TracBrowser for help on using the repository browser.