- Timestamp:
- 06/03/05 21:37:54 (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestDISP.C
r6028 r7137 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): Abelardo Moralejo 1/2005 <mailto:moralejo@pd.infn.it>! 19 ! 20 ! Copyright: MAGIC Software Development, 2000-2003 21 ! 22 ! 23 \* ======================================================================== */ 24 25 ////////////////////////////////////////////////////////////////////////////////// 26 // 27 // macro RanForestDISP.C 28 // 29 // Version of Random Forest macro, intended to run on real data without assuming 30 // any particular source position in the camera FOV. So we do not use the DIST 31 // parameter, but rather (one version of) DISP method. We do not use either the 32 // sign of the third longitudinal moment with respect to the center of the camera 33 // (since the source can be somewhere else like in the case of wobble mode). 34 // 35 ////////////////////////////////////////////////////////////////////////////////// 36 37 38 #include "MMcEvt.hxx" 39 40 void RanForestDISP(Int_t minsize=1000) 41 { 42 Int_t mincorepixels = 6; // minimum number of core pixels. 43 // Int_t mincorepixels = 0; 44 45 TString skipevents; 46 // skipevents = "(MMcEvt.fImpact>15000.)"; 47 // skipevents = "(MNewImagePar.fLeakage1>0.001)"; 48 // skipevents = "(MHillas.fSize>1000)"; 49 50 // Skip calibration events, leaking events and flash events: 51 skipevents = "(MNewImagePar.fLeakage1>0.1) || ({1.5+log10(MNewImagePar.fConc)*2.33333}-{6.5-log10(MHillas.fSize)}>0.)"; 52 53 // 54 // Create a empty Parameter List and an empty Task List 55 // The tasklist is identified in the eventloop by its name 56 // 57 MParList plist; 58 59 MTaskList tlist; 60 plist.AddToList(&tlist); 61 62 63 MReadMarsFile read("Events", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Calib_6slices/tc_3.84_2.56/star_gamma_train_new.root"); 64 Float_t numgammas = read.GetEntries(); 65 66 read.AddFile("star_train_20050110_CrabNebulaW1.root"); 67 Float_t numhadrons = read.GetEntries() - numgammas; 68 69 // Fraction of gammas to be used in training: 70 // Float_t gamma_frac = 1.5*(numhadrons / numgammas); 71 Float_t gamma_frac = 1.; 72 73 // Total hadron number to be processed in training: 74 // If you run RF over different data runs, using always the same # of training gammas, 75 // you should also use always a fixed number of hadrons for training (although the data 76 // runs may be of different duration). Otherwise the meaning of the output hadronness 77 // would be different for the different subsamples, since its absolute values depend on 78 // the statistics of the training sample. The number set here is the number BEFORE the 79 // "a priori cuts" set above through "skipevents". 80 81 Float_t select_n_hadrons = 37000.; 82 Float_t hadron_frac = select_n_hadrons / numhadrons; 83 84 read.DisableAutoScheme(); 85 86 tlist.AddToList(&read); 87 88 MFParticleId fgamma("MMcEvt", '=', MMcEvt::kGAMMA); 89 tlist.AddToList(&fgamma); 90 91 MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA); 92 tlist.AddToList(&fhadrons); 93 94 95 MHMatrix matrixg("MatrixGammas"); 96 97 // TEST TEST TEST 98 // matrixg.AddColumn("MMcEvt.fImpact"); 99 // matrixg.AddColumn("MMcEvt.fMuonCphFraction"); 100 101 // matrixg.AddColumn("MNewImagePar.fLeakage1"); 102 // matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))"); 103 104 // matrixg.AddColumn("MMcEvt.fTelescopePhi"); 105 // matrixg.AddColumn("MSigmabar.fSigmabar"); 106 107 108 matrixg.AddColumn("MHillas.fSize"); 109 110 // We do not use DIST, since it depends on knowing the source position on camera. 111 // In wobble mode the position changes with time. Besides, we intend to do a 2-d map 112 // with the DISP method, and so we cannot use DIST in the cuts (since it would amde the 113 // camera center a "provoleged" position). 114 115 matrixg.AddColumn("(1.-(MHillas.fWidth/MHillas.fLength))*(0.9776+(0.101062*log10(MHillas.fSize)))"); 116 117 // matrixg.AddColumn("MHillas.fWidth"); 118 // matrixg.AddColumn("MHillas.fLength"); 119 120 121 // Scaling from real data? 122 // matrixg.AddColumn("MHillas.fWidth/(-8.37847+(9.58826*log10(MHillas.fSize)))"); 123 // matrixg.AddColumn("MHillas.fLength/(-40.1409+(29.3044*log10(MHillas.fSize)))"); 124 125 // Scaling from MC: 126 matrixg.AddColumn("(MHillas.fWidth*3.37e-3)/(-0.028235+(0.03231*log10(MHillas.fSize)))"); 127 matrixg.AddColumn("(MHillas.fLength*3.37e-3)/(-0.13527+(0.09876*log10(MHillas.fSize)))"); 128 129 matrixg.AddColumn("MNewImagePar.fConc"); 130 131 // TEST TEST TEST TEST 132 // matrixg.AddColumn("MConcentration.Conc7"); 133 134 // matrixg.AddColumn("MNewImagePar.fConc1"); 135 // matrixg.AddColumn("MNewImagePar.fAsym"); 136 // matrixg.AddColumn("MHillasExt.fM3Trans"); 137 // matrixg.AddColumn("abs(MHillasSrc.fHeadTail)"); 138 139 140 // We cannot use the 3rd moment with a sign referred to the nominal source position, if we 141 // assume we do not know it a priori. We cannot use either the sign as referred to the 142 // reconstructed source position through the DISP method, since that method uses precisely 143 // fM3Long to solve the "head-tail" asimetry and hence all showers would be "gamma-like" from 144 // that point of view. We add however the 3rd moment with no sign to the separation parameters. 145 146 matrixg.AddColumn("abs(MHillasExt.fM3Long)"); 147 148 // matrixg.AddColumn("MHillasSrc.fAlpha"); 149 150 TString dummy("(MHillas.fSize<"); 151 dummy += minsize; 152 153 // AM TEST 154 dummy += ") || (MNewImagePar.fNumCorePixels<"; 155 dummy += mincorepixels; 156 157 // AM, useful FOR High Energy only: 158 // dummy += ") || (MImagePar.fNumSatPixelsLG>0"; 159 dummy += ")"; 160 161 MContinue SizeCut(dummy); 162 tlist.AddToList(&SizeCut); 163 164 165 if (skipevents != "") 166 { 167 MContinue SCut(skipevents); 168 tlist.AddToList(&SCut); 169 } 170 171 plist.AddToList(&matrixg); 172 173 MHMatrix matrixh("MatrixHadrons"); 174 matrixh.AddColumns(matrixg.GetColumns()); 175 plist.AddToList(&matrixh); 176 177 MFEventSelector reduce_training_hadrons; 178 reduce_training_hadrons.SetSelectionRatio(hadron_frac); 179 MFilterList hadfilter; 180 hadfilter.AddToList(&reduce_training_hadrons); 181 hadfilter.AddToList(&fhadrons); 182 tlist.AddToList(&hadfilter); 183 184 MFillH fillmath("MatrixHadrons"); 185 fillmath.SetFilter(&hadfilter); 186 tlist.AddToList(&fillmath); 187 188 MContinue skiphad(&fhadrons); 189 tlist.AddToList(&skiphad); 190 191 MFEventSelector reduce_training_gammas; 192 reduce_training_gammas.SetSelectionRatio(gamma_frac); 193 tlist.AddToList(&reduce_training_gammas); 194 195 MFillH fillmatg("MatrixGammas"); 196 fillmatg.SetFilter(&reduce_training_gammas); 197 tlist.AddToList(&fillmatg); 198 199 // 200 // Create and setup the eventloop 201 // 202 MEvtLoop evtloop; 203 evtloop.SetParList(&plist); 204 205 MProgressBar *bar = new MProgressBar; 206 evtloop.SetProgressBar(bar); 207 // 208 // Execute your analysis 209 // 210 if (!evtloop.Eventloop()) 211 return; 212 213 tlist.PrintStatistics(); 214 215 216 // --------------------------------------------------------------- 217 // --------------------------------------------------------------- 218 // second event loop: the trees of the random forest are grown, 219 // the event loop is now actually a tree loop (loop of training 220 // process) 221 // --------------------------------------------------------------- 222 // --------------------------------------------------------------- 223 224 MTaskList tlist2; 225 plist.Replace(&tlist2); 226 227 // 228 // matrixh.ShuffleRows(9999); 229 // matrixg.ShuffleRows(1111); 230 // 231 232 MRanForestGrow rfgrow2; 233 rfgrow2.SetNumTrees(100); 234 rfgrow2.SetNumTry(3); 235 rfgrow2.SetNdSize(10); 236 237 MWriteRootFile rfwrite2("RF.root"); 238 rfwrite2.AddContainer("MRanTree","Tree"); //save all trees 239 MFillH fillh2("MHRanForestGini"); 240 241 tlist2.AddToList(&rfgrow2); 242 tlist2.AddToList(&rfwrite2); 243 tlist2.AddToList(&fillh2); 244 245 // gRandom is accessed from MRanForest (-> bootstrap aggregating) 246 // and MRanTree (-> random split selection) and should be initialized 247 // here if you want to set a certain random number generator 248 if(gRandom) 249 delete gRandom; 250 // gRandom = new TRandom3(0); // Takes seed from the computer clock 251 gRandom = new TRandom3(); // Uses always same seed 252 // 253 // Execute tree growing (now the eventloop is actually a treeloop) 254 // 255 256 evtloop.SetProgressBar(bar); 257 if (!evtloop.Eventloop()) 258 return; 259 260 tlist2.PrintStatistics(); 261 262 plist.FindObject("MHRanForestGini")->DrawClone(); 263 264 // --------------------------------------------------------------- 265 // --------------------------------------------------------------- 266 // third event loop: the control sample (star2.root) is processed 267 // through the previously grown random forest, 268 // 269 // the histograms MHHadronness (quality of g/h-separation) and 270 // MHRanForest are created and displayed. 271 // MHRanForest shows the convergence of the classification error 272 // as function of the number of grown (and combined) trees 273 // and tells the user how many trees are actually needed in later 274 // classification tasks. 275 // --------------------------------------------------------------- 276 // --------------------------------------------------------------- 277 278 MTaskList tlist3; 279 280 plist.Replace(&tlist3); 281 282 283 MReadMarsFile read3("Events", "star_20050110_CrabNebulaW1.root"); 284 // MReadMarsFile read3("Events", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Calib_6slices/tc_3.84_2.56/star_gamma_test.root"); 285 286 read3.DisableAutoScheme(); 287 tlist3.AddToList(&read3); 288 tlist3.AddToList(&SizeCut); 289 290 if (skipevents != "") 291 tlist3.AddToList(&SCut); 292 293 MRanForestCalc calc; 294 tlist3.AddToList(&calc); 295 296 297 TString outfname = "star_S"; 298 outfname += minsize; 299 outfname += "_20050110_CrabNebulaW1.root"; 300 301 // outfname += "_mcgammas.root"; 302 303 // parameter containers you want to save 304 // 305 306 MWriteRootFile write(outfname, "recreate"); 307 308 // write.AddContainer("MSigmabar", "Events"); 309 310 write.AddContainer("MRawEvtHeader", "Events"); 311 write.AddContainer("MMcEvt", "Events", kFALSE); 312 write.AddContainer("MHillas", "Events"); 313 write.AddContainer("MHillasExt", "Events"); 314 write.AddContainer("MImagePar", "Events"); 315 write.AddContainer("MNewImagePar", "Events"); 316 write.AddContainer("MHillasSrc", "Events"); 317 write.AddContainer("MHadronness", "Events"); 318 write.AddContainer("MPointingPos", "Events"); 319 write.AddContainer("MTime", "Events", kFALSE); 320 write.AddContainer("MConcentration", "Events", kFALSE); 321 322 tlist3.AddToList(&write); 323 324 evtloop.SetProgressBar(bar); 325 326 // 327 // Execute your analysis 328 // 329 if (!evtloop.Eventloop()) 330 return; 331 332 tlist3.PrintStatistics(); 333 334 // bar->DestroyWindow(); 335 336 return; 337 } 1 matrixg.AddColumn("(1.-(MHillas.fWidth/MHillas.fLength))*(0.9776+(0.101062*log10(MHillas.fSize)))");
Note:
See TracChangeset
for help on using the changeset viewer.