Changeset 3550 for trunk/MagicSoft/Mars

Timestamp:

03/19/04 16:07:28 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHFalseSource.cc (modified) (8 diffs)
• MHFalseSource.h (modified) (3 diffs)

trunk/MagicSoft/Mars/mhist/MHFalseSource.cc

@@ -233 +233 @@
 // --------------------------------------------------------------------------
 //
+// 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 +275 @@
 
         MBinning b;
-        b.SetEdges(40, -r, r);
+        b.SetEdges(100, -r, r);
         SetBinning(&fHist, &b, &b, &binsa);
     }
@@ -455 +468 @@
                 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 +598 @@
     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 +663 @@
     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 +682 @@
             if (h->GetBinContent(1)==0)
                 continue;
-
-            func.SetParLimits(3, 0, 1);
 
             // First fit a polynom in the off region
@@ -659 +690 @@
             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 +779 @@
     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;
+    }
+}
+

trunk/MagicSoft/Mars/mhist/MHFalseSource.h

@@ -30 +30 @@
 
     Float_t fAlphaCut;          // Alpha cut
-
     Float_t fBgMean;            // Background mean
 
@@ -52 +51 @@
     TH1 *GetHistByName(const TString name) { return &fHist; }
 
+    void FitSignificance(); //*MENU*
+
     void SetAlphaCut(Float_t alpha); //*MENU*
     void SetAlphaPlus5()  { SetAlphaCut(fAlphaCut+5); } //*MENU*
@@ -60 +61 @@
     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