Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc
===================================================================
--- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc	(revision 3557)
+++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc	(revision 3568)
@@ -95,4 +95,8 @@
 //  - a more clever (and faster) algorithm to fill the histogram, eg by
 //    calculating alpha once and fill the two coils around the mean
+//  - more drawing options...
+//  - Move Significance() to a more 'general' place and implement
+//    also different algorithms like (Li/Ma)
+//  - implement fit for best alpha distribution -- online (Paint)
 //
 //////////////////////////////////////////////////////////////////////////////
@@ -101,7 +105,9 @@
 #include <TF1.h>
 #include <TH2.h>
+#include <TGraph.h>
 #include <TStyle.h>
 #include <TCanvas.h>
 #include <TPaveText.h>
+#include <TStopwatch.h>
 
 #include "MGeomCam.h"
@@ -236,4 +242,7 @@
 // Calculate Significance as
 // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
+//
+// s: total number of events in signal region
+// b: number of background events in signal region
 // 
 Double_t MHFalseSource::Significance(Double_t s, Double_t b)
@@ -245,4 +254,37 @@
 }
 
+// --------------------------------------------------------------------------
+//
+//  calculates the significance according to Li & Ma
+//  ApJ 272 (1983) 317
+//
+//    s: total number of events in signal region
+//
+//    Double_t fGamma;   // Nbg = Non - gamma * Noff
+//  - the effective number of background events (fNoff), and fGamma :
+//
+//    fGamma = b / fNoff;
+//    fGamma = fdNbg / sqrt(fNoff);
+//    fGamma = fdNbg*fdNbg / fNbg;
+//    fNoff  =  b*b  / (fdNbg*fdNbg);
+//     Double_t fNbg;     // number of background events in signal region
+//    b = fNOff *fGamma;
+/*
+Double_t MHFindSignificance::SignificanceLiMa(Double_t non, Double_t noff, Double_t gamma)
+{
+    if (gamma <= 0.0  ||  non <= 0.0  ||  noff <= 0.0)
+    {
+        *siglima = 0.0;
+        return kFALSE;
+    }
+
+    Double_t help1 = TMath::Log( (gamma+1)*s  / (gamma*(s+noff)) );
+    Double_t help2 = TMath::Log( (gamma+1)*noff / (       s+noff ) );
+
+    Double_t siglima = TMath::Sqrt((s*help1+noff*help2)*2);
+
+    return non<gamma*noff ? -siglima : siglima;
+}
+*/
 // --------------------------------------------------------------------------
 //
@@ -276,5 +318,5 @@
 
         MBinning b;
-        b.SetEdges(100, -r, r);
+        b.SetEdges(20, -r, r);
         SetBinning(&fHist, &b, &b, &binsa);
     }
@@ -311,10 +353,7 @@
     Double_t rho = 0;
     if (fTime && fObservatory && fPointPos)
-    {
-        const Double_t ra  = fPointPos->GetRa();
-        const Double_t dec = fPointPos->GetDec();
-
-        rho = MAstroSky2Local(*fTime, *fObservatory).RotationAngle(ra, dec);
-    }
+        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
+    //if (fPointPos)
+    //    rho = fPointPos->RotationAngle(*fObservatory);
 
     MSrcPosCam src;
@@ -554,5 +593,4 @@
     h4->Draw("colz");
     h4->SetBit(kCanDelete);
-
 }
 
@@ -587,5 +625,5 @@
     gPad->cd();
 }
-
+/*
 Double_t fcn(Double_t *arg, Double_t *p)
 {
@@ -595,5 +633,5 @@
 
     const Double_t f1 = p[0]*TMath::Exp(-0.5*dx*dx);
-    const Double_t f2 = p[3]*x*x + p[4];
+    const Double_t f2 = p[3] + p[5]*x*x;
 
     return f1 + f2;
@@ -602,5 +640,5 @@
 Double_t FcnI1(Double_t x, Double_t *p)
 {
-    return (p[3]*x*x/3+p[4])*x;
+    return (p[5]*x*x/3+p[3])*x;
 }
 Double_t FcnI2(Double_t x, Double_t *p)
@@ -615,54 +653,212 @@
     return f2;
 }
-
-
-void MHFalseSource::FitSignificance(Float_t sigmax, Float_t bgmin, Float_t bgmax)
-{
-    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} (red, green) / significance (black)",       50,  0, 35);
-    TH1D h5a("prob1",  "Fit probability",  50,  0, 1.1);
-    TH1D h4b("chisq2", "\\chi^{2} (red, green) / significance (black)",       50,  0, 35);
-    TH1D h5b("prob",  "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);
+*/
+/*
+class MHSignificance : public MH
+{
+private:
+    TH1D fHist;
+
+    MParameterD *fParam;
+
+public:
+    MHSignificance() : fParam(0)
+    {
+        fHist.SetName("Alpha");
+        fHist.SetTitle("Distribution of \\alpha");
+        fHist.SetXTitle("\\alpha [\\circ]");
+        fHist.SetYTitle("Counts");
+    }
+    Int_t SetupFill(MParList *p)
+    {
+        fHist.Reset();
+
+        fParam = (MParameterD*)p->FindCreateObj("Significance", "MParameterD");
+        if (fParam)
+            return kFALSE;
+
+        return kTRUE;
+    }
+    Int_t Process(MParContainer *p, Double_t w=1)
+    {
+        MHillasSrc *hil = dynamic_cast<MHillasSrc*>(p);
+        if (!hil)
+        {
+            *fLog << err << dbginf << "Got no MHillasSrc as argument of Fill()..." << endl;
+            return kFALSE;
+        }
+
+        fHist->Fill(hil->GetAlpha(), w);
+
+        return kTRUE;
+    }
+    Int_t Finalize()
+    {
+        if (fHist.GetEntries()==0)
+        {
+            *fLog << err << "Histogram empty." << endl;
+            return kFALSE;
+        }
+
+        Float_t sigmax=15;
+        Float_t bgmin =45;
+        Float_t bgmax =80;
+
+        fHist.SetNameTitle("Significance",
+                           Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
+                                sigmax, bgmin, bgmax));
+
+        // Implementing the function yourself is not faster at all!
+        TF1 func("gaus(0) + pol2(3)", fcn, 0, 90, 6);
+        TArrayD maxpar(func.GetNpar());
+
+        func.FixParameter(1, 0);
+        func.FixParameter(4, 0);
+        func.SetParLimits(3, -1, 1);
+
+        const Double_t alpha0 = fHist.GetBinContent(1);
+
+        // First fit a polynom in the off region
+        func.FixParameter(0, 0);
+        func.FixParameter(2, 1);
+        func.ReleaseParameter(3);
+        func.ReleaseParameter(5);
+
+        h->Fit(&func, "N0Q", "", bgmin, bgmax);
+
+        // Now fit a gaus in the on region on top of the polynom
+        func.SetParameter(0, alpha0-func.GetParameter(4));
+        func.SetParameter(2, sigmax*0.75);
+
+        func.ReleaseParameter(0);
+        func.ReleaseParameter(2);
+        func.FixParameter(3, func.GetParameter(3));
+        func.FixParameter(5, func.GetParameter(5));
+
+        func.SetParLimits(2, 0, 80);
+        h->Fit(&func, "N0Q", "", 0, sigmax);
+
+        TArrayD p(func.GetNpar(), func.GetParameters());
+
+        Double_t sig=0;
+
+        const Int_t n = hist->GetBin(ix+1, iy+1);
+        if (!(func.GetParameter(0)>alpha0*2 ||
+              func.GetParameter(2)<2.5      ||
+              func.GetParameter(2)>70))
+        {
+            // Implementing the integral as analytical function
+            // gives the same result in the order of 10e-5
+            // and it is not faster at all...
+            const Double_t s = func.Integral(0, 15);
+
+            func.SetParameter(0, 0);
+            func.SetParameter(2, 1);
+
+            const Double_t b = func.Integral(0, 15);
+
+            sig = Significance(s, b);
+        }
+
+        fParam->SetValue(sig);
+        fParam->SetReadyToSave();
+        return kTRUE;
+    }
+    ClassDef(MHSignificance, 0)
+};
+*/
+
+
+// --------------------------------------------------------------------------
+//
+// This is a preliminary implementation of a alpha-fit procedure for
+// all possible source positions. It will be moved into its own
+// more powerfull class soon.
+//
+// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
+//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
+// or
+//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
+//
+// Parameter [1] is fixed to 0 while the alpha peak should be
+// symmetric around alpha=0.
+//
+// Parameter [4] is fixed to 0 because the first derivative at
+// alpha=0 should be 0, too.
+//
+// In a first step the background is fitted between bgmin and bgmax,
+// while the parameters [0]=0 and [2]=1 are fixed.
+//
+// In a second step the signal region (alpha<sigmax) is fittet using
+// the whole function with parameters [1], [3], [4] and [5] fixed.
+//
+// The number of excess and background events are calculated as
+//   s = int(0, sigint, gaus(0)+pol2(3))
+//   b = int(0, sigint,         pol2(3))
+//
+// The Significance is calculated using the Significance() member
+// function.
+//
+void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
+{
+    TH1D h0a("A",          "", 50,   0, 4000);
+    TH1D h4a("chisq1",     "", 50,   0,   35);
+    //TH1D h5a("prob1",      "", 50,   0,  1.1);
+    TH1D h6("signifcance", "", 50, -20,   20);
+
+    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 h0b("a",         "Parameter a (red), A (blue)", 50, 0, 4000);
+    TH1D h4b("\\chi^{2}", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35);
+    //TH1D h5b("prob",      "Fit probability: Bg(red), F(blue)", 50, 0, 1.1);
 
     h0a.SetLineColor(kBlue);
     h4a.SetLineColor(kBlue);
-    h5a.SetLineColor(kBlue);
+    //h5a.SetLineColor(kBlue);
     h0b.SetLineColor(kRed);
     h4b.SetLineColor(kRed);
-    h5b.SetLineColor(kRed);
-
-    TH1 *hist = fHist.Project3D("xy_fit");
+    //h5b.SetLineColor(kRed);
+
+    TH1 *hist  = fHist.Project3D("xy_fit");
+    hist->SetDirectory(0);
+    TH1 *hists = fHist.Project3D("xy_fit");
+    hists->SetDirectory(0);
+    TH1 *histb = fHist.Project3D("xy_fit");
+    histb->SetDirectory(0);
     hist->Reset();
+    hists->Reset();
+    histb->Reset();
     hist->SetNameTitle("Significance",
                        Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
                             sigmax, bgmin, bgmax));
+    hists->SetNameTitle("Excess",     Form("Number of excess events for \\alpha<%.0f\\circ", sigint));
+    histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint));
     hist->SetXTitle(fHist.GetXaxis()->GetTitle());
+    hists->SetXTitle(fHist.GetXaxis()->GetTitle());
+    histb->SetXTitle(fHist.GetXaxis()->GetTitle());
     hist->SetYTitle(fHist.GetXaxis()->GetTitle());
+    hists->SetYTitle(fHist.GetXaxis()->GetTitle());
+    histb->SetYTitle(fHist.GetXaxis()->GetTitle());
 
     //                      xmin, xmax, npar
-    TF1 func("MyFunc", fcn, 0,    90,   5);
-    TArrayD par(5);
-
-    func.SetParName(0, "A");
-    func.SetParName(1, "mu");
-    func.SetParName(2, "sigma");
-    func.SetParName(3, "a");
-    func.SetParName(4, "b");
-
+    //TF1 func("MyFunc", fcn, 0, 90, 6);
+    // Implementing the function yourself is only about 5% faster
+    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
+    TArrayD maxpar(func.GetNpar());
+
+    /*
+     func.SetParName(0, "A");
+     func.SetParName(1, "mu");
+     func.SetParName(2, "sigma");
+     func.SetParName(3, "a");
+     func.SetParName(4, "b");
+     func.SetParName(5, "c");
+    */
+
+    func.FixParameter(1, 0);
+    func.FixParameter(4, 0);
+    func.SetParLimits(2, 0, 90);
     func.SetParLimits(3, -1, 1);
 
@@ -676,11 +872,22 @@
     Int_t maxy=0;
 
+    TStopwatch clk;
+    clk.Start();
+
+    *fLog << inf;
+    *fLog << "Signal fit:     alpha < " << sigmax << endl;
+    *fLog << "Integration:    alpha < " << sigint << endl;
+    *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
+    *fLog << "Polynom order:  " << (int)polynom << endl;
+    *fLog << "Fitting False Source Plot..." << flush;
+
     TH1 *h=0;
-    for (int ix=0; ix<nx; ix++)
-        for (int iy=0; iy<ny; iy++)
+    for (int ix=1; ix<=nx; ix++)
+        for (int iy=1; iy<=ny; iy++)
         {
-            h = fHist.ProjectionZ("AlphaFit", ix+1, ix+1, iy+1, iy+1);
+            h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
 
             const Double_t alpha0 = h->GetBinContent(1);
+            const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);
 
             // Check for the regios which is not filled...
@@ -693,8 +900,9 @@
             // First fit a polynom in the off region
             func.FixParameter(0, 0);
-            func.FixParameter(1, 0);
             func.FixParameter(2, 1);
             func.ReleaseParameter(3);
-            func.ReleaseParameter(4);
+
+            for (int i=5; i<func.GetNpar(); i++)
+                func.ReleaseParameter(i);
 
             h->Fit(&func, "N0Q", "", bgmin, bgmax);
@@ -702,9 +910,5 @@
 
             h4a.Fill(func.GetChisquare());
-            h5a.Fill(func.GetProb());
-
-            // Now fit a gaus in the on region on top of the polynom
-            func.SetParameter(0, alpha0-func.GetParameter(4));
-            func.SetParameter(2, sigmax*0.75);
+            //h5a.Fill(func.GetProb());
 
             //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
@@ -715,34 +919,57 @@
             func.ReleaseParameter(2);
             func.FixParameter(3, func.GetParameter(3));
-            func.FixParameter(4, func.GetParameter(4));
-
-            func.SetParLimits(2, 0, 80);
+            for (int i=5; i<func.GetNpar(); i++)
+                func.FixParameter(i, func.GetParameter(i));
+
+            // Do not allow signals smaller than the background
+            const Double_t A  = alpha0-func.GetParameter(3);
+            const Double_t dA = TMath::Abs(A);
+            func.SetParLimits(0, -dA*4, dA*4);
+            func.SetParLimits(2, 0, 90);
+
+            // Now fit a gaus in the on region on top of the polynom
+            func.SetParameter(0, A);
+            func.SetParameter(2, sigmax*0.75);
+
             h->Fit(&func, "N0Q", "", 0, sigmax);
             //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
 
+            TArrayD p(func.GetNpar(), func.GetParameters());
+
             // Fill results into some histograms
-            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));
+            h0a.Fill(p[0]);
+            h0b.Fill(p[3]);
+            h1.Fill(p[1]);
+            h2.Fill(p[2]);
+            if (polynom>1)
+                h3.Fill(p[5]);
             h4b.Fill(func.GetChisquare());
-            h5b.Fill(func.GetProb());
-
-            Double_t sig=0;
-
-            const Int_t n = hist->GetBin(ix+1, iy+1);
-            if (!(func.GetParameter(0)>alpha0*2 ||
-                  func.GetParameter(2)<2.5      ||
-                  func.GetParameter(2)>70
-                /*func.GetProb()<0.005 ||*/))
+            //h5b.Fill(func.GetProb());
+
+            // Implementing the integral as analytical function
+            // gives the same result in the order of 10e-5
+            // and it is not faster at all...
+            //const Bool_t ok = /*p[0]>=0 && /*p[0]<alpha0*2 &&*/ p[2]>1.75;// && p[2]<88.5;
+            /*
+            if (p[0]<-fabs(alpha0-p[3])*7 && p[2]<alphaw/3)
             {
-                const Double_t b = FcnI1(15, func.GetParameters());
-                const Double_t s = FcnI2(15, func.GetParameters());
-
-                sig = Significance(s+b, b);
+                func.SetParameter(0, alpha0-p[3]);
+                cout << p[2] << " " << p[0] << " " << alpha0-p[3] << endl;
             }
-
-            hist->SetBinContent(n, sig==0?1e-15:sig);
+            */
+
+            const Double_t s = func.Integral(0, sigint);
+
+            func.SetParameter(0, 0);
+            func.SetParameter(2, 1);
+
+            const Double_t b   = func.Integral(0, sigint);
+            const Double_t sig = Significance(s, b);
+
+            const Int_t n = hist->GetBin(ix, iy);
+            hists->SetBinContent(n, s-b);
+            histb->SetBinContent(n, b);
+
+            hist->SetBinContent(n, sig);
             if (sig!=0)
                 h6.Fill(sig);
@@ -751,14 +978,21 @@
             {
                 maxs = sig;
-                maxx = ix+1;
-                maxy = iy+1;
-                for(int i=0; i<par.GetSize(); i++)
-                    par[i] = func.GetParameter(i);
+                maxx = ix;
+                maxy = iy;
+                maxpar = p;
             }
         }
 
+    *fLog << inf << "done." << endl;
+
     h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
     h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
 
+    //hists->SetMinimum(GetMinimumGT(*hists));
+    histb->SetMinimum(GetMinimumGT(*histb));
+
+    clk.Stop();
+    clk.Print();
+
     TCanvas *c=new TCanvas;
 
@@ -766,5 +1000,16 @@
     c->cd(1);
     gPad->SetBorderMode(0);
-    gPad->Divide(1,2, 0, 0);
+    hists->Draw("colz");
+    hists->SetBit(kCanDelete);
+    c->cd(2);
+    gPad->SetBorderMode(0);
+    hist->Draw("colz");
+    hist->SetBit(kCanDelete);
+    c->cd(3);
+    gPad->SetBorderMode(0);
+    histb->Draw("colz");
+    histb->SetBit(kCanDelete);
+    c->cd(4);
+    gPad->Divide(1,3, 0, 0);
     TVirtualPad *p=gPad;
     p->SetBorderMode(0);
@@ -776,22 +1021,21 @@
     gPad->SetBorderMode(0);
     h3.DrawCopy();
-    c->cd(2);
-    gPad->SetBorderMode(0);
-    hist->Draw("colz");
-    hist->SetDirectory(NULL);
-    hist->SetBit(kCanDelete);
-    c->cd(3);
+    p->cd(3);
+    gPad->SetBorderMode(0);
+    h2.DrawCopy();
+    c->cd(6);
+    gPad->Divide(1,2, 0, 0);
+    TVirtualPad *q=gPad;
+    q->SetBorderMode(0);
+    q->cd(1);
+    gPad->SetBorderMode(0);
     gPad->SetBorderMode(0);
     h4b.DrawCopy();
     h4a.DrawCopy("same");
     h6.DrawCopy("same");
-    c->cd(4);
-    gPad->SetBorderMode(0);
-    h2.DrawCopy();
-    c->cd(6);
-    gPad->SetBorderMode(0);
-    h5b.DrawCopy();
-    h5a.DrawCopy("same");
-
+    q->cd(2);
+    gPad->SetBorderMode(0);
+    //h5b.DrawCopy();
+    //h5a.DrawCopy("same");
     c->cd(5);
     gPad->SetBorderMode(0);
@@ -801,5 +1045,7 @@
                                  fHist.GetXaxis()->GetBinCenter(maxx),
                                  fHist.GetYaxis()->GetBinCenter(maxy), maxs);
+
         TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
+        result->SetDirectory(NULL);
         result->SetNameTitle("Result \\alpha", title);
         result->SetBit(kCanDelete);
@@ -808,8 +1054,8 @@
         result->Draw();
 
-        TF1 f1("MyFunc1", fcn, 0, 90, 5);
-        TF1 f2("MyFunc2", fcn, 0, 90, 5);
-        f1.SetParameters(par.GetArray());
-        f2.SetParameters(par.GetArray());
+        TF1 f1("", func.GetExpFormula(), 0, 90);
+        TF1 f2("", func.GetExpFormula(), 0, 90);
+        f1.SetParameters(maxpar.GetArray());
+        f2.SetParameters(maxpar.GetArray());
         f2.FixParameter(0, 0);
         f2.FixParameter(1, 0);
@@ -824,14 +1070,52 @@
         leg->SetBorderSize(1);
         leg->SetTextSize(0.04);
-        leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
+        leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
+        //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
         leg->AddLine(0, 0.65, 0, 0.65);
-        leg->AddText(0.06, 0.54, Form("A=%.2f", par[0]))->SetTextAlign(11);
-        leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", par[2]))->SetTextAlign(11);
-        leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fixed)", par[1]))->SetTextAlign(11);
-        leg->AddText(0.60, 0.54, Form("a=%.2f", par[3]))->SetTextAlign(11);
-        leg->AddText(0.60, 0.34, Form("b=%.2f", par[4]))->SetTextAlign(11);
+        leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
+        leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
+        leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
+        leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
+        leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
+        if (polynom>1)
+            leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
         leg->SetBit(kCanDelete);
         leg->Draw();
-    }
-
-}
+
+        q->cd(2);
+
+        TGraph *g = new TGraph;
+        g->SetPoint(0, 0, 0);
+
+        Int_t max=0;
+        Float_t maxsig=0;
+        for (int i=1; i<89; i++)
+        {
+            const Double_t s = f1.Integral(0, (float)i);
+            const Double_t b = f2.Integral(0, (float)i);
+
+            const Double_t sig = Significance(s, b);
+
+            g->SetPoint(g->GetN(), i, sig);
+
+            if (sig>maxsig)
+            {
+                max = i;
+                maxsig = sig;
+            }
+        }
+
+        g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
+        g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
+        g->GetHistogram()->SetYTitle("Significance");
+        g->SetBit(kCanDelete);
+        g->Draw("AP");
+
+        leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
+        leg->SetBorderSize(1);
+        leg->SetTextSize(0.1);
+        leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max));
+        leg->SetBit(kCanDelete);
+        leg->Draw();
+    }
+}
Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h
===================================================================
--- trunk/MagicSoft/Mars/mhist/MHFalseSource.h	(revision 3557)
+++ trunk/MagicSoft/Mars/mhist/MHFalseSource.h	(revision 3568)
@@ -51,5 +51,5 @@
     TH1 *GetHistByName(const TString name) { return &fHist; }
 
-    void FitSignificance(Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80); //*MENU*
+    void FitSignificance(Float_t sigint=15, Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80, Byte_t polynom=2); //*MENU*
     void FitSignificanceStd() { FitSignificance(); } //*MENU*
 
