source: trunk/MagicSoft/Mars/macros/calibrate_data.C@ 3749

Last change on this file since 3749 was 3567, checked in by gaug, 22 years ago
*** empty log message ***
File size: 23.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): Hendrik Bartko, 03/2004 <mailto:hbartko@mppmu.mpg.de>
19! Markus Gaug, 03/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25#include "MAGIC.h"
26
27const TString defpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/";
28const TString defrout = "output_test.root";
29
30const Int_t defpedr [] = {12461};
31const Int_t defcalr [] = {12526,12527,12528,12529};
32const Int_t defdatar[] = {12517,12518};
33
34void calibrate_data(const TString inpath=defpath,
35 const Int_t psize=1, const Int_t pedruns[]=defpedr,
36 const Int_t csize=4, const Int_t calruns[]=defcalr,
37 const Int_t dsize=2, const Int_t dataruns[]=defdatar,
38 const TString resname=defrout)
39
40{
41
42 MRunIter pruns;
43 MRunIter cruns;
44 MRunIter druns;
45
46 for (Int_t i=0;i<psize;i++) {
47 cout << "Adding pedestal run: " << pedruns[i] << endl;
48 pruns.AddRun(pedruns[i],inpath);
49 }
50 for (Int_t i=0;i<csize;i++) {
51 cout << "Adding calibration run: " << calruns[i] << endl;
52 cruns.AddRun(calruns[i],inpath);
53 }
54 for (Int_t i=0;i<dsize;i++) {
55 cout << "Adding data run: " << dataruns[i] << endl;
56 druns.AddRun(dataruns[i],inpath);
57 }
58
59 MStatusDisplay *display = new MStatusDisplay;
60 display->SetUpdateTime(3000);
61 display->Resize(850,700);
62
63 gStyle->SetOptStat(1111);
64 gStyle->SetOptFit();
65
66 /************************************/
67 /* FIRST LOOP: PEDESTAL COMPUTATION */
68 /************************************/
69
70 MParList plist1;
71 MTaskList tlist1;
72 plist1.AddToList(&tlist1);
73
74 // containers
75 MPedestalCam pedcam;
76 MBadPixelsCam badcam;
77
78 //
79 // for excluding pixels from the beginning:
80 //
81 // badcam.AsciiRead("badpixels.dat");
82
83 plist1.AddToList(&pedcam);
84 plist1.AddToList(&badcam);
85
86 //tasks
87 MReadMarsFile read("Events");
88 read.DisableAutoScheme();
89 static_cast<MRead&>(read).AddFiles(pruns);
90
91 MGeomApply geomapl;
92 MPedCalcPedRun pedcalc;
93 MGeomCamMagic geomcam;
94
95 tlist1.AddToList(&read);
96 tlist1.AddToList(&geomapl);
97 tlist1.AddToList(&pedcalc);
98
99 // Create and execute the event looper
100 MEvtLoop pedloop;
101 pedloop.SetParList(&plist1);
102 pedloop.SetDisplay(display);
103
104 cout << "*************************" << endl;
105 cout << "** COMPUTING PEDESTALS **" << endl;
106 cout << "*************************" << endl;
107
108 if (!pedloop.Eventloop())
109 return;
110
111 tlist1.PrintStatistics();
112
113 //
114 // Now the short version:
115 //
116/*
117 MJCalibration calloop;
118 calloop.SetInput(&cruns);
119 calloop.SetDisplay(display);
120 if (!calloop.Process(pedloop.GetPedestalCam()))
121 return;
122#if 0
123*/
124 //
125 // The longer version:
126 //
127
128 //
129 // Create a empty Parameter List and an empty Task List
130 //
131 MParList plist2;
132 MTaskList tlist2;
133 plist2.AddToList(&tlist2);
134 plist2.AddToList(&pedcam);
135 plist2.AddToList(&badcam);
136
137 gLog << endl;;
138 gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
139 gLog << endl;
140
141 MReadMarsFile read2("Events");
142 read2.DisableAutoScheme();
143 static_cast<MRead&>(read2).AddFiles(cruns);
144
145 MGeomCamMagic geomcam;
146 MExtractedSignalCam sigcam;
147 MArrivalTimeCam timecam;
148 MCalibrationChargeCam calcam;
149 // MCalibrationChargePINDiode pindiode;
150 // MCalibrationChargeBlindPix blindpix;
151
152 MHCalibrationRelTimeCam histtime;
153 MHCalibrationChargeCam histcharge;
154 // MHCalibrationChargePINDiode histpin;
155 // MHCalibrationChargeBlindPix histblind;
156 //
157 //
158 // Get the previously created MPedestalCam into the new Parameter List
159 //
160 plist2.AddToList(&geomcam);
161 plist2.AddToList(&sigcam);
162 plist2.AddToList(&timecam);
163 plist2.AddToList(&calcam);
164 plist2.AddToList(&histtime);
165 plist2.AddToList(&histcharge);
166// plist2.AddToList(&histpin);
167// plist2.AddToList(&histblind);
168
169 //
170 // We saw that the signal jumps between slices,
171 // thus take the sliding window
172 //
173 MExtractSignal2 sigcalc;
174 MExtractPINDiode pincalc;
175 MExtractBlindPixel blindcalc;
176
177 MArrivalTimeCalc2 timecalc;
178 MCalibrationChargeCalc calcalc;
179 MGeomApply geomapl;
180
181 MFillH filltime( "MHCalibrationRelTimeCam" , "MArrivalTimeCam");
182// MFillH fillpin ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
183// MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
184 MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam");
185
186 //
187 // Skip the HiGain vs. LoGain calibration
188 //
189 calcalc.SkipHiLoGainCalibration();
190
191 //
192 // Apply a filter against cosmics
193 // (was directly in MCalibrationCalc in earlier versions)
194 //
195 MFCosmics cosmics;
196 MContinue cont(&cosmics);
197
198 tlist2.AddToList(&read2);
199 tlist2.AddToList(&geomapl);
200 tlist2.AddToList(&sigcalc);
201 // tlist2.AddToList(&blindcalc);
202 // tlist2.AddToList(&pincalc);
203 //
204 // In case, you want to skip the cosmics rejection,
205 // uncomment the next line
206 //
207 tlist2.AddToList(&cont);
208 //
209 // In case, you want to skip the somewhat lengthy calculation
210 // of the arrival times using a spline, uncomment the next two lines
211 //
212 tlist2.AddToList(&timecalc);
213 tlist2.AddToList(&filltime);
214// tlist2.AddToList(&fillpin);
215// tlist2.AddToList(&fillblind);
216 tlist2.AddToList(&fillcam);
217 //
218 tlist2.AddToList(&calcalc);
219 //
220 // Create and setup the eventloop
221 //
222 MEvtLoop evtloop;
223 evtloop.SetParList(&plist2);
224 evtloop.SetDisplay(display);
225
226 cout << "***************************" << endl;
227 cout << "** COMPUTING CALIBRATION **" << endl;
228 cout << "***************************" << endl;
229
230 //
231 // Execute second analysis
232 //
233 if (!evtloop.Eventloop())
234 return;
235
236 tlist2.PrintStatistics();
237
238 //
239 // print the most important results of all pixels to a file
240 //
241 // MLog gauglog;
242 // gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
243 // calcam.SetLogStream(&gauglog);
244 // calcam.Print();
245 // calcam.SetLogStream(&gLog);
246 //
247 // just one example how to get the plots of individual pixels
248 //
249 // histblind.DrawClone("all");
250 // histcharge[5].DrawClone("all");
251 // histcharge(5).DrawClone("all");
252 // histtime[5].DrawClone("fourierevents");
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;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)");
264 MHCamera disp10 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel 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;NotFitted", "Pixels that could not be fitted");
269 MHCamera disp15 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
270 MHCamera disp16 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels HI Gain");
271 MHCamera disp17 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO 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 // Blind Pixel Method
315 disp9.SetCamContent(calcam, 16);
316 disp9.SetCamError( calcam, 17);
317 disp10.SetCamContent(calcam,18);
318 disp10.SetCamError( calcam,19);
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(badcam,7);
329 disp15.SetCamContent(badcam,9);
330 disp16.SetCamContent(badcam,16);
331 disp17.SetCamContent(badcam,15);
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("Charge [FADC units]");
364 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
365 disp3.SetYTitle("P_{Charge} [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("Conversion Factor [Phot/FADC Count]");
375 disp10.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [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, 3);
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,3);
418
419 CamDraw(c2, disp3,calcam,1, 1 , 4);
420
421 // Reduced Sigmas
422 TCanvas &c3 = display->AddTab("Red.Sigma");
423 c3.Divide(2,3);
424
425 CamDraw(c3, disp4,calcam,1, 2 , 2);
426 CamDraw(c3, disp5,calcam,2, 2 , 2);
427
428 // F-Factor Method
429 TCanvas &c4 = display->AddTab("F-Factor");
430 c4.Divide(3,3);
431
432 CamDraw(c4, disp6,calcam,1, 3 , 2);
433 CamDraw(c4, disp7,calcam,2, 3 , 2);
434 CamDraw(c4, disp8,calcam,3, 3 , 2);
435
436 // Blind Pixel Method
437 TCanvas &c5 = display->AddTab("BlindPix");
438 c5.Divide(2,3);
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,3);
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
468 // Valid flags
469 TCanvas &c9 = display->AddTab("Validity");
470 c9.Divide(4,2);
471
472 CamDraw(c9, disp20,calcam,1,4,0);
473 CamDraw(c9, disp21,calcam,2,4,0);
474 CamDraw(c9, disp22,calcam,3,4,0);
475 CamDraw(c9, disp23,calcam,4,4,0);
476
477
478 // Pedestals
479 TCanvas &c10 = display->AddTab("Pedestals");
480 c10.Divide(2,3);
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,3);
488
489 CamDraw(c11,disp26,calcam,1,3,2);
490 CamDraw(c11,disp27,calcam,2,3,2);
491 CamDraw(c11,disp38,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,3);
503
504 CamDraw(c13,disp31,calcam,1,2,2);
505 CamDraw(c13,disp32,calcam,2,2,2);
506
507 /************************************************************************/
508 /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */
509 /************************************************************************/
510
511 // Create an empty Parameter List and an empty Task List
512 MParList plist3;
513 MTaskList tlist3;
514 plist3.AddToList(&tlist3);
515
516 // containers
517 MCerPhotEvt photevt;
518 MPedPhotCam pedphotcam;
519 MSrcPosCam srccam;
520 MRawRunHeader runhead;
521
522 plist3.AddToList(&geomcam );
523 plist3.AddToList(&pedcam );
524 plist3.AddToList(&calcam );
525 plist3.AddToList(&badcam );
526 plist3.AddToList(&timecam );
527 plist3.AddToList(&sigcam );
528 plist3.AddToList(&histtime);
529 plist3.AddToList(&photevt);
530 plist3.AddToList(&pedphotcam);
531 plist3.AddToList(&srccam);
532 plist3.AddToList(&runhead);
533
534 //tasks
535 MReadMarsFile read3("Events");
536 read3.DisableAutoScheme();
537 static_cast<MRead&>(read3).AddFiles(druns);
538
539 MCalibrateData photcalc;
540 photcalc.SetCalibrationMode(MCalibrateData::kFfactor); // !!! was only MCalibrate
541 // MPedPhotCalc pedphotcalc; // already done by MCalibrate Data
542 // MCerPhotCalc cerphotcalc; // already done by MCalibrate Data
543
544 tlist3.AddToList(&read3);
545 tlist3.AddToList(&geomapl);
546 tlist3.AddToList(&sigcalc);
547 tlist3.AddToList(&timecalc);
548 // tlist3.AddToList(&cerphotcalc); // already done by MCalibrate Data
549 tlist3.AddToList(&photcalc);
550 // tlist3.AddToList(&pedphotcalc); // already done by MCalibrate Data
551
552 MWriteRootFile write(resname);
553
554 write.AddContainer("MGeomCam" , "RunHeaders");
555 write.AddContainer("MRawRunHeader" , "RunHeaders");
556 write.AddContainer("MSrcPosCam" , "RunHeaders");
557 write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
558 // write.AddContainer("MPedPhotCam","RunHeaders"); // Attention, was in Events - Tree!!
559 write.AddContainer("MPedestalCam" , "RunHeaders");
560 write.AddContainer("MHCalibrationRelTimeCam","RunHeaders");
561
562 write.AddContainer("MCerPhotEvt" , "Events");
563 write.AddContainer("MRawEvtHeader" , "Events");
564 write.AddContainer("MBadPixelsCam" , "Events");
565 write.AddContainer("MPedPhotCam" , "Events");
566
567 tlist3.AddToList(&write);
568
569 // Create and execute eventloop
570 MEvtLoop evtloop3;
571 evtloop3.SetParList(&plist3);
572
573 cout << "*************************************************************" << endl;
574 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
575 cout << "*************************************************************" << endl;
576
577 if (!evtloop3.Eventloop())
578 return;
579 tlist3.PrintStatistics();
580
581}
582
583void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
584{
585
586 c.cd(i);
587 gPad->SetBorderMode(0);
588 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
589 // obj1->AddNotify(evt);
590
591 c.cd(i+j);
592 gPad->SetBorderMode(0);
593 obj1->Draw();
594 ((MHCamera*)obj1)->SetPrettyPalette();
595
596 if (fit != 0)
597 {
598 c.cd(i+2*j);
599 gPad->SetBorderMode(0);
600 TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
601
602// obj2->Sumw2();
603 obj2->Draw();
604 obj2->SetBit(kCanDelete);
605
606 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
607 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
608 const Double_t integ = obj2->Integral("width")/2.5066283;
609 const Double_t mean = obj2->GetMean();
610 const Double_t rms = obj2->GetRMS();
611 const Double_t width = max-min;
612
613 if (rms == 0. || width == 0. )
614 return;
615
616 switch (fit)
617 {
618 case 1:
619 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
620 sgaus->SetBit(kCanDelete);
621 sgaus->SetParNames("Area","#mu","#sigma");
622 sgaus->SetParameters(integ/rms,mean,rms);
623 sgaus->SetParLimits(0,0.,integ);
624 sgaus->SetParLimits(1,min,max);
625 sgaus->SetParLimits(2,0,width/1.5);
626 obj2->Fit("sgaus","QLR");
627 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
628 break;
629
630 case 2:
631 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
632 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
633 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
634 dgaus->SetBit(kCanDelete);
635 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
636 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
637 integ/width/2.,(max+mean)/2.,width/4.);
638 // The left-sided Gauss
639 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
640 dgaus->SetParLimits(1,min+(width/10.),mean);
641 dgaus->SetParLimits(2,0,width/2.);
642 // The right-sided Gauss
643 dgaus->SetParLimits(3,0,integ);
644 dgaus->SetParLimits(4,mean,max-(width/10.));
645 dgaus->SetParLimits(5,0,width/2.);
646 obj2->Fit("dgaus","QLRM");
647 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
648 break;
649
650 case 3:
651 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
652 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
653 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
654 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
655 tgaus->SetBit(kCanDelete);
656 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
657 "A_{2}","#mu_{2}","#sigma_{2}",
658 "A_{3}","#mu_{3}","#sigma_{3}");
659 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
660 integ/width/3.,(max+mean)/2.,width/4.,
661 integ/width/3.,mean,width/2.);
662 // The left-sided Gauss
663 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
664 tgaus->SetParLimits(1,min+(width/10.),mean);
665 tgaus->SetParLimits(2,width/15.,width/2.);
666 // The right-sided Gauss
667 tgaus->SetParLimits(3,0.,integ);
668 tgaus->SetParLimits(4,mean,max-(width/10.));
669 tgaus->SetParLimits(5,width/15.,width/2.);
670 // The Gauss describing the outliers
671 tgaus->SetParLimits(6,0.,integ);
672 tgaus->SetParLimits(7,min,max);
673 tgaus->SetParLimits(8,width/4.,width/1.5);
674 obj2->Fit("tgaus","QLRM");
675 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
676 break;
677 case 4:
678 obj2->Fit("pol0","Q");
679 obj2->GetFunction("pol0")->SetLineColor(kYellow);
680 break;
681 case 9:
682 break;
683 default:
684 obj2->Fit("gaus","Q");
685 obj2->GetFunction("gaus")->SetLineColor(kYellow);
686 break;
687 }
688
689 TArrayI s0(3);
690 s0[0] = 6;
691 s0[1] = 1;
692 s0[2] = 2;
693
694 TArrayI s1(3);
695 s1[0] = 3;
696 s1[1] = 4;
697 s1[2] = 5;
698
699 TArrayI inner(1);
700 inner[0] = 0;
701
702 TArrayI outer(1);
703 outer[0] = 1;
704
705 // Just to get the right (maximum) binning
706 TH1D *half[4];
707 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
708 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
709 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
710 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
711
712 for (int i=0; i<4; i++)
713 {
714 half[i]->SetLineColor(kRed+i);
715 half[i]->SetDirectory(0);
716 half[i]->SetBit(kCanDelete);
717 half[i]->Draw("same");
718 }
719
720 gPad->Modified();
721 gPad->Update();
722
723 }
724}
725
726
Note: See TracBrowser for help on using the repository browser.