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

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