Ignore:
Timestamp:
06/11/08 12:32:36 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhflux
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

    r8907 r8935  
    2020!   Author(s): Thomas Bretz 1/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
    2121!
    22 !   Copyright: MAGIC Software Development, 2000-2005
     22!   Copyright: MAGIC Software Development, 2000-2008
    2323!
    2424!
     
    3131//  calculates the migration matrix E-est vs. E-true
    3232//  for different bins in Theta
     33//
     34//  Class Version 2:
     35//  - fHResolution
     36//  + fHResolutionEst
     37//  + fHResolutionMC
    3338//
    3439//////////////////////////////////////////////////////////////////////////////
     
    8287    fHEnergy.SetZTitle("\\Theta [\\circ]");
    8388
    84     fHResolution.SetDirectory(NULL);
    85     fHResolution.SetName("EnergyRes");
    86     fHResolution.SetTitle("Histogram in \\Delta E/E vs E_{est} and E_{mc}");
    87     fHResolution.SetXTitle("E_{est} [GeV]");
    88     fHResolution.SetYTitle("E_{mc} [GeV]");
    89     fHResolution.SetZTitle("E_{est}/E_{mc}-1");
     89    fHResolutionEst.SetDirectory(NULL);
     90    fHResolutionEst.SetName("ResEnergyEst");
     91    fHResolutionEst.SetTitle("Histogram in \\Delta E/E_{est} vs E_{est}");
     92    fHResolutionEst.SetXTitle("E_{est} [GeV]");
     93    fHResolutionEst.SetYTitle("1-E_{mc}/E_{est}");
     94
     95    fHResolutionMC.SetDirectory(NULL);
     96    fHResolutionMC.SetName("ResEnergyMC");
     97    fHResolutionMC.SetTitle("Histogram in \\Delta E/E_{mc} vs E_{mc}");
     98    fHResolutionMC.SetXTitle("E_{mc} [GeV]");
     99    fHResolutionMC.SetYTitle("E_{est}/E_{mc}-1");
    90100
    91101    fHImpact.SetDirectory(NULL);
     
    105115    SetBinning(&fHEnergy,     &binse, &binse, &binst);
    106116    SetBinning(&fHImpact,     &binsi, &binsr);
    107     SetBinning(&fHResolution, &binse, &binse, &binsr);
     117
     118    SetBinning(&fHResolutionEst, &binse, &binsr);
     119    SetBinning(&fHResolutionMC,  &binse, &binsr);
    108120
    109121    // For some unknown reasons this must be called after
    110122    // the binning has been initialized at least once
    111123    fHEnergy.Sumw2();
    112     fHResolution.Sumw2();
    113124    fHImpact.Sumw2();
     125    fHResolutionEst.Sumw2();
     126    fHResolutionMC.Sumw2();
    114127}
    115128
     
    154167    SetBinning(&fHEnergy,     &binse, &binse, &binst);
    155168    SetBinning(&fHImpact,     &binsi, &binsr);
    156     SetBinning(&fHResolution, &binse, &binse, &binsr);
     169
     170    SetBinning(&fHResolutionEst, &binse, &binsr);
     171    SetBinning(&fHResolutionMC,  &binse, &binsr);
    157172
    158173    fChisq = 0;
     
    161176    fHEnergy.Reset();
    162177    fHImpact.Reset();
    163     fHResolution.Reset();
     178
     179    fHResolutionEst.Reset();
     180    fHResolutionMC.Reset();
    164181
    165182    return kTRUE;
     
    176193    const Double_t imp   = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
    177194    const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
    178     const Double_t resE  = (eest-etru)/etru;
     195
     196    const Double_t resEst  = (eest-etru)/eest;
     197    const Double_t resMC   = (eest-etru)/etru;
    179198
    180199    fHEnergy.Fill(eest, etru, theta, w);
    181     fHResolution.Fill(eest, etru, resE, w);
    182     fHImpact.Fill(imp, resE, w);
     200    fHImpact.Fill(imp, resEst, w);
     201
     202    fHResolutionEst.Fill(eest, resEst, w);
     203    fHResolutionMC.Fill(etru, resMC, w);
    183204
    184205    // For the fit we use a different quantity
     
    223244    }
    224245
    225     TH1D *h = (TH1D*)fHResolution.ProjectionZ("Resolution");
     246    TH1D *h = (TH1D*)fHResolutionEst.ProjectionY("Dummy", -1, -1, "s");
    226247    h->Fit("gaus", "Q0", "", -1.0, 0.25);
    227248
     
    270291    pad->cd(1);
    271292
    272     TH1D *hx=0;
    273     TH1D *hy=0;
     293    TH1 *hx=0;
     294    TH1 *hy=0;
    274295
    275296    if (pad->GetPad(1))
     
    277298        pad->GetPad(1)->cd(1);
    278299
    279         if ((hx=(TH1D*)gPad->FindObject("EnergyEst_ex")))
    280         {
    281             TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ex");
    282             hx->Reset();
    283             hx->Add(h2);
    284             delete h2;
    285         }
    286 
    287         if ((hy=(TH1D*)gPad->FindObject("EnergyEst_ey")))
    288         {
    289             TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ey");
    290             hy->Reset();
    291             hy->Add(h2);
    292             delete h2;
    293         }
     300        if (gPad->FindObject("EnergyEst_ex"))
     301            hx = fHEnergy.Project3D("ex");
     302
     303        if (gPad->FindObject("EnergyEst_ey"))
     304            hy = fHEnergy.Project3D("ey");
    294305
    295306        if (hx && hy)
    296307        {
     308            hx->SetLineColor(kBlue);
     309            hx->SetMarkerColor(kBlue);
    297310            hy->SetMaximum();
    298311            hy->SetMaximum(TMath::Max(hx->GetMaximum(), hy->GetMaximum())*1.2);
     
    305318            pad->GetPad(1)->GetPad(2)->cd(1);
    306319            if ((hx=(TH1D*)gPad->FindObject("EnergyEst_ez")))
    307             {
    308                 TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ez");
    309                 hx->Reset();
    310                 hx->Add(h2);
    311                 delete h2;
    312             }
     320                fHEnergy.Project3D("ez");
    313321
    314322            pad->GetPad(1)->GetPad(2)->cd(2);
    315             hx = (TH1D*)fHResolution.ProjectionZ("Resolution");
     323            hx = (TH1D*)fHResolutionEst.ProjectionY("Resolution", -1, -1, "e");
    316324            TPaveStats *stats = dynamic_cast<TPaveStats*>(hx->FindObject("stats"));
    317325            if (stats)
     
    323331
    324332            hx->Fit("gaus", "Q", "", -1.0, 0.25);
     333            hx->GetFunction("gaus")->SetLineColor(kBlue);
     334            hx->GetFunction("gaus")->SetLineWidth(2);
    325335            gPad=NULL;
    326336            gStyle->SetOptFit(101);
     
    331341    {
    332342        pad->GetPad(2)->cd(1);
    333         UpdatePlot(fHEnergy, "yx", kTRUE);
     343        if (gPad->FindObject("EnergyEst_yx"))
     344        {
     345            TH2D *hyx = static_cast<TH2D*>(fHEnergy.Project3D("yx"));
     346            UpdateProf(*hyx, kTRUE);
     347        }
    334348
    335349        TLine *l = (TLine*)gPad->FindObject("TLine");
     
    346360
    347361        pad->GetPad(2)->cd(2);
    348         UpdatePlot(fHResolution, "zy");
     362        UpdateProf(fHResolutionEst, kFALSE);
    349363
    350364        pad->GetPad(2)->cd(3);
    351         UpdatePlot(fHResolution, "zx");
     365        UpdateProf(fHResolutionMC, kFALSE);
    352366    }
    353367}
    354368
    355 void MHEnergyEst::UpdatePlot(TH3 &h, const char *how, Bool_t logy)
    356 {
    357     TH2D *hyx=0;
    358     if (!(hyx=(TH2D*)gPad->FindObject(MString::Format("%s_%s", h.GetName(), how))))
     369void MHEnergyEst::UpdateProf(TH2 &h, Bool_t logy)
     370{
     371    const TString pname = MString::Format("Prof%s", h.GetName());
     372
     373    if (!gPad->FindObject(pname))
    359374        return;
    360375
    361     TH2D *h2 = (TH2D*)h.Project3D(MString::Format("dum_%s", how));
    362     hyx->Reset();
    363     hyx->Add(h2);
    364     delete h2;
    365 
    366     TH1D *hx = 0;
    367     if ((hx=(TH1D*)gPad->FindObject(MString::Format("Prof%s", h.GetName()))))
    368     {
    369         hx = hyx->ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s");
    370 
    371         if (logy && hx->GetMaximum()>0)
    372             gPad->SetLogy();
    373     }
    374 }
    375 
    376 TH1 *MHEnergyEst::MakePlot(TH3 &h, const char *how)
    377 {
    378     gPad->SetBorderMode(0);
     376    TH1D *hx = h.ProfileX(pname, -1, -1, "s");
     377    hx->SetLineColor(kBlue);
     378    hx->SetMarkerColor(kBlue);
     379
     380    if (logy && hx->GetMaximum()>0)
     381        gPad->SetLogy();
     382}
     383
     384TH1 *MHEnergyEst::MakeProj(const char *how)
     385{
     386    TH1 *p = fHEnergy.Project3D(how);
     387    p->SetDirectory(NULL);
     388    p->SetBit(kCanDelete);
     389    p->SetBit(TH1::kNoStats);
     390    p->SetMarkerStyle(kFullDotMedium);
     391    p->SetLineColor(kBlue);
     392
     393    return p;
     394}
     395
     396TH1 *MHEnergyEst::MakeProf(TH2 &h)
     397{
     398    TH1 *p = h.ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s");
     399    p->SetDirectory(NULL);
     400    p->SetBit(kCanDelete);
     401    p->SetLineWidth(2);
     402    p->SetLineColor(kBlue);
     403    p->SetFillStyle(4000);
     404    p->SetStats(kFALSE);
     405
     406    return p;
     407}
     408
     409// --------------------------------------------------------------------------
     410//
     411// Draw the histogram
     412//
     413void MHEnergyEst::Draw(Option_t *opt)
     414{
     415    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
     416
     417    // Do the projection before painting the histograms into
     418    // the individual pads
     419    AppendPad("");
     420
     421    pad->SetBorderMode(0);
     422    pad->Divide(2, 1, 1e-10, 1e-10);
     423
     424    TH1 *h;
     425
     426    // ----------------------------------------
     427
     428    pad->cd(1);
     429    gPad->SetBorderMode(0);
     430
     431    gPad->Divide(1, 2, 1e-10, 1e-10);
     432
     433    TVirtualPad *pad2 = gPad;
     434
     435    // ----------------------------------------
     436
     437    pad2->cd(1);
     438    gPad->SetBorderMode(0);
     439
     440    gPad->SetGridx();
     441    gPad->SetGridy();
     442    gPad->SetLogx();
     443
     444    h = MakeProj("ey");
     445    h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)");
     446    h->SetXTitle("E [GeV]"); // E_mc
     447    h->SetYTitle("Counts");
     448    h->Draw();
     449
     450    h = MakeProj("ex");
     451    h->SetLineColor(kBlue);
     452    h->SetMarkerColor(kBlue);
     453    h->Draw("same");
     454
     455    // ----------------------------------------
     456
     457    pad2->cd(2);
     458    gPad->SetBorderMode(0);
     459
     460    TVirtualPad *pad3 = gPad;
     461    pad3->Divide(2, 1, 1e-10, 1e-10);
     462    pad3->cd(1);
     463    gPad->SetBorderMode(0);
     464    gPad->SetGridx();
     465    gPad->SetGridy();
     466
     467    h = MakeProj("ez");
     468    h->SetTitle("Zenith Angle Distribution");
     469    h->GetXaxis()->SetMoreLogLabels();
     470    h->GetXaxis()->SetNoExponent();
     471    h->Draw();
     472
     473    // ----------------------------------------
     474
     475    pad3->cd(2);
     476    gPad->SetBorderMode(0);
     477    gPad->SetGridx();
     478    gPad->SetGridy();
     479
     480    h = fHResolutionEst.ProjectionY("_py");
     481    h->SetTitle("Distribution of \\Delta E/E_{est}");
     482    h->SetDirectory(NULL);
     483    h->SetBit(kCanDelete);
     484    h->GetXaxis()->SetRangeUser(-1.3, 1.3);
     485    h->Draw();
     486    // ----------------------------------------
     487
     488    pad->cd(2);
     489    gPad->SetBorderMode(0);
     490
     491    gPad->Divide(1, 3, 1e-10, 1e-10);
     492    pad2 = gPad;
     493
     494    // ----------------------------------------
     495
     496    pad2->cd(1);
     497    gPad->SetBorderMode(0);
     498    gPad->SetLogy();
    379499    gPad->SetLogx();
    380500    gPad->SetGridx();
     
    384504    //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
    385505
    386     TH2D *h2 = (TH2D*)h.Project3D(how);
     506    TH2D *h2 = (TH2D*)fHEnergy.Project3D("yx");
    387507    h2->SetDirectory(NULL);
    388508    h2->SetBit(kCanDelete);
    389509    h2->SetFillColor(kBlue);
    390     h2->SetLineColor(kRed);
    391 
    392     TH1D *h1 = h2->ProfileX(Form("Prof%s", h.GetName()), -1, -1, "s");
    393     h1->SetDirectory(NULL);
    394     h1->SetBit(kCanDelete);
    395     h1->SetLineWidth(2);
    396     h1->SetLineColor(kBlue);
    397     h1->SetFillStyle(4000);
    398     h1->SetStats(kFALSE);
     510    h2->SetTitle("Estimated energy E_{est} vs Monte Carlo energy E_{mc}");
     511
     512    TH1 *h1 = MakeProf(*h2);
    399513
    400514    h2->Draw("");
    401     h1->Draw("E0 hist C same");
    402 //    h1->Draw("Chistsame");
    403 
    404     return h2;
    405 }
    406 
    407 
    408 // --------------------------------------------------------------------------
    409 //
    410 // Draw the histogram
    411 //
    412 void MHEnergyEst::Draw(Option_t *opt)
    413 {
    414     TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    415 
    416     // Do the projection before painting the histograms into
    417     // the individual pads
    418     AppendPad("");
    419 
    420     pad->SetBorderMode(0);
    421     pad->Divide(2, 1, 1e-10, 1e-10);
    422 
    423     TH1 *h;
    424 
    425     pad->cd(1);
    426     gPad->SetBorderMode(0);
    427 
    428     gPad->Divide(1, 2, 1e-10, 1e-10);
    429 
    430     TVirtualPad *pad2 = gPad;
    431 
    432     pad2->cd(1);
    433     gPad->SetBorderMode(0);
    434 
    435     gPad->SetGridx();
    436     gPad->SetGridy();
    437     gPad->SetLogx();
    438     h = (TH1D*)fHEnergy.Project3D("ey");
    439     h->SetBit(TH1::kNoStats);
    440     h->SetTitle("Energy disribution: Monte Carlo E_{mc} (black), Estimated E_{est} (blue)");
    441     h->SetXTitle("E [GeV]"); // E_mc
    442     h->SetYTitle("Counts");
    443     h->SetBit(kCanDelete);
    444     h->SetDirectory(NULL);
    445     h->SetMarkerStyle(kFullDotMedium);
    446     h->Draw();
    447 
    448     h = (TH1D*)fHEnergy.Project3D("ex");
    449     h->SetBit(TH1::kNoTitle|TH1::kNoStats);
    450     h->SetXTitle("E [GeV]"); // E_est
    451     h->SetYTitle("Counts");
    452     h->SetBit(kCanDelete);
    453     h->SetDirectory(NULL);
    454     h->SetMarkerStyle(kFullDotMedium);
    455     h->SetLineColor(kBlue);
    456     h->SetMarkerColor(kBlue);
    457     h->Draw("same");
    458 
    459     // FIXME: LEGEND
    460 
    461     pad2->cd(2);
    462     gPad->SetBorderMode(0);
    463 
    464     TVirtualPad *pad3 = gPad;
    465     pad3->Divide(2, 1, 1e-10, 1e-10);
    466     pad3->cd(1);
    467     gPad->SetBorderMode(0);
    468     gPad->SetGridx();
    469     gPad->SetGridy();
    470     h = fHEnergy.Project3D("ez");
    471     h->SetTitle("Zenith Angle Distribution");
    472     h->SetBit(TH1::kNoStats);
    473     h->SetDirectory(NULL);
    474     h->SetXTitle("\\Theta [\\circ]");
    475     h->SetBit(kCanDelete);
    476     h->Draw();
    477 
    478     pad3->cd(2);
    479     gPad->SetBorderMode(0);
    480     gPad->SetGridx();
    481     gPad->SetGridy();
    482     /*
    483     h = fHImpact.ProjectionX("Impact", -1, -1, "e");
    484     h->SetBit(TH1::kNoStats);
    485     h->SetTitle("Distribution of Impact");
    486     h->SetDirectory(NULL);
    487     h->SetXTitle("Impact [m]");
    488     h->SetBit(kCanDelete);
    489     h->Draw();*/
    490     // ----------------------------------------
    491     h = fHResolution.ProjectionZ("Resolution");
    492     //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
    493     h->SetYTitle("Counts");
    494     h->SetTitle("Distribution of \\Delta E/E");
    495     h->SetDirectory(NULL);
    496     h->SetBit(kCanDelete);
    497     h->GetXaxis()->SetRangeUser(-1.3, 1.3);
    498     h->Draw("");
    499     //h->Fit("gaus");
    500     // ----------------------------------------
    501 
    502     pad->cd(2);
    503     gPad->SetBorderMode(0);
    504 
    505     gPad->Divide(1, 3, 1e-10, 1e-10);
    506     pad2 = gPad;
    507 
    508     pad2->cd(1);
    509     gPad->SetLogy();
    510     h = MakePlot(fHEnergy, "xy");
    511     h->SetXTitle("E_{mc} [GeV]");
    512     h->SetYTitle("E_{est} [GeV]");
    513     h->SetTitle("Estimated energy E_{est} vs Monte Carlo energy E_{mc}");
     515    h1->Draw("E0 same");
    514516
    515517    TLine line;
     
    520522    line.SetLineStyle(kDashed);
    521523
     524    // ----------------------------------------
     525
    522526    pad2->cd(2);
    523     h = MakePlot(fHResolution, "zy");
    524     h->SetXTitle("E_{mc} [GeV]");
    525     h->SetYTitle("(E_{est}/E_{mc}-1");
    526     h->SetTitle("Energy resolution \\Delta E/E vs Monte Carlo Energy E_{mc}");
    527     h->SetMinimum(-1.3);
    528     h->SetMaximum(1.3);
    529 
    530     line.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
     527    gPad->SetBorderMode(0);
     528    gPad->SetLogx();
     529    gPad->SetGridx();
     530    gPad->SetGridy();
     531    fHResolutionEst.Draw();
     532    MakeProf(fHResolutionEst)->Draw("E0 same");
     533
     534    fHResolutionEst.GetXaxis()->SetMoreLogLabels();
     535    fHResolutionEst.GetXaxis()->SetNoExponent();
     536
     537    line.DrawLine(fHResolutionEst.GetXaxis()->GetXmin(), 0, fHResolutionEst.GetXaxis()->GetXmax(), 0);
     538
     539    // ----------------------------------------
    531540
    532541    pad2->cd(3);
    533     h = MakePlot(fHResolution, "zx");
    534     h->SetXTitle("E_{est} [GeV]");
    535     h->SetYTitle("(E_{est}/E_{mc}-1");
    536     h->SetTitle("Energy resolution \\Delta E/E vs estimated Energy E_{est}");
    537     h->SetMinimum(-1.3);
    538     h->SetMaximum(1.3);
    539 
    540     line.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
     542    gPad->SetBorderMode(0);
     543    gPad->SetLogx();
     544    gPad->SetGridx();
     545    gPad->SetGridy();
     546    fHResolutionMC.Draw();
     547    MakeProf(fHResolutionMC)->Draw("E0 same");
     548
     549    fHResolutionMC.GetXaxis()->SetMoreLogLabels();
     550    fHResolutionMC.GetXaxis()->SetNoExponent();
     551
     552    line.DrawLine(fHResolutionMC.GetXaxis()->GetXmin(), 0, fHResolutionMC.GetXaxis()->GetXmax(), 0);
    541553}
    542554
  • trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h

    r7692 r8935  
    3333
    3434    TH3D fHEnergy;
    35     TH3D fHResolution;
     35    TH2D fHResolutionEst;
     36    TH2D fHResolutionMC;
    3637    TH2D fHImpact;
    3738
     
    3940    Double_t fBias;
    4041
    41     TH1 *MakePlot(TH3 &h, const char *how);
    42     void UpdatePlot(TH3 &h, const char *how, Bool_t logy=kFALSE);
     42    TH1 *MakeProj(const char *how);
     43    TH1 *MakeProf(TH2 &h);
     44    void UpdateProf(TH2 &h, Bool_t logy);
    4345
    4446    Double_t GetVal(Int_t i) const;
     
    6264    void Print(Option_t *o="") const;
    6365
    64     //ClassDef(MHEnergyEst, 2) //
    65     ClassDef(MHEnergyEst, 1) // Histogram for the result of the energy reconstruction
     66    ClassDef(MHEnergyEst, 2) // Histogram for the result of the energy reconstruction
    6667};
    6768
Note: See TracChangeset for help on using the changeset viewer.