source: trunk/Mars/macros/train/traindisp.C@ 13389

Last change on this file since 13389 was 8729, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 7.2 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, 2000-2007
21!
22!
23\* ======================================================================== *//////////////////////////////////////////////////////////////////////////////
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// traindisp.C
28// ===========
29//
30// Trains a random forest for disp estimation.
31//
32// The additional code below shows how the MagicCuts can be used in
33// training and testing. There is also code which allows to change the
34// slope of the input spectrum. If a min- or max-break energy is given
35// the new slope is only applied in this region. Further possibilities
36// for additional cuts are shown.
37//
38// SequencesOn: Monte Carlo Sequences used for training
39// SequencesOff: Monte Carlo Sequences used for testing
40//
41// Use:
42// opt.AddPreCut to use cut for test and training
43// opt.AddTestCut to use cut for test only
44// opt.AddTrainCut to use cut for train only
45// opt.AddPreTask a task executed before filling the matrix/evaluating the RF
46// opt.AddPostTask a task executed before training, after evaluating the RF
47//
48/////////////////////////////////////////////////////////////////////////////
49
50void traindisp()
51{
52 MDataSet set("mcsponde/mcdataset-for-training.txt");
53 set.SetNumAnalysis(1); // Necessary
54
55 // alternatively use:
56 //MDataSet set("mcctesttrain.txt", "/magic/montacrlo/sequences", "/magic/montecarlo/star");
57
58 MJTrainDisp opt;
59 //opt.SetDebug();
60
61 // These are control parameters how to train the random forest (use with care!)
62 //opt.SetNumTrees(100);
63 //opt.SetNdSize(5);
64 //opt.SetNumTry(0);
65
66 // --- Change the theta-cut for the displayed cut-effieincy ---
67 //opt.SetThetaCut(0.16);
68
69 // ------- Parameters to train Random Forest --------
70 opt.AddParameter("MHillas.fLength");
71 opt.AddParameter("MHillas.fWidth");
72 opt.AddParameter("MHillas.fSize");
73 opt.AddParameter("MNewImagePar.fLeakage1");
74 opt.AddParameter("MPointingPos.fZd");
75 opt.AddParameter("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)");
76 //opt.AddParameter("MNewImagePar.fLeakage2");
77 //opt.AddParameter("MNewImagePar.fConc");
78 //opt.AddParameter("MNewImagePar.fConc1");
79 //opt.AddParameter("MNewImagePar.fUsedArea");
80 //opt.AddParameter("MNewImagePar.fCoreArea");
81 //opt.AddParameter("MNewImagePar.fNumUsedPixels");
82 //opt.AddParameter("MNewImagePar.fNumCorePixels");
83 //opt.AddParameter("MImagePar.fSizeSinglePixels");
84 //opt.AddParameter("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)");
85 //opt.AddParameter("MHillasExt.fAsym*sign(MHillasSrc.fCosDeltaAlpha)");
86
87 // -------------------- Run ----------------------------
88
89 MStatusDisplay *d = new MStatusDisplay;
90 opt.SetDisplay(d);
91
92 /*
93 // -------------------- Magic-Cuts ----------------------
94 MFMagicCuts cuts;
95 cuts.SetHadronnessCut(MFMagicCuts::kArea);
96 cuts.SetThetaCut(MFMagicCuts::kNone);
97
98 TArrayD arr(13);
99 arr[0] = 1.15136;
100 arr[1] = 0.215;
101 arr[2] = 0.215468;
102 arr[3] = 5.63973;
103 arr[4] = 0.0836169;
104 arr[5] = -0.07;
105 arr[6] = 7.2;
106 arr[7] = 0.5;
107 arr[8] = 0.0681437;
108 arr[9] = 2.62932;
109 arr[10] = 1.51279;
110 arr[11] = 0.0507821;
111
112 cuts.SetVariables(arr);
113
114 opt.AddTestCut(&cuts);
115
116 // -------------------- Other cuts ----------------------
117 opt.AddPreCut("MNewImagePar.fLeakage1<0.0001");
118 opt.AddPreCut("MHillas.fSize>200");
119 opt.AddPreCut("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg<0.3");
120 opt.AddPreCut("MPointingPos.fZd<25");
121 opt.AddTestCut("MHillasExt.fM3Long *sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg>-0.07");
122 opt.AddTestCut("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg<(MHillasSrc.fDist*MGeomCam.fConvMm2Deg-0.5)*7.2");
123
124 // ---------------------- Histogram --------------------
125 MHn hist("MyHist", "Dist Residual (Disp-Dist)");
126
127 MBinning bins(50, -0.5, 0.5);
128
129 hist.AddHist("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*3.37e-3", "Disp.fVal-MHillasSrc.fDist*3.37e-3");
130 hist.InitName("ResM3l;M3Long;ResidualDist");
131 hist.InitTitle(";M3l [\\circ];Disp-Dist [\\circ];");
132 hist.SetDrawOption("colz profx");
133 hist.SetBinnings(&bins);
134
135 hist.AddHist("MHillasExt.fAsym*sign(MHillasSrc.fCosDeltaAlpha)*3.37e-3", "Disp.fVal-MHillasSrc.fDist*3.37e-3");
136 hist.InitName("ResAsym;M3Long;ResidualDist");
137 hist.InitTitle(";Asym [\\circ];Disp-Dist [\\circ];");
138 hist.SetDrawOption("colz profx");
139 hist.SetBinnings(&bins);
140
141 MFillH fill(&hist, "", "FillMyHist");
142 //fill.SetWeight(); // Enable weights to be used to fill this histogram
143 opt.AddTestTask(&fill);
144
145 // -------------------- Energy Slope --------------------
146 // Note, that weight normally doesn't improve anything here.
147 // This is a way to throw away events to a different slope
148 MFEnergySlope slope(-4.0); // New slope for mc spectrum
149 opt.AddPreCut(&slope); // throw away events to change slope
150
151 // This is a way to weight the events to a different spectrum
152 MMcSpectrumWeight weight;
153 weight.SetFormula("pow(X/300, -2.31-0.26*log10(X/300))");
154 opt.SetWeights(&weight);
155
156 // ------------------ Zd distribution -------------------
157 TFile file("ganymed00001111.root");
158
159 MStatusArray arr;
160 if (arr.Read()<=0)
161 return;
162 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
163 if (!vstime)
164 return -1;
165
166 MMcSpectrumWeight weight;
167 weight.SetWeightsZd(vstime);
168 opt.AddPreTask(&weight);
169
170 // ------- A trick for testing the parametrization -------
171 // This will overwrite the random forest calculated disp.
172 // Therefore you can test our disp-parametrization.
173 MParameterCalc calc("(1.15136 +"
174 "0.0681437*MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3 +"
175 "0.0507821*(log10(MHillas.fSize)-1.51279)^2*(log10(MHillas.fSize)>1.51279) +"
176 "2.62932*MNewImagePar.fLeakage1)*"
177 "(1-MHillas.fWidth/MHillas.fLength)");
178 calc.SetNameParameter("Disp");
179 opt.AddPostTask(&calc);
180
181 // ------------------------------------------------------
182 */
183
184 // To allow overwrite of existing files
185 // opt.SetOverwrite();
186
187 // The number is the number of events read from the file
188 // randomly and is the number of events before cuts
189 opt.Train("rf-disp.root", set, 30000);
190}
Note: See TracBrowser for help on using the repository browser.