Index: trunk/MagicSoft/Mars/mbase/MMath.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MMath.cc (revision 3594) +++ trunk/MagicSoft/Mars/mbase/MMath.cc (revision 3595) @@ -59,4 +59,6 @@ // alpha = t_on/t_off; // t: observation time // +// Returns -1 if calculation failed... +// Double_t MMath::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha) { @@ -64,5 +66,5 @@ if (sum<=0 || alpha<=0) - return 0; + return -1; const Double_t l = s*TMath::Log(s/sum*(alpha+1)/alpha); Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3594) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3595) @@ -263,5 +263,6 @@ Double_t MHFalseSource::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha) { - return MMath::SignificanceLiMa(s, b); + Double_t lima = MMath::SignificanceLiMa(s, b); + return lima<0 ? MMath::Significance(s, b) : lima; } @@ -340,11 +341,10 @@ hsrc.SetSrcPos(&src); - // This is because a 3D histogram x and y are exchanged... - const Int_t nx = fHist.GetNbinsY(); - const Int_t ny = fHist.GetNbinsX(); + const Int_t nx = fHist.GetNbinsX(); + const Int_t ny = fHist.GetNbinsY(); Axis_t cx[nx]; Axis_t cy[ny]; - fHist.GetYaxis()->GetCenter(cx); - fHist.GetXaxis()->GetCenter(cy); + fHist.GetXaxis()->GetCenter(cx); + fHist.GetYaxis()->GetCenter(cy); for (int ix=0; ixcd(3); - if ((h3 = (TH2D*)gPad->FindObject("Alpha_xy_on"))) + if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) ProjectOn(h3); padsave->cd(4); - if ((h2 = (TH2D*)gPad->FindObject("Alpha_xy_off"))) + if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) ProjectOff(h2); padsave->cd(2); - if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_xy_sig"))) + if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig"))) { const Int_t nx = h4->GetXaxis()->GetNbins(); @@ -489,5 +488,5 @@ const Double_t b = h2->GetBinContent(n); - const Double_t sig = Significance(s, b); + const Double_t sig = SignificanceLiMa(s, b); h4->SetBinContent(n, sig); @@ -510,5 +509,5 @@ // This is because a 3D histogram x and y are vice versa // Than for their projections - TH1 *h = fHist.ProjectionZ("Alpha", maxy, maxy, maxx, maxx); + TH1 *h = fHist.ProjectionZ("Alpha", maxx, maxx, maxy, maxy); h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s)); } @@ -545,8 +544,8 @@ gPad->SetBorderMode(0); fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); - TH1 *h2 = fHist.Project3D("xy_off"); + TH1 *h2 = fHist.Project3D("yx_off"); h2->SetDirectory(NULL); - h2->SetXTitle(fHist.GetYaxis()->GetTitle()); - h2->SetYTitle(fHist.GetXaxis()->GetTitle()); + h2->SetXTitle(fHist.GetXaxis()->GetTitle()); + h2->SetYTitle(fHist.GetYaxis()->GetTitle()); h2->Draw("colz"); h2->SetBit(kCanDelete); @@ -555,9 +554,9 @@ gPad->SetBorderMode(0); fHist.GetZaxis()->SetRangeUser(0,fAlphaCut); - TH1 *h3 = fHist.Project3D("xy_on"); + TH1 *h3 = fHist.Project3D("yx_on"); fHist.GetZaxis()->SetRange(0,9999); h3->SetDirectory(NULL); - h3->SetXTitle(fHist.GetYaxis()->GetTitle()); - h3->SetYTitle(fHist.GetXaxis()->GetTitle()); + h3->SetXTitle(fHist.GetXaxis()->GetTitle()); + h3->SetYTitle(fHist.GetYaxis()->GetTitle()); h3->Draw("colz"); h3->SetBit(kCanDelete); @@ -566,10 +565,10 @@ gPad->SetBorderMode(0); fHist.GetZaxis()->SetRange(0,0); - TH1 *h4 = fHist.Project3D("xy_sig"); // Do this to get the correct binning.... + TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning.... fHist.GetZaxis()->SetRange(0,9999); h4->SetTitle("Significance"); h4->SetDirectory(NULL); - h4->SetXTitle(fHist.GetYaxis()->GetTitle()); - h4->SetYTitle(fHist.GetXaxis()->GetTitle()); + h4->SetXTitle(fHist.GetXaxis()->GetTitle()); + h4->SetYTitle(fHist.GetYaxis()->GetTitle()); h4->Draw("colz"); h4->SetBit(kCanDelete); @@ -659,9 +658,9 @@ //h5b.SetLineColor(kRed); - TH1 *hist = fHist.Project3D("xy_fit"); + TH1 *hist = fHist.Project3D("yx_fit"); hist->SetDirectory(0); - TH1 *hists = fHist.Project3D("xy_fit"); + TH1 *hists = fHist.Project3D("yx_fit"); hists->SetDirectory(0); - TH1 *histb = fHist.Project3D("xy_fit"); + TH1 *histb = fHist.Project3D("yx_fit"); histb->SetDirectory(0); hist->Reset(); @@ -673,10 +672,10 @@ 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.GetYaxis()->GetTitle()); - hists->SetXTitle(fHist.GetYaxis()->GetTitle()); - histb->SetXTitle(fHist.GetYaxis()->GetTitle()); - hist->SetYTitle(fHist.GetXaxis()->GetTitle()); - hists->SetYTitle(fHist.GetXaxis()->GetTitle()); - histb->SetYTitle(fHist.GetXaxis()->GetTitle()); + hist->SetXTitle(fHist.GetXaxis()->GetTitle()); + hists->SetXTitle(fHist.GetXaxis()->GetTitle()); + histb->SetXTitle(fHist.GetXaxis()->GetTitle()); + hist->SetYTitle(fHist.GetYaxis()->GetTitle()); + hists->SetYTitle(fHist.GetYaxis()->GetTitle()); + histb->SetYTitle(fHist.GetYaxis()->GetTitle()); const Double_t w = fHist.GetZaxis()->GetBinWidth(1); @@ -728,5 +727,5 @@ // This is because a 3D histogram x and y are vice versa // Than for their projections - h = fHist.ProjectionZ("AlphaFit", iy, iy, ix, ix); + h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy); const Double_t alpha0 = h->GetBinContent(1); @@ -810,7 +809,6 @@ const Double_t b = func.Integral(0, sigint)/w; - const Double_t sig = Significance(s, b); - - // This is because a 3D histogram x and y are exchanged... + const Double_t sig = SignificanceLiMa(s, b); + const Int_t n = hist->GetBin(ix, iy); hists->SetBinContent(n, s-b); @@ -894,5 +892,5 @@ hist->GetYaxis()->GetBinCenter(maxy), maxs); - TH1 *result = fHist.ProjectionZ("AlphaFit", maxy, maxy, maxx, maxx); + TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); result->SetDirectory(NULL); result->SetNameTitle("Result \\alpha", title);