source: trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibrate_data_ped_standard.C@ 6232

Last change on this file since 6232 was 4998, checked in by hbartko, 21 years ago
*** empty log message ***
File size: 22.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): Hendrik Bartko, 03/2004 <mailto:hbartko@mppmu.mpg.de>
19! Markus Gaug, 03/2004 <mailto:markus@ifae.es>
20! Wolfgang Wittek 07/2004 <mailto:wittek@mppmu.mpg.de>
21! Daniel Mazin 08/2004 <mailto:mazin@mppmu.mpg.de>
22!
23! Copyright: MAGIC Software Development, 2000-2004
24!
25!
26\* ======================================================================== */
27
28//-------------------------------------
29//
30// changes by W. Wittek :
31// - call SetRange() and SetWindows() for
32// MExtractedFixedWindowPeakSearch
33// - call SetRange() and SetWindowSize() for MPedCalcPedRun
34// - call MPedCalcFromData in 3rd loop (data calibration)
35//
36//-------------------------------------
37
38
39
40//const TString defpath = "/local/rootdata/2004_08_23/";
41
42// const TString defpath = "/mnt/data21a/Crab_August/";
43
44 const TString defpath = "/home/pcmagic04/pratik/analysis/2004_08_14/";
45
46//const TString defcontinuoslightfile =
47// "/.magic/magicserv01/MAGIC/rootdata2/2004_04_22/20040813_32834_C_3EG1727+04-100CL_E.root";
48const TString defcontinuoslightfile =
49 "/mnt/data21a/Crab_August/20040820_34060_C_CrabNebula-100CL_E.root";
50
51// ON data Mkn 421
52
53const TString defrout = "34994_test.root";
54const Int_t defpedr = 10060; // {17191,17192};
55const Int_t defcalr = 10051; // {17193,17198,17200};
56const Int_t defdatar = 23210; // {17206,17207,17208};
57
58
59//void calibrate_data_pedestal(const TString inpath=defpath,
60// const Int_t pedrun=defpedr,
61// const Int_t calrun=defcalr,
62// const Int_t datarun=defdatar,
63// const TString continuoslightfile=defcontinuoslightfile,
64// const TString resname=defrout)
65
66void calibrate_data_ped_standard(const TString inpath = defpath,
67 const Int_t pedrun = defpedr,
68 const Int_t calrun = defcalr,
69 const Int_t datarun = defdatar,
70 const TString continuoslightfile = defcontinuoslightfile,
71 const TString resname = defrout)
72
73{
74 cout << " inpath : " << inpath << endl;
75 cout << " pedrun : " << pedrun << endl;
76 cout << " calrun : " << calrun << endl;
77 cout << " datarun : " << datarun << endl;
78 cout << " continuoslightfile : " << continuoslightfile << endl;
79 cout << " FILE to write : " << resname << endl;
80
81 gLog.SetNoColors();
82 // the parameters should be identical to the ones of MPedCalcPedRun
83
84 Int_t WindowSize = 6; // size of the chosen integration window for signal and time extraction
85
86 MExtractFixedWindowPeakSearch extractor;
87 extractor.SetRange(0, 14, 0, 14);
88 extractor.SetWindows( WindowSize, WindowSize, 4);
89
90
91
92 MExtractTimeHighestIntegral timeext;
93 timeext.SetRange(0, 14, 0, 14);
94 timeext.SetWindowSize(WindowSize,WindowSize);
95
96
97 MRunIter pruns;
98 MRunIter cruns;
99 MRunIter druns;
100
101/*
102 for (Int_t i=0;i<psize;i++) {
103 cout << "Adding pedestal run: " << pedruns[i] << endl;
104 pruns.AddRun(pedruns[i],inpath);
105 }
106 for (Int_t i=0;i<csize;i++) {
107 cout << "Adding calibration run: " << calruns[i] << endl;
108 cruns.AddRun(calruns[i],inpath);
109 }
110 for (Int_t i=0;i<dsize;i++) {
111 cout << "Adding data run: " << dataruns[i] << endl;
112 druns.AddRun(dataruns[i],inpath);
113 }
114*/
115
116 cout << "Adding pedestal run: " << pedrun << endl;
117 pruns.AddRun(pedrun,inpath);
118 cout << "Adding calibration run: " << calrun << endl;
119 cruns.AddRun(calrun,inpath);
120 cruns.AddRun(10001,inpath);
121cruns.AddRun(10002,inpath);
122cruns.AddRun(10003,inpath);
123cruns.AddRun(10004,inpath);
124cruns.AddRun(10005,inpath);
125cruns.AddRun(10006,inpath);
126cruns.AddRun(10007,inpath);
127cruns.AddRun(10008,inpath);
128cruns.AddRun(10009,inpath);
129cruns.AddRun(10010,inpath);
130cruns.AddRun(10011,inpath);
131cruns.AddRun(10012,inpath);
132cruns.AddRun(10013,inpath);
133cruns.AddRun(10014,inpath);
134cruns.AddRun(10015,inpath);
135cruns.AddRun(10016,inpath);
136cruns.AddRun(10017,inpath);
137cruns.AddRun(10018,inpath);
138cruns.AddRun(10019,inpath);
139cruns.AddRun(10020,inpath);
140cruns.AddRun(10021,inpath);
141cruns.AddRun(10022,inpath);
142cruns.AddRun(10023,inpath);
143cruns.AddRun(10024,inpath);
144cruns.AddRun(10025,inpath);
145cruns.AddRun(10026,inpath);
146cruns.AddRun(10027,inpath);
147cruns.AddRun(10028,inpath);
148cruns.AddRun(10029,inpath);
149cruns.AddRun(10030,inpath);
150cruns.AddRun(10031,inpath);
151cruns.AddRun(10032,inpath);
152cruns.AddRun(10033,inpath);
153cruns.AddRun(10034,inpath);
154cruns.AddRun(10035,inpath);
155cruns.AddRun(10036,inpath);
156cruns.AddRun(10037,inpath);
157cruns.AddRun(10038,inpath);
158cruns.AddRun(10039,inpath);
159cruns.AddRun(10040,inpath);
160cruns.AddRun(10041,inpath);
161cruns.AddRun(10042,inpath);
162cruns.AddRun(10043,inpath);
163cruns.AddRun(10044,inpath);
164cruns.AddRun(10045,inpath);
165cruns.AddRun(10046,inpath);
166cruns.AddRun(10047,inpath);
167cruns.AddRun(10048,inpath);
168cruns.AddRun(10049,inpath);
169cruns.AddRun(10050,inpath);
170
171
172 cout << "Adding data run: " << datarun << endl;
173 druns.AddRun(datarun,inpath);
174
175
176 MStatusDisplay *display = new MStatusDisplay;
177 display->SetUpdateTime(3000);
178 display->Resize(850,700);
179
180 gStyle->SetOptStat(1111);
181 gStyle->SetOptFit();
182
183 /************************************/
184 /* FIRST LOOP: PEDESTAL COMPUTATION */
185 /************************************/
186
187 MParList plist1;
188 MTaskList tlist1;
189 plist1.AddToList(&tlist1);
190
191 // containers
192 MPedestalCam pedcam;
193 MBadPixelsCam badcam;
194 //
195 // for excluding pixels from the beginning:
196 //
197 // badcam.AsciiRead("badpixels.dat");
198
199
200 plist1.AddToList(&pedcam);
201 plist1.AddToList(&badcam);
202
203 //tasks
204 MReadMarsFile read("Events");
205 read.DisableAutoScheme();
206 static_cast<MRead&>(read).AddFiles(pruns);
207
208 MGeomApply geomapl;
209
210 // the parameters should be identical to the ones of the extractor
211 MPedCalcPedRun pedcalc;
212 pedcalc.SetRange(0,14,0,14);
213 pedcalc.SetWindowSize( WindowSize, WindowSize);
214
215 MGeomCamMagic geomcam;
216
217 tlist1.AddToList(&read);
218 tlist1.AddToList(&geomapl);
219 tlist1.AddToList(&pedcalc);
220
221 // Create and execute the event looper
222 MEvtLoop pedloop;
223 pedloop.SetParList(&plist1);
224 pedloop.SetDisplay(display);
225
226 cout << "*************************" << endl;
227 cout << "** COMPUTING PEDESTALS **" << endl;
228 cout << "*************************" << endl;
229
230 if (!pedloop.Eventloop())
231 return;
232
233 tlist1.PrintStatistics();
234
235 //
236 // Now the short version:
237 //
238 //
239 // Now setup the new tasks for the calibration:
240 // ---------------------------------------------------
241 //
242 MJCalibration calloop;
243 calloop.SetInput(&cruns);
244
245
246// const MCalibrationCam::PulserColor_t color=MCalibrationCam::kGREEN;
247//cout << "COLOR = " << color << endl;
248
249
250 //calloop.SetColor(color); // neu
251 // calloop.SetFullDisplay();
252 //
253 calloop.SetExtractor(&extractor);
254 //
255 // Set the corr. cams:
256 //
257 calloop.SetBadPixels(badcam);
258 //
259 // The next two commands are for the display:
260 //
261 calloop.SetDisplay(display);
262// calloop.SetFullDisplay(); //new
263
264 //
265 // Apply rel. time calibration:
266 //
267 //calloop.SetRelTimeCalibration();
268 //calloop.SetTimeExtractor(&timeext);
269 //
270 // Do the event-loop:
271 //
272 cout << "***************************" << endl;
273 cout << "** COMPUTING CALIBRATION **" << endl;
274 cout << "***************************" << endl;
275
276 if (!calloop.Process(pedcam))
277 return;
278cout << " end of Process " << endl;
279
280 badcam.Print();
281
282 MBadPixelsCam &badbad = calloop.GetBadPixels();
283 MCalibrationChargeCam &calcam = calloop.GetCalibrationCam();
284 MCalibrationRelTimeCam &timecam = calloop.GetRelTimeCam();
285 MCalibrationQECam &qecam = calloop.GetQECam();
286
287 badbad.Print();
288
289 /************************************************************************/
290 /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */
291 /************************************************************************/
292
293
294 // First: a small loop to calculate the pedestals from the data for the first some events
295
296
297 // Create an empty Parameter List and an empty Task List
298
299
300 MParList plist3a;
301 MTaskList tlist3a;
302 plist3a.AddToList(&tlist3a);
303
304 MPedestalCam pedcamdata;
305 pedcamdata.SetName("MPedestalCamFromData");
306
307
308 plist3a.AddToList(&pedcamdata);
309
310
311 MRawRunHeader runhead3a;
312 plist3a.AddToList(&geomcam );
313 plist3a.AddToList(&runhead3a);
314
315
316
317 //tasks
318
319 MPedCalcFromLoGain peddatacalc;
320
321
322 peddatacalc.SetPedContainerName("MPedestalCamFromData");
323 peddatacalc.SetPedestalUpdate(kFALSE);
324 peddatacalc.SetRange( 0, 11, 1, 14);
325 peddatacalc.SetWindowSize( 12, 14);
326 peddatacalc.SetNumEventsDump(500);
327 peddatacalc.SetMaxHiGainVar(40);
328
329
330 MReadMarsFile read3a("Events");
331 read3a.DisableAutoScheme();
332 static_cast<MRead&>(read3a).AddFiles(druns);
333
334 tlist3a.AddToList(&read3a);
335 tlist3a.AddToList(&geomapl);
336 tlist3a.AddToList(&peddatacalc);
337
338
339 MEvtLoop evtloop3a;
340 evtloop3a.SetParList(&plist3a);
341
342 cout << " event loop over the first 500 data events for pedestal calculation" << endl;
343
344 if ( !evtloop3a.Eventloop(500) )
345 return;
346
347// if (!evtloop3a.PreProcess())
348// return;
349
350// for (int i=0; i<500; i++){
351// if (!tlist3a.Process()) break;
352// }
353
354// evtloop3a.PostProcess();
355
356
357 // Now the main third loop for the cosmics and pedestal calibration
358
359
360 MParList plist3;
361 MTaskList tlist3;
362 plist3.AddToList(&tlist3);
363
364
365 // containers
366
367 MCerPhotEvt photevt;
368 MPedPhotCam pedphotcam;
369
370 MPedPhotCam pedphotcamdata;
371 pedphotcamdata.SetName("MPedPhotCamFromData");
372
373 MSrcPosCam srccam;
374 MRawRunHeader runhead;
375 MExtractedSignalCam sigcam;
376// MSrcPosCam pntcam; //new!
377// pntcam.SetName("MPntPosCam");
378
379// MCameraDC dccam;
380// MStarCam starcam;
381
382// changes by wolfgang:
383 MObservatory obs;
384// MStarCam sourcecam;
385// sourcecam.SetName("MSourceCam");
386// MHTelAxisFromStars htelaxis;
387 plist3.AddToList(&obs);
388// plist3.AddToList(&htelaxis);
389// plist3.AddToList(&sourcecam);
390// plist3.AddToList(&pntcam); //new!
391// end changes
392
393 plist3.AddToList(&geomcam );
394 plist3.AddToList(&pedcam);
395 plist3.AddToList(&pedcamdata);
396 plist3.AddToList(&calcam );
397 plist3.AddToList(&qecam );
398 plist3.AddToList(&badbad );
399 plist3.AddToList(&timecam );
400 plist3.AddToList(&sigcam );
401 plist3.AddToList(&photevt);
402 plist3.AddToList(&pedphotcam);
403 plist3.AddToList(&pedphotcamdata);
404 plist3.AddToList(&srccam);
405 plist3.AddToList(&runhead);
406// plist3.AddToList(&dccam);
407// plist3.AddToList(&starcam);
408
409
410 //tasks
411// MReadMarsFile read3("Events");
412
413 MReadReports read3;
414 read3.AddTree("Currents");
415 read3.AddTree("Drive");
416 read3.AddTree("Events","MTime.",kTRUE);
417
418// read3.DisableAutoScheme();
419 static_cast<MRead&>(read3).AddFiles(druns);
420
421 read3.AddToBranchList("MReportCurrents.*");
422 read3.AddToBranchList("MReportDrive.*");
423 read3.AddToBranchList("MRawCrateArray.*");
424 read3.AddToBranchList("MRawEvtData.*");
425 read3.AddToBranchList("MRawEvtHeader.*");
426
427
428 Float_t mindc = 0.9; //[uA]
429 MCalibrateDC dccal;
430 dccal.SetFileName(continuoslightfile);
431 dccal.SetMinDCAllowed(mindc);
432
433 const Int_t numblind = 5;
434 const Short_t x[numblind] = { 47, 124, 470, 475, 571};
435 const TArrayS blindpixels(numblind,(Short_t*)x);
436 Float_t ringinterest = 100; //[mm]
437 Float_t tailcut = 2.5;
438 UInt_t integratedevents = 1;
439
440 // We need the MAGIC mirror geometry from a MC header:
441// TString geometryfile = "../Mars/Gamma_zbin9_90_7_1480to1489_w0.root";
442
443 TString geometryfile = "/.magic/magicserv01/MAGIC/mcdata/Period016/spot_1cm/standard/gamma/Gamma_zbin9_90_7_1740to1749_w0.root";
444
445 // We need the Bright Star Catalog:
446
447 // TString catalogfile = "mtemp/mmpi/macros/bsc5.dat";
448
449 TString catalogfile = "bsc5.dat";
450
451 MFindStars findstars;
452 // findstars.SetBlindPixels(blindpixels); ?? WHY?
453 findstars.SetRingInterest(ringinterest);
454 findstars.SetDCTailCut(tailcut);
455 findstars.SetNumIntegratedEvents(integratedevents);
456 findstars.SetMinuitPrintOutLevel(-1);
457 findstars.SetGeometryFile(geometryfile);
458 findstars.SetBSCFile(catalogfile);
459 const Double_t ra = MAstro::Hms2Rad(11, 4, 26);
460 const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
461// const Double_t ra = MAstro::Hms2Rad(20, 0, 0);
462// const Double_t dec = MAstro::Dms2Rad(20, 48, 0);
463 findstars.SetRaDec(ra,dec);
464 findstars.SetLimMag(8);
465 findstars.SetRadiusFOV(1.5);
466
467
468 //$$$$$$$$$$$$$$$$ ww
469/*
470 MSourceDirections sdirs;
471 sdirs.SetGeometryFile(geometryfile);
472 const Double_t ra = MAstro::Hms2Rad(11, 4, 26);
473 const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
474 sdirs.SetRaDec(ra,dec);
475 sdirs.AddDirection(ra,dec,1,"Mkn 421");
476 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
477 const Double_t dec = MAstro::Dms2Rad(38, 14, 29);
478 sdirs.AddDirection(ra,dec,1,"My_UMa 51");
479 sdirs.SetRadiusFOV(3);
480
481 Int_t InputType=1;
482 MTelAxisFromStars telaxis;
483 telaxis.SetInputType(InputType);
484 MFillH fillhisto("MHTelAxisFromStars[MHTelAxisFromStars]","");
485 htelaxis.SetInputType(InputType);
486*/
487 //$$$$$$$$$$$$$$$$ ww
488
489
490 peddatacalc.SetPedestalUpdate(kTRUE);
491 MArrivalTimeCalc2 timecalc;
492
493
494 MCalibrateData photcalc;
495 photcalc.SetCalibrationMode(MCalibrateData::kFfactor);
496 photcalc.EnablePedestalType(MCalibrateData::kEvent);
497 photcalc.EnablePedestalType(MCalibrateData::kRun);
498
499
500 MBadPixelsCalc badcalc;
501
502 MPointingPosCalc pntposcalc;
503// MPointingPosCalc mpntposcalc; // not necessary, wolfgang calculates by himself
504
505 // for the interpolation of the bad pixels comment the following lines in
506 // done in the next step (Analysis.C)
507 /*
508 MBadPixelsTreat badtreat;
509 badtreat.SetUseInterpolation();
510 badtreat.SetProcessRMS();
511 badtreat.SetNumMinNeighbors(int);
512 */
513
514
515 // set the names of the pedestal containers to be used:
516// extractor.SetNamePedContainer("MPedestalCamFromData");
517 extractor.SetNamePedestalCam("MPedestalCamFromData");
518 timeext.SetNamePedestalCam("MPedestalCamFromData");
519// photcalc.SetNamePedADCRunContainer("MPedestalCam");
520// photcalc.SetNamePedADCEventContainer("MPedestalCamFromData");
521// photcalc.SetNamePedPhotRunContainer("MPedPhotCam");
522// photcalc.SetNamePedPhotEventContainer("MPedPhotCamFromData");
523 photcalc.SetNamePedADCContainer("MPedestalCamFromData");
524 photcalc.SetNamePedPhotContainer("MPedPhotCamFromData");
525
526 badcalc.SetNamePedPhotContainer("MPedPhotCamFromData");
527
528
529
530 tlist3.AddToList(&read3);
531 tlist3.AddToList(&geomapl);
532// tlist3.AddToList(&dccal, "Currents");
533// tlist3.AddToList(&findstars, "Currents");
534// tlist3.AddToList(&sdirs);
535// tlist3.AddToList(&telaxis);
536// tlist3.AddToList(&mpntposcalc);
537 tlist3.AddToList(&pntposcalc);
538// tlist3.AddToList(&fillhisto); // changed by wolfgang
539
540 tlist3.AddToList(&peddatacalc, "Events");
541 tlist3.AddToList(&extractor, "Events");
542 tlist3.AddToList(&timeext, "Events");
543 tlist3.AddToList(&photcalc, "Events");
544 tlist3.AddToList(&badcalc, "Events");
545
546
547 MWriteRootFile write(resname);
548
549 write.AddContainer("MGeomCam" , "RunHeaders");
550 write.AddContainer("MRawRunHeader" , "RunHeaders");
551// write.AddContainer("MSrcPosCam" , "RunHeaders");
552 write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
553 write.AddContainer("MCalibrationQECam" , "RunHeaders");
554 write.AddContainer("MCalibrationRelTimeCam", "RunHeaders");
555
556 write.AddContainer("MTime" , "Events");
557 write.AddContainer("MPedestalCam" , "Events");
558 write.AddContainer("MPedestalCamFromData","Events");
559 write.AddContainer("MCerPhotEvt" , "Events");
560 write.AddContainer("MRawEvtHeader" , "Events");
561 write.AddContainer("MBadPixelsCam" , "Events");
562 write.AddContainer("MPedPhotCamFromData", "Events");
563 write.AddContainer("MArrivalTimeCam", "Events");
564
565// write.AddContainer("MStarCam", "Events");
566 write.AddContainer("MSrcPosCam", "Events");
567
568 write.AddContainer("MReportDrive", "Events"); //new
569// write.AddContainer("MPntPosCam", "Events");
570 write.AddContainer("MPointingPos", "Events");
571
572 tlist3.AddToList(&write);
573
574 // Create and execute eventloop
575 MEvtLoop evtloop3;
576 evtloop3.SetParList(&plist3);
577
578 cout << "*************************************************************" << endl;
579 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
580 cout << "*************************************************************" << endl;
581
582 int maxevents = 3000;
583
584 if (!evtloop3.Eventloop(maxevents))
585 return;
586 tlist3.PrintStatistics();
587
588// changes by wolfgang
589// TObject *obj = plist3.FindObject("MHTelAxisFromStars");
590// obj->DrawClone();
591
592}
593
594void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
595{
596
597 c.cd(i);
598 gPad->SetBorderMode(0);
599 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
600 // obj1->AddNotify(evt);
601
602 c.cd(i+j);
603 gPad->SetBorderMode(0);
604 obj1->Draw();
605 ((MHCamera*)obj1)->SetPrettyPalette();
606
607 if (fit != 0)
608 {
609 c.cd(i+2*j);
610 gPad->SetBorderMode(0);
611 TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
612
613// obj2->Sumw2();
614 obj2->Draw();
615 obj2->SetBit(kCanDelete);
616
617 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
618 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
619 const Double_t integ = obj2->Integral("width")/2.5066283;
620 const Double_t mean = obj2->GetMean();
621 const Double_t rms = obj2->GetRMS();
622 const Double_t width = max-min;
623
624 if (rms == 0. || width == 0. )
625 return;
626
627 switch (fit)
628 {
629 case 1:
630 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
631 sgaus->SetBit(kCanDelete);
632 sgaus->SetParNames("Area","#mu","#sigma");
633 sgaus->SetParameters(integ/rms,mean,rms);
634 sgaus->SetParLimits(0,0.,integ);
635 sgaus->SetParLimits(1,min,max);
636 sgaus->SetParLimits(2,0,width/1.5);
637 obj2->Fit("sgaus","QLR");
638 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
639 break;
640
641 case 2:
642 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
643 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
644 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
645 dgaus->SetBit(kCanDelete);
646 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
647 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
648 integ/width/2.,(max+mean)/2.,width/4.);
649 // The left-sided Gauss
650 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
651 dgaus->SetParLimits(1,min+(width/10.),mean);
652 dgaus->SetParLimits(2,0,width/2.);
653 // The right-sided Gauss
654 dgaus->SetParLimits(3,0,integ);
655 dgaus->SetParLimits(4,mean,max-(width/10.));
656 dgaus->SetParLimits(5,0,width/2.);
657 obj2->Fit("dgaus","QLRM");
658 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
659 break;
660
661 case 3:
662 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
663 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
664 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
665 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
666 tgaus->SetBit(kCanDelete);
667 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
668 "A_{2}","#mu_{2}","#sigma_{2}",
669 "A_{3}","#mu_{3}","#sigma_{3}");
670 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
671 integ/width/3.,(max+mean)/2.,width/4.,
672 integ/width/3.,mean,width/2.);
673 // The left-sided Gauss
674 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
675 tgaus->SetParLimits(1,min+(width/10.),mean);
676 tgaus->SetParLimits(2,width/15.,width/2.);
677 // The right-sided Gauss
678 tgaus->SetParLimits(3,0.,integ);
679 tgaus->SetParLimits(4,mean,max-(width/10.));
680 tgaus->SetParLimits(5,width/15.,width/2.);
681 // The Gauss describing the outliers
682 tgaus->SetParLimits(6,0.,integ);
683 tgaus->SetParLimits(7,min,max);
684 tgaus->SetParLimits(8,width/4.,width/1.5);
685 obj2->Fit("tgaus","QLRM");
686 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
687 break;
688 case 4:
689 obj2->Fit("pol0","Q");
690 obj2->GetFunction("pol0")->SetLineColor(kYellow);
691 break;
692 case 9:
693 break;
694 default:
695 obj2->Fit("gaus","Q");
696 obj2->GetFunction("gaus")->SetLineColor(kYellow);
697 break;
698 }
699
700 TArrayI s0(3);
701 s0[0] = 6;
702 s0[1] = 1;
703 s0[2] = 2;
704
705 TArrayI s1(3);
706 s1[0] = 3;
707 s1[1] = 4;
708 s1[2] = 5;
709
710 TArrayI inner(1);
711 inner[0] = 0;
712
713 TArrayI outer(1);
714 outer[0] = 1;
715
716 // Just to get the right (maximum) binning
717 TH1D *half[4];
718 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
719 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
720 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
721 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
722
723 for (int i=0; i<4; i++)
724 {
725 half[i]->SetLineColor(kRed+i);
726 half[i]->SetDirectory(0);
727 half[i]->SetBit(kCanDelete);
728 half[i]->Draw("same");
729 }
730
731 gPad->Modified();
732 gPad->Update();
733
734 }
735}
736
737
Note: See TracBrowser for help on using the repository browser.