Ignore:
Timestamp:
08/24/07 09:37:31 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mjtrain
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.cc

    r8657 r8704  
    8787#include "MBinning.h"
    8888#include "MH3.h"
     89#include "MHn.h"
    8990#include "MHThetaSq.h"
    9091
    … …  
    115116    line.SetLineWidth(1);
    116117
    117     c.cd(i);
     118    c.cd(i+4);
     119    gPad->SetBorderMode(0);
     120    gPad->SetFrameBorderMode(0);
     121    //gPad->SetFillColor(kWhite);
     122    gPad->SetLogx();
     123    gPad->SetGridx();
     124    gPad->SetGridy();
     125    //gPad->SetLeftMargin(0.12);
     126    //gPad->SetRightMargin(0.12);
     127
     128    const Float_t cutval = hist.GetYaxis()->GetBinLowEdge(4);
     129
     130    TH1D heff;
     131    heff.SetName(Form("Eff%s", hist.GetName()));
     132    heff.SetTitle(Form("Cut efficiency vs. %s for \\vartheta<%.3f", hist.GetName(), TMath::Sqrt(cutval)));
     133    heff.SetDirectory(0);
     134    heff.SetXTitle(hist.GetXaxis()->GetTitle());
     135    heff.SetYTitle("Efficiency");
     136
     137    MH::SetBinning(&heff, hist.GetXaxis());
     138
     139
     140    for (int x=0; x<=hist.GetNbinsX()+1; x++)
     141    {
     142        const Double_t n0 = hist.Integral(x, x, -1, -1);
     143        if (n0>0)
     144            heff.SetBinContent(x, hist.Integral(x, x, -1, 3)/n0);
     145    }
     146
     147    heff.SetMinimum(0);
     148    heff.SetMaximum(1);
     149    heff.DrawCopy();
     150
     151    line.DrawLine(10, 0.5, 31623, 0.5);
     152
     153    c.cd(i+0);
    118154    gPad->SetBorderMode(0);
    119155    gPad->SetFrameBorderMode(0);
    … …  
    179215{
    180216    TCanvas &c = fDisplay->AddTab("Disp");
    181     c.Divide(2,2);
     217    c.Divide(2,3);
    182218
    183219    DisplayHist(c, 1, hsize);
    … …  
    215251
    216252    if (!set.AddFilesOn(readtrn))
    217         return kFALSE;
     253    {
     254        *fLog << err << "ERROR - Adding SequencesOn." << endl;
     255        return kFALSE;
     256    }
    218257    if (!set.AddFilesOff(readtst))
    219         return kFALSE;
     258    {
     259        *fLog << err << "ERROR - Adding SequencesOff." << endl;
     260        return kFALSE;
     261    }
    220262
    221263    // ----------------------- Setup Matrix ------------------
    … …  
    277319
    278320    MParameterD par("ThetaSquaredCut");
    279     par.SetVal(0.2);
     321    par.SetVal(0.215*0.215);
    280322    plist.AddToList(&par);
    281323
    … …  
    297339    cont.SetInverted();
    298340
     341    const char *rule = "(MHillasSrc.fDist*MGeomCam.fConvMm2Deg)^2 + (Disp.fVal)^2 - (2*MHillasSrc.fDist*MGeomCam.fConvMm2Deg*Disp.fVal*cos(MHillasSrc.fAlpha*kDeg2Rad))";
     342
     343    MParameterCalc calcthetasq(rule, "MThetaSqCalc");
     344    calcthetasq.SetNameParameter("ThetaSquared");
     345
     346    MChisqEval eval;
     347    eval.SetY1("sqrt(ThetaSquared.fVal)");
     348
     349    // ----------- Setup binnings ----------------
     350    MBinning binsS(50,  10,    100000, "BinningSize",         "log");
     351    MBinning binsE(70,  10,    31623,  "BinningEnergy",       "log");
     352    MBinning binsG(50, -10,    10,     "BinningSlope",        "lin");
     353    MBinning binsX(50, -1,     1,      "BinningResidualDist", "lin");
     354    MBinning binsL(50,  0,     0.3,    "BinningLeakage",      "lin");
     355    MBinning binsT(51, -0.005, 0.505,  "BinningTheta",        "asin");
     356    MBinning binsC(50,  1e-2,  1,      "BinningConc",         "log");
     357    MBinning binsW(50,  0,     0.5,    "BinningLength",       "lin");
     358    MBinning binsM(50,  0,     0.3,    "BinningWidth",        "lin");
     359    MBinning binsV(75,  0, par.GetVal()*25, "BinningThetaSq", "lin");
     360
     361    plist.AddToList(&binsG);
     362    plist.AddToList(&binsS);
     363    plist.AddToList(&binsX);
     364    plist.AddToList(&binsE);
     365    plist.AddToList(&binsL);
     366    plist.AddToList(&binsT);
     367    plist.AddToList(&binsC);
     368    plist.AddToList(&binsV);
     369    plist.AddToList(&binsW);
     370    plist.AddToList(&binsM);
     371
     372    // ----------- Setup some histograms ----------------
     373
    299374    MHThetaSq hist;
    300375    hist.SkipHistTime();
    … …  
    302377    hist.SkipHistEnergy();
    303378
    304     MFillH fillh(&hist, "", "FillThetaSq");
    305 
    306     // 0 =  disp^2 - 2*disp*dist*cos(alpha) + dist^2
    307 
    308     // cos^2 -1 = - sin^2
    309 
    310     // disp = +dist* (cos(alpha) +/- sqrt(cos^2(alpha) - 1) )
    311 
    312     const char *rule = "(MHillasSrc.fDist*MGeomCam.fConvMm2Deg)^2 + (Disp.fVal)^2 - (2*MHillasSrc.fDist*MGeomCam.fConvMm2Deg*Disp.fVal*cos(MHillasSrc.fAlpha*kDeg2Rad))";
    313 
    314     MParameterCalc calcthetasq(rule, "MThetaSqCalc");
    315     calcthetasq.SetNameParameter("ThetaSquared");
    316 
    317     MChisqEval eval;
    318     eval.SetY1("sqrt(ThetaSquared.fVal)");
     379    // To speed it up we could precalculate it.
     380    const char *res = "Disp.fVal-MHillasSrc.fDist*3.37e-3";
     381
     382    MHn hres1("Disp1", "Xi Residual (Dist/Disp)");
     383    hres1.AddHist("MHillas.fSize", res);
     384    hres1.InitName("ResSize;Size;ResidualDist");
     385    hres1.InitTitle(";S [phe];Disp-Dist [\\circ];");
     386    hres1.SetDrawOption("colz profx");
     387    hres1.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3", res);
     388    hres1.InitName("ResSlope;Slope;ResidualDist");
     389    hres1.InitTitle(";Slope;Disp-Dist [\\circ];");
     390    hres1.SetDrawOption("colz profx");
     391    hres1.AddHist("MNewImagePar.fLeakage1", res);
     392    hres1.InitName("ResLeak;Leakage;ResidualDist");
     393    hres1.InitTitle(";Leak;Disp-Dist [\\circ];");
     394    hres1.SetDrawOption("colz profx");
     395    hres1.AddHist("MPointingPos.fZd", res);
     396    hres1.InitName("ResTheta;Theta;ResidualDist");
     397    hres1.InitTitle(";Zd [\\circ];Disp-Dist [\\circ];");
     398    hres1.SetDrawOption("colz profx");
     399
     400    MHn hres2("Disp2", "Dist Residual (Disp-Dist)");
     401    hres2.AddHist("MHillas.fLength*3.37e-3", res);
     402    hres2.InitName("ResLength;Length;ResidualDist");
     403    hres2.InitTitle(";L [\\circ];Disp-Dist [\\circ];");
     404    hres2.SetDrawOption("colz profx");
     405    hres2.AddHist("MNewImagePar.fConc1", res);
     406    hres2.InitName("ResConc1;Conc;ResidualDist");
     407    hres2.InitTitle(";C;Disp-Dist [\\circ];");
     408    hres2.SetDrawOption("colz profx");
     409    hres2.AddHist("MHillas.fWidth*3.37e-3", res);
     410    hres2.InitName("ResWidth;Width;ResidualDist");
     411    hres2.InitTitle(";W [\\circ];Disp-Dist [\\circ];");
     412    hres2.SetDrawOption("colz profx");
     413    hres2.AddHist("MMcEvt.fEnergy", res);
     414    hres2.InitName("ResEmc;Energy;ResidualDist");
     415    hres2.InitTitle(";E_{mc} [GeV];Disp-Dist [\\circ];");
     416    hres2.SetDrawOption("colz profx");
    319417
    320418    MH3 hdisp1("MHillas.fSize",  "ThetaSquared.fVal");
    321419    MH3 hdisp2("MMcEvt.fEnergy", "ThetaSquared.fVal");
    322     hdisp1.SetTitle("\\vartheta distribution vs. Size:Size [phe]:\\vartheta [\\circ]");
    323     hdisp2.SetTitle("\\vartheta distribution vs. Energy:Enerhy [GeV]:\\vartheta [\\circ]");
    324 
    325     MBinning binsx( 70, 10, 31623, "BinningMH3X", "log");
    326     MBinning binsy( 75, 0,  0.3,   "BinningMH3Y", "lin");
    327 
    328     plist.AddToList(&binsx);
    329     plist.AddToList(&binsy);
    330 
    331     MFillH fillh2a(&hdisp1, "", "FillSize");
    332     MFillH fillh2b(&hdisp2, "", "FillEnergy");
    333     fillh2a.SetDrawOption("colz profx");
    334     fillh2b.SetDrawOption("colz profx");
    335     fillh2a.SetNameTab("Size");
    336     fillh2b.SetNameTab("Energy");
     420    hdisp1.SetName("Size;Size;ThetaSq");
     421    hdisp2.SetName("Energy;Energy;ThetaSq");
     422    hdisp1.SetTitle("\\vartheta distribution vs. Size:Size [phe]:\\vartheta^2 [\\circ]");
     423    hdisp2.SetTitle("\\vartheta distribution vs. Energy:Energy [GeV]:\\vartheta^2 [\\circ]");
     424
     425    // -------------- Setup fill tasks ----------------
     426
     427    MFillH fillh(&hist, "", "FillThetaSq");
     428    MFillH fillh2a(&hres1, "", "FillResiduals1");
     429    MFillH fillh2b(&hres2, "", "FillResiduals2");
     430    MFillH fillh2c(&hdisp1, "", "FillSize");
     431    MFillH fillh2d(&hdisp2, "", "FillEnergy");
     432    fillh2c.SetBit(MFillH::kDoNotDisplay);
     433    fillh2d.SetBit(MFillH::kDoNotDisplay);
     434
     435    // --------------- Setup weighting -------------------
     436
     437    if (fEnableWeights)
     438    {
     439        fillh.SetWeight();
     440        fillh2a.SetWeight();
     441        fillh2b.SetWeight();
     442        fillh2c.SetWeight();
     443        fillh2d.SetWeight();
     444        eval.SetNameWeight();
     445    }
     446
     447    // --------------- Setup tasklist -------------------
    337448
    338449    tlist.AddToList(&readtst);
     450    tlist.AddToList(fPreTasks);
    339451    tlist.AddToList(&cont);
    340452    tlist.AddToList(&rf);
    341453    tlist.AddToList(&calcthetasq);
     454    tlist.AddToList(fPostTasks);
    342455    tlist.AddToList(&fillh);
    343456    tlist.AddToList(&fillh2a);
    344457    tlist.AddToList(&fillh2b);
     458    tlist.AddToList(&fillh2c);
     459    tlist.AddToList(&fillh2d);
     460    tlist.AddToList(fTestTasks);
    345461    tlist.AddToList(&eval);
     462
     463    // ------------- Setup/run eventloop -----------------
    346464
    347465    MEvtLoop loop(fTitle);
    … …  
    355473        return kFALSE;
    356474
     475    // ---------------- Prepare result -------------------
     476
    357477    // Print the result
    358478    *fLog << inf;
    … …  
    361481    hist.GetAlphaFitter().Print("result");
    362482
     483    // The user has closed the display
     484    if (!fDisplay)
     485        return kTRUE;
     486
    363487    DisplayResult(hdisp1, hdisp2);
    364488
    365489    SetPathOut(out);
    366     if (!WriteDisplay(0, "UPDATE"))
    367         return kFALSE;
    368 
    369     return kTRUE;
     490    return WriteDisplay(0, "UPDATE");
    370491}
    371492

  • trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc

    r8646 r8704  
    7676
    7777// histograms
     78#include "MHn.h"
     79#include "MBinning.h"
    7880#include "MHEnergyEst.h"
    7981
    … …  
    128130    train.AddColumn("MMcEvt.fImpact/100");
    129131    train.AddColumn("MMcEvt.fTelescopeTheta*TMath::RadToDeg()");
    130     train.AddColumn("MMcEvt.fEnergy");
     132    train.AddColumn(fTrainParameter);
    131133
    132134
    … …  
    156158    rf.SetDebug(fDebug>1);
    157159    rf.SetNameOutput("MEnergyEst");
     160    rf.SetFunction(fResultFunction);
    158161
    159162    /*
    … …  
    188191    cont.SetInverted();
    189192
     193    // -------------------------------------------------------------
     194    MBinning binsS(50, 10,     100000, "BinningSize",           "log");
     195    MBinning binsE(70, 10,     31623,  "BinningEnergy",         "log");
     196    MBinning binsG(50,-10,     10,     "BinningSlope",          "lin");
     197    MBinning binsR(50, -1,     1,      "BinningEnergyResidual", "lin");
     198    MBinning binsL(50,  0,     0.3,    "BinningLeakage",        "lin");
     199    MBinning binsT(51, -0.005, 0.505,  "BinningTheta",          "asin");
     200    MBinning binsD(50,  0,     1.6,    "BinningDist",           "lin");
     201    MBinning binsC(50,  1e-2,  1,      "BinningConc",           "log");
     202
     203    plist.AddToList(&binsG);
     204    plist.AddToList(&binsS);
     205    plist.AddToList(&binsR);
     206    plist.AddToList(&binsE);
     207    plist.AddToList(&binsL);
     208    plist.AddToList(&binsT);
     209    plist.AddToList(&binsD);
     210    plist.AddToList(&binsC);
     211
    190212    MHEnergyEst hist;
     213
     214    // To speed it up we could precalculate it.
     215    const char *res = "log10(MEnergyEst.fVal)-log10(MMcEvt.fEnergy)";
     216
     217    MHn hres1("Energy1", "Energy Residual (lg E_{est} - lg E_{mc})");
     218    hres1.AddHist("MHillas.fSize", res);
     219    hres1.InitName("ResSize;Size;EnergyResidual");
     220    hres1.InitTitle(";S [phe];\\Delta lg E;");
     221    hres1.SetDrawOption("colz profx");
     222    hres1.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3", res);
     223    hres1.InitName("ResSlope;Slope;EnergyResidual");
     224    hres1.InitTitle(";Slope;\\Delta lg E;");
     225    hres1.SetDrawOption("colz profx");
     226    hres1.AddHist("MNewImagePar.fLeakage1", res);
     227    hres1.InitName("ResLeak;Leakage;EnergyResidual");
     228    hres1.InitTitle(";Leak;\\Delta lg E;");
     229    hres1.SetDrawOption("colz profx");
     230    hres1.AddHist("MHillasSrc.fDist*3.37e-3", res);
     231    hres1.InitName("ResDist;Dist;EnergyResidual");
     232    hres1.InitTitle(";D [\\circ];\\Delta lg E;");
     233    hres1.SetDrawOption("colz profx");
     234
     235    MHn hres2("Energy2", "Energy Residual (lg E_{est} - lg E_{mc})");
     236    hres2.AddHist("MMcEvt.fEnergy", res);
     237    hres2.InitName("ResEmc;Energy;EnergyResidual");
     238    hres2.InitTitle(";E_{mc} [GeV];\\Delta lg E;");
     239    hres2.SetDrawOption("colz profx");
     240    hres2.SetAutoRange(kFALSE, kFALSE, kFALSE);
     241    hres2.AddHist("MPointingPos.fZd", res);
     242    hres2.InitName("ResTheta;Theta;EnergyResidual");
     243    hres2.InitTitle(";Zd [\\circ];\\Delta lg E;");
     244    hres2.SetDrawOption("colz profx");
     245    hres2.AddHist("MEnergyEst.fVal", res);
     246    hres2.InitName("ResEest;Energy;EnergyResidual");
     247    hres2.InitTitle(";E_{est} [GeV];\\Delta lg E;");
     248    hres2.SetDrawOption("colz profx");
     249    hres2.SetAutoRange(kFALSE, kFALSE, kFALSE);
     250    hres2.AddHist("MNewImagePar.fConc1", res);
     251    hres2.InitName("ResConc1;Conc;EnergyResidual");
     252    hres2.InitTitle(";C;\\Delta lg E;");
     253    hres2.SetDrawOption("colz profx");
     254
    191255    MFillH fillh(&hist);
     256    MFillH fillh1(&hres1, "", "FillResiduals1");
     257    MFillH fillh2(&hres2, "", "FillResiduals2");
     258
    192259    if (fEnableWeights)
     260    {
    193261        fillh.SetWeight();
     262        fillh1.SetWeight();
     263        fillh2.SetWeight();
     264    }
    194265
    195266    tlist.AddToList(&readtst);
    … …  
    199270    tlist.AddToList(fPostTasks);
    200271    tlist.AddToList(&fillh);
     272    tlist.AddToList(&fillh1);
     273    tlist.AddToList(&fillh2);
     274    tlist.AddToList(fTestTasks);
    201275
    202276    MEvtLoop loop(fTitle);

  • trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.h

    r7412 r8704  
    1010class MJTrainEnergy : public MJTrainRanForest
    1111{
     12private:
     13    TString fTrainParameter;
     14    TString fResultFunction;
     15
    1216public:
    13     MJTrainEnergy() { }
     17    MJTrainEnergy() { SetTrainLin(); }
     18
     19    void SetTrainLog() { SetTrainFunc("log10(MMcEvt.fEnergy)", "pow(10, x)"); }
     20    void SetTrainLin() { SetTrainFunc("MMcEvt.fEnergy", "x"); }
     21
     22    void SetTrainFunc(const char *par, const char *res)
     23    {
     24        fTrainParameter = par;
     25        fResultFunction = res;
     26    }
     27
    1428    Bool_t Train(const char *out, const MDataSet &set, Int_t num);
    1529
Note: See TracChangeset for help on using the changeset viewer.