source: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C@ 2445

Last change on this file since 2445 was 2437, checked in by wittek, 23 years ago
*** empty log message ***
File size: 72.2 KB
Line 
1
2#include "CT1EgyEst.C"
3
4
5void InitBinnings(MParList *plist)
6{
7 gLog << "InitBinnings" << endl;
8
9 MBinning *binse = new MBinning("BinningE");
10 //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
11 binse->SetEdges(30, 2, 5);
12 plist->AddToList(binse);
13
14 MBinning *binssize = new MBinning("BinningSize");
15 binssize->SetEdgesLog(50, 10, 1.0e5);
16 plist->AddToList(binssize);
17
18 MBinning *binsdistc = new MBinning("BinningDist");
19 binsdistc->SetEdges(50, 0, 1.4);
20 plist->AddToList(binsdistc);
21
22 MBinning *binswidth = new MBinning("BinningWidth");
23 binswidth->SetEdges(50, 0, 1.0);
24 plist->AddToList(binswidth);
25
26 MBinning *binslength = new MBinning("BinningLength");
27 binslength->SetEdges(50, 0, 1.0);
28 plist->AddToList(binslength);
29
30 MBinning *binsalpha = new MBinning("BinningAlpha");
31 binsalpha->SetEdges(100, -100, 100);
32 plist->AddToList(binsalpha);
33
34 MBinning *binsasym = new MBinning("BinningAsym");
35 binsasym->SetEdges(50, -1.5, 1.5);
36 plist->AddToList(binsasym);
37
38 MBinning *binsm3l = new MBinning("BinningM3Long");
39 binsm3l->SetEdges(50, -1.5, 1.5);
40 plist->AddToList(binsm3l);
41
42 MBinning *binsm3t = new MBinning("BinningM3Trans");
43 binsm3t->SetEdges(50, -1.5, 1.5);
44 plist->AddToList(binsm3t);
45
46
47 //.....
48 MBinning *binsb = new MBinning("BinningSigmabar");
49 binsb->SetEdges( 100, 0.0, 5.0);
50 plist->AddToList(binsb);
51
52 MBinning *binth = new MBinning("BinningTheta");
53 Double_t yedge[9] =
54 {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
55 TArrayD yed(9,yedge);
56 binth->SetEdges(yed);
57 plist->AddToList(binth);
58
59 MBinning *bincosth = new MBinning("BinningCosTheta");
60 Double_t zedge[9];
61 for (Int_t i=0; i<9; i++)
62 {
63 zedge[8-i] = cos(yedge[i]/kRad2Deg);
64 }
65 TArrayD zed(9,zedge);
66 bincosth->SetEdges(zed);
67 plist->AddToList(bincosth);
68
69 MBinning *binsdiff = new MBinning("BinningDiffsigma2");
70 binsdiff->SetEdges(100, -5.0, 20.0);
71 plist->AddToList(binsdiff);
72}
73
74void DeleteBinnings(MParList *plist)
75{
76 gLog << "DeleteBinnings" << endl;
77
78 if (!plist)
79 {
80 gLog << "Deletebinnins : MParlist no longer existing" << endl;
81 return;
82 }
83
84 TObject *bin;
85
86 bin = plist->FindObject("BinningSize");
87 if (bin) delete bin;
88
89 bin = plist->FindObject("BinningDist");
90 if (bin) delete bin;
91
92 bin = plist->FindObject("BinningWidth");
93 if (bin) delete bin;
94
95 bin = plist->FindObject("BinningLength");
96 if (bin) delete bin;
97
98 bin = plist->FindObject("BinningAlpha");
99 if (bin) delete bin;
100
101 bin = plist->FindObject("BinningAsym");
102 if (bin) delete bin;
103
104 bin = plist->FindObject("BinningM3Long");
105 if (bin) delete bin;
106
107 bin = plist->FindObject("BinningM3Trans");
108 if (bin) delete bin;
109
110 //.....
111 bin = plist->FindObject("BinningSigmabar");
112 if (bin) delete bin;
113
114 bin = plist->FindObject("BinningTheta");
115 if (bin) delete bin;
116
117 bin = plist->FindObject("BinningCosTheta");
118 if (bin) delete bin;
119
120 bin = plist->FindObject("BinningDiffsigma2");
121 if (bin) delete bin;
122
123 gLog << "exit DeleteBinnings" << endl;
124}
125
126//************************************************************************
127void ONOFFCT1Analysis()
128{
129 gLog.SetNoColors();
130
131 if (gRandom)
132 delete gRandom;
133 gRandom = new TRandom3(0);
134
135 //-----------------------------------------------
136 const char *offfile = "~magican/ct1test/wittek/offdata.preproc";
137
138 //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc";
139 const char *onfile = "~magican/ct1test/wittek/mkn421_00-01";
140
141 const char *mcfile = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc";
142 //-----------------------------------------------
143
144 // path for input for Mars
145 TString inPath = "~magican/ct1test/wittek/marsoutput/optionD/";
146
147 // path for output from Mars
148 TString outPath = "~magican/ct1test/wittek/marsoutput/optionD/";
149
150 //-----------------------------------------------
151
152 //TEnv env("macros/CT1env.rc");
153 //Bool_t printEnv = kFALSE;
154
155 //************************************************************************
156
157 // Job A :
158 // - produce MHSigmaTheta plots for ON and OFF data
159 // - write out (or read in) these MHSigmaTheta plots
160 // - read ON (or OFF or MC) data
161 // - pad the events;
162 // - write root file for ON (or OFF or MC) data (ON1.root, ...);
163
164 Bool_t JobA = kTRUE;
165 Bool_t WPad = kTRUE; // write out padding histograms ?
166 Bool_t RPad = kFALSE; // read in padding histograms ?
167 Bool_t Wout = kFALSE; // write out root file ON1.root
168 // (or OFF1.root or MC1.root)?
169
170
171 // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file
172 // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
173 // - if CTrainRF = TRUE : create matrices of training events
174 // - if RTree = TRUE : read in trees, otherwise create trees
175 // - calculate hadroness for method of RANDOM FOREST
176 // - update the input files with the hadronesses (ON2.root, OFF2.root
177 // or MC2.root)
178
179 Bool_t JobB_RF_UP = kFALSE;
180 Bool_t RTrainRF = kFALSE; // read in matrices of training events
181 Bool_t CTrainRF = kFALSE; // create matrices of training events
182 Bool_t RTree = kFALSE; // read in trees
183 Bool_t WRF = kFALSE; // update input root file ?
184
185
186
187 // Job C:
188 // - read ON1.root and MC1.root files
189 // which should have been updated to contain the hadronnesses
190 // for the method of
191 // NEAREST NEIGHBORS
192 // SUPERCUTS and
193 // RF
194 // - produce Neyman-Pearson plots
195
196 Bool_t JobC = kFALSE;
197
198
199 // Job D :
200 // - select g/h separation method XX
201 // - read ON2 (or MC2) root file
202 // - apply cuts in hadronness
203 // - make plots
204
205 Bool_t JobD = kFALSE;
206
207
208
209 // Job E_EST_UP :
210 // - read MC1.root file
211 // - select g/h separation method XX
212 // - optimize energy estimation for events passing the final cuts
213 // - write parameters of energy estimator onto file
214 // - update ON1.root, OFF1.root and MC1.root files with estimated energy
215 // (ON_XX1.root, OFF_XX1.root and MC_XX1.root)
216
217 Bool_t JobE_EST_UP = kFALSE;
218 Bool_t WESTUP = kFALSE; // update root files ?
219
220
221
222 // Job F_XX :
223 // - select g/h separation method XX
224 // - read MC_XX2.root file
225   // - calculate eff. collection area
226 // - read ON_XX2.root file
227 // - apply final cuts
228 // - calculate flux
229 // - write root file for ON data after final cuts (ON3.root))
230
231 Bool_t JobF_XX = kFALSE;
232 Bool_t WXX = kFALSE; // write out root file ON3.root ?
233
234
235 //************************************************************************
236
237
238 //---------------------------------------------------------------------
239 // Job A
240 //=========
241 // read ON data file
242
243 // - produce the 2D-histogram "sigmabar versus Theta"
244 // (SigmaTheta_ON.root) for ON data
245 // (to be used for the padding of the MC gamma data)
246
247 // - write a file of ON events (ON1.root)
248 // (after the standard cuts, before the g/h separation)
249 // (to be used together with the corresponding MC gamma file (MC1.root)
250 // for the optimization of the g/h separation)
251
252
253 if (JobA)
254 {
255 gLog << "=====================================================" << endl;
256 gLog << "Macro CT1Analysis : Start of Job A" << endl;
257 gLog << "" << endl;
258 gLog << "Macro CT1Analysis : JobA, WPad, RPad, Wout = "
259 << (JobA ? "kTRUE" : "kFALSE") << ", "
260 << (WPad ? "kTRUE" : "kFALSE") << ", "
261 << (RPad ? "kTRUE" : "kFALSE") << ", "
262 << (Wout ? "kTRUE" : "kFALSE") << endl;
263
264
265 //--------------------------------------------------
266 // names of ON and OFF files to be read
267 // for generating the histograms to be used in the padding
268 TString fileON = onfile;
269 TString fileOFF = offfile;
270 gLog << "fileON, fileOFF = " << fileON << ", " << fileOFF << endl;
271
272 // name of file to conatin the histograms for the padding
273 TString outNameSigTh = outPath;
274 outNameSigTh += "SigmaTheta";
275 outNameSigTh += ".root";
276
277 //--------------------------------------------------
278 // type of data to be padded
279 TString typeInput = "ON";
280 //TString typeInput = "OFF";
281 //TString typeInput = "MC";
282 gLog << "typeInput = " << typeInput << endl;
283
284
285 // name of input root file
286 if (typeInput == "ON")
287 TString filenamein(onfile);
288 else if (typeInput == "OFF")
289 TString filenamein(offfile);
290 else if (typeInput == "MC")
291 TString filenamein(mcfile);
292 gLog << "data to be padded : " << filenamein << endl;
293
294 // name of output root file
295 TString outNameImage = outPath;
296 outNameImage += typeInput;
297 outNameImage += "1.root";
298 gLog << "padded data to be written onto : " << outNameImage << endl;
299
300 //--------------------------------------------------
301
302 //************************************************************
303 // generate histograms to be used in the padding
304 //
305 // read ON and OFF data files
306 // generate (or read in) the padding histograms for ON and OFF data
307 // and merge these histograms
308
309 MCT1PadONOFF pad;
310 pad.SetName("MCT1PadONOFF");
311 pad.SetPadFlag(1);
312 pad.SetDataType(typeInput);
313
314 // generate the padding histograms
315 if (!RPad)
316 {
317 gLog << "=====================================================" << endl;
318 gLog << "Start generating the padding histograms" << endl;
319 //-----------------------------------------
320 // ON events
321
322 gLog << "-----------" << endl;
323 gLog << "ON events :" << endl;
324 gLog << "-----------" << endl;
325
326 MTaskList tliston;
327 MParList pliston;
328
329 MCT1ReadPreProc readON(fileON);
330
331 //MFCT1SelBasic selthetaon;
332 //selthetaon.SetCuts(-100.0, 29.5, 35.5);
333 //MContinue contthetaon(&selthetaon);
334
335 MBlindPixelCalc blindon;
336 blindon.SetUseBlindPixels();
337
338 MFCT1SelBasic selbasicon;
339 MContinue contbasicon(&selbasicon);
340
341 MHBlindPixels blindON("BlindPixelsON");
342 MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels");
343 fillblindON.SetName("FillBlind");
344
345 MSigmabarCalc sigbarcalcon;
346
347 MHSigmaTheta sigthON("SigmaThetaON");
348 MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt");
349 fillsigthetaON.SetName("FillSigTheta");
350
351 //*****************************
352 // entries in MParList
353
354 pliston.AddToList(&tliston);
355 InitBinnings(&pliston);
356 pliston.AddToList(&blindON);
357 pliston.AddToList(&sigthON);
358
359
360 //*****************************
361 // entries in MTaskList
362
363 tliston.AddToList(&readON);
364 //tliston.AddToList(&contthetaon);
365
366 tliston.AddToList(&blindon);
367
368 tliston.AddToList(&contbasicon);
369 tliston.AddToList(&fillblindON);
370 tliston.AddToList(&sigbarcalcon);
371 tliston.AddToList(&fillsigthetaON);
372
373 MProgressBar baron;
374 MEvtLoop evtloopon;
375 evtloopon.SetParList(&pliston);
376 evtloopon.SetProgressBar(&baron);
377
378 Int_t maxeventson = -1;
379 //Int_t maxeventson = 10000;
380 if ( !evtloopon.Eventloop(maxeventson) )
381 return;
382
383 tliston.PrintStatistics(0, kTRUE);
384
385 blindON.DrawClone();
386 sigthON.DrawClone();
387
388 // save the histograms for the padding
389 TH2D *sigthon = sigthON.GetSigmaTheta();
390 TH3D *sigpixthon = sigthON.GetSigmaPixTheta();
391 TH3D *diffpixthon = sigthON.GetDiffPixTheta();
392
393 TH2D *blindidthon = blindON.GetBlindId();
394 TH2D *blindnthon = blindON.GetBlindN();
395
396 //-----------------------------------------
397 // OFF events
398
399 gLog << "------------" << endl;
400 gLog << "OFF events :" << endl;
401 gLog << "------------" << endl;
402
403 MTaskList tlistoff;
404 MParList plistoff;
405
406 MCT1ReadPreProc readOFF(fileOFF);
407
408 MFCT1SelBasic selthetaoff;
409 selthetaoff.SetCuts(-100.0, 29.5, 35.5);
410 MContinue contthetaoff(&selthetaoff);
411
412 MBlindPixelCalc blindoff;
413 blindoff.SetUseBlindPixels();
414
415 MFCT1SelBasic selbasicoff;
416 MContinue contbasicoff(&selbasicoff);
417
418 MHBlindPixels blindOFF("BlindPixelsOFF");
419 MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels");
420 fillblindOFF.SetName("FillBlindOFF");
421
422 MSigmabarCalc sigbarcalcoff;
423
424 MHSigmaTheta sigthOFF("SigmaThetaOFF");
425 MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt");
426 fillsigthetaOFF.SetName("FillSigThetaOFF");
427
428 //*****************************
429 // entries in MParList
430
431 plistoff.AddToList(&tlistoff);
432 InitBinnings(&plistoff);
433 plistoff.AddToList(&blindOFF);
434 plistoff.AddToList(&sigthOFF);
435
436
437 //*****************************
438 // entries in MTaskList
439
440 tlistoff.AddToList(&readOFF);
441 //tlistoff.AddToList(&contthetaoff);
442
443 tlistoff.AddToList(&blindoff);
444
445 tlistoff.AddToList(&contbasicoff);
446 tlistoff.AddToList(&fillblindOFF);
447 tlistoff.AddToList(&sigbarcalcoff);
448 tlistoff.AddToList(&fillsigthetaOFF);
449
450 MProgressBar baroff;
451 MEvtLoop evtloopoff;
452 evtloopoff.SetParList(&plistoff);
453 evtloopoff.SetProgressBar(&baroff);
454
455 Int_t maxeventsoff = -1;
456 //Int_t maxeventsoff = 20000;
457 if ( !evtloopoff.Eventloop(maxeventsoff) )
458 return;
459
460 tlistoff.PrintStatistics(0, kTRUE);
461
462 blindOFF.DrawClone();
463 sigthOFF.DrawClone();
464
465 // save the histograms for the padding
466 TH2D *sigthoff = sigthOFF.GetSigmaTheta();
467 TH3D *sigpixthoff = sigthOFF.GetSigmaPixTheta();
468 TH3D *diffpixthoff = sigthOFF.GetDiffPixTheta();
469
470 TH2D *blindidthoff = blindOFF.GetBlindId();
471 TH2D *blindnthoff = blindOFF.GetBlindN();
472
473
474 //-----------------------------------------
475
476 gLog << "End of generating the padding histograms" << endl;
477 gLog << "=====================================================" << endl;
478
479 pad.MergeHistograms(sigthon, sigthoff,
480 sigpixthon, sigpixthoff,
481 diffpixthon, diffpixthoff,
482 blindidthon, blindidthoff,
483 blindnthon, blindnthoff);
484
485 if (WPad)
486 {
487 // write the padding histograms onto a file ---------
488 pad.WriteTargetDist(outNameSigTh);
489 }
490 }
491
492 // read the padding histograms ---------------------------
493 if (RPad)
494 {
495 pad.ReadTargetDist(outNameSigTh);
496 }
497
498
499 //************************************************************
500
501 if (Wout)
502 {
503 gLog << "=====================================================" << endl;
504 gLog << "Start the padding" << endl;
505
506 //-----------------------------------------------------------
507 MTaskList tliston;
508 MParList pliston;
509
510 char *sourceName = "MSrcPosCam";
511 MSrcPosCam source(sourceName);
512
513 // geometry is needed in MHHillas... classes
514 MGeomCam *fGeom =
515 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
516
517 //-------------------------------------------
518 // create the tasks which should be executed
519 //
520
521 MCT1ReadPreProc read(filenamein);
522
523 MFCT1SelBasic seltheta;
524 seltheta.SetCuts(-100.0, 29.5, 35.5);
525 MContinue conttheta(&seltheta);
526
527 if (typeInput == "ON")
528 {
529 MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc",
530 "MCT1PointingCorrCalc");
531 }
532
533 MBlindPixelCalc blindbeforepad;
534 blindbeforepad.SetUseBlindPixels();
535 blindbeforepad.SetName("BlindBeforePadding");
536
537 MBlindPixelCalc blind;
538 blind.SetUseBlindPixels();
539 blind.SetName("BlindAfterPadding");
540
541 MFCT1SelBasic selbasic;
542 MContinue contbasic(&selbasic);
543 contbasic.SetName("SelBasic");
544
545 MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels");
546 fillblind.SetName("HBlind");
547
548 MSigmabarCalc sigbarcalc;
549
550 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
551 fillsigtheta.SetName("HSigmaTheta");
552
553 MImgCleanStd clean;
554
555
556 // calculation of image parameters ---------------------
557 TString fHilName = "MHillas";
558 TString fHilNameExt = "MHillasExt";
559 TString fHilNameSrc = "MHillasSrc";
560 TString fImgParName = "MNewImagePar";
561
562 MHillasCalc hcalc;
563 hcalc.SetNameHillas(fHilName);
564 hcalc.SetNameHillasExt(fHilNameExt);
565 hcalc.SetNameNewImgPar(fImgParName);
566
567 MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc);
568 hsrccalc.SetInput(fHilName);
569
570 MFillH hfill1("MHHillas", fHilName);
571 hfill1.SetName("HHillas");
572
573 MFillH hfill2("MHStarMap", fHilName);
574 hfill2.SetName("HStarMap");
575
576 MFillH hfill3("MHHillasExt", fHilNameSrc);
577 hfill3.SetName("HHillasExt");
578
579 MFillH hfill4("MHHillasSrc", fHilNameSrc);
580 hfill4.SetName("HHillasSrc");
581
582 MFillH hfill5("MHNewImagePar", fImgParName);
583 hfill5.SetName("HNewImagePar");
584 // --------------------------------------------------
585
586 MFCT1SelStandard selstandard(fHilNameSrc);
587 selstandard.SetHillasName(fHilName);
588 selstandard.SetImgParName(fImgParName);
589 selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0);
590 MContinue contstandard(&selstandard);
591 contstandard.SetName("SelStandard");
592
593
594 MWriteRootFile write(outNameImage);
595
596 write.AddContainer("MRawRunHeader", "RunHeaders");
597 write.AddContainer("MTime", "Events");
598 write.AddContainer("MMcEvt", "Events");
599 write.AddContainer("ThetaOrig", "Events");
600 write.AddContainer("MSrcPosCam", "Events");
601 write.AddContainer("MSigmabar", "Events");
602 write.AddContainer("MHillas", "Events");
603 write.AddContainer("MHillasExt", "Events");
604 write.AddContainer("MHillasSrc", "Events");
605 write.AddContainer("MNewImagePar", "Events");
606
607
608 //*****************************
609 // entries in MParList
610
611 pliston.AddToList(&tliston);
612 InitBinnings(&pliston);
613
614 pliston.AddToList(&source);
615
616
617 //*****************************
618 // entries in MTaskList
619
620 tliston.AddToList(&read);
621 //tliston.AddToList(&conttheta);
622
623 tliston.AddToList(&blindbeforepad);
624 tliston.AddToList(&pad);
625 if (typeInput == "ON")
626 tliston.AddToList(&pointcorr);
627 tliston.AddToList(&blind);
628
629 tliston.AddToList(&contbasic);
630 tliston.AddToList(&fillblind);
631 tliston.AddToList(&sigbarcalc);
632 tliston.AddToList(&fillsigtheta);
633 tliston.AddToList(&clean);
634
635 tliston.AddToList(&hcalc);
636 tliston.AddToList(&hsrccalc);
637
638 tliston.AddToList(&hfill1);
639 tliston.AddToList(&hfill2);
640 tliston.AddToList(&hfill3);
641 tliston.AddToList(&hfill4);
642 tliston.AddToList(&hfill5);
643
644 tliston.AddToList(&contstandard);
645 tliston.AddToList(&write);
646
647 //*****************************
648
649 //-------------------------------------------
650 // Execute event loop
651 //
652 MProgressBar bar;
653 MEvtLoop evtloop;
654 evtloop.SetParList(&pliston);
655 //evtloop.ReadEnv(env, "", printEnv);
656 evtloop.SetProgressBar(&bar);
657 // evtloop.Write();
658
659 Int_t maxevents = -1;
660 //Int_t maxevents = 1000;
661 if ( !evtloop.Eventloop(maxevents) )
662 return;
663
664 tliston.PrintStatistics(0, kTRUE);
665
666
667 //-------------------------------------------
668 // Display the histograms
669
670 pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone();
671 pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone();
672
673 pliston.FindObject("MHHillas")->DrawClone();
674 pliston.FindObject("MHHillasExt")->DrawClone();
675 pliston.FindObject("MHHillasSrc")->DrawClone();
676 pliston.FindObject("MHNewImagePar")->DrawClone();
677 pliston.FindObject("MHStarMap")->DrawClone();
678
679 DeleteBinnings(&pliston);
680
681 gLog << "End of padding" << endl;
682 gLog << "=====================================================" << endl;
683 }
684
685
686 gLog << "Macro CT1Analysis : End of Job A" << endl;
687 gLog << "===================================================" << endl;
688 }
689
690
691
692
693
694 //---------------------------------------------------------------------
695 // Job B_RF_UP
696 //============
697
698
699 // - create (or read in) the matrices of training events for gammas
700 // and hadrons
701 // - create (or read in) the trees
702 // - then read ON1.root (or MC1.root) file
703 // - calculate the hadroness for the method of RANDOM FOREST
704 // - update input root file with the hadroness
705
706
707 if (JobB_RF_UP)
708 {
709 gLog << "=====================================================" << endl;
710 gLog << "Macro CT1Analysis : Start of Job B_RF_UP" << endl;
711
712 gLog << "" << endl;
713 gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
714 << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", "
715 << (RTrainRF ? "kTRUE" : "kFALSE") << ", "
716 << (CTrainRF ? "kTRUE" : "kFALSE") << ", "
717 << (RTree ? "kTRUE" : "kFALSE") << ", "
718 << (WRF ? "kTRUE" : "kFALSE") << endl;
719
720
721 //--------------------------------------------
722 // parameters for the random forest
723 Int_t NumTrees = 100;
724 Int_t NumTry = 3;
725 Int_t NdSize = 1;
726
727
728 TString hadRFName = "HadRF";
729 Float_t maxhadronness = 0.23;
730 Float_t maxalpha = 20.0;
731 Float_t maxdist = 10.0;
732
733 TString fHilName = "MHillas";
734 TString fHilNameExt = "MHillasExt";
735 TString fHilNameSrc = "MHillasSrc";
736 TString fImgParName = "MNewImagePar";
737
738
739 //--------------------------------------------
740 // file to be updated (ON, OFF or MC)
741
742 TString typeInput = "ON";
743 //TString typeInput = "OFF";
744 //TString typeInput = "MC";
745 gLog << "typeInput = " << typeInput << endl;
746
747 // name of input root file
748 TString filenameData = outPath;
749 filenameData += typeInput;
750 filenameData += "1.root";
751 gLog << "filenameData = " << filenameData << endl;
752
753 // name of output root file
754 TString outNameImage = outPath;
755 outNameImage += typeInput;
756 outNameImage += "2.root";
757 //TString outNameImage = filenameData;
758
759 gLog << "outNameImage = " << outNameImage << endl;
760
761 //--------------------------------------------
762 // files to be read for generating the matrices of training events
763 // "hadrons" :
764 TString filenameHad = outPath;
765 filenameHad += "OFF";
766 filenameHad += "1.root";
767 Int_t howManyHadrons = 1000000;
768 gLog << "filenameHad = " << filenameHad << ", howManyHadrons = "
769 << howManyHadrons << endl;
770
771
772 // "gammas" :
773 TString filenameMC = outPath;
774 filenameMC += "MC";
775 filenameMC += "1.root";
776 Int_t howManyGammas = 50000;
777 gLog << "filenameMC = " << filenameMC << ", howManyGammas = "
778 << howManyGammas << endl;
779
780 //--------------------------------------------
781 // files of training events
782
783 TString outNameGammas = outPath;
784 outNameGammas += "RFmatrix_gammas_";
785 outNameGammas += "MC";
786 outNameGammas += ".root";
787
788 TString typeMatrixHadrons = "OFF";
789 gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
790
791 TString outNameHadrons = outPath;
792 outNameHadrons += "RFmatrix_hadrons_";
793 outNameHadrons += typeMatrixHadrons;
794 outNameHadrons += ".root";
795
796
797 MHMatrix matrixg("MatrixGammas");
798 matrixg.EnableGraphicalOutput();
799
800 matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
801 matrixg.AddColumn("MSigmabar.fSigmabar");
802 matrixg.AddColumn("log10(MHillas.fSize)");
803 matrixg.AddColumn("MHillasSrc.fDist");
804 matrixg.AddColumn("MHillas.fWidth");
805 matrixg.AddColumn("MHillas.fLength");
806 matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
807 matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
808 matrixg.AddColumn("MNewImagePar.fConc");
809 matrixg.AddColumn("MNewImagePar.fLeakage1");
810
811 MHMatrix matrixh("MatrixHadrons");
812 matrixh.EnableGraphicalOutput();
813
814 matrixh.AddColumns(matrixg.GetColumns());
815
816 //--------------------------------------------
817 // file of trees of the random forest
818
819 TString outRF = outPath;
820 outRF += "RF.root";
821
822
823 //*************************************************************************
824 // read in matrices of training events
825if (RTrainRF)
826 {
827 const char* mtxName = "MatrixGammas";
828
829 gLog << "" << endl;
830 gLog << "========================================================" << endl;
831 gLog << "Get matrix for (gammas)" << endl;
832 gLog << "matrix name = " << mtxName << endl;
833 gLog << "name of root file = " << outNameGammas << endl;
834 gLog << "" << endl;
835
836
837 // read in the object with the name 'mtxName' from file 'outNameGammas'
838 //
839 TFile fileg(outNameGammas);
840
841 matrixg.Read(mtxName);
842 matrixg.Print("SizeCols");
843
844
845 //*****************************************************************
846
847 const char* mtxName = "MatrixHadrons";
848
849 gLog << "" << endl;
850 gLog << "========================================================" << endl;
851 gLog << " Get matrix for (hadrons)" << endl;
852 gLog << "matrix name = " << mtxName << endl;
853 gLog << "name of root file = " << outNameHadrons << endl;
854 gLog << "" << endl;
855
856
857 // read in the object with the name 'mtxName' from file 'outNameHadrons'
858 //
859 TFile fileh(outNameHadrons);
860
861 matrixh.Read(mtxName);
862 matrixh.Print("SizeCols");
863 }
864
865
866 //*************************************************************************
867 // create matrices of training events
868if (CTrainRF)
869 {
870 gLog << "" << endl;
871 gLog << "========================================================" << endl;
872 gLog << " Create matrices of training events" << endl;
873 gLog << " Gammas :" << endl;
874
875
876 MParList plistg;
877 MTaskList tlistg;
878 MFilterList flistg;
879
880 MParList plisth;
881 MTaskList tlisth;
882 MFilterList flisth;
883
884 MReadMarsFile readg("Events", filenameMC);
885 readg.DisableAutoScheme();
886
887 MReadMarsFile readh("Events", filenameHad);
888 readh.DisableAutoScheme();
889
890 MFParticleId fgamma("MMcEvt", '=', kGAMMA);
891 fgamma.SetName("gammaID");
892 MFParticleId fhadrons("MMcEvt", '!', kGAMMA);
893 fhadrons.SetName("hadronID)");
894
895 TString mgname("costhg");
896 MBinning bing("Binning"+mgname);
897 bing.SetEdges(10, 0., 1.0);
898
899 MH3 gref("cos(MMcEvt.fTelescopeTheta)");
900 gref.SetName(mgname);
901 MH::SetBinning(&gref.GetHist(), &bing);
902 for (Int_t i=1; i<=gref.GetNbins(); i++)
903 gref.GetHist().SetBinContent(i, 1.0);
904
905 MFEventSelector2 selectorg(gref);
906 selectorg.SetNumMax(howManyGammas);
907 selectorg.SetName("selectGammas");
908
909 MFillH fillmatg("MatrixGammas");
910 fillmatg.SetFilter(&flistg);
911 fillmatg.SetName("fillGammas");
912
913
914 //***************************** fill gammas ***
915 // entries in MFilterList
916
917 flistg.AddToList(&fgamma);
918 flistg.AddToList(&selectorg);
919
920 //*****************************
921 // entries in MParList
922
923 plistg.AddToList(&tlistg);
924 InitBinnings(&plistg);
925
926 plistg.AddToList(&matrixg);
927
928 //*****************************
929 // entries in MTaskList
930
931 tlistg.AddToList(&readg);
932 tlistg.AddToList(&flistg);
933 tlistg.AddToList(&fillmatg);
934
935 //*****************************
936
937 MProgressBar matrixbar;
938 MEvtLoop evtloopg;
939 evtloopg.SetParList(&plistg);
940 //evtloopg.ReadEnv(env, "", printEnv);
941 evtloopg.SetProgressBar(&matrixbar);
942
943 Int_t maxevents = -1;
944 if (!evtloopg.Eventloop(maxevents))
945 return;
946
947 tlistg.PrintStatistics(0, kTRUE);
948
949
950 //***************************** fill hadrons ***
951
952 gLog << " Hadrons :" << endl;
953
954 TString mhname("costhh");
955 MBinning binh("Binning"+mhname);
956 binh.SetEdges(10, 0., 1.0);
957
958 MH3 href("cos(MMcEvt.fTelescopeTheta)");
959 href.SetName(mhname);
960 MH::SetBinning(&href.GetHist(), &binh);
961 for (Int_t i=1; i<=href.GetNbins(); i++)
962 href.GetHist().SetBinContent(i, 1.0);
963
964 MFEventSelector2 selectorh(href);
965 //selectorh.SetNumMax(howManyHadrons);
966 // select as many hadrons as gammas
967 selectorh.SetNumMax(matrixg.GetM().GetNrows());
968 selectorh.SetName("selectHadrons");
969
970 MFillH fillmath("MatrixHadrons");
971 fillmath.SetFilter(&flisth);
972 fillmath.SetName("fillHadrons");
973
974
975 // entries in MFilterList
976
977 flisth.AddToList(&fhadrons);
978 flisth.AddToList(&selectorh);
979
980 //*****************************
981 // entries in MParList
982
983 plisth.AddToList(&tlisth);
984 InitBinnings(&plisth);
985
986 plisth.AddToList(&matrixh);
987
988 //*****************************
989 // entries in MTaskList
990
991 tlisth.AddToList(&readh);
992 tlisth.AddToList(&flisth);
993 tlisth.AddToList(&fillmath);
994
995 //*****************************
996
997 MProgressBar matrixbar;
998 MEvtLoop evtlooph;
999 evtlooph.SetParList(&plisth);
1000 //evtlooph.ReadEnv(env, "", printEnv);
1001 evtlooph.SetProgressBar(&matrixbar);
1002
1003 Int_t maxevents = -1;
1004 if (!evtlooph.Eventloop(maxevents))
1005 return;
1006
1007 tlisth.PrintStatistics(0, kTRUE);
1008 //*****************************************************
1009
1010
1011 // write out matrices of training events events
1012
1013 gLog << "" << endl;
1014 gLog << "========================================================" << endl;
1015 gLog << "Write out matrices of training events" << endl;
1016
1017
1018 //-------------------------------------------
1019 // "gammas"
1020 gLog << "Gammas :" << endl;
1021 matrixg.Print("SizeCols");
1022
1023 TFile writeg(outNameGammas, "RECREATE", "");
1024 matrixg.Write();
1025
1026 gLog << "" << endl;
1027 gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file "
1028 << outNameGammas << endl;
1029
1030 //-------------------------------------------
1031 // "hadrons"
1032 gLog << "Hadrons :" << endl;
1033 matrixh.Print("SizeCols");
1034
1035 TFile writeh(outNameHadrons, "RECREATE", "");
1036 matrixh.Write();
1037
1038 gLog << "" << endl;
1039 gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file "
1040 << outNameHadrons << endl;
1041
1042 }
1043 //********** end of creating matrices of training events ***********
1044
1045
1046 MRanForest *fRanForest;
1047 MRanTree *fRanTree;
1048 //-----------------------------------------------------------------
1049 // read in the trees of the random forest
1050 if (RTree)
1051 {
1052 MParList plisttr;
1053 MTaskList tlisttr;
1054 plisttr.AddToList(&tlisttr);
1055
1056 MReadTree readtr("TREE", outRF);
1057 readtr.DisableAutoScheme();
1058
1059 MRanForestFill rffill;
1060 rffill.SetNumTrees(NumTrees);
1061
1062 // list of tasks for the loop over the trees
1063
1064 tlisttr.AddToList(&readtr);
1065 tlisttr.AddToList(&rffill);
1066
1067 //-------------------
1068 // Execute tree loop
1069 //
1070 MEvtLoop evtlooptr;
1071 evtlooptr.SetName("ReadRFTrees");
1072 evtlooptr.SetParList(&plisttr);
1073 if (!evtlooptr.Eventloop())
1074 return;
1075
1076 tlisttr.PrintStatistics(0, kTRUE);
1077
1078
1079 // get adresses of objects which are used in the next eventloop
1080 fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
1081 if (!fRanForest)
1082 {
1083 gLog << err << dbginf << "MRanForest not found... aborting." << endl;
1084 return kFALSE;
1085 }
1086
1087 fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
1088 if (!fRanTree)
1089 {
1090 gLog << err << dbginf << "MRanTree not found... aborting." << endl;
1091 return kFALSE;
1092 }
1093
1094 }
1095
1096 //-----------------------------------------------------------------
1097 // grow the trees of the random forest (event loop = tree loop)
1098
1099 if (!RTree)
1100 {
1101
1102 gLog << "" << endl;
1103 gLog << "========================================================" << endl;
1104 gLog << "Macro CT1Analysis : start growing trees" << endl;
1105
1106 MTaskList tlist2;
1107 MParList plist2;
1108 plist2.AddToList(&tlist2);
1109
1110 plist2.AddToList(&matrixg);
1111 plist2.AddToList(&matrixh);
1112
1113 MRanForestGrow rfgrow2;
1114 rfgrow2.SetNumTrees(NumTrees);
1115 rfgrow2.SetNumTry(NumTry);
1116 rfgrow2.SetNdSize(NdSize);
1117
1118 MWriteRootFile rfwrite2(outRF);
1119 rfwrite2.AddContainer("MRanTree", "TREE");
1120
1121 // list of tasks for the loop over the trees
1122
1123 tlist2.AddToList(&rfgrow2);
1124 tlist2.AddToList(&rfwrite2);
1125
1126 //-------------------
1127 // Execute tree loop
1128 //
1129 MEvtLoop treeloop;
1130 treeloop.SetName("GrowRFTrees");
1131 treeloop.SetParList(&plist2);
1132
1133 if ( !treeloop.Eventloop() )
1134 return;
1135
1136 tlist2.PrintStatistics(0, kTRUE);
1137
1138 // get adresses of objects which are used in the next eventloop
1139 fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
1140 if (!fRanForest)
1141 {
1142 gLog << err << dbginf << "MRanForest not found... aborting." << endl;
1143 return kFALSE;
1144 }
1145
1146 fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
1147 if (!fRanTree)
1148 {
1149 gLog << err << dbginf << "MRanTree not found... aborting." << endl;
1150 return kFALSE;
1151 }
1152
1153 }
1154 // end of growing the trees of the random forest
1155 //-----------------------------------------------------------------
1156
1157
1158
1159
1160 //-----------------------------------------------------------------
1161 // Update the input files with the RF hadronness
1162 //
1163
1164 if (WRF)
1165 {
1166 gLog << "" << endl;
1167 gLog << "========================================================" << endl;
1168 gLog << "Update input file '" << filenameData
1169 << "' with the RF hadronness" << endl;
1170
1171 MTaskList tliston;
1172 MParList pliston;
1173
1174
1175 // geometry is needed in MHHillas... classes
1176 MGeomCam *fGeom =
1177 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1178
1179 //-------------------------------------------
1180 // create the tasks which should be executed
1181 //
1182
1183 MReadMarsFile read("Events", filenameData);
1184 read.DisableAutoScheme();
1185
1186
1187 //.......................................................................
1188 // calculate hadronnes for method of RANDOM FOREST
1189
1190
1191 MRanForestCalc rfcalc;
1192 rfcalc.SetHadronnessName(hadRFName);
1193
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(hadRFName, "Events");
1212
1213 //-----------------------------------------------------------------
1214 // geometry is needed in MHHillas... classes
1215 MGeomCam *fGeom =
1216 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1217
1218
1219
1220 MFCT1SelFinal selfinalgh(fHilNameSrc);
1221 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1222 selfinalgh.SetHadronnessName(hadRFName);
1223 selfinalgh.SetName("SelFinalgh");
1224 MContinue contfinalgh(&selfinalgh);
1225 contfinalgh.SetName("ContSelFinalgh");
1226
1227 MFillH fillranfor("MHRanForest");
1228 fillranfor.SetName("HRanForest");
1229
1230 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
1231 fillhadrf.SetName("HhadRF");
1232
1233 MFCT1SelFinal selfinal(fHilNameSrc);
1234 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1235 selfinal.SetHadronnessName(hadRFName);
1236 selfinal.SetName("SelFinal");
1237 MContinue contfinal(&selfinal);
1238 contfinal.SetName("ContSelFinal");
1239
1240 TString mh3name = "abs(Alpha)";
1241 MBinning binsalphaabs("Binning"+mh3name);
1242 binsalphaabs.SetEdges(50, -2.0, 98.0);
1243
1244 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
1245 alphaabs.SetName(mh3name);
1246 MFillH alpha(&alphaabs);
1247 alpha.SetName("FillAlphaAbs");
1248
1249
1250 MFillH hfill1("MHHillas", fHilName);
1251 hfill1.SetName("HHillas");
1252
1253 MFillH hfill2("MHStarMap", fHilName);
1254 hfill2.SetName("HStarMap");
1255
1256 MFillH hfill3("MHHillasExt", fHilNameSrc);
1257 hfill3.SetName("HHillasExt");
1258
1259 MFillH hfill4("MHHillasSrc", fHilNameSrc);
1260 hfill4.SetName("HHillasSrc");
1261
1262 MFillH hfill5("MHNewImagePar", fImgParName);
1263 hfill5.SetName("HNewImagePar");
1264
1265 //*****************************
1266 // entries in MParList
1267
1268 pliston.AddToList(&tliston);
1269 InitBinnings(&pliston);
1270
1271 pliston.AddToList(fRanForest);
1272 pliston.AddToList(fRanTree);
1273
1274 pliston.AddToList(&binsalphaabs);
1275 pliston.AddToList(&alphaabs);
1276
1277
1278 //*****************************
1279 // entries in MTaskList
1280
1281 tliston.AddToList(&read);
1282
1283 tliston.AddToList(&rfcalc);
1284 tliston.AddToList(&fillranfor);
1285 tliston.AddToList(&fillhadrf);
1286
1287 tliston.AddToList(&write);
1288 tliston.AddToList(&contfinalgh);
1289
1290 tliston.AddToList(&alpha);
1291 tliston.AddToList(&hfill1);
1292 tliston.AddToList(&hfill2);
1293 tliston.AddToList(&hfill3);
1294 tliston.AddToList(&hfill4);
1295 tliston.AddToList(&hfill5);
1296
1297 tliston.AddToList(&contfinal);
1298
1299 //*****************************
1300
1301 //-------------------------------------------
1302 // Execute event loop
1303 //
1304 MProgressBar bar;
1305 MEvtLoop evtloop;
1306 evtloop.SetParList(&pliston);
1307 evtloop.SetProgressBar(&bar);
1308
1309 Int_t maxevents = -1;
1310 if ( !evtloop.Eventloop(maxevents) )
1311 return;
1312
1313 tliston.PrintStatistics(0, kTRUE);
1314
1315
1316 //-------------------------------------------
1317 // Display the histograms
1318 //
1319 pliston.FindObject("MHRanForest")->DrawClone();
1320 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
1321
1322 pliston.FindObject("MHHillas")->DrawClone();
1323 pliston.FindObject("MHHillasExt")->DrawClone();
1324 pliston.FindObject("MHHillasSrc")->DrawClone();
1325 pliston.FindObject("MHNewImagePar")->DrawClone();
1326 pliston.FindObject("MHStarMap")->DrawClone();
1327
1328
1329 //-------------------------------------------
1330 // fit alpha distribution to get the number of excess events and
1331 // calculate significance of gamma signal in the alpha plot
1332
1333 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
1334 TH1 &alphaHist = absalpha->GetHist();
1335 alphaHist.SetXTitle("|alpha| [\\circ]");
1336 alphaHist.SetName("alpha-macro");
1337
1338 Double_t alphasig = 13.1;
1339 Double_t alphamin = 30.0;
1340 Double_t alphamax = 90.0;
1341 Int_t degree = 2;
1342 Double_t significance = -99.0;
1343 Bool_t drawpoly = kTRUE;
1344 Bool_t fitgauss = kTRUE;
1345 Bool_t print = kTRUE;
1346
1347 MHFindSignificance findsig;
1348 findsig.SetRebin(kTRUE);
1349 findsig.SetReduceDegree(kFALSE);
1350
1351 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
1352 alphasig, drawpoly, fitgauss, print);
1353 significance = findsig.GetSignificance();
1354 Float_t alphasi = findsig.GetAlphasi();
1355
1356 gLog << "For file '" << filenameData << "' : " << endl;
1357 gLog << "Significance of gamma signal after supercuts : "
1358 << significance << " (for |alpha| < " << alphasi << " degrees)"
1359 << endl;
1360
1361 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
1362
1363 //-------------------------------------------
1364
1365
1366 DeleteBinnings(&pliston);
1367 }
1368
1369 gLog << "Macro CT1Analysis : End of Job B_RF_UP" << endl;
1370 gLog << "=======================================================" << endl;
1371 }
1372 //---------------------------------------------------------------------
1373
1374
1375
1376 //---------------------------------------------------------------------
1377 // Job C
1378 //======
1379
1380 // - read ON1 and MC1 data files
1381 // which should have been updated to contain the hadronnesses
1382 // for the method of NEAREST NEIGHBORS and for the SUOERCUTS
1383 // - produce Neyman-Pearson plots
1384
1385 if (JobC)
1386 {
1387 gLog << "=====================================================" << endl;
1388 gLog << "Macro CT1Analysis : Start of Job C" << endl;
1389
1390 gLog << "" << endl;
1391 gLog << "Macro CT1Analysis : JobC = " << JobC << endl;
1392
1393
1394 // name of input data file
1395 TString filenameData = outPath;
1396 filenameData += "OFF";
1397 filenameData += "2.root";
1398 gLog << "filenameData = " << filenameData << endl;
1399
1400 // name of input MC file
1401 TString filenameMC = outPath;
1402 filenameMC += "MC";
1403 filenameMC += "2.root";
1404 gLog << "filenameMC = " << filenameMC << endl;
1405
1406
1407 //-----------------------------------------------------------------
1408
1409 MTaskList tliston;
1410 MParList pliston;
1411
1412
1413 // geometry is needed in MHHillas... classes
1414 MGeomCam *fGeom =
1415 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1416
1417 //-------------------------------------------
1418 // create the tasks which should be executed
1419 //
1420
1421 MReadMarsFile read("Events", filenameMC);
1422 read.AddFile(filenameData);
1423 read.DisableAutoScheme();
1424
1425
1426 //.......................................................................
1427 // names of hadronness containers
1428
1429 TString hadNNName = "HadNN";
1430 TString hadSCName = "HadSC";
1431 TString hadRFName = "HadRF";
1432
1433 //.......................................................................
1434
1435
1436 TString fHilName = "MHillas";
1437 TString fHilNameExt = "MHillasExt";
1438 TString fHilNameSrc = "MHillasSrc";
1439 TString fImgParName = "MNewImagePar";
1440
1441 Float_t maxhadronness = 0.40;
1442 Float_t maxalpha = 20.0;
1443 Float_t maxdist = 10.0;
1444
1445 MFCT1SelFinal selfinalgh(fHilNameSrc);
1446 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1447 selfinalgh.SetHadronnessName(hadRFName);
1448 selfinalgh.SetName("SelFinalgh");
1449 MContinue contfinalgh(&selfinalgh);
1450 contfinalgh.SetName("ContSelFinalgh");
1451
1452 //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
1453 //fillhadnn.SetName("HhadNN");
1454 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
1455 fillhadsc.SetName("HhadSC");
1456 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
1457 fillhadrf.SetName("HhadRF");
1458
1459 MFCT1SelFinal selfinal(fHilNameSrc);
1460 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1461 selfinal.SetHadronnessName(hadRFName);
1462 selfinal.SetName("SelFinal");
1463 MContinue contfinal(&selfinal);
1464 contfinal.SetName("ContSelFinal");
1465
1466
1467 MFillH hfill1("MHHillas", fHilName);
1468 hfill1.SetName("HHillas");
1469
1470 MFillH hfill2("MHStarMap", fHilName);
1471 hfill2.SetName("HStarMap");
1472
1473 MFillH hfill3("MHHillasExt", fHilNameSrc);
1474 hfill3.SetName("HHillasExt");
1475
1476 MFillH hfill4("MHHillasSrc", fHilNameSrc);
1477 hfill4.SetName("HHillasSrc");
1478
1479 MFillH hfill5("MHNewImagePar", fImgParName);
1480 hfill5.SetName("HNewImagePar");
1481
1482
1483 //*****************************
1484 // entries in MParList
1485
1486 pliston.AddToList(&tliston);
1487 InitBinnings(&pliston);
1488
1489
1490 //*****************************
1491 // entries in MTaskList
1492
1493 tliston.AddToList(&read);
1494
1495 //tliston.AddToList(&fillhadnn);
1496 tliston.AddToList(&fillhadsc);
1497 tliston.AddToList(&fillhadrf);
1498
1499 tliston.AddToList(&contfinalgh);
1500 tliston.AddToList(&hfill1);
1501 tliston.AddToList(&hfill2);
1502 tliston.AddToList(&hfill3);
1503 tliston.AddToList(&hfill4);
1504 tliston.AddToList(&hfill5);
1505
1506 tliston.AddToList(&contfinal);
1507
1508 //*****************************
1509
1510 //-------------------------------------------
1511 // Execute event loop
1512 //
1513 MProgressBar bar;
1514 MEvtLoop evtloop;
1515 evtloop.SetParList(&pliston);
1516 evtloop.SetProgressBar(&bar);
1517
1518 Int_t maxevents = -1;
1519 //Int_t maxevents = 35000;
1520 if ( !evtloop.Eventloop(maxevents) )
1521 return;
1522
1523 tliston.PrintStatistics(0, kTRUE);
1524
1525
1526 //-------------------------------------------
1527 // Display the histograms
1528 //
1529
1530 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
1531 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
1532 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
1533
1534 pliston.FindObject("MHHillas")->DrawClone();
1535 pliston.FindObject("MHHillasExt")->DrawClone();
1536 pliston.FindObject("MHHillasSrc")->DrawClone();
1537 pliston.FindObject("MHNewImagePar")->DrawClone();
1538 pliston.FindObject("MHStarMap")->DrawClone();
1539
1540 DeleteBinnings(&pliston);
1541
1542 gLog << "Macro CT1Analysis : End of Job C" << endl;
1543 gLog << "===================================================" << endl;
1544 }
1545
1546
1547 //---------------------------------------------------------------------
1548 // Job D
1549 //======
1550
1551 // - select g/h separation method XX
1552 // - read ON2 (or MC2) root file
1553 // - apply cuts in hadronness
1554 // - make plots
1555
1556
1557 if (JobD)
1558 {
1559 gLog << "=====================================================" << endl;
1560 gLog << "Macro CT1Analysis : Start of Job D" << endl;
1561
1562 gLog << "" << endl;
1563 gLog << "Macro CT1Analysis : JobD = "
1564 << JobD << endl;
1565
1566 // type of data to be analysed
1567 TString typeData = "ON";
1568 //TString typeData = "OFF";
1569 //TString typeData = "MC";
1570 gLog << "typeData = " << typeData << endl;
1571
1572 TString ext = "2.root";
1573
1574
1575 //------------------------------
1576 // selection of g/h separation method
1577 // and definition of final selections
1578
1579 //TString XX("NN");
1580 //TString XX("SC");
1581 TString XX("RF");
1582 TString fhadronnessName("Had");
1583 fhadronnessName += XX;
1584 gLog << "fhadronnessName = " << fhadronnessName << endl;
1585
1586 // maximum values of the hadronness, |ALPHA| and DIST
1587 Float_t maxhadronness = 0.30;
1588 Float_t maxalpha = 20.0;
1589 Float_t maxdist = 10.0;
1590 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
1591 << maxhadronness << ", " << maxalpha << ", "
1592 << maxdist << endl;
1593
1594
1595 //------------------------------
1596 // name of data file to be analysed
1597 TString filenameData(outPath);
1598 filenameData += typeData;
1599 filenameData += ext;
1600 gLog << "filenameData = " << filenameData << endl;
1601
1602
1603
1604 //*************************************************************************
1605 //
1606 // Analyse the data
1607 //
1608
1609 MTaskList tliston;
1610 MParList pliston;
1611
1612 // geometry is needed in MHHillas... classes
1613 MGeomCam *fGeom =
1614 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1615
1616
1617 TString fHilName = "MHillas";
1618 TString fHilNameExt = "MHillasExt";
1619 TString fHilNameSrc = "MHillasSrc";
1620 TString fImgParName = "MNewImagePar";
1621
1622 //-------------------------------------------
1623 // create the tasks which should be executed
1624 //
1625
1626 MReadMarsFile read("Events", filenameData);
1627 read.DisableAutoScheme();
1628
1629
1630 //-----------------------------------------------------------------
1631 // geometry is needed in MHHillas... classes
1632 MGeomCam *fGeom =
1633 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1634
1635 MFCT1SelFinal selfinalgh(fHilNameSrc);
1636 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1637 selfinalgh.SetHadronnessName(fhadronnessName);
1638 selfinalgh.SetName("SelFinalgh");
1639 MContinue contfinalgh(&selfinalgh);
1640 contfinalgh.SetName("ContSelFinalgh");
1641
1642 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
1643 //fillhadnn.SetName("HhadNN");
1644 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
1645 fillhadsc.SetName("HhadSC");
1646 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
1647 fillhadrf.SetName("HhadRF");
1648
1649
1650 MFillH hfill1("MHHillas", fHilName);
1651 hfill1.SetName("HHillas");
1652
1653 MFillH hfill2("MHStarMap", fHilName);
1654 hfill2.SetName("HStarMap");
1655
1656 MFillH hfill3("MHHillasExt", fHilNameSrc);
1657 hfill3.SetName("HHillasExt");
1658
1659 MFillH hfill4("MHHillasSrc", fHilNameSrc);
1660 hfill4.SetName("HHillasSrc");
1661
1662 MFillH hfill5("MHNewImagePar", fImgParName);
1663 hfill5.SetName("HNewImagePar");
1664
1665 MFCT1SelFinal selfinal(fHilNameSrc);
1666 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1667 selfinal.SetHadronnessName(fhadronnessName);
1668 selfinal.SetName("SelFinal");
1669 MContinue contfinal(&selfinal);
1670 contfinal.SetName("ContSelFinal");
1671
1672
1673 //*****************************
1674 // entries in MParList
1675
1676 pliston.AddToList(&tliston);
1677 InitBinnings(&pliston);
1678
1679
1680 //*****************************
1681 // entries in MTaskList
1682
1683 tliston.AddToList(&read);
1684
1685 tliston.AddToList(&contfinalgh);
1686
1687 //tliston.AddToList(&fillhadnn);
1688 tliston.AddToList(&fillhadsc);
1689 tliston.AddToList(&fillhadrf);
1690
1691 tliston.AddToList(&hfill1);
1692 tliston.AddToList(&hfill2);
1693 tliston.AddToList(&hfill3);
1694 tliston.AddToList(&hfill4);
1695 tliston.AddToList(&hfill5);
1696
1697 tliston.AddToList(&contfinal);
1698
1699 //*****************************
1700
1701 //-------------------------------------------
1702 // Execute event loop
1703 //
1704 MProgressBar bar;
1705 MEvtLoop evtloop;
1706 evtloop.SetParList(&pliston);
1707 evtloop.SetProgressBar(&bar);
1708
1709 Int_t maxevents = -1;
1710 //Int_t maxevents = 10000;
1711 if ( !evtloop.Eventloop(maxevents) )
1712 return;
1713
1714 tliston.PrintStatistics(0, kTRUE);
1715
1716
1717
1718 //-------------------------------------------
1719 // Display the histograms
1720 //
1721
1722 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
1723 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
1724 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
1725
1726 pliston.FindObject("MHHillas")->DrawClone();
1727 pliston.FindObject("MHHillasExt")->DrawClone();
1728 pliston.FindObject("MHHillasSrc")->DrawClone();
1729 pliston.FindObject("MHNewImagePar")->DrawClone();
1730 pliston.FindObject("MHStarMap")->DrawClone();
1731
1732 //-------------------------------------------
1733 // fit alpha distribution to get the number of excess events
1734 //
1735
1736 MHHillasSrc* hillasSrc =
1737 (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
1738 TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
1739
1740 MHOnSubtraction onsub;
1741 onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
1742 //-------------------------------------------
1743
1744 DeleteBinnings(&pliston);
1745
1746 gLog << "Macro CT1Analysis : End of Job D" << endl;
1747 gLog << "=======================================================" << endl;
1748 }
1749 //---------------------------------------------------------------------
1750
1751
1752 //---------------------------------------------------------------------
1753 // Job E_EST_UP
1754 //============
1755
1756 // - read MC2.root file
1757 // - select g/h separation method XX
1758 // - optimize energy estimator for events passing final cuts
1759 // - write parameters of energy estimator onto file "energyest_XX.root"
1760 //
1761 // - read ON2.root and MC2.root files
1762 // - update input root file with the estimated energy
1763 // (ON_XX2.root, MC_XX2.root)
1764
1765
1766 if (JobE_EST_UP)
1767 {
1768 gLog << "=====================================================" << endl;
1769 gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl;
1770
1771 gLog << "" << endl;
1772 gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = "
1773 << JobE_EST_UP << ", " << WESTUP << endl;
1774
1775
1776 TString typeON = "ON";
1777 TString typeOFF = "OFF";
1778 TString typeMC = "MC";
1779 TString ext = "2.root";
1780 TString extout = "3.root";
1781
1782 //------------------------------
1783 // name of MC file to be used for optimizing the energy estimator
1784 TString filenameOpt(outPath);
1785 filenameOpt += typeMC;
1786 filenameOpt += ext;
1787 gLog << "filenameOpt = " << filenameOpt << endl;
1788
1789 //------------------------------
1790 // selection of g/h separation method
1791 // and definition of final selections
1792
1793 //TString XX("NN");
1794 //TString XX("SC");
1795 TString XX("RF");
1796 TString fhadronnessName("Had");
1797 fhadronnessName += XX;
1798 gLog << "fhadronnessName = " << fhadronnessName << endl;
1799
1800 // maximum values of the hadronness, |alpha| and dist
1801 Float_t maxhadronness = 0.40;
1802 Float_t maxalpha = 20.0;
1803 Float_t maxdist = 10.0;
1804 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
1805 << maxhadronness << ", " << maxalpha << ", "
1806 << maxdist << endl;
1807
1808 // name of file containing the parameters of the energy estimator
1809 TString energyParName(outPath);
1810 energyParName += "energyest_";
1811 energyParName += XX;
1812 energyParName += ".root";
1813 gLog << "energyParName = " << energyParName << endl;
1814
1815
1816 //------------------------------
1817 // name of ON file to be updated
1818 TString filenameON(outPath);
1819 filenameON += typeON;
1820 filenameON += ext;
1821 gLog << "filenameON = " << filenameON << endl;
1822
1823 // name of OFF file to be updated
1824 TString filenameOFF(outPath);
1825 filenameOFF += typeOFF;
1826 filenameOFF += ext;
1827 gLog << "filenameOFF = " << filenameOFF << endl;
1828
1829 // name of MC file to be updated
1830 TString filenameMC(outPath);
1831 filenameMC += typeMC;
1832 filenameMC += ext;
1833 gLog << "filenameMC = " << filenameMC << endl;
1834
1835 //------------------------------
1836 // name of updated ON file
1837 TString filenameONup(outPath);
1838 filenameONup += typeON;
1839 filenameONup += "_";
1840 filenameONup += XX;
1841 filenameONup += extout;
1842 gLog << "filenameONup = " << filenameONup << endl;
1843
1844 // name of updated OFF file
1845 TString filenameOFFup(outPath);
1846 filenameOFFup += typeOFF;
1847 filenameOFFup += "_";
1848 filenameOFFup += XX;
1849 filenameOFFup += extout;
1850 gLog << "filenameOFFup = " << filenameOFFup << endl;
1851
1852 // name of updated MC file
1853 TString filenameMCup(outPath);
1854 filenameMCup += typeMC;
1855 filenameMCup += "_";
1856 filenameMCup += XX;
1857 filenameMCup += extout;
1858 gLog << "filenameMCup = " << filenameMCup << endl;
1859
1860 //-----------------------------------------------------------
1861
1862 TString fHilName = "MHillas";
1863 TString fHilNameExt = "MHillasExt";
1864 TString fHilNameSrc = "MHillasSrc";
1865 TString fImgParName = "MNewImagePar";
1866
1867 //===========================================================
1868 //
1869 // Optimization of energy estimator
1870 //
1871
1872 TString inpath("");
1873 TString outpath("");
1874 Int_t howMany = 2000;
1875 CT1EEst(inpath, filenameOpt, outpath, energyParName,
1876 fHilName, fHilNameSrc, fhadronnessName,
1877 howMany, maxhadronness, maxalpha, maxdist);
1878
1879 //-----------------------------------------------------------
1880 //
1881 // Read in parameters of energy estimator
1882 //
1883 gLog << "================================================================"
1884 << endl;
1885 gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '"
1886 << energyParName << "'" << endl;
1887 TFile enparam(energyParName);
1888 MMcEnergyEst mcest("MMcEnergyEst");
1889 mcest.Read("MMcEnergyEst");
1890 enparam.Close();
1891
1892 TArrayD parA(5);
1893 TArrayD parB(7);
1894 for (Int_t i=0; i<parA.GetSize(); i++)
1895 parA[i] = mcest.GetCoeff(i);
1896 for (Int_t i=0; i<parB.GetSize(); i++)
1897 parB[i] = mcest.GetCoeff( i+parA.GetSize() );
1898
1899
1900 if (WESTUP)
1901 {
1902 //========== start update ============================================
1903 //
1904 // Update ON, OFF and MC root files with the estimated energy
1905
1906 //---------------------------------------------------
1907 // Update OFF data
1908 //
1909 gLog << "============================================================"
1910 << endl;
1911 gLog << "Macro CT1Analysis.C : update file '" << filenameOFF
1912 << "'" << endl;
1913
1914 MTaskList tlistoff;
1915 MParList plistoff;
1916
1917
1918 // geometry is needed in MHHillas... classes
1919 MGeomCam *fGeom =
1920 (MGeomCam*)plistoff->FindCreateObj("MGeomCamCT1", "MGeomCam");
1921
1922 //-------------------------------------------
1923 // create the tasks which should be executed
1924 //
1925
1926 MReadMarsFile read("Events", filenameOFF);
1927 read.DisableAutoScheme();
1928
1929 //---------------------------
1930 // calculate estimated energy
1931
1932 MEnergyEstParam eest2(fHilName);
1933 eest2.Add(fHilNameSrc);
1934
1935 eest2.SetCoeffA(parA);
1936 eest2.SetCoeffB(parB);
1937
1938 //.......................................................................
1939
1940 MWriteRootFile write(filenameOFFup);
1941
1942 write.AddContainer("MRawRunHeader", "RunHeaders");
1943 write.AddContainer("MTime", "Events");
1944 write.AddContainer("MMcEvt", "Events");
1945 write.AddContainer("ThetaOrig", "Events");
1946 write.AddContainer("MSrcPosCam", "Events");
1947 write.AddContainer("MSigmabar", "Events");
1948 write.AddContainer("MHillas", "Events");
1949 write.AddContainer("MHillasExt", "Events");
1950 write.AddContainer("MHillasSrc", "Events");
1951 write.AddContainer("MNewImagePar", "Events");
1952
1953 //write.AddContainer("HadNN", "Events");
1954 write.AddContainer("HadSC", "Events");
1955 write.AddContainer("HadRF", "Events");
1956
1957 write.AddContainer("MEnergyEst", "Events");
1958
1959 //-----------------------------------------------------------------
1960
1961 MFCT1SelFinal selfinal(fHilNameSrc);
1962 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1963 selfinal.SetHadronnessName(fhadronnessName);
1964 MContinue contfinal(&selfinal);
1965
1966
1967 //*****************************
1968 // entries in MParList
1969
1970 plistoff.AddToList(&tlistoff);
1971 InitBinnings(&plistoff);
1972
1973
1974 //*****************************
1975 // entries in MTaskList
1976
1977 tlistoff.AddToList(&read);
1978 tlistoff.AddToList(&eest2);
1979 tlistoff.AddToList(&write);
1980 tlistoff.AddToList(&contfinal);
1981
1982 //*****************************
1983
1984 //-------------------------------------------
1985 // Execute event loop
1986 //
1987 MProgressBar bar;
1988 MEvtLoop evtloop;
1989 evtloop.SetParList(&plistoff);
1990 evtloop.SetProgressBar(&bar);
1991
1992 Int_t maxevents = -1;
1993 //Int_t maxevents = 1000;
1994 if ( !evtloop.Eventloop(maxevents) )
1995 return;
1996
1997 tlistoff.PrintStatistics(0, kTRUE);
1998 DeleteBinnings(&plistoff);
1999
2000 //---------------------------------------------------
2001
2002 //---------------------------------------------------
2003 // Update ON data
2004 //
2005 gLog << "============================================================"
2006 << endl;
2007 gLog << "Macro CT1Analysis.C : update file '" << filenameON
2008 << "'" << endl;
2009
2010 MTaskList tliston;
2011 MParList pliston;
2012
2013
2014 // geometry is needed in MHHillas... classes
2015 MGeomCam *fGeom =
2016 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2017
2018 //-------------------------------------------
2019 // create the tasks which should be executed
2020 //
2021
2022 MReadMarsFile read("Events", filenameON);
2023 read.DisableAutoScheme();
2024
2025 //---------------------------
2026 // calculate estimated energy
2027
2028 MEnergyEstParam eest2(fHilName);
2029 eest2.Add(fHilNameSrc);
2030
2031 eest2.SetCoeffA(parA);
2032 eest2.SetCoeffB(parB);
2033
2034 //.......................................................................
2035
2036 MWriteRootFile write(filenameONup);
2037
2038 write.AddContainer("MRawRunHeader", "RunHeaders");
2039 write.AddContainer("MTime", "Events");
2040 write.AddContainer("MMcEvt", "Events");
2041 write.AddContainer("ThetaOrig", "Events");
2042 write.AddContainer("MSrcPosCam", "Events");
2043 write.AddContainer("MSigmabar", "Events");
2044 write.AddContainer("MHillas", "Events");
2045 write.AddContainer("MHillasExt", "Events");
2046 write.AddContainer("MHillasSrc", "Events");
2047 write.AddContainer("MNewImagePar", "Events");
2048
2049 //write.AddContainer("HadNN", "Events");
2050 write.AddContainer("HadSC", "Events");
2051 write.AddContainer("HadRF", "Events");
2052
2053 write.AddContainer("MEnergyEst", "Events");
2054
2055 //-----------------------------------------------------------------
2056
2057 MFCT1SelFinal selfinal(fHilNameSrc);
2058 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2059 selfinal.SetHadronnessName(fhadronnessName);
2060 MContinue contfinal(&selfinal);
2061
2062
2063 //*****************************
2064 // entries in MParList
2065
2066 pliston.AddToList(&tliston);
2067 InitBinnings(&pliston);
2068
2069
2070 //*****************************
2071 // entries in MTaskList
2072
2073 tliston.AddToList(&read);
2074 tliston.AddToList(&eest2);
2075 tliston.AddToList(&write);
2076 tliston.AddToList(&contfinal);
2077
2078 //*****************************
2079
2080 //-------------------------------------------
2081 // Execute event loop
2082 //
2083 MProgressBar bar;
2084 MEvtLoop evtloop;
2085 evtloop.SetParList(&pliston);
2086 evtloop.SetProgressBar(&bar);
2087
2088 Int_t maxevents = -1;
2089 //Int_t maxevents = 1000;
2090 if ( !evtloop.Eventloop(maxevents) )
2091 return;
2092
2093 tliston.PrintStatistics(0, kTRUE);
2094 DeleteBinnings(&pliston);
2095
2096 //---------------------------------------------------
2097
2098 //---------------------------------------------------
2099 // Update MC data
2100 //
2101 gLog << "============================================================"
2102 << endl;
2103 gLog << "Macro CT1Analysis.C : update file '" << filenameMC
2104 << "'" << endl;
2105
2106 MTaskList tlistmc;
2107 MParList plistmc;
2108
2109 //-------------------------------------------
2110 // create the tasks which should be executed
2111 //
2112
2113 MReadMarsFile read("Events", filenameMC);
2114 read.DisableAutoScheme();
2115
2116 //---------------------------
2117 // calculate estimated energy
2118
2119 MEnergyEstParam eest2(fHilName);
2120 eest2.Add(fHilNameSrc);
2121
2122 eest2.SetCoeffA(parA);
2123 eest2.SetCoeffB(parB);
2124
2125 //.......................................................................
2126
2127 MWriteRootFile write(filenameMCup);
2128
2129 write.AddContainer("MRawRunHeader", "RunHeaders");
2130 write.AddContainer("MTime", "Events");
2131 write.AddContainer("MMcEvt", "Events");
2132 write.AddContainer("ThetaOrig", "Events");
2133 write.AddContainer("MSrcPosCam", "Events");
2134 write.AddContainer("MSigmabar", "Events");
2135 write.AddContainer("MHillas", "Events");
2136 write.AddContainer("MHillasExt", "Events");
2137 write.AddContainer("MHillasSrc", "Events");
2138 write.AddContainer("MNewImagePar", "Events");
2139
2140 //write.AddContainer("HadNN", "Events");
2141 write.AddContainer("HadSC", "Events");
2142 write.AddContainer("HadRF", "Events");
2143
2144 write.AddContainer("MEnergyEst", "Events");
2145
2146 //-----------------------------------------------------------------
2147
2148 MFCT1SelFinal selfinal(fHilNameSrc);
2149 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2150 selfinal.SetHadronnessName(fhadronnessName);
2151 MContinue contfinal(&selfinal);
2152
2153
2154 //*****************************
2155 // entries in MParList
2156
2157 plistmc.AddToList(&tlistmc);
2158 InitBinnings(&plistmc);
2159
2160
2161 //*****************************
2162 // entries in MTaskList
2163
2164 tlistmc.AddToList(&read);
2165 tlistmc.AddToList(&eest2);
2166 tlistmc.AddToList(&write);
2167 tlistmc.AddToList(&contfinal);
2168
2169 //*****************************
2170
2171 //-------------------------------------------
2172 // Execute event loop
2173 //
2174 MProgressBar bar;
2175 MEvtLoop evtloop;
2176 evtloop.SetParList(&plistmc);
2177 evtloop.SetProgressBar(&bar);
2178
2179 Int_t maxevents = -1;
2180 //Int_t maxevents = 1000;
2181 if ( !evtloop.Eventloop(maxevents) )
2182 return;
2183
2184 tlistmc.PrintStatistics(0, kTRUE);
2185 DeleteBinnings(&plistmc);
2186
2187
2188 //========== end update ============================================
2189 }
2190
2191 enparam.Close();
2192
2193 gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl;
2194 gLog << "=======================================================" << endl;
2195 }
2196 //---------------------------------------------------------------------
2197
2198
2199 //---------------------------------------------------------------------
2200 // Job F_XX
2201 //=========
2202
2203 // - select g/h separation method XX
2204 // - read MC_XX2.root file
2205   // - calculate eff. collection area
2206 // - read ON_XX2.root file
2207 // - apply final cuts
2208 // - calculate flux
2209 // - write root file for ON data after final cuts (ON_XX3.root))
2210
2211
2212 if (JobF_XX)
2213 {
2214 gLog << "=====================================================" << endl;
2215 gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
2216
2217 gLog << "" << endl;
2218 gLog << "Macro CT1Analysis : JobF_XX, WXX = "
2219 << JobF_XX << ", " << WXX << endl;
2220
2221 // type of data to be analysed
2222 //TString typeData = "ON";
2223 //TString typeData = "OFF";
2224 TString typeData = "MC";
2225 gLog << "typeData = " << typeData << endl;
2226
2227 TString typeMC = "MC";
2228 TString ext = "3.root";
2229 TString extout = "4.root";
2230
2231 //------------------------------
2232 // selection of g/h separation method
2233 // and definition of final selections
2234
2235 //TString XX("NN");
2236 //TString XX("SC");
2237 TString XX("RF");
2238 TString fhadronnessName("Had");
2239 fhadronnessName += XX;
2240 gLog << "fhadronnessName = " << fhadronnessName << endl;
2241
2242 // maximum values of the hadronness, |ALPHA| and DIST
2243 Float_t maxhadronness = 0.40;
2244 Float_t maxalpha = 20.0;
2245 Float_t maxdist = 10.0;
2246 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2247 << maxhadronness << ", " << maxalpha << ", "
2248 << maxdist << endl;
2249
2250
2251 //------------------------------
2252 // name of MC file to be used for calculating the eff. collection areas
2253 TString filenameArea(outPath);
2254 filenameArea += typeMC;
2255 filenameArea += "_";
2256 filenameArea += XX;
2257 filenameArea += ext;
2258 gLog << "filenameArea = " << filenameArea << endl;
2259
2260 //------------------------------
2261 // name of file containing the eff. collection areas
2262 TString collareaName(outPath);
2263 collareaName += "area_";
2264 collareaName += XX;
2265 collareaName += ".root";
2266 gLog << "collareaName = " << collareaName << endl;
2267
2268 //------------------------------
2269 // name of data file to be analysed
2270 TString filenameData(outPath);
2271 filenameData += typeData;
2272 filenameData += "_";
2273 filenameData += XX;
2274 filenameData += ext;
2275 gLog << "filenameData = " << filenameData << endl;
2276
2277 //------------------------------
2278 // name of output data file (after the final cuts)
2279 TString filenameDataout(outPath);
2280 filenameDataout += typeData;
2281 filenameDataout += "_";
2282 filenameDataout += XX;
2283 filenameDataout += extout;
2284 gLog << "filenameDataout = " << filenameDataout << endl;
2285
2286
2287 //====================================================================
2288 gLog << "-----------------------------------------------" << endl;
2289 gLog << "Start calculation of effective collection areas" << endl;
2290 MParList parlist;
2291 MTaskList tasklist;
2292
2293 //---------------------------------------
2294 // Setup the tasks to be executed
2295 //
2296 MReadMarsFile reader("Events", filenameArea);
2297 reader.DisableAutoScheme();
2298
2299 MFCT1SelFinal cuthadrons;
2300 cuthadrons.SetHadronnessName(fhadronnessName);
2301 cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
2302
2303 MContinue conthadrons(&cuthadrons);
2304
2305 //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
2306 //MHMcCT1CollectionArea* collarea;
2307
2308 MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
2309 filler.SetName("CollectionArea");
2310
2311 //********************************
2312 // entries in MParList
2313
2314 parlist.AddToList(&tasklist);
2315 InitBinnings(&parlist);
2316 //parlist.AddToList(collarea);
2317
2318 //********************************
2319 // entries in MTaskList
2320
2321 tasklist.AddToList(&reader);
2322 tasklist.AddToList(&conthadrons);
2323 tasklist.AddToList(&filler);
2324
2325 //********************************
2326
2327 //-----------------------------------------
2328 // Execute event loop
2329 //
2330 MEvtLoop magic;
2331 magic.SetParList(&parlist);
2332
2333 MProgressBar bar;
2334 magic.SetProgressBar(&bar);
2335 if (!magic.Eventloop())
2336 return;
2337
2338 tasklist.PrintStatistics(0, kTRUE);
2339
2340 // Calculate effective collection areas
2341 // and display the histograms
2342 //
2343
2344 MHMcCT1CollectionArea *collarea =
2345 (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
2346
2347 collarea->CalcEfficiency();
2348 collarea->DrawClone("lego");
2349
2350 //---------------------------------------------
2351 // Write histograms to a file
2352 //
2353
2354 TFile f(collareaName, "RECREATE");
2355 collarea->GetHist()->Write();
2356 collarea->GetHAll()->Write();
2357 collarea->GetHSel()->Write();
2358 f.Close();
2359
2360 //delete collarea;
2361
2362 gLog << "Calculation of effective collection areas done" << endl;
2363 gLog << "-----------------------------------------------" << endl;
2364 //------------------------------------------------------------------
2365
2366
2367 //*************************************************************************
2368 //
2369 // Analyse the data
2370 //
2371
2372 MTaskList tliston;
2373 MParList pliston;
2374
2375 // geometry is needed in MHHillas... classes
2376 MGeomCam *fGeom =
2377 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2378
2379 TString fHilName = "MHillas";
2380 TString fHilNameExt = "MHillasExt";
2381 TString fHilNameSrc = "MHillasSrc";
2382 TString fImgParName = "MNewImagePar";
2383
2384 //-------------------------------------------
2385 // create the tasks which should be executed
2386 //
2387
2388 MReadMarsFile read("Events", filenameData);
2389 read.DisableAutoScheme();
2390
2391 //.......................................................................
2392
2393
2394 MWriteRootFile write(filenameDataout);
2395
2396 write.AddContainer("MRawRunHeader", "RunHeaders");
2397 write.AddContainer("MTime", "Events");
2398 write.AddContainer("MMcEvt", "Events");
2399 write.AddContainer("ThetaOrig", "Events");
2400 write.AddContainer("MSrcPosCam", "Events");
2401 write.AddContainer("MSigmabar", "Events");
2402 write.AddContainer("MHillas", "Events");
2403 write.AddContainer("MHillasExt", "Events");
2404 write.AddContainer("MHillasSrc", "Events");
2405 write.AddContainer("MNewImagePar", "Events");
2406
2407 //write.AddContainer("HadNN", "Events");
2408 write.AddContainer("HadSC", "Events");
2409 write.AddContainer("HadRF", "Events");
2410
2411 write.AddContainer("MEnergyEst", "Events");
2412
2413
2414 //-----------------------------------------------------------------
2415 // geometry is needed in MHHillas... classes
2416 MGeomCam *fGeom =
2417 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2418
2419 MFCT1SelFinal selfinalgh(fHilNameSrc);
2420 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2421 selfinalgh.SetHadronnessName(fhadronnessName);
2422 selfinalgh.SetName("SelFinalgh");
2423 MContinue contfinalgh(&selfinalgh);
2424 contfinalgh.SetName("ContSelFinalgh");
2425
2426 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
2427 //fillhadnn.SetName("HhadNN");
2428 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
2429 fillhadsc.SetName("HhadSC");
2430 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
2431 fillhadrf.SetName("HhadRF");
2432
2433
2434 MFillH hfill1("MHHillas", fHilName);
2435 hfill1.SetName("HHillas");
2436
2437 MFillH hfill2("MHStarMap", fHilName);
2438 hfill2.SetName("HStarMap");
2439
2440 MFillH hfill3("MHHillasExt", fHilNameSrc);
2441 hfill3.SetName("HHillasExt");
2442
2443 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2444 hfill4.SetName("HHillasSrc");
2445
2446 MFillH hfill5("MHNewImagePar", fImgParName);
2447 hfill5.SetName("HNewImagePar");
2448
2449 MFCT1SelFinal selfinal(fHilNameSrc);
2450 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2451 selfinal.SetHadronnessName(fhadronnessName);
2452 selfinal.SetName("SelFinal");
2453 MContinue contfinal(&selfinal);
2454 contfinal.SetName("ContSelFinal");
2455
2456
2457 //*****************************
2458 // entries in MParList
2459
2460 pliston.AddToList(&tliston);
2461 InitBinnings(&pliston);
2462
2463
2464 //*****************************
2465 // entries in MTaskList
2466
2467 tliston.AddToList(&read);
2468
2469 tliston.AddToList(&contfinalgh);
2470
2471 if (WXX)
2472 tliston.AddToList(&write);
2473
2474 //tliston.AddToList(&fillhadnn);
2475 tliston.AddToList(&fillhadsc);
2476 tliston.AddToList(&fillhadrf);
2477
2478 tliston.AddToList(&hfill1);
2479 tliston.AddToList(&hfill2);
2480 tliston.AddToList(&hfill3);
2481 tliston.AddToList(&hfill4);
2482 tliston.AddToList(&hfill5);
2483
2484 tliston.AddToList(&contfinal);
2485
2486 //*****************************
2487
2488 //-------------------------------------------
2489 // Execute event loop
2490 //
2491 MProgressBar bar;
2492 MEvtLoop evtloop;
2493 evtloop.SetParList(&pliston);
2494 evtloop.SetProgressBar(&bar);
2495
2496 Int_t maxevents = -1;
2497 //Int_t maxevents = 10000;
2498 if ( !evtloop.Eventloop(maxevents) )
2499 return;
2500
2501 tliston.PrintStatistics(0, kTRUE);
2502
2503
2504 //-------------------------------------------
2505 // Display the histograms
2506 //
2507
2508 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
2509 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2510 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2511
2512 pliston.FindObject("MHHillas")->DrawClone();
2513 pliston.FindObject("MHHillasExt")->DrawClone();
2514 pliston.FindObject("MHHillasSrc")->DrawClone();
2515 pliston.FindObject("MHNewImagePar")->DrawClone();
2516 pliston.FindObject("MHStarMap")->DrawClone();
2517
2518 DeleteBinnings(&pliston);
2519
2520 gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
2521 gLog << "=======================================================" << endl;
2522 }
2523 //---------------------------------------------------------------------
2524
2525}
2526
2527
2528
Note: See TracBrowser for help on using the repository browser.