Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 5830) +++ trunk/MagicSoft/Mars/Changelog (revision 5831) @@ -21,9 +21,15 @@ -*-*- END OF LINE -*-*- + + 2005/01/14 Abelardo Moralejo + + * mtemp/mpadova/macros/trainsubsample.C, RanForestDISP.C + - added. Two macros used in the analysis of real wobble mode data. + + 2005/01/14 Daniela Dorner * macros/sql/filldotrun.C - added new arehucas-verions and changed code accordingly - Index: trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestDISP.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestDISP.C (revision 5831) +++ trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestDISP.C (revision 5831) @@ -0,0 +1,313 @@ +/* ======================================================================== *\ + ! + ! * + ! * This file is part of MARS, the MAGIC Analysis and Reconstruction + ! * Software. It is distributed to you in the hope that it can be a useful + ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. + ! * It is distributed WITHOUT ANY WARRANTY. + ! * + ! * Permission to use, copy, modify and distribute this software and its + ! * documentation for any purpose is hereby granted without fee, + ! * provided that the above copyright notice appear in all copies and + ! * that both that copyright notice and this permission notice appear + ! * in supporting documentation. It is provided "as is" without express + ! * or implied warranty. + ! * + ! + ! + ! Author(s): Abelardo Moralejo 1/2005 ! + ! + ! Copyright: MAGIC Software Development, 2000-2003 + ! + ! + \* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////////// +// +// macro RanForestDISP.C +// +// Version of Random Forest macro, intended to run on real data without assuming +// any particular source position in the camera FOV. So we do not use the DIST +// parameter, but rather (one version of) DISP method. We do not use either the +// sign of the third longitudinal moment with respect to the center of the camera +// (since the source can be somewhere else like in the case of wobble mode). +// +////////////////////////////////////////////////////////////////////////////////// + + +#include "MMcEvt.hxx" + +void RanForestDISP(Int_t minsize=1000) +{ + Int_t mincorepixels = 6; // minimum number of core pixels. +// Int_t mincorepixels = 0; + + TString skipevents; + // skipevents = "(MMcEvt.fImpact>15000.)"; + // skipevents = "(MNewImagePar.fLeakage1>0.001)"; + // skipevents = "(MHillas.fSize>1000)"; + + // Skip calibration events, leaking events and flash events: + skipevents = "(MNewImagePar.fLeakage1>0.1) || ({1.5+log10(MNewImagePar.fConc)*2.33333}-{6.5-log10(MHillas.fSize)}>0.)"; + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + + MReadMarsFile read("Events", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Star_6slices_tc32/star_gamma_train_new.root"); + Float_t numgammas = read.GetEntries(); + + read.AddFile("/data2/magic/data/rootdata/Crab/Period021/20040915/star_train_20040915_CrabNebulaW1.root"); + Float_t numhadrons = read.GetEntries() - numgammas; + + // Fraction of gammas to be used in training: + // Float_t gamma_frac = 1.5*(numhadrons / numgammas); + + Float_t gamma_frac = 1.; + + read.DisableAutoScheme(); + + tlist.AddToList(&read); + + MFParticleId fgamma("MMcEvt", '=', MMcEvt::kGAMMA); + tlist.AddToList(&fgamma); + + MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA); + tlist.AddToList(&fhadrons); + + + MHMatrix matrixg("MatrixGammas"); + + // TEST TEST TEST +// matrixg.AddColumn("MMcEvt.fImpact"); +// matrixg.AddColumn("MMcEvt.fMuonCphFraction"); + + // matrixg.AddColumn("MNewImagePar.fLeakage1"); + // matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))"); + + // matrixg.AddColumn("MMcEvt.fTelescopePhi"); + // matrixg.AddColumn("MSigmabar.fSigmabar"); + + + matrixg.AddColumn("MHillas.fSize"); + + // We do not use DIST, since it depends on knowing the source position on camera. + // In wobble mode the position changes with time. Besides, we intend to do a 2-d map + // with the DISP method, and so we cannot use DIST in the cuts (since it would amde the + // camera center a "provoleged" position). + + matrixg.AddColumn("(1.-(MHillas.fWidth/MHillas.fLength))*(0.9776+0.101062*log10(MHillas.fSize))"); + +// matrixg.AddColumn("MHillas.fWidth"); +// matrixg.AddColumn("MHillas.fLength"); + + + // Scaling from real data? + // matrixg.AddColumn("MHillas.fWidth/(-8.37847+(9.58826*log10(MHillas.fSize)))"); + // matrixg.AddColumn("MHillas.fLength/(-40.1409+(29.3044*log10(MHillas.fSize)))"); + + // Scaling from MC: + matrixg.AddColumn("(MHillas.fWidth*3.37e-3)/(-0.028235+(0.03231*log10(MHillas.fSize)))"); + matrixg.AddColumn("(MHillas.fLength*3.37e-3)/(-0.13527+(0.09876*log10(MHillas.fSize)))"); + + matrixg.AddColumn("MNewImagePar.fConc"); + + // TEST TEST TEST TEST + // matrixg.AddColumn("MConcentration.Conc7"); + + // matrixg.AddColumn("MNewImagePar.fConc1"); + // matrixg.AddColumn("MNewImagePar.fAsym"); + // matrixg.AddColumn("MHillasExt.fM3Trans"); + // matrixg.AddColumn("abs(MHillasSrc.fHeadTail)"); + + + // We cannot use the 3rd moment with a sign referred to the nominal source position, if we + // assume we do not know it a priori. We cannot use either the sign as referred to the + // reconstructed source position through the DISP method, since that method uses precisely + // fM3Long to solve the "head-tail" asimetry and hence all showers would be "gamma-like" from + // that point of view. We add however the 3rd moment with no sign to the separation parameters. + + matrixg.AddColumn("abs(MHillasExt.fM3Long)"); + + // matrixg.AddColumn("MHillasSrc.fAlpha"); + + TString dummy("(MHillas.fSize<"); + dummy += minsize; + + // AM TEST + dummy += ") || (MNewImagePar.fNumCorePixels<"; + dummy += mincorepixels; + + // AM, useful FOR High Energy only: + // dummy += ") || (MImagePar.fNumSatPixelsLG>0"; + dummy += ")"; + + MContinue SizeCut(dummy); + tlist.AddToList(&SizeCut); + + + if (skipevents != "") + { + MContinue SCut(skipevents); + tlist.AddToList(&SCut); + } + + plist.AddToList(&matrixg); + + MHMatrix matrixh("MatrixHadrons"); + matrixh.AddColumns(matrixg.GetColumns()); + plist.AddToList(&matrixh); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&fhadrons); + tlist.AddToList(&fillmath); + + MFEventSelector reduce_training_gammas; + reduce_training_gammas.SetSelectionRatio(gamma_frac); + tlist.AddToList(&reduce_training_gammas); + MContinue skiphad(&fhadrons); + tlist.AddToList(&skiphad); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&reduce_training_gammas); + tlist.AddToList(&fillmatg); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + MProgressBar *bar = new MProgressBar; + evtloop.SetProgressBar(bar); + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // second event loop: the trees of the random forest are grown, + // the event loop is now actually a tree loop (loop of training + // process) + // --------------------------------------------------------------- + // --------------------------------------------------------------- + + MTaskList tlist2; + plist.Replace(&tlist2); + + // + // matrixh.ShuffleRows(9999); + // matrixg.ShuffleRows(1111); + // + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(100); + rfgrow2.SetNumTry(3); + rfgrow2.SetNdSize(10); + + MWriteRootFile rfwrite2("RF.root"); + rfwrite2.AddContainer("MRanTree","Tree"); //save all trees + MFillH fillh2("MHRanForestGini"); + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); + + // gRandom is accessed from MRanForest (-> bootstrap aggregating) + // and MRanTree (-> random split selection) and should be initialized + // here if you want to set a certain random number generator + if(gRandom) + delete gRandom; + // gRandom = new TRandom3(0); // Takes seed from the computer clock + gRandom = new TRandom3(); // Uses always same seed + // + // Execute tree growing (now the eventloop is actually a treeloop) + // + + evtloop.SetProgressBar(bar); + if (!evtloop.Eventloop()) + return; + + tlist2.PrintStatistics(); + + plist.FindObject("MHRanForestGini")->DrawClone(); + + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // third event loop: the control sample (star2.root) is processed + // through the previously grown random forest, + // + // the histograms MHHadronness (quality of g/h-separation) and + // MHRanForest are created and displayed. + // MHRanForest shows the convergence of the classification error + // as function of the number of grown (and combined) trees + // and tells the user how many trees are actually needed in later + // classification tasks. + // --------------------------------------------------------------- + // --------------------------------------------------------------- + + MTaskList tlist3; + + plist.Replace(&tlist3); + + + MReadMarsFile read3("Events", "/data2/magic/data/rootdata/Crab/Period021/20040915/star_20040915_CrabNebulaW1.root"); + + read3.DisableAutoScheme(); + tlist3.AddToList(&read3); + tlist3.AddToList(&SizeCut); + + if (skipevents != "") + tlist3.AddToList(&SCut); + + MRanForestCalc calc; + tlist3.AddToList(&calc); + + // parameter containers you want to save + // + + TString outfname = "star_S"; + outfname += minsize; + outfname += "_20040915_CrabNebulaW1.root"; + MWriteRootFile write(outfname, "recreate"); + + // write.AddContainer("MSigmabar", "Events"); + +// write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MHadronness", "Events"); + write.AddContainer("MPointingPos", "Events"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MConcentration", "Events"); + + tlist3.AddToList(&write); + + evtloop.SetProgressBar(bar); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist3.PrintStatistics(); + + // bar->DestroyWindow(); + + return; +} Index: trunk/MagicSoft/Mars/mtemp/mpadova/macros/trainsubsample.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpadova/macros/trainsubsample.C (revision 5831) +++ trunk/MagicSoft/Mars/mtemp/mpadova/macros/trainsubsample.C (revision 5831) @@ -0,0 +1,100 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! Author(s): Abelardo Moralejo 1/2005 ! +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////////////////////// +// +// Macro trainsubsample.C +// +// Reads a star file from real data, and makes a subsample of it, adding a "fake" +// MMcEvt container with MMcEvt.fPartId=14. This intended to produce, from real data, +// a "hadron training file" to be used in random forest (or other) g/h separation. The +// MMcEvt container is added so that whatever the method will consider these (real) +// showers to be all protons (before cuts, most will be, even in the case of strong +// source as Crab). +// +///////////////////////////////////////////////////////////////////////////////////// + +void trainsubsample() +{ + MParList plist; + MTaskList tlist; + + plist.AddToList(&tlist); + + MReadMarsFile read("Events", "star_20040915_CrabNebulaW1.root"); + + MFEventSelector evsel; + evsel.SetSelectionRatio(0.05); // Fraction of events in output + + read.DisableAutoScheme(); + tlist.AddToList(&read); + + MMcEvt* mcevt = new MMcEvt; + mcevt->SetPartId(14); + plist.AddToList(mcevt); + + MWriteRootFile write("star_train_20040915_CrabNebulaW1.root", "recreate"); + + write.AddContainer("MMcEvt", "Events", kFALSE); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MConcentration", "Events", kFALSE); + write.AddContainer("MPointingPos", "Events", kFALSE); + + + write.AddContainer("MGeomCam", "RunHeaders", kFALSE); + write.AddContainer("MMcConfigRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcFadcHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcTrigHeader", "RunHeaders", kFALSE); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + + write.SetFilter(&evsel); + + tlist.AddToList(&evsel); + tlist.AddToList(&write); + + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + MProgressBar *bar = new MProgressBar; + evtloop.SetProgressBar(bar); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + bar->DestroyWindow(); + + return; +}