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

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