Index: trunk/MagicSoft/Mars/macros/optim/optimdisp.C =================================================================== --- trunk/MagicSoft/Mars/macros/optim/optimdisp.C (revision 8655) +++ trunk/MagicSoft/Mars/macros/optim/optimdisp.C (revision 8656) @@ -8,12 +8,27 @@ // -------------------- Setup ---------------------------- - opt.AddParameter("1-(MHillas.fWidth/MHillas.fLength)"); // M[0] - opt.AddParameter("log10(MNewImagePar.fLeakage1+1)"); // M[1] - opt.SetParameter(0, 1.30871); // Setup [0] - opt.SetParameter(1, 5.81119); // Setup [1] - opt.SetParameter(2, 0.763486); // Setup [2] + opt.AddParameter("1-MHillas.fWidth/MHillas.fLength"); // M[0] + opt.AddParameter("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg"); // M[1] + opt.AddParameter("MNewImagePar.fLeakage1"); // M[2] + opt.AddParameter("log10(MHillas.fSize)"); // M[3] - char *r = "([0]+([1]*pow(M[1], [2])))*M[0]"; // Rule to calc Disp + // -------------- Parametrization -------------------- + + char *r = "M[0]*([0] + [1]*M[1] + [2]*M[2] + (M[3]>[3])*[4]*(M[3]-[3])^2)"; + + opt.FixParameter(0, 1.266195); + opt.FixParameter(1, 0.100577); + opt.FixParameter(2, 1.80309); + opt.FixParameter(3, 2.87177); + opt.FixParameter(4, 0.616823); + + opt.AddPreCut("MNewImagePar.fLeakage1>0"); + opt.AddPreCut("log10(MHillas.fSize)<3.2"); + + //opt.AddPreCut("abs(MHillasSrc.fDCA*MGeomCam.fConvMm2Deg)<0.2"); + //opt.AddPreCut("(MHillasSrc.fDist*MGeomCam.fConvMm2Deg-0.5)*7.2>MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg"); + opt.AddPreCut("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg>-0.07"); + opt.AddPreCut("DataType.fVal>0.5"); // -------------------- Run ---------------------------- @@ -26,15 +41,17 @@ MFMagicCuts cuts; cuts.SetHadronnessCut(MFMagicCuts::kArea); - cuts.SetThetaCut(MFMagicCuts::kOn); + cuts.SetThetaCut(MFMagicCuts::kNone); 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; + arr[2]= 0.0836169; + arr[3]= 5.63973; + arr[4]= 0.0844195; + arr[5]= -0.07; + arr[6]= 13.2; + arr[7]= 1.0; + arr[8]= 7.2; + arr[9]= 0.5; cuts.SetVariables(arr); @@ -47,26 +64,17 @@ -------------------- Other cuts ---------------------- - opt.AddPreCut("MNewImagePar.fLeakage1<0.0001"); - */ + */ - opt.RunDisp("ganymedmcpart.root", r); + opt.RunDisp("ganymed00000001-summary.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()) - */ + + Par | 0 1 2 3 4 | Cut + -------+---------------------------------+----------------------- + Fit 1 | 1.3 0.8 fix=0 fix=0 fix=0 | Leak1==0 lgSize<2.5 + Fit 2 | fix fix fix=0 2.4 0.3 | Leak1==0 + Fit 2 | fix fix 1.8 fix fix | Leak1>0 + Fit 3 | free free fix fix fix | -/- + -------+---------------------------------+----------------------- +*/ Index: trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.cc (revision 8655) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.cc (revision 8656) @@ -189,5 +189,5 @@ // Run Disp optimization // -Bool_t MJTrainDisp::Train(const char *out, const MDataSet &set, Int_t num) +Bool_t MJTrainDisp::TrainDisp(const char *out, const MDataSet &set, Int_t num) { SetTitle(Form("TrainDisp: %s", out)); @@ -224,5 +224,5 @@ if (fEnableWeights) train.AddColumn("MWeight.fVal"); - train.AddColumn(fgTrainParameter); + train.AddColumn(fTrainParameter); // ----------------------- Fill Matrix RF ---------------------- @@ -240,5 +240,5 @@ // ------------------------ Train RF -------------------------- - MRanForestCalc rf(fTitle); + MRanForestCalc rf("TrainDisp", fTitle); rf.SetNumTrees(fNumTrees); rf.SetNdSize(fNdSize); @@ -268,4 +268,6 @@ gLog.Separator("Test"); + MH::SetPalette("pretty"); + MParList plist; MTaskList tlist; @@ -329,6 +331,6 @@ MFillH fillh2a(&hdisp1, "", "FillSize"); MFillH fillh2b(&hdisp2, "", "FillEnergy"); - fillh2a.SetDrawOption("blue profx"); - fillh2b.SetDrawOption("blue profx"); + fillh2a.SetDrawOption("colz profx"); + fillh2b.SetDrawOption("colz profx"); fillh2a.SetNameTab("Size"); fillh2b.SetNameTab("Energy"); @@ -354,5 +356,7 @@ // Print the result - *fLog << inf << "Rule: " << rule << endl; + *fLog << inf; + *fLog << "Rule: " << rule << endl; + *fLog << "Disp: " << fTrainParameter << endl; hist.GetAlphaFitter().Print("result"); @@ -366,2 +370,162 @@ } +/* +#include "MParameterCalc.h" +#include "MHillasCalc.h" +#include "../mpointing/MSrcPosRndm.h" + +Bool_t MJTrainDisp::TrainGhostbuster(const char *out, const MDataSet &set, Int_t num) +{ + SetTitle(Form("TrainGhostbuster: %s", out)); + + if (fDisplay) + fDisplay->SetTitle(fTitle); + + if (!set.IsValid()) + { + *fLog << err << "ERROR - DataSet invalid!" << endl; + return kFALSE; + } + + if (!HasWritePermission(out)) + return kFALSE; + + *fLog << inf; + fLog->Separator(GetDescriptor()); + + // --------------------- Setup files -------------------- + MReadMarsFile readtrn("Events"); + MReadMarsFile readtst("Events"); + readtrn.DisableAutoScheme(); + readtst.DisableAutoScheme(); + + if (!set.AddFilesOn(readtrn)) + return kFALSE; + if (!set.AddFilesOff(readtst)) + return kFALSE; + + // ----------------------- Setup Matrix ------------------ + MHMatrix train("Train"); + train.AddColumns(fRules); + if (fEnableWeights) + train.AddColumn("MWeight.fVal"); + train.AddColumn("sign(MHillasSrc.fCosDeltaAlpha)==sign(SignStore.fVal)"); + + MParameterCalc calc("MHillasSrc.fCosDeltaAlpha", "SignStore"); + calc.SetNameParameter("SignStore"); + + MSrcPosRndm rndm; + rndm.SetRule(fTrainParameter); + //rndm.SetDistOfSource(120*3.37e-3); + + MHillasCalc hcalc; + hcalc.SetFlags(MHillasCalc::kCalcHillasSrc); + + // ----------------------- Fill Matrix RF ---------------------- + MTFillMatrix fill(fTitle); + fill.SetDisplay(fDisplay); + fill.SetLogStream(fLog); + fill.SetDestMatrix1(&train, num); + fill.SetReader(&readtrn); + fill.AddPreCuts(fPreCuts); + fill.AddPreCuts(fTrainCuts); + fill.AddPreTasks(fPreTasks); + fill.AddPostTasks(fPostTasks); + fill.AddPostTask(&calc); + fill.AddPostTask(&rndm); + fill.AddPostTask(&hcalc); + if (!fill.Process()) + return kFALSE; + + // ------------------------ Train RF -------------------------- + MRanForestCalc rf("TrainGhostbuster", fTitle); + rf.SetNumTrees(fNumTrees); + rf.SetNdSize(fNdSize); + rf.SetNumTry(fNumTry); + rf.SetNumObsoleteVariables(1); + rf.SetLastDataColumnHasWeights(fEnableWeights); + rf.SetDisplay(fDisplay); + rf.SetLogStream(fLog); + rf.SetFileName(out); + rf.SetDebug(fDebug>1); + rf.SetNameOutput("Sign"); + + if (!rf.TrainRegression(train)) // regression (best choice) + return kFALSE; + + // --------------------- Display result ---------------------- + + gLog.Separator("Test"); + + MH::SetPalette("pretty"); + + MParList plist; + MTaskList tlist; + plist.AddToList(this); + plist.AddToList(&tlist); + //plist.AddToList(&b); + + MFilterList list; + if (!list.AddToList(fPreCuts)) + *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl; + if (!list.AddToList(fTestCuts)) + *fLog << err << "ERROR - Calling MFilterList::AddToList for fTestCuts failed!" << endl; + + MContinue cont(&list); + cont.SetInverted(); + + const char *rule = "abs(Sign.fVal-(sign(MHillasSrc.fCosDeltaAlpha)==sign(SignStore.fVal)))"; + + //MChisqEval eval; + //eval.SetY1("sqrt(ThetaSquared.fVal)"); + + MH3 hsign1("MHillas.fSize", rule); + MH3 hsign2("MMcEvt.fEnergy", rule); + hsign1.SetTitle("Was ist das? vs. Size:Size [phe]:XXX"); + hsign2.SetTitle("Was ist das? vs. Energy:Enerhy [GeV]:XXXX"); + + MBinning binsx( 70, 10, 31623, "BinningMH3X", "log"); + MBinning binsy( 51, 0, 1, "BinningMH3Y", "lin"); + + plist.AddToList(&binsx); + plist.AddToList(&binsy); + + MFillH fillh2a(&hsign1, "", "FillSize"); + MFillH fillh2b(&hsign2, "", "FillEnergy"); + fillh2a.SetDrawOption("profx colz"); + fillh2b.SetDrawOption("profx colz"); + fillh2a.SetNameTab("Size"); + fillh2b.SetNameTab("Energy"); + + tlist.AddToList(&readtst); + tlist.AddToList(&cont); + tlist.AddToList(&calc); + tlist.AddToList(&rndm); + tlist.AddToList(&hcalc); + tlist.AddToList(&rf); + tlist.AddToList(&fillh2a); + tlist.AddToList(&fillh2b); + //tlist.AddToList(&eval); + + MEvtLoop loop(fTitle); + loop.SetLogStream(fLog); + loop.SetDisplay(fDisplay); + loop.SetParList(&plist); + //if (!SetupEnv(loop)) + // return kFALSE; + + if (!loop.Eventloop()) + return kFALSE; + + // Print the result + *fLog << inf << "Rule: " << rule << endl; + + //DisplayResult(hdisp1, hdisp2); + + SetPathOut(out); + if (!WriteDisplay(0, "UPDATE")) + return kFALSE; + + return kTRUE; +} +*/ Index: trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.h =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.h (revision 8655) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.h (revision 8656) @@ -27,4 +27,5 @@ Bool_t Train(const char *out, const MDataSet &set, Int_t num); + //Bool_t TrainGhostbuster(const char *out, const MDataSet &set, Int_t num); ClassDef(MJTrainDisp, 0)//Class to train Random Forest disp estimator