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

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