Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3556) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3557) @@ -103,4 +103,5 @@ #include #include +#include #include "MGeomCam.h" @@ -490,5 +491,5 @@ 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,16 +617,16 @@ -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,13 +648,14 @@ TH1 *hist = fHist.Project3D("xy_fit"); + hist->Reset(); + hist->SetNameTitle("Significance", + Form("Fit Region: Signal<%.1f\\circ, %.1f\\circSetXTitle(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,8 +670,9 @@ 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,7 +682,12 @@ 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,5 +698,5 @@ 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,6 +705,6 @@ // 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,5 +718,5 @@ 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,79 +730,86 @@ 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; iSetMinimum(-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,6 +821,17 @@ 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(); + } + +} Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3556) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3557) @@ -51,5 +51,6 @@ 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*