source: trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C@ 4571

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