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

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