source: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C@ 2551

Last change on this file since 2551 was 2543, checked in by wittek, 22 years ago
*** empty log message ***
File size: 95.0 KB
Line 
1
2
3#include "CT1EgyEst.C"
4
5
6void InitBinnings(MParList *plist)
7{
8 gLog << "InitBinnings" << endl;
9
10 //--------------------------------------------
11 MBinning *binse = new MBinning("BinningE");
12 //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
13
14 //This is Daniel's binning in energy:
15 binse->SetEdgesLog(14, 296.296, 86497.6);
16 plist->AddToList(binse);
17
18 //--------------------------------------------
19
20 MBinning *binssize = new MBinning("BinningSize");
21 binssize->SetEdgesLog(50, 10, 1.0e5);
22 plist->AddToList(binssize);
23
24 MBinning *binsdistc = new MBinning("BinningDist");
25 binsdistc->SetEdges(50, 0, 1.4);
26 plist->AddToList(binsdistc);
27
28 MBinning *binswidth = new MBinning("BinningWidth");
29 binswidth->SetEdges(50, 0, 1.0);
30 plist->AddToList(binswidth);
31
32 MBinning *binslength = new MBinning("BinningLength");
33 binslength->SetEdges(50, 0, 1.0);
34 plist->AddToList(binslength);
35
36 MBinning *binsalpha = new MBinning("BinningAlpha");
37 binsalpha->SetEdges(100, -100, 100);
38 plist->AddToList(binsalpha);
39
40 MBinning *binsasym = new MBinning("BinningAsym");
41 binsasym->SetEdges(50, -1.5, 1.5);
42 plist->AddToList(binsasym);
43
44 MBinning *binsm3l = new MBinning("BinningM3Long");
45 binsm3l->SetEdges(50, -1.5, 1.5);
46 plist->AddToList(binsm3l);
47
48 MBinning *binsm3t = new MBinning("BinningM3Trans");
49 binsm3t->SetEdges(50, -1.5, 1.5);
50 plist->AddToList(binsm3t);
51
52
53 //.....
54 MBinning *binsb = new MBinning("BinningSigmabar");
55 binsb->SetEdges( 100, 0.0, 5.0);
56 plist->AddToList(binsb);
57
58 MBinning *binth = new MBinning("BinningTheta");
59 // this is Daniel's binning in theta
60 //Double_t yedge[8] =
61 // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
62 // this is our binning
63 Double_t yedge[9] =
64 {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
65 TArrayD yed;
66 yed.Set(9,yedge);
67 binth->SetEdges(yed);
68 plist->AddToList(binth);
69
70 MBinning *bincosth = new MBinning("BinningCosTheta");
71 Double_t zedge[9];
72 for (Int_t i=0; i<9; i++)
73 {
74 zedge[8-i] = cos(yedge[i]/kRad2Deg);
75 }
76 TArrayD zed;
77 zed.Set(9,zedge);
78 bincosth->SetEdges(zed);
79 plist->AddToList(bincosth);
80
81 MBinning *binsdiff = new MBinning("BinningDiffsigma2");
82 binsdiff->SetEdges(100, -5.0, 20.0);
83 plist->AddToList(binsdiff);
84
85 // robert ----------------------------------------------
86 MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
87 binsalphaf->SetEdges(100, -100, 100);
88 plist->AddToList(binsalphaf);
89
90 MBinning *binsdifftime = new MBinning("BinningTimeDiff");
91 binsdifftime->SetEdges(50, 0, 10);
92 plist->AddToList(binsdifftime);
93
94 MBinning *binstime = new MBinning("BinningTime");
95 binstime->SetEdges(50, 44500, 61000);
96 plist->AddToList(binstime);
97}
98
99
100void DeleteBinnings(MParList *plist)
101{
102 gLog << "DeleteBinnings" << endl;
103
104 TObject *bin;
105
106 //--------------------------------------------
107 bin = plist->FindObject("BinningE");
108 if (bin) delete bin;
109
110 //--------------------------------------------
111
112 bin = plist->FindObject("BinningSize");
113 if (bin) delete bin;
114
115 bin = plist->FindObject("BinningDist");
116 if (bin) delete bin;
117
118 bin = plist->FindObject("BinningWidth");
119 if (bin) delete bin;
120
121 bin = plist->FindObject("BinningLength");
122 if (bin) delete bin;
123
124 bin = plist->FindObject("BinningAlpha");
125 if (bin) delete bin;
126
127 bin = plist->FindObject("BinningAsym");
128 if (bin) delete bin;
129
130 bin = plist->FindObject("BinningM3Long");
131 if (bin) delete bin;
132
133 bin = plist->FindObject("BinningM3Trans");
134 if (bin) delete bin;
135
136 //.....
137 bin = plist->FindObject("BinningSigmabar");
138 if (bin) delete bin;
139
140 bin = plist->FindObject("BinningTheta");
141 if (bin) delete bin;
142
143 bin = plist->FindObject("BinningCosTheta");
144 if (bin) delete bin;
145
146 bin = plist->FindObject("BinningDiffsigma2");
147 if (bin) delete bin;
148
149
150 // robert ----------------------------------------------
151 bin = plist->FindObject("BinningAlphaFlux");
152 if (bin) delete bin;
153
154 bin = plist->FindObject("BinningTimeDiff");
155 if (bin) delete bin;
156
157 bin = plist->FindObject("BinningTime");
158 if (bin) delete bin;
159}
160
161
162
163//************************************************************************
164void ONOFFCT1Analysis()
165{
166 gLog.SetNoColors();
167
168 if (gRandom)
169 delete gRandom;
170 gRandom = new TRandom3(0);
171
172 //-----------------------------------------------
173 const char *offfile = "~magican/ct1test/wittek/offdata.preproc";
174
175 //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc";
176 const char *onfile = "~magican/ct1test/wittek/mkn421_00-01";
177
178 const char *mcfile = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc";
179 //-----------------------------------------------
180
181 // path for input for Mars
182 TString inPath = "~magican/ct1test/wittek/marsoutput/optionD/";
183
184 // path for output from Mars
185 TString outPath = "~magican/ct1test/wittek/marsoutput/optionD/";
186
187 //-----------------------------------------------
188
189 //TEnv env("macros/CT1env.rc");
190 //Bool_t printEnv = kFALSE;
191
192 //************************************************************************
193
194 // Job A :
195 // - produce MHSigmaTheta plots for ON and OFF data
196 // - write out (or read in) these MHSigmaTheta plots
197 // - read ON (or OFF or MC) data
198 // - pad the events;
199 // - write root file for ON (or OFF or MC) data (ON1.root, ...);
200
201 Bool_t JobA = kFALSE;
202 Bool_t WPad = kFALSE; // write out padding histograms ?
203 Bool_t RPad = kFALSE; // read in padding histograms ?
204 Bool_t Wout = kFALSE; // write out root file ON1.root
205 // (or OFF1.root or MC1.root)?
206
207
208 // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file
209 // - if CTrainRF = TRUE : create matrices of training events
210 // and root files of training and test events
211 // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
212 // - if RTree = TRUE : read in trees, otherwise create trees
213 // - calculate hadroness for method of RANDOM FOREST
214 // - update the input files with the hadronesses (ON2.root, OFF2.root
215 // or MC2.root)
216
217 Bool_t JobB_RF_UP = kTRUE;
218 Bool_t CTrainRF = kTRUE; // create matrices of training events
219 // and root files of training and test events
220 Bool_t RTrainRF = kFALSE; // read in matrices of training events
221 Bool_t RTree = kFALSE; // read in trees (otherwise grow them)
222 Bool_t WRF = kFALSE; // update input root file ?
223
224
225 // Job B_SC_UP : read ON2.root (or MC2.root) file
226 // - depending on WParSC : create (or read in) supercuts parameter values
227 // - calculate hadroness for the SUPERCUTS
228 // - update the input files with the hadroness (==>ON3.root or MC3.root)
229
230 Bool_t JobB_SC_UP = kFALSE;
231 Bool_t CMatrix = kFALSE; // create training and test matrices
232 Bool_t RMatrix = kFALSE; // read training and test matrices from file
233 Bool_t WOptimize = kFALSE; // do optimization using the training sample
234 // and write supercuts parameter values
235 // onto the file parSCfile
236 Bool_t RTest = kFALSE; // test the supercuts using the test matrix
237 Bool_t WSC = kFALSE; // update input root file ?
238
239
240 // Job C:
241 // - read ON3.root and MC3.root files
242 // which should have been updated to contain the hadronnesses
243 // for the method of
244 // RF
245 // SUPERCUTS and
246 // - produce Neyman-Pearson plots
247
248 Bool_t JobC = kFALSE;
249
250
251 // Job D :
252 // - select g/h separation method XX
253 // - read ON3 (or MC3) root file
254 // - apply cuts in hadronness
255 // - make plots
256
257 Bool_t JobD = kFALSE;
258
259
260
261 // Job E_XX : extended version of E_XX (including flux plots)
262 // - select g/h separation method XX
263 // - read MC root file
264   // - calculate eff. collection area
265 // - optimize energy estimator
266 // - read ON root file
267 // - apply final cuts
268 // - calculate flux
269 // - write root file for ON data after final cuts
270
271
272 Bool_t JobE_XX = kFALSE;
273 Bool_t CCollArea= kFALSE; // calculate eff. collection areas
274 Bool_t OEEst = kFALSE; // optimize energy estimator
275 Bool_t WEX = kFALSE; // update root file ?
276 Bool_t WRobert = kFALSE; // write out Robert's file ?
277
278
279
280 //************************************************************************
281
282
283 //---------------------------------------------------------------------
284 // Job A
285 //=========
286 // read ON data file
287
288 // - produce the 2D-histogram "sigmabar versus Theta"
289 // (SigmaTheta_ON.root) for ON data
290 // (to be used for the padding of the MC gamma data)
291
292 // - write a file of ON events (ON1.root)
293 // (after the standard cuts, before the g/h separation)
294 // (to be used together with the corresponding MC gamma file (MC1.root)
295 // for the optimization of the g/h separation)
296
297
298 if (JobA)
299 {
300 gLog << "=====================================================" << endl;
301 gLog << "Macro CT1Analysis : Start of Job A" << endl;
302 gLog << "" << endl;
303 gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = "
304 << (JobA ? "kTRUE" : "kFALSE") << ", "
305 << (WPad ? "kTRUE" : "kFALSE") << ", "
306 << (RPad ? "kTRUE" : "kFALSE") << ", "
307 << (Wout ? "kTRUE" : "kFALSE") << endl;
308
309
310 //--------------------------------------------------
311 // names of ON and OFF files to be read
312 // for generating the histograms to be used in the padding
313 TString fileON = onfile;
314 TString fileOFF = offfile;
315 gLog << "fileON, fileOFF = " << fileON << ", " << fileOFF << endl;
316
317 // name of file to conatin the histograms for the padding
318 TString outNameSigTh = outPath;
319 outNameSigTh += "SigmaTheta";
320 outNameSigTh += ".root";
321
322 //--------------------------------------------------
323 // type of data to be padded
324 //TString typeInput = "ON";
325 //TString typeInput = "OFF";
326 TString typeInput = "MC";
327 gLog << "typeInput = " << typeInput << endl;
328
329
330 // name of input root file
331 if (typeInput == "ON")
332 TString filenamein(onfile);
333 else if (typeInput == "OFF")
334 TString filenamein(offfile);
335 else if (typeInput == "MC")
336 TString filenamein(mcfile);
337 gLog << "data to be padded : " << filenamein << endl;
338
339 // name of output root file
340 TString outNameImage = outPath;
341 outNameImage += typeInput;
342 outNameImage += "1.root";
343 gLog << "padded data to be written onto : " << outNameImage << endl;
344
345 //--------------------------------------------------
346
347 //************************************************************
348 // generate histograms to be used in the padding
349 //
350 // read ON and OFF data files
351 // generate (or read in) the padding histograms for ON and OFF data
352 // and merge these histograms
353
354 MCT1PadONOFF pad;
355 pad.SetName("MCT1PadONOFF");
356 pad.SetPadFlag(1);
357 pad.SetDataType(typeInput);
358
359 // generate the padding histograms
360 if (!RPad)
361 {
362 gLog << "=====================================================" << endl;
363 gLog << "Start generating the padding histograms" << endl;
364 //-----------------------------------------
365 // ON events
366
367 gLog << "-----------" << endl;
368 gLog << "ON events :" << endl;
369 gLog << "-----------" << endl;
370
371 MTaskList tliston;
372 MParList pliston;
373
374 MCT1ReadPreProc readON(fileON);
375
376 //MFCT1SelBasic selthetaon;
377 //selthetaon.SetCuts(-100.0, 29.5, 35.5);
378 //MContinue contthetaon(&selthetaon);
379
380 MBlindPixelCalc blindon;
381 blindon.SetUseBlindPixels();
382
383 MFCT1SelBasic selbasicon;
384 MContinue contbasicon(&selbasicon);
385
386 MHBlindPixels blindON("BlindPixelsON");
387 MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels");
388 fillblindON.SetName("FillBlind");
389
390 MSigmabarCalc sigbarcalcon;
391
392 MHSigmaTheta sigthON("SigmaThetaON");
393 MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt");
394 fillsigthetaON.SetName("FillSigTheta");
395
396 //*****************************
397 // entries in MParList
398
399 pliston.AddToList(&tliston);
400 InitBinnings(&pliston);
401 pliston.AddToList(&blindON);
402 pliston.AddToList(&sigthON);
403
404
405 //*****************************
406 // entries in MTaskList
407
408 tliston.AddToList(&readON);
409 //tliston.AddToList(&contthetaon);
410
411 tliston.AddToList(&blindon);
412
413 tliston.AddToList(&contbasicon);
414 tliston.AddToList(&fillblindON);
415 tliston.AddToList(&sigbarcalcon);
416 tliston.AddToList(&fillsigthetaON);
417
418 MProgressBar baron;
419 MEvtLoop evtloopon;
420 evtloopon.SetParList(&pliston);
421 evtloopon.SetProgressBar(&baron);
422
423 Int_t maxeventson = -1;
424 //Int_t maxeventson = 10000;
425 if ( !evtloopon.Eventloop(maxeventson) )
426 return;
427
428 tliston.PrintStatistics(0, kTRUE);
429
430 blindON.DrawClone();
431 sigthON.DrawClone();
432
433 // save the histograms for the padding
434 TH2D *sigthon = sigthON.GetSigmaTheta();
435 TH3D *sigpixthon = sigthON.GetSigmaPixTheta();
436 TH3D *diffpixthon = sigthON.GetDiffPixTheta();
437
438 TH2D *blindidthon = blindON.GetBlindId();
439 TH2D *blindnthon = blindON.GetBlindN();
440
441 //-----------------------------------------
442 // OFF events
443
444 gLog << "------------" << endl;
445 gLog << "OFF events :" << endl;
446 gLog << "------------" << endl;
447
448 MTaskList tlistoff;
449 MParList plistoff;
450
451 MCT1ReadPreProc readOFF(fileOFF);
452
453 MFCT1SelBasic selthetaoff;
454 selthetaoff.SetCuts(-100.0, 29.5, 35.5);
455 MContinue contthetaoff(&selthetaoff);
456
457 MBlindPixelCalc blindoff;
458 blindoff.SetUseBlindPixels();
459
460 MFCT1SelBasic selbasicoff;
461 MContinue contbasicoff(&selbasicoff);
462
463 MHBlindPixels blindOFF("BlindPixelsOFF");
464 MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels");
465 fillblindOFF.SetName("FillBlindOFF");
466
467 MSigmabarCalc sigbarcalcoff;
468
469 MHSigmaTheta sigthOFF("SigmaThetaOFF");
470 MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt");
471 fillsigthetaOFF.SetName("FillSigThetaOFF");
472
473 //*****************************
474 // entries in MParList
475
476 plistoff.AddToList(&tlistoff);
477 InitBinnings(&plistoff);
478 plistoff.AddToList(&blindOFF);
479 plistoff.AddToList(&sigthOFF);
480
481
482 //*****************************
483 // entries in MTaskList
484
485 tlistoff.AddToList(&readOFF);
486 //tlistoff.AddToList(&contthetaoff);
487
488 tlistoff.AddToList(&blindoff);
489
490 tlistoff.AddToList(&contbasicoff);
491 tlistoff.AddToList(&fillblindOFF);
492 tlistoff.AddToList(&sigbarcalcoff);
493 tlistoff.AddToList(&fillsigthetaOFF);
494
495 MProgressBar baroff;
496 MEvtLoop evtloopoff;
497 evtloopoff.SetParList(&plistoff);
498 evtloopoff.SetProgressBar(&baroff);
499
500 Int_t maxeventsoff = -1;
501 //Int_t maxeventsoff = 20000;
502 if ( !evtloopoff.Eventloop(maxeventsoff) )
503 return;
504
505 tlistoff.PrintStatistics(0, kTRUE);
506
507 blindOFF.DrawClone();
508 sigthOFF.DrawClone();
509
510 // save the histograms for the padding
511 TH2D *sigthoff = sigthOFF.GetSigmaTheta();
512 TH3D *sigpixthoff = sigthOFF.GetSigmaPixTheta();
513 TH3D *diffpixthoff = sigthOFF.GetDiffPixTheta();
514
515 TH2D *blindidthoff = blindOFF.GetBlindId();
516 TH2D *blindnthoff = blindOFF.GetBlindN();
517
518
519 //-----------------------------------------
520
521 gLog << "End of generating the padding histograms" << endl;
522 gLog << "=====================================================" << endl;
523
524 pad.MergeHistograms(sigthon, sigthoff,
525 sigpixthon, sigpixthoff,
526 diffpixthon, diffpixthoff,
527 blindidthon, blindidthoff,
528 blindnthon, blindnthoff);
529
530 if (WPad)
531 {
532 // write the padding histograms onto a file ---------
533 pad.WriteTargetDist(outNameSigTh);
534 }
535 }
536
537 // read the padding histograms ---------------------------
538 if (RPad)
539 {
540 pad.ReadTargetDist(outNameSigTh);
541 }
542
543
544 //************************************************************
545
546 if (Wout)
547 {
548 gLog << "=====================================================" << endl;
549 gLog << "Start the padding" << endl;
550
551 //-----------------------------------------------------------
552 MTaskList tliston;
553 MParList pliston;
554
555 char *sourceName = "MSrcPosCam";
556 MSrcPosCam source(sourceName);
557
558 // geometry is needed in MHHillas... classes
559 MGeomCam *fGeom =
560 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
561
562 //-------------------------------------------
563 // create the tasks which should be executed
564 //
565
566 MCT1ReadPreProc read(filenamein);
567
568 MFCT1SelBasic seltheta;
569 seltheta.SetCuts(-100.0, 29.5, 35.5);
570 MContinue conttheta(&seltheta);
571
572 if (typeInput == "ON")
573 {
574 MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc",
575 "MCT1PointingCorrCalc");
576 }
577
578 MBlindPixelCalc blindbeforepad;
579 blindbeforepad.SetUseBlindPixels();
580 blindbeforepad.SetName("BlindBeforePadding");
581
582 MBlindPixelCalc blind;
583 blind.SetUseBlindPixels();
584 blind.SetName("BlindAfterPadding");
585
586 MFCT1SelBasic selbasic;
587 MContinue contbasic(&selbasic);
588 contbasic.SetName("SelBasic");
589
590 MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels");
591 fillblind.SetName("HBlind");
592
593 MSigmabarCalc sigbarcalc;
594
595 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
596 fillsigtheta.SetName("HSigmaTheta");
597
598 MImgCleanStd clean;
599
600
601 // calculation of image parameters ---------------------
602 TString fHilName = "MHillas";
603 TString fHilNameExt = "MHillasExt";
604 TString fHilNameSrc = "MHillasSrc";
605 TString fImgParName = "MNewImagePar";
606
607 MHillasCalc hcalc;
608 hcalc.SetNameHillas(fHilName);
609 hcalc.SetNameHillasExt(fHilNameExt);
610 hcalc.SetNameNewImgPar(fImgParName);
611
612 MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc);
613 hsrccalc.SetInput(fHilName);
614
615 MFillH hfill1("MHHillas", fHilName);
616 hfill1.SetName("HHillas");
617
618 MFillH hfill2("MHStarMap", fHilName);
619 hfill2.SetName("HStarMap");
620
621 MFillH hfill3("MHHillasExt", fHilNameSrc);
622 hfill3.SetName("HHillasExt");
623
624 MFillH hfill4("MHHillasSrc", fHilNameSrc);
625 hfill4.SetName("HHillasSrc");
626
627 MFillH hfill5("MHNewImagePar", fImgParName);
628 hfill5.SetName("HNewImagePar");
629 // --------------------------------------------------
630
631 MFCT1SelStandard selstandard(fHilNameSrc);
632 selstandard.SetHillasName(fHilName);
633 selstandard.SetImgParName(fImgParName);
634 selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0);
635 MContinue contstandard(&selstandard);
636 contstandard.SetName("SelStandard");
637
638
639 MWriteRootFile write(outNameImage);
640
641 write.AddContainer("MRawRunHeader", "RunHeaders");
642 write.AddContainer("MTime", "Events");
643 write.AddContainer("MMcEvt", "Events");
644 write.AddContainer("ThetaOrig", "Events");
645 write.AddContainer("MSrcPosCam", "Events");
646 write.AddContainer("MSigmabar", "Events");
647 write.AddContainer("MHillas", "Events");
648 write.AddContainer("MHillasExt", "Events");
649 write.AddContainer("MHillasSrc", "Events");
650 write.AddContainer("MNewImagePar", "Events");
651
652
653 //*****************************
654 // entries in MParList
655
656 pliston.AddToList(&tliston);
657 InitBinnings(&pliston);
658
659 pliston.AddToList(&source);
660
661
662 //*****************************
663 // entries in MTaskList
664
665 tliston.AddToList(&read);
666 //tliston.AddToList(&conttheta);
667
668 tliston.AddToList(&blindbeforepad);
669 tliston.AddToList(&pad);
670 if (typeInput == "ON")
671 tliston.AddToList(&pointcorr);
672 tliston.AddToList(&blind);
673
674 tliston.AddToList(&contbasic);
675 tliston.AddToList(&fillblind);
676 tliston.AddToList(&sigbarcalc);
677 tliston.AddToList(&fillsigtheta);
678 tliston.AddToList(&clean);
679
680 tliston.AddToList(&hcalc);
681 tliston.AddToList(&hsrccalc);
682
683 tliston.AddToList(&hfill1);
684 tliston.AddToList(&hfill2);
685 tliston.AddToList(&hfill3);
686 tliston.AddToList(&hfill4);
687 tliston.AddToList(&hfill5);
688
689 tliston.AddToList(&contstandard);
690 tliston.AddToList(&write);
691
692 //*****************************
693
694 //-------------------------------------------
695 // Execute event loop
696 //
697 MProgressBar bar;
698 MEvtLoop evtloop;
699 evtloop.SetParList(&pliston);
700 //evtloop.ReadEnv(env, "", printEnv);
701 evtloop.SetProgressBar(&bar);
702 // evtloop.Write();
703
704 Int_t maxevents = -1;
705 //Int_t maxevents = 1000;
706 if ( !evtloop.Eventloop(maxevents) )
707 return;
708
709 tliston.PrintStatistics(0, kTRUE);
710
711
712 //-------------------------------------------
713 // Display the histograms
714
715 pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone();
716 pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone();
717
718 pliston.FindObject("MHHillas")->DrawClone();
719 pliston.FindObject("MHHillasExt")->DrawClone();
720 pliston.FindObject("MHHillasSrc")->DrawClone();
721 pliston.FindObject("MHNewImagePar")->DrawClone();
722 pliston.FindObject("MHStarMap")->DrawClone();
723
724 DeleteBinnings(&pliston);
725
726 gLog << "End of padding" << endl;
727 gLog << "=====================================================" << endl;
728 }
729
730
731 gLog << "Macro CT1Analysis : End of Job A" << endl;
732 gLog << "===================================================" << endl;
733 }
734
735
736
737
738
739 //---------------------------------------------------------------------
740 // Job B_RF_UP
741 //============
742
743
744 // - create (or read in) the matrices of training events for gammas
745 // and hadrons
746 // - create (or read in) the trees
747 // - then read ON1.root (or MC1.root) file
748 // - calculate the hadroness for the method of RANDOM FOREST
749 // - update input root file with the hadroness
750
751
752 if (JobB_RF_UP)
753 {
754 gLog << "=====================================================" << endl;
755 gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
756
757 gLog << "" << endl;
758 gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
759 << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", "
760 << (RTrainRF ? "kTRUE" : "kFALSE") << ", "
761 << (CTrainRF ? "kTRUE" : "kFALSE") << ", "
762 << (RTree ? "kTRUE" : "kFALSE") << ", "
763 << (WRF ? "kTRUE" : "kFALSE") << endl;
764
765
766 //--------------------------------------------
767 // parameters for the random forest
768 Int_t NumTrees = 100;
769 Int_t NumTry = 3;
770 Int_t NdSize = 1;
771
772
773 TString hadRFName = "HadRF";
774 Float_t maxhadronness = 0.23;
775 Float_t maxalpha = 20.0;
776 Float_t maxdist = 10.0;
777
778 TString fHilName = "MHillas";
779 TString fHilNameExt = "MHillasExt";
780 TString fHilNameSrc = "MHillasSrc";
781 TString fImgParName = "MNewImagePar";
782
783
784 //--------------------------------------------
785 // file to be updated (ON, OFF or MC)
786
787 //TString typeInput = "ON";
788 //TString typeInput = "OFF";
789 TString typeInput = "MC";
790 gLog << "typeInput = " << typeInput << endl;
791
792 // name of input root file
793 TString filenameData = outPath;
794 filenameData += typeInput;
795 filenameData += "1.root";
796 gLog << "filenameData = " << filenameData << endl;
797
798 // name of output root file
799 TString outNameImage = outPath;
800 outNameImage += typeInput;
801 outNameImage += "2.root";
802 //TString outNameImage = filenameData;
803
804 gLog << "outNameImage = " << outNameImage << endl;
805
806 //--------------------------------------------
807 // files to be read for generating
808 // - the matrices of training events
809 // - and the root files of training and test events
810
811 TString typeMatrixHadrons = "OFF";
812 gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
813
814 // "hadrons" :
815 TString filenameHad = outPath;
816 filenameHad += typeMatrixHadrons;
817 filenameHad += "1.root";
818 Int_t howManyHadronsTrain = 8000;
819 Int_t howManyHadronsTest = 8000;
820 gLog << "filenameHad = " << filenameHad << ", howManyHadronsTrain = "
821 << howManyHadronsTrain << ", howManyHadronsTest = "
822 << howManyHadronsTest << endl;
823
824
825 // "gammas" :
826 TString filenameMC = outPath;
827 filenameMC += "MC";
828 filenameMC += "1.root";
829 Int_t howManyGammasTrain = 8000;
830 Int_t howManyGammasTest = 8000;
831 gLog << "filenameMC = " << filenameMC << ", howManyGammasTrain = "
832 << howManyGammasTrain << ", howManyGammasTest = "
833 << howManyGammasTest << endl;
834
835 //--------------------------------------------
836 // files for the matrices of training events
837
838 TString outNameGammas = outPath;
839 outNameGammas += "RFmatrix_gammas_Train_";
840 outNameGammas += "MC";
841 outNameGammas += ".root";
842
843 TString outNameHadrons = outPath;
844 outNameHadrons += "RFmatrix_hadrons_Train_";
845 outNameHadrons += typeMatrixHadrons;
846 outNameHadrons += ".root";
847
848
849 //--------------------------------------------
850 // root files for the training events
851
852 TString outNameGammasTrain = outPath;
853 outNameGammasTrain += "RF_gammas_Train_";
854 outNameGammasTrain += "MC";
855 outNameGammasTrain += ".root";
856
857 TString outNameHadronsTrain = outPath;
858 outNameHadronsTrain += "RF_hadrons_Train_";
859 outNameHadronsTrain += typeMatrixHadrons;
860 outNameHadronsTrain += ".root";
861
862
863 //--------------------------------------------
864 // root files for the test events
865
866 TString outNameGammasTest = outPath;
867 outNameGammasTest += "RF_gammas_Test_";
868 outNameGammasTest += "MC";
869 outNameGammasTest += ".root";
870
871 TString outNameHadronsTest = outPath;
872 outNameHadronsTest += "RF_hadrons_Test_";
873 outNameHadronsTest += typeMatrixHadrons;
874 outNameHadronsTest += ".root";
875
876 //--------------------------------------------------------------------
877
878
879 MHMatrix matrixg("MatrixGammas");
880 matrixg.EnableGraphicalOutput();
881
882 matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
883 matrixg.AddColumn("MSigmabar.fSigmabar");
884 matrixg.AddColumn("log10(MHillas.fSize)");
885 matrixg.AddColumn("MHillasSrc.fDist");
886 matrixg.AddColumn("MHillas.fWidth");
887 matrixg.AddColumn("MHillas.fLength");
888 matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
889 matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
890 matrixg.AddColumn("MNewImagePar.fConc");
891 matrixg.AddColumn("MNewImagePar.fLeakage1");
892
893 MHMatrix matrixh("MatrixHadrons");
894 matrixh.EnableGraphicalOutput();
895
896 matrixh.AddColumns(matrixg.GetColumns());
897
898 //--------------------------------------------
899 // file of trees of the random forest
900
901 TString outRF = outPath;
902 outRF += "RF.root";
903
904
905 //*************************************************************************
906 // read in matrices of training events
907if (RTrainRF)
908 {
909 const char* mtxName = "MatrixGammas";
910
911 gLog << "" << endl;
912 gLog << "========================================================" << endl;
913 gLog << "Get matrix for (gammas)" << endl;
914 gLog << "matrix name = " << mtxName << endl;
915 gLog << "name of root file = " << outNameGammas << endl;
916 gLog << "" << endl;
917
918
919 // read in the object with the name 'mtxName' from file 'outNameGammas'
920 //
921 TFile fileg(outNameGammas);
922
923 matrixg.Read(mtxName);
924 matrixg.Print("SizeCols");
925
926
927 //*****************************************************************
928
929 const char* mtxName = "MatrixHadrons";
930
931 gLog << "" << endl;
932 gLog << "========================================================" << endl;
933 gLog << " Get matrix for (hadrons)" << endl;
934 gLog << "matrix name = " << mtxName << endl;
935 gLog << "name of root file = " << outNameHadrons << endl;
936 gLog << "" << endl;
937
938
939 // read in the object with the name 'mtxName' from file 'outNameHadrons'
940 //
941 TFile fileh(outNameHadrons);
942
943 matrixh.Read(mtxName);
944 matrixh.Print("SizeCols");
945 }
946
947
948 //*************************************************************************
949 // create matrices of training events
950 // and root files of training and test events
951
952if (CTrainRF)
953 {
954 gLog << "" << endl;
955 gLog << "========================================================" << endl;
956 gLog << " Create matrices of training events and root files of training and test events"
957 << endl;
958 gLog << " Gammas :" << endl;
959 gLog << "---------" << endl;
960
961 MParList plistg;
962 MTaskList tlistg;
963
964 MFilterList flistgtrain;
965 flistgtrain.SetName("GammaFilterTrain");
966 MFilterList flistgtest;
967 flistgtest.SetName("GammaFilterTest");
968
969 MReadMarsFile readg("Events", filenameMC);
970 readg.DisableAutoScheme();
971
972
973 MFParticleId fgamma("MMcEvt", '=', kGAMMA);
974 fgamma.SetName("gammaID");
975
976 TString mgname("costhg");
977 MBinning bing("Binning"+mgname);
978 bing.SetEdges(10, 0., 1.0);
979
980 MH3 gref("cos(MMcEvt.fTelescopeTheta)");
981 gref.SetName(mgname);
982 MH::SetBinning(&gref.GetHist(), &bing);
983 for (Int_t i=1; i<=gref.GetNbins(); i++)
984 gref.GetHist().SetBinContent(i, 1.0);
985
986 MFEventSelector2 selectorg(gref);
987 selectorg.SetNumMax(howManyGammasTrain+howManyGammasTest);
988 selectorg.SetName("selectGammasTrainTest");
989
990 TString strgtrain = "gRandom->Rndm()*(";
991 strgtrain += howManyGammasTrain;
992 strgtrain += "+";
993 strgtrain += howManyGammasTest;
994 strgtrain += ")<";
995 strgtrain += howManyGammasTrain;
996 MF selsplitg(strgtrain);
997
998 MFillH fillmatg("MatrixGammas");
999 fillmatg.SetFilter(&flistgtrain);
1000 fillmatg.SetName("fillGammas");
1001
1002 //-----------------------
1003 // for writing the root files of training and test events
1004 // for gammas
1005
1006 MWriteRootFile writetraing(outNameGammasTrain, "RECREATE");
1007 writetraing.SetName("WriteGammasTrain");
1008 writetraing.SetFilter(&flistgtrain);
1009
1010 writetraing.AddContainer("MRawRunHeader", "RunHeaders");
1011 writetraing.AddContainer("MTime", "Events");
1012 writetraing.AddContainer("MMcEvt", "Events");
1013 writetraing.AddContainer("ThetaOrig", "Events");
1014 writetraing.AddContainer("MSrcPosCam", "Events");
1015 writetraing.AddContainer("MSigmabar", "Events");
1016 writetraing.AddContainer("MHillas", "Events");
1017 writetraing.AddContainer("MHillasExt", "Events");
1018 writetraing.AddContainer("MHillasSrc", "Events");
1019 writetraing.AddContainer("MNewImagePar", "Events");
1020
1021 MContinue contg(&flistgtrain);
1022
1023 MWriteRootFile writetestg(outNameGammasTest, "RECREATE");
1024 writetestg.SetName("WriteGammasTest");
1025 writetestg.SetFilter(&flistgtest);
1026
1027 writetestg.AddContainer("MRawRunHeader", "RunHeaders");
1028 writetestg.AddContainer("MTime", "Events");
1029 writetestg.AddContainer("MMcEvt", "Events");
1030 writetestg.AddContainer("ThetaOrig", "Events");
1031 writetestg.AddContainer("MSrcPosCam", "Events");
1032 writetestg.AddContainer("MSigmabar", "Events");
1033 writetestg.AddContainer("MHillas", "Events");
1034 writetestg.AddContainer("MHillasExt", "Events");
1035 writetestg.AddContainer("MHillasSrc", "Events");
1036 writetestg.AddContainer("MNewImagePar", "Events");
1037
1038 //-----------------------
1039
1040 //***************************** fill gammas ***
1041 // entries in MFilterList
1042
1043 flistgtrain.AddToList(&fgamma);
1044 flistgtrain.AddToList(&selectorg);
1045 flistgtrain.AddToList(&selsplitg);
1046
1047 flistgtest.AddToList(&fgamma);
1048 flistgtest.AddToList(&selectorg);
1049 flistgtest.AddToList(&selsplitg);
1050
1051 //*****************************
1052 // entries in MParList
1053
1054 plistg.AddToList(&tlistg);
1055 InitBinnings(&plistg);
1056
1057 plistg.AddToList(&matrixg);
1058
1059 //*****************************
1060 // entries in MTaskList
1061
1062 tlistg.AddToList(&readg);
1063
1064 tlistg.AddToList(&flistgtrain);
1065 tlistg.AddToList(&fillmatg);
1066 tlistg.AddToList(&writetraing);
1067
1068 tlistg.AddToList(&contg);
1069
1070 tlistg.AddToList(&flistgtest);
1071 tlistg.AddToList(&writetestg);
1072
1073 //*****************************
1074
1075 MProgressBar matrixbar;
1076 MEvtLoop evtloopg;
1077 evtloopg.SetName("FillGammaMatrix");
1078 evtloopg.SetParList(&plistg);
1079 //evtloopg.ReadEnv(env, "", printEnv);
1080 evtloopg.SetProgressBar(&matrixbar);
1081
1082 Int_t maxevents = -1;
1083 if (!evtloopg.Eventloop(maxevents))
1084 return;
1085
1086 tlistg.PrintStatistics(0, kTRUE);
1087
1088 matrixg.Print("SizeCols");
1089 Int_t generatedgTrain = matrixg.GetM().GetNrows();
1090 if ( fabs(generatedgTrain-howManyGammasTrain) >
1091 3.0*sqrt(howManyGammasTrain) )
1092 {
1093 gLog << "ONOFFCT1Analysis.C : no.of generated gamma training events ("
1094 << generatedgTrain << ") is incompatible with the no.of requested events ("
1095 << howManyGammasTrain << ")" << endl;
1096 }
1097
1098
1099 Int_t generatedgTest = writetestg.GetNumExecutions();
1100 if ( fabs(generatedgTest-howManyGammasTest) >
1101 3.0*sqrt(howManyGammasTest) )
1102 {
1103 gLog << "ONOFFCT1Analysis.C : no.of generated gamma test events ("
1104 << generatedgTest << ") is incompatible with the no.of requested events ("
1105 << howManyGammasTest << ")" << endl;
1106 }
1107
1108 //***************************** fill hadrons ***
1109
1110 gLog << " Hadrons :" << endl;
1111 gLog << "----------" << endl;
1112
1113 MParList plisth;
1114 MTaskList tlisth;
1115
1116 MFilterList flisthtrain;
1117 flisthtrain.SetName("HadronFilterTrain");
1118 MFilterList flisthtest;
1119 flisthtest.SetName("HadronFilterTest");
1120
1121 MReadMarsFile readh("Events", filenameHad);
1122 readh.DisableAutoScheme();
1123
1124 MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
1125 fhadrons.SetName("hadronID)");
1126
1127 TString mhname("costhh");
1128 MBinning binh("Binning"+mhname);
1129 binh.SetEdges(10, 0., 1.0);
1130
1131 MH3 href("cos(MMcEvt.fTelescopeTheta)");
1132 href.SetName(mhname);
1133 MH::SetBinning(&href.GetHist(), &binh);
1134 for (Int_t i=1; i<=href.GetNbins(); i++)
1135 href.GetHist().SetBinContent(i, 1.0);
1136
1137 MFEventSelector2 selectorh(href);
1138 selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest);
1139 selectorh.SetName("selectHadronsTrainTest");
1140
1141 TString strhtrain = "gRandom->Rndm()*(";
1142 strhtrain += howManyHadronsTrain;
1143 strhtrain += "+";
1144 strhtrain += howManyHadronsTest;
1145 strhtrain += ")<";
1146 strhtrain += howManyHadronsTrain;
1147 MF selsplith(strhtrain);
1148
1149 MFillH fillmath("MatrixHadrons");
1150 fillmath.SetFilter(&flisthtrain);
1151 fillmath.SetName("fillHadrons");
1152
1153 //-----------------------
1154 // for writing the root files of training and test events
1155 // for hadrons
1156
1157 MWriteRootFile writetrainh(outNameHadronsTrain, "RECREATE");
1158 writetrainh.SetName("WriteHadronsTrain");
1159 writetrainh.SetFilter(&flisthtrain);
1160
1161 writetrainh.AddContainer("MRawRunHeader", "RunHeaders");
1162 writetrainh.AddContainer("MTime", "Events");
1163 writetrainh.AddContainer("MMcEvt", "Events");
1164 writetrainh.AddContainer("ThetaOrig", "Events");
1165 writetrainh.AddContainer("MSrcPosCam", "Events");
1166 writetrainh.AddContainer("MSigmabar", "Events");
1167 writetrainh.AddContainer("MHillas", "Events");
1168 writetrainh.AddContainer("MHillasExt", "Events");
1169 writetrainh.AddContainer("MHillasSrc", "Events");
1170 writetrainh.AddContainer("MNewImagePar", "Events");
1171
1172 MContinue conth(&flisthtrain);
1173
1174 MWriteRootFile writetesth(outNameHadronsTest, "RECREATE");
1175 writetesth.SetName("WriteHadronsTest");
1176 writetesth.SetFilter(&flisthtest);
1177
1178 writetesth.AddContainer("MRawRunHeader", "RunHeaders");
1179 writetesth.AddContainer("MTime", "Events");
1180 writetesth.AddContainer("MMcEvt", "Events");
1181 writetesth.AddContainer("ThetaOrig", "Events");
1182 writetesth.AddContainer("MSrcPosCam", "Events");
1183 writetesth.AddContainer("MSigmabar", "Events");
1184 writetesth.AddContainer("MHillas", "Events");
1185 writetesth.AddContainer("MHillasExt", "Events");
1186 writetesth.AddContainer("MHillasSrc", "Events");
1187 writetesth.AddContainer("MNewImagePar", "Events");
1188
1189 //-----------------------
1190
1191 // entries in MFilterList
1192
1193 flisthtrain.AddToList(&fhadrons);
1194 flisthtrain.AddToList(&selectorh);
1195 flisthtrain.AddToList(&selsplith);
1196
1197 flisthtest.AddToList(&fhadrons);
1198 flisthtest.AddToList(&selectorh);
1199 flisthtest.AddToList(&selsplith);
1200
1201 //*****************************
1202 // entries in MParList
1203
1204 plisth.AddToList(&tlisth);
1205 InitBinnings(&plisth);
1206
1207 plisth.AddToList(&matrixh);
1208
1209 //*****************************
1210 // entries in MTaskList
1211
1212 tlisth.AddToList(&readh);
1213
1214 tlisth.AddToList(&flisthtrain);
1215 tlisth.AddToList(&fillmath);
1216 tlisth.AddToList(&writetrainh);
1217
1218 tlisth.AddToList(&conth);
1219
1220 tlisth.AddToList(&flisthtest);
1221 tlisth.AddToList(&writetesth);
1222
1223 //*****************************
1224
1225 MProgressBar matrixbar;
1226 MEvtLoop evtlooph;
1227 evtlooph.SetName("FillHadronMatrix");
1228 evtlooph.SetParList(&plisth);
1229 //evtlooph.ReadEnv(env, "", printEnv);
1230 evtlooph.SetProgressBar(&matrixbar);
1231
1232 Int_t maxevents = -1;
1233 if (!evtlooph.Eventloop(maxevents))
1234 return;
1235
1236 tlisth.PrintStatistics(0, kTRUE);
1237
1238 matrixh.Print("SizeCols");
1239 Int_t generatedhTrain = matrixh.GetM().GetNrows();
1240 if ( fabs(generatedhTrain-howManyHadronsTrain) >
1241 3.0*sqrt(howManyHadronsTrain) )
1242 {
1243 gLog << "ONOFFCT1Analysis.C : no.of generated hadron training events ("
1244 << generatedhTrain << ") is incompatible with the no.of requested events ("
1245 << howManyHadronsTrain << ")" << endl;
1246 }
1247
1248
1249 Int_t generatedhTest = writetesth.GetNumExecutions();
1250 if ( fabs(generatedhTest-howManyHadronsTest) >
1251 3.0*sqrt(howManyHadronsTest) )
1252 {
1253 gLog << "ONOFFCT1Analysis.C : no.of generated gamma test events ("
1254 << generatedhTest << ") is incompatible with the no.of requested events ("
1255 << howManyHadronsTest << ")" << endl;
1256 }
1257
1258
1259 //*****************************************************
1260
1261
1262 // write out matrices of training events
1263
1264 gLog << "" << endl;
1265 gLog << "========================================================" << endl;
1266 gLog << "Write out matrices of training events" << endl;
1267
1268
1269 //-------------------------------------------
1270 // "gammas"
1271 gLog << "Gammas :" << endl;
1272 matrixg.Print("SizeCols");
1273
1274 TFile writeg(outNameGammas, "RECREATE", "");
1275 matrixg.Write();
1276
1277 gLog << "" << endl;
1278 gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file "
1279 << outNameGammas << endl;
1280
1281 //-------------------------------------------
1282 // "hadrons"
1283 gLog << "Hadrons :" << endl;
1284 matrixh.Print("SizeCols");
1285
1286 TFile writeh(outNameHadrons, "RECREATE", "");
1287 matrixh.Write();
1288
1289 gLog << "" << endl;
1290 gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file "
1291 << outNameHadrons << endl;
1292
1293 }
1294 //********** end of creating matrices of training events ***********
1295
1296
1297 MRanForest *fRanForest;
1298 MRanTree *fRanTree;
1299 //-----------------------------------------------------------------
1300 // read in the trees of the random forest
1301 if (RTree)
1302 {
1303 MParList plisttr;
1304 MTaskList tlisttr;
1305 plisttr.AddToList(&tlisttr);
1306
1307 MReadTree readtr("TREE", outRF);
1308 readtr.DisableAutoScheme();
1309
1310 MRanForestFill rffill;
1311 rffill.SetNumTrees(NumTrees);
1312
1313 // list of tasks for the loop over the trees
1314
1315 tlisttr.AddToList(&readtr);
1316 tlisttr.AddToList(&rffill);
1317
1318 //-------------------
1319 // Execute tree loop
1320 //
1321 MEvtLoop evtlooptr;
1322 evtlooptr.SetName("ReadRFTrees");
1323 evtlooptr.SetParList(&plisttr);
1324 if (!evtlooptr.Eventloop())
1325 return;
1326
1327 tlisttr.PrintStatistics(0, kTRUE);
1328
1329
1330 // get adresses of objects which are used in the next eventloop
1331 fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
1332 if (!fRanForest)
1333 {
1334 gLog << err << dbginf << "MRanForest not found... aborting." << endl;
1335 return kFALSE;
1336 }
1337
1338 fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
1339 if (!fRanTree)
1340 {
1341 gLog << err << dbginf << "MRanTree not found... aborting." << endl;
1342 return kFALSE;
1343 }
1344
1345 }
1346
1347 //-----------------------------------------------------------------
1348 // grow the trees of the random forest (event loop = tree loop)
1349
1350 if (!RTree)
1351 {
1352
1353 gLog << "" << endl;
1354 gLog << "========================================================" << endl;
1355 gLog << "Macro CT1Analysis : start growing trees" << endl;
1356
1357 MTaskList tlist2;
1358 MParList plist2;
1359 plist2.AddToList(&tlist2);
1360
1361 plist2.AddToList(&matrixg);
1362 plist2.AddToList(&matrixh);
1363
1364 MRanForestGrow rfgrow2;
1365 rfgrow2.SetNumTrees(NumTrees);
1366 rfgrow2.SetNumTry(NumTry);
1367 rfgrow2.SetNdSize(NdSize);
1368
1369 MWriteRootFile rfwrite2(outRF);
1370 rfwrite2.AddContainer("MRanTree", "TREE");
1371
1372 MFillH fillh2("MHRanForestGini");
1373
1374 // list of tasks for the loop over the trees
1375
1376 tlist2.AddToList(&rfgrow2);
1377 tlist2.AddToList(&rfwrite2);
1378 tlist2.AddToList(&fillh2);
1379
1380 //-------------------
1381 // Execute tree loop
1382 //
1383 MEvtLoop treeloop;
1384 treeloop.SetName("GrowRFTrees");
1385 treeloop.SetParList(&plist2);
1386
1387 if ( !treeloop.Eventloop() )
1388 return;
1389
1390 tlist2.PrintStatistics(0, kTRUE);
1391
1392 plist2.FindObject("MHRanForestGini")->DrawClone();
1393
1394
1395 // get adresses of objects which are used in the next eventloop
1396 fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
1397 if (!fRanForest)
1398 {
1399 gLog << err << dbginf << "MRanForest not found... aborting." << endl;
1400 return kFALSE;
1401 }
1402
1403 fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
1404 if (!fRanTree)
1405 {
1406 gLog << err << dbginf << "MRanTree not found... aborting." << endl;
1407 return kFALSE;
1408 }
1409
1410 }
1411 // end of growing the trees of the random forest
1412 //-----------------------------------------------------------------
1413
1414
1415 //-----------------------------------------------------------------
1416 // Update the input files with the RF hadronness
1417 //
1418
1419 if (WRF)
1420 {
1421 gLog << "" << endl;
1422 gLog << "========================================================" << endl;
1423 gLog << "Update input file '" << filenameData
1424 << "' with the RF hadronness" << endl;
1425
1426 MTaskList tliston;
1427 MParList pliston;
1428
1429
1430 // geometry is needed in MHHillas... classes
1431 MGeomCam *fGeom =
1432 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1433
1434 //-------------------------------------------
1435 // create the tasks which should be executed
1436 //
1437
1438 MReadMarsFile read("Events", filenameData);
1439 read.DisableAutoScheme();
1440
1441
1442 //.......................................................................
1443 // calculate hadronnes for method of RANDOM FOREST
1444
1445
1446 MRanForestCalc rfcalc;
1447 rfcalc.SetHadronnessName(hadRFName);
1448
1449
1450 //.......................................................................
1451
1452 //MWriteRootFile write(outNameImage, "UPDATE");
1453 MWriteRootFile write(outNameImage, "RECREATE");
1454
1455 write.AddContainer("MRawRunHeader", "RunHeaders");
1456 write.AddContainer("MTime", "Events");
1457 write.AddContainer("MMcEvt", "Events");
1458 write.AddContainer("ThetaOrig", "Events");
1459 write.AddContainer("MSrcPosCam", "Events");
1460 write.AddContainer("MSigmabar", "Events");
1461 write.AddContainer("MHillas", "Events");
1462 write.AddContainer("MHillasExt", "Events");
1463 write.AddContainer("MHillasSrc", "Events");
1464 write.AddContainer("MNewImagePar", "Events");
1465
1466 write.AddContainer(hadRFName, "Events");
1467
1468 //-----------------------------------------------------------------
1469
1470
1471 MFCT1SelFinal selfinalgh(fHilNameSrc);
1472 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1473 selfinalgh.SetHadronnessName(hadRFName);
1474 selfinalgh.SetName("SelFinalgh");
1475 MContinue contfinalgh(&selfinalgh);
1476 contfinalgh.SetName("ContSelFinalgh");
1477
1478 MFillH fillranfor("MHRanForest");
1479 fillranfor.SetName("HRanForest");
1480
1481 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
1482 fillhadrf.SetName("HhadRF");
1483
1484 MFCT1SelFinal selfinal(fHilNameSrc);
1485 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1486 selfinal.SetHadronnessName(hadRFName);
1487 selfinal.SetName("SelFinal");
1488 MContinue contfinal(&selfinal);
1489 contfinal.SetName("ContSelFinal");
1490
1491 TString mh3name = "abs(Alpha)";
1492 MBinning binsalphaabs("Binning"+mh3name);
1493 binsalphaabs.SetEdges(50, -2.0, 98.0);
1494
1495 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
1496 alphaabs.SetName(mh3name);
1497 MFillH alpha(&alphaabs);
1498 alpha.SetName("FillAlphaAbs");
1499
1500
1501 MFillH hfill1("MHHillas", fHilName);
1502 hfill1.SetName("HHillas");
1503
1504 MFillH hfill2("MHStarMap", fHilName);
1505 hfill2.SetName("HStarMap");
1506
1507 MFillH hfill3("MHHillasExt", fHilNameSrc);
1508 hfill3.SetName("HHillasExt");
1509
1510 MFillH hfill4("MHHillasSrc", fHilNameSrc);
1511 hfill4.SetName("HHillasSrc");
1512
1513 MFillH hfill5("MHNewImagePar", fImgParName);
1514 hfill5.SetName("HNewImagePar");
1515
1516 //*****************************
1517 // entries in MParList
1518
1519 pliston.AddToList(&tliston);
1520 InitBinnings(&pliston);
1521
1522 pliston.AddToList(fRanForest);
1523 pliston.AddToList(fRanTree);
1524
1525 pliston.AddToList(&binsalphaabs);
1526 pliston.AddToList(&alphaabs);
1527
1528
1529 //*****************************
1530 // entries in MTaskList
1531
1532 tliston.AddToList(&read);
1533
1534 tliston.AddToList(&rfcalc);
1535 tliston.AddToList(&fillranfor);
1536 tliston.AddToList(&fillhadrf);
1537
1538 tliston.AddToList(&write);
1539 tliston.AddToList(&contfinalgh);
1540
1541 tliston.AddToList(&alpha);
1542 tliston.AddToList(&hfill1);
1543 tliston.AddToList(&hfill2);
1544 tliston.AddToList(&hfill3);
1545 tliston.AddToList(&hfill4);
1546 tliston.AddToList(&hfill5);
1547
1548 tliston.AddToList(&contfinal);
1549
1550 //*****************************
1551
1552 //-------------------------------------------
1553 // Execute event loop
1554 //
1555 MProgressBar bar;
1556 MEvtLoop evtloop;
1557 evtloop.SetName("UpdateRootFile");
1558 evtloop.SetParList(&pliston);
1559 evtloop.SetProgressBar(&bar);
1560
1561 Int_t maxevents = -1;
1562 if ( !evtloop.Eventloop(maxevents) )
1563 return;
1564
1565 tliston.PrintStatistics(0, kTRUE);
1566
1567
1568 //-------------------------------------------
1569 // Display the histograms
1570 //
1571 pliston.FindObject("MHRanForest")->DrawClone();
1572 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
1573
1574 pliston.FindObject("MHHillas")->DrawClone();
1575 pliston.FindObject("MHHillasExt")->DrawClone();
1576 pliston.FindObject("MHHillasSrc")->DrawClone();
1577 pliston.FindObject("MHNewImagePar")->DrawClone();
1578 pliston.FindObject("MHStarMap")->DrawClone();
1579
1580
1581 //-------------------------------------------
1582 // fit alpha distribution to get the number of excess events and
1583 // calculate significance of gamma signal in the alpha plot
1584
1585 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
1586 TH1 &alphaHist = absalpha->GetHist();
1587 alphaHist.SetXTitle("|alpha| [\\circ]");
1588 alphaHist.SetName("alpha-macro");
1589
1590 Double_t alphasig = 13.1;
1591 Double_t alphamin = 30.0;
1592 Double_t alphamax = 90.0;
1593 Int_t degree = 2;
1594 Double_t significance = -99.0;
1595 Bool_t drawpoly = kTRUE;
1596 Bool_t fitgauss = kTRUE;
1597 Bool_t print = kTRUE;
1598
1599 MHFindSignificance findsig;
1600 findsig.SetRebin(kTRUE);
1601 findsig.SetReduceDegree(kFALSE);
1602
1603 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
1604 alphasig, drawpoly, fitgauss, print);
1605 significance = findsig.GetSignificance();
1606 Float_t alphasi = findsig.GetAlphasi();
1607
1608 gLog << "For file '" << filenameData << "' : " << endl;
1609 gLog << "Significance of gamma signal after supercuts : "
1610 << significance << " (for |alpha| < " << alphasi << " degrees)"
1611 << endl;
1612
1613 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
1614
1615 //-------------------------------------------
1616
1617
1618 DeleteBinnings(&pliston);
1619 }
1620
1621 gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
1622 gLog << "=======================================================" << endl;
1623 }
1624 //---------------------------------------------------------------------
1625
1626
1627 //---------------------------------------------------------------------
1628 // Job B_SC_UP
1629 //============
1630
1631 // - create (or read in) optimum supercuts parameter values
1632 //
1633 // - calculate the hadroness for the supercuts
1634 //
1635 // - update input root file, including the hadroness
1636
1637
1638 if (JobB_SC_UP)
1639 {
1640 gLog << "=====================================================" << endl;
1641 gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl;
1642
1643 gLog << "" << endl;
1644 gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
1645 << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", "
1646 << (CMatrix ? "kTRUE" : "kFALSE") << ", "
1647 << (RMatrix ? "kTRUE" : "kFALSE") << ", "
1648 << (WOptimize ? "kTRUE" : "kFALSE") << ", "
1649 << (RTest ? "kTRUE" : "kFALSE") << ", "
1650 << (WSC ? "kTRUE" : "kFALSE") << endl;
1651
1652
1653 //--------------------------------------------
1654 // file which contains the initial parameter values for the supercuts
1655 // if parSCinit ="" the initial values are taken from the constructor of
1656 // MCT1Supercuts
1657
1658 TString parSCinit = outPath;
1659 //parSCinit += "parSC_1709d";
1660 parSCinit = "";
1661
1662 gLog << "parSCinit = " << parSCinit << endl;
1663
1664 //---------------
1665 // file onto which the optimal parameter values for the supercuts
1666 // are written
1667
1668 TString parSCfile = outPath;
1669 parSCfile += "parSC_2310a";
1670
1671 gLog << "parSCfile = " << parSCfile << endl;
1672
1673 //--------------------------------------------
1674 // file to be updated (either ON or MC)
1675
1676 //TString typeInput = "ON";
1677 //TString typeInput = "OFF";
1678 TString typeInput = "MC";
1679 gLog << "typeInput = " << typeInput << endl;
1680
1681 // name of input root file
1682 TString filenameData = outPath;
1683 filenameData += typeInput;
1684 filenameData += "2.root";
1685 gLog << "filenameData = " << filenameData << endl;
1686
1687 // name of output root file
1688 TString outNameImage = outPath;
1689 outNameImage += typeInput;
1690 outNameImage += "3.root";
1691
1692
1693 //TString outNameImage = filenameData;
1694
1695 gLog << "outNameImage = " << outNameImage << endl;
1696
1697 //--------------------------------------------
1698 // files to be read for optimizing the supercuts
1699 //
1700 // for the training
1701 TString filenameTrain = outPath;
1702 filenameTrain += "ON";
1703 filenameTrain += "1.root";
1704 Int_t howManyTrain = 800000;
1705 gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = "
1706 << howManyTrain << endl;
1707
1708 // for testing
1709 TString filenameTest = outPath;
1710 filenameTest += "ON";
1711 filenameTest += "1.root";
1712 Int_t howManyTest = 800000;
1713
1714 gLog << "filenameTest = " << filenameTest << ", howManyTest = "
1715 << howManyTest << endl;
1716
1717
1718 //--------------------------------------------
1719 // files to contain the matrices (generated from filenameTrain and
1720 // filenameTest)
1721 //
1722 // for the training
1723 TString fileMatrixTrain = outPath;
1724 fileMatrixTrain += "MatrixTrainSC";
1725 fileMatrixTrain += ".root";
1726 gLog << "fileMatrixTrain = " << fileMatrixTrain << endl;
1727
1728 // for testing
1729 TString fileMatrixTest = outPath;
1730 fileMatrixTest += "MatrixTestSC";
1731 fileMatrixTest += ".root";
1732 gLog << "fileMatrixTest = " << fileMatrixTest << endl;
1733
1734
1735
1736 //---------------------------------------------------------------------
1737 // Training and test matrices :
1738 // - either create them and write them onto a file
1739 // - or read them from a file
1740
1741
1742 MCT1FindSupercuts findsuper;
1743 findsuper.SetFilenameParam(parSCfile);
1744 findsuper.SetHadronnessName("HadSC");
1745
1746
1747 //--------------------------
1748 // create matrices and write them onto files
1749 if (CMatrix)
1750 {
1751 MH3 &mh3 = *(new MH3("MHillas.fSize"));
1752 mh3.SetName("Target distribution for SIZE");
1753
1754 if (filenameTrain == filenameTest)
1755 {
1756 if ( !findsuper.DefineTrainTestMatrix(filenameTrain,
1757 howManyTrain, mh3, howManyTest, mh3,
1758 fileMatrixTrain, fileMatrixTest) )
1759 {
1760 *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
1761 return;
1762 }
1763
1764 }
1765 else
1766 {
1767 if ( !findsuper.DefineTrainMatrix(filenameTrain,
1768 howManyTrain, mh3, fileMatrixTrain) )
1769 {
1770 *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
1771 return;
1772 }
1773
1774 if ( !findsuper.DefineTestMatrix( filenameTest,
1775 howManyTest, mh3, fileMatrixTest) )
1776 {
1777 *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
1778 return;
1779 }
1780 }
1781 }
1782
1783 //--------------------------
1784 // read matrices from files
1785 //
1786
1787 if (RMatrix)
1788 findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
1789 //--------------------------
1790
1791
1792
1793 //---------------------------------------------------------------------
1794 // optimize supercuts using the training sample
1795 //
1796 // the initial values are taken
1797 // - from the file parSCinit (if != "")
1798 // - or from the arrays params and steps (if their sizes are != 0)
1799 // - or from the MCT1Supercuts constructor
1800
1801
1802if (WOptimize)
1803 {
1804 gLog << "CT1Analysis.C : optimize the supercuts using the training matrix"
1805 << endl;
1806
1807 TArrayD params(0);
1808 TArrayD steps(0);
1809
1810 if (parSCinit == "")
1811 {
1812 Double_t vparams[104] = {
1813 // LengthUp
1814 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
1815 0.007388, -0.013463,
1816 // LengthLo
1817 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993,
1818 0.000080, -0.007157,
1819 // WidthUp
1820 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353,
1821 0.020711, -0.016703,
1822 // WidthLo
1823 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
1824 0.006126, -0.002849,
1825 // DistUp
1826 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909,
1827 0.189697, 0.0,
1828 // DistLo
1829 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374,
1830 -0.001750, 0.0,
1831 // AsymUp
1832 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1833 0.0, 0.0,
1834 // AsymLo
1835 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1836 0.0, 0.0,
1837 // ConcUp
1838 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1839 0.0, 0.0,
1840 // ConcLo
1841 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1842 0.0, 0.0,
1843 // Leakage1Up
1844 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1845 0.0, 0.0,
1846 // Leakage1Lo
1847 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1848 0.0, 0.0,
1849 // AlphaUp
1850 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0,
1851 0.0, 0.0 };
1852
1853 Double_t vsteps[104] = {
1854 // LengthUp
1855 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002,
1856 0.0008, 0.002,
1857 // LengthLo
1858 0.02, 0.003, 0.05, 0.006, 0.002, 0.3,
1859 0.0001, 0.0008,
1860 // WidthUp
1861 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006,
1862 0.002, 0.002,
1863 // WidthLo
1864 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002,
1865 0.0007, 0.003,
1866 // DistUp
1867 0.2, 0.0, 0.3, 0.02, 0.0, 0.03,
1868 0.02, 0.0
1869 // DistLo
1870 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005,
1871 0.0002, 0.0
1872 // AsymUp
1873 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1874 0.0, 0.0,
1875 // AsymLo
1876 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1877 0.0, 0.0,
1878 // ConcUp
1879 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1880 0.0, 0.0,
1881 // ConcLo
1882 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1883 0.0, 0.0,
1884 // Leakage1Up
1885 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1886 0.0, 0.0,
1887 // Leakage1Lo
1888 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1889 0.0, 0.0,
1890 // AlphaUp
1891 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1892 0.0, 0.0 };
1893
1894 params.Set(104, vparams);
1895 steps.Set (104, vsteps );
1896 }
1897
1898 Bool_t rf;
1899 rf = findsuper.FindParams(parSCinit, params, steps);
1900
1901 if (!rf)
1902 {
1903 gLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
1904 return;
1905 }
1906 }
1907
1908 //--------------------------------------
1909 // test the supercuts on the test sample
1910 //
1911
1912 if (RTest)
1913 {
1914 gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl;
1915 Bool_t rt = findsuper.TestParams();
1916 if (!rt)
1917 {
1918 gLog << "CT1Analysis.C : test of supercuts on the test matrix failed"
1919 << endl;
1920 }
1921
1922 }
1923
1924
1925 //-----------------------------------------------------------------
1926 // Update the input files with the SC hadronness
1927 //
1928
1929 if (WSC)
1930 {
1931 gLog << "" << endl;
1932 gLog << "========================================================" << endl;
1933 gLog << "Update input file '" << filenameData
1934 << "' with the SC hadronness" << endl;
1935
1936
1937 //----------------------------------------------------
1938 // read in optimum parameter values for the supercuts
1939
1940 TFile inparam(parSCfile);
1941 MCT1Supercuts scin;
1942 scin.Read("MCT1Supercuts");
1943 inparam.Close();
1944
1945 gLog << "Parameter values for supercuts were read in from file '"
1946 << parSCfile << "'" << endl;
1947
1948 TArrayD supercutsPar;
1949 supercutsPar = scin.GetParameters();
1950
1951 TArrayD supercutsStep;
1952 supercutsStep = scin.GetStepsizes();
1953
1954 gLog << "Parameter values for supercuts : " << endl;
1955 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
1956 {
1957 gLog << supercutsPar[i] << ", ";
1958 }
1959 gLog << endl;
1960
1961 gLog << "Step values for supercuts : " << endl;
1962 for (Int_t i=0; i<supercutsStep.GetSize(); i++)
1963 {
1964 gLog << supercutsStep[i] << ", ";
1965 }
1966 gLog << endl;
1967
1968
1969 //----------------------------------------------------
1970 MTaskList tliston;
1971 MParList pliston;
1972
1973 // set the parameters of the supercuts
1974 MCT1Supercuts supercuts;
1975 supercuts.SetParameters(supercutsPar);
1976 gLog << "parameter values for the supercuts used for updating the input file ' "
1977 << filenameData << "'" << endl;
1978 supercutsPar = supercuts.GetParameters();
1979 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
1980 {
1981 gLog << supercutsPar[i] << ", ";
1982 }
1983 gLog << endl;
1984
1985
1986 // geometry is needed in MHHillas... classes
1987 MGeomCam *fGeom =
1988 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1989
1990 //-------------------------------------------
1991 // create the tasks which should be executed
1992 //
1993
1994 MReadMarsFile read("Events", filenameData);
1995 read.DisableAutoScheme();
1996
1997 TString fHilName = "MHillas";
1998 TString fHilNameExt = "MHillasExt";
1999 TString fHilNameSrc = "MHillasSrc";
2000 TString fImgParName = "MNewImagePar";
2001
2002
2003 //.......................................................................
2004 // calculation of hadroness for the supercuts
2005 // (=0.25 if fullfilled, =0.75 otherwise)
2006
2007 TString hadSCName = "HadSC";
2008 MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
2009 sccalc.SetHadronnessName(hadSCName);
2010
2011
2012 //.......................................................................
2013
2014
2015 //MWriteRootFile write(outNameImage, "UPDATE");
2016 //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
2017
2018
2019 MWriteRootFile write(outNameImage, "RECREATE");
2020
2021 write.AddContainer("MRawRunHeader", "RunHeaders");
2022 write.AddContainer("MTime", "Events");
2023 write.AddContainer("MMcEvt", "Events");
2024 write.AddContainer("ThetaOrig", "Events");
2025 write.AddContainer("MSrcPosCam", "Events");
2026 write.AddContainer("MSigmabar", "Events");
2027 write.AddContainer("MHillas", "Events");
2028 write.AddContainer("MHillasExt", "Events");
2029 write.AddContainer("MHillasSrc", "Events");
2030 write.AddContainer("MNewImagePar", "Events");
2031
2032 write.AddContainer("HadRF", "Events");
2033 write.AddContainer(hadSCName, "Events");
2034
2035
2036 //-----------------------------------------------------------------
2037 // geometry is needed in MHHillas... classes
2038 MGeomCam *fGeom =
2039 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2040
2041 Float_t maxhadronness = 0.40;
2042 Float_t maxalpha = 20.0;
2043 Float_t maxdist = 10.0;
2044
2045 MFCT1SelFinal selfinalgh(fHilNameSrc);
2046 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2047 selfinalgh.SetHadronnessName(hadSCName);
2048 selfinalgh.SetName("SelFinalgh");
2049 MContinue contfinalgh(&selfinalgh);
2050 contfinalgh.SetName("ContSelFinalgh");
2051
2052 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
2053 fillhadsc.SetName("HhadSC");
2054
2055 MFCT1SelFinal selfinal(fHilNameSrc);
2056 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2057 selfinal.SetHadronnessName(hadSCName);
2058 selfinal.SetName("SelFinal");
2059 MContinue contfinal(&selfinal);
2060 contfinal.SetName("ContSelFinal");
2061
2062 TString mh3name = "abs(Alpha)";
2063 MBinning binsalphaabs("Binning"+mh3name);
2064 binsalphaabs.SetEdges(50, -2.0, 98.0);
2065
2066 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
2067 alphaabs.SetName(mh3name);
2068
2069 TH1 &alphahist = alphaabs->GetHist();
2070
2071 MFillH alpha(&alphaabs);
2072 alpha.SetName("FillAlphaAbs");
2073
2074 MFillH hfill1("MHHillas", fHilName);
2075 hfill1.SetName("HHillas");
2076
2077 MFillH hfill2("MHStarMap", fHilName);
2078 hfill2.SetName("HStarMap");
2079
2080 MFillH hfill3("MHHillasExt", fHilNameSrc);
2081 hfill3.SetName("HHillasExt");
2082
2083 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2084 hfill4.SetName("HHillasSrc");
2085
2086 MFillH hfill5("MHNewImagePar", fImgParName);
2087 hfill5.SetName("HNewImagePar");
2088
2089 //*****************************
2090 // entries in MParList
2091
2092 pliston.AddToList(&tliston);
2093 InitBinnings(&pliston);
2094
2095 pliston.AddToList(&supercuts);
2096
2097 pliston.AddToList(&binsalphaabs);
2098 pliston.AddToList(&alphaabs);
2099
2100 //*****************************
2101 // entries in MTaskList
2102
2103 tliston.AddToList(&read);
2104
2105 tliston.AddToList(&sccalc);
2106 tliston.AddToList(&fillhadsc);
2107
2108 tliston.AddToList(&write);
2109 tliston.AddToList(&contfinalgh);
2110
2111 tliston.AddToList(&alpha);
2112 tliston.AddToList(&hfill1);
2113 tliston.AddToList(&hfill2);
2114 tliston.AddToList(&hfill3);
2115 tliston.AddToList(&hfill4);
2116 tliston.AddToList(&hfill5);
2117
2118 tliston.AddToList(&contfinal);
2119
2120 //*****************************
2121
2122 //-------------------------------------------
2123 // Execute event loop
2124 //
2125 MProgressBar bar;
2126 MEvtLoop evtloop;
2127 evtloop.SetParList(&pliston);
2128 evtloop.SetProgressBar(&bar);
2129
2130 Int_t maxevents = -1;
2131 if ( !evtloop.Eventloop(maxevents) )
2132 return;
2133
2134 tliston.PrintStatistics(0, kTRUE);
2135
2136
2137 //-------------------------------------------
2138 // Display the histograms
2139 //
2140 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2141
2142 pliston.FindObject("MHHillas")->DrawClone();
2143 pliston.FindObject("MHHillasExt")->DrawClone();
2144 pliston.FindObject("MHHillasSrc")->DrawClone();
2145 pliston.FindObject("MHNewImagePar")->DrawClone();
2146 pliston.FindObject("MHStarMap")->DrawClone();
2147
2148 //-------------------------------------------
2149 // fit alpha distribution to get the number of excess events and
2150 // calculate significance of gamma signal in the alpha plot
2151
2152 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
2153 TH1 &alphaHist = absalpha->GetHist();
2154 alphaHist.SetXTitle("|alpha| [\\circ]");
2155 alphaHist.SetName("alpha-macro");
2156
2157 Double_t alphasig = 13.1;
2158 Double_t alphamin = 30.0;
2159 Double_t alphamax = 90.0;
2160 Int_t degree = 2;
2161 Double_t significance = -99.0;
2162 Bool_t drawpoly = kTRUE;
2163 Bool_t fitgauss = kTRUE;
2164 Bool_t print = kTRUE;
2165
2166 MHFindSignificance findsig;
2167 findsig.SetRebin(kTRUE);
2168 findsig.SetReduceDegree(kFALSE);
2169
2170 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
2171 alphasig, drawpoly, fitgauss, print);
2172 significance = findsig.GetSignificance();
2173 Float_t alphasi = findsig.GetAlphasi();
2174
2175 gLog << "For file '" << filenameData << "' : " << endl;
2176 gLog << "Significance of gamma signal after supercuts : "
2177 << significance << " (for |alpha| < " << alphasi << " degrees)"
2178 << endl;
2179
2180 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
2181
2182 //-------------------------------------------
2183
2184 DeleteBinnings(&pliston);
2185 }
2186
2187
2188 gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
2189 gLog << "=======================================================" << endl;
2190 }
2191 //---------------------------------------------------------------------
2192
2193
2194
2195 //---------------------------------------------------------------------
2196 // Job C
2197 //======
2198
2199 // - read ON1 and MC1 data files
2200 // which should have been updated to contain the hadronnesses
2201 // for the method of Random Forest and for the SUPERCUTS
2202 // - produce Neyman-Pearson plots
2203
2204 if (JobC)
2205 {
2206 gLog << "=====================================================" << endl;
2207 gLog << "Macro CT1Analysis : Start of Job C" << endl;
2208
2209 gLog << "" << endl;
2210 gLog << "Macro CT1Analysis : JobC = "
2211 << (JobC ? "kTRUE" : "kFALSE") << endl;
2212
2213
2214 // name of input data file
2215 TString filenameData = outPath;
2216 filenameData += "OFF";
2217 filenameData += "3.root";
2218 gLog << "filenameData = " << filenameData << endl;
2219
2220 // name of input MC file
2221 TString filenameMC = outPath;
2222 filenameMC += "MC";
2223 filenameMC += "3.root";
2224 gLog << "filenameMC = " << filenameMC << endl;
2225
2226
2227 //-----------------------------------------------------------------
2228
2229 MTaskList tliston;
2230 MParList pliston;
2231
2232
2233 // geometry is needed in MHHillas... classes
2234 MGeomCam *fGeom =
2235 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2236
2237 //-------------------------------------------
2238 // create the tasks which should be executed
2239 //
2240
2241 MReadMarsFile read("Events", filenameMC);
2242 read.AddFile(filenameData);
2243 read.DisableAutoScheme();
2244
2245
2246 //.......................................................................
2247 // names of hadronness containers
2248
2249 TString hadSCName = "HadSC";
2250 TString hadRFName = "HadRF";
2251
2252 //.......................................................................
2253
2254
2255 TString fHilName = "MHillas";
2256 TString fHilNameExt = "MHillasExt";
2257 TString fHilNameSrc = "MHillasSrc";
2258 TString fImgParName = "MNewImagePar";
2259
2260 Float_t maxhadronness = 0.40;
2261 Float_t maxalpha = 20.0;
2262 Float_t maxdist = 10.0;
2263
2264 MFCT1SelFinal selfinalgh(fHilNameSrc);
2265 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2266 selfinalgh.SetHadronnessName(hadRFName);
2267 selfinalgh.SetName("SelFinalgh");
2268 MContinue contfinalgh(&selfinalgh);
2269 contfinalgh.SetName("ContSelFinalgh");
2270
2271
2272 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
2273 fillhadsc.SetName("HhadSC");
2274 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
2275 fillhadrf.SetName("HhadRF");
2276
2277 MFCT1SelFinal selfinal(fHilNameSrc);
2278 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2279 selfinal.SetHadronnessName(hadRFName);
2280 selfinal.SetName("SelFinal");
2281 MContinue contfinal(&selfinal);
2282 contfinal.SetName("ContSelFinal");
2283
2284
2285 MFillH hfill1("MHHillas", fHilName);
2286 hfill1.SetName("HHillas");
2287
2288 MFillH hfill2("MHStarMap", fHilName);
2289 hfill2.SetName("HStarMap");
2290
2291 MFillH hfill3("MHHillasExt", fHilNameSrc);
2292 hfill3.SetName("HHillasExt");
2293
2294 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2295 hfill4.SetName("HHillasSrc");
2296
2297 MFillH hfill5("MHNewImagePar", fImgParName);
2298 hfill5.SetName("HNewImagePar");
2299
2300
2301 //*****************************
2302 // entries in MParList
2303
2304 pliston.AddToList(&tliston);
2305 InitBinnings(&pliston);
2306
2307
2308 //*****************************
2309 // entries in MTaskList
2310
2311 tliston.AddToList(&read);
2312
2313 tliston.AddToList(&fillhadsc);
2314 tliston.AddToList(&fillhadrf);
2315
2316 tliston.AddToList(&contfinalgh);
2317 tliston.AddToList(&hfill1);
2318 tliston.AddToList(&hfill2);
2319 tliston.AddToList(&hfill3);
2320 tliston.AddToList(&hfill4);
2321 tliston.AddToList(&hfill5);
2322
2323 tliston.AddToList(&contfinal);
2324
2325 //*****************************
2326
2327 //-------------------------------------------
2328 // Execute event loop
2329 //
2330 MProgressBar bar;
2331 MEvtLoop evtloop;
2332 evtloop.SetParList(&pliston);
2333 evtloop.SetProgressBar(&bar);
2334
2335 Int_t maxevents = -1;
2336 //Int_t maxevents = 35000;
2337 if ( !evtloop.Eventloop(maxevents) )
2338 return;
2339
2340 tliston.PrintStatistics(0, kTRUE);
2341
2342
2343 //-------------------------------------------
2344 // Display the histograms
2345 //
2346
2347 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2348 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2349
2350 pliston.FindObject("MHHillas")->DrawClone();
2351 pliston.FindObject("MHHillasExt")->DrawClone();
2352 pliston.FindObject("MHHillasSrc")->DrawClone();
2353 pliston.FindObject("MHNewImagePar")->DrawClone();
2354 pliston.FindObject("MHStarMap")->DrawClone();
2355
2356 DeleteBinnings(&pliston);
2357
2358 gLog << "Macro CT1Analysis : End of Job C" << endl;
2359 gLog << "===================================================" << endl;
2360 }
2361
2362
2363
2364 //---------------------------------------------------------------------
2365 // Job D
2366 //======
2367
2368 // - select g/h separation method XX
2369 // - read ON2 (or MC2) root file
2370 // - apply cuts in hadronness
2371 // - make plots
2372
2373
2374 if (JobD)
2375 {
2376 gLog << "=====================================================" << endl;
2377 gLog << "Macro CT1Analysis : Start of Job D" << endl;
2378
2379 gLog << "" << endl;
2380 gLog << "Macro CT1Analysis : JobD = "
2381 << (JobD ? "kTRUE" : "kFALSE") << endl;
2382
2383
2384 // type of data to be analysed
2385 TString typeData = "ON";
2386 //TString typeData = "OFF";
2387 //TString typeData = "MC";
2388 gLog << "typeData = " << typeData << endl;
2389
2390 TString ext = "3.root";
2391
2392
2393 //------------------------------
2394 // selection of g/h separation method
2395 // and definition of final selections
2396
2397 //TString XX("SC");
2398 TString XX("RF");
2399 TString fhadronnessName("Had");
2400 fhadronnessName += XX;
2401 gLog << "fhadronnessName = " << fhadronnessName << endl;
2402
2403 // maximum values of the hadronness, |ALPHA| and DIST
2404 Float_t maxhadronness = 0.233;
2405 Float_t maxalpha = 20.0;
2406 Float_t maxdist = 10.0;
2407 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2408 << maxhadronness << ", " << maxalpha << ", "
2409 << maxdist << endl;
2410
2411
2412 //------------------------------
2413 // name of data file to be analysed
2414 TString filenameData(outPath);
2415 filenameData += typeData;
2416 filenameData += ext;
2417 gLog << "filenameData = " << filenameData << endl;
2418
2419
2420
2421 //*************************************************************************
2422 //
2423 // Analyse the data
2424 //
2425
2426 MTaskList tliston;
2427 MParList pliston;
2428
2429 // geometry is needed in MHHillas... classes
2430 MGeomCam *fGeom =
2431 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2432
2433
2434 TString fHilName = "MHillas";
2435 TString fHilNameExt = "MHillasExt";
2436 TString fHilNameSrc = "MHillasSrc";
2437 TString fImgParName = "MNewImagePar";
2438
2439 //-------------------------------------------
2440 // create the tasks which should be executed
2441 //
2442
2443 MReadMarsFile read("Events", filenameData);
2444 read.DisableAutoScheme();
2445
2446
2447 //-----------------------------------------------------------------
2448 // geometry is needed in MHHillas... classes
2449 MGeomCam *fGeom =
2450 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2451
2452 MFCT1SelFinal selfinalgh(fHilNameSrc);
2453 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2454 selfinalgh.SetHadronnessName(fhadronnessName);
2455 selfinalgh.SetName("SelFinalgh");
2456 MContinue contfinalgh(&selfinalgh);
2457 contfinalgh.SetName("ContSelFinalgh");
2458
2459 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
2460 fillhadsc.SetName("HhadSC");
2461 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
2462 fillhadrf.SetName("HhadRF");
2463
2464 TString mh3name = "abs(Alpha)";
2465 MBinning binsalphaabs("Binning"+mh3name);
2466 binsalphaabs.SetEdges(50, -2.0, 98.0);
2467
2468 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
2469 alphaabs.SetName(mh3name);
2470
2471 TH1 &alphahist = alphaabs->GetHist();
2472
2473 MFillH alpha(&alphaabs);
2474 alpha.SetName("FillAlphaAbs");
2475
2476 MFillH hfill1("MHHillas", fHilName);
2477 hfill1.SetName("HHillas");
2478
2479 MFillH hfill2("MHStarMap", fHilName);
2480 hfill2.SetName("HStarMap");
2481
2482 MFillH hfill3("MHHillasExt", fHilNameSrc);
2483 hfill3.SetName("HHillasExt");
2484
2485 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2486 hfill4.SetName("HHillasSrc");
2487
2488 MFillH hfill5("MHNewImagePar", fImgParName);
2489 hfill5.SetName("HNewImagePar");
2490
2491 MFCT1SelFinal selfinal(fHilNameSrc);
2492 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2493 selfinal.SetHadronnessName(fhadronnessName);
2494 selfinal.SetName("SelFinal");
2495 MContinue contfinal(&selfinal);
2496 contfinal.SetName("ContSelFinal");
2497
2498
2499 //*****************************
2500 // entries in MParList
2501
2502 pliston.AddToList(&tliston);
2503 InitBinnings(&pliston);
2504 pliston.AddToList(&binsalphaabs);
2505 pliston.AddToList(&alphaabs);
2506
2507 //*****************************
2508 // entries in MTaskList
2509
2510 tliston.AddToList(&read);
2511
2512 tliston.AddToList(&contfinalgh);
2513
2514 tliston.AddToList(&fillhadsc);
2515 tliston.AddToList(&fillhadrf);
2516
2517 tliston.AddToList(&alpha);
2518 tliston.AddToList(&hfill1);
2519 tliston.AddToList(&hfill2);
2520 tliston.AddToList(&hfill3);
2521 tliston.AddToList(&hfill4);
2522 tliston.AddToList(&hfill5);
2523
2524 tliston.AddToList(&contfinal);
2525
2526 //*****************************
2527
2528 //-------------------------------------------
2529 // Execute event loop
2530 //
2531 MProgressBar bar;
2532 MEvtLoop evtloop;
2533 evtloop.SetParList(&pliston);
2534 evtloop.SetProgressBar(&bar);
2535
2536 Int_t maxevents = -1;
2537 //Int_t maxevents = 10000;
2538 if ( !evtloop.Eventloop(maxevents) )
2539 return;
2540
2541 tliston.PrintStatistics(0, kTRUE);
2542
2543
2544 //-------------------------------------------
2545 // Display the histograms
2546 //
2547
2548 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2549 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2550
2551 pliston.FindObject("MHHillas")->DrawClone();
2552 pliston.FindObject("MHHillasExt")->DrawClone();
2553 pliston.FindObject("MHHillasSrc")->DrawClone();
2554 pliston.FindObject("MHNewImagePar")->DrawClone();
2555 pliston.FindObject("MHStarMap")->DrawClone();
2556
2557
2558 //-------------------------------------------
2559
2560 // fit alpha distribution to get the number of excess events and
2561 // calculate significance of gamma signal in the alpha plot
2562
2563 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
2564 TH1 &alphaHist = absalpha->GetHist();
2565 alphaHist.SetXTitle("|alpha| [\\circ]");
2566 alphaHist.SetName("alpha-JobD");
2567
2568 Double_t alphasig = 13.1;
2569 Double_t alphamin = 30.0;
2570 Double_t alphamax = 90.0;
2571 Int_t degree = 2;
2572 Double_t significance = -99.0;
2573 Bool_t drawpoly = kTRUE;
2574 Bool_t fitgauss = kTRUE;
2575 Bool_t print = kTRUE;
2576
2577 MHFindSignificance findsig;
2578 findsig.SetRebin(kTRUE);
2579 findsig.SetReduceDegree(kFALSE);
2580
2581 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
2582 alphasig, drawpoly, fitgauss, print);
2583 significance = findsig.GetSignificance();
2584 Float_t alphasi = findsig.GetAlphasi();
2585
2586 gLog << "For file '" << filenameData << "' : " << endl;
2587 gLog << "Significance of gamma signal after supercuts : "
2588 << significance << " (for |alpha| < " << alphasi << " degrees)"
2589 << endl;
2590
2591 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
2592
2593 //-------------------------------------------
2594
2595
2596 DeleteBinnings(&pliston);
2597
2598 gLog << "Macro CT1Analysis : End of Job D" << endl;
2599 gLog << "=======================================================" << endl;
2600 }
2601 //---------------------------------------------------------------------
2602
2603
2604
2605
2606
2607 //---------------------------------------------------------------------
2608 // Job E_XX
2609 //=========
2610
2611 // - select g/h separation method XX
2612 // - read MC_XX2.root file
2613   // - calculate eff. collection area
2614 // - read ON_XX2.root file
2615 // - apply final cuts
2616 // - calculate flux
2617 // - write root file for ON data after final cuts (ON_XX3.root))
2618
2619
2620 if (JobE_XX)
2621 {
2622 gLog << "=====================================================" << endl;
2623 gLog << "Macro CT1Analysis : Start of Job E_XX" << endl;
2624
2625 gLog << "" << endl;
2626 gLog << "Macro CT1Analysis : JobE_XX, CCollArea, OEEst, WEX = "
2627 << (JobE_XX ? "kTRUE" : "kFALSE") << ", "
2628 << (CCollArea?"kTRUE" : "kFALSE") << ", "
2629 << (OEEst ? "kTRUE" : "kFALSE") << ", "
2630 << (WEX ? "kTRUE" : "kFALSE") << endl;
2631
2632
2633 // type of data to be analysed
2634 //TString typeData = "ON";
2635 //TString typeData = "OFF";
2636 TString typeData = "MC";
2637 gLog << "typeData = " << typeData << endl;
2638
2639 TString typeMC = "MC";
2640 TString ext = "3.root";
2641 TString extout = "4.root";
2642
2643 //------------------------------
2644 // selection of g/h separation method
2645 // and definition of final selections
2646
2647 //TString XX("SC");
2648 TString XX("RF");
2649 TString fhadronnessName("Had");
2650 fhadronnessName += XX;
2651 gLog << "fhadronnessName = " << fhadronnessName << endl;
2652
2653 // maximum values of the hadronness, |ALPHA| and DIST
2654 Float_t maxhadronness = 0.23;
2655 Float_t maxalpha = 20.0;
2656 Float_t maxdist = 10.0;
2657 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2658 << maxhadronness << ", " << maxalpha << ", "
2659 << maxdist << endl;
2660
2661 //------------------------------
2662 // name of MC file to be used for optimizing the energy estimator
2663 TString filenameOpt(outPath);
2664 filenameOpt += typeMC;
2665 filenameOpt += ext;
2666 gLog << "filenameOpt = " << filenameOpt << endl;
2667
2668 //------------------------------
2669 // name of file containing the parameters of the energy estimator
2670 TString energyParName(outPath);
2671 energyParName += "energyest_";
2672 energyParName += XX;
2673 energyParName += ".root";
2674 gLog << "energyParName = " << energyParName << endl;
2675
2676 //------------------------------
2677 // name of MC file to be used for calculating the eff. collection areas
2678 TString filenameArea(outPath);
2679 filenameArea += typeMC;
2680 filenameArea += ext;
2681 gLog << "filenameArea = " << filenameArea << endl;
2682
2683 //------------------------------
2684 // name of file containing the eff. collection areas
2685 TString collareaName(outPath);
2686 collareaName += "area_";
2687 collareaName += XX;
2688 collareaName += ".root";
2689 gLog << "collareaName = " << collareaName << endl;
2690
2691 //------------------------------
2692 // name of data file to be analysed
2693 TString filenameData(outPath);
2694 filenameData += typeData;
2695 filenameData += ext;
2696 gLog << "filenameData = " << filenameData << endl;
2697
2698 //------------------------------
2699 // name of output data file (after the final cuts)
2700 TString filenameDataout(outPath);
2701 filenameDataout += typeData;
2702 filenameDataout += "_";
2703 filenameDataout += XX;
2704 filenameDataout += extout;
2705 gLog << "filenameDataout = " << filenameDataout << endl;
2706
2707 //------------------------------
2708 // name of file containing histograms for flux calculastion
2709 TString filenameResults(outPath);
2710 filenameResults += typeData;
2711 filenameResults += "Results_";
2712 filenameResults += XX;
2713 filenameResults += extout;
2714 gLog << "filenameResults = " << filenameResults << endl;
2715
2716
2717 //====================================================================
2718
2719 MHMcCT1CollectionArea collarea;
2720 collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
2721
2722 MParList parlist;
2723 InitBinnings(&parlist);
2724
2725 if (CCollArea)
2726 {
2727 gLog << "-----------------------------------------------" << endl;
2728 gLog << "Start calculation of effective collection areas" << endl;
2729
2730
2731 MTaskList tasklist;
2732
2733 //---------------------------------------
2734 // Setup the tasks to be executed
2735 //
2736 MReadMarsFile reader("Events", filenameArea);
2737 reader.DisableAutoScheme();
2738
2739 MFCT1SelFinal cuthadrons;
2740 cuthadrons.SetHadronnessName(fhadronnessName);
2741 cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
2742
2743 MContinue conthadrons(&cuthadrons);
2744
2745
2746 MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
2747 filler.SetName("CollectionArea");
2748
2749 //********************************
2750 // entries in MParList
2751
2752 parlist.AddToList(&tasklist);
2753
2754 parlist.AddToList(&collarea);
2755
2756 //********************************
2757 // entries in MTaskList
2758
2759 tasklist.AddToList(&reader);
2760 tasklist.AddToList(&conthadrons);
2761 tasklist.AddToList(&filler);
2762
2763 //********************************
2764
2765 //-----------------------------------------
2766 // Execute event loop
2767 //
2768 MEvtLoop magic;
2769 magic.SetParList(&parlist);
2770
2771 MProgressBar bar;
2772 magic.SetProgressBar(&bar);
2773 if (!magic.Eventloop())
2774 return;
2775
2776 tasklist.PrintStatistics(0, kTRUE);
2777
2778 // Calculate effective collection areas
2779 // and display the histograms
2780 //
2781 //MHMcCT1CollectionArea *collarea =
2782 // (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
2783 collarea.CalcEfficiency();
2784 collarea.DrawClone();
2785
2786
2787
2788 //---------------------------------------------
2789 // Write histograms to a file
2790 //
2791
2792 TFile f(collareaName, "RECREATE");
2793 //collarea.GetHist()->Write();
2794 //collarea.GetHAll()->Write();
2795 //collarea.GetHSel()->Write();
2796 collarea.Write();
2797
2798 f.Close();
2799
2800 gLog << "Collection area plots written onto file " << collareaName << endl;
2801
2802 gLog << "Calculation of effective collection areas done" << endl;
2803 gLog << "-----------------------------------------------" << endl;
2804 //------------------------------------------------------------------
2805 }
2806
2807 if (!CCollArea)
2808 {
2809 gLog << "-----------------------------------------------" << endl;
2810 gLog << "Read in effective collection areas from file "
2811 << collareaName << endl;
2812
2813 TFile collfile(collareaName);
2814 collfile.ls();
2815 collarea.Read("MHMcCT1CollectionArea");
2816 collarea.DrawClone();
2817
2818 gLog << "Effective collection areas were read in from file "
2819 << collareaName << endl;
2820 gLog << "-----------------------------------------------" << endl;
2821 }
2822
2823
2824 // save binnings for call to CT1EEst
2825 MBinning *binsE = (MBinning*)parlist.FindObject("BinningE");
2826 if (!binsE)
2827 {
2828 gLog << "Object 'BinningE' not found in MParList" << endl;
2829 return;
2830 }
2831 MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
2832 if (!binsTheta)
2833 {
2834 gLog << "Object 'BinningTheta' not found in MParList" << endl;
2835 return;
2836 }
2837
2838 //-------------------------------------
2839 TString fHilName = "MHillas";
2840 TString fHilNameExt = "MHillasExt";
2841 TString fHilNameSrc = "MHillasSrc";
2842 TString fImgParName = "MNewImagePar";
2843
2844
2845 if (OEEst)
2846 {
2847 //===========================================================
2848 //
2849 // Optimization of energy estimator
2850 //
2851 gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
2852
2853 TString inpath("");
2854 TString outpath("");
2855 Int_t howMany = 2000;
2856 CT1EEst(inpath, filenameOpt, outpath, energyParName,
2857 fHilName, fHilNameSrc, fhadronnessName,
2858 howMany, maxhadronness, maxalpha, maxdist,
2859 binsE, binsTheta);
2860 gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
2861 }
2862
2863 if (WEX)
2864 {
2865 //-----------------------------------------------------------
2866 //
2867 // Read in parameters of energy estimator ("MMcEnergyEst")
2868 // and migration matrix ("MHMcEnergyMigration")
2869 //
2870 gLog << "================================================================"
2871 << endl;
2872 gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
2873 << energyParName << "'" << endl;
2874 TFile enparam(energyParName);
2875 enparam.ls();
2876 MMcEnergyEst mcest("MMcEnergyEst");
2877 mcest.Read("MMcEnergyEst");
2878
2879 //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
2880 gLog << "Parameters of energy estimator were read in" << endl;
2881
2882
2883 gLog << "Read in Migration matrix" << endl;
2884
2885 MHMcEnergyMigration mighiston("MHMcEnergyMigration");
2886 mighiston.Read("MHMcEnergyMigration");
2887 //MHMcEnergyMigration &mighiston =
2888 // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
2889
2890 gLog << "Migration matrix was read in" << endl;
2891
2892
2893 TArrayD parA(mcest.GetNumCoeffA());
2894 TArrayD parB(mcest.GetNumCoeffB());
2895 for (Int_t i=0; i<parA.GetSize(); i++)
2896 parA[i] = mcest.GetCoeff(i);
2897 for (Int_t i=0; i<parB.GetSize(); i++)
2898 parB[i] = mcest.GetCoeff( i+parA.GetSize() );
2899
2900 //*************************************************************************
2901 //
2902 // Analyse the data
2903 //
2904 gLog << "============================================================"
2905 << endl;
2906 gLog << "Analyse the data" << endl;
2907
2908 MTaskList tliston;
2909 MParList pliston;
2910
2911 // geometry is needed in MHHillas... classes
2912 MGeomCam *fGeom =
2913 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2914
2915
2916 //-------------------------------------------
2917 // create the tasks which should be executed
2918 //
2919
2920 MReadMarsFile read("Events", filenameData);
2921 read.DisableAutoScheme();
2922
2923 //.......................................................................
2924
2925 gLog << "CT1Analysis.C : write root file '" << filenameDataout
2926 << "'" << endl;
2927
2928 //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
2929
2930
2931 MWriteRootFile write(filenameDataout, "RECREATE");
2932
2933 write.AddContainer("MRawRunHeader", "RunHeaders");
2934 write.AddContainer("MTime", "Events");
2935 write.AddContainer("MMcEvt", "Events");
2936 write.AddContainer("ThetaOrig", "Events");
2937 write.AddContainer("MSrcPosCam", "Events");
2938 write.AddContainer("MSigmabar", "Events");
2939 write.AddContainer("MHillas", "Events");
2940 write.AddContainer("MHillasExt", "Events");
2941 write.AddContainer("MHillasSrc", "Events");
2942 write.AddContainer("MNewImagePar", "Events");
2943
2944 //write.AddContainer("HadNN", "Events");
2945 write.AddContainer("HadSC", "Events");
2946 write.AddContainer("HadRF", "Events");
2947
2948 write.AddContainer("MEnergyEst", "Events");
2949
2950
2951 //-----------------------------------------------------------------
2952 // geometry is needed in MHHillas... classes
2953 MGeomCam *fGeom =
2954 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2955
2956 MFCT1SelFinal selfinalgh(fHilNameSrc);
2957 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2958 selfinalgh.SetHadronnessName(fhadronnessName);
2959 selfinalgh.SetName("SelFinalgh");
2960 MContinue contfinalgh(&selfinalgh);
2961 contfinalgh.SetName("ContSelFinalgh");
2962
2963 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
2964 //fillhadnn.SetName("HhadNN");
2965 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
2966 fillhadsc.SetName("HhadSC");
2967 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
2968 fillhadrf.SetName("HhadRF");
2969
2970 //---------------------------
2971 // calculate estimated energy
2972
2973 MEnergyEstParam eeston(fHilName);
2974 eeston.Add(fHilNameSrc);
2975
2976 eeston.SetCoeffA(parA);
2977 eeston.SetCoeffB(parB);
2978
2979 //---------------------------
2980 // calculate estimated energy using Daniel's parameters
2981
2982 //MEnergyEstParamDanielMkn421 eeston(fHilName);
2983 //eeston.Add(fHilNameSrc);
2984 //eeston.SetCoeffA(parA);
2985 //eeston.SetCoeffB(parB);
2986
2987
2988 //---------------------------
2989
2990
2991 MFillH hfill1("MHHillas", fHilName);
2992 hfill1.SetName("HHillas");
2993
2994 MFillH hfill2("MHStarMap", fHilName);
2995 hfill2.SetName("HStarMap");
2996
2997 MFillH hfill3("MHHillasExt", fHilNameSrc);
2998 hfill3.SetName("HHillasExt");
2999
3000 MFillH hfill4("MHHillasSrc", fHilNameSrc);
3001 hfill4.SetName("HHillasSrc");
3002
3003 MFillH hfill5("MHNewImagePar", fImgParName);
3004 hfill5.SetName("HNewImagePar");
3005
3006 //---------------------------
3007 // new from Robert
3008
3009 MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
3010 hfill6.SetName("HTimeDiffTheta");
3011
3012 MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
3013 hfill6a.SetName("HTimeDiffTime");
3014
3015 MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
3016 hfill7.SetName("HAlphaEnergyTheta");
3017
3018 MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
3019 hfill7a.SetName("HAlphaEnergyTime");
3020
3021 MFillH hfill7b("MHThetabarTime", fHilNameSrc);
3022 hfill7b.SetName("HThetabarTime");
3023
3024 MFillH hfill7c("MHEnergyTime", "MMcEvt");
3025 hfill7c.SetName("HEnergyTime");
3026
3027
3028 //---------------------------
3029
3030 MFCT1SelFinal selfinal(fHilNameSrc);
3031 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
3032 selfinal.SetHadronnessName(fhadronnessName);
3033 selfinal.SetName("SelFinal");
3034 MContinue contfinal(&selfinal);
3035 contfinal.SetName("ContSelFinal");
3036
3037
3038 //*****************************
3039 // entries in MParList
3040
3041 pliston.AddToList(&tliston);
3042 InitBinnings(&pliston);
3043
3044
3045 //*****************************
3046 // entries in MTaskList
3047
3048 tliston.AddToList(&read);
3049
3050 // robert
3051 tliston.AddToList(&hfill6); //timediff
3052 tliston.AddToList(&hfill6a); //timediff
3053
3054 tliston.AddToList(&contfinalgh);
3055 tliston.AddToList(&eeston);
3056
3057 tliston.AddToList(&write);
3058
3059 //tliston.AddToList(&fillhadnn);
3060 tliston.AddToList(&fillhadsc);
3061 tliston.AddToList(&fillhadrf);
3062
3063 tliston.AddToList(&hfill1);
3064 tliston.AddToList(&hfill2);
3065 tliston.AddToList(&hfill3);
3066 tliston.AddToList(&hfill4);
3067 tliston.AddToList(&hfill5);
3068
3069 //robert
3070 tliston.AddToList(&hfill7);
3071 tliston.AddToList(&hfill7a);
3072 tliston.AddToList(&hfill7b);
3073 tliston.AddToList(&hfill7c);
3074
3075 tliston.AddToList(&contfinal);
3076
3077 //*****************************
3078
3079 //-------------------------------------------
3080 // Execute event loop
3081 //
3082 MProgressBar bar;
3083 MEvtLoop evtloop;
3084 evtloop.SetParList(&pliston);
3085 evtloop.SetProgressBar(&bar);
3086
3087 Int_t maxevents = -1;
3088 if ( !evtloop.Eventloop(maxevents) )
3089 return;
3090
3091 tliston.PrintStatistics(0, kTRUE);
3092
3093
3094 //-------------------------------------------
3095 // Display the histograms
3096 //
3097
3098 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
3099
3100 gLog << "before hadRF" << endl;
3101 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
3102
3103 gLog << "before hadSC" << endl;
3104 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
3105
3106 gLog << "before MHHillas" << endl;
3107 pliston.FindObject("MHHillas")->DrawClone();
3108
3109 gLog << "before MHHillasExt" << endl;
3110 pliston.FindObject("MHHillasExt")->DrawClone();
3111
3112 gLog << "before MHHillasSrc" << endl;
3113 pliston.FindObject("MHHillasSrc")->DrawClone();
3114
3115 gLog << "before MHNewImagePar" << endl;
3116 pliston.FindObject("MHNewImagePar")->DrawClone();
3117
3118 gLog << "before MHStarMap" << endl;
3119 pliston.FindObject("MHStarMap")->DrawClone();
3120
3121 gLog << "before DeleteBinnings" << endl;
3122
3123 DeleteBinnings(&pliston);
3124
3125 gLog << "before Robert's code" << endl;
3126
3127
3128//rwagner write all relevant histograms onto a file
3129
3130 if (WRobert)
3131 {
3132 gLog << "=======================================================" << endl;
3133 gLog << "Write results onto file '" << filenameResults << "'" << endl;
3134
3135 TFile outfile(filenameResults,"recreate");
3136
3137 MHHillasSrc* hillasSrc =
3138 (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
3139 TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
3140 alphaHist->Write();
3141 gLog << "Alpha plot has been written out" << endl;
3142
3143
3144 MHAlphaEnergyTheta* aetH =
3145 (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
3146 TH3D* aetHist = (TH3D*)(aetH->GetHist());
3147 aetHist->SetName("aetHist");
3148 aetHist->Write();
3149 gLog << "AlphaEnergyTheta plot has been written out" << endl;
3150
3151 MHAlphaEnergyTime* aetH2 =
3152 (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
3153 TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
3154 aetHist2->SetName("aetimeHist");
3155// aetHist2->DrawClone();
3156 aetHist2->Write();
3157 gLog << "AlphaEnergyTime plot has been written out" << endl;
3158
3159 MHThetabarTime* tbt =
3160 (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
3161 TProfile* tbtHist = (TProfile*)(tbt->GetHist());
3162 tbtHist->SetName("tbtHist");
3163 tbtHist->Write();
3164 gLog << "ThetabarTime plot has been written out" << endl;
3165
3166 MHEnergyTime* ent =
3167 (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
3168 TH2D* entHist = (TH2D*)(ent->GetHist());
3169 entHist->SetName("entHist");
3170 entHist->Write();
3171 gLog << "EnergyTime plot has been written out" << endl;
3172
3173 MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
3174 TH2D* timeHist = (TH2D*)(time->GetHist());
3175 timeHist->SetName("MHTimeDiffTheta");
3176 timeHist->SetTitle("Time diffs");
3177 timeHist->Write();
3178 gLog << "TimeDiffTheta plot has been written out" << endl;
3179
3180
3181 MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
3182 TH2D* timeHist2 = (TH2D*)(time2->GetHist());
3183 timeHist2->SetName("MHTimeDiffTime");
3184 timeHist2->SetTitle("Time diffs");
3185 timeHist2->Write();
3186 gLog << "TimeDiffTime plot has been written out" << endl;
3187
3188//rwagner write also collareas to same file
3189 collarea->GetHist()->Write();
3190 collarea->GetHAll()->Write();
3191 collarea->GetHSel()->Write();
3192 gLog << "Effective collection areas have been written out" << endl;
3193
3194//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
3195//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
3196
3197 gLog << "before closing outfile" << endl;
3198
3199 //outfile.Close();
3200 gLog << "Results were written onto file '" << filenameResults
3201 << "'" << endl;
3202 gLog << "=======================================================" << endl;
3203 }
3204
3205 }
3206
3207 gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
3208 gLog << "=======================================================" << endl;
3209 }
3210 //---------------------------------------------------------------------
3211
3212}
3213
3214
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227
3228
3229
3230
3231
3232
3233
3234
3235
3236
3237
3238
3239
3240
3241
Note: See TracBrowser for help on using the repository browser.