Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3667) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3682) @@ -46,4 +46,10 @@ // pads) you can change the AlphaCut 'online' // +// Each Z-Projection (Alpha-histogram) is scaled such, that the number +// of entries fits the maximum number of entries in all Z-Projections. +// This should correct for the different acceptance in the center +// and at the vorder of the camera. This, however, produces more noise +// at the border. +// // Here is a slightly simplified version of the algorithm: // ------------------------------------------------------ @@ -139,6 +145,6 @@ // MHFalseSource::MHFalseSource(const char *name, const char *title) - : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), - fAlphaCut(12.5), fBgMean(55), fDistMin(-1), fDistMax(-1) + : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5), + fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1) { // @@ -155,4 +161,15 @@ fHist.SetYTitle("y [\\circ]"); fHist.SetZTitle("\\alpha [\\circ]"); +} + +void MHFalseSource::MakeSymmetric(TH1 *h) +{ + const Float_t min = TMath::Abs(h->GetMinimum()); + const Float_t max = TMath::Abs(h->GetMaximum()); + + const Float_t absmax = TMath::Max(min, max)*1.002; + + h->SetMaximum( absmax); + h->SetMinimum(-absmax); } @@ -337,8 +354,12 @@ // and/or MTaskList // Source dependant distance cut - if (fDistMin>0 && hsrc.GetDist()*fMm2Deg0 && hsrc.GetDist()*fMm2Deg0 && hil->GetLength()>fDistMax*hsrc.GetDist()) + if (fMaxDist>0 && hsrc.GetDist()*fMm2Deg>fMaxDist) + continue; + + if (fMaxLD>0 && hil->GetLength()>fMaxLD*hsrc.GetDist()) + continue; + if (fMinLD>0 && hil->GetLength()Reset(); - h2->Add(p); + h2->Add(p, all->GetMaximum()); + h2->Divide(all); // Delete p @@ -393,5 +415,5 @@ // range (0, fAlphaCut) // -void MHFalseSource::ProjectOn(TH2D *h3) +void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all) { TAxis &axe = *fHist.GetZaxis(); @@ -410,5 +432,8 @@ // Move contents from projection to h3 h3->Reset(); - h3->Add(p); + h3->Add(p, all->GetMaximum()); + h3->Divide(all); + + // Delete p delete p; @@ -454,22 +479,30 @@ TH2D* h3; TH2D* h4; + TH2D* h5; // Update projection of all-events - padsave->GetPad(1)->cd(3); + padsave->GetPad(2)->cd(3); if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all"))) ProjectAll(h0); + // Update projection of on-events + padsave->GetPad(1)->cd(1); + if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) + ProjectOn(h3, h0); + // Update projection of off-events - padsave->GetPad(2)->cd(1); + padsave->GetPad(1)->cd(3); if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) - ProjectOff(h2); - - // Update projection of on-events + ProjectOff(h2, h0); + padsave->GetPad(2)->cd(2); - if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) - ProjectOn(h3); + if ((h5 = (TH2D*)gPad->FindObject("Alpha_yx_diff"))) + { + h5->Add(h3, h2, -1); + MakeSymmetric(h5); + } // Update projection of significance - padsave->GetPad(2)->cd(3); + padsave->GetPad(1)->cd(2); if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig"))) { @@ -483,8 +516,8 @@ Int_t max = h4->GetBin(nx, ny); - for (int ix=0; ixGetBin(ix+1, iy+1); + const Int_t n = h4->GetBin(ix, iy); const Double_t s = h3->GetBinContent(n); @@ -503,6 +536,8 @@ } + MakeSymmetric(h4); + // Update projection of 'the best alpha-plot' - padsave->GetPad(1)->cd(1); + padsave->GetPad(2)->cd(1); if ((h1 = (TH1D*)gPad->FindObject("Alpha")) && max>0) { @@ -560,49 +595,11 @@ TObject *catalog = GetCatalog(); - // draw the 2D histogram Sigmabar versus Theta + // Initialize upper part pad->cd(1); gPad->SetBorderMode(0); gPad->Divide(3, 1); - delete pad->GetPad(1)->GetPad(2); - - pad->GetPad(1)->cd(3); - gPad->SetBorderMode(0); - TH1 *h0 = fHist.Project3D("yx_all"); - h0->SetDirectory(NULL); - h0->SetXTitle(fHist.GetXaxis()->GetTitle()); - h0->SetYTitle(fHist.GetYaxis()->GetTitle()); - h0->Draw("colz"); - h0->SetBit(kCanDelete); - catalog->Draw("mirror same"); - + + // PAD #1 pad->GetPad(1)->cd(1); - gPad->SetBorderMode(0); - - TH1 *h1 = fHist.ProjectionZ("Alpha"); - h1->SetDirectory(NULL); - h1->SetTitle("Distribution of \\alpha"); - h1->SetXTitle(fHist.GetZaxis()->GetTitle()); - h1->SetYTitle("Counts"); - h1->Draw(opt); - h1->SetBit(kCanDelete); - - pad->cd(2); - gPad->SetBorderMode(0); - gPad->Divide(3, 1); - - pad = gPad; - - pad->cd(1); - gPad->SetBorderMode(0); - fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); - TH1 *h2 = fHist.Project3D("yx_off"); - h2->SetDirectory(NULL); - h2->SetXTitle(fHist.GetXaxis()->GetTitle()); - h2->SetYTitle(fHist.GetYaxis()->GetTitle()); - h2->Draw("colz"); - h2->SetBit(kCanDelete); - catalog->Draw("mirror same"); - - pad->cd(2); gPad->SetBorderMode(0); fHist.GetZaxis()->SetRangeUser(0,fAlphaCut); @@ -616,5 +613,6 @@ catalog->Draw("mirror same"); - pad->cd(3); + // PAD #2 + pad->GetPad(1)->cd(2); gPad->SetBorderMode(0); fHist.GetZaxis()->SetRange(0,0); @@ -627,4 +625,54 @@ h4->Draw("colz"); h4->SetBit(kCanDelete); + catalog->Draw("mirror same"); + + // PAD #3 + pad->GetPad(1)->cd(3); + gPad->SetBorderMode(0); + fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); + TH1 *h2 = fHist.Project3D("yx_off"); + h2->SetDirectory(NULL); + h2->SetXTitle(fHist.GetXaxis()->GetTitle()); + h2->SetYTitle(fHist.GetYaxis()->GetTitle()); + h2->Draw("colz"); + h2->SetBit(kCanDelete); + catalog->Draw("mirror same"); + + // Initialize lower part + pad->cd(2); + gPad->SetBorderMode(0); + gPad->Divide(3, 1); + + // PAD #4 + pad->GetPad(2)->cd(1); + gPad->SetBorderMode(0); + TH1 *h1 = fHist.ProjectionZ("Alpha"); + h1->SetDirectory(NULL); + h1->SetTitle("Distribution of \\alpha"); + h1->SetXTitle(fHist.GetZaxis()->GetTitle()); + h1->SetYTitle("Counts"); + h1->Draw(opt); + h1->SetBit(kCanDelete); + + // PAD #5 + pad->GetPad(2)->cd(2); + gPad->SetBorderMode(0); + TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff"); + h5->Add(h2, -1); + h5->SetTitle("Difference of on- and off-distribution"); + h5->SetDirectory(NULL); + h5->Draw("colz"); + h5->SetBit(kCanDelete); + catalog->Draw("mirror same"); + + // PAD #6 + pad->GetPad(2)->cd(3); + gPad->SetBorderMode(0); + TH1 *h0 = fHist.Project3D("yx_all"); + h0->SetDirectory(NULL); + h0->SetXTitle(fHist.GetXaxis()->GetTitle()); + h0->SetYTitle(fHist.GetYaxis()->GetTitle()); + h0->Draw("colz"); + h0->SetBit(kCanDelete); catalog->Draw("mirror same"); } @@ -729,6 +777,6 @@ Form("Fit Region: Signal<%.1f\\circ, %.1f\\circSetNameTitle("Excess", Form("Number of excess events for \\alpha<%.0f\\circ", sigint)); - histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint)); + hists->SetName("Excess"); + histb->SetName("Background"); hist->SetXTitle(fHist.GetXaxis()->GetTitle()); hists->SetXTitle(fHist.GetXaxis()->GetTitle()); @@ -776,4 +824,8 @@ *fLog << "Polynom order: " << (int)polynom << endl; *fLog << "Fitting False Source Plot..." << flush; + + TH1 *h0 = fHist.Project3D("yx_entries"); + Float_t entries = h0->GetMaximum(); + delete h0; TH1 *h=0; @@ -790,4 +842,6 @@ if (alpha0==0) continue; + + h->Scale(entries/h->GetEntries()); if (alpha0>maxalpha0) @@ -887,6 +941,12 @@ h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5); + hists->SetTitle(Form("Excess events for \\alpha<%.0f\\circ (N_{max}=%d)", sigint, (int)hists->GetMaximum())); + histb->SetTitle(Form("Background events for \\alpha<%.0f\\circ", sigint)); + //hists->SetMinimum(GetMinimumGT(*hists)); histb->SetMinimum(GetMinimumGT(*histb)); + + MakeSymmetric(hists); + MakeSymmetric(hist); clk.Stop(); @@ -950,4 +1010,6 @@ TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); + result->Scale(entries/h->GetEntries()); + result->SetDirectory(NULL); result->SetNameTitle("Result \\alpha", title); Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3667) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3682) @@ -29,6 +29,9 @@ Float_t fBgMean; // Background mean - Float_t fDistMin; // Min dist - Float_t fDistMax; // Maximum distance in percent of dist + Float_t fMinDist; // Min dist + Float_t fMaxDist; // Max dist + + Float_t fMinLD; // Minimum distance in percent of dist + Float_t fMaxLD; // Maximum distance in percent of dist TH3D fHist; // Alpha vs. x and y @@ -40,9 +43,11 @@ void Modified(); - void ProjectOff(TH2D *h); - void ProjectOn(TH2D *h); void ProjectAll(TH2D *h); + void ProjectOff(TH2D *h, TH2D *all); + void ProjectOn(TH2D *h, TH2D *all); TObject *GetCatalog(); + + void MakeSymmetric(TH1 *h); public: @@ -57,6 +62,8 @@ void FitSignificanceStd() { FitSignificance(); } //*MENU* - void SetDistMin(Float_t dist) { fDistMin = dist; } // Absolute minimum distance - void SetDistMax(Float_t ratio) { fDistMax = ratio; } // Maximum ratio between length/dist + void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance + void SetMaxDist(Float_t dist) { fMaxDist = dist; } // Absolute maximum distance + void SetMinLD(Float_t ratio) { fMinLD = ratio; } // Minimum ratio between length/dist + void SetMaxLD(Float_t ratio) { fMaxLD = ratio; } // Maximum ratio between length/dist void SetAlphaCut(Float_t alpha); //*MENU*