source: trunk/MagicSoft/Mars/macros/CT1Analysis.C@ 3990

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