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

Last change on this file since 2550 was 2543, checked in by wittek, 22 years ago
*** empty log message ***
File size: 90.6 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 // Job B_RF_UP
868 //============
869
870
871 // - create (or read in) the matrices of training events for gammas
872 // and hadrons
873 // - create (or read in) the trees
874 // - then read ON1.root (or MC1.root) file
875 // - calculate the hadroness for the method of RANDOM FOREST
876 // - update input root file with the hadroness
877
878
879 if (JobB_RF_UP)
880 {
881 gLog << "=====================================================" << endl;
882 gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
883
884 gLog << "" << endl;
885 gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
886 << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", "
887 << (RTrainRF ? "kTRUE" : "kFALSE") << ", "
888 << (CTrainRF ? "kTRUE" : "kFALSE") << ", "
889 << (RTree ? "kTRUE" : "kFALSE") << ", "
890 << (WRF ? "kTRUE" : "kFALSE") << endl;
891
892
893 //--------------------------------------------
894 // parameters for the random forest
895 Int_t NumTrees = 100;
896 Int_t NumTry = 3;
897 Int_t NdSize = 1;
898
899 // cut in RF hadronness
900 TString hadRFName = "HadRF";
901 Float_t maxhadronness = 0.23;
902 Float_t maxalpha = 20.0;
903 Float_t maxdist = 10.0;
904
905 TString fHilName = "MHillas";
906 TString fHilNameExt = "MHillasExt";
907 TString fHilNameSrc = "MHillasSrc";
908 TString fImgParName = "MNewImagePar";
909
910
911 //--------------------------------------------
912 // file to be updated (either ON or MC)
913
914 //TString typeInput = "ON";
915 TString typeInput = "MC";
916 gLog << "typeInput = " << typeInput << endl;
917
918 // name of input root file
919 TString filenameData = outPath;
920 filenameData += typeInput;
921 filenameData += "1.root";
922 gLog << "filenameData = " << filenameData << endl;
923
924 // name of output root file
925 TString outNameImage = outPath;
926 outNameImage += typeInput;
927 outNameImage += "2.root";
928 //TString outNameImage = filenameData;
929
930 gLog << "outNameImage = " << outNameImage << endl;
931
932 //--------------------------------------------
933 // files to be read for generating the matrices of training events
934 // "hadrons" :
935 TString filenameON = outPath;
936 filenameON += "ON";
937 filenameON += "1.root";
938 Int_t howManyHadrons = 35000;
939 gLog << "filenameON = " << filenameON << ", howManyHadrons = "
940 << howManyHadrons << endl;
941
942
943 // "gammas" :
944 TString filenameMC = outPath;
945 filenameMC += "MC";
946 filenameMC += "1.root";
947 Int_t howManyGammas = 35000;
948 gLog << "filenameMC = " << filenameMC << ", howManyGammas = "
949 << howManyGammas << endl;
950
951 //--------------------------------------------
952 // files of training events
953
954 TString outNameGammas = outPath;
955 outNameGammas += "RFmatrix_gammas_";
956 outNameGammas += "MC";
957 outNameGammas += ".root";
958
959 TString typeMatrixHadrons = "ON";
960 gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
961
962 TString outNameHadrons = outPath;
963 outNameHadrons += "RFmatrix_hadrons_";
964 outNameHadrons += typeMatrixHadrons;
965 outNameHadrons += ".root";
966
967
968 MHMatrix matrixg("MatrixGammas");
969 matrixg.EnableGraphicalOutput();
970
971 matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
972 matrixg.AddColumn("MSigmabar.fSigmabar");
973 matrixg.AddColumn("log10(MHillas.fSize)");
974 matrixg.AddColumn("MHillasSrc.fDist");
975 matrixg.AddColumn("MHillas.fWidth");
976 matrixg.AddColumn("MHillas.fLength");
977 matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
978 matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
979 matrixg.AddColumn("MNewImagePar.fConc");
980 matrixg.AddColumn("MNewImagePar.fLeakage1");
981
982 MHMatrix matrixh("MatrixHadrons");
983 matrixh.EnableGraphicalOutput();
984
985 matrixh.AddColumns(matrixg.GetColumns());
986
987 //--------------------------------------------
988 // file of trees of the random forest
989
990 TString outRF = outPath;
991 outRF += "RF.root";
992
993
994 //*************************************************************************
995 // read in matrices of training events
996if (RTrainRF)
997 {
998 const char* mtxName = "MatrixGammas";
999
1000 gLog << "" << endl;
1001 gLog << "========================================================" << endl;
1002 gLog << "Get matrix for (gammas)" << endl;
1003 gLog << "matrix name = " << mtxName << endl;
1004 gLog << "name of root file = " << outNameGammas << endl;
1005 gLog << "" << endl;
1006
1007
1008 // read in the object with the name 'mtxName' from file 'outNameGammas'
1009 //
1010 TFile fileg(outNameGammas);
1011
1012 matrixg.Read(mtxName);
1013 matrixg.Print("SizeCols");
1014
1015
1016 //*****************************************************************
1017
1018 const char* mtxName = "MatrixHadrons";
1019
1020 gLog << "" << endl;
1021 gLog << "========================================================" << endl;
1022 gLog << " Get matrix for (hadrons)" << endl;
1023 gLog << "matrix name = " << mtxName << endl;
1024 gLog << "name of root file = " << outNameHadrons << endl;
1025 gLog << "" << endl;
1026
1027
1028 // read in the object with the name 'mtxName' from file 'outNameHadrons'
1029 //
1030 TFile fileh(outNameHadrons);
1031
1032 matrixh.Read(mtxName);
1033 matrixh.Print("SizeCols");
1034 }
1035
1036
1037 //*************************************************************************
1038 // create matrices of training events
1039if (CTrainRF)
1040 {
1041 gLog << "" << endl;
1042 gLog << "========================================================" << endl;
1043 gLog << " Create matrices of training events" << endl;
1044 gLog << " Gammas :" << endl;
1045
1046
1047 MParList plistg;
1048 MTaskList tlistg;
1049 MFilterList flistg;
1050
1051 MParList plisth;
1052 MTaskList tlisth;
1053 MFilterList flisth;
1054
1055 MReadMarsFile readg("Events", filenameMC);
1056 readg.DisableAutoScheme();
1057
1058 MReadMarsFile readh("Events", filenameON);
1059 readh.DisableAutoScheme();
1060
1061 MFParticleId fgamma("MMcEvt", '=', kGAMMA);
1062 fgamma.SetName("gammaID");
1063
1064 MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
1065 fhadrons.SetName("hadronID)");
1066
1067 TString mgname("costhg");
1068 MBinning bing("Binning"+mgname);
1069 bing.SetEdges(10, 0., 1.0);
1070
1071 MH3 gref("cos(MMcEvt.fTelescopeTheta)");
1072 gref.SetName(mgname);
1073 MH::SetBinning(&gref.GetHist(), &bing);
1074 for (Int_t i=1; i<=gref.GetNbins(); i++)
1075 gref.GetHist().SetBinContent(i, 1.0);
1076
1077
1078 MFEventSelector2 selectorg(gref);
1079 selectorg.SetNumMax(howManyGammas);
1080 selectorg.SetName("selectGammas");
1081
1082
1083 MFillH fillmatg("MatrixGammas");
1084 fillmatg.SetFilter(&flistg);
1085 fillmatg.SetName("fillGammas");
1086
1087
1088
1089 //***************************** fill gammas ***
1090 // entries in MFilterList
1091
1092 flistg.AddToList(&fgamma);
1093 flistg.AddToList(&selectorg);
1094
1095 //*****************************
1096 // entries in MParList
1097
1098 plistg.AddToList(&tlistg);
1099 InitBinnings(&plistg);
1100
1101 plistg.AddToList(&matrixg);
1102
1103 //*****************************
1104 // entries in MTaskList
1105
1106 tlistg.AddToList(&readg);
1107 tlistg.AddToList(&flistg);
1108 tlistg.AddToList(&fillmatg);
1109
1110 //*****************************
1111
1112 MProgressBar matrixbar;
1113 MEvtLoop evtloopg;
1114 evtloopg.SetParList(&plistg);
1115 //evtloopg.ReadEnv(env, "", printEnv);
1116 evtloopg.SetProgressBar(&matrixbar);
1117
1118 Int_t maxevents = -1;
1119 if (!evtloopg.Eventloop(maxevents))
1120 return;
1121
1122 tlistg.PrintStatistics(0, kTRUE);
1123
1124
1125 //***************************** fill hadrons ***
1126
1127 gLog << " Hadrons :" << endl;
1128
1129 TString mhname("costhh");
1130 MBinning binh("Binning"+mhname);
1131 binh.SetEdges(10, 0., 1.0);
1132
1133 MH3 href("cos(MMcEvt.fTelescopeTheta)");
1134 href.SetName(mhname);
1135 MH::SetBinning(&href.GetHist(), &binh);
1136 for (Int_t i=1; i<=href.GetNbins(); i++)
1137 href.GetHist().SetBinContent(i, 1.0);
1138
1139 MFEventSelector2 selectorh(href);
1140 //selectorh.SetNumMax(howManyHadrons);
1141 // select as many hadrons as gammas
1142 selectorh.SetNumMax(matrixg.GetM().GetNrows());
1143 selectorh.SetName("selectHadrons");
1144
1145 MFillH fillmath("MatrixHadrons");
1146 fillmath.SetFilter(&flisth);
1147 fillmath.SetName("fillHadrons");
1148
1149
1150 // entries in MFilterList
1151
1152 flisth.AddToList(&fhadrons);
1153 flisth.AddToList(&selectorh);
1154
1155 //*****************************
1156 // entries in MParList
1157
1158 plisth.AddToList(&tlisth);
1159 InitBinnings(&plisth);
1160
1161 plisth.AddToList(&matrixh);
1162
1163 //*****************************
1164 // entries in MTaskList
1165
1166 tlisth.AddToList(&readh);
1167 tlisth.AddToList(&flisth);
1168 tlisth.AddToList(&fillmath);
1169
1170 //*****************************
1171
1172 MProgressBar matrixbar;
1173 MEvtLoop evtlooph;
1174 evtlooph.SetParList(&plisth);
1175 //evtlooph.ReadEnv(env, "", printEnv);
1176 evtlooph.SetProgressBar(&matrixbar);
1177
1178 Int_t maxevents = -1;
1179 if (!evtlooph.Eventloop(maxevents))
1180 return;
1181
1182 tlisth.PrintStatistics(0, kTRUE);
1183
1184
1185
1186 //*****************************************************
1187
1188
1189 // write out matrices of training events
1190
1191 gLog << "" << endl;
1192 gLog << "========================================================" << endl;
1193 gLog << "Write out matrices of training events" << endl;
1194
1195
1196 //-------------------------------------------
1197 // "gammas"
1198 gLog << "Gammas :" << endl;
1199 matrixg.Print("SizeCols");
1200
1201 TFile writeg(outNameGammas, "RECREATE", "");
1202 matrixg.Write();
1203
1204 gLog << "" << endl;
1205 gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file "
1206 << outNameGammas << endl;
1207
1208 //-------------------------------------------
1209 // "hadrons"
1210 gLog << "Hadrons :" << endl;
1211 matrixh.Print("SizeCols");
1212
1213 TFile writeh(outNameHadrons, "RECREATE", "");
1214 matrixh.Write();
1215
1216 gLog << "" << endl;
1217 gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file "
1218 << outNameHadrons << endl;
1219
1220 }
1221 //********** end of creating matrices of training events ***********
1222
1223
1224 MRanForest *fRanForest;
1225 MRanTree *fRanTree;
1226 //-----------------------------------------------------------------
1227 // read in the trees of the random forest
1228 if (RTree)
1229 {
1230 MParList plisttr;
1231 MTaskList tlisttr;
1232 plisttr.AddToList(&tlisttr);
1233
1234 MReadTree readtr("TREE", outRF);
1235 readtr.DisableAutoScheme();
1236
1237 MRanForestFill rffill;
1238 rffill.SetNumTrees(NumTrees);
1239
1240 // list of tasks for the loop over the trees
1241
1242 tlisttr.AddToList(&readtr);
1243 tlisttr.AddToList(&rffill);
1244
1245 //-------------------
1246 // Execute tree loop
1247 //
1248 MEvtLoop evtlooptr;
1249 evtlooptr.SetName("ReadRFTrees");
1250 evtlooptr.SetParList(&plisttr);
1251 if (!evtlooptr.Eventloop())
1252 return;
1253
1254 tlisttr.PrintStatistics(0, kTRUE);
1255
1256
1257 // get adresses of objects which are used in the next eventloop
1258 fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
1259 if (!fRanForest)
1260 {
1261 *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
1262 return kFALSE;
1263 }
1264
1265 fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
1266 if (!fRanTree)
1267 {
1268 *fLog << err << dbginf << "MRanTree not found... aborting." << endl;
1269 return kFALSE;
1270 }
1271
1272 }
1273
1274 //-----------------------------------------------------------------
1275 // grow the trees of the random forest (event loop = tree loop)
1276
1277 if (!RTree)
1278 {
1279
1280 gLog << "" << endl;
1281 gLog << "========================================================" << endl;
1282 gLog << "Macro CT1Analysis : start growing trees" << endl;
1283
1284 MTaskList tlist2;
1285 MParList plist2;
1286 plist2.AddToList(&tlist2);
1287
1288 plist2.AddToList(&matrixg);
1289 plist2.AddToList(&matrixh);
1290
1291 MRanForestGrow rfgrow2;
1292 rfgrow2.SetNumTrees(NumTrees);
1293 rfgrow2.SetNumTry(NumTry);
1294 rfgrow2.SetNdSize(NdSize);
1295
1296 MWriteRootFile rfwrite2(outRF);
1297 rfwrite2.AddContainer("MRanTree", "TREE");
1298
1299 MFillH fillh2("MHRanForestGini");
1300
1301 // list of tasks for the loop over the trees
1302
1303 tlist2.AddToList(&rfgrow2);
1304 tlist2.AddToList(&rfwrite2);
1305 tlist2.AddToList(&fillh2);
1306
1307 //-------------------
1308 // Execute tree loop
1309 //
1310 MEvtLoop treeloop;
1311 treeloop.SetName("GrowRFTrees");
1312 treeloop.SetParList(&plist2);
1313
1314 if ( !treeloop.Eventloop() )
1315 return;
1316
1317 tlist2.PrintStatistics(0, kTRUE);
1318
1319 plist2.FindObject("MHRanForestGini")->DrawClone();
1320
1321
1322 // get adresses of objects which are used in the next eventloop
1323 fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
1324 if (!fRanForest)
1325 {
1326 *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
1327 return kFALSE;
1328 }
1329
1330 fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
1331 if (!fRanTree)
1332 {
1333 *fLog << err << dbginf << "MRanTree not found... aborting." << endl;
1334 return kFALSE;
1335 }
1336
1337 }
1338 // end of growing the trees of the random forest
1339 //-----------------------------------------------------------------
1340
1341
1342
1343 //-----------------------------------------------------------------
1344 // Update the input files with the RF hadronness
1345 //
1346 if (WRF)
1347 {
1348 gLog << "" << endl;
1349 gLog << "========================================================" << endl;
1350 gLog << "Update input file '" << filenameData
1351 << "' with the RF hadronness" << endl;
1352
1353 MTaskList tliston;
1354 MParList pliston;
1355
1356
1357 // geometry is needed in MHHillas... classes
1358 MGeomCam *fGeom =
1359 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1360
1361 //-------------------------------------------
1362 // create the tasks which should be executed
1363 //
1364
1365 MReadMarsFile read("Events", filenameData);
1366 read.DisableAutoScheme();
1367
1368
1369 //.......................................................................
1370 // calculate hadronnes for method of RANDOM FOREST
1371
1372
1373 MRanForestCalc rfcalc;
1374 rfcalc.SetHadronnessName(hadRFName);
1375
1376 //.......................................................................
1377
1378 //MWriteRootFile write(outNameImage, "UPDATE");
1379 MWriteRootFile write(outNameImage, "RECREATE");
1380
1381 write.AddContainer("MRawRunHeader", "RunHeaders");
1382 write.AddContainer("MTime", "Events");
1383 write.AddContainer("MMcEvt", "Events");
1384 write.AddContainer("ThetaOrig", "Events");
1385 write.AddContainer("MSrcPosCam", "Events");
1386 write.AddContainer("MSigmabar", "Events");
1387 write.AddContainer("MHillas", "Events");
1388 write.AddContainer("MHillasExt", "Events");
1389 write.AddContainer("MHillasSrc", "Events");
1390 write.AddContainer("MNewImagePar", "Events");
1391
1392 write.AddContainer(hadRFName, "Events");
1393
1394
1395 //-----------------------------------------------------------------
1396
1397
1398 MFCT1SelFinal selfinalgh(fHilNameSrc);
1399 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1400 selfinalgh.SetHadronnessName(hadRFName);
1401 selfinalgh.SetName("SelFinalgh");
1402 MContinue contfinalgh(&selfinalgh);
1403 contfinalgh.SetName("ContSelFinalgh");
1404
1405 MFillH fillranfor("MHRanForest");
1406 fillranfor.SetName("HRanForest");
1407
1408 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
1409 fillhadrf.SetName("HhadRF");
1410
1411 MFCT1SelFinal selfinal(fHilNameSrc);
1412 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1413 selfinal.SetHadronnessName(hadRFName);
1414 selfinal.SetName("SelFinal");
1415 MContinue contfinal(&selfinal);
1416 contfinal.SetName("ContSelFinal");
1417
1418 TString mh3name = "abs(Alpha)";
1419 MBinning binsalphaabs("Binning"+mh3name);
1420 binsalphaabs.SetEdges(50, -2.0, 98.0);
1421
1422 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
1423 alphaabs.SetName(mh3name);
1424 MFillH alpha(&alphaabs);
1425 alpha.SetName("FillAlphaAbs");
1426
1427 MFillH hfill1("MHHillas", fHilName);
1428 hfill1.SetName("HHillas");
1429
1430 MFillH hfill2("MHStarMap", fHilName);
1431 hfill2.SetName("HStarMap");
1432
1433 MFillH hfill3("MHHillasExt", fHilNameSrc);
1434 hfill3.SetName("HHillasExt");
1435
1436 MFillH hfill4("MHHillasSrc", fHilNameSrc);
1437 hfill4.SetName("HHillasSrc");
1438
1439 MFillH hfill5("MHNewImagePar", fImgParName);
1440 hfill5.SetName("HNewImagePar");
1441
1442 //*****************************
1443 // entries in MParList
1444
1445 pliston.AddToList(&tliston);
1446 InitBinnings(&pliston);
1447
1448 pliston.AddToList(fRanForest);
1449 pliston.AddToList(fRanTree);
1450
1451 pliston.AddToList(&binsalphaabs);
1452 pliston.AddToList(&alphaabs);
1453
1454 //*****************************
1455 // entries in MTaskList
1456
1457 tliston.AddToList(&read);
1458
1459 tliston.AddToList(&rfcalc);
1460 tliston.AddToList(&fillranfor);
1461 tliston.AddToList(&fillhadrf);
1462
1463 tliston.AddToList(&write);
1464 tliston.AddToList(&contfinalgh);
1465
1466 tliston.AddToList(&alpha);
1467 tliston.AddToList(&hfill1);
1468 tliston.AddToList(&hfill2);
1469 tliston.AddToList(&hfill3);
1470 tliston.AddToList(&hfill4);
1471 tliston.AddToList(&hfill5);
1472
1473 tliston.AddToList(&contfinal);
1474
1475 //*****************************
1476
1477 //-------------------------------------------
1478 // Execute event loop
1479 //
1480 MProgressBar bar;
1481 MEvtLoop evtloop;
1482 evtloop.SetParList(&pliston);
1483 evtloop.SetProgressBar(&bar);
1484
1485 Int_t maxevents = -1;
1486 if ( !evtloop.Eventloop(maxevents) )
1487 return;
1488
1489 tliston.PrintStatistics(0, kTRUE);
1490
1491
1492 //-------------------------------------------
1493 // Display the histograms
1494 //
1495 pliston.FindObject("MHRanForest")->DrawClone();
1496 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
1497
1498 pliston.FindObject("MHHillas")->DrawClone();
1499 pliston.FindObject("MHHillasExt")->DrawClone();
1500 pliston.FindObject("MHHillasSrc")->DrawClone();
1501 pliston.FindObject("MHNewImagePar")->DrawClone();
1502 pliston.FindObject("MHStarMap")->DrawClone();
1503
1504 //-------------------------------------------
1505 // fit alpha distribution to get the number of excess events and
1506 // calculate significance of gamma signal in the alpha plot
1507
1508 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
1509 TH1 &alphaHist = absalpha->GetHist();
1510 alphaHist.SetXTitle("|alpha| [\\circ]");
1511 alphaHist.SetName("alpha-macro");
1512
1513 Double_t alphasig = 13.1;
1514 Double_t alphamin = 30.0;
1515 Double_t alphamax = 90.0;
1516 Int_t degree = 2;
1517 Double_t significance = -99.0;
1518 Bool_t drawpoly = kTRUE;
1519 Bool_t fitgauss = kTRUE;
1520 Bool_t print = kTRUE;
1521
1522 MHFindSignificance findsig;
1523 findsig.SetRebin(kTRUE);
1524 findsig.SetReduceDegree(kFALSE);
1525
1526 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
1527 alphasig, drawpoly, fitgauss, print);
1528 significance = findsig.GetSignificance();
1529 Float_t alphasi = findsig.GetAlphasi();
1530
1531 gLog << "For file '" << filenameData << "' : " << endl;
1532 gLog << "Significance of gamma signal after supercuts : "
1533 << significance << " (for |alpha| < " << alphasi << " degrees)"
1534 << endl;
1535
1536 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
1537
1538 //-------------------------------------------
1539
1540
1541 DeleteBinnings(&pliston);
1542 }
1543
1544 gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
1545 gLog << "=======================================================" << endl;
1546 }
1547 //---------------------------------------------------------------------
1548
1549
1550
1551
1552 //---------------------------------------------------------------------
1553 // Job B_SC_UP
1554 //============
1555
1556 // - create (or read in) optimum supercuts parameter values
1557 //
1558 // - calculate the hadroness for the supercuts
1559 //
1560 // - update input root file, including the hadroness
1561
1562
1563 if (JobB_SC_UP)
1564 {
1565 gLog << "=====================================================" << endl;
1566 gLog << "Macro CT1Analysis : Start of Job B_SC_UP" << endl;
1567
1568 gLog << "" << endl;
1569 gLog << "Macro CT1Analysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
1570 << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", "
1571 << (CMatrix ? "kTRUE" : "kFALSE") << ", "
1572 << (RMatrix ? "kTRUE" : "kFALSE") << ", "
1573 << (WOptimize ? "kTRUE" : "kFALSE") << ", "
1574 << (RTest ? "kTRUE" : "kFALSE") << ", "
1575 << (WSC ? "kTRUE" : "kFALSE") << endl;
1576
1577
1578 //--------------------------------------------
1579 // file which contains the initial parameter values for the supercuts
1580 // if parSCinit ="" the initial values are taken from the constructor of
1581 // MCT1Supercuts
1582
1583 TString parSCinit = outPath;
1584 //parSCinit += "parSC_1709d";
1585 parSCinit = "";
1586
1587 gLog << "parSCinit = " << parSCinit << endl;
1588
1589 //---------------
1590 // file onto which the optimal parameter values for the supercuts
1591 // are written
1592
1593 TString parSCfile = outPath;
1594 parSCfile += "parSC_2310a";
1595
1596 gLog << "parSCfile = " << parSCfile << endl;
1597
1598 //--------------------------------------------
1599 // file to be updated (either ON or MC)
1600
1601 //TString typeInput = "ON";
1602 TString typeInput = "MC";
1603 gLog << "typeInput = " << typeInput << endl;
1604
1605 // name of input root file
1606 TString filenameData = outPath;
1607 filenameData += typeInput;
1608 filenameData += "2.root";
1609 gLog << "filenameData = " << filenameData << endl;
1610
1611 // name of output root file
1612 TString outNameImage = outPath;
1613 outNameImage += typeInput;
1614 outNameImage += "3.root";
1615
1616
1617 //TString outNameImage = filenameData;
1618
1619 gLog << "outNameImage = " << outNameImage << endl;
1620
1621 //--------------------------------------------
1622 // files to be read for optimizing the supercuts
1623 //
1624 // for the training
1625 TString filenameTrain = outPath;
1626 filenameTrain += "ON";
1627 filenameTrain += "1.root";
1628 Int_t howManyTrain = 800000;
1629 gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = "
1630 << howManyTrain << endl;
1631
1632 // for testing
1633 TString filenameTest = outPath;
1634 filenameTest += "ON";
1635 filenameTest += "1.root";
1636 Int_t howManyTest = 800000;
1637
1638 gLog << "filenameTest = " << filenameTest << ", howManyTest = "
1639 << howManyTest << endl;
1640
1641
1642 //--------------------------------------------
1643 // files to contain the matrices (generated from filenameTrain and
1644 // filenameTest)
1645 //
1646 // for the training
1647 TString fileMatrixTrain = outPath;
1648 fileMatrixTrain += "MatrixTrainSC";
1649 fileMatrixTrain += ".root";
1650 gLog << "fileMatrixTrain = " << fileMatrixTrain << endl;
1651
1652 // for testing
1653 TString fileMatrixTest = outPath;
1654 fileMatrixTest += "MatrixTestSC";
1655 fileMatrixTest += ".root";
1656 gLog << "fileMatrixTest = " << fileMatrixTest << endl;
1657
1658
1659
1660 //---------------------------------------------------------------------
1661 // Training and test matrices :
1662 // - either create them and write them onto a file
1663 // - or read them from a file
1664
1665
1666 MCT1FindSupercuts findsuper;
1667 findsuper.SetFilenameParam(parSCfile);
1668 findsuper.SetHadronnessName("HadSC");
1669
1670
1671 //--------------------------
1672 // create matrices and write them onto files
1673 if (CMatrix)
1674 {
1675 MH3 &mh3 = *(new MH3("MHillas.fSize"));
1676 mh3.SetName("Target distribution for SIZE");
1677
1678 if (filenameTrain == filenameTest)
1679 {
1680 if ( !findsuper.DefineTrainTestMatrix(filenameTrain,
1681 howManyTrain, mh3, howManyTest, mh3,
1682 fileMatrixTrain, fileMatrixTest) )
1683 {
1684 *fLog << "CT1Analysis.C : DefineTrainTestMatrix failed" << endl;
1685 return;
1686 }
1687
1688 }
1689 else
1690 {
1691 if ( !findsuper.DefineTrainMatrix(filenameTrain,
1692 howManyTrain, mh3, fileMatrixTrain) )
1693 {
1694 *fLog << "CT1Analysis.C : DefineTrainMatrix failed" << endl;
1695 return;
1696 }
1697
1698 if ( !findsuper.DefineTestMatrix( filenameTest,
1699 howManyTest, mh3, fileMatrixTest) )
1700 {
1701 *fLog << "CT1Analysis.C : DefineTestMatrix failed" << endl;
1702 return;
1703 }
1704 }
1705 }
1706
1707 //--------------------------
1708 // read matrices from files
1709 //
1710
1711 if (RMatrix)
1712 findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
1713 //--------------------------
1714
1715
1716
1717 //---------------------------------------------------------------------
1718 // optimize supercuts using the training sample
1719 //
1720 // the initial values are taken
1721 // - from the file parSCinit (if != "")
1722 // - or from the arrays params and steps (if their sizes are != 0)
1723 // - or from the MCT1Supercuts constructor
1724
1725
1726if (WOptimize)
1727 {
1728 gLog << "CT1Analysis.C : optimize the supercuts using the training matrix"
1729 << endl;
1730
1731 TArrayD params(0);
1732 TArrayD steps(0);
1733
1734 if (parSCinit == "")
1735 {
1736 Double_t vparams[104] = {
1737 // LengthUp
1738 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
1739 0.007388, -0.013463,
1740 // LengthLo
1741 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993,
1742 0.000080, -0.007157,
1743 // WidthUp
1744 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353,
1745 0.020711, -0.016703,
1746 // WidthLo
1747 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
1748 0.006126, -0.002849,
1749 // DistUp
1750 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909,
1751 0.189697, 0.0,
1752 // DistLo
1753 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374,
1754 -0.001750, 0.0,
1755 // AsymUp
1756 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1757 0.0, 0.0,
1758 // AsymLo
1759 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1760 0.0, 0.0,
1761 // ConcUp
1762 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1763 0.0, 0.0,
1764 // ConcLo
1765 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1766 0.0, 0.0,
1767 // Leakage1Up
1768 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1769 0.0, 0.0,
1770 // Leakage1Lo
1771 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1772 0.0, 0.0,
1773 // AlphaUp
1774 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0,
1775 0.0, 0.0 };
1776
1777 Double_t vsteps[104] = {
1778 // LengthUp
1779 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002,
1780 0.0008, 0.002,
1781 // LengthLo
1782 0.02, 0.003, 0.05, 0.006, 0.002, 0.3,
1783 0.0001, 0.0008,
1784 // WidthUp
1785 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006,
1786 0.002, 0.002,
1787 // WidthLo
1788 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002,
1789 0.0007, 0.003,
1790 // DistUp
1791 0.2, 0.0, 0.3, 0.02, 0.0, 0.03,
1792 0.02, 0.0
1793 // DistLo
1794 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005,
1795 0.0002, 0.0
1796 // AsymUp
1797 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1798 0.0, 0.0,
1799 // AsymLo
1800 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1801 0.0, 0.0,
1802 // ConcUp
1803 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1804 0.0, 0.0,
1805 // ConcLo
1806 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1807 0.0, 0.0,
1808 // Leakage1Up
1809 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1810 0.0, 0.0,
1811 // Leakage1Lo
1812 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1813 0.0, 0.0,
1814 // AlphaUp
1815 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1816 0.0, 0.0 };
1817
1818 params.Set(104, vparams);
1819 steps.Set (104, vsteps );
1820 }
1821
1822 Bool_t rf;
1823 rf = findsuper.FindParams(parSCinit, params, steps);
1824
1825 if (!rf)
1826 {
1827 gLog << "CT1Analysis.C : optimization of supercuts failed" << endl;
1828 return;
1829 }
1830 }
1831
1832 //--------------------------------------
1833 // test the supercuts on the test sample
1834 //
1835
1836 if (RTest)
1837 {
1838 gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl;
1839 Bool_t rt = findsuper.TestParams();
1840 if (!rt)
1841 {
1842 gLog << "CT1Analysis.C : test of supercuts on the test matrix failed"
1843 << endl;
1844 }
1845
1846 }
1847
1848
1849 //-----------------------------------------------------------------
1850 // Update the input files with the SC hadronness
1851 //
1852
1853 if (WSC)
1854 {
1855 gLog << "" << endl;
1856 gLog << "========================================================" << endl;
1857 gLog << "Update input file '" << filenameData
1858 << "' with the SC hadronness" << endl;
1859
1860
1861 //----------------------------------------------------
1862 // read in optimum parameter values for the supercuts
1863
1864 TFile inparam(parSCfile);
1865 MCT1Supercuts scin;
1866 scin.Read("MCT1Supercuts");
1867 inparam.Close();
1868
1869 gLog << "Parameter values for supercuts were read in from file '"
1870 << parSCfile << "'" << endl;
1871
1872 TArrayD supercutsPar;
1873 supercutsPar = scin.GetParameters();
1874
1875 TArrayD supercutsStep;
1876 supercutsStep = scin.GetStepsizes();
1877
1878 gLog << "Parameter values for supercuts : " << endl;
1879 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
1880 {
1881 gLog << supercutsPar[i] << ", ";
1882 }
1883 gLog << endl;
1884
1885 gLog << "Step values for supercuts : " << endl;
1886 for (Int_t i=0; i<supercutsStep.GetSize(); i++)
1887 {
1888 gLog << supercutsStep[i] << ", ";
1889 }
1890 gLog << endl;
1891
1892
1893 //----------------------------------------------------
1894 MTaskList tliston;
1895 MParList pliston;
1896
1897 // set the parameters of the supercuts
1898 MCT1Supercuts supercuts;
1899 supercuts.SetParameters(supercutsPar);
1900 gLog << "parameter values for the supercuts used for updating the input file ' "
1901 << filenameData << "'" << endl;
1902 supercutsPar = supercuts.GetParameters();
1903 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
1904 {
1905 gLog << supercutsPar[i] << ", ";
1906 }
1907 gLog << endl;
1908
1909
1910 // geometry is needed in MHHillas... classes
1911 MGeomCam *fGeom =
1912 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1913
1914 //-------------------------------------------
1915 // create the tasks which should be executed
1916 //
1917
1918 MReadMarsFile read("Events", filenameData);
1919 read.DisableAutoScheme();
1920
1921 TString fHilName = "MHillas";
1922 TString fHilNameExt = "MHillasExt";
1923 TString fHilNameSrc = "MHillasSrc";
1924 TString fImgParName = "MNewImagePar";
1925
1926
1927 //.......................................................................
1928 // calculation of hadroness for the supercuts
1929 // (=0.25 if fullfilled, =0.75 otherwise)
1930
1931 TString hadSCName = "HadSC";
1932 MCT1SupercutsCalc sccalc(fHilName, fHilNameSrc);
1933 sccalc.SetHadronnessName(hadSCName);
1934
1935
1936 //.......................................................................
1937
1938
1939 //MWriteRootFile write(outNameImage, "UPDATE");
1940 //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
1941
1942
1943 MWriteRootFile write(outNameImage, "RECREATE");
1944
1945 write.AddContainer("MRawRunHeader", "RunHeaders");
1946 write.AddContainer("MTime", "Events");
1947 write.AddContainer("MMcEvt", "Events");
1948 write.AddContainer("ThetaOrig", "Events");
1949 write.AddContainer("MSrcPosCam", "Events");
1950 write.AddContainer("MSigmabar", "Events");
1951 write.AddContainer("MHillas", "Events");
1952 write.AddContainer("MHillasExt", "Events");
1953 write.AddContainer("MHillasSrc", "Events");
1954 write.AddContainer("MNewImagePar", "Events");
1955
1956 write.AddContainer("HadRF", "Events");
1957 write.AddContainer(hadSCName, "Events");
1958
1959
1960 //-----------------------------------------------------------------
1961 // geometry is needed in MHHillas... classes
1962 MGeomCam *fGeom =
1963 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1964
1965 Float_t maxhadronness = 0.40;
1966 Float_t maxalpha = 20.0;
1967 Float_t maxdist = 10.0;
1968
1969 MFCT1SelFinal selfinalgh(fHilNameSrc);
1970 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1971 selfinalgh.SetHadronnessName(hadSCName);
1972 selfinalgh.SetName("SelFinalgh");
1973 MContinue contfinalgh(&selfinalgh);
1974 contfinalgh.SetName("ContSelFinalgh");
1975
1976 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
1977 fillhadsc.SetName("HhadSC");
1978
1979 MFCT1SelFinal selfinal(fHilNameSrc);
1980 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1981 selfinal.SetHadronnessName(hadSCName);
1982 selfinal.SetName("SelFinal");
1983 MContinue contfinal(&selfinal);
1984 contfinal.SetName("ContSelFinal");
1985
1986 TString mh3name = "abs(Alpha)";
1987 MBinning binsalphaabs("Binning"+mh3name);
1988 binsalphaabs.SetEdges(50, -2.0, 98.0);
1989
1990 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
1991 alphaabs.SetName(mh3name);
1992
1993 TH1 &alphahist = alphaabs->GetHist();
1994
1995 MFillH alpha(&alphaabs);
1996 alpha.SetName("FillAlphaAbs");
1997
1998 MFillH hfill1("MHHillas", fHilName);
1999 hfill1.SetName("HHillas");
2000
2001 MFillH hfill2("MHStarMap", fHilName);
2002 hfill2.SetName("HStarMap");
2003
2004 MFillH hfill3("MHHillasExt", fHilNameSrc);
2005 hfill3.SetName("HHillasExt");
2006
2007 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2008 hfill4.SetName("HHillasSrc");
2009
2010 MFillH hfill5("MHNewImagePar", fImgParName);
2011 hfill5.SetName("HNewImagePar");
2012
2013 //*****************************
2014 // entries in MParList
2015
2016 pliston.AddToList(&tliston);
2017 InitBinnings(&pliston);
2018
2019 pliston.AddToList(&supercuts);
2020
2021 pliston.AddToList(&binsalphaabs);
2022 pliston.AddToList(&alphaabs);
2023
2024 //*****************************
2025 // entries in MTaskList
2026
2027 tliston.AddToList(&read);
2028
2029 tliston.AddToList(&sccalc);
2030 tliston.AddToList(&fillhadsc);
2031
2032 tliston.AddToList(&write);
2033 tliston.AddToList(&contfinalgh);
2034
2035 tliston.AddToList(&alpha);
2036 tliston.AddToList(&hfill1);
2037 tliston.AddToList(&hfill2);
2038 tliston.AddToList(&hfill3);
2039 tliston.AddToList(&hfill4);
2040 tliston.AddToList(&hfill5);
2041
2042 tliston.AddToList(&contfinal);
2043
2044 //*****************************
2045
2046 //-------------------------------------------
2047 // Execute event loop
2048 //
2049 MProgressBar bar;
2050 MEvtLoop evtloop;
2051 evtloop.SetParList(&pliston);
2052 evtloop.SetProgressBar(&bar);
2053
2054 Int_t maxevents = -1;
2055 if ( !evtloop.Eventloop(maxevents) )
2056 return;
2057
2058 tliston.PrintStatistics(0, kTRUE);
2059
2060
2061 //-------------------------------------------
2062 // Display the histograms
2063 //
2064 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2065
2066 pliston.FindObject("MHHillas")->DrawClone();
2067 pliston.FindObject("MHHillasExt")->DrawClone();
2068 pliston.FindObject("MHHillasSrc")->DrawClone();
2069 pliston.FindObject("MHNewImagePar")->DrawClone();
2070 pliston.FindObject("MHStarMap")->DrawClone();
2071
2072 //-------------------------------------------
2073 // fit alpha distribution to get the number of excess events and
2074 // calculate significance of gamma signal in the alpha plot
2075
2076 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
2077 TH1 &alphaHist = absalpha->GetHist();
2078 alphaHist.SetXTitle("|alpha| [\\circ]");
2079 alphaHist.SetName("alpha-macro");
2080
2081 Double_t alphasig = 13.1;
2082 Double_t alphamin = 30.0;
2083 Double_t alphamax = 90.0;
2084 Int_t degree = 2;
2085 Double_t significance = -99.0;
2086 Bool_t drawpoly = kTRUE;
2087 Bool_t fitgauss = kTRUE;
2088 Bool_t print = kTRUE;
2089
2090 MHFindSignificance findsig;
2091 findsig.SetRebin(kTRUE);
2092 findsig.SetReduceDegree(kFALSE);
2093
2094 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
2095 alphasig, drawpoly, fitgauss, print);
2096 significance = findsig.GetSignificance();
2097 Float_t alphasi = findsig.GetAlphasi();
2098
2099 gLog << "For file '" << filenameData << "' : " << endl;
2100 gLog << "Significance of gamma signal after supercuts : "
2101 << significance << " (for |alpha| < " << alphasi << " degrees)"
2102 << endl;
2103
2104 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
2105
2106 //-------------------------------------------
2107
2108 DeleteBinnings(&pliston);
2109 }
2110
2111
2112 gLog << "Macro CT1Analysis : End of Job B_SC_UP" << endl;
2113 gLog << "=======================================================" << endl;
2114 }
2115 //---------------------------------------------------------------------
2116
2117
2118 //---------------------------------------------------------------------
2119 // Job C
2120 //======
2121
2122 // - read ON1 and MC1 data files
2123 // which should have been updated to contain the hadronnesses
2124 // for the method of Random Forest and for the SUPERCUTS
2125 // - produce Neyman-Pearson plots
2126
2127 if (JobC)
2128 {
2129 gLog << "=====================================================" << endl;
2130 gLog << "Macro CT1Analysis : Start of Job C" << endl;
2131
2132 gLog << "" << endl;
2133 gLog << "Macro CT1Analysis : JobC = "
2134 << (JobC ? "kTRUE" : "kFALSE") << endl;
2135
2136
2137 // name of input data file
2138 TString filenameData = outPath;
2139 filenameData += "ON";
2140 filenameData += "3.root";
2141 gLog << "filenameData = " << filenameData << endl;
2142
2143 // name of input MC file
2144 TString filenameMC = outPath;
2145 filenameMC += "MC";
2146 filenameMC += "3.root";
2147 gLog << "filenameMC = " << filenameMC << endl;
2148
2149
2150 //-----------------------------------------------------------------
2151
2152 MTaskList tliston;
2153 MParList pliston;
2154
2155
2156 // geometry is needed in MHHillas... classes
2157 MGeomCam *fGeom =
2158 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2159
2160 //-------------------------------------------
2161 // create the tasks which should be executed
2162 //
2163
2164 MReadMarsFile read("Events", filenameMC);
2165 read.AddFile(filenameData);
2166 read.DisableAutoScheme();
2167
2168
2169 //.......................................................................
2170 // names of hadronness containers
2171
2172 //TString hadNNName = "HadNN";
2173 TString hadSCName = "HadSC";
2174 TString hadRFName = "HadRF";
2175
2176 //.......................................................................
2177
2178
2179 TString fHilName = "MHillas";
2180 TString fHilNameExt = "MHillasExt";
2181 TString fHilNameSrc = "MHillasSrc";
2182 TString fImgParName = "MNewImagePar";
2183
2184 Float_t maxhadronness = 0.40;
2185 Float_t maxalpha = 20.0;
2186 Float_t maxdist = 10.0;
2187
2188 MFCT1SelFinal selfinalgh(fHilNameSrc);
2189 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2190 selfinalgh.SetHadronnessName(hadSCName);
2191 selfinalgh.SetName("SelFinalgh");
2192 MContinue contfinalgh(&selfinalgh);
2193 contfinalgh.SetName("ContSelFinalgh");
2194
2195 //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
2196 //fillhadnn.SetName("HhadNN");
2197 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
2198 fillhadsc.SetName("HhadSC");
2199 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
2200 fillhadrf.SetName("HhadRF");
2201
2202 MFCT1SelFinal selfinal(fHilNameSrc);
2203 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2204 selfinal.SetHadronnessName(hadSCName);
2205 selfinal.SetName("SelFinal");
2206 MContinue contfinal(&selfinal);
2207 contfinal.SetName("ContSelFinal");
2208
2209
2210 MFillH hfill1("MHHillas", fHilName);
2211 hfill1.SetName("HHillas");
2212
2213 MFillH hfill2("MHStarMap", fHilName);
2214 hfill2.SetName("HStarMap");
2215
2216 MFillH hfill3("MHHillasExt", fHilNameSrc);
2217 hfill3.SetName("HHillasExt");
2218
2219 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2220 hfill4.SetName("HHillasSrc");
2221
2222 MFillH hfill5("MHNewImagePar", fImgParName);
2223 hfill5.SetName("HNewImagePar");
2224
2225
2226 //*****************************
2227 // entries in MParList
2228
2229 pliston.AddToList(&tliston);
2230 InitBinnings(&pliston);
2231
2232
2233 //*****************************
2234 // entries in MTaskList
2235
2236 tliston.AddToList(&read);
2237
2238 //tliston.AddToList(&fillhadnn);
2239 tliston.AddToList(&fillhadsc);
2240 tliston.AddToList(&fillhadrf);
2241
2242 tliston.AddToList(&contfinalgh);
2243 tliston.AddToList(&hfill1);
2244 tliston.AddToList(&hfill2);
2245 tliston.AddToList(&hfill3);
2246 tliston.AddToList(&hfill4);
2247 tliston.AddToList(&hfill5);
2248
2249 tliston.AddToList(&contfinal);
2250
2251 //*****************************
2252
2253 //-------------------------------------------
2254 // Execute event loop
2255 //
2256 MProgressBar bar;
2257 MEvtLoop evtloop;
2258 evtloop.SetParList(&pliston);
2259 evtloop.SetProgressBar(&bar);
2260
2261 Int_t maxevents = -1;
2262 //Int_t maxevents = 35000;
2263 if ( !evtloop.Eventloop(maxevents) )
2264 return;
2265
2266 tliston.PrintStatistics(0, kTRUE);
2267
2268
2269 //-------------------------------------------
2270 // Display the histograms
2271 //
2272
2273 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
2274 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2275 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2276
2277 pliston.FindObject("MHHillas")->DrawClone();
2278 pliston.FindObject("MHHillasExt")->DrawClone();
2279 pliston.FindObject("MHHillasSrc")->DrawClone();
2280 pliston.FindObject("MHNewImagePar")->DrawClone();
2281 pliston.FindObject("MHStarMap")->DrawClone();
2282
2283 DeleteBinnings(&pliston);
2284
2285 gLog << "Macro CT1Analysis : End of Job C" << endl;
2286 gLog << "===================================================" << endl;
2287 }
2288
2289
2290 //---------------------------------------------------------------------
2291 // Job D
2292 //======
2293
2294 // - select g/h separation method XX
2295 // - read ON2 (or MC2) root file
2296 // - apply cuts in hadronness
2297 // - make plots
2298
2299
2300 if (JobD)
2301 {
2302 gLog << "=====================================================" << endl;
2303 gLog << "Macro CT1Analysis : Start of Job D" << endl;
2304
2305 gLog << "" << endl;
2306 gLog << "Macro CT1Analysis : JobD = "
2307 << (JobD ? "kTRUE" : "kFALSE") << endl;
2308
2309
2310 // type of data to be analysed
2311 TString typeData = "ON";
2312 //TString typeData = "OFF";
2313 //TString typeData = "MC";
2314 gLog << "typeData = " << typeData << endl;
2315
2316 TString ext = "3.root";
2317
2318
2319 //------------------------------
2320 // selection of g/h separation method
2321 // and definition of final selections
2322
2323 //TString XX("SC");
2324 TString XX("RF");
2325 TString fhadronnessName("Had");
2326 fhadronnessName += XX;
2327 gLog << "fhadronnessName = " << fhadronnessName << endl;
2328
2329 // maximum values of the hadronness, |ALPHA| and DIST
2330 Float_t maxhadronness = 0.233;
2331 Float_t maxalpha = 20.0;
2332 Float_t maxdist = 10.0;
2333 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2334 << maxhadronness << ", " << maxalpha << ", "
2335 << maxdist << endl;
2336
2337
2338 //------------------------------
2339 // name of data file to be analysed
2340 TString filenameData(outPath);
2341 filenameData += typeData;
2342 filenameData += ext;
2343 gLog << "filenameData = " << filenameData << endl;
2344
2345
2346
2347 //*************************************************************************
2348 //
2349 // Analyse the data
2350 //
2351
2352 MTaskList tliston;
2353 MParList pliston;
2354
2355 // geometry is needed in MHHillas... classes
2356 MGeomCam *fGeom =
2357 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2358
2359
2360 TString fHilName = "MHillas";
2361 TString fHilNameExt = "MHillasExt";
2362 TString fHilNameSrc = "MHillasSrc";
2363 TString fImgParName = "MNewImagePar";
2364
2365 //-------------------------------------------
2366 // create the tasks which should be executed
2367 //
2368
2369 MReadMarsFile read("Events", filenameData);
2370 read.DisableAutoScheme();
2371
2372
2373 //-----------------------------------------------------------------
2374 // geometry is needed in MHHillas... classes
2375 MGeomCam *fGeom =
2376 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2377
2378 MFCT1SelFinal selfinalgh(fHilNameSrc);
2379 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2380 selfinalgh.SetHadronnessName(fhadronnessName);
2381 selfinalgh.SetName("SelFinalgh");
2382 MContinue contfinalgh(&selfinalgh);
2383 contfinalgh.SetName("ContSelFinalgh");
2384
2385 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
2386 fillhadsc.SetName("HhadSC");
2387 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
2388 fillhadrf.SetName("HhadRF");
2389
2390 TString mh3name = "abs(Alpha)";
2391 MBinning binsalphaabs("Binning"+mh3name);
2392 binsalphaabs.SetEdges(50, -2.0, 98.0);
2393
2394 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
2395 alphaabs.SetName(mh3name);
2396
2397 TH1 &alphahist = alphaabs->GetHist();
2398
2399 MFillH alpha(&alphaabs);
2400 alpha.SetName("FillAlphaAbs");
2401
2402 MFillH hfill1("MHHillas", fHilName);
2403 hfill1.SetName("HHillas");
2404
2405 MFillH hfill2("MHStarMap", fHilName);
2406 hfill2.SetName("HStarMap");
2407
2408 MFillH hfill3("MHHillasExt", fHilNameSrc);
2409 hfill3.SetName("HHillasExt");
2410
2411 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2412 hfill4.SetName("HHillasSrc");
2413
2414 MFillH hfill5("MHNewImagePar", fImgParName);
2415 hfill5.SetName("HNewImagePar");
2416
2417 MFCT1SelFinal selfinal(fHilNameSrc);
2418 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2419 selfinal.SetHadronnessName(fhadronnessName);
2420 selfinal.SetName("SelFinal");
2421 MContinue contfinal(&selfinal);
2422 contfinal.SetName("ContSelFinal");
2423
2424
2425 //*****************************
2426 // entries in MParList
2427
2428 pliston.AddToList(&tliston);
2429 InitBinnings(&pliston);
2430 pliston.AddToList(&binsalphaabs);
2431 pliston.AddToList(&alphaabs);
2432
2433 //*****************************
2434 // entries in MTaskList
2435
2436 tliston.AddToList(&read);
2437
2438 tliston.AddToList(&contfinalgh);
2439
2440 tliston.AddToList(&fillhadsc);
2441 tliston.AddToList(&fillhadrf);
2442
2443 tliston.AddToList(&alpha);
2444 tliston.AddToList(&hfill1);
2445 tliston.AddToList(&hfill2);
2446 tliston.AddToList(&hfill3);
2447 tliston.AddToList(&hfill4);
2448 tliston.AddToList(&hfill5);
2449
2450 tliston.AddToList(&contfinal);
2451
2452 //*****************************
2453
2454 //-------------------------------------------
2455 // Execute event loop
2456 //
2457 MProgressBar bar;
2458 MEvtLoop evtloop;
2459 evtloop.SetParList(&pliston);
2460 evtloop.SetProgressBar(&bar);
2461
2462 Int_t maxevents = -1;
2463 //Int_t maxevents = 10000;
2464 if ( !evtloop.Eventloop(maxevents) )
2465 return;
2466
2467 tliston.PrintStatistics(0, kTRUE);
2468
2469
2470 //-------------------------------------------
2471 // Display the histograms
2472 //
2473
2474 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2475 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2476
2477 pliston.FindObject("MHHillas")->DrawClone();
2478 pliston.FindObject("MHHillasExt")->DrawClone();
2479 pliston.FindObject("MHHillasSrc")->DrawClone();
2480 pliston.FindObject("MHNewImagePar")->DrawClone();
2481 pliston.FindObject("MHStarMap")->DrawClone();
2482
2483
2484 //-------------------------------------------
2485
2486 // fit alpha distribution to get the number of excess events and
2487 // calculate significance of gamma signal in the alpha plot
2488
2489 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
2490 TH1 &alphaHist = absalpha->GetHist();
2491 alphaHist.SetXTitle("|alpha| [\\circ]");
2492 alphaHist.SetName("alpha-JobD");
2493
2494 Double_t alphasig = 13.1;
2495 Double_t alphamin = 30.0;
2496 Double_t alphamax = 90.0;
2497 Int_t degree = 2;
2498 Double_t significance = -99.0;
2499 Bool_t drawpoly = kTRUE;
2500 Bool_t fitgauss = kTRUE;
2501 Bool_t print = kTRUE;
2502
2503 MHFindSignificance findsig;
2504 findsig.SetRebin(kTRUE);
2505 findsig.SetReduceDegree(kFALSE);
2506
2507 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
2508 alphasig, drawpoly, fitgauss, print);
2509 significance = findsig.GetSignificance();
2510 Float_t alphasi = findsig.GetAlphasi();
2511
2512 gLog << "For file '" << filenameData << "' : " << endl;
2513 gLog << "Significance of gamma signal after supercuts : "
2514 << significance << " (for |alpha| < " << alphasi << " degrees)"
2515 << endl;
2516
2517 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
2518
2519 //-------------------------------------------
2520
2521
2522 DeleteBinnings(&pliston);
2523
2524 gLog << "Macro CT1Analysis : End of Job D" << endl;
2525 gLog << "=======================================================" << endl;
2526 }
2527 //---------------------------------------------------------------------
2528
2529
2530
2531
2532 //---------------------------------------------------------------------
2533 // Job E_XX
2534 //=========
2535
2536 // - select g/h separation method XX
2537 // - read MC_XX2.root file
2538   // - calculate eff. collection area
2539 // - read ON_XX2.root file
2540 // - apply final cuts
2541 // - calculate flux
2542 // - write root file for ON data after final cuts (ON_XX3.root))
2543
2544
2545 if (JobE_XX)
2546 {
2547 gLog << "=====================================================" << endl;
2548 gLog << "Macro CT1Analysis : Start of Job E_XX" << endl;
2549
2550 gLog << "" << endl;
2551 gLog << "Macro CT1Analysis : JobE_XX, CCollArea, OEEst, WEX = "
2552 << (JobE_XX ? "kTRUE" : "kFALSE") << ", "
2553 << (CCollArea?"kTRUE" : "kFALSE") << ", "
2554 << (OEEst ? "kTRUE" : "kFALSE") << ", "
2555 << (WEX ? "kTRUE" : "kFALSE") << endl;
2556
2557
2558 // type of data to be analysed
2559 //TString typeData = "ON";
2560 //TString typeData = "OFF";
2561 TString typeData = "MC";
2562 gLog << "typeData = " << typeData << endl;
2563
2564 TString typeMC = "MC";
2565 TString ext = "3.root";
2566 TString extout = "4.root";
2567
2568 //------------------------------
2569 // selection of g/h separation method
2570 // and definition of final selections
2571
2572 //TString XX("SC");
2573 TString XX("RF");
2574 TString fhadronnessName("Had");
2575 fhadronnessName += XX;
2576 gLog << "fhadronnessName = " << fhadronnessName << endl;
2577
2578 // maximum values of the hadronness, |ALPHA| and DIST
2579 Float_t maxhadronness = 0.23;
2580 Float_t maxalpha = 20.0;
2581 Float_t maxdist = 10.0;
2582 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2583 << maxhadronness << ", " << maxalpha << ", "
2584 << maxdist << endl;
2585
2586 //------------------------------
2587 // name of MC file to be used for optimizing the energy estimator
2588 TString filenameOpt(outPath);
2589 filenameOpt += typeMC;
2590 filenameOpt += ext;
2591 gLog << "filenameOpt = " << filenameOpt << endl;
2592
2593 //------------------------------
2594 // name of file containing the parameters of the energy estimator
2595 TString energyParName(outPath);
2596 energyParName += "energyest_";
2597 energyParName += XX;
2598 energyParName += ".root";
2599 gLog << "energyParName = " << energyParName << endl;
2600
2601 //------------------------------
2602 // name of MC file to be used for calculating the eff. collection areas
2603 TString filenameArea(outPath);
2604 filenameArea += typeMC;
2605 filenameArea += ext;
2606 gLog << "filenameArea = " << filenameArea << endl;
2607
2608 //------------------------------
2609 // name of file containing the eff. collection areas
2610 TString collareaName(outPath);
2611 collareaName += "area_";
2612 collareaName += XX;
2613 collareaName += ".root";
2614 gLog << "collareaName = " << collareaName << endl;
2615
2616 //------------------------------
2617 // name of data file to be analysed
2618 TString filenameData(outPath);
2619 filenameData += typeData;
2620 filenameData += ext;
2621 gLog << "filenameData = " << filenameData << endl;
2622
2623 //------------------------------
2624 // name of output data file (after the final cuts)
2625 TString filenameDataout(outPath);
2626 filenameDataout += typeData;
2627 filenameDataout += "_";
2628 filenameDataout += XX;
2629 filenameDataout += extout;
2630 gLog << "filenameDataout = " << filenameDataout << endl;
2631
2632 //------------------------------
2633 // name of file containing histograms for flux calculastion
2634 TString filenameResults(outPath);
2635 filenameResults += typeData;
2636 filenameResults += "Results_";
2637 filenameResults += XX;
2638 filenameResults += extout;
2639 gLog << "filenameResults = " << filenameResults << endl;
2640
2641
2642 //====================================================================
2643
2644 MHMcCT1CollectionArea collarea;
2645 collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
2646
2647 MParList parlist;
2648 InitBinnings(&parlist);
2649
2650 if (CCollArea)
2651 {
2652 gLog << "-----------------------------------------------" << endl;
2653 gLog << "Start calculation of effective collection areas" << endl;
2654
2655
2656 MTaskList tasklist;
2657
2658 //---------------------------------------
2659 // Setup the tasks to be executed
2660 //
2661 MReadMarsFile reader("Events", filenameArea);
2662 reader.DisableAutoScheme();
2663
2664 MFCT1SelFinal cuthadrons;
2665 cuthadrons.SetHadronnessName(fhadronnessName);
2666 cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
2667
2668 MContinue conthadrons(&cuthadrons);
2669
2670
2671 MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
2672 filler.SetName("CollectionArea");
2673
2674 //********************************
2675 // entries in MParList
2676
2677 parlist.AddToList(&tasklist);
2678
2679 parlist.AddToList(&collarea);
2680
2681 //********************************
2682 // entries in MTaskList
2683
2684 tasklist.AddToList(&reader);
2685 tasklist.AddToList(&conthadrons);
2686 tasklist.AddToList(&filler);
2687
2688 //********************************
2689
2690 //-----------------------------------------
2691 // Execute event loop
2692 //
2693 MEvtLoop magic;
2694 magic.SetParList(&parlist);
2695
2696 MProgressBar bar;
2697 magic.SetProgressBar(&bar);
2698 if (!magic.Eventloop())
2699 return;
2700
2701 tasklist.PrintStatistics(0, kTRUE);
2702
2703 // Calculate effective collection areas
2704 // and display the histograms
2705 //
2706 //MHMcCT1CollectionArea *collarea =
2707 // (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
2708 collarea.CalcEfficiency();
2709 collarea.DrawClone();
2710
2711
2712
2713 //---------------------------------------------
2714 // Write histograms to a file
2715 //
2716
2717 TFile f(collareaName, "RECREATE");
2718 //collarea.GetHist()->Write();
2719 //collarea.GetHAll()->Write();
2720 //collarea.GetHSel()->Write();
2721 collarea.Write();
2722
2723 f.Close();
2724
2725 gLog << "Collection area plots written onto file " << collareaName << endl;
2726
2727 gLog << "Calculation of effective collection areas done" << endl;
2728 gLog << "-----------------------------------------------" << endl;
2729 //------------------------------------------------------------------
2730 }
2731
2732 if (!CCollArea)
2733 {
2734 gLog << "-----------------------------------------------" << endl;
2735 gLog << "Read in effective collection areas from file "
2736 << collareaName << endl;
2737
2738 TFile collfile(collareaName);
2739 collfile.ls();
2740 collarea.Read("MHMcCT1CollectionArea");
2741 collarea.DrawClone();
2742
2743 gLog << "Effective collection areas were read in from file "
2744 << collareaName << endl;
2745 gLog << "-----------------------------------------------" << endl;
2746 }
2747
2748
2749 // save binnings for call to CT1EEst
2750 MBinning *binsE = (MBinning*)parlist.FindObject("BinningE");
2751 if (!binsE)
2752 {
2753 gLog << "Object 'BinningE' not found in MParList" << endl;
2754 return;
2755 }
2756 MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
2757 if (!binsTheta)
2758 {
2759 gLog << "Object 'BinningTheta' not found in MParList" << endl;
2760 return;
2761 }
2762
2763 //-------------------------------------
2764 TString fHilName = "MHillas";
2765 TString fHilNameExt = "MHillasExt";
2766 TString fHilNameSrc = "MHillasSrc";
2767 TString fImgParName = "MNewImagePar";
2768
2769
2770 if (OEEst)
2771 {
2772 //===========================================================
2773 //
2774 // Optimization of energy estimator
2775 //
2776 gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
2777
2778 TString inpath("");
2779 TString outpath("");
2780 Int_t howMany = 2000;
2781 CT1EEst(inpath, filenameOpt, outpath, energyParName,
2782 fHilName, fHilNameSrc, fhadronnessName,
2783 howMany, maxhadronness, maxalpha, maxdist,
2784 binsE, binsTheta);
2785 gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
2786 }
2787
2788 if (WEX)
2789 {
2790 //-----------------------------------------------------------
2791 //
2792 // Read in parameters of energy estimator ("MMcEnergyEst")
2793 // and migration matrix ("MHMcEnergyMigration")
2794 //
2795 gLog << "================================================================"
2796 << endl;
2797 gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
2798 << energyParName << "'" << endl;
2799 TFile enparam(energyParName);
2800 enparam.ls();
2801 MMcEnergyEst mcest("MMcEnergyEst");
2802 mcest.Read("MMcEnergyEst");
2803
2804 //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
2805 gLog << "Parameters of energy estimator were read in" << endl;
2806
2807
2808 gLog << "Read in Migration matrix" << endl;
2809
2810 MHMcEnergyMigration mighiston("MHMcEnergyMigration");
2811 mighiston.Read("MHMcEnergyMigration");
2812 //MHMcEnergyMigration &mighiston =
2813 // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
2814
2815 gLog << "Migration matrix was read in" << endl;
2816
2817
2818 TArrayD parA(mcest.GetNumCoeffA());
2819 TArrayD parB(mcest.GetNumCoeffB());
2820 for (Int_t i=0; i<parA.GetSize(); i++)
2821 parA[i] = mcest.GetCoeff(i);
2822 for (Int_t i=0; i<parB.GetSize(); i++)
2823 parB[i] = mcest.GetCoeff( i+parA.GetSize() );
2824
2825 //*************************************************************************
2826 //
2827 // Analyse the data
2828 //
2829 gLog << "============================================================"
2830 << endl;
2831 gLog << "Analyse the data" << endl;
2832
2833 MTaskList tliston;
2834 MParList pliston;
2835
2836 // geometry is needed in MHHillas... classes
2837 MGeomCam *fGeom =
2838 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2839
2840
2841 //-------------------------------------------
2842 // create the tasks which should be executed
2843 //
2844
2845 MReadMarsFile read("Events", filenameData);
2846 read.DisableAutoScheme();
2847
2848 //.......................................................................
2849
2850 gLog << "CT1Analysis.C : write root file '" << filenameDataout
2851 << "'" << endl;
2852
2853 //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
2854
2855
2856 MWriteRootFile write(filenameDataout, "RECREATE");
2857
2858 write.AddContainer("MRawRunHeader", "RunHeaders");
2859 write.AddContainer("MTime", "Events");
2860 write.AddContainer("MMcEvt", "Events");
2861 write.AddContainer("ThetaOrig", "Events");
2862 write.AddContainer("MSrcPosCam", "Events");
2863 write.AddContainer("MSigmabar", "Events");
2864 write.AddContainer("MHillas", "Events");
2865 write.AddContainer("MHillasExt", "Events");
2866 write.AddContainer("MHillasSrc", "Events");
2867 write.AddContainer("MNewImagePar", "Events");
2868
2869 //write.AddContainer("HadNN", "Events");
2870 write.AddContainer("HadSC", "Events");
2871 write.AddContainer("HadRF", "Events");
2872
2873 write.AddContainer("MEnergyEst", "Events");
2874
2875
2876 //-----------------------------------------------------------------
2877 // geometry is needed in MHHillas... classes
2878 MGeomCam *fGeom =
2879 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2880
2881 MFCT1SelFinal selfinalgh(fHilNameSrc);
2882 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2883 selfinalgh.SetHadronnessName(fhadronnessName);
2884 selfinalgh.SetName("SelFinalgh");
2885 MContinue contfinalgh(&selfinalgh);
2886 contfinalgh.SetName("ContSelFinalgh");
2887
2888 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
2889 //fillhadnn.SetName("HhadNN");
2890 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
2891 fillhadsc.SetName("HhadSC");
2892 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
2893 fillhadrf.SetName("HhadRF");
2894
2895 //---------------------------
2896 // calculate estimated energy
2897
2898 MEnergyEstParam eeston(fHilName);
2899 eeston.Add(fHilNameSrc);
2900
2901 eeston.SetCoeffA(parA);
2902 eeston.SetCoeffB(parB);
2903
2904 //---------------------------
2905 // calculate estimated energy using Daniel's parameters
2906
2907 //MEnergyEstParamDanielMkn421 eeston(fHilName);
2908 //eeston.Add(fHilNameSrc);
2909 //eeston.SetCoeffA(parA);
2910 //eeston.SetCoeffB(parB);
2911
2912
2913 //---------------------------
2914
2915
2916 MFillH hfill1("MHHillas", fHilName);
2917 hfill1.SetName("HHillas");
2918
2919 MFillH hfill2("MHStarMap", fHilName);
2920 hfill2.SetName("HStarMap");
2921
2922 MFillH hfill3("MHHillasExt", fHilNameSrc);
2923 hfill3.SetName("HHillasExt");
2924
2925 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2926 hfill4.SetName("HHillasSrc");
2927
2928 MFillH hfill5("MHNewImagePar", fImgParName);
2929 hfill5.SetName("HNewImagePar");
2930
2931 //---------------------------
2932 // new from Robert
2933
2934 MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
2935 hfill6.SetName("HTimeDiffTheta");
2936
2937 MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
2938 hfill6a.SetName("HTimeDiffTime");
2939
2940 MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
2941 hfill7.SetName("HAlphaEnergyTheta");
2942
2943 MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
2944 hfill7a.SetName("HAlphaEnergyTime");
2945
2946 MFillH hfill7b("MHThetabarTime", fHilNameSrc);
2947 hfill7b.SetName("HThetabarTime");
2948
2949 MFillH hfill7c("MHEnergyTime", "MMcEvt");
2950 hfill7c.SetName("HEnergyTime");
2951
2952
2953 //---------------------------
2954
2955 MFCT1SelFinal selfinal(fHilNameSrc);
2956 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2957 selfinal.SetHadronnessName(fhadronnessName);
2958 selfinal.SetName("SelFinal");
2959 MContinue contfinal(&selfinal);
2960 contfinal.SetName("ContSelFinal");
2961
2962
2963 //*****************************
2964 // entries in MParList
2965
2966 pliston.AddToList(&tliston);
2967 InitBinnings(&pliston);
2968
2969
2970 //*****************************
2971 // entries in MTaskList
2972
2973 tliston.AddToList(&read);
2974
2975 // robert
2976 tliston.AddToList(&hfill6); //timediff
2977 tliston.AddToList(&hfill6a); //timediff
2978
2979 tliston.AddToList(&contfinalgh);
2980 tliston.AddToList(&eeston);
2981
2982 tliston.AddToList(&write);
2983
2984 //tliston.AddToList(&fillhadnn);
2985 tliston.AddToList(&fillhadsc);
2986 tliston.AddToList(&fillhadrf);
2987
2988 tliston.AddToList(&hfill1);
2989 tliston.AddToList(&hfill2);
2990 tliston.AddToList(&hfill3);
2991 tliston.AddToList(&hfill4);
2992 tliston.AddToList(&hfill5);
2993
2994 //robert
2995 tliston.AddToList(&hfill7);
2996 tliston.AddToList(&hfill7a);
2997 tliston.AddToList(&hfill7b);
2998 tliston.AddToList(&hfill7c);
2999
3000 tliston.AddToList(&contfinal);
3001
3002 //*****************************
3003
3004 //-------------------------------------------
3005 // Execute event loop
3006 //
3007 MProgressBar bar;
3008 MEvtLoop evtloop;
3009 evtloop.SetParList(&pliston);
3010 evtloop.SetProgressBar(&bar);
3011
3012 Int_t maxevents = -1;
3013 if ( !evtloop.Eventloop(maxevents) )
3014 return;
3015
3016 tliston.PrintStatistics(0, kTRUE);
3017
3018
3019 //-------------------------------------------
3020 // Display the histograms
3021 //
3022
3023 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
3024
3025 gLog << "before hadRF" << endl;
3026 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
3027
3028 gLog << "before hadSC" << endl;
3029 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
3030
3031 gLog << "before MHHillas" << endl;
3032 pliston.FindObject("MHHillas")->DrawClone();
3033
3034 gLog << "before MHHillasExt" << endl;
3035 pliston.FindObject("MHHillasExt")->DrawClone();
3036
3037 gLog << "before MHHillasSrc" << endl;
3038 pliston.FindObject("MHHillasSrc")->DrawClone();
3039
3040 gLog << "before MHNewImagePar" << endl;
3041 pliston.FindObject("MHNewImagePar")->DrawClone();
3042
3043 gLog << "before MHStarMap" << endl;
3044 pliston.FindObject("MHStarMap")->DrawClone();
3045
3046 gLog << "before DeleteBinnings" << endl;
3047
3048 DeleteBinnings(&pliston);
3049
3050 gLog << "before Robert's code" << endl;
3051
3052
3053//rwagner write all relevant histograms onto a file
3054
3055 if (WRobert)
3056 {
3057 gLog << "=======================================================" << endl;
3058 gLog << "Write results onto file '" << filenameResults << "'" << endl;
3059
3060 TFile outfile(filenameResults,"recreate");
3061
3062 MHHillasSrc* hillasSrc =
3063 (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
3064 TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
3065 alphaHist->Write();
3066 gLog << "Alpha plot has been written out" << endl;
3067
3068
3069 MHAlphaEnergyTheta* aetH =
3070 (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
3071 TH3D* aetHist = (TH3D*)(aetH->GetHist());
3072 aetHist->SetName("aetHist");
3073 aetHist->Write();
3074 gLog << "AlphaEnergyTheta plot has been written out" << endl;
3075
3076 MHAlphaEnergyTime* aetH2 =
3077 (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
3078 TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
3079 aetHist2->SetName("aetimeHist");
3080// aetHist2->DrawClone();
3081 aetHist2->Write();
3082 gLog << "AlphaEnergyTime plot has been written out" << endl;
3083
3084 MHThetabarTime* tbt =
3085 (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
3086 TProfile* tbtHist = (TProfile*)(tbt->GetHist());
3087 tbtHist->SetName("tbtHist");
3088 tbtHist->Write();
3089 gLog << "ThetabarTime plot has been written out" << endl;
3090
3091 MHEnergyTime* ent =
3092 (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
3093 TH2D* entHist = (TH2D*)(ent->GetHist());
3094 entHist->SetName("entHist");
3095 entHist->Write();
3096 gLog << "EnergyTime plot has been written out" << endl;
3097
3098 MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
3099 TH2D* timeHist = (TH2D*)(time->GetHist());
3100 timeHist->SetName("MHTimeDiffTheta");
3101 timeHist->SetTitle("Time diffs");
3102 timeHist->Write();
3103 gLog << "TimeDiffTheta plot has been written out" << endl;
3104
3105
3106 MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
3107 TH2D* timeHist2 = (TH2D*)(time2->GetHist());
3108 timeHist2->SetName("MHTimeDiffTime");
3109 timeHist2->SetTitle("Time diffs");
3110 timeHist2->Write();
3111 gLog << "TimeDiffTime plot has been written out" << endl;
3112
3113//rwagner write also collareas to same file
3114 collarea->GetHist()->Write();
3115 collarea->GetHAll()->Write();
3116 collarea->GetHSel()->Write();
3117 gLog << "Effective collection areas have been written out" << endl;
3118
3119//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
3120//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
3121
3122 gLog << "before closing outfile" << endl;
3123
3124 //outfile.Close();
3125 gLog << "Results were written onto file '" << filenameResults
3126 << "'" << endl;
3127 gLog << "=======================================================" << endl;
3128 }
3129
3130 }
3131
3132 gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
3133 gLog << "=======================================================" << endl;
3134 }
3135 //---------------------------------------------------------------------
3136
3137}
3138
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
Note: See TracBrowser for help on using the repository browser.