Ignore:
Timestamp:
03/22/04 09:25:49 (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

    r3557 r3568  
    9595//  - a more clever (and faster) algorithm to fill the histogram, eg by
    9696//    calculating alpha once and fill the two coils around the mean
     97//  - more drawing options...
     98//  - Move Significance() to a more 'general' place and implement
     99//    also different algorithms like (Li/Ma)
     100//  - implement fit for best alpha distribution -- online (Paint)
    97101//
    98102//////////////////////////////////////////////////////////////////////////////
    … …  
    101105#include <TF1.h>
    102106#include <TH2.h>
     107#include <TGraph.h>
    103108#include <TStyle.h>
    104109#include <TCanvas.h>
    105110#include <TPaveText.h>
     111#include <TStopwatch.h>
    106112
    107113#include "MGeomCam.h"
    … …  
    236242// Calculate Significance as
    237243// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
     244//
     245// s: total number of events in signal region
     246// b: number of background events in signal region
    238247//
    239248Double_t MHFalseSource::Significance(Double_t s, Double_t b)
    … …  
    245254}
    246255
     256// --------------------------------------------------------------------------
     257//
     258//  calculates the significance according to Li & Ma
     259//  ApJ 272 (1983) 317
     260//
     261//    s: total number of events in signal region
     262//
     263//    Double_t fGamma;   // Nbg = Non - gamma * Noff
     264//  - the effective number of background events (fNoff), and fGamma :
     265//
     266//    fGamma = b / fNoff;
     267//    fGamma = fdNbg / sqrt(fNoff);
     268//    fGamma = fdNbg*fdNbg / fNbg;
     269//    fNoff  =  b*b  / (fdNbg*fdNbg);
     270//     Double_t fNbg;     // number of background events in signal region
     271//    b = fNOff *fGamma;
     272/*
     273Double_t MHFindSignificance::SignificanceLiMa(Double_t non, Double_t noff, Double_t gamma)
     274{
     275    if (gamma <= 0.0  ||  non <= 0.0  ||  noff <= 0.0)
     276    {
     277        *siglima = 0.0;
     278        return kFALSE;
     279    }
     280
     281    Double_t help1 = TMath::Log( (gamma+1)*s  / (gamma*(s+noff)) );
     282    Double_t help2 = TMath::Log( (gamma+1)*noff / (       s+noff ) );
     283
     284    Double_t siglima = TMath::Sqrt((s*help1+noff*help2)*2);
     285
     286    return non<gamma*noff ? -siglima : siglima;
     287}
     288*/
    247289// --------------------------------------------------------------------------
    248290//
    … …  
    276318
    277319        MBinning b;
    278         b.SetEdges(100, -r, r);
     320        b.SetEdges(20, -r, r);
    279321        SetBinning(&fHist, &b, &b, &binsa);
    280322    }
    … …  
    311353    Double_t rho = 0;
    312354    if (fTime && fObservatory && fPointPos)
    313     {
    314         const Double_t ra  = fPointPos->GetRa();
    315         const Double_t dec = fPointPos->GetDec();
    316 
    317         rho = MAstroSky2Local(*fTime, *fObservatory).RotationAngle(ra, dec);
    318     }
     355        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
     356    //if (fPointPos)
     357    //    rho = fPointPos->RotationAngle(*fObservatory);
    319358
    320359    MSrcPosCam src;
    … …  
    554593    h4->Draw("colz");
    555594    h4->SetBit(kCanDelete);
    556 
    557595}
    558596
    … …  
    587625    gPad->cd();
    588626}
    589 
     627/*
    590628Double_t fcn(Double_t *arg, Double_t *p)
    591629{
    … …  
    595633
    596634    const Double_t f1 = p[0]*TMath::Exp(-0.5*dx*dx);
    597     const Double_t f2 = p[3]*x*x + p[4];
     635    const Double_t f2 = p[3] + p[5]*x*x;
    598636
    599637    return f1 + f2;
    … …  
    602640Double_t FcnI1(Double_t x, Double_t *p)
    603641{
    604     return (p[3]*x*x/3+p[4])*x;
     642    return (p[5]*x*x/3+p[3])*x;
    605643}
    606644Double_t FcnI2(Double_t x, Double_t *p)
    … …  
    615653    return f2;
    616654}
    617 
    618 
    619 void MHFalseSource::FitSignificance(Float_t sigmax, Float_t bgmin, Float_t bgmax)
    620 {
    621     TH1D h0a("A",     "Parameter A",       50,  0, 4000);
    622     TH1D h0b("a",     "Parameter a",       50,  0, 4000);
    623     TH1D h1("mu",     "Parameter \\mu",    50, -1, 1);
    624     TH1D h2("sigma",  "Parameter \\sigma", 50,  0, 90);
    625     TH1D h3("b",      "Parameter b",       50, -0.1, 0.1);
    626     TH1D h4a("chisq1", "\\chi^{2} (red, green) / significance (black)",       50,  0, 35);
    627     TH1D h5a("prob1",  "Fit probability",  50,  0, 1.1);
    628     TH1D h4b("chisq2", "\\chi^{2} (red, green) / significance (black)",       50,  0, 35);
    629     TH1D h5b("prob",  "Fit probability",  50,  0, 1.1);
    630     TH1D h6("Signif",  "Significance",     50,  -20, 20);
    631     h0a.SetDirectory(NULL);
    632     h0b.SetDirectory(NULL);
    633     h1.SetDirectory(NULL);
    634     h2.SetDirectory(NULL);
    635     h3.SetDirectory(NULL);
    636     h4a.SetDirectory(NULL);
    637     h5b.SetDirectory(NULL);
    638     h4a.SetDirectory(NULL);
    639     h5b.SetDirectory(NULL);
    640     h6.SetDirectory(NULL);
     655*/
     656/*
     657class MHSignificance : public MH
     658{
     659private:
     660    TH1D fHist;
     661
     662    MParameterD *fParam;
     663
     664public:
     665    MHSignificance() : fParam(0)
     666    {
     667        fHist.SetName("Alpha");
     668        fHist.SetTitle("Distribution of \\alpha");
     669        fHist.SetXTitle("\\alpha [\\circ]");
     670        fHist.SetYTitle("Counts");
     671    }
     672    Int_t SetupFill(MParList *p)
     673    {
     674        fHist.Reset();
     675
     676        fParam = (MParameterD*)p->FindCreateObj("Significance", "MParameterD");
     677        if (fParam)
     678            return kFALSE;
     679
     680        return kTRUE;
     681    }
     682    Int_t Process(MParContainer *p, Double_t w=1)
     683    {
     684        MHillasSrc *hil = dynamic_cast<MHillasSrc*>(p);
     685        if (!hil)
     686        {
     687            *fLog << err << dbginf << "Got no MHillasSrc as argument of Fill()..." << endl;
     688            return kFALSE;
     689        }
     690
     691        fHist->Fill(hil->GetAlpha(), w);
     692
     693        return kTRUE;
     694    }
     695    Int_t Finalize()
     696    {
     697        if (fHist.GetEntries()==0)
     698        {
     699            *fLog << err << "Histogram empty." << endl;
     700            return kFALSE;
     701        }
     702
     703        Float_t sigmax=15;
     704        Float_t bgmin =45;
     705        Float_t bgmax =80;
     706
     707        fHist.SetNameTitle("Significance",
     708                           Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
     709                                sigmax, bgmin, bgmax));
     710
     711        // Implementing the function yourself is not faster at all!
     712        TF1 func("gaus(0) + pol2(3)", fcn, 0, 90, 6);
     713        TArrayD maxpar(func.GetNpar());
     714
     715        func.FixParameter(1, 0);
     716        func.FixParameter(4, 0);
     717        func.SetParLimits(3, -1, 1);
     718
     719        const Double_t alpha0 = fHist.GetBinContent(1);
     720
     721        // First fit a polynom in the off region
     722        func.FixParameter(0, 0);
     723        func.FixParameter(2, 1);
     724        func.ReleaseParameter(3);
     725        func.ReleaseParameter(5);
     726
     727        h->Fit(&func, "N0Q", "", bgmin, bgmax);
     728
     729        // Now fit a gaus in the on region on top of the polynom
     730        func.SetParameter(0, alpha0-func.GetParameter(4));
     731        func.SetParameter(2, sigmax*0.75);
     732
     733        func.ReleaseParameter(0);
     734        func.ReleaseParameter(2);
     735        func.FixParameter(3, func.GetParameter(3));
     736        func.FixParameter(5, func.GetParameter(5));
     737
     738        func.SetParLimits(2, 0, 80);
     739        h->Fit(&func, "N0Q", "", 0, sigmax);
     740
     741        TArrayD p(func.GetNpar(), func.GetParameters());
     742
     743        Double_t sig=0;
     744
     745        const Int_t n = hist->GetBin(ix+1, iy+1);
     746        if (!(func.GetParameter(0)>alpha0*2 ||
     747              func.GetParameter(2)<2.5      ||
     748              func.GetParameter(2)>70))
     749        {
     750            // Implementing the integral as analytical function
     751            // gives the same result in the order of 10e-5
     752            // and it is not faster at all...
     753            const Double_t s = func.Integral(0, 15);
     754
     755            func.SetParameter(0, 0);
     756            func.SetParameter(2, 1);
     757
     758            const Double_t b = func.Integral(0, 15);
     759
     760            sig = Significance(s, b);
     761        }
     762
     763        fParam->SetValue(sig);
     764        fParam->SetReadyToSave();
     765        return kTRUE;
     766    }
     767    ClassDef(MHSignificance, 0)
     768};
     769*/
     770
     771
     772// --------------------------------------------------------------------------
     773//
     774// This is a preliminary implementation of a alpha-fit procedure for
     775// all possible source positions. It will be moved into its own
     776// more powerfull class soon.
     777//
     778// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
     779//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
     780// or
     781//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
     782//
     783// Parameter [1] is fixed to 0 while the alpha peak should be
     784// symmetric around alpha=0.
     785//
     786// Parameter [4] is fixed to 0 because the first derivative at
     787// alpha=0 should be 0, too.
     788//
     789// In a first step the background is fitted between bgmin and bgmax,
     790// while the parameters [0]=0 and [2]=1 are fixed.
     791//
     792// In a second step the signal region (alpha<sigmax) is fittet using
     793// the whole function with parameters [1], [3], [4] and [5] fixed.
     794//
     795// The number of excess and background events are calculated as
     796//   s = int(0, sigint, gaus(0)+pol2(3))
     797//   b = int(0, sigint,         pol2(3))
     798//
     799// The Significance is calculated using the Significance() member
     800// function.
     801//
     802void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
     803{
     804    TH1D h0a("A",          "", 50,   0, 4000);
     805    TH1D h4a("chisq1",     "", 50,   0,   35);
     806    //TH1D h5a("prob1",      "", 50,   0,  1.1);
     807    TH1D h6("signifcance", "", 50, -20,   20);
     808
     809    TH1D h1("mu",    "Parameter \\mu",    50,   -1,    1);
     810    TH1D h2("sigma", "Parameter \\sigma", 50,    0,   90);
     811    TH1D h3("b",     "Parameter b",       50, -0.1,  0.1);
     812
     813    TH1D h0b("a",         "Parameter a (red), A (blue)", 50, 0, 4000);
     814    TH1D h4b("\\chi^{2}", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35);
     815    //TH1D h5b("prob",      "Fit probability: Bg(red), F(blue)", 50, 0, 1.1);
    641816
    642817    h0a.SetLineColor(kBlue);
    643818    h4a.SetLineColor(kBlue);
    644     h5a.SetLineColor(kBlue);
     819    //h5a.SetLineColor(kBlue);
    645820    h0b.SetLineColor(kRed);
    646821    h4b.SetLineColor(kRed);
    647     h5b.SetLineColor(kRed);
    648 
    649     TH1 *hist = fHist.Project3D("xy_fit");
     822    //h5b.SetLineColor(kRed);
     823
     824    TH1 *hist  = fHist.Project3D("xy_fit");
     825    hist->SetDirectory(0);
     826    TH1 *hists = fHist.Project3D("xy_fit");
     827    hists->SetDirectory(0);
     828    TH1 *histb = fHist.Project3D("xy_fit");
     829    histb->SetDirectory(0);
    650830    hist->Reset();
     831    hists->Reset();
     832    histb->Reset();
    651833    hist->SetNameTitle("Significance",
    652834                       Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
    653835                            sigmax, bgmin, bgmax));
     836    hists->SetNameTitle("Excess",     Form("Number of excess events for \\alpha<%.0f\\circ", sigint));
     837    histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint));
    654838    hist->SetXTitle(fHist.GetXaxis()->GetTitle());
     839    hists->SetXTitle(fHist.GetXaxis()->GetTitle());
     840    histb->SetXTitle(fHist.GetXaxis()->GetTitle());
    655841    hist->SetYTitle(fHist.GetXaxis()->GetTitle());
     842    hists->SetYTitle(fHist.GetXaxis()->GetTitle());
     843    histb->SetYTitle(fHist.GetXaxis()->GetTitle());
    656844
    657845    //                      xmin, xmax, npar
    658     TF1 func("MyFunc", fcn, 0,    90,   5);
    659     TArrayD par(5);
    660 
    661     func.SetParName(0, "A");
    662     func.SetParName(1, "mu");
    663     func.SetParName(2, "sigma");
    664     func.SetParName(3, "a");
    665     func.SetParName(4, "b");
    666 
     846    //TF1 func("MyFunc", fcn, 0, 90, 6);
     847    // Implementing the function yourself is only about 5% faster
     848    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
     849    TArrayD maxpar(func.GetNpar());
     850
     851    /*
     852     func.SetParName(0, "A");
     853     func.SetParName(1, "mu");
     854     func.SetParName(2, "sigma");
     855     func.SetParName(3, "a");
     856     func.SetParName(4, "b");
     857     func.SetParName(5, "c");
     858    */
     859
     860    func.FixParameter(1, 0);
     861    func.FixParameter(4, 0);
     862    func.SetParLimits(2, 0, 90);
    667863    func.SetParLimits(3, -1, 1);
    668864
    … …  
    676872    Int_t maxy=0;
    677873
     874    TStopwatch clk;
     875    clk.Start();
     876
     877    *fLog << inf;
     878    *fLog << "Signal fit:     alpha < " << sigmax << endl;
     879    *fLog << "Integration:    alpha < " << sigint << endl;
     880    *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
     881    *fLog << "Polynom order:  " << (int)polynom << endl;
     882    *fLog << "Fitting False Source Plot..." << flush;
     883
    678884    TH1 *h=0;
    679     for (int ix=0; ix<nx; ix++)
    680         for (int iy=0; iy<ny; iy++)
     885    for (int ix=1; ix<=nx; ix++)
     886        for (int iy=1; iy<=ny; iy++)
    681887        {
    682             h = fHist.ProjectionZ("AlphaFit", ix+1, ix+1, iy+1, iy+1);
     888            h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
    683889
    684890            const Double_t alpha0 = h->GetBinContent(1);
     891            const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);
    685892
    686893            // Check for the regios which is not filled...
    … …  
    693900            // First fit a polynom in the off region
    694901            func.FixParameter(0, 0);
    695             func.FixParameter(1, 0);
    696902            func.FixParameter(2, 1);
    697903            func.ReleaseParameter(3);
    698             func.ReleaseParameter(4);
     904
     905            for (int i=5; i<func.GetNpar(); i++)
     906                func.ReleaseParameter(i);
    699907
    700908            h->Fit(&func, "N0Q", "", bgmin, bgmax);
    … …  
    702910
    703911            h4a.Fill(func.GetChisquare());
    704             h5a.Fill(func.GetProb());
    705 
    706             // Now fit a gaus in the on region on top of the polynom
    707             func.SetParameter(0, alpha0-func.GetParameter(4));
    708             func.SetParameter(2, sigmax*0.75);
     912            //h5a.Fill(func.GetProb());
    709913
    710914            //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
    … …  
    715919            func.ReleaseParameter(2);
    716920            func.FixParameter(3, func.GetParameter(3));
    717             func.FixParameter(4, func.GetParameter(4));
    718 
    719             func.SetParLimits(2, 0, 80);
     921            for (int i=5; i<func.GetNpar(); i++)
     922                func.FixParameter(i, func.GetParameter(i));
     923
     924            // Do not allow signals smaller than the background
     925            const Double_t A  = alpha0-func.GetParameter(3);
     926            const Double_t dA = TMath::Abs(A);
     927            func.SetParLimits(0, -dA*4, dA*4);
     928            func.SetParLimits(2, 0, 90);
     929
     930            // Now fit a gaus in the on region on top of the polynom
     931            func.SetParameter(0, A);
     932            func.SetParameter(2, sigmax*0.75);
     933
    720934            h->Fit(&func, "N0Q", "", 0, sigmax);
    721935            //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
    722936
     937            TArrayD p(func.GetNpar(), func.GetParameters());
     938
    723939            // Fill results into some histograms
    724             h0a.Fill(func.GetParameter(0));
    725             h0b.Fill(func.GetParameter(4));
    726             h1.Fill(func.GetParameter(1));
    727             h2.Fill(func.GetParameter(2));
    728             h3.Fill(func.GetParameter(3));
     940            h0a.Fill(p[0]);
     941            h0b.Fill(p[3]);
     942            h1.Fill(p[1]);
     943            h2.Fill(p[2]);
     944            if (polynom>1)
     945                h3.Fill(p[5]);
    729946            h4b.Fill(func.GetChisquare());
    730             h5b.Fill(func.GetProb());
    731 
    732             Double_t sig=0;
    733 
    734             const Int_t n = hist->GetBin(ix+1, iy+1);
    735             if (!(func.GetParameter(0)>alpha0*2 ||
    736                   func.GetParameter(2)<2.5      ||
    737                   func.GetParameter(2)>70
    738                 /*func.GetProb()<0.005 ||*/))
     947            //h5b.Fill(func.GetProb());
     948
     949            // Implementing the integral as analytical function
     950            // gives the same result in the order of 10e-5
     951            // and it is not faster at all...
     952            //const Bool_t ok = /*p[0]>=0 && /*p[0]<alpha0*2 &&*/ p[2]>1.75;// && p[2]<88.5;
     953            /*
     954            if (p[0]<-fabs(alpha0-p[3])*7 && p[2]<alphaw/3)
    739955            {
    740                 const Double_t b = FcnI1(15, func.GetParameters());
    741                 const Double_t s = FcnI2(15, func.GetParameters());
    742 
    743                 sig = Significance(s+b, b);
     956                func.SetParameter(0, alpha0-p[3]);
     957                cout << p[2] << " " << p[0] << " " << alpha0-p[3] << endl;
    744958            }
    745 
    746             hist->SetBinContent(n, sig==0?1e-15:sig);
     959            */
     960
     961            const Double_t s = func.Integral(0, sigint);
     962
     963            func.SetParameter(0, 0);
     964            func.SetParameter(2, 1);
     965
     966            const Double_t b   = func.Integral(0, sigint);
     967            const Double_t sig = Significance(s, b);
     968
     969            const Int_t n = hist->GetBin(ix, iy);
     970            hists->SetBinContent(n, s-b);
     971            histb->SetBinContent(n, b);
     972
     973            hist->SetBinContent(n, sig);
    747974            if (sig!=0)
    748975                h6.Fill(sig);
    … …  
    751978            {
    752979                maxs = sig;
    753                 maxx = ix+1;
    754                 maxy = iy+1;
    755                 for(int i=0; i<par.GetSize(); i++)
    756                     par[i] = func.GetParameter(i);
     980                maxx = ix;
     981                maxy = iy;
     982                maxpar = p;
    757983            }
    758984        }
    759985
     986    *fLog << inf << "done." << endl;
     987
    760988    h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    761989    h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    762990
     991    //hists->SetMinimum(GetMinimumGT(*hists));
     992    histb->SetMinimum(GetMinimumGT(*histb));
     993
     994    clk.Stop();
     995    clk.Print();
     996
    763997    TCanvas *c=new TCanvas;
    764998
    … …  
    7661000    c->cd(1);
    7671001    gPad->SetBorderMode(0);
    768     gPad->Divide(1,2, 0, 0);
     1002    hists->Draw("colz");
     1003    hists->SetBit(kCanDelete);
     1004    c->cd(2);
     1005    gPad->SetBorderMode(0);
     1006    hist->Draw("colz");
     1007    hist->SetBit(kCanDelete);
     1008    c->cd(3);
     1009    gPad->SetBorderMode(0);
     1010    histb->Draw("colz");
     1011    histb->SetBit(kCanDelete);
     1012    c->cd(4);
     1013    gPad->Divide(1,3, 0, 0);
    7691014    TVirtualPad *p=gPad;
    7701015    p->SetBorderMode(0);
    … …  
    7761021    gPad->SetBorderMode(0);
    7771022    h3.DrawCopy();
    778     c->cd(2);
    779     gPad->SetBorderMode(0);
    780     hist->Draw("colz");
    781     hist->SetDirectory(NULL);
    782     hist->SetBit(kCanDelete);
    783     c->cd(3);
     1023    p->cd(3);
     1024    gPad->SetBorderMode(0);
     1025    h2.DrawCopy();
     1026    c->cd(6);
     1027    gPad->Divide(1,2, 0, 0);
     1028    TVirtualPad *q=gPad;
     1029    q->SetBorderMode(0);
     1030    q->cd(1);
     1031    gPad->SetBorderMode(0);
    7841032    gPad->SetBorderMode(0);
    7851033    h4b.DrawCopy();
    7861034    h4a.DrawCopy("same");
    7871035    h6.DrawCopy("same");
    788     c->cd(4);
    789     gPad->SetBorderMode(0);
    790     h2.DrawCopy();
    791     c->cd(6);
    792     gPad->SetBorderMode(0);
    793     h5b.DrawCopy();
    794     h5a.DrawCopy("same");
    795 
     1036    q->cd(2);
     1037    gPad->SetBorderMode(0);
     1038    //h5b.DrawCopy();
     1039    //h5a.DrawCopy("same");
    7961040    c->cd(5);
    7971041    gPad->SetBorderMode(0);
    … …  
    8011045                                 fHist.GetXaxis()->GetBinCenter(maxx),
    8021046                                 fHist.GetYaxis()->GetBinCenter(maxy), maxs);
     1047
    8031048        TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
     1049        result->SetDirectory(NULL);
    8041050        result->SetNameTitle("Result \\alpha", title);
    8051051        result->SetBit(kCanDelete);
    … …  
    8081054        result->Draw();
    8091055
    810         TF1 f1("MyFunc1", fcn, 0, 90, 5);
    811         TF1 f2("MyFunc2", fcn, 0, 90, 5);
    812         f1.SetParameters(par.GetArray());
    813         f2.SetParameters(par.GetArray());
     1056        TF1 f1("", func.GetExpFormula(), 0, 90);
     1057        TF1 f2("", func.GetExpFormula(), 0, 90);
     1058        f1.SetParameters(maxpar.GetArray());
     1059        f2.SetParameters(maxpar.GetArray());
    8141060        f2.FixParameter(0, 0);
    8151061        f2.FixParameter(1, 0);
    … …  
    8241070        leg->SetBorderSize(1);
    8251071        leg->SetTextSize(0.04);
    826         leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
     1072        leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
     1073        //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
    8271074        leg->AddLine(0, 0.65, 0, 0.65);
    828         leg->AddText(0.06, 0.54, Form("A=%.2f", par[0]))->SetTextAlign(11);
    829         leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", par[2]))->SetTextAlign(11);
    830         leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fixed)", par[1]))->SetTextAlign(11);
    831         leg->AddText(0.60, 0.54, Form("a=%.2f", par[3]))->SetTextAlign(11);
    832         leg->AddText(0.60, 0.34, Form("b=%.2f", par[4]))->SetTextAlign(11);
     1075        leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
     1076        leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
     1077        leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
     1078        leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
     1079        leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
     1080        if (polynom>1)
     1081            leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
    8331082        leg->SetBit(kCanDelete);
    8341083        leg->Draw();
    835     }
    836 
    837 }
     1084
     1085        q->cd(2);
     1086
     1087        TGraph *g = new TGraph;
     1088        g->SetPoint(0, 0, 0);
     1089
     1090        Int_t max=0;
     1091        Float_t maxsig=0;
     1092        for (int i=1; i<89; i++)
     1093        {
     1094            const Double_t s = f1.Integral(0, (float)i);
     1095            const Double_t b = f2.Integral(0, (float)i);
     1096
     1097            const Double_t sig = Significance(s, b);
     1098
     1099            g->SetPoint(g->GetN(), i, sig);
     1100
     1101            if (sig>maxsig)
     1102            {
     1103                max = i;
     1104                maxsig = sig;
     1105            }
     1106        }
     1107
     1108        g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
     1109        g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
     1110        g->GetHistogram()->SetYTitle("Significance");
     1111        g->SetBit(kCanDelete);
     1112        g->Draw("AP");
     1113
     1114        leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
     1115        leg->SetBorderSize(1);
     1116        leg->SetTextSize(0.1);
     1117        leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max));
     1118        leg->SetBit(kCanDelete);
     1119        leg->Draw();
     1120    }
     1121}

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

    r3557 r3568  
    5151    TH1 *GetHistByName(const TString name) { return &fHist; }
    5252
    53     void FitSignificance(Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80); //*MENU*
     53    void FitSignificance(Float_t sigint=15, Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80, Byte_t polynom=2); //*MENU*
    5454    void FitSignificanceStd() { FitSignificance(); } //*MENU*
    5555
Note: See TracChangeset for help on using the changeset viewer.