Changeset 3666 for trunk/MagicSoft/Mars/mhist

Timestamp:

04/06/04 13:41:56 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHFalseSource.cc (modified) (33 diffs)
• MHFalseSource.h (modified) (4 diffs)
• MHStarMap.cc (modified) (1 diff)

trunk/MagicSoft/Mars/mhist/MHFalseSource.cc

@@ -99 +99 @@
 //    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 +119 @@
 #include "MObservatory.h"
 #include "MPointingPos.h"
+#include "MAstroCatalog.h"
 #include "MAstroSky2Local.h"
 #include "MStatusDisplay.h"
@@ -136 +139 @@
 //
 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 +190 @@
 // --------------------------------------------------------------------------
 //
-// 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 +198 @@
 Double_t MHFalseSource::Significance(Double_t s, Double_t b)
 {
-    return MMath::Significance(s, b);
+    return MMath::SignificanceSym(s, b);
 }
 
@@ -263 +212 @@
 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 +229 @@
 {
     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 +243 @@
     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 +262 @@
     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 +279 @@
 Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w)
 {
-    MHillas *hil = (MHillas*)par;
-
+    const MHillas *hil = dynamic_cast<const MHillas*>(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 +296 @@
     //    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 +313 @@
         for (int iy=0; iy<ny; iy++)
         {
+            // check distance... to get a circle plot
             if (TMath::Hypot(cx[ix], cy[iy])>maxr)
                 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 +334 @@
             }
 
+            // FIXME: This should be replaced by an external MFilter
+            //        and/or MTaskList
+            // Source dependant distance cut
+            if (fDistMin>0 && hsrc.GetDist()*fMm2Deg<fDistMin)
+                continue;
+
+            if (fDistMax>0 && 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 +419 @@
 // --------------------------------------------------------------------------
 //
-// 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 +450 @@
 
     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 +495 @@
                 h4->SetBinContent(n, sig);
 
-                if (sig>h4->GetBinContent(max) && sig>0 && ix*ix+iy*iy<nr*nr/9)
+                if (sig>h4->GetBinContent(max) && sig>0/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)<nr*nr/9*/)
                 {
                     max = n;
@@ -501 +503 @@
             }
 
-        padsave->cd(1);
+        // Update projection of 'the best alpha-plot'
+        padsave->GetPad(1)->cd(1);
         if ((h1 = (TH1D*)gPad->FindObject("Alpha")) && max>0)
         {
@@ -520 +523 @@
 // --------------------------------------------------------------------------
 //
+// 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 +556 @@
     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 +587 @@
     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 +602 @@
     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 +614 @@
     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 +627 @@
     h4->Draw("colz");
     h4->SetBit(kCanDelete);
+    catalog->Draw("mirror same");
 }
 
@@ -596 +650 @@
     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 +695 @@
 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 +744 @@
     // 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 +759 @@
     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 +803 @@
 
             h->Fit(&func, "N0Q", "", bgmin, bgmax);
-            //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
 
             h4a.Fill(func.GetChisquare());
@@ -774 +828 @@
 
             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 +873 @@
                 h6.Fill(sig);
 
-            if (sig>maxs && ix*ix+iy*iy<nr*nr/9)
+            if (sig>maxs/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)<nr*nr/9*/)
             {
                 maxs = sig;
@@ -838 +891 @@
 
     clk.Stop();
-    clk.Print();
+    clk.Print("m");
 
     TCanvas *c=new TCanvas;
@@ -849 +902 @@
     hists->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 +998 @@
             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);

trunk/MagicSoft/Mars/mhist/MHFalseSource.h

@@ -12 +12 @@
 class TH2D;
 
-class MHillasSrc;
-class MEnergyEst;
 class MParList;
 class MTime;
@@ -22 +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 +42 @@
     void ProjectOff(TH2D *h);
     void ProjectOn(TH2D *h);
+    void ProjectAll(TH2D *h);
+
+    TObject *GetCatalog();
 
 public:
@@ -46 +52 @@
     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*

trunk/MagicSoft/Mars/mhist/MHStarMap.cc

@@ -34 +34 @@
 // 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"