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

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