Index: /trunk/MagicSoft/Mars/macros/optim/optimonoff.C =================================================================== --- /trunk/MagicSoft/Mars/macros/optim/optimonoff.C (revision 8705) +++ /trunk/MagicSoft/Mars/macros/optim/optimonoff.C (revision 8706) @@ -6,18 +6,24 @@ // ------------------- Xi ----------------------- - opt.FixParameter(0, 1.3245); - opt.FixParameter(6, 8.2957); - opt.FixParameter(7, 0.8677); - + opt.FixParameter( 0, 1.15136); + opt.FixParameter( 8, 0.0681437); + opt.FixParameter( 9, 2.62932); + opt.FixParameter(10, 1.51279); + opt.FixParameter(11, 0.0507821); + // --------------- Theta Sq --------------------- - opt.SetParameter(1, 0.204); + opt.FixParameter( 1, 0.215); // ----------------- Area ----------------------- - opt.SetParameter(2, 0.215); - opt.SetParameter(3, 5.486); - opt.SetParameter(4, 0.0897); + opt.SetParameter( 2, 0.21); + opt.SetParameter( 3, 5.6); + opt.SetParameter( 4, 0.083); // ----------------- M3L ------------------------ - opt.FixParameter(5, -0.07); + opt.FixParameter( 5, -0.07); + + // ---------------- Slope ----------------------- + opt.FixParameter( 6, 7.2); + opt.FixParameter( 7, 0.5); // ----------------- Fit ------------------------ Index: /trunk/MagicSoft/Mars/macros/optim/optimwobble.C =================================================================== --- /trunk/MagicSoft/Mars/macros/optim/optimwobble.C (revision 8705) +++ /trunk/MagicSoft/Mars/macros/optim/optimwobble.C (revision 8706) @@ -6,18 +6,24 @@ // ------------------- Xi ----------------------- - opt.FixParameter(0, 1.221); - opt.FixParameter(6, 13.425); - opt.FixParameter(7, 1.04343); + opt.FixParameter( 0, 1.15136); + opt.FixParameter( 8, 0.0681437); + opt.FixParameter( 9, 2.62932); + opt.FixParameter(10, 1.51279); + opt.FixParameter(11, 0.0507821); // --------------- Theta Sq --------------------- - opt.SetParameter(1, 0.23); + opt.FixParameter( 1, 0.215); // ----------------- Area ----------------------- - opt.SetParameter(2, 0.25); - opt.SetParameter(3, 5.21); - opt.SetParameter(4, 0.11); + opt.SetParameter( 2, 0.21); + opt.SetParameter( 3, 5.6); + opt.SetParameter( 4, 0.083); // ----------------- M3L ------------------------ - opt.SetParameter(5, -0.089); + opt.FixParameter( 5, -0.07); + + // ---------------- Slope ----------------------- + opt.FixParameter( 6, 7.2); + opt.FixParameter( 7, 0.5); // ----------------- Fit ------------------------ @@ -28,8 +34,12 @@ fit.SetMinimizationStrategy(MAlphaFitter::kSignificanceLogExcess); + // For wobble mode with three off-regions + // (sets also ScaleMode to kUser) + fit.SetScaleUser(1./3); + // ----------------- Cuts ----------------------- MFMagicCuts cuts; - cuts.SetThetaCut(MFMagicCuts::kWobble); + cuts.SetThetaCut(MFMagicCuts::kOn); /* You can also try (switch on scaling): @@ -47,5 +57,4 @@ opt.SetDisplay(d); - // opt.SetPathOut("optimwobble.root"); - opt.RunOnOff("wobblelza-abs.root", &cuts, &fit); + opt.RunOnOff("ganymed00000001-summary.root", &cuts, &fit); } Index: /trunk/MagicSoft/Mars/macros/train/traindisp.C =================================================================== --- /trunk/MagicSoft/Mars/macros/train/traindisp.C (revision 8705) +++ /trunk/MagicSoft/Mars/macros/train/traindisp.C (revision 8706) @@ -1,5 +1,54 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz, 11/2005 +! +! Copyright: MAGIC Software Development, 2000-2007 +! +! +\* ======================================================================== *////////////////////////////////////////////////////////////////////////////// + +///////////////////////////////////////////////////////////////////////////// +// +// traindisp.C +// =========== +// +// Trains a random forest for disp estimation. +// +// The additional code below shows how the MagicCuts can be used in +// training and testing. There is also code which allows to change the +// slope of the input spectrum. If a min- or max-break energy is given +// the new slope is only applied in this region. Further possibilities +// for additional cuts are shown. +// +// SequencesOn: Monte Carlo Sequences used for training +// SequencesOff: Monte Carlo Sequences used for testing +// +// Use: +// opt.AddPreCut to use cut for test and training +// opt.AddTestCut to use cut for test only +// opt.AddTrainCut to use cut for train only +// opt.AddPreTask a task executed before filling the matrix/evaluating the RF +// opt.AddPostTask a task executed before training, after evaluating the RF +// +///////////////////////////////////////////////////////////////////////////// + void traindisp() { - MDataSet set("mctesttrain.txt"); + MDataSet set("mcsponde/mcdataset-for-training.txt"); set.SetNumAnalysis(1); // Necessary @@ -10,4 +59,9 @@ //opt.SetDebug(); + // These are control parameters how to train the random forest (use with care!) + //opt.SetNumTrees(100); + //opt.SetNdSize(5); + //opt.SetNumTry(0); + // ------- Parameters to train Random Forest -------- opt.AddParameter("MHillas.fLength"); @@ -15,15 +69,16 @@ opt.AddParameter("MHillas.fSize"); opt.AddParameter("MNewImagePar.fLeakage1"); - opt.AddParameter("MNewImagePar.fLeakage2"); - opt.AddParameter("MNewImagePar.fConc"); - opt.AddParameter("MNewImagePar.fConc1"); opt.AddParameter("MPointingPos.fZd"); - opt.AddParameter("MNewImagePar.fUsedArea"); - opt.AddParameter("MNewImagePar.fCoreArea"); - opt.AddParameter("MNewImagePar.fNumUsedPixels"); - opt.AddParameter("MNewImagePar.fNumCorePixels"); - opt.AddParameter("MImagePar.fSizeSinglePixels"); - opt.AddParameter("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)"); - opt.AddParameter("MHillasExt.fAsym*sign(MHillasSrc.fCosDeltaAlpha)"); + opt.AddParameter("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)"); + //opt.AddParameter("MNewImagePar.fLeakage2"); + //opt.AddParameter("MNewImagePar.fConc"); + //opt.AddParameter("MNewImagePar.fConc1"); + //opt.AddParameter("MNewImagePar.fUsedArea"); + //opt.AddParameter("MNewImagePar.fCoreArea"); + //opt.AddParameter("MNewImagePar.fNumUsedPixels"); + //opt.AddParameter("MNewImagePar.fNumCorePixels"); + //opt.AddParameter("MImagePar.fSizeSinglePixels"); + //opt.AddParameter("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)"); + //opt.AddParameter("MHillasExt.fAsym*sign(MHillasSrc.fCosDeltaAlpha)"); // -------------------- Run ---------------------------- @@ -38,13 +93,17 @@ 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; + TArrayD arr(13); + arr[0] = 1.15136; + arr[1] = 0.215; + arr[2] = 0.215468; + arr[3] = 5.63973; + arr[4] = 0.0836169; + arr[5] = -0.07; + arr[6] = 7.2; + arr[7] = 0.5; + arr[8] = 0.0681437; + arr[9] = 2.62932; + arr[10] = 1.51279; + arr[11] = 0.0507821; cuts.SetVariables(arr); @@ -58,8 +117,46 @@ -------------------- Other cuts ---------------------- opt.AddPreCut("MNewImagePar.fLeakage1<0.0001"); - */ - opt.AddPreCut("MHillas.fSize>200"); opt.AddPreCut("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg<0.3"); + opt.AddPreCut("MPointingPos.fZd<25"); + + ------------------ Zd distribution ------------------- + TFile file("ganymed00001111.root"); + + MStatusArray arr; + if (arr.Read()<=0) + return; + TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime"); + if (!vstime) + return; + + MMcSpectrumWeight weight; + weight.SetWeightsZd(vstime); + opt.AddPreTask(&weight); + + ---------------------- Histogram -------------------- + + MHn hist("MyHist", "Energy Residual (lg E_{est} - lg E_{mc})"); + + MBinning bins(50, -0.5, 0.5); + + hist.AddHist("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*3.37e-3", "Disp.fVal-MHillasSrc.fDist*3.37e-3"); + hist.InitName("ResM3l;M3Long;ResidualDist"); + hist.InitTitle(";M3l [\\circ];Disp-Dist [\\circ];"); + hist.SetDrawOption("colz profx"); + hist.SetBinnings(&bins); + + hist.AddHist("MHillasExt.fAsym*sign(MHillasSrc.fCosDeltaAlpha)*3.37e-3", "Disp.fVal-MHillasSrc.fDist*3.37e-3"); + hist.InitName("ResAsym;M3Long;ResidualDist"); + hist.InitTitle(";Asym [\\circ];Disp-Dist [\\circ];"); + hist.SetDrawOption("colz profx"); + hist.SetBinnings(&bins); + + MFillH fill(&hist, "", "FillMyHist"); + //fill.SetWeight(); // Enable weights to be used to fill this histogram + opt.AddTestTask(&fill); + + ------------------------------------------------------ + */ // The number is the number of events read from the file @@ -67,11 +164,2 @@ opt.Train("rf-disp.root", set, 30000); } -/* - // SequencesOn: Monte Carlo Sequences used for training - // SequencesOff: Monte Carlo Sequences used for testing - - // Use: - // opt.AddPreCut to use cut for test and training - // opt.AddTestCut to use cut for test only - // opt.AddTrainCut to use cut for train only - */ Index: /trunk/MagicSoft/Mars/macros/train/trainenergy.C =================================================================== --- /trunk/MagicSoft/Mars/macros/train/trainenergy.C (revision 8705) +++ /trunk/MagicSoft/Mars/macros/train/trainenergy.C (revision 8706) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 11/2005 ! -! Copyright: MAGIC Software Development, 2000-2006 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -36,8 +36,18 @@ // for additional cuts are shown. // +// SequencesOn: Monte Carlo Sequences used for training +// SequencesOff: Monte Carlo Sequences used for testing +// +// Use: +// opt.AddPreCut to use cut for test and training +// opt.AddTestCut to use cut for test only +// opt.AddTrainCut to use cut for train only +// opt.AddPreTask a task executed before filling the matrix/evaluating the RF +// opt.AddPostTask a task executed before training, after evaluating the RF +// ///////////////////////////////////////////////////////////////////////////// void trainenergy() { - MDataSet set("mcdataset.txt"); + MDataSet set("mcsponde/mcdataset-for-training.txt"); set.SetNumAnalysis(1); // Necessary @@ -48,15 +58,30 @@ //opt.SetDebug(); + // These are control parameters how to train the random forest (use with care!) + //opt.SetNumTrees(100); + //opt.SetNdSize(5); + //opt.SetNumTry(0); + // ------- Parameters to train Random Forest -------- + // The following parameters are the best parameters opt.AddParameter("MHillas.fSize"); opt.AddParameter("MHillasSrc.fDist"); opt.AddParameter("MPointingPos.fZd"); - opt.AddParameter("MHillas.GetArea"); - opt.AddParameter("MNewImagePar.fUsedArea"); - opt.AddParameter("MNewImagePar.fCoreArea"); opt.AddParameter("MNewImagePar.fLeakage1"); - opt.AddParameter("MNewImagePar.fLeakage2"); - opt.AddParameter("MNewImagePar.fConc"); - opt.AddParameter("MNewImagePar.fConc1"); + opt.AddParameter("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)"); + // The influence of these parameters is unclear + //opt.AddParameter("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)"); + //opt.AddParameter("MNewImagePar.fConcCOG"); + //opt.AddParameter("MHillas.GetArea"); + //opt.AddParameter("MNewImagePar.fConc1"); + //opt.AddParameter("MNewImagePar.fConc"); + //opt.AddParameter("MNewImagePar.fUsedArea"); + //opt.AddParameter("MNewImagePar.fCoreArea"); + //opt.AddParameter("MNewImagePar.fLeakage2"); + + // Setup how to train the RF. This one gives best results so far + opt.SetTrainFunc("log(MMcEvt.fEnergy)/log(MHillas.fSize)", "MHillas.fSize^x"); + // opt.SetTrainLog(); // Train in log-energy + // opt.SetTrainLin(); // Train just in energy // -------------------- Run ---------------------------- @@ -65,50 +90,78 @@ opt.SetDisplay(d); + // -------------------- Magic-Cuts ---------------------- + // It is recommended to test with your cuts, but train with all events + + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kArea); + cuts.SetThetaCut(MFMagicCuts::kOn); + + TArrayD arr(13); + arr[0] = 1.15136; + arr[1] = 0.215; + arr[2] = 0.215468; + arr[3] = 5.63973; + arr[4] = 0.0836169; + arr[5] = -0.07; + arr[6] = 7.2; + arr[7] = 0.5; + arr[8] = 0.0681437; + arr[9] = 2.62932; + arr[10] = 1.51279; + arr[11] = 0.0507821; + + cuts.SetVariables(arr); + + // opt.AddPreCut(&cuts); // add cut for training and testing + // opt.AddTrainCut(&cuts); // add cut for training only + opt.AddTestCut(&cuts); // add cut for testing only + /* - -------------------- Magic-Cuts ---------------------- - MFMagicCuts cuts; - cuts.SetHadronnessCut(MFMagicCuts::kArea); - cuts.SetThetaCut(MFMagicCuts::kOn); + -------------------- Energy Slope -------------------- + // This is a way to throw away events to a different slope + MFEnergySlope slope(-4.0); // New slope for mc spectrum + opt.AddPreCut(&slope); // throw away events to change slope - 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); - - -------------------- Energy Slope -------------------- - MFEnergySlope slope(-2.8); // New slope for mc spectrum - slope.SetMcMinEnergy(80); // Set break energy from -2.6 to -2.8 - opt.AddPreCut(&slope); // throw away events to change slope + // This is a way to weight the events to a different spectrum + MMcSpectrumWeight weight; + weight.SetFormula("pow(X/300, -2.31-0.26*log10(X/300))"); + opt.SetWeights(&weight); -------------------- Other cuts ---------------------- opt.AddPreCut("MHillasSrc.fDist*MGeomCam.fConvMm2Deg<1.0"); opt.AddPreCut("MHillas.fSize>200"); + + ------------------ Zd distribution ------------------- + TFile file("ganymed00001111.root"); + + MStatusArray arr; + if (arr.Read()<=0) + return; + TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime"); + if (!vstime) + return -1; + + MMcSpectrumWeight weight; + weight.SetWeightsZd(vstime); + opt.AddPreTask(&weight); + + ---------------------- Histogram -------------------- + MHn hist("MyHist", "Energy Residual (lg E_{est} - lg E_{mc})"); + + hist.AddHist("MNewImagePar.fConcCOG", "log10(MEnergyEst.fVal)-log10(MMcEvt.fEnergy)"); + hist.InitName("ResConc;Conc;EnergyResidual"); + hist.InitTitle(";C;\\Delta lg E;"); + hist.SetDrawOption("colz profx"); + + MBinning bins(50, 0, 1); + hist.SetBinnings(&bins); + + MFillH fill(&hist, "", "FillMyHist"); + //fill.SetWeight(); // Enable weights to be used to fill this histogram + opt.AddTestTask(&fill); + + ------------------------------------------------------ */ - - // Things which could be implemented in a future version - // PreCut in Energy oder Size? - // PreCut in Max Dist - // PreCut in Theta^2 (0.25 oder 0.3) - // New spectral slope - // Zenith angle distribution in OnTime opt.Train("rf-energy.root", set, 30000); } -/* - // SequencesOn: Monte Carlo Sequences used for training - // SequencesOff: Monte Carlo Sequences used for testing - - // Use: - // opt.AddPreCut to use cut for test and training - // opt.AddTestCut to use cut for test only - // opt.AddTrainCut to use cut for train only - */ Index: /trunk/MagicSoft/Mars/mfbase/MFDataPhrase.cc =================================================================== --- /trunk/MagicSoft/Mars/mfbase/MFDataPhrase.cc (revision 8705) +++ /trunk/MagicSoft/Mars/mfbase/MFDataPhrase.cc (revision 8706) @@ -107,8 +107,11 @@ // the class description above. // -MFDataPhrase::MFDataPhrase(const char *text, const char *name, const char *title) : fAllowEmpty(kFALSE) +MFDataPhrase::MFDataPhrase(const char *text, const char *name, const char *title) : fData(NULL), fAllowEmpty(kFALSE) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); + + if (!text) + return; *fLog << inf << "Parsing filter rule... " << flush; @@ -162,6 +165,5 @@ Int_t MFDataPhrase::Process() { - if (fData) - fResult = fData->GetValue()>0.5; + fResult = fData ? fData->GetValue()>0.5 : fDefault; return kTRUE; Index: /trunk/MagicSoft/Mars/mfbase/MFDataPhrase.h =================================================================== --- /trunk/MagicSoft/Mars/mfbase/MFDataPhrase.h (revision 8705) +++ /trunk/MagicSoft/Mars/mfbase/MFDataPhrase.h (revision 8706) @@ -18,4 +18,5 @@ MData *fData; Bool_t fAllowEmpty; // Not a Bit to be stored + Bool_t fDefault; // Default for empty case Bool_t fResult; //! @@ -38,4 +39,7 @@ void SetAllowEmpty(Bool_t b=kTRUE) { fAllowEmpty = b; } Bool_t IsAllowEmpty() const { return fAllowEmpty; } + + void SetDefault(Bool_t b) { fDefault = b; } + Bool_t GetDefault() const { return fDefault; } // MFilter Index: /trunk/MagicSoft/Mars/mjobs/MJCut.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 8705) +++ /trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 8706) @@ -666,6 +666,5 @@ // Filter for VsSize - MFDataPhrase ftheta("", "CutT"); - ftheta.SetAllowEmpty(); + MFDataPhrase ftheta(0, "CutT"); // ------------- Loop Off Data -------------------- Index: /trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 8705) +++ /trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 8706) @@ -53,5 +53,5 @@ void PrintSetup(const MAlphaFitter &fit) const; Bool_t DisplayResult(const TH2D &mh1) const; - Bool_t IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const; + //Bool_t IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const; TArrayD FitSpectrum(TH1D &spectrum) const; TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const; Index: /trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 8705) +++ /trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 8706) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MRanForestCalc.cc,v 1.29 2007-08-24 08:33:48 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MRanForestCalc.cc,v 1.30 2007-08-24 12:58:49 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -167,12 +167,17 @@ // In the case weights should be used initialize the // corresponding array + Double_t sum = 0; + TArrayF weights(nrows); if (fLastDataColumnHasWeights) + { for (Int_t j=0; j