Changeset 3557 for trunk

Timestamp:

03/19/04 20:19:03 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHFalseSource.cc (modified) (11 diffs)
• MHFalseSource.h (modified) (1 diff)

trunk/MagicSoft/Mars/mhist/MHFalseSource.cc

@@ -103 +103 @@
 #include <TStyle.h>
 #include <TCanvas.h>
+#include <TPaveText.h>
 
 #include "MGeomCam.h"
@@ -490 +491 @@
 
             TH1 *h = fHist.ProjectionZ("Alpha", maxx, maxx, maxy, maxy);
-            h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma=%.1f)", x, y, s));
+            h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s));
         }
     }
@@ -616 +617 @@
 
 
-void MHFalseSource::FitSignificance()
-{
-    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);
+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);
@@ -647 +648 @@
 
     TH1 *hist = fHist.Project3D("xy_fit");
+    hist->Reset();
+    hist->SetNameTitle("Significance",
+                       Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
+                            sigmax, bgmin, bgmax));
+    hist->SetXTitle(fHist.GetXaxis()->GetTitle());
+    hist->SetYTitle(fHist.GetXaxis()->GetTitle());
 
     //                      xmin, xmax, npar
     TF1 func("MyFunc", fcn, 0,    90,   5);
-
-    TArrayD samx(50);
-    TArrayD samw(50);
-
-    // func.CalcGaussLegendreSamplingPoints(50, samx.GetArray(), samw.GetArray());
-    // func.IntegralFast(50, samx.GetArray(), samw.GetArray(), 0, 15);
+    TArrayD par(5);
 
     func.SetParName(0, "A");
@@ -668 +670 @@
     const Int_t ny = fHist.GetYaxis()->GetNbins();
 
+    Double_t maxalpha0=0;
     Double_t maxs=3;
-    TH1 *result=0;
-    TF1 *fres=0;
-
+
+    Int_t maxx=0;
+    Int_t maxy=0;
 
     TH1 *h=0;
@@ -679 +682 @@
             h = fHist.ProjectionZ("AlphaFit", ix+1, ix+1, iy+1, iy+1);
 
+            const Double_t alpha0 = h->GetBinContent(1);
+
             // Check for the regios which is not filled...
-            if (h->GetBinContent(1)==0)
+            if (alpha0==0)
                 continue;
+
+            if (alpha0>maxalpha0)
+                maxalpha0=alpha0;
 
             // First fit a polynom in the off region
@@ -690 +698 @@
             func.ReleaseParameter(4);
 
-            h->Fit(&func, "N0Q", "", 35, 75);
+            h->Fit(&func, "N0Q", "", bgmin, bgmax);
             //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
 
@@ -697 +705 @@
 
             // 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.SetParameter(0, alpha0-func.GetParameter(4));
+            func.SetParameter(2, sigmax*0.75);
 
             //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
@@ -710 +718 @@
 
             func.SetParLimits(2, 0, 80);
-            h->Fit(&func, "N0Q", "", 0, 35);
+            h->Fit(&func, "N0Q", "", 0, sigmax);
             //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
 
@@ -722 +730 @@
             h5b.Fill(func.GetProb());
 
+            Double_t sig=0;
+
             const Int_t n = hist->GetBin(ix+1, iy+1);
-            if (func.GetParameter(0)>h->GetBinContent(1)*2 ||
-                func.GetParameter(2)<2.5 || func.GetParameter(2)>70
-                /*func.GetProb()<0.005 ||*/)
+            if (!(func.GetParameter(0)>alpha0*2 ||
+                  func.GetParameter(2)<2.5      ||
+                  func.GetParameter(2)>70
+                /*func.GetProb()<0.005 ||*/))
             {
-                hist->SetBinContent(n, 0);
-                continue;
+                const Double_t b = FcnI1(15, func.GetParameters());
+                const Double_t s = FcnI2(15, func.GetParameters());
+
+                sig = Significance(s+b, b);
             }
 
-            //*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);
+            hist->SetBinContent(n, sig==0?1e-15:sig);
+            if (sig!=0)
+                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");
+                maxx = ix+1;
+                maxy = iy+1;
+                for(int i=0; i<par.GetSize(); i++)
+                    par[i] = func.GetParameter(i);
             }
         }
 
-    hist->SetMinimum(-10);
-    hist->SetMaximum(maxs);
+    h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
+    h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
 
     TCanvas *c=new TCanvas;
 
-    c->Divide(3,2);
+    c->Divide(3,2, 0, 0);
     c->cd(1);
+    gPad->SetBorderMode(0);
+    gPad->Divide(1,2, 0, 0);
+    TVirtualPad *p=gPad;
+    p->SetBorderMode(0);
+    p->cd(1);
+    gPad->SetBorderMode(0);
     h0b.DrawCopy();
     h0a.DrawCopy("same");
+    p->cd(2);
+    gPad->SetBorderMode(0);
+    h3.DrawCopy();
     c->cd(2);
+    gPad->SetBorderMode(0);
     hist->Draw("colz");
     hist->SetDirectory(NULL);
     hist->SetBit(kCanDelete);
     c->cd(3);
+    gPad->SetBorderMode(0);
+    h4b.DrawCopy();
+    h4a.DrawCopy("same");
+    h6.DrawCopy("same");
+    c->cd(4);
+    gPad->SetBorderMode(0);
     h2.DrawCopy();
-    c->cd(4);
-    h3.DrawCopy();
+    c->cd(6);
+    gPad->SetBorderMode(0);
+    h5b.DrawCopy();
+    h5a.DrawCopy("same");
+
     c->cd(5);
-    h4a.DrawCopy();
-    h4b.DrawCopy("same");
-    h6.DrawCopy("same");
-    c->cd(6);
-    h5a.DrawCopy();
-    h5b.DrawCopy("same");
-
-    if (result)
-    {
-        c=new TCanvas;
+    gPad->SetBorderMode(0);
+    if (maxx>0 && maxy>0)
+    {
+        const char *title = Form(" \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f) ",
+                                 fHist.GetXaxis()->GetBinCenter(maxx),
+                                 fHist.GetYaxis()->GetBinCenter(maxy), maxs);
+        TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
+        result->SetNameTitle("Result \\alpha", title);
         result->SetBit(kCanDelete);
+        result->SetXTitle("\\alpha [\\circ]");
+        result->SetYTitle("Counts");
         result->Draw();
 
         TF1 f1("MyFunc1", fcn, 0, 90, 5);
         TF1 f2("MyFunc2", fcn, 0, 90, 5);
-        f1.SetParameters(fres->GetParameters());
-        f2.SetParameters(fres->GetParameters());
+        f1.SetParameters(par.GetArray());
+        f2.SetParameters(par.GetArray());
         f2.FixParameter(0, 0);
         f2.FixParameter(1, 0);
@@ -806 +821 @@
         f1.DrawCopy("same");
 
-        delete fres;
-    }
-}
-
+        TPaveText *leg = new TPaveText(0.35, 0.10, 0.90, 0.35, "brNDC");
+        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->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->SetBit(kCanDelete);
+        leg->Draw();
+    }
+
+}

trunk/MagicSoft/Mars/mhist/MHFalseSource.h

@@ -51 +51 @@
     TH1 *GetHistByName(const TString name) { return &fHist; }
 
-    void FitSignificance(); //*MENU*
+    void FitSignificance(Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80); //*MENU*
+    void FitSignificanceStd() { FitSignificance(); } //*MENU*
 
     void SetAlphaCut(Float_t alpha); //*MENU*