| 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): Thomas Bretz, 6/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
|
|---|
| 19 | !
|
|---|
| 20 | ! Copyright: MAGIC Software Development, 2000-2005
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 25 | //
|
|---|
| 26 | // MJOptimizeDisp
|
|---|
| 27 | //
|
|---|
| 28 | // Class for otimizing the disp parameters. For more details see
|
|---|
| 29 | // MJOptimize.
|
|---|
| 30 | //
|
|---|
| 31 | // MJOptimizeDisp minimzes the width of the Theta^2 distribution by
|
|---|
| 32 | // minimizing <Theta^2> (the RMS^2 of Theta^2), while Theta^2 is
|
|---|
| 33 | // calculated as:
|
|---|
| 34 | // Theta^2 = d^2 + p^2 - 2*d*p*cos(a)
|
|---|
| 35 | // with:
|
|---|
| 36 | // d: MHillasSrc.fDist [deg]
|
|---|
| 37 | // p: Disp as calculated by the given rule, eg: [0] (a constant)
|
|---|
| 38 | // a: MHillasSrc.fAlpha [rad]
|
|---|
| 39 | //
|
|---|
| 40 | // Example:
|
|---|
| 41 | // --------
|
|---|
| 42 | // MJOptimizeDisp opt;
|
|---|
| 43 | // opt.SetDebug(2);
|
|---|
| 44 | // opt.SetOptimizer(MJOptimize::kMigrad);
|
|---|
| 45 | // opt.SetNumEvents(20000);
|
|---|
| 46 | // opt.EnableTestTrain();
|
|---|
| 47 | // opt.AddParameter("1-(MHillas.fWidth/MHillas.fLength)");
|
|---|
| 48 | // opt.SetParameter(0, 1, 0, 2);
|
|---|
| 49 | // char *r = "[0]*M[0]"; //Rule to calculate disp
|
|---|
| 50 | // MStatusDisplay *d = new MStatusDisplay;
|
|---|
| 51 | // opt.SetDisplay(d);
|
|---|
| 52 | // opt.RunDisp("ganymed-summary.root", r);
|
|---|
| 53 | //
|
|---|
| 54 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 55 | #include "MJOptimizeDisp.h"
|
|---|
| 56 |
|
|---|
| 57 | #include "MHMatrix.h"
|
|---|
| 58 |
|
|---|
| 59 | // environment
|
|---|
| 60 | #include "MLog.h"
|
|---|
| 61 | #include "MLogManip.h"
|
|---|
| 62 | #include "MStatusDisplay.h"
|
|---|
| 63 |
|
|---|
| 64 | // eventloop
|
|---|
| 65 | #include "MParList.h"
|
|---|
| 66 | #include "MTaskList.h"
|
|---|
| 67 | #include "MEvtLoop.h"
|
|---|
| 68 |
|
|---|
| 69 | // parameters
|
|---|
| 70 | #include "MParameters.h"
|
|---|
| 71 | #include "MGeomCamFACT.h"
|
|---|
| 72 |
|
|---|
| 73 | // histograms
|
|---|
| 74 | #include "MH3.h"
|
|---|
| 75 | #include "MBinning.h"
|
|---|
| 76 | #include "../mhflux/MAlphaFitter.h"
|
|---|
| 77 | #include "../mhflux/MHThetaSq.h"
|
|---|
| 78 | #include "../mtools/MChisqEval.h"
|
|---|
| 79 |
|
|---|
| 80 | // tasks
|
|---|
| 81 | #include "MReadTree.h"
|
|---|
| 82 | #include "MMatrixLoop.h"
|
|---|
| 83 | #include "MParameterCalc.h"
|
|---|
| 84 | #include "MFillH.h"
|
|---|
| 85 | #include "MContinue.h"
|
|---|
| 86 | #include "MWriteRootFile.h"
|
|---|
| 87 |
|
|---|
| 88 | // filters
|
|---|
| 89 | #include "MFilterList.h"
|
|---|
| 90 | #include "MFDataMember.h"
|
|---|
| 91 |
|
|---|
| 92 | using namespace std;
|
|---|
| 93 |
|
|---|
| 94 | ClassImp(MJOptimizeDisp);
|
|---|
| 95 |
|
|---|
| 96 | //------------------------------------------------------------------------
|
|---|
| 97 | //
|
|---|
| 98 | // Read all events from file which do match rules and optimize
|
|---|
| 99 | // disp estimator.
|
|---|
| 100 | //
|
|---|
| 101 | Bool_t MJOptimizeDisp::RunDisp(const char *fname, const char *rule, MTask *weights)
|
|---|
| 102 | {
|
|---|
| 103 | SetTitle(Form("OptimizeDisp: %s", fname));
|
|---|
| 104 |
|
|---|
| 105 | if (fDisplay)
|
|---|
| 106 | fDisplay->SetTitle(fTitle);
|
|---|
| 107 |
|
|---|
| 108 | fLog->Separator("Preparing Disp optimization");
|
|---|
| 109 |
|
|---|
| 110 | MParList parlist;
|
|---|
| 111 |
|
|---|
| 112 | MParameterI par1("DataType");
|
|---|
| 113 | par1.SetVal(1);
|
|---|
| 114 | parlist.AddToList(&par1);
|
|---|
| 115 |
|
|---|
| 116 | MParameterD par2("ThetaSquaredCut");
|
|---|
| 117 | par2.SetVal(0.2);
|
|---|
| 118 | parlist.AddToList(&par2);
|
|---|
| 119 |
|
|---|
| 120 | MAlphaFitter fit;
|
|---|
| 121 | fit.SetPolynomOrder(0);
|
|---|
| 122 | fit.SetSignalFitMax(0.8);
|
|---|
| 123 | fit.EnableBackgroundFit(kFALSE);
|
|---|
| 124 | fit.SetSignalFunction(MAlphaFitter::kThetaSq);
|
|---|
| 125 | fit.SetMinimizationStrategy(MAlphaFitter::kGaussSigma);
|
|---|
| 126 | parlist.AddToList(&fit);
|
|---|
| 127 |
|
|---|
| 128 | // To avoid this hard-coded we have to switch to MReadMarsFile
|
|---|
| 129 | MGeomCamFACT geom; // For GetConvMm2Deg
|
|---|
| 130 | parlist.AddToList(&geom);
|
|---|
| 131 |
|
|---|
| 132 | MHMatrix m("M");
|
|---|
| 133 | AddRulesToMatrix(m);
|
|---|
| 134 | parlist.AddToList(&m);
|
|---|
| 135 |
|
|---|
| 136 | const Int_t num1 = m.AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
|
|---|
| 137 | const Int_t num2 = m.AddColumn("MHillasSrc.fAlpha*TMath::DegToRad()");
|
|---|
| 138 | const Int_t num3 = m.AddColumn("MHillas.fSize");
|
|---|
| 139 |
|
|---|
| 140 | MHThetaSq hist;
|
|---|
| 141 | hist.SkipHistTime();
|
|---|
| 142 | hist.SkipHistTheta();
|
|---|
| 143 | //hist.SkipHistEnergy();
|
|---|
| 144 | //hist.ForceUsingSize();
|
|---|
| 145 | hist.InitMapping(&m, 1);
|
|---|
| 146 |
|
|---|
| 147 | MFDataMember filter("DataType.fVal", '>', 0.5);
|
|---|
| 148 | fPreCuts.Add(&filter);
|
|---|
| 149 |
|
|---|
| 150 | MParameterCalc calc1(rule, "MParameters");
|
|---|
| 151 | calc1.SetNameParameter("Disp");
|
|---|
| 152 | parlist.AddToList(&calc1);
|
|---|
| 153 |
|
|---|
| 154 | const char *disp = "Disp.fVal";
|
|---|
| 155 | const char *n1 = Form("M[%d]", num1);
|
|---|
| 156 |
|
|---|
| 157 | const char *rule2 = Form("(%s*%s) + (%s*%s) - (2*%s*%s*cos(M[%d]))",
|
|---|
| 158 | n1, n1, disp, disp, n1, disp, num2);
|
|---|
| 159 |
|
|---|
| 160 | MParameterCalc calc2(rule2, "MThetaSqCalc");
|
|---|
| 161 | calc2.SetNameParameter("ThetaSquared");
|
|---|
| 162 |
|
|---|
| 163 | MReadTree read("Events");
|
|---|
| 164 | // NECESSARY BECAUSE OF MDataFormula GetRules missing
|
|---|
| 165 | read.DisableAutoScheme();
|
|---|
| 166 | if (fname)
|
|---|
| 167 | read.AddFile(fname);
|
|---|
| 168 | else
|
|---|
| 169 | AddSequences(read, fNamesOn);
|
|---|
| 170 | if (!FillMatrix(read, parlist, kTRUE))
|
|---|
| 171 | return kFALSE;
|
|---|
| 172 |
|
|---|
| 173 | fPreCuts.Remove(&filter);
|
|---|
| 174 |
|
|---|
| 175 | MTaskList tasklist;
|
|---|
| 176 | parlist.Replace(&tasklist);
|
|---|
| 177 |
|
|---|
| 178 | MFillH fill(&hist);
|
|---|
| 179 | if (weights)
|
|---|
| 180 | fill.SetWeight();
|
|---|
| 181 |
|
|---|
| 182 | MChisqEval eval;
|
|---|
| 183 | eval.SetY1(fUseThetaSq?"ThetaSquared.fVal":"sqrt(ThetaSquared.fVal)");
|
|---|
| 184 | if (weights)
|
|---|
| 185 | eval.SetNameWeight();
|
|---|
| 186 |
|
|---|
| 187 | MMatrixLoop loop(&m);
|
|---|
| 188 |
|
|---|
| 189 | const char *n3 = Form("M[%d]", num3);
|
|---|
| 190 | MH3 hdisp(n3, "sqrt(ThetaSquared.fVal)");
|
|---|
| 191 | hdisp.SetTitle("\\vartheta^{2} distribution vs. Size:Size [phe]:\\vartheta^{2} [\\circ^{2}]");
|
|---|
| 192 |
|
|---|
| 193 | MBinning binsx(100, 10, 100000, "BinningMH3X", "log");
|
|---|
| 194 | MBinning binsy(100, 0, 2, "BinningMH3Y", "lin");
|
|---|
| 195 |
|
|---|
| 196 | parlist.AddToList(&binsx);
|
|---|
| 197 | parlist.AddToList(&binsy);
|
|---|
| 198 |
|
|---|
| 199 | MFillH fillh2(&hdisp);
|
|---|
| 200 | fillh2.SetDrawOption("blue profx");
|
|---|
| 201 |
|
|---|
| 202 | tasklist.AddToList(&loop);
|
|---|
| 203 | tasklist.AddToList(&calc1);
|
|---|
| 204 | tasklist.AddToList(&calc2);
|
|---|
| 205 | if (weights)
|
|---|
| 206 | tasklist.AddToList(weights);
|
|---|
| 207 | //tasklist.AddToList(&fill);
|
|---|
| 208 | //tasklist.AddToList(&fillh2);
|
|---|
| 209 | tasklist.AddToList(&eval);
|
|---|
| 210 |
|
|---|
| 211 | // Optimize with the tasklist in this parameterlist
|
|---|
| 212 | if (!Optimize(parlist))
|
|---|
| 213 | return kFALSE;
|
|---|
| 214 |
|
|---|
| 215 | // Print the result
|
|---|
| 216 | *fLog << inf;
|
|---|
| 217 | *fLog << "Finished processing of " << fname << endl;
|
|---|
| 218 | *fLog << "With Rule: " << rule << endl;
|
|---|
| 219 | //hist.GetAlphaFitter().Print("result");
|
|---|
| 220 |
|
|---|
| 221 | if (fNameOut.IsNull())
|
|---|
| 222 | return kTRUE;
|
|---|
| 223 |
|
|---|
| 224 | fLog->Separator("Writing result");
|
|---|
| 225 |
|
|---|
| 226 | MTaskList tlist;
|
|---|
| 227 | MParList plist;
|
|---|
| 228 | plist.AddToList(&tlist);
|
|---|
| 229 | plist.AddToList(&geom);
|
|---|
| 230 |
|
|---|
| 231 | MReadTree read2("Events");
|
|---|
| 232 | read2.DisableAutoScheme();
|
|---|
| 233 | read2.AddFiles(read);
|
|---|
| 234 | tlist.AddToList(&read2);
|
|---|
| 235 |
|
|---|
| 236 | MFilterList flist;
|
|---|
| 237 | fPreCuts.Add(&filter);
|
|---|
| 238 | flist.AddToList(fPreCuts);
|
|---|
| 239 | flist.SetName("PreCuts");
|
|---|
| 240 | flist.SetInverted(kTRUE);
|
|---|
| 241 |
|
|---|
| 242 | MContinue cont(&flist);
|
|---|
| 243 | tlist.AddToList(&cont);
|
|---|
| 244 |
|
|---|
| 245 | MParameterCalc calcx("0", "CountEvents");
|
|---|
| 246 | tlist.AddToList(&calcx);
|
|---|
| 247 |
|
|---|
| 248 | // FIXME: Handle fTestTrain<0
|
|---|
| 249 | MContinue contt("TMath::Odd(CountEvents.fNumExecutions)>0.5", "ContTrain");
|
|---|
| 250 | plist.AddToList(&calcx);
|
|---|
| 251 | tlist.AddToList(&contt);
|
|---|
| 252 |
|
|---|
| 253 | TString r(rule);
|
|---|
| 254 |
|
|---|
| 255 | TIter Next(&fRules);
|
|---|
| 256 | TObject *o=0;
|
|---|
| 257 | Int_t num = 0;
|
|---|
| 258 | while ((o=Next()))
|
|---|
| 259 | {
|
|---|
| 260 | r.ReplaceAll(Form("M[%d]", num), Form("M%d.fVal", num));
|
|---|
| 261 |
|
|---|
| 262 | MParameterCalc *c = new MParameterCalc(o->GetName(), Form("CalcM%d", num));
|
|---|
| 263 | c->SetNameParameter(Form("M%d", num++));
|
|---|
| 264 | c->SetBit(kCanDelete);
|
|---|
| 265 | tlist.AddToList(c);
|
|---|
| 266 | }
|
|---|
| 267 |
|
|---|
| 268 | MParameterCalc calcr(r, "CalcDisp");
|
|---|
| 269 | calcr.SetVariables(GetParameters());
|
|---|
| 270 | calcr.SetNameParameter("Disp");
|
|---|
| 271 |
|
|---|
| 272 | MParameterCalc calct("MHillasSrc.fDist^2*MGeomCam.fConvMm2Deg^2 + Disp.fVal^2 - 2*MHillasSrc.fDist*MGeomCam.fConvMm2Deg*Disp.fVal*cos(MHillasSrc.fAlpha*TMath::DegToRad())", "CalcTheta");
|
|---|
| 273 | calct.SetNameParameter("ThetaSquared");
|
|---|
| 274 |
|
|---|
| 275 | tlist.AddToList(&calcr);
|
|---|
| 276 | tlist.AddToList(&calct);
|
|---|
| 277 |
|
|---|
| 278 | tlist.AddToList(fTestTasks);
|
|---|
| 279 |
|
|---|
| 280 | MWriteRootFile write(fNameOut);
|
|---|
| 281 | write.AddContainer("MHillas", "Events");
|
|---|
| 282 | write.AddContainer("MHillasSrc", "Events");
|
|---|
| 283 | write.AddContainer("MHillasExt", "Events");
|
|---|
| 284 | write.AddContainer("MImagePar", "Events");
|
|---|
| 285 | write.AddContainer("MNewImagePar", "Events");
|
|---|
| 286 | write.AddContainer("Disp", "Events");
|
|---|
| 287 | write.AddContainer("ThetaSquared", "Events");
|
|---|
| 288 | write.AddContainer("MMcEvt", "Events");
|
|---|
| 289 | write.AddContainer("DataType", "Events");
|
|---|
| 290 | write.AddContainer("Weight", "Events");
|
|---|
| 291 | write.AddContainer("FileId", "Events");
|
|---|
| 292 | write.AddContainer("EvtNumber", "Events");
|
|---|
| 293 | tlist.AddToList(&write);
|
|---|
| 294 |
|
|---|
| 295 | MEvtLoop loop2;
|
|---|
| 296 | loop2.SetDisplay(fDisplay);
|
|---|
| 297 | loop2.SetParList(&plist);
|
|---|
| 298 | if (!loop2.Eventloop())
|
|---|
| 299 | return kFALSE;
|
|---|
| 300 |
|
|---|
| 301 | // Store result if requested
|
|---|
| 302 | TObjArray arr;
|
|---|
| 303 | if (fDisplay)
|
|---|
| 304 | arr.Add(fDisplay);
|
|---|
| 305 | arr.Add(&calc1);
|
|---|
| 306 | return WriteContainer(arr, fNameOut);
|
|---|
| 307 | }
|
|---|
