| 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 Hengstebeck 3/2003 <mailto:hengsteb@alwa02.physik.uni-siegen.de>!
|
|---|
| 19 | !
|
|---|
| 20 | ! Copyright: MAGIC Software Development, 2000-2003
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 |
|
|---|
| 25 | #include "MMcEvt.hxx"
|
|---|
| 26 |
|
|---|
| 27 | void RanForestPD(Int_t minsize=0)
|
|---|
| 28 | {
|
|---|
| 29 | //
|
|---|
| 30 | // Create a empty Parameter List and an empty Task List
|
|---|
| 31 | // The tasklist is identified in the eventloop by its name
|
|---|
| 32 | //
|
|---|
| 33 | MParList plist;
|
|---|
| 34 |
|
|---|
| 35 | MTaskList tlist;
|
|---|
| 36 | plist.AddToList(&tlist);
|
|---|
| 37 |
|
|---|
| 38 | MReadMarsFile read("Events", "star_gamma_train.root");
|
|---|
| 39 | Float_t numgammas = read.GetEntries();
|
|---|
| 40 |
|
|---|
| 41 | read.AddFile("star_proton_train.root");
|
|---|
| 42 | read.AddFile("star_helium_train.root");
|
|---|
| 43 |
|
|---|
| 44 | Float_t numhadrons = read.GetEntries() - numgammas;
|
|---|
| 45 |
|
|---|
| 46 | // Fraction of gammas to be used in training:
|
|---|
| 47 | // (less than total sample, to speed up execution)
|
|---|
| 48 | //
|
|---|
| 49 | Float_t gamma_frac = 1.5*(numhadrons / numgammas);
|
|---|
| 50 |
|
|---|
| 51 |
|
|---|
| 52 | read.DisableAutoScheme();
|
|---|
| 53 |
|
|---|
| 54 | tlist.AddToList(&read);
|
|---|
| 55 |
|
|---|
| 56 | MFParticleId fgamma("MMcEvt", '=', MMcEvt::kGAMMA);
|
|---|
| 57 | tlist.AddToList(&fgamma);
|
|---|
| 58 |
|
|---|
| 59 | MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA);
|
|---|
| 60 | tlist.AddToList(&fhadrons);
|
|---|
| 61 |
|
|---|
| 62 |
|
|---|
| 63 | MHMatrix matrixg("MatrixGammas");
|
|---|
| 64 |
|
|---|
| 65 | // matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))");
|
|---|
| 66 |
|
|---|
| 67 | // matrixg.AddColumn("MMcEvt.fTelescopePhi");
|
|---|
| 68 | // matrixg.AddColumn("MSigmabar.fSigmabar");
|
|---|
| 69 |
|
|---|
| 70 | matrixg.AddColumn("MHillas.fSize");
|
|---|
| 71 | matrixg.AddColumn("MHillasSrc.fDist");
|
|---|
| 72 |
|
|---|
| 73 | // matrixg.AddColumn("MHillas.fWidth");
|
|---|
| 74 | // matrixg.AddColumn("MHillas.fLength");
|
|---|
| 75 |
|
|---|
| 76 | // Scaled Width and Length:
|
|---|
| 77 | //
|
|---|
| 78 | matrixg.AddColumn("MHillas.fWidth/(-13.1618+10.0492*log10(MHillas.fSize))");
|
|---|
| 79 | matrixg.AddColumn("MHillas.fLength/(-57.3784+32.6131*log10(MHillas.fSize))");
|
|---|
| 80 |
|
|---|
| 81 |
|
|---|
| 82 | // matrixg.AddColumn("(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
|
|---|
| 83 | matrixg.AddColumn("MNewImagePar.fConc");
|
|---|
| 84 |
|
|---|
| 85 | // matrixg.AddColumn("MNewImagePar.fConc1");
|
|---|
| 86 | // matrixg.AddColumn("MNewImagePar.fAsym");
|
|---|
| 87 | // matrixg.AddColumn("MHillasExt.fM3Trans");
|
|---|
| 88 | // matrixg.AddColumn("abs(MHillasSrc.fHeadTail)");
|
|---|
| 89 | // matrixg.AddColumn("MNewImagePar.fLeakage1");
|
|---|
| 90 |
|
|---|
| 91 | //
|
|---|
| 92 | // Third moment, with sign referred to source position:
|
|---|
| 93 | //
|
|---|
| 94 | matrixg.AddColumn("MHillasExt.fM3Long*MHillasSrc.fCosDeltaAlpha/abs(MHillasSrc.fCosDeltaAlpha)");
|
|---|
| 95 | matrixg.AddColumn("MHillasSrc.fAlpha");
|
|---|
| 96 |
|
|---|
| 97 |
|
|---|
| 98 | //
|
|---|
| 99 | // SIZE cut:
|
|---|
| 100 | //
|
|---|
| 101 | TString dummy("(MHillas.fSize<");
|
|---|
| 102 | dummy += minsize;
|
|---|
| 103 |
|
|---|
| 104 | // AM, useful FOR High Energy only:
|
|---|
| 105 | // dummy += ") || (MImagePar.fNumSatPixelsLG>0";
|
|---|
| 106 | dummy += ")";
|
|---|
| 107 |
|
|---|
| 108 | MContinue SizeCut(dummy);
|
|---|
| 109 | tlist.AddToList(&SizeCut);
|
|---|
| 110 |
|
|---|
| 111 | plist.AddToList(&matrixg);
|
|---|
| 112 |
|
|---|
| 113 | MHMatrix matrixh("MatrixHadrons");
|
|---|
| 114 | matrixh.AddColumns(matrixg.GetColumns());
|
|---|
| 115 | plist.AddToList(&matrixh);
|
|---|
| 116 |
|
|---|
| 117 | MFillH fillmath("MatrixHadrons");
|
|---|
| 118 | fillmath.SetFilter(&fhadrons);
|
|---|
| 119 | tlist.AddToList(&fillmath);
|
|---|
| 120 |
|
|---|
| 121 |
|
|---|
| 122 | MFEventSelector reduce_training_gammas;
|
|---|
| 123 | reduce_training_gammas.SetSelectionRatio(gamma_frac);
|
|---|
| 124 | tlist.AddToList(&reduce_training_gammas);
|
|---|
| 125 |
|
|---|
| 126 | MFillH fillmatg("MatrixGammas");
|
|---|
| 127 | fillmatg.SetFilter(&fgamma);
|
|---|
| 128 | fillmatg.SetFilter(&reduce_training_gammas);
|
|---|
| 129 | tlist.AddToList(&fillmatg);
|
|---|
| 130 |
|
|---|
| 131 | //
|
|---|
| 132 | // Create and setup the eventloop
|
|---|
| 133 | //
|
|---|
| 134 | MEvtLoop evtloop;
|
|---|
| 135 | evtloop.SetParList(&plist);
|
|---|
| 136 |
|
|---|
| 137 | MProgressBar bar;
|
|---|
| 138 | evtloop.SetProgressBar(&bar);
|
|---|
| 139 | //
|
|---|
| 140 | // Execute your analysis
|
|---|
| 141 | //
|
|---|
| 142 | if (!evtloop.Eventloop())
|
|---|
| 143 | return;
|
|---|
| 144 |
|
|---|
| 145 | tlist.PrintStatistics();
|
|---|
| 146 |
|
|---|
| 147 | // ---------------------------------------------------------------
|
|---|
| 148 | // ---------------------------------------------------------------
|
|---|
| 149 | // second event loop: the trees of the random forest are grown,
|
|---|
| 150 | // the event loop is now actually a tree loop (loop of training
|
|---|
| 151 | // process)
|
|---|
| 152 | // ---------------------------------------------------------------
|
|---|
| 153 | // ---------------------------------------------------------------
|
|---|
| 154 |
|
|---|
| 155 | MTaskList tlist2;
|
|---|
| 156 | plist.Replace(&tlist2);
|
|---|
| 157 |
|
|---|
| 158 | //
|
|---|
| 159 | // Just a check:
|
|---|
| 160 | // matrixh.ShuffleRows(9999);
|
|---|
| 161 | // matrixg.ShuffleRows(1111);
|
|---|
| 162 | //
|
|---|
| 163 |
|
|---|
| 164 | MRanForestGrow rfgrow2;
|
|---|
| 165 | rfgrow2.SetNumTrees(100);
|
|---|
| 166 | rfgrow2.SetNumTry(3);
|
|---|
| 167 | rfgrow2.SetNdSize(10);
|
|---|
| 168 |
|
|---|
| 169 | MWriteRootFile rfwrite2("RF.root");
|
|---|
| 170 | rfwrite2.AddContainer("MRanTree","Tree"); //save all trees
|
|---|
| 171 | MFillH fillh2("MHRanForestGini");
|
|---|
| 172 |
|
|---|
| 173 | tlist2.AddToList(&rfgrow2);
|
|---|
| 174 | tlist2.AddToList(&rfwrite2);
|
|---|
| 175 | tlist2.AddToList(&fillh2);
|
|---|
| 176 |
|
|---|
| 177 | // gRandom is accessed from MRanForest (-> bootstrap aggregating)
|
|---|
| 178 | // and MRanTree (-> random split selection) and should be initialized
|
|---|
| 179 | // here if you want to set a certain random number generator
|
|---|
| 180 | if(gRandom)
|
|---|
| 181 | delete gRandom;
|
|---|
| 182 | // gRandom = new TRandom3(0); // Takes seed from the computer clock
|
|---|
| 183 | gRandom = new TRandom3(); // Uses always same seed
|
|---|
| 184 | //
|
|---|
| 185 | // Execute tree growing (now the eventloop is actually a treeloop)
|
|---|
| 186 | //
|
|---|
| 187 |
|
|---|
| 188 | evtloop.SetProgressBar(&bar);
|
|---|
| 189 | if (!evtloop.Eventloop())
|
|---|
| 190 | return;
|
|---|
| 191 |
|
|---|
| 192 | tlist2.PrintStatistics();
|
|---|
| 193 |
|
|---|
| 194 | plist.FindObject("MHRanForestGini")->DrawClone();
|
|---|
| 195 |
|
|---|
| 196 | // ---------------------------------------------------------------
|
|---|
| 197 | // ---------------------------------------------------------------
|
|---|
| 198 | // third event loop: the control sample (star2.root) is processed
|
|---|
| 199 | // through the previously grown random forest,
|
|---|
| 200 | //
|
|---|
| 201 | // the histograms MHHadronness (quality of g/h-separation) and
|
|---|
| 202 | // MHRanForest are created and displayed.
|
|---|
| 203 | // MHRanForest shows the convergence of the classification error
|
|---|
| 204 | // as function of the number of grown (and combined) trees
|
|---|
| 205 | // and tells the user how many trees are actually needed in later
|
|---|
| 206 | // classification tasks.
|
|---|
| 207 | // ---------------------------------------------------------------
|
|---|
| 208 | // ---------------------------------------------------------------
|
|---|
| 209 |
|
|---|
| 210 | MTaskList tlist3;
|
|---|
| 211 |
|
|---|
| 212 | plist.Replace(&tlist3);
|
|---|
| 213 |
|
|---|
| 214 | MReadMarsFile read3("Events", "star_gamma_test.root");
|
|---|
| 215 |
|
|---|
| 216 | read3.AddFile("star_proton_test.root");
|
|---|
| 217 | read3.AddFile("star_helium_test.root");
|
|---|
| 218 |
|
|---|
| 219 |
|
|---|
| 220 | /*
|
|---|
| 221 | MReadMarsFile read3("Events", "star_gammas_test_extended.root");
|
|---|
| 222 | read3.AddFile("star_protons_test_extended.root");
|
|---|
| 223 | read3.AddFile("star_helium_test_extended.root");
|
|---|
| 224 | */
|
|---|
| 225 |
|
|---|
| 226 | read3.DisableAutoScheme();
|
|---|
| 227 | tlist3.AddToList(&read3);
|
|---|
| 228 | tlist3.AddToList(&SizeCut);
|
|---|
| 229 |
|
|---|
| 230 | MRanForestCalc calc;
|
|---|
| 231 | tlist3.AddToList(&calc);
|
|---|
| 232 |
|
|---|
| 233 | // parameter containers you want to save
|
|---|
| 234 | //
|
|---|
| 235 |
|
|---|
| 236 | TString outfname = "star_S";
|
|---|
| 237 | outfname += minsize;
|
|---|
| 238 | outfname += "_all.root";
|
|---|
| 239 | MWriteRootFile write(outfname, "recreate");
|
|---|
| 240 |
|
|---|
| 241 |
|
|---|
| 242 | // MWriteRootFile write("star_surviving_hadrons.root", "recreate");
|
|---|
| 243 | // MWriteRootFile write("star_surviving_gammas.root", "recreate");
|
|---|
| 244 |
|
|---|
| 245 | // write.AddContainer("MSigmabar", "Events");
|
|---|
| 246 |
|
|---|
| 247 | write.AddContainer("MMcEvt", "Events");
|
|---|
| 248 | write.AddContainer("MHillas", "Events");
|
|---|
| 249 | write.AddContainer("MHillasExt", "Events");
|
|---|
| 250 | write.AddContainer("MImagePar", "Events");
|
|---|
| 251 | write.AddContainer("MNewImagePar", "Events");
|
|---|
| 252 | write.AddContainer("MHillasSrc", "Events");
|
|---|
| 253 | write.AddContainer("MHadronness", "Events");
|
|---|
| 254 |
|
|---|
| 255 | write.AddContainer("MRawRunHeader", "RunHeaders");
|
|---|
| 256 | write.AddContainer("MSrcPosCam", "RunHeaders");
|
|---|
| 257 | write.AddContainer("MMcRunHeader", "RunHeaders");
|
|---|
| 258 |
|
|---|
| 259 | /*
|
|---|
| 260 | write.AddContainer("MMcTrig", "Events");
|
|---|
| 261 | write.AddContainer("MRawEvtData", "Events");
|
|---|
| 262 | write.AddContainer("MRawEvtHeader", "Events");
|
|---|
| 263 | */
|
|---|
| 264 |
|
|---|
| 265 |
|
|---|
| 266 | MFillH fillh3a("MHHadronness");
|
|---|
| 267 | MFillH fillh3b("MHRanForest");
|
|---|
| 268 |
|
|---|
| 269 | tlist3.AddToList(&fillh3a);
|
|---|
| 270 | tlist3.AddToList(&fillh3b);
|
|---|
| 271 |
|
|---|
| 272 |
|
|---|
| 273 | tlist3.AddToList(&write);
|
|---|
| 274 |
|
|---|
| 275 | evtloop.SetProgressBar(&bar);
|
|---|
| 276 |
|
|---|
| 277 | //
|
|---|
| 278 | // Execute your analysis
|
|---|
| 279 | //
|
|---|
| 280 | if (!evtloop.Eventloop())
|
|---|
| 281 | return;
|
|---|
| 282 |
|
|---|
| 283 | tlist3.PrintStatistics();
|
|---|
| 284 |
|
|---|
| 285 | plist.FindObject("MHRanForest")->DrawClone();
|
|---|
| 286 |
|
|---|
| 287 | TCanvas *had = new TCanvas;
|
|---|
| 288 | had->cd();
|
|---|
| 289 | plist.FindObject("MHHadronness")->Draw();
|
|---|
| 290 |
|
|---|
| 291 | TString gifname = "had_S";
|
|---|
| 292 | gifname += minsize;
|
|---|
| 293 | gifname += ".gif";
|
|---|
| 294 |
|
|---|
| 295 | had->SaveAs(gifname);
|
|---|
| 296 | had->Close();
|
|---|
| 297 |
|
|---|
| 298 | plist.FindObject("MHHadronness")->DrawClone();
|
|---|
| 299 | plist.FindObject("MHHadronness")->Print();//*/
|
|---|
| 300 |
|
|---|
| 301 | return;
|
|---|
| 302 | }
|
|---|
