Ignore:
Timestamp:
03/19/04 16:07:28 (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

    r3548 r3550  
    233233// --------------------------------------------------------------------------
    234234//
     235// Calculate Significance as
     236// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
     237//
     238Double_t MHFalseSource::Significance(Double_t s, Double_t b)
     239{
     240    const Double_t k = b==0 ? 0 : s/b;
     241    const Double_t f = s+k*k*b;
     242
     243    return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
     244}
     245
     246// --------------------------------------------------------------------------
     247//
    235248// Set binnings (takes BinningFalseSource) and prepare filling of the
    236249// histogram.
    … …  
    262275
    263276        MBinning b;
    264         b.SetEdges(40, -r, r);
     277        b.SetEdges(100, -r, r);
    265278        SetBinning(&fHist, &b, &b, &binsa);
    266279    }
    … …  
    455468                const Double_t b = h2->GetBinContent(n);
    456469
    457                 const Double_t k = b==0 ? 0 : s/b;
    458                 const Double_t f = s+k*k*b;
    459 
    460                 const Double_t sig = f==0 ? 0 : (s-b)/TMath::Sqrt(f);
    461 
    462                 //if (b>0 && s>0)
    463                 //    *fLog << dbg << b << " " << s << endl;
    464 
    465                 // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
     470                const Double_t sig = Significance(s, b);
     471
    466472                h4->SetBinContent(n, sig);
    467473
    … …  
    592598    return f1 + f2;
    593599}
    594 /*
    595 Double_t FcnIntegral(Double_t *arg, Double_t *p)
    596 {
    597     static const Double_t sqrt2   = TMath::Sqrt(2);
     600
     601Double_t FcnI1(Double_t x, Double_t *p)
     602{
     603    return (p[3]*x*x/3+p[4])*x;
     604}
     605Double_t FcnI2(Double_t x, Double_t *p)
     606{
     607    static const Double_t sqrt2   = TMath::Sqrt(2.);
    598608    static const Double_t sqrt2pi = TMath::Sqrt(TMath::TwoPi());
    599609
    600     const Double_t x = arg[0];
    601 
    602610    const Double_t dx = (x-p[1])/p[2];
    603611
    604     const Double_t f1 = p[0]*p[2]*sqrt2pi*TMath::Erf(dx/sqrt2)/2;
    605     const Double_t f2 = p[3]*x*2;
    606 
    607     return f1+f2;
    608 }
    609 */
    610 /*
     612    const Double_t f2 = p[0]*p[2]*sqrt2pi*TMath::Erf(dx/sqrt2)/2;
     613
     614    return f2;
     615}
     616
     617
    611618void MHFalseSource::FitSignificance()
    612619{
    613     TH1D h0("A",     "Parameter A",     50, 0, 10000);
    614     TH1D h1("mu",    "Parameter mu",    50, -1, 1);
    615     TH1D h2("sigma", "Parameter sigma", 50, 0, 20);
    616     TH1D h3("b",     "Parameter b",     50, 0.001, 0.01);
    617     h0.SetDirectory(NULL);
     620    TH1D h0a("A",     "Parameter A",      50,  0, 4000);
     621    TH1D h0b("a",     "Parameter a",      50,  0, 4000);
     622    TH1D h1("mu",    "Parameter mu",      50, -1, 1);
     623    TH1D h2("sigma", "Parameter sigma",   50,  0, 90);
     624    TH1D h3("b",     "Parameter b",       50, -0.1, 0.1);
     625    TH1D h4a("chisq1", "\\chi^{2}",       50,  0, 35);
     626    TH1D h5a("prob1",  "Fit probability", 50,  0, 1.1);
     627    TH1D h4b("chisq2", "\\chi^{2}",       50,  0, 35);
     628    TH1D h5b("prob2",  "Fit probability", 50,  0, 1.1);
     629    TH1D h6("Signif",  "Significance", 50,  -20, 20);
     630    h0a.SetDirectory(NULL);
     631    h0b.SetDirectory(NULL);
    618632    h1.SetDirectory(NULL);
    619633    h2.SetDirectory(NULL);
    620634    h3.SetDirectory(NULL);
     635    h4a.SetDirectory(NULL);
     636    h5b.SetDirectory(NULL);
     637    h4a.SetDirectory(NULL);
     638    h5b.SetDirectory(NULL);
     639    h6.SetDirectory(NULL);
     640
     641    h0a.SetLineColor(kBlue);
     642    h4a.SetLineColor(kBlue);
     643    h5a.SetLineColor(kBlue);
     644    h0b.SetLineColor(kRed);
     645    h4b.SetLineColor(kRed);
     646    h5b.SetLineColor(kRed);
    621647
    622648    TH1 *hist = fHist.Project3D("xy_fit");
    … …  
    637663    func.SetParName(4, "b");
    638664
     665    func.SetParLimits(3, -1, 1);
     666
    639667    const Int_t nx = fHist.GetXaxis()->GetNbins();
    640668    const Int_t ny = fHist.GetYaxis()->GetNbins();
     669
     670    Double_t maxs=3;
     671    TH1 *result=0;
     672    TF1 *fres=0;
     673
    641674
    642675    TH1 *h=0;
    … …  
    649682            if (h->GetBinContent(1)==0)
    650683                continue;
    651 
    652             func.SetParLimits(3, 0, 1);
    653684
    654685            // First fit a polynom in the off region
    … …  
    659690            func.ReleaseParameter(4);
    660691
    661             h->Fit(&func, "N0Q", "", 35, 80);
    662             *fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
     692            h->Fit(&func, "N0Q", "", 35, 75);
     693            //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
     694
     695            h4a.Fill(func.GetChisquare());
     696            h5a.Fill(func.GetProb());
    663697
    664698            // Now fit a gaus in the on region on top of the polynom
     699            func.SetParameter(0, h->GetBinContent(1)-func.GetParameter(4));
     700            func.SetParameter(2, 80);
     701
     702            //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
     703            //func.SetParLimits(2, 5, 90);
     704
    665705            func.ReleaseParameter(0);
    666706            //func.ReleaseParameter(1);
    667707            func.ReleaseParameter(2);
    668 
    669             func.SetParameter(0, h->GetBinContent(1));
    670             func.SetParameter(2, 10);
    671 
    672             func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
    673             func.SetParLimits(2, 0, 90);
    674 
    675708            func.FixParameter(3, func.GetParameter(3));
    676709            func.FixParameter(4, func.GetParameter(4));
    677710
    678             h->Fit(&func, "N0Q", "", 0, 20);
    679             *fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
     711            func.SetParLimits(2, 0, 80);
     712            h->Fit(&func, "N0Q", "", 0, 35);
     713            //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
    680714
    681715            // Fill results into some histograms
    682             h0.Fill(func.GetParameter(0));
     716            h0a.Fill(func.GetParameter(0));
     717            h0b.Fill(func.GetParameter(4));
    683718            h1.Fill(func.GetParameter(1));
    684719            h2.Fill(func.GetParameter(2));
    685720            h3.Fill(func.GetParameter(3));
    686 
    687             *fLog << dbg << "     " << func.GetChisquare() << " " << func.GetNDF() << " " << func.GetNpx() << " " << func.GetNumberFreeParameters() << " " << func.GetNumberFitPoints() << endl;
     721            h4b.Fill(func.GetChisquare());
     722            h5b.Fill(func.GetProb());
    688723
    689724            const Int_t n = hist->GetBin(ix+1, iy+1);
    690             hist->SetBinContent(n, func.GetParameter(0));
     725            if (func.GetParameter(0)>h->GetBinContent(1)*2 ||
     726                func.GetParameter(2)<2.5 || func.GetParameter(2)>70
     727                /*func.GetProb()<0.005 ||*/)
     728            {
     729                hist->SetBinContent(n, 0);
     730                continue;
     731            }
     732
     733            //*fLog << dbg << "     Chisq=" << func.GetChisquare() << " NDF=" << func.GetNDF()
     734            //    << " Prob=" << func.GetProb() << " FitP=" << func.GetNumberFitPoints() << endl;
     735
     736            Double_t b = FcnI1(15, func.GetParameters());
     737            Double_t s = FcnI2(15, func.GetParameters());
     738
     739            //*fLog << dbg << func.GetParameter(3) << ": " << b << "    ";
     740            //*fLog << dbg << func.GetParameter(0) << ": " << s;
     741            //*fLog << endl;
     742
     743            if (b<0)
     744                continue;
     745
     746            const Double_t sig = Significance(s+b, b);
     747
     748            hist->SetBinContent(n, sig);
     749            h6.Fill(sig);
     750
     751            if (sig>maxs)
     752            {
     753                maxs = sig;
     754                if (result)
     755                {
     756                    delete result;
     757                    delete fres;
     758                }
     759                result=(TH1*)h->Clone("Result \\alpha");
     760                fres  =(TF1*)func.Clone("Result Func");
     761            }
    691762        }
    692763
    693     hist->SetMinimum(0);
    694     hist->SetMaximum(10000);
     764    hist->SetMinimum(-10);
     765    hist->SetMaximum(maxs);
    695766
    696767    TCanvas *c=new TCanvas;
    697768
    698     c->Divide(2,2);
     769    c->Divide(3,2);
    699770    c->cd(1);
    700     h0.DrawCopy();
     771    h0b.DrawCopy();
     772    h0a.DrawCopy("same");
    701773    c->cd(2);
    702774    hist->Draw("colz");
    … …  
    707779    c->cd(4);
    708780    h3.DrawCopy();
    709 }
    710 */
     781    c->cd(5);
     782    h4a.DrawCopy();
     783    h4b.DrawCopy("same");
     784    h6.DrawCopy("same");
     785    c->cd(6);
     786    h5a.DrawCopy();
     787    h5b.DrawCopy("same");
     788
     789    if (result)
     790    {
     791        c=new TCanvas;
     792        result->SetBit(kCanDelete);
     793        result->Draw();
     794
     795        TF1 f1("MyFunc1", fcn, 0, 90, 5);
     796        TF1 f2("MyFunc2", fcn, 0, 90, 5);
     797        f1.SetParameters(fres->GetParameters());
     798        f2.SetParameters(fres->GetParameters());
     799        f2.FixParameter(0, 0);
     800        f2.FixParameter(1, 0);
     801        f2.FixParameter(2, 1);
     802        f1.SetLineColor(kGreen);
     803        f2.SetLineColor(kRed);
     804
     805        f2.DrawCopy("same");
     806        f1.DrawCopy("same");
     807
     808        delete fres;
     809    }
     810}
     811

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

    r3548 r3550  
    3030
    3131    Float_t fAlphaCut;          // Alpha cut
    32 
    3332    Float_t fBgMean;            // Background mean
    3433
    … …  
    5251    TH1 *GetHistByName(const TString name) { return &fHist; }
    5352
     53    void FitSignificance(); //*MENU*
     54
    5455    void SetAlphaCut(Float_t alpha); //*MENU*
    5556    void SetAlphaPlus5()  { SetAlphaCut(fAlphaCut+5); } //*MENU*
    … …  
    6061    void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU*
    6162
    62     //void FitSignificance(); //*MENU*
    63 
    6463    void Paint(Option_t *opt="");
    6564    void Draw(Option_t *option="");
     65
     66    static Double_t Significance(Double_t s, Double_t b);
    6667
    6768    ClassDef(MHFalseSource, 1) //3D-histogram in alpha, x and y
Note: See TracChangeset for help on using the changeset viewer.