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

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