Ignore:
Timestamp:
04/08/04 19:52:13 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhist
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mhist/MHFalseSource.cc

    r3666 r3682  
    4646// pads) you can change the AlphaCut 'online'
    4747//
     48// Each Z-Projection (Alpha-histogram) is scaled such, that the number
     49// of entries fits the maximum number of entries in all Z-Projections.
     50// This should correct for the different acceptance in the center
     51// and at the vorder of the camera. This, however, produces more noise
     52// at the border.
     53//
    4854// Here is a slightly simplified version of the algorithm:
    4955// ------------------------------------------------------
    … …  
    139145//
    140146MHFalseSource::MHFalseSource(const char *name, const char *title)
    141     : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1),
    142       fAlphaCut(12.5), fBgMean(55), fDistMin(-1), fDistMax(-1)
     147    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5),
     148    fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
    143149{
    144150    //
    … …  
    155161    fHist.SetYTitle("y [\\circ]");
    156162    fHist.SetZTitle("\\alpha [\\circ]");
     163}
     164
     165void MHFalseSource::MakeSymmetric(TH1 *h)
     166{
     167    const Float_t min = TMath::Abs(h->GetMinimum());
     168    const Float_t max = TMath::Abs(h->GetMaximum());
     169
     170    const Float_t absmax = TMath::Max(min, max)*1.002;
     171
     172    h->SetMaximum( absmax);
     173    h->SetMinimum(-absmax);
    157174}
    158175
    … …  
    337354            //        and/or MTaskList
    338355            // Source dependant distance cut
    339             if (fDistMin>0 && hsrc.GetDist()*fMm2Deg<fDistMin)
     356            if (fMinDist>0 && hsrc.GetDist()*fMm2Deg<fMinDist)
    340357                continue;
    341 
    342             if (fDistMax>0 && hil->GetLength()>fDistMax*hsrc.GetDist())
     358            if (fMaxDist>0 && hsrc.GetDist()*fMm2Deg>fMaxDist)
     359                continue;
     360
     361            if (fMaxLD>0 && hil->GetLength()>fMaxLD*hsrc.GetDist())
     362                continue;
     363            if (fMinLD>0 && hil->GetLength()<fMinLD*hsrc.GetDist())
    343364                continue;
    344365
    … …  
    358379// the same number of bins than for on-data
    359380//
    360 void MHFalseSource::ProjectOff(TH2D *h2)
     381void MHFalseSource::ProjectOff(TH2D *h2, TH2D *all)
    361382{
    362383    TAxis &axe = *fHist.GetZaxis();
    … …  
    379400    // Move contents from projection to h2
    380401    h2->Reset();
    381     h2->Add(p);
     402    h2->Add(p, all->GetMaximum());
     403    h2->Divide(all);
    382404
    383405    // Delete p
    … …  
    393415// range (0, fAlphaCut)
    394416//
    395 void MHFalseSource::ProjectOn(TH2D *h3)
     417void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all)
    396418{
    397419    TAxis &axe = *fHist.GetZaxis();
    … …  
    410432    // Move contents from projection to h3
    411433    h3->Reset();
    412     h3->Add(p);
     434    h3->Add(p, all->GetMaximum());
     435    h3->Divide(all);
     436
     437    // Delete p
    413438    delete p;
    414439
    … …  
    454479    TH2D* h3;
    455480    TH2D* h4;
     481    TH2D* h5;
    456482
    457483    // Update projection of all-events
    458     padsave->GetPad(1)->cd(3);
     484    padsave->GetPad(2)->cd(3);
    459485    if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all")))
    460486        ProjectAll(h0);
    461487
     488    // Update projection of on-events
     489    padsave->GetPad(1)->cd(1);
     490    if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
     491        ProjectOn(h3, h0);
     492
    462493    // Update projection of off-events
    463     padsave->GetPad(2)->cd(1);
     494    padsave->GetPad(1)->cd(3);
    464495    if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
    465         ProjectOff(h2);
    466 
    467     // Update projection of on-events
     496        ProjectOff(h2, h0);
     497
    468498    padsave->GetPad(2)->cd(2);
    469     if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
    470         ProjectOn(h3);
     499    if ((h5 = (TH2D*)gPad->FindObject("Alpha_yx_diff")))
     500    {
     501        h5->Add(h3, h2, -1);
     502        MakeSymmetric(h5);
     503    }
    471504
    472505    // Update projection of significance
    473     padsave->GetPad(2)->cd(3);
     506    padsave->GetPad(1)->cd(2);
    474507    if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig")))
    475508    {
    … …  
    483516        Int_t max = h4->GetBin(nx, ny);
    484517
    485         for (int ix=0; ix<nx; ix++)
    486             for (int iy=0; iy<ny; iy++)
     518        for (int ix=1; ix<=nx; ix++)
     519            for (int iy=1; iy<=ny; iy++)
    487520            {
    488                 const Int_t n = h4->GetBin(ix+1, iy+1);
     521                const Int_t n = h4->GetBin(ix, iy);
    489522
    490523                const Double_t s = h3->GetBinContent(n);
    … …  
    503536            }
    504537
     538        MakeSymmetric(h4);
     539
    505540        // Update projection of 'the best alpha-plot'
    506         padsave->GetPad(1)->cd(1);
     541        padsave->GetPad(2)->cd(1);
    507542        if ((h1 = (TH1D*)gPad->FindObject("Alpha")) && max>0)
    508543        {
    … …  
    560595    TObject *catalog = GetCatalog();
    561596
    562     // draw the 2D histogram Sigmabar versus Theta
     597    // Initialize upper part
    563598    pad->cd(1);
    564599    gPad->SetBorderMode(0);
    565600    gPad->Divide(3, 1);
    566     delete pad->GetPad(1)->GetPad(2);
    567 
    568     pad->GetPad(1)->cd(3);
    569     gPad->SetBorderMode(0);
    570     TH1 *h0 = fHist.Project3D("yx_all");
    571     h0->SetDirectory(NULL);
    572     h0->SetXTitle(fHist.GetXaxis()->GetTitle());
    573     h0->SetYTitle(fHist.GetYaxis()->GetTitle());
    574     h0->Draw("colz");
    575     h0->SetBit(kCanDelete);
    576     catalog->Draw("mirror same");
    577 
     601
     602    // PAD #1
    578603    pad->GetPad(1)->cd(1);
    579     gPad->SetBorderMode(0);
    580 
    581     TH1 *h1 = fHist.ProjectionZ("Alpha");
    582     h1->SetDirectory(NULL);
    583     h1->SetTitle("Distribution of \\alpha");
    584     h1->SetXTitle(fHist.GetZaxis()->GetTitle());
    585     h1->SetYTitle("Counts");
    586     h1->Draw(opt);
    587     h1->SetBit(kCanDelete);
    588 
    589     pad->cd(2);
    590     gPad->SetBorderMode(0);
    591     gPad->Divide(3, 1);
    592 
    593     pad = gPad;
    594 
    595     pad->cd(1);
    596     gPad->SetBorderMode(0);
    597     fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
    598     TH1 *h2 = fHist.Project3D("yx_off");
    599     h2->SetDirectory(NULL);
    600     h2->SetXTitle(fHist.GetXaxis()->GetTitle());
    601     h2->SetYTitle(fHist.GetYaxis()->GetTitle());
    602     h2->Draw("colz");
    603     h2->SetBit(kCanDelete);
    604     catalog->Draw("mirror same");
    605 
    606     pad->cd(2);
    607604    gPad->SetBorderMode(0);
    608605    fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
    … …  
    616613    catalog->Draw("mirror same");
    617614
    618     pad->cd(3);
     615    // PAD #2
     616    pad->GetPad(1)->cd(2);
    619617    gPad->SetBorderMode(0);
    620618    fHist.GetZaxis()->SetRange(0,0);
    … …  
    627625    h4->Draw("colz");
    628626    h4->SetBit(kCanDelete);
     627    catalog->Draw("mirror same");
     628
     629    // PAD #3
     630    pad->GetPad(1)->cd(3);
     631    gPad->SetBorderMode(0);
     632    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
     633    TH1 *h2 = fHist.Project3D("yx_off");
     634    h2->SetDirectory(NULL);
     635    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
     636    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
     637    h2->Draw("colz");
     638    h2->SetBit(kCanDelete);
     639    catalog->Draw("mirror same");
     640
     641    // Initialize lower part
     642    pad->cd(2);
     643    gPad->SetBorderMode(0);
     644    gPad->Divide(3, 1);
     645
     646    // PAD #4
     647    pad->GetPad(2)->cd(1);
     648    gPad->SetBorderMode(0);
     649    TH1 *h1 = fHist.ProjectionZ("Alpha");
     650    h1->SetDirectory(NULL);
     651    h1->SetTitle("Distribution of \\alpha");
     652    h1->SetXTitle(fHist.GetZaxis()->GetTitle());
     653    h1->SetYTitle("Counts");
     654    h1->Draw(opt);
     655    h1->SetBit(kCanDelete);
     656
     657    // PAD #5
     658    pad->GetPad(2)->cd(2);
     659    gPad->SetBorderMode(0);
     660    TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff");
     661    h5->Add(h2, -1);
     662    h5->SetTitle("Difference of on- and off-distribution");
     663    h5->SetDirectory(NULL);
     664    h5->Draw("colz");
     665    h5->SetBit(kCanDelete);
     666    catalog->Draw("mirror same");
     667
     668    // PAD #6
     669    pad->GetPad(2)->cd(3);
     670    gPad->SetBorderMode(0);
     671    TH1 *h0 = fHist.Project3D("yx_all");
     672    h0->SetDirectory(NULL);
     673    h0->SetXTitle(fHist.GetXaxis()->GetTitle());
     674    h0->SetYTitle(fHist.GetYaxis()->GetTitle());
     675    h0->Draw("colz");
     676    h0->SetBit(kCanDelete);
    629677    catalog->Draw("mirror same");
    630678}
    … …  
    729777                       Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
    730778                            sigmax, bgmin, bgmax));
    731     hists->SetNameTitle("Excess",     Form("Number of excess events for \\alpha<%.0f\\circ", sigint));
    732     histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint));
     779    hists->SetName("Excess");
     780    histb->SetName("Background");
    733781    hist->SetXTitle(fHist.GetXaxis()->GetTitle());
    734782    hists->SetXTitle(fHist.GetXaxis()->GetTitle());
    … …  
    776824    *fLog << "Polynom order:  " << (int)polynom << endl;
    777825    *fLog << "Fitting False Source Plot..." << flush;
     826
     827    TH1 *h0 = fHist.Project3D("yx_entries");
     828    Float_t entries = h0->GetMaximum();
     829    delete h0;
    778830
    779831    TH1 *h=0;
    … …  
    790842            if (alpha0==0)
    791843                continue;
     844
     845            h->Scale(entries/h->GetEntries());
    792846
    793847            if (alpha0>maxalpha0)
    … …  
    887941    h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    888942
     943    hists->SetTitle(Form("Excess events for \\alpha<%.0f\\circ (N_{max}=%d)", sigint, (int)hists->GetMaximum()));
     944    histb->SetTitle(Form("Background events for \\alpha<%.0f\\circ", sigint));
     945
    889946    //hists->SetMinimum(GetMinimumGT(*hists));
    890947    histb->SetMinimum(GetMinimumGT(*histb));
     948
     949    MakeSymmetric(hists);
     950    MakeSymmetric(hist);
    891951
    892952    clk.Stop();
    … …  
    9501010
    9511011        TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
     1012        result->Scale(entries/h->GetEntries());
     1013
    9521014        result->SetDirectory(NULL);
    9531015        result->SetNameTitle("Result \\alpha", title);

  • trunk/MagicSoft/Mars/mhist/MHFalseSource.h

    r3666 r3682  
    2929    Float_t fBgMean;             // Background mean
    3030
    31     Float_t fDistMin;            // Min dist
    32     Float_t fDistMax;            // Maximum distance in percent of dist
     31    Float_t fMinDist;            // Min dist
     32    Float_t fMaxDist;            // Max dist
     33
     34    Float_t fMinLD;              // Minimum distance in percent of dist
     35    Float_t fMaxLD;              // Maximum distance in percent of dist
    3336
    3437    TH3D    fHist;               // Alpha vs. x and y
    … …  
    4043    void Modified();
    4144
    42     void ProjectOff(TH2D *h);
    43     void ProjectOn(TH2D *h);
    4445    void ProjectAll(TH2D *h);
     46    void ProjectOff(TH2D *h, TH2D *all);
     47    void ProjectOn(TH2D *h, TH2D *all);
    4548
    4649    TObject *GetCatalog();
     50
     51    void MakeSymmetric(TH1 *h);
    4752
    4853public:
    … …  
    5762    void FitSignificanceStd() { FitSignificance(); } //*MENU*
    5863
    59     void SetDistMin(Float_t dist)  { fDistMin = dist;  } // Absolute minimum distance
    60     void SetDistMax(Float_t ratio) { fDistMax = ratio; } // Maximum ratio between length/dist
     64    void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance
     65    void SetMaxDist(Float_t dist) { fMaxDist = dist; } // Absolute maximum distance
     66    void SetMinLD(Float_t ratio)  { fMinLD = ratio; }  // Minimum ratio between length/dist
     67    void SetMaxLD(Float_t ratio)  { fMaxLD = ratio; }  // Maximum ratio between length/dist
    6168
    6269    void SetAlphaCut(Float_t alpha); //*MENU*
Note: See TracChangeset for help on using the changeset viewer.