Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3549) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3550) @@ -233,4 +233,17 @@ // -------------------------------------------------------------------------- // +// Calculate Significance as +// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b +// +Double_t MHFalseSource::Significance(Double_t s, Double_t b) +{ + const Double_t k = b==0 ? 0 : s/b; + const Double_t f = s+k*k*b; + + return f==0 ? 0 : (s-b)/TMath::Sqrt(f); +} + +// -------------------------------------------------------------------------- +// // Set binnings (takes BinningFalseSource) and prepare filling of the // histogram. @@ -262,5 +275,5 @@ MBinning b; - b.SetEdges(40, -r, r); + b.SetEdges(100, -r, r); SetBinning(&fHist, &b, &b, &binsa); } @@ -455,13 +468,6 @@ const Double_t b = h2->GetBinContent(n); - const Double_t k = b==0 ? 0 : s/b; - const Double_t f = s+k*k*b; - - const Double_t sig = f==0 ? 0 : (s-b)/TMath::Sqrt(f); - - //if (b>0 && s>0) - // *fLog << dbg << b << " " << s << endl; - - // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b + const Double_t sig = Significance(s, b); + h4->SetBinContent(n, sig); @@ -592,31 +598,51 @@ return f1 + f2; } -/* -Double_t FcnIntegral(Double_t *arg, Double_t *p) -{ - static const Double_t sqrt2 = TMath::Sqrt(2); + +Double_t FcnI1(Double_t x, Double_t *p) +{ + return (p[3]*x*x/3+p[4])*x; +} +Double_t FcnI2(Double_t x, Double_t *p) +{ + static const Double_t sqrt2 = TMath::Sqrt(2.); static const Double_t sqrt2pi = TMath::Sqrt(TMath::TwoPi()); - const Double_t x = arg[0]; - const Double_t dx = (x-p[1])/p[2]; - const Double_t f1 = p[0]*p[2]*sqrt2pi*TMath::Erf(dx/sqrt2)/2; - const Double_t f2 = p[3]*x*2; - - return f1+f2; -} -*/ -/* + const Double_t f2 = p[0]*p[2]*sqrt2pi*TMath::Erf(dx/sqrt2)/2; + + return f2; +} + + void MHFalseSource::FitSignificance() { - TH1D h0("A", "Parameter A", 50, 0, 10000); - TH1D h1("mu", "Parameter mu", 50, -1, 1); - TH1D h2("sigma", "Parameter sigma", 50, 0, 20); - TH1D h3("b", "Parameter b", 50, 0.001, 0.01); - h0.SetDirectory(NULL); + TH1D h0a("A", "Parameter A", 50, 0, 4000); + TH1D h0b("a", "Parameter a", 50, 0, 4000); + TH1D h1("mu", "Parameter mu", 50, -1, 1); + TH1D h2("sigma", "Parameter sigma", 50, 0, 90); + TH1D h3("b", "Parameter b", 50, -0.1, 0.1); + TH1D h4a("chisq1", "\\chi^{2}", 50, 0, 35); + TH1D h5a("prob1", "Fit probability", 50, 0, 1.1); + TH1D h4b("chisq2", "\\chi^{2}", 50, 0, 35); + TH1D h5b("prob2", "Fit probability", 50, 0, 1.1); + TH1D h6("Signif", "Significance", 50, -20, 20); + h0a.SetDirectory(NULL); + h0b.SetDirectory(NULL); h1.SetDirectory(NULL); h2.SetDirectory(NULL); h3.SetDirectory(NULL); + h4a.SetDirectory(NULL); + h5b.SetDirectory(NULL); + h4a.SetDirectory(NULL); + h5b.SetDirectory(NULL); + h6.SetDirectory(NULL); + + h0a.SetLineColor(kBlue); + h4a.SetLineColor(kBlue); + h5a.SetLineColor(kBlue); + h0b.SetLineColor(kRed); + h4b.SetLineColor(kRed); + h5b.SetLineColor(kRed); TH1 *hist = fHist.Project3D("xy_fit"); @@ -637,6 +663,13 @@ func.SetParName(4, "b"); + func.SetParLimits(3, -1, 1); + const Int_t nx = fHist.GetXaxis()->GetNbins(); const Int_t ny = fHist.GetYaxis()->GetNbins(); + + Double_t maxs=3; + TH1 *result=0; + TF1 *fres=0; + TH1 *h=0; @@ -649,6 +682,4 @@ if (h->GetBinContent(1)==0) continue; - - func.SetParLimits(3, 0, 1); // First fit a polynom in the off region @@ -659,44 +690,85 @@ func.ReleaseParameter(4); - h->Fit(&func, "N0Q", "", 35, 80); - *fLog << dbg << ix << "/" << iy << ": " << func.GetParameter(3) << " " << func.GetParameter(4) << endl; + h->Fit(&func, "N0Q", "", 35, 75); + //*fLog << dbg << ix << "/" << iy << ": " << func.GetParameter(3) << " " << func.GetParameter(4) << endl; + + h4a.Fill(func.GetChisquare()); + h5a.Fill(func.GetProb()); // Now fit a gaus in the on region on top of the polynom + func.SetParameter(0, h->GetBinContent(1)-func.GetParameter(4)); + func.SetParameter(2, 80); + + //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1)); + //func.SetParLimits(2, 5, 90); + func.ReleaseParameter(0); //func.ReleaseParameter(1); func.ReleaseParameter(2); - - func.SetParameter(0, h->GetBinContent(1)); - func.SetParameter(2, 10); - - func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1)); - func.SetParLimits(2, 0, 90); - func.FixParameter(3, func.GetParameter(3)); func.FixParameter(4, func.GetParameter(4)); - h->Fit(&func, "N0Q", "", 0, 20); - *fLog << dbg << " " << func.GetParameter(0) << " " << func.GetParameter(1) << " " << func.GetParameter(2) << endl; + func.SetParLimits(2, 0, 80); + h->Fit(&func, "N0Q", "", 0, 35); + //*fLog << dbg << " " << func.GetParameter(0) << " " << func.GetParameter(1) << " " << func.GetParameter(2) << endl; // Fill results into some histograms - h0.Fill(func.GetParameter(0)); + h0a.Fill(func.GetParameter(0)); + h0b.Fill(func.GetParameter(4)); h1.Fill(func.GetParameter(1)); h2.Fill(func.GetParameter(2)); h3.Fill(func.GetParameter(3)); - - *fLog << dbg << " " << func.GetChisquare() << " " << func.GetNDF() << " " << func.GetNpx() << " " << func.GetNumberFreeParameters() << " " << func.GetNumberFitPoints() << endl; + h4b.Fill(func.GetChisquare()); + h5b.Fill(func.GetProb()); const Int_t n = hist->GetBin(ix+1, iy+1); - hist->SetBinContent(n, func.GetParameter(0)); + if (func.GetParameter(0)>h->GetBinContent(1)*2 || + func.GetParameter(2)<2.5 || func.GetParameter(2)>70 + /*func.GetProb()<0.005 ||*/) + { + hist->SetBinContent(n, 0); + continue; + } + + //*fLog << dbg << " Chisq=" << func.GetChisquare() << " NDF=" << func.GetNDF() + // << " Prob=" << func.GetProb() << " FitP=" << func.GetNumberFitPoints() << endl; + + Double_t b = FcnI1(15, func.GetParameters()); + Double_t s = FcnI2(15, func.GetParameters()); + + //*fLog << dbg << func.GetParameter(3) << ": " << b << " "; + //*fLog << dbg << func.GetParameter(0) << ": " << s; + //*fLog << endl; + + if (b<0) + continue; + + const Double_t sig = Significance(s+b, b); + + hist->SetBinContent(n, sig); + h6.Fill(sig); + + if (sig>maxs) + { + maxs = sig; + if (result) + { + delete result; + delete fres; + } + result=(TH1*)h->Clone("Result \\alpha"); + fres =(TF1*)func.Clone("Result Func"); + } } - hist->SetMinimum(0); - hist->SetMaximum(10000); + hist->SetMinimum(-10); + hist->SetMaximum(maxs); TCanvas *c=new TCanvas; - c->Divide(2,2); + c->Divide(3,2); c->cd(1); - h0.DrawCopy(); + h0b.DrawCopy(); + h0a.DrawCopy("same"); c->cd(2); hist->Draw("colz"); @@ -707,4 +779,33 @@ c->cd(4); h3.DrawCopy(); -} -*/ + c->cd(5); + h4a.DrawCopy(); + h4b.DrawCopy("same"); + h6.DrawCopy("same"); + c->cd(6); + h5a.DrawCopy(); + h5b.DrawCopy("same"); + + if (result) + { + c=new TCanvas; + result->SetBit(kCanDelete); + result->Draw(); + + TF1 f1("MyFunc1", fcn, 0, 90, 5); + TF1 f2("MyFunc2", fcn, 0, 90, 5); + f1.SetParameters(fres->GetParameters()); + f2.SetParameters(fres->GetParameters()); + f2.FixParameter(0, 0); + f2.FixParameter(1, 0); + f2.FixParameter(2, 1); + f1.SetLineColor(kGreen); + f2.SetLineColor(kRed); + + f2.DrawCopy("same"); + f1.DrawCopy("same"); + + delete fres; + } +} + Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3549) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3550) @@ -30,5 +30,4 @@ Float_t fAlphaCut; // Alpha cut - Float_t fBgMean; // Background mean @@ -52,4 +51,6 @@ TH1 *GetHistByName(const TString name) { return &fHist; } + void FitSignificance(); //*MENU* + void SetAlphaCut(Float_t alpha); //*MENU* void SetAlphaPlus5() { SetAlphaCut(fAlphaCut+5); } //*MENU* @@ -60,8 +61,8 @@ void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU* - //void FitSignificance(); //*MENU* - void Paint(Option_t *opt=""); void Draw(Option_t *option=""); + + static Double_t Significance(Double_t s, Double_t b); ClassDef(MHFalseSource, 1) //3D-histogram in alpha, x and y