source: trunk/MagicSoft/Mars/mtemp/mmpi/hit.cc@ 6985

Last change on this file since 6985 was 6703, checked in by mazin, 21 years ago
*** empty log message ***
File size: 25.2 KB
Line 
1#include <TROOT.h>
2#include <TClass.h>
3#include <TSystem.h>
4#include <TGClient.h>
5#include <TApplication.h>
6#include <TObjectTable.h>
7
8#include "MHMatrix.h"
9#include "MParList.h"
10#include "MTaskList.h"
11#include "MWriteRootFile.h"
12#include "MReadMarsFile.h"
13#include "MReadReports.h"
14#include "MApplySupercuts.h"
15#include "MHFindSignificance.h"
16#include "MEvtLoop.h"
17#include "MProgressBar.h"
18#include "MContinue.h"
19#include "MGeomApply.h"
20#include "MHadronness.h"
21#include "MFilterList.h"
22#include "MF.h"
23#include "MFindSupercutsONOFFThetaLoop.h"
24#include "MRanForestFill.h"
25
26#include "MArgs.h"
27#include "MArray.h"
28#include "MDirIter.h"
29
30#include "MStatusDisplay.h"
31
32#include "MSequence.h"
33#include "MJStar.h"
34#include "MH3.h"
35#include "MRFEnergyEst.h"
36
37#include "MRanForestCalc.h"
38#include "MTaskInteractive.h"
39#include "MHillas.h"
40//#include "MFHadAlpha.h"
41
42#include "MLog.h"
43#include "MLogManip.h"
44#include "MPrint.h"
45
46using namespace std;
47
48
49Int_t PreProcess(MParList *plist);
50Int_t Process();
51
52static void StartUpMessage()
53{
54 gLog << all << endl;
55
56 // 1 2 3 4 5 6
57 // 123456789012345678901234567890123456789012345678901234567890
58 gLog << "========================================================" << endl;
59 gLog << " Melibea - MARS V" << MARSVER << endl;
60 gLog << " MARS -- Supercuts Hillas Image Treatment" << endl;
61 gLog << " Compiled on <" << __DATE__ << ">" << endl;
62 gLog << " Using ROOT v" << ROOTVER << endl;
63 gLog << "========================================================" << endl;
64
65 gLog << endl;
66}
67
68static void Usage()
69{
70 // 1 2 3 4 5 6 7 8
71 // 12345678901234567890123456789012345678901234567890123456789012345678901234567890
72 gLog << all << endl;
73 gLog << "Sorry the usage is:" << endl;
74
75
76}
77
78int main(int argc, char **argv)
79{
80 StartUpMessage();
81
82 //
83 // Evaluate arguments
84 //
85 MArgs arg(argc, argv, kTRUE);
86
87 if (arg.HasOnly("-V") || arg.HasOnly("--version"))
88 return 0;
89
90 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
91 {
92 Usage();
93 return -1;
94 }
95
96 gLog.Setup(arg);
97
98 const Bool_t kOptimize = arg.HasOnlyAndRemove("--optimize");
99 const Bool_t kSupercuts = arg.HasOnlyAndRemove("--supercuts");
100 const Bool_t kOverwrite = arg.HasOnlyAndRemove("-f");
101 const Bool_t kPrintSeq = arg.HasOnlyAndRemove("--print-seq");
102 const Bool_t kBatch = arg.HasOnlyAndRemove("-b");
103 const Bool_t kPlot = arg.HasOnlyAndRemove("-p");
104 const Bool_t kMc = arg.HasOnlyAndRemove("--mc");
105 const Bool_t kUseSeq = !arg.HasOnlyAndRemove("--nouseseq");
106
107 // Variable that decides whether the optimized cuts are used
108 // in the test sample.
109 const Bool_t TestParams = !arg.HasOnlyAndRemove("--opttest");
110 const Bool_t CombineCosThetaBinsForTrainSample = arg.HasOnlyAndRemove("--combcosthtrain");
111 const Bool_t CombineCosThetaBinsForTestSample = arg.HasOnlyAndRemove("--combcosthtest");
112 const Double_t whichfractiontrain = arg.GetFloatAndRemove("--ontrainfrac=",0.5);
113 const Double_t whichfractiontest = arg.GetFloatAndRemove("--ontestfrac=",0.5);
114 const Double_t whichfractiontrainOFF = arg.GetFloatAndRemove("--offtrainfrac=",0.5);
115 const Double_t whichfractiontestOFF = arg.GetFloatAndRemove("--offtestfrac=",0.5);
116 const Double_t gammaeff = arg.GetFloatAndRemove("--geff=",0.6);
117 const Double_t alphasig = arg.GetFloatAndRemove("--alphasig=",6);
118 const Double_t alphabkgmin = arg.GetFloatAndRemove("--alphabgmin=",20);
119 const Double_t alphabkgmax = arg.GetFloatAndRemove("--alphabgmax=",80);
120 const Double_t SizeLow = arg.GetFloatAndRemove("--sizelow=",60);
121 const Double_t SizeUp = arg.GetFloatAndRemove("--sizeup=",1000000);
122 const Double_t LeakageMax = arg.GetFloatAndRemove("--leakagemax=",0.25);
123 const Double_t DistMax = arg.GetFloatAndRemove("--distmax=",1.0);
124 const Double_t DistMin = arg.GetFloatAndRemove("--distmin=",0.2);
125 const Double_t fThetaMax = arg.GetFloatAndRemove("--zdmax=",30.);
126 const Double_t fThetaMin = arg.GetFloatAndRemove("--zdmin=",0.);
127 const Int_t NAlphaBins = arg.GetIntAndRemove("--nalphabin=",30);
128 const Double_t AlphaBinLow = arg.GetFloatAndRemove("--alphabinlow=",0);
129 const Double_t AlphaBinUp = arg.GetFloatAndRemove("--alphabinup=",90);
130 const Bool_t NormFactorFromAlphaBkg = !arg.HasOnlyAndRemove("--normbgmeth");
131 const Bool_t UseFittedQuantities = !arg.HasOnlyAndRemove("--usefittedq");
132 const Bool_t NotUseTheta = !arg.HasOnlyAndRemove("--notuseth");
133 const Bool_t UseStaticCuts = !arg.HasOnlyAndRemove("--staticcuts");
134 const Bool_t ReadInitParamsFromAsciiFile = !arg.HasOnlyAndRemove("--ascii");
135 const Bool_t kPrintFiles = !arg.HasOnlyAndRemove("--printseq");
136 const Int_t NInitSCPar = arg.GetIntAndRemove("--ninitscpar=",104);
137 const Int_t degree = arg.GetIntAndRemove("--degree=",2);
138 const TString InitSCParamAsciiFile = arg.GetStringAndRemove("--inscparascii=","mtemp/mmpi/asciifiles/StartingValuesForSmallSizes1.txt");
139 const TString ONDataFilename = arg.GetStringAndRemove("--inondataopt=","./*CrabNebula*.root");
140 const TString OFFDataFilename = arg.GetStringAndRemove("--inoffdataopt=","./*OffCrab*.root");
141 const TString PathForFiles = arg.GetStringAndRemove("--outopt=","./");
142
143 const TString kInpath = arg.GetStringAndRemove("--ind=", "");
144 const TString kOutpath = arg.GetStringAndRemove("--out=", "");
145 const TString kSCfile = arg.GetStringAndRemove("--scf=", "");
146 const TString kSource = arg.GetStringAndRemove("--src=", "");
147 const TString kDate = arg.GetStringAndRemove("--date=", "");
148 const TString kRFTreeFile = arg.GetStringAndRemove("--rftree=", "");
149 const Double_t kAlphaMin = arg.GetFloatAndRemove("--alphamin=",-1);
150 const Double_t kAlphaMax = arg.GetFloatAndRemove("--alphamax=",-1);
151 const Double_t kAlphaSig = arg.GetFloatAndRemove("--alphasignal=",-1);
152 const Int_t kDegree = arg.GetIntAndRemove("--poldeg=",-1);
153 const Int_t kNumTrees = arg.GetIntAndRemove("--ntrees=",100);
154
155
156 if (arg.GetNumOptions()>0)
157 {
158 gLog << warn << "WARNING - Unknown commandline options..." << endl;
159 arg.Print("options");
160 gLog << endl;
161 return -1;
162 }
163
164 if (arg.GetNumArguments()!=1)
165 {
166 Usage();
167 return -1;
168 }
169 //
170 // Setup sequence file and check for its existance
171 //
172 const TString kSequence = arg.GetArgumentStr(0);
173
174 if (gSystem->AccessPathName(kSequence, kFileExists))
175 {
176 gLog << err << "Sorry, sequence file '" << kSequence << "' doesn't exist." << endl;
177 return -1;
178 }
179
180 //
181 // Setup sequence and check its validity
182 //
183 MSequence seq(kSequence);
184 if (kPrintSeq)
185 {
186 gLog << all;
187 gLog.Separator(kSequence);
188 seq.Print();
189 gLog << endl;
190 }
191 if (!seq.IsValid())
192 {
193 gLog << err << "Sequence read but not valid!" << endl << endl;
194 return -1;
195 }
196
197 //
198 // Process print options
199 //
200
201 MDirIter iter;
202
203 if(kUseSeq) {
204 gLog << inf;
205 gLog << "Calculate image parameters from sequence ";
206 gLog << seq.GetName() << endl;
207 gLog << endl;
208
209
210 const Int_t n0 = seq.SetupDatRuns(iter, kInpath, "I");
211 const Int_t n1 = seq.GetNumDatRuns();
212
213 if (n0==0)
214 {
215 gLog << err << "ERROR - No input files of sequence found!" << endl;
216 return kFALSE;
217 }
218 if (n0!=n1)
219 {
220 gLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl;
221 return kFALSE;
222 }
223
224 }
225
226
227
228 //
229 // Initialize root
230 //
231 MArray::Class()->IgnoreTObjectStreamer();
232 MParContainer::Class()->IgnoreTObjectStreamer();
233
234 TApplication app("Star", &argc, argv);
235 if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch)
236 {
237 gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl;
238 return 1;
239 }
240
241
242 // **************************************************** //
243 // Setup flags for optimization //
244 // **************************************************** //
245
246
247 // Name of the root file where alpha distributions, TTree objects
248 // with info about the events and cuts applied and info support histograms
249 // will be stored.
250 // Write only the name of the file. The Path
251 // is the one defined previously
252
253 TString RootFilename = kOutpath + "RootFileDynCuts.root";
254 TString ONFiles = ONDataFilename + "*_I_*.root";
255 TString OFFFiles = OFFDataFilename + "*_I_*.root";
256 // Vector containing the theta bins in which data (ON/OFF train/test)
257 // will be divided. Actually this vector contains the cosinus of
258 // these theta bins. The dimension of the vector is N+1, where
259 // N is the number of theta bins intended. The first component of the
260 // vector is the low bin edge of the first bin, and the last
261 // vector component the upper bin edge of the last bin.
262
263 TArrayD CosThetaRangeVector(2);
264 CosThetaRangeVector[1] = 0.866; //30
265 CosThetaRangeVector[0] = 0.5; // 60
266
267 // Object of MCT1FindSupercutsONOFFThetaLoop created, data that was specified
268 // above is introduced and ... and the party starts.
269
270 MFindSupercutsONOFFThetaLoop FindSupercuts("MFindSupercutsONOFFThetaLoop",
271 "Optimizer for the supercuts");
272
273
274 // ****************************************************** //
275 // Calculate and/or apply dynamical cuts //
276 // ****************************************************** //
277
278 // Hillas data file
279 //TString datafile = kInpath + "*_I_CrabNebula_E.root";
280
281 TString datafile;
282 if (kUseSeq)
283 datafile = kInpath + kDate + "*_I_" + kSource + "*.root";
284 else
285 datafile = kInpath + "*" + kSource + "*.root";
286
287 MReadMarsFile read("Events");
288 read.DisableAutoScheme();
289 if (kUseSeq)
290 read.AddFiles(iter);
291 else
292 read.AddFile(datafile);
293
294 cout << " datafile is : " << datafile << endl;
295
296 MReadReports readreal;
297 readreal.AddTree("Events", "MTime.", kTRUE);
298 readreal.AddTree("Drive");
299 readreal.AddTree("EffectiveOnTime");
300 if (kUseSeq)
301 readreal.AddFiles(iter);
302 else
303 readreal.AddFile(datafile);
304
305// readreal.AddToBranchList("MReportDrive.*");
306// readreal.AddToBranchList("MRawEvtHeader.*");
307// readreal.AddToBranchList("MEffectiveOnTime.*");
308
309
310 MGeomApply geomapl;
311
312 MApplySupercuts appcuts;
313 appcuts.SetSCFilename(kSCfile);
314
315 MFilterList flist;
316 MF hadrf("MHadronness.fHadronness<0.5");
317 MF sizef("MHillas.fSize>300");
318 MF leakmax("MNewImagePar.fLeakage1<0.25");
319 MF distmax("MHillasSrc.fDist<300");
320 MF distmin("MHillasSrc.fDist>60");
321 MF widthmin("MHillas.fWidth>15");
322 MF widthmax("MHillas.fWidth<120");
323
324 TString ThetaCutMinString ("MPointingPos.fZd>"); // new! // for magic
325 ThetaCutMinString += fThetaMin;
326 MContinue ThetaCutMin(ThetaCutMinString);
327 ThetaCutMin.SetInverted();
328
329 TString ThetaCutMaxString ("MPointingPos.fZd<"); // new! // for magic
330 ThetaCutMaxString += fThetaMax;
331 MContinue ThetaCutMax(ThetaCutMaxString);
332 ThetaCutMax.SetInverted();
333
334 MContinue cont("(MHillasSrc.fDist<300) && (MHillasSrc.fDist>60) && (MNewImagePar.fLeakage1<0.25) && (MHillas.fSize>60.)");
335 cont.SetInverted();
336
337
338 MRFEnergyEst rfs;
339 //rfs.SetFile("/datm3/data/RFEnergyEst/EForestsBeforeCutsLZA.root");
340 //rfs.SetFile("/datm3/data/RFEnergyEst/EForests19990101_10003_I_MCGamTrainHZA_E_10_5.root");
341 rfs.SetFile("/home/pcmagic16/mazin/mars/files/EForestsHZA.root");
342
343
344 //TString outname = Form("%s{s/_I_/_Q_}", kOutpath.Data());
345 MWriteRootFile &write=*new MWriteRootFile(2, Form("%s{s/_I_/_Q_}", kOutpath.Data()), "RECREATE");
346 // Write the Events tree
347 write.AddContainer("MHillas", "Events");
348 write.AddContainer("MHillasExt", "Events");
349 write.AddContainer("MHillasSrc", "Events");
350 write.AddContainer("MImagePar", "Events");
351 write.AddContainer("MNewImagePar", "Events");
352 write.AddContainer("MRawEvtHeader", "Events", kFALSE);
353 write.AddContainer("MPointingPos", "Events");
354 write.AddContainer("MHadronness","Events", kFALSE);
355 write.AddContainer("MEnergyEst","Events", kFALSE);
356
357 if(kMc)
358 {
359 // Monte Carlo
360 write.AddContainer("MMcEvt", "Events");
361 write.AddContainer("MMcTrig", "Events");
362 // Monte Carlo Headers
363 write.AddContainer("MMcTrigHeader", "RunHeaders");
364 write.AddContainer("MMcConfigRunHeader", "RunHeaders");
365 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
366 write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
367 }
368 else
369 {
370 write.AddContainer("MTime", "Events");
371 // Run Header
372 write.AddContainer("MRawRunHeader", "RunHeaders");
373 write.AddContainer("MBadPixelsCam", "RunHeaders", kFALSE);
374 write.AddContainer("MGeomCam", "RunHeaders", kFALSE);
375 //write.AddContainer("MObservatory", "RunHeaders");
376 // Drive
377 //write.AddContainer("MSrcPosCam", "Drive");
378 write.AddContainer("MReportDrive", "Drive", kFALSE);
379 write.AddContainer("MTimeDrive", "Drive", kFALSE);
380 // Effective On Time
381 write.AddContainer("MEffectiveOnTime", "EffectiveOnTime", kFALSE);
382 write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime", kFALSE);
383 }
384
385
386/*
387 write.AddContainer("MHillas","Events");
388 write.AddContainer("MHillasExt","Events");
389 write.AddContainer("MHillasSrc","Events");
390 write.AddContainer("MImagePar","Events");
391 write.AddContainer("MNewImagePar","Events");
392 write.AddContainer("MPointingPos","Events");
393 write.AddContainer("MRawEvtHeader","Events");
394 write.AddContainer("MTime","Events", kFALSE);
395 write.AddContainer("MHadronness","Events");
396 write.AddContainer("MEnergyEst","Events");
397 // Write the RunHeader tree
398 write.AddContainer("MRawRunHeader","RunHeaders", kFALSE);
399 write.AddContainer("MBadPixelsCam","RunHeaders",kFALSE);
400 write.AddContainer("MGeomCam","RunHeaders",kFALSE);
401 // Write the Drive tree
402 write.AddContainer("MReportDrive","Drive",kFALSE);
403 write.AddContainer("MTimeDrive","Drive",kFALSE);
404 // Write the EffectiveOnTime tree
405 write.AddContainer("MEffectiveOnTime","EffectiveOnTime",kFALSE);
406 write.AddContainer("MTimeEffectiveOnTime","EffectiveOnTime",kFALSE);
407
408
409 write.AddContainer("MMcConfigRunHeader", "RunHeaders", kFALSE);
410 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
411 write.AddContainer("MMcFadcHeader", "RunHeaders", kFALSE);
412 write.AddContainer("MMcTrigHeader", "RunHeaders", kFALSE);
413 write.AddContainer("MMcEvt", "Events", kFALSE);
414 write.AddContainer("MMcTrig", "Events", kFALSE);
415*/
416 cout << datafile << endl;
417 MParList plist;
418 MTaskList tlist;
419 plist.AddToList(&tlist);
420
421 MProgressBar bar;
422
423 switch(kSupercuts)
424 {
425 case kTRUE:
426 switch(kOptimize){
427 case kTRUE:
428 //cout << ONFiles << endl;
429 FindSupercuts.SetPathForFiles(PathForFiles);
430 FindSupercuts.SetDataONOFFRootFilenames(ONFiles, OFFFiles);
431 FindSupercuts.SetFractionTrainTestOnOffEvents(whichfractiontrain,
432 whichfractiontest,
433 whichfractiontrainOFF,
434 whichfractiontestOFF);
435 FindSupercuts.SetGammaEfficiency(gammaeff);
436 FindSupercuts.SetAlphaSig(alphasig);
437 // Bkg alpha region is set
438 FindSupercuts.SetAlphaBkgMin(alphabkgmin);
439 FindSupercuts.SetAlphaBkgMax(alphabkgmax);
440 // alpha bkg and signal region set in object FindSupercuts
441 // are re-checked in order to be sure that make sense
442 FindSupercuts.CheckAlphaSigBkg();
443 // Degree of the polynomials used to fit the ON OFF data is set
444 FindSupercuts.SetDegreeON(degree);
445 FindSupercuts.SetDegreeOFF(degree);
446 // binning for alpha plots is defined
447 FindSupercuts.SetAlphaPlotBinining(NAlphaBins, AlphaBinLow, AlphaBinUp);
448 // Size range is defined
449 FindSupercuts.SetSizeRange(SizeLow, SizeUp);
450 FindSupercuts.SetFilters(LeakageMax, DistMax, DistMin);
451 FindSupercuts.SetNormFactorFromAlphaBkg(NormFactorFromAlphaBkg);
452 FindSupercuts.SetUseFittedQuantities(UseFittedQuantities);
453 FindSupercuts.SetVariableUseStaticCuts(UseStaticCuts);
454 FindSupercuts.SetVariableNotUseTheta(NotUseTheta);
455 FindSupercuts.SetReadMatricesFromFile(kFALSE);// normal case kFALSE
456 FindSupercuts.SetTrainParameters(kTRUE);
457 FindSupercuts.SetSkipOptimization(kFALSE);
458
459 FindSupercuts.SetUseInitialSCParams(kTRUE);
460 FindSupercuts.SetTestParameters(TestParams);
461 FindSupercuts.SetHadronnessName("MHadSC");
462 FindSupercuts.SetHadronnessNameOFF("MHadOFFSC");
463 FindSupercuts.SetAlphaDistributionsRootFilename (RootFilename);
464 FindSupercuts.SetCosThetaRangeVector (CosThetaRangeVector);
465
466// energy estimation from Thomas Hengstebeck
467
468
469 // Names for all root files (matrices, alpha distributions...) are created
470 FindSupercuts.SetALLNames();
471
472 if(ReadInitParamsFromAsciiFile)
473 {
474 // Initial SC Parameters and steps are retrieved from Ascii file
475 if(!FindSupercuts.ReadSCParamsFromAsciiFile(InitSCParamAsciiFile,NInitSCPar))
476 {
477 cout << "Initial SC Parameters could not be read from Ascii file "
478 << InitSCParamAsciiFile << endl
479 << "Aborting execution of macro... " << endl;
480 return -1;
481 }
482 }
483
484 // Finally loop over all theta bins defined is executed
485 if (!FindSupercuts.LoopOverThetaRanges())
486 {
487 cout << "Function MFindSupercutsONOFFThetaLoop::LoopOverThetaRanges()" << endl
488 << "could not be performed" << endl;
489 }
490
491 // Nex and Significance are computed vs alphasig
492 if (!FindSupercuts.ComputeNexSignificanceVSAlphaSig())
493 {
494 cout << "Function MFindSupercutsONOFFThetaLoop::ComputeNexSignificanceVSAlphaSig()" << endl
495 << "could not be performed" << endl;
496 }
497
498 // Several theta bins are combined to produced a single alpha plot (for train and test)
499 // with single Nex and significances
500 if (CombineCosThetaBinsForTrainSample || CombineCosThetaBinsForTestSample)
501 {
502 if(!FindSupercuts.ComputeOverallSignificance(CombineCosThetaBinsForTrainSample,
503 CombineCosThetaBinsForTestSample))
504 {
505 cout << "Function MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance" << endl
506 << "could not be performed" << endl;
507 }
508 }
509 break;
510 case kFALSE:
511 if(kAlphaMin!=-1)
512 appcuts.SetBGAlphaMin(kAlphaMin);
513 if(kAlphaMax!=-1)
514 appcuts.SetBGAlphaMax(kAlphaMax);
515 if(kAlphaSig!=-1)
516 appcuts.SetAlphaSig(kAlphaSig);
517 if(kDegree!=-1)
518 appcuts.SetPolDegree(kDegree);
519 if(kPlot)
520 appcuts.SetPlotON();
521 appcuts.SetPlotSizeLow(SizeLow);
522 appcuts.SetPlotSizeUp(SizeUp);
523 cout << "alphasignal=" << kAlphaSig << endl;
524 //flist.AddToList(&hadrf);
525 flist.AddToList(&leakmax);
526 flist.AddToList(&distmax);
527 flist.AddToList(&distmin);
528 flist.AddToList(&widthmin);
529 flist.AddToList(&widthmax);
530 flist.AddToList(&sizef);
531 appcuts.SetFilter(&flist);
532 //write.SetFilter(&flist);
533
534 tlist.AddToList(kMc ? (MTask*)&read : (MTask*)&readreal);
535 cout << " is montecarlo ? " << kMc << endl;
536
537 tlist.AddToList(&geomapl);
538 tlist.AddToList(&ThetaCutMin, "Events");
539 tlist.AddToList(&ThetaCutMax, "Events");
540 tlist.AddToList(&flist, "Events");
541 tlist.AddToList(&appcuts, "Events");
542 tlist.AddToList(&rfs, "Events");
543 tlist.AddToList(&cont, "Events");
544
545 //tlist.AddToList(&hsigma);
546MPrint print1("MTimeDrive");
547MPrint print2("MEffectiveOnTime");
548// tlist.AddToList(&print1, "Drive");
549// tlist.AddToList(&print2, "EffectiveOnTime");
550
551 tlist.AddToList(&write);
552 MEvtLoop shitloop;
553 shitloop.SetParList(&plist);
554 //
555 // Start to loop over all events
556 //
557 MProgressBar bar;
558 shitloop.SetProgressBar(&bar);
559
560 if (!shitloop.Eventloop())
561 return -1;
562
563 /*
564
565 if (!shitloop.PreProcess())
566 return -1;
567
568 while (tlist.Process())
569 {
570 //MHadronness *fH = (MHadronness *)plist.FindObject("MHadronness");
571 //cout<< "Hadronness " << fH->GetHadronness() << endl;
572 }
573
574 if (!shitloop.PostProcess())
575 return -1;
576 */
577 tlist.PrintStatistics();
578
579
580 break;
581 }
582 break;
583 case kFALSE:
584 // PUT THINGS FOR THE RANDOM FOREST
585
586
587 gLog << endl;
588 gLog.Separator("Applying RF to data");
589
590
591 //
592 // Create a empty Parameter List and an empty Task List
593 // The tasklist is identified in the eventloop by its name
594 //
595 //MParList plist;
596
597 //MTaskList tlist;
598 //plist.AddToList(&tlist);
599 //
600 // ---------------------------------------------------------------
601 // ---------------------------------------------------------------
602 // first event loop: the trees of the random forest are read in
603 // ---------------------------------------------------------------
604 // ---------------------------------------------------------------
605 //
606 MReadTree readrf("Tree",kRFTreeFile);
607 readrf.DisableAutoScheme();
608
609 MRanForestFill rffill;
610 rffill.SetNumTrees(kNumTrees);
611
612 tlist.AddToList(&readrf);
613 tlist.AddToList(&rffill);
614
615 //
616 // Eventloop
617 //
618 MEvtLoop evtloop;
619 evtloop.SetParList(&plist);
620 evtloop.SetProgressBar(&bar);
621
622 if (!evtloop.Eventloop())
623 return 1;
624
625
626 // ---------------------------------------------------------------
627 //
628 // Apply RF to data
629 // ---------------------------------------------------------------
630 // ---------------------------------------------------------------
631 MTaskList tlist3;
632 plist.Replace(&tlist3);
633
634 // Calc RF
635 MRanForestCalc calc;
636
637 // Fill Hist (Only for MC file with Gammas & Hadrons mixed)
638 //MFillH fillh3a("MHHadronness");
639 //MFillH fillh3b("MHRanForest");
640
641
642 // My intercative task to apply cuts
643 //MTaskInteractive mytask;
644 //mytask.SetPreProcess(PreProcess);
645 //mytask.SetProcess(Process);
646
647 //
648 // Hadronness cut
649 //
650// MFHadAlpha cut;
651// cut.SetCutType(MFHadAlpha::kHadCut );
652// cut.SetInverted();
653// MContinue contHad(&cut);
654
655
656
657
658 // Add to TaskList
659 tlist3.AddToList(kMc ? (MTask*)&read : (MTask*)&readreal);
660 //tlist3.AddToList(&SizeCut);
661 tlist3.AddToList(&calc);
662 tlist3.AddToList(&rfs);
663 // tlist3.AddToList(&fillh3a);
664// tlist3.AddToList(&fillh3b);
665
666 //tlist.AddToList(&mytask);
667 //tlist.AddToList(&contHad);
668
669
670// tlist3.AddToList(&flist);
671 tlist3.AddToList(&ThetaCutMax);
672 tlist3.AddToList(&ThetaCutMin);
673 tlist3.AddToList(&write);
674
675
676
677
678 //
679 // Execute your analysis
680 //
681 // MProgressBar bar;
682 evtloop.SetProgressBar(&bar);
683 if (!evtloop.Eventloop())
684 return 1;
685
686 tlist3.PrintStatistics();
687
688 // plist.FindObject("MHRanForest")->DrawClone();
689// plist.FindObject("MHHadronness")->DrawClone();
690
691
692
693
694// END OF PUTTING THINGS FOR THE RANDOM FOREST
695 break;
696 }
697
698//abelardo
699 ////////////////////////////////////////////////////////////////////
700 //
701 // Second event loop: simply copy the tree MMcEvtBasic in the tree
702 // "OriginalMC" from the input file to the output file:
703 delete &write;
704 if(kMc)
705 //if(0==1)
706 {
707 MParList plist2;
708 MTaskList tlist2;
709
710 plist2.AddToList(&tlist2);
711
712 MReadMarsFile read2("OriginalMC",datafile);
713// read2.AddFile(AnalysisFilename);
714 read2.DisableAutoScheme();
715
716 tlist2.AddToList(&read2);
717
718 MWriteRootFile writeOrig(2, Form("%s{s/_I_/_Q_}", kOutpath.Data()),"UPDATE");
719 writeOrig.AddContainer("MMcEvtBasic", "OriginalMC", kFALSE);
720
721 tlist2.AddToList(&writeOrig);
722
723 MEvtLoop evtloop2;
724 bar.SetWindowName("Copying OriginalMC tree...");
725 evtloop2.SetProgressBar(&bar);
726 evtloop2.SetParList(&plist2);
727
728 if (!evtloop2.Eventloop())
729 return -1;
730
731 tlist2.PrintStatistics();
732 }
733 // end ab
734 return 0;
735}
736
737Int_t OptimizeCuts()
738{
739 return 0;
740}
741
742
743
744
745
746
747// ---------------------------------------------------------------------
748//
749// Interactive Task
750//
751//
752// Data members of the InteractiveTask
753MParList *fParList;
754MTaskList *fTaskList;
755MHadronness *fHadronness;
756MHillas *fHillas;
757
758
759Int_t PreProcess(MParList *plist)
760{
761 // Get parameter and tasklist, see Process
762 fParList = plist;
763 fTaskList = (MTaskList*)plist->FindObject("MTaskList");
764
765 // Search in the PreProcess for the objects
766 fHadronness = (MHadronness*)fParList->FindObject("MHadronness");
767 fHillas = (MHillas*)fParList->FindObject("MHillas");
768
769 return kTRUE;
770}
771
772// --------------------------------------------------------------
773//
774Int_t Process()
775{
776 Float_t size = fHillas->GetSize()/0.18; //photons!
777 Float_t had = fHadronness->GetHadronness();
778
779 if(size<500)
780 return kCONTINUE;
781
782 if(size>500 && size <1000 && had>0.32)
783 return kCONTINUE;
784
785 if(size>1000 && size <2000 && had>0.28)
786 return kCONTINUE;
787
788 if(size>2000 && size <4000 && had>0.12)
789 return kCONTINUE;
790
791 if(size>4000 && size <8000 && had>0.21)
792 return kCONTINUE;
793
794 if(size>8000 && had>0.49)
795 return kCONTINUE;
796
797 return kTRUE;
798}
799
800
801
Note: See TracBrowser for help on using the repository browser.