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

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