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

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