source: trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc@ 7562

Last change on this file since 7562 was 7552, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 14.6 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 11/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2006
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJTrainSeparation
28//
29////////////////////////////////////////////////////////////////////////////
30#include "MJTrainSeparation.h"
31
32#include <TF1.h>
33#include <TH2.h>
34#include <TChain.h>
35#include <TGraph.h>
36#include <TVirtualPad.h>
37
38#include "MHMatrix.h"
39
40#include "MLog.h"
41#include "MLogManip.h"
42
43// tools
44#include "MMath.h"
45#include "MDataSet.h"
46#include "MTFillMatrix.h"
47#include "MChisqEval.h"
48#include "MStatusDisplay.h"
49
50// eventloop
51#include "MParList.h"
52#include "MTaskList.h"
53#include "MEvtLoop.h"
54
55// tasks
56#include "MReadMarsFile.h"
57#include "MContinue.h"
58#include "MFillH.h"
59#include "MRanForestCalc.h"
60#include "MParameterCalc.h"
61
62// container
63#include "MMcEvt.hxx"
64#include "MParameters.h"
65
66// histograms
67#include "MBinning.h"
68#include "MH3.h"
69#include "MHHadronness.h"
70
71// filter
72#include "MF.h"
73#include "MFEventSelector.h"
74#include "MFilterList.h"
75
76ClassImp(MJTrainSeparation);
77
78using namespace std;
79
80void MJTrainSeparation::DisplayResult(MH3 &h31, MH3 &h32)
81{
82 TH2 &g = (TH2&)h32.GetHist();
83 TH2 &h = (TH2&)h31.GetHist();
84
85 h.SetMarkerColor(kRed);
86 g.SetMarkerColor(kGreen);
87
88 const Int_t nx = h.GetNbinsX();
89 const Int_t ny = h.GetNbinsY();
90
91 gROOT->SetSelectedPad(NULL);
92
93 TGraph gr1;
94 TGraph gr2;
95 for (int x=0; x<nx; x++)
96 {
97 TH1 *hx = h.ProjectionY("H_py", x+1, x+1);
98 TH1 *gx = g.ProjectionY("G_py", x+1, x+1);
99
100 Double_t max1 = -1;
101 Double_t max2 = -1;
102 Int_t maxy1 = 0;
103 Int_t maxy2 = 0;
104 for (int y=0; y<ny; y++)
105 {
106 const Float_t s = gx->Integral(1, y+1);
107 const Float_t b = hx->Integral(1, y+1);
108 const Float_t sig1 = MMath::SignificanceLiMa(s+b, b);
109 const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s);
110 if (sig1>max1)
111 {
112 maxy1 = y;
113 max1 = sig1;
114 }
115 if (sig2>max2)
116 {
117 maxy2 = y;
118 max2 = sig2;
119 }
120 }
121
122 gr1.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy1+1));
123 gr2.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy2+1));
124
125 delete hx;
126 delete gx;
127 }
128
129 fDisplay->AddTab("OptCut");
130 gPad->SetLogx();
131 h.DrawCopy();
132 g.DrawCopy("same");
133 gr1.SetMarkerStyle(kFullDotMedium);
134 gr1.DrawClone("LP")->SetBit(kCanDelete);
135 gr2.SetLineColor(kBlue);
136 gr2.SetMarkerStyle(kFullDotMedium);
137 gr2.DrawClone("LP")->SetBit(kCanDelete);
138}
139
140/*
141Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
142{
143 fLog->Separator("Initialize energy weighting");
144
145 if (!CheckEnv(w))
146 {
147 *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
148 return kFALSE;
149 }
150
151 TChain chain("RunHeaders");
152 set.AddFilesOn(chain);
153
154 MMcCorsikaRunHeader *h=0;
155 chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
156 chain.GetEntry(1);
157
158 if (!h)
159 {
160 *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
161 return kFALSE;
162 }
163
164 if (!w.Set(*h))
165 {
166 *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
167 return kFALSE;
168 }
169
170 w.Print();
171 return kTRUE;
172}
173
174Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
175{
176 // Some debug output
177 fLog->Separator("Compiling original MC distribution");
178
179 weight.SetNameMcEvt("MMcEvtBasic");
180 const TString w(weight.GetFormulaWeights());
181 weight.SetNameMcEvt();
182
183 *fLog << inf << "Using weights: " << w << endl;
184 *fLog << "Please stand by, this may take a while..." << flush;
185
186 if (fDisplay)
187 fDisplay->SetStatusLine1("Compiling MC distribution...");
188
189 // Create chain
190 TChain chain("OriginalMC");
191 set.AddFilesOn(chain);
192
193 // Prepare histogram
194 h.Reset();
195
196 // Fill histogram from chain
197 h.SetDirectory(gROOT);
198 if (h.InheritsFrom(TH2::Class()))
199 {
200 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
201 h.SetXTitle("\\Theta [\\circ]");
202 h.SetYTitle("E [GeV]");
203 h.SetZTitle("Counts");
204 chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
205 }
206 else
207 {
208 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
209 h.SetXTitle("\\Theta [\\circ]");
210 h.SetYTitle("Counts");
211 chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
212 }
213 h.SetDirectory(0);
214
215 *fLog << "done." << endl;
216 if (fDisplay)
217 fDisplay->SetStatusLine2("done.");
218
219 if (h.GetEntries()>0)
220 return kTRUE;
221
222 *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
223
224 return h.GetEntries()>0;
225}
226*/
227
228Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
229{
230 TChain chain("OriginalMC");
231 set.AddFilesOn(chain);
232
233 min = chain.GetMinimum("MMcEvtBasic.fEnergy");
234 max = chain.GetMaximum("MMcEvtBasic.fEnergy");
235
236 num = chain.GetEntries();
237
238 if (num<100)
239 *fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;
240
241 return num>=100;
242}
243
244Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const
245{
246 TChain chain1("Events");
247 set.AddFilesOff(chain1);
248
249 const Double_t num = chain1.GetEntries();
250 if (num<100)
251 {
252 *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
253 return -1;
254 }
255
256 TChain chain("EffectiveOnTime");
257 set.AddFilesOff(chain);
258
259 TH1F h;
260 h.SetDirectory(gROOT);
261 h.SetNameTitle("OnTime", "Effective on-time");
262 chain.Draw("MEffectiveOnTime.fVal>>OnTime", "", "goff");
263 h.SetDirectory(0);
264
265 if (h.Integral()<1)
266 {
267 *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
268 return -1;
269 }
270
271 *fLog << inf << "Found " << num << " events in " << h.Integral();
272 *fLog << "s (" << num/h.Integral() << "Hz)" << endl;
273
274 return num/h.Integral();
275}
276
277/*
278 Scale:
279
280
281 TF1 fold("old", "x^(-2.6)", emin, emax);
282 TF1 fnew("new", "x^(-4.0)", emin, emax);
283
284 TF1 q("q", "new/old", emin, emax);
285
286 Double_t scale = 1./q.GetMaximum(emin, emax);
287
288 // Anzahl produzierter Events vor MFEnergySlope:
289 Double_t nold = fold.Integral(emin, emax);
290
291 // Anzahl produzierter Events nach MFEnergySlope:
292 Double_t nnew = fnew.Integral(emin, emax)*scale;
293
294 class MFSpectrum : MMcSpectrumWeight
295 {
296 Double_t fScale;
297 Bool_t fResult;
298
299 MFSpectrum::MFSpectrum(const char *name, const char *title)
300 {
301 fName = name ? name : "MMcSpectrumWeight";
302 fTitle = title ? title : "Task to calculate weights to change the energy spectrum";
303
304 Init(fName, fTitle);
305
306 }
307
308 Int_t PreProcess(MParList *pList)
309 {
310 Int_t rc = MFSpectrumWeight::PreProcess(pList);
311 if (rc!=kTRUE)
312 return rc;
313
314 fScale = fEval->GetMaximum(fEnergyMin, fEnergyMax);
315
316 return kTRUE;
317 }
318
319 Int_t Process()
320 {
321 const Double_t e = fMcEvt->GetEnergy();
322
323 Double_t prob = fFunc->Eval(e)/fScale;
324
325 const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope);
326 const Float_t Nrnd = ;
327
328 fResult = Nexp >= gRandom->Uniform();
329 }
330
331 }
332
333
334 */
335
336Bool_t MJTrainSeparation::AutoTrain()
337{
338 Double_t num, min, max;
339 if (!GetEventsProduced(fDataSetTrain, num, min, max))
340 return kFALSE;
341
342 *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;
343
344 // Target spectrum
345 TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max);
346 flx.SetParameter(0, 1e-5);
347
348 // Number n0 of events this spectrum would produce per s and m^2
349 const Double_t n0 = flx.Integral(min, max); //[#]
350
351 // Area produced in MC
352 const Double_t A = TMath::Pi()*300*300; //[m²]
353
354 // Rate R of events this spectrum would produce per s
355 const Double_t R = n0*A; //[Hz]
356
357 // Number N of events produced (in trainings sample)
358 const Double_t N = num; //[#]
359
360 // This correponds to an observation time T [s]
361 const Double_t T = N/R; //[s]
362
363 // With an average data rate after star of
364 const Double_t r = GetDataRate(fDataSetTrain); //[Hz]
365
366 // this yields a number of n events to be read for training
367 const Double_t n = r*T; //[#]
368
369 if (r<0)
370 return kFALSE;
371
372 *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl;
373 *fLog << "For a data rate of " << r << "Hz this corresponds to " << n << " data events." << endl;
374
375 fNumTrainOn = (UInt_t)-1;
376 fNumTrainOff = TMath::Nint(n);
377
378 /*
379 An event rate dependent selection?
380 ----------------------------------
381 Total average data rate: R
382 Goal number of events: N
383 Number of data events: N0
384 Rate assigned to single evt: r
385
386 Selection probability: N/N0 * r/R
387
388 f := N/N0 * r
389
390 MF f("f * MEventRate.fRate < rand");
391 */
392
393
394 return kTRUE;
395}
396
397Bool_t MJTrainSeparation::Train(const char *out)
398{
399 if (!fDataSetTrain.IsValid())
400 {
401 *fLog << err << "ERROR - DataSet for training invalid!" << endl;
402 return kFALSE;
403 }
404 if (!fDataSetTest.IsValid())
405 {
406 *fLog << err << "ERROR - DataSet for testing invalid!" << endl;
407 return kFALSE;
408 }
409
410 // ----------------------- Auto Train? ----------------------
411
412 if (fAutoTrain)
413 if (!AutoTrain())
414 return kFALSE;
415
416 // --------------------- Setup files --------------------
417 MReadMarsFile read1("Events");
418 MReadMarsFile read2("Events");
419 MReadMarsFile read3("Events");
420 MReadMarsFile read4("Events");
421 read1.DisableAutoScheme();
422 read2.DisableAutoScheme();
423 read3.DisableAutoScheme();
424 read4.DisableAutoScheme();
425
426 fDataSetTrain.AddFilesOn(read1);
427 fDataSetTrain.AddFilesOff(read3);
428
429 fDataSetTest.AddFilesOff(read2);
430 fDataSetTest.AddFilesOn(read4);
431
432 // ----------------------- Setup RF ----------------------
433 MHMatrix train("Train");
434 train.AddColumns(fRules);
435 train.AddColumn("MHadronness.fVal");
436
437 // ----------------------- Fill Matrix RF ----------------------
438
439 MParameterD had("MHadronness");
440
441 MParList plistx;
442 plistx.AddToList(&had);
443 plistx.AddToList(this);
444
445 MTFillMatrix fill;
446 fill.SetLogStream(fLog);
447 fill.SetDisplay(fDisplay);
448 fill.AddPreCuts(fPreCuts);
449 fill.AddPreCuts(fTrainCuts);
450
451 // Set classifier for gammas
452 had.SetVal(0);
453 fill.SetName("FillGammas");
454 fill.SetDestMatrix1(&train, fNumTrainOn);
455 fill.SetReader(&read1);
456 if (!fill.Process(plistx))
457 return kFALSE;
458
459 // Set classifier for hadrons
460 had.SetVal(1);
461 fill.SetName("FillBackground");
462 fill.SetDestMatrix1(&train, fNumTrainOff);
463 fill.SetReader(&read3);
464 if (!fill.Process(plistx))
465 return kFALSE;
466
467 // ------------------------ Train RF --------------------------
468
469 MRanForestCalc rf;
470 rf.SetNumTrees(fNumTrees);
471 rf.SetNdSize(fNdSize);
472 rf.SetNumTry(fNumTry);
473 rf.SetNumObsoleteVariables(1);
474 rf.SetDebug(fDebug);
475 rf.SetDisplay(fDisplay);
476 rf.SetLogStream(fLog);
477 rf.SetFileName(out);
478 rf.SetNameOutput("MHadronness");
479
480 //MBinning b(2, -0.5, 1.5, "BinningHadronness", "lin");
481
482 //if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification
483 // return;
484
485 //if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification
486 // return;
487
488 if (!rf.TrainSingleRF(train)) // regression
489 return kFALSE;
490
491 //fDisplay = rf.GetDisplay();
492
493 // --------------------- Display result ----------------------
494 gLog.Separator("Test");
495
496 MParList plist;
497 MTaskList tlist;
498 plist.AddToList(this);
499 plist.AddToList(&tlist);
500
501 MMcEvt mcevt;
502 plist.AddToList(&mcevt);
503
504 // ----- Setup histograms -----
505 MBinning binsy(100, 0 , 1, "BinningMH3Y", "lin");
506 MBinning binsx( 50, 10, 100000, "BinningMH3X", "log");
507
508 plist.AddToList(&binsx);
509 plist.AddToList(&binsy);
510
511 MH3 h31("MHillas.fSize", "MHadronness.fVal");
512 MH3 h32("MHillas.fSize", "MHadronness.fVal");
513 h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness");
514 h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness");
515
516 MHHadronness hist;
517
518 // ----- Setup tasks -----
519 MFillH fillh0(&hist, "", "FillHadronness");
520 MFillH fillh1(&h31);
521 MFillH fillh2(&h32);
522 fillh1.SetNameTab("Background");
523 fillh2.SetNameTab("Gammas");
524 fillh0.SetBit(MFillH::kDoNotDisplay);
525
526 // ----- Setup filter -----
527 MFilterList precuts;
528 precuts.AddToList(fPreCuts);
529 precuts.AddToList(fTestCuts);
530
531 MContinue c0(&precuts);
532 c0.SetName("PreCuts");
533 c0.SetInverted();
534
535 MFEventSelector sel;
536 sel.SetNumSelectEvts(fNumTestOn);
537
538 MContinue c1(&sel);
539 c1.SetInverted();
540
541 // ----- Setup tasklist -----
542 tlist.AddToList(&read2);
543 tlist.AddToList(&c0);
544 tlist.AddToList(&c1);
545 tlist.AddToList(&rf);
546 tlist.AddToList(&fillh0);
547 tlist.AddToList(&fillh1);
548
549 // ----- Run eventloop on gammas -----
550 MEvtLoop loop;
551 loop.SetDisplay(fDisplay);
552 loop.SetLogStream(fLog);
553 loop.SetParList(&plist);
554
555 if (!loop.Eventloop())
556 return kFALSE;
557
558 // ----- Setup and run eventloop on background -----
559 sel.SetNumSelectEvts(fNumTestOff);
560 fillh0.ResetBit(MFillH::kDoNotDisplay);
561
562 tlist.Replace(&read4);
563 tlist.Replace(&fillh2);
564
565 if (!loop.Eventloop())
566 return kFALSE;
567
568 DisplayResult(h31, h32);
569
570 if (!WriteDisplay(out))
571 return kFALSE;
572
573 return kTRUE;
574}
575
Note: See TracBrowser for help on using the repository browser.