Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7400) +++ trunk/MagicSoft/Mars/Changelog (revision 7401) @@ -18,4 +18,16 @@ -*-*- END OF LINE -*-*- + 2005/11/15 Thomas Bretz + + * macros/optim/optimdisp.C, macros/optim/optimenergy.C, + macros/optim/rfenergyest.C: + - added some more comments + - updated + + * mtools/MTFillMatrix.[h,cc]: + - added the possibility to use PreCuts like in MJOptim + + + 2005/11/14 Daniela Dorner Index: trunk/MagicSoft/Mars/macros/optim/optimdisp.C =================================================================== --- trunk/MagicSoft/Mars/macros/optim/optimdisp.C (revision 7400) +++ trunk/MagicSoft/Mars/macros/optim/optimdisp.C (revision 7401) @@ -1,27 +1,19 @@ +/* More explanations at the end of the file */ void optimdisp() { - MJOptimizeDisp opt; - opt.SetDebug(2); - opt.SetOptimizer(MJOptimize::kSimplex); - opt.EnableTestTrain(); + MJOptimizeDisp opt; // Initialize optimization class + opt.SetDebug(2); // Choose the level of output + opt.SetOptimizer(MJOptimize::kMigrad); // Choose the fit algorithm (see MJOptimize) + opt.EnableTestTrain(); // Split sample into test and train // -------------------- Setup ---------------------------- - //opt.AddParameter("1-(MHillas.fWidth/MHillas.fLength)"); - opt.AddParameter("atan(1./(1-(MHillas.fWidth/MHillas.fLength)))"); - //opt.AddParameter("log10(MNewImagePar.fLeakage1+1)"); + opt.AddParameter("1-(MHillas.fWidth/MHillas.fLength)"); // M[0] + opt.AddParameter("log10(MNewImagePar.fLeakage1+1)"); // M[1] - opt.FixParameter(0, 0.743871, 0, 2); //0.562676 - opt.FixParameter(1, 0.329471, 0, 2); //0.317315 - opt.FixParameter(2, 1.500002, 0, 10); //1.55754 + opt.SetParameter(0, 1.30871); // Setup [0] + opt.SetParameter(1, 5.81119); // Setup [1] + opt.SetParameter(2, 0.763486); // Setup [2] - // Leakage==0: RMS = 0.030 - // Leakage>0: RMS = 0.010 - - // all: 0.029 (0.605991, 0.332266, 1.50825) - - //char *r = "([0]+(pow(M[1],[1])*[2]))/M[0]"; - //char *r = "([0]+(pow(M[1],[1])*[2]))*M[0]"; - //char *r = "[0]*M[0]"; - char *r = "[0]/M[0]"; + char *r = "([0]+([1]*pow(M[1], [2])))*M[0]"; // Rule to calc Disp // -------------------- Run ---------------------------- @@ -30,5 +22,47 @@ opt.SetDisplay(d); - opt.AddPreCut("MNewImagePar.fLeakage1<0.0001"); + /* + -------------------- Magic-Cuts ---------------------- + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kArea); + cuts.SetThetaCut(MFMagicCuts::kOn); + + TArrayD arr(10); + arr[0]= 1.3245; + arr[1]= 0.208700; + arr[2]= 0.229200; + arr[3]= 5.305200; + arr[4]= 0.098930; + arr[5]= -0.082950; + arr[6]= 8.2957; + arr[7]= 0.8677; + + cuts.SetVariables(arr); + + opt.AddPreCut(&cuts); + + -------------------- Other cuts ---------------------- + opt.AddPreCut("MNewImagePar.fLeakage1<0.0001"); + */ + opt.RunDisp("ganymedmcpart.root", r); } + +/* ------------------ Good strategy ------------------- + 1) first fix parameters: + opt.FixParameter(1, 0); + opt.FixParameter(2, 1); + and process only showers without leakage + opt.AddPreCut("MNewImagePar.fLeakage1<0.0001"); + 2) release parameters 1 and 2 and fix 0 to the result of 1) + opt.FixParameter(0, 0.8362); + opt.SetParameter(1, 2.0); + opt.SetParameter(2, 0.8); + and process only showers with leakage + opt.AddPreCut("MNewImagePar.fLeakage1>0.0001"); + 3) release all parameters and start with the result of 1) and 2) + opt.SetParameter(0, 0.8362); + opt.SetParameter(1, 5.84); + opt.SetParameter(2, 0.76); + and process all showers. (Comment out opt.PreCuts()) + */ Index: trunk/MagicSoft/Mars/macros/optim/optimenergy.C =================================================================== --- trunk/MagicSoft/Mars/macros/optim/optimenergy.C (revision 7400) +++ trunk/MagicSoft/Mars/macros/optim/optimenergy.C (revision 7401) @@ -18,5 +18,28 @@ opt.SetDisplay(d); - //opt.AddPreCut("MNewImagePar.fLeakage1<0.0001"); + /* + -------------------- Magic-Cuts ---------------------- + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kArea); + cuts.SetThetaCut(MFMagicCuts::kOn); + + TArrayD arr(10); + arr[0]= 1.3245; + arr[1]= 0.208700; + arr[2]= 0.229200; + arr[3]= 5.305200; + arr[4]= 0.098930; + arr[5]= -0.082950; + arr[6]= 8.2957; + arr[7]= 0.8677; + + cuts.SetVariables(arr); + + opt.AddPreCut(&cuts); + + -------------------- Other cuts ---------------------- + opt.AddPreCut("MPointingPos.fZd<7"); + */ + opt.RunEnergy("ganymedmcpart.root", r); } Index: trunk/MagicSoft/Mars/macros/optim/rfenergyest.C =================================================================== --- trunk/MagicSoft/Mars/macros/optim/rfenergyest.C (revision 7400) +++ trunk/MagicSoft/Mars/macros/optim/rfenergyest.C (revision 7401) @@ -1,6 +1,6 @@ void rfenergyest() { - MSequence seqtst("~/Software/Mars/mranforest/sequencemc-test.txt"); - MSequence seqtrn("~/Software/Mars/mranforest/sequencemc-train.txt"); + MSequence seqtst("~/Software/mc/sequencemc-test.txt"); + MSequence seqtrn("~/Software/mc/sequencemc-train.txt"); if (!seqtst.IsValid()) @@ -22,6 +22,6 @@ MDirIter iter, iter2; - seqtrn.SetupDatRuns(iter, MSequence::kImages, "~/Software/mc"); - seqtst.SetupDatRuns(iter2, MSequence::kImages, "~/Software/mc"); + seqtrn.SetupDatRuns(iter, MSequence::kImages, "~/Software/mc/img-abs"); + seqtst.SetupDatRuns(iter2, MSequence::kImages, "~/Software/mc/img-abs"); read.AddFiles(iter); @@ -60,4 +60,27 @@ fill.SetReader(&read); + /* ---------- It doesn't seem to improve anything ---------- + MFMagicCuts cuts; + + cuts.SetHadronnessCut(MFMagicCuts::kArea); + cuts.SetThetaCut(MFMagicCuts::kOn); + + TArrayD arr(10); + arr[0]= 1.3245; + arr[1]= 0.208700; + arr[2]= 0.229200; + arr[3]= 5.305200; + arr[4]= 0.098930; + arr[5]= -0.082950; + arr[6]= 8.2957; + arr[7]= 0.8677; + + cuts.SetVariables(arr); + + fill.AddPreCut(&cuts); + + --------------- Use the cuts also in test-loop ---------------- + */ + if (!fill.Process()) return; @@ -71,11 +94,11 @@ MBinning b(32, 10, 100000, "BinningEnergyEst", "log"); /* - if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification + if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification with one tree per bin return; - if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification + if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification into different bins return; */ - if (!rf.TrainSingleRF(train)) // regression + if (!rf.TrainSingleRF(train)) // regression (best choice) return; @@ -90,7 +113,10 @@ MHEnergyEst hist; + //MContinue cont(&cuts); + //cont.SetInverted(); MFillH fillh(&hist); tlist.AddToList(&read2); + //tlist.AddToList(&cont); tlist.AddToList(&rf); tlist.AddToList(&fillh); Index: trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc (revision 7400) +++ trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc (revision 7401) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 12/2003 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -81,4 +81,5 @@ #include +// environment #include "MHMatrix.h" @@ -86,11 +87,16 @@ #include "MLogManip.h" +// eventloop #include "MParList.h" #include "MTaskList.h" #include "MEvtLoop.h" +// tasks #include "MRead.h" #include "MFillH.h" #include "MContinue.h" + +// filters +#include "MF.h" #include "MFilterList.h" #include "MFEventSelector.h" @@ -174,4 +180,28 @@ delete fReference; } + +//------------------------------------------------------------------------ +// +// Add a cut which is used to fill the matrix, eg "MMcEvt.fOartId<1.5" +// (The rule is applied, nit inverted: The matrix is filled with +// the events fullfilling the condition) +// +void MTFillMatrix::AddPreCut(const char *rule) +{ + MFilter *f = new MF(rule); + f->SetBit(kCanDelete); + AddPreCut(f); +} + +//------------------------------------------------------------------------ +// +// Add a cut which is used to fill the matrix. If kCanDelete is set +// MJOptimize takes the ownership. +// +void MTFillMatrix::AddPreCut(MFilter *f) +{ + fPreCuts.Add(f); +} + // -------------------------------------------------------------------------- @@ -236,4 +266,15 @@ // + // Continue for PreCuts + // + MFilterList list; + list.SetName("PreCuts"); + if (!list.AddToList(fPreCuts)) + *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl; + + MContinue cont0(&list); + cont0.SetInverted(); + + // // Continue for all events which are not (SetInverted()) // selected by the 'selector' @@ -268,4 +309,6 @@ // entries in MTaskList tlist.AddToList(fReader); // Read events + if (fPreCuts.GetEntries()>0) + tlist.AddToList(&cont0); // PreCuts if (fReference && selector) tlist.AddToList(&cont); // select a sample of events