Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3597) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3666) @@ -99,4 +99,6 @@ // also different algorithms like (Li/Ma) // - implement fit for best alpha distribution -- online (Paint) +// - currently a constant pointing position is assumed in Fill() +// - the center of rotation need not to be the camera center // ////////////////////////////////////////////////////////////////////////////// @@ -117,4 +119,5 @@ #include "MObservatory.h" #include "MPointingPos.h" +#include "MAstroCatalog.h" #include "MAstroSky2Local.h" #include "MStatusDisplay.h" @@ -136,5 +139,6 @@ // MHFalseSource::MHFalseSource(const char *name, const char *title) - : fMm2Deg(-1), fUseMmScale(kTRUE), fAlphaCut(12.5), fBgMean(55) + : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), + fAlphaCut(12.5), fBgMean(55), fDistMin(-1), fDistMax(-1) { // @@ -186,59 +190,4 @@ // -------------------------------------------------------------------------- // -// Use this function to setup your own conversion factor between degrees -// and millimeters. The conversion factor should be the one calculated in -// MGeomCam. Use this function with Caution: You could create wrong values -// by setting up your own scale factor. -// -void MHFalseSource::SetMm2Deg(Float_t mmdeg) -{ - if (mmdeg<0) - { - *fLog << warn << dbginf << "Warning - Conversion factor < 0 - nonsense. Ignored." << endl; - return; - } - - if (fMm2Deg>=0) - *fLog << warn << dbginf << "Warning - Conversion factor already set. Overwriting" << endl; - - fMm2Deg = mmdeg; -} - -// -------------------------------------------------------------------------- -// -// With this function you can convert the histogram ('on the fly') between -// degrees and millimeters. -// -void MHFalseSource::SetMmScale(Bool_t mmscale) -{ - if (fUseMmScale == mmscale) - return; - - if (fMm2Deg<0) - { - *fLog << warn << dbginf << "Warning - Sorry, no conversion factor for conversion available." << endl; - return; - } - - if (fUseMmScale) - { - fHist.SetXTitle("x [mm]"); - fHist.SetYTitle("y [mm]"); - - fHist.Scale(1./fMm2Deg); - } - else - { - fHist.SetXTitle("x [\\circ]"); - fHist.SetYTitle("y [\\circ]"); - - fHist.Scale(1./fMm2Deg); - } - - fUseMmScale = mmscale; -} - -// -------------------------------------------------------------------------- -// // Calculate Significance as // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b @@ -249,5 +198,5 @@ Double_t MHFalseSource::Significance(Double_t s, Double_t b) { - return MMath::Significance(s, b); + return MMath::SignificanceSym(s, b); } @@ -263,6 +212,9 @@ Double_t MHFalseSource::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha) { - const Double_t lima = MMath::SignificanceLiMa(s, b); - return lima<0 ? 0 : lima; + return MMath::SignificanceLiMaSigned(s, b); + /* + const Double_t lima = MMath::SignificanceLiMa(s, b); + return lima<0 ? 0 : lima; + */ } @@ -277,12 +229,11 @@ { const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam"); - if (geom) + if (!geom) { - fMm2Deg = geom->GetConvMm2Deg(); - fUseMmScale = kFALSE; - - fHist.SetXTitle("x [\\circ]"); - fHist.SetYTitle("y [\\circ]"); + *fLog << err << "MGeomCam not found... aborting." << endl; + return kFALSE; } + + fMm2Deg = geom->GetConvMm2Deg(); MBinning binsa; @@ -292,8 +243,5 @@ if (!bins) { - Float_t r = geom ? geom->GetMaxRadius() : 600; - r /= 3; - if (!fUseMmScale) - r *= fMm2Deg; + const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg; MBinning b; @@ -314,6 +262,11 @@ fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory")); if (!fObservatory) - *fLog << err << "MObservatory not found... no derotation." << endl; - + *fLog << warn << "MObservatory not found... no derotation." << endl; + + // FIXME: Because the pointing position could change we must check + // for the current pointing position and add a offset in the + // Fill function! + fRa = fPointPos->GetRa(); + fDec = fPointPos->GetDec(); return kTRUE; @@ -326,6 +279,12 @@ Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w) { - MHillas *hil = (MHillas*)par; - + const MHillas *hil = dynamic_cast(par); + if (!hil) + { + *fLog << err << "MHFalseSource::Fill: No container specified!" << endl; + return kFALSE; + } + + // Get max radius... const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1)); @@ -337,8 +296,10 @@ // rho = fPointPos->RotationAngle(*fObservatory); + // Create necessary containers for calculation MSrcPosCam src; MHillasSrc hsrc; hsrc.SetSrcPos(&src); + // Get number of bins and bin-centers const Int_t nx = fHist.GetNbinsX(); const Int_t ny = fHist.GetNbinsY(); @@ -352,16 +313,19 @@ for (int iy=0; iymaxr) continue; + // rotate center of bin + // precalculation of sin/cos doesn't accelerate TVector2 v(cx[ix], cy[iy]); if (rho!=0) v=v.Rotate(rho); - if (!fUseMmScale) - v *= 1./fMm2Deg; - + // convert degrees to millimeters + v *= 1./fMm2Deg; src.SetXY(v); + // Source dependant hillas parameters if (!hsrc.Calc(hil)) { @@ -370,6 +334,15 @@ } + // FIXME: This should be replaced by an external MFilter + // and/or MTaskList + // Source dependant distance cut + if (fDistMin>0 && hsrc.GetDist()*fMm2Deg0 && hil->GetLength()>fDistMax*hsrc.GetDist()) + continue; + + // Fill histogram const Double_t alpha = hsrc.GetAlpha(); - fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w); } @@ -446,10 +419,30 @@ // -------------------------------------------------------------------------- // -// Update the projections +// Create projection for all data, taking sum of bin contents of +// range (0, 90) - corresponding to the number of entries in this slice. +// +void MHFalseSource::ProjectAll(TH2D *h3) +{ + h3->SetTitle("Number of entries"); + + // Get projection for range + TH2D *p = (TH2D*)fHist.Project3D("yx_all"); + + // Move contents from projection to h3 + h3->Reset(); + h3->Add(p); + delete p; + + // Set Minimum as minimum value Greater Than 0 + h3->SetMinimum(GetMinimumGT(*h3)); +} + +// -------------------------------------------------------------------------- +// +// Update the projections and paint them // void MHFalseSource::Paint(Option_t *opt) { - // sigma = (s-b)/sqrt(s+k*k*b) mit k=s/b - + // Set pretty color palette gStyle->SetPalette(1, 0); @@ -457,22 +450,31 @@ TH1D* h1; + TH2D* h0; TH2D* h2; TH2D* h3; TH2D* h4; - padsave->cd(3); + // Update projection of all-events + padsave->GetPad(1)->cd(3); + if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all"))) + ProjectAll(h0); + + // Update projection of off-events + padsave->GetPad(2)->cd(1); + if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) + ProjectOff(h2); + + // Update projection of on-events + padsave->GetPad(2)->cd(2); if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) ProjectOn(h3); - padsave->cd(4); - if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) - ProjectOff(h2); - - padsave->cd(2); + // Update projection of significance + padsave->GetPad(2)->cd(3); if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig"))) { const Int_t nx = h4->GetXaxis()->GetNbins(); const Int_t ny = h4->GetYaxis()->GetNbins(); - const Int_t nr = nx*nx + ny*ny; + //const Int_t nr = nx*nx + ny*ny; Int_t maxx=nx/2; @@ -493,5 +495,5 @@ h4->SetBinContent(n, sig); - if (sig>h4->GetBinContent(max) && sig>0 && ix*ix+iy*iyh4->GetBinContent(max) && sig>0/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)cd(1); + // Update projection of 'the best alpha-plot' + padsave->GetPad(1)->cd(1); if ((h1 = (TH1D*)gPad->FindObject("Alpha")) && max>0) { @@ -520,4 +523,28 @@ // -------------------------------------------------------------------------- // +// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude +// is set to 9, while the fov is equal to the current fov of the false +// source plot. +// +TObject *MHFalseSource::GetCatalog() +{ + const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1)); + + // Create catalog... + MAstroCatalog stars; + stars.SetLimMag(9); + stars.SetGuiActive(kFALSE); + stars.SetRadiusFOV(maxr); + stars.SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad()); + stars.ReadBSC("bsc5.dat"); + + TObject *o = (MAstroCatalog*)stars.Clone(); + o->SetBit(kCanDelete); + + return o; +} + +// -------------------------------------------------------------------------- +// // Draw the histogram // @@ -529,9 +556,27 @@ AppendPad(""); - pad->Divide(2, 2); + pad->Divide(1, 2, 0, 0.03); + + TObject *catalog = GetCatalog(); // draw the 2D histogram Sigmabar versus Theta 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->GetPad(1)->cd(1); + gPad->SetBorderMode(0); + TH1 *h1 = fHist.ProjectionZ("Alpha"); h1->SetDirectory(NULL); @@ -542,5 +587,11 @@ h1->SetBit(kCanDelete); - pad->cd(4); + 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); @@ -551,6 +602,7 @@ h2->Draw("colz"); h2->SetBit(kCanDelete); - - pad->cd(3); + catalog->Draw("mirror same"); + + pad->cd(2); gPad->SetBorderMode(0); fHist.GetZaxis()->SetRangeUser(0,fAlphaCut); @@ -562,6 +614,7 @@ h3->Draw("colz"); h3->SetBit(kCanDelete); - - pad->cd(2); + catalog->Draw("mirror same"); + + pad->cd(3); gPad->SetBorderMode(0); fHist.GetZaxis()->SetRange(0,0); @@ -574,4 +627,5 @@ h4->Draw("colz"); h4->SetBit(kCanDelete); + catalog->Draw("mirror same"); } @@ -596,11 +650,13 @@ TVirtualPad *padsave = gPad; padsave->Modified(); - padsave->cd(1); + padsave->GetPad(1)->cd(1); gPad->Modified(); - padsave->cd(2); + padsave->GetPad(1)->cd(3); gPad->Modified(); - padsave->cd(3); + padsave->GetPad(2)->cd(1); gPad->Modified(); - padsave->cd(4); + padsave->GetPad(2)->cd(2); + gPad->Modified(); + padsave->GetPad(2)->cd(3); gPad->Modified(); gPad->cd(); @@ -639,4 +695,6 @@ void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom) { + TObject *catalog = GetCatalog(); + TH1D h0a("A", "", 50, 0, 4000); TH1D h4a("chisq1", "", 50, 0, 35); @@ -686,13 +744,10 @@ // Implementing the function yourself is only about 5% faster TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90); + //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90); TArrayD maxpar(func.GetNpar()); - /* - func.SetParName(0, "A"); - func.SetParName(1, "mu"); - func.SetParName(2, "sigma"); - func.SetParName(3, "a"); - func.SetParName(4, "b"); - func.SetParName(5, "c"); + /* func.SetParName(0, "A"); + * func.SetParName(1, "mu"); + * func.SetParName(2, "sigma"); */ @@ -704,5 +759,5 @@ const Int_t nx = hist->GetXaxis()->GetNbins(); const Int_t ny = hist->GetYaxis()->GetNbins(); - const Int_t nr = nx*nx+ny*ny; + //const Int_t nr = nx*nx+ny*ny; Double_t maxalpha0=0; @@ -748,5 +803,4 @@ h->Fit(&func, "N0Q", "", bgmin, bgmax); - //*fLog << dbg << ix << "/" << iy << ": " << func.GetParameter(3) << " " << func.GetParameter(4) << endl; h4a.Fill(func.GetChisquare()); @@ -774,5 +828,4 @@ h->Fit(&func, "N0Q", "", 0, sigmax); - //*fLog << dbg << " " << func.GetParameter(0) << " " << func.GetParameter(1) << " " << func.GetParameter(2) << endl; TArrayD p(func.GetNpar(), func.GetParameters()); @@ -820,5 +873,5 @@ h6.Fill(sig); - if (sig>maxs && ix*ix+iy*iymaxs/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)Draw("colz"); hists->SetBit(kCanDelete); + catalog->Draw("mirror same"); c->cd(2); gPad->SetBorderMode(0); hist->Draw("colz"); hist->SetBit(kCanDelete); + catalog->Draw("mirror same"); c->cd(3); gPad->SetBorderMode(0); histb->Draw("colz"); histb->SetBit(kCanDelete); + catalog->Draw("mirror same"); c->cd(4); gPad->Divide(1,3, 0, 0); @@ -942,5 +998,5 @@ const Double_t b = f2.Integral(0, (float)i)/w; - const Double_t sig = Significance(s, b); + const Double_t sig = SignificanceLiMa(s, b); g->SetPoint(g->GetN(), i, sig); Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3597) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3666) @@ -12,6 +12,4 @@ class TH2D; -class MHillasSrc; -class MEnergyEst; class MParList; class MTime; @@ -22,15 +20,20 @@ { private: - MTime *fTime; //! container to take the event time from - MPointingPos *fPointPos; //! container to take pointing position from - MObservatory *fObservatory; //! conteiner to take observatory location from + MTime *fTime; //! container to take the event time from + MPointingPos *fPointPos; //! container to take pointing position from + MObservatory *fObservatory; //! conteiner to take observatory location from - Float_t fMm2Deg; // conversion factor for display in degrees - Bool_t fUseMmScale; // which scale to use? + Float_t fMm2Deg; // conversion factor for display in degrees - Float_t fAlphaCut; // Alpha cut - Float_t fBgMean; // Background mean + Float_t fAlphaCut; // Alpha cut + Float_t fBgMean; // Background mean - TH3D fHist; // Alpha vs. x and y + Float_t fDistMin; // Min dist + Float_t fDistMax; // Maximum distance in percent of dist + + TH3D fHist; // Alpha vs. x and y + + Double_t fRa; + Double_t fDec; Int_t DistancetoPrimitive(Int_t px, Int_t py); @@ -39,4 +42,7 @@ void ProjectOff(TH2D *h); void ProjectOn(TH2D *h); + void ProjectAll(TH2D *h); + + TObject *GetCatalog(); public: @@ -46,11 +52,11 @@ Bool_t Fill(const MParContainer *par, const Stat_t w=1); - void SetMmScale(Bool_t mmscale=kTRUE); - void SetMm2Deg(Float_t mmdeg); - TH1 *GetHistByName(const TString name) { return &fHist; } void FitSignificance(Float_t sigint=15, Float_t sigmax=70, Float_t bgmin=40, Float_t bgmax=70, Byte_t polynom=1); //*MENU* 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 SetAlphaCut(Float_t alpha); //*MENU* Index: trunk/MagicSoft/Mars/mhist/MHStarMap.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHStarMap.cc (revision 3597) +++ trunk/MagicSoft/Mars/mhist/MHStarMap.cc (revision 3666) @@ -34,11 +34,15 @@ // For a given a shower, a series of points along its main axis are filled // into the 2-dim histogram (x, y) of the camera plane. -// +// // Before filling a point (x, y) into the histogram it is // - shifted by (xSrc, ySrc) (the expected source position) // - and rotated in order to compensate the rotation of the // sky image in the camera -// -// +// +// To be done: +// - simplification of the algorithm, by doing all translation before +// calculation of the line +// - the center of rotation need not to be the camera center +// ///////////////////////////////////////////////////////////////////////////// #include "MHStarMap.h"