Changeset 3595 for trunk/MagicSoft
- Timestamp:
- 03/24/04 12:01:41 (21 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mbase/MMath.cc
r3581 r3595 59 59 // alpha = t_on/t_off; // t: observation time 60 60 // 61 // Returns -1 if calculation failed... 62 // 61 63 Double_t MMath::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha) 62 64 { … … 64 66 65 67 if (sum<=0 || alpha<=0) 66 return 0;68 return -1; 67 69 68 70 const Double_t l = s*TMath::Log(s/sum*(alpha+1)/alpha); -
trunk/MagicSoft/Mars/mhist/MHFalseSource.cc
r3593 r3595 263 263 Double_t MHFalseSource::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha) 264 264 { 265 return MMath::SignificanceLiMa(s, b); 265 Double_t lima = MMath::SignificanceLiMa(s, b); 266 return lima<0 ? MMath::Significance(s, b) : lima; 266 267 } 267 268 … … 340 341 hsrc.SetSrcPos(&src); 341 342 342 // This is because a 3D histogram x and y are exchanged... 343 const Int_t nx = fHist.GetNbinsY(); 344 const Int_t ny = fHist.GetNbinsX(); 343 const Int_t nx = fHist.GetNbinsX(); 344 const Int_t ny = fHist.GetNbinsY(); 345 345 Axis_t cx[nx]; 346 346 Axis_t cy[ny]; 347 fHist.Get Yaxis()->GetCenter(cx);348 fHist.Get Xaxis()->GetCenter(cy);347 fHist.GetXaxis()->GetCenter(cx); 348 fHist.GetYaxis()->GetCenter(cy); 349 349 350 350 for (int ix=0; ix<nx; ix++) … … 372 372 const Double_t alpha = hsrc.GetAlpha(); 373 373 374 // This is because a 3D histogram x and y are exchanged... 375 fHist.Fill(cy[iy], cx[ix], TMath::Abs(alpha), w); 374 fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w); 376 375 } 377 376 } … … 400 399 401 400 // Get projection for range 402 TH2D *p = (TH2D*)fHist.Project3D(" xy_off");401 TH2D *p = (TH2D*)fHist.Project3D("yx_off"); 403 402 404 403 // Reset range … … 431 430 432 431 // Get projection for range 433 TH2D *p = (TH2D*)fHist.Project3D(" xy_on");432 TH2D *p = (TH2D*)fHist.Project3D("yx_on"); 434 433 435 434 // Reset range … … 463 462 464 463 padsave->cd(3); 465 if ((h3 = (TH2D*)gPad->FindObject("Alpha_ xy_on")))464 if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) 466 465 ProjectOn(h3); 467 466 468 467 padsave->cd(4); 469 if ((h2 = (TH2D*)gPad->FindObject("Alpha_ xy_off")))468 if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) 470 469 ProjectOff(h2); 471 470 472 471 padsave->cd(2); 473 if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_ xy_sig")))472 if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig"))) 474 473 { 475 474 const Int_t nx = h4->GetXaxis()->GetNbins(); … … 489 488 const Double_t b = h2->GetBinContent(n); 490 489 491 const Double_t sig = Significance (s, b);490 const Double_t sig = SignificanceLiMa(s, b); 492 491 493 492 h4->SetBinContent(n, sig); … … 510 509 // This is because a 3D histogram x and y are vice versa 511 510 // Than for their projections 512 TH1 *h = fHist.ProjectionZ("Alpha", max y, maxy, maxx, maxx);511 TH1 *h = fHist.ProjectionZ("Alpha", maxx, maxx, maxy, maxy); 513 512 h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s)); 514 513 } … … 545 544 gPad->SetBorderMode(0); 546 545 fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); 547 TH1 *h2 = fHist.Project3D(" xy_off");546 TH1 *h2 = fHist.Project3D("yx_off"); 548 547 h2->SetDirectory(NULL); 549 h2->SetXTitle(fHist.Get Yaxis()->GetTitle());550 h2->SetYTitle(fHist.Get Xaxis()->GetTitle());548 h2->SetXTitle(fHist.GetXaxis()->GetTitle()); 549 h2->SetYTitle(fHist.GetYaxis()->GetTitle()); 551 550 h2->Draw("colz"); 552 551 h2->SetBit(kCanDelete); … … 555 554 gPad->SetBorderMode(0); 556 555 fHist.GetZaxis()->SetRangeUser(0,fAlphaCut); 557 TH1 *h3 = fHist.Project3D(" xy_on");556 TH1 *h3 = fHist.Project3D("yx_on"); 558 557 fHist.GetZaxis()->SetRange(0,9999); 559 558 h3->SetDirectory(NULL); 560 h3->SetXTitle(fHist.Get Yaxis()->GetTitle());561 h3->SetYTitle(fHist.Get Xaxis()->GetTitle());559 h3->SetXTitle(fHist.GetXaxis()->GetTitle()); 560 h3->SetYTitle(fHist.GetYaxis()->GetTitle()); 562 561 h3->Draw("colz"); 563 562 h3->SetBit(kCanDelete); … … 566 565 gPad->SetBorderMode(0); 567 566 fHist.GetZaxis()->SetRange(0,0); 568 TH1 *h4 = fHist.Project3D(" xy_sig"); // Do this to get the correct binning....567 TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning.... 569 568 fHist.GetZaxis()->SetRange(0,9999); 570 569 h4->SetTitle("Significance"); 571 570 h4->SetDirectory(NULL); 572 h4->SetXTitle(fHist.Get Yaxis()->GetTitle());573 h4->SetYTitle(fHist.Get Xaxis()->GetTitle());571 h4->SetXTitle(fHist.GetXaxis()->GetTitle()); 572 h4->SetYTitle(fHist.GetYaxis()->GetTitle()); 574 573 h4->Draw("colz"); 575 574 h4->SetBit(kCanDelete); … … 659 658 //h5b.SetLineColor(kRed); 660 659 661 TH1 *hist = fHist.Project3D(" xy_fit");660 TH1 *hist = fHist.Project3D("yx_fit"); 662 661 hist->SetDirectory(0); 663 TH1 *hists = fHist.Project3D(" xy_fit");662 TH1 *hists = fHist.Project3D("yx_fit"); 664 663 hists->SetDirectory(0); 665 TH1 *histb = fHist.Project3D(" xy_fit");664 TH1 *histb = fHist.Project3D("yx_fit"); 666 665 histb->SetDirectory(0); 667 666 hist->Reset(); … … 673 672 hists->SetNameTitle("Excess", Form("Number of excess events for \\alpha<%.0f\\circ", sigint)); 674 673 histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint)); 675 hist->SetXTitle(fHist.Get Yaxis()->GetTitle());676 hists->SetXTitle(fHist.Get Yaxis()->GetTitle());677 histb->SetXTitle(fHist.Get Yaxis()->GetTitle());678 hist->SetYTitle(fHist.Get Xaxis()->GetTitle());679 hists->SetYTitle(fHist.Get Xaxis()->GetTitle());680 histb->SetYTitle(fHist.Get Xaxis()->GetTitle());674 hist->SetXTitle(fHist.GetXaxis()->GetTitle()); 675 hists->SetXTitle(fHist.GetXaxis()->GetTitle()); 676 histb->SetXTitle(fHist.GetXaxis()->GetTitle()); 677 hist->SetYTitle(fHist.GetYaxis()->GetTitle()); 678 hists->SetYTitle(fHist.GetYaxis()->GetTitle()); 679 histb->SetYTitle(fHist.GetYaxis()->GetTitle()); 681 680 682 681 const Double_t w = fHist.GetZaxis()->GetBinWidth(1); … … 728 727 // This is because a 3D histogram x and y are vice versa 729 728 // Than for their projections 730 h = fHist.ProjectionZ("AlphaFit", i y, iy, ix, ix);729 h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy); 731 730 732 731 const Double_t alpha0 = h->GetBinContent(1); … … 810 809 811 810 const Double_t b = func.Integral(0, sigint)/w; 812 const Double_t sig = Significance(s, b); 813 814 // This is because a 3D histogram x and y are exchanged... 811 const Double_t sig = SignificanceLiMa(s, b); 812 815 813 const Int_t n = hist->GetBin(ix, iy); 816 814 hists->SetBinContent(n, s-b); … … 894 892 hist->GetYaxis()->GetBinCenter(maxy), maxs); 895 893 896 TH1 *result = fHist.ProjectionZ("AlphaFit", max y, maxy, maxx, maxx);894 TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); 897 895 result->SetDirectory(NULL); 898 896 result->SetNameTitle("Result \\alpha", title);
Note:
See TracChangeset
for help on using the changeset viewer.