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

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