Index: /trunk/MagicSoft/Mars/mhist/MHFalseSource.cc
===================================================================
--- /trunk/MagicSoft/Mars/mhist/MHFalseSource.cc	(revision 3548)
+++ /trunk/MagicSoft/Mars/mhist/MHFalseSource.cc	(revision 3548)
@@ -0,0 +1,710 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// MHFalseSource
+//
+// Create a false source plot. For the Binning in x,y the object MBinning
+// "BinningFalseSource" is taken from the paremeter list.
+//
+// The binning in alpha is currently fixed as 18bins from 0 to 90deg.
+//
+// If MTime, MObservatory and MPointingPos is available in the paremeter
+// list each image is derotated.
+//
+// MHFalseSource fills a 3D histogram with alpha distribution for
+// each (derotated) x and y position.
+//
+// Only a radius of 1.2deg around the camera center is taken into account.
+//
+// The displayed histogram is the projection of alpha<fAlphaCut into
+// the x,y plain and the projection of alpha>90-fAlphaCut
+//
+// By using the context menu (find it in a small region between the single
+// pads) you can change the AlphaCut 'online'
+//
+// Here is a slightly simplified version of the algorithm:
+// ------------------------------------------------------
+//    MHillas hil; // Taken as second argument in MFillH
+//
+//    MSrcPosCam src;
+//    MHillasSrc hsrc;
+//    hsrc.SetSrcPos(&src);
+//
+//    for (int ix=0; ix<nx; ix++)
+//        for (int iy=0; iy<ny; iy++)
+//        {
+//            TVector2 v(cx[ix], cy[iy]); //cx center of x-bin ix
+//            if (rho!=0)                 //cy center of y-bin iy
+//                v.Rotate(-rho);         //image rotation angle
+//
+//            src.SetXY(v);               //source position in the camera
+//
+//            if (!hsrc.Calc(hil))        //calc source dependant hillas
+//                return;
+//
+//            //fill absolute alpha into histogram
+//            const Double_t alpha = hsrc.GetAlpha();
+//            fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
+//        }
+//    }
+//
+// Use MHFalse Source like this:
+// -----------------------------
+//    MFillH fill("MHFalseSource", "MHillas");
+// or
+//    MHFalseSource hist;
+//    hist.SetAlphaCut(12.5);  // The default value
+//    hist.SetBgmean(55);      // The default value
+//    MFillH fill(&hist, "MHillas");
+//
+// To be implemented:
+// ------------------
+//  - a switch to use alpha or |alpha|
+//  - taking the binning for alpha from the parlist (binning is
+//    currently fixed)
+//  - a possibility to change the fit-function (eg using a different TF1)
+//  - a possibility to change the fit algorithm (eg which paremeters
+//    are fixed at which time)
+//  - use a different varaible than alpha
+//  - a possiblity to change the algorithm which is used to calculate
+//    alpha (or another parameter) is missing (this could be done using
+//    a tasklist somewhere...)
+//  - a more clever (and faster) algorithm to fill the histogram, eg by
+//    calculating alpha once and fill the two coils around the mean
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MHFalseSource.h"
+
+#include <TF1.h>
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+#include "MGeomCam.h"
+#include "MSrcPosCam.h"
+#include "MHillasSrc.h"
+#include "MTime.h"
+#include "MObservatory.h"
+#include "MPointingPos.h"
+#include "MAstroSky2Local.h"
+#include "MStatusDisplay.h"
+
+#include "MBinning.h"
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MHFalseSource);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default Constructor
+//
+MHFalseSource::MHFalseSource(const char *name, const char *title)
+    : fMm2Deg(-1), fUseMmScale(kTRUE), fAlphaCut(12.5), fBgMean(55)
+{
+    //
+    //   set the name and title of this object
+    //
+    fName  = name  ? name  : "MHFalseSource";
+    fTitle = title ? title : "3D-plot of Alpha vs x, y";
+
+    fHist.SetDirectory(NULL);
+
+    fHist.SetName("Alpha");
+    fHist.SetTitle("3D-plot of Alpha vs x, y");
+    fHist.SetXTitle("x [\\circ]");
+    fHist.SetYTitle("y [\\circ]");
+    fHist.SetZTitle("\\alpha [\\circ]");
+}
+
+// --------------------------------------------------------------------------
+//
+// Set the alpha cut (|alpha|<fAlphaCut) which is use to estimate the
+// number of excess events
+//
+void MHFalseSource::SetAlphaCut(Float_t alpha)
+{
+    if (alpha<0)
+        *fLog << warn << "Alpha<0... taking absolute value." << endl;
+
+    fAlphaCut = TMath::Abs(alpha);
+
+    Modified();
+}
+
+// --------------------------------------------------------------------------
+//
+// Set mean alpha around which the off sample is determined
+// (fBgMean-fAlphaCut/2<|fAlpha|<fBgMean+fAlphaCut/2) which is use
+// to estimate the number of off events
+//
+void MHFalseSource::SetBgMean(Float_t alpha)
+{
+    if (alpha<0)
+        *fLog << warn << "Alpha<0... taking absolute value." << endl;
+
+    fBgMean = TMath::Abs(alpha);
+
+    Modified();
+}
+
+// --------------------------------------------------------------------------
+//
+// 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;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set binnings (takes BinningFalseSource) and prepare filling of the
+// histogram.
+//
+// Also search for MTime, MObservatory and MPointingPos
+// 
+Bool_t MHFalseSource::SetupFill(const MParList *plist)
+{
+    const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
+    if (geom)
+    {
+        fMm2Deg = geom->GetConvMm2Deg();
+        fUseMmScale = kFALSE;
+
+        fHist.SetXTitle("x [\\circ]");
+        fHist.SetYTitle("y [\\circ]");
+    }
+
+    MBinning binsa;
+    binsa.SetEdges(18, 0, 90);
+
+    const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
+    if (!bins)
+    {
+        Float_t r = geom ? geom->GetMaxRadius() : 600;
+        r /= 3;
+        if (!fUseMmScale)
+            r *= fMm2Deg;
+
+        MBinning b;
+        b.SetEdges(40, -r, r);
+        SetBinning(&fHist, &b, &b, &binsa);
+    }
+    else
+        SetBinning(&fHist, bins, bins, &binsa);
+
+    fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
+    if (!fPointPos)
+        *fLog << warn << "MPointingPos not found... no derotation." << endl;
+
+    fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
+    if (!fTime)
+        *fLog << warn << "MTime not found... no derotation." << endl;
+
+    fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
+    if (!fObservatory)
+        *fLog << err << "MObservatory not found...  no derotation." << endl;
+
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Fill the histogram. For details see the code or the class description
+// 
+Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w)
+{
+    MHillas *hil = (MHillas*)par;
+
+    const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
+
+    // Get camera rotation angle
+    Double_t rho = 0;
+    if (fTime && fObservatory && fPointPos)
+    {
+        const Double_t ra  = fPointPos->GetRa();
+        const Double_t dec = fPointPos->GetDec();
+
+        rho = MAstroSky2Local(*fTime, *fObservatory).RotationAngle(ra, dec);
+    }
+
+    MSrcPosCam src;
+    MHillasSrc hsrc;
+    hsrc.SetSrcPos(&src);
+
+    const Int_t nx = fHist.GetNbinsX();
+    const Int_t ny = fHist.GetNbinsY();
+
+    Axis_t cx[nx];
+    Axis_t cy[ny];
+    fHist.GetXaxis()->GetCenter(cx);
+    fHist.GetYaxis()->GetCenter(cy);
+
+    for (int ix=0; ix<nx; ix++)
+    {
+        for (int iy=0; iy<ny; iy++)
+        {
+            if (TMath::Hypot(cx[ix], cy[iy])>maxr)
+                continue;
+
+            TVector2 v(cx[ix], cy[iy]);
+            if (rho!=0)
+                v.Rotate(-rho);
+
+            if (!fUseMmScale)
+                v *= 1./fMm2Deg;
+
+            src.SetXY(v);
+
+            if (!hsrc.Calc(hil))
+            {
+                *fLog << warn << "Calculation of MHillasSrc failed for x=" << cx[ix] << " y=" << cy[iy] << endl;
+                return kFALSE;
+            }
+
+            const Double_t alpha = hsrc.GetAlpha();
+            fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
+        }
+    }
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Create projection for off data, taking sum of bin contents of
+// range (fBgMean-fAlphaCut/2, fBgMean+fAlphaCut) Making sure to take
+// the same number of bins than for on-data
+//
+void MHFalseSource::ProjectOff(TH2D *h2)
+{
+    TAxis &axe = *fHist.GetZaxis();
+
+    // Find range to cut (left edge and width)
+    const Int_t f = axe.FindBin(fBgMean-fAlphaCut/2);
+    const Int_t l = axe.FindBin(fAlphaCut)+f-1;
+
+    axe.SetRange(f, l);
+    const Float_t cut1 = axe.GetBinLowEdge(f);
+    const Float_t cut2 = axe.GetBinUpEdge(l);
+    h2->SetTitle(Form("Distribution of %.1f\\circ<|\\alpha|<%.1f\\circ in x,y", cut1, cut2));
+
+    // Get projection for range
+    TH2D *p = (TH2D*)fHist.Project3D("xy_off");
+
+    // Reset range
+    axe.SetRange(0,9999);
+
+    // Move contents from projection to h2
+    h2->Reset();
+    h2->Add(p);
+
+    // Delete p
+    delete p;
+
+    // Set Minimum as minimum value Greater Than 0
+    h2->SetMinimum(GetMinimumGT(*h2));
+}
+
+// --------------------------------------------------------------------------
+//
+// Create projection for on data, taking sum of bin contents of
+// range (0, fAlphaCut)
+//
+void MHFalseSource::ProjectOn(TH2D *h3)
+{
+    TAxis &axe = *fHist.GetZaxis();
+
+    // Find range to cut
+    axe.SetRangeUser(0, fAlphaCut);
+    const Float_t cut = axe.GetBinUpEdge(axe.GetLast());
+    h3->SetTitle(Form("Distribution of |\\alpha|<%.1f\\circ in x,y", cut));
+
+    // Get projection for range
+    TH2D *p = (TH2D*)fHist.Project3D("xy_on");
+
+    // Reset range
+    axe.SetRange(0,9999);
+
+    // 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
+//
+void MHFalseSource::Paint(Option_t *opt)
+{
+    // sigma = (s-b)/sqrt(s+k*k*b) mit k=s/b
+
+    gStyle->SetPalette(1, 0);
+
+    TVirtualPad *padsave = gPad;
+
+    TH1D* h1;
+    TH2D* h2;
+    TH2D* h3;
+    TH2D* h4;
+
+    padsave->cd(3);
+    if ((h3 = (TH2D*)gPad->FindObject("Alpha_xy_on")))
+        ProjectOn(h3);
+
+    padsave->cd(4);
+    if ((h2 = (TH2D*)gPad->FindObject("Alpha_xy_off")))
+        ProjectOff(h2);
+
+    padsave->cd(2);
+    if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_xy_sig")))
+    {
+        const Int_t nx = h4->GetXaxis()->GetNbins();
+        const Int_t ny = h4->GetYaxis()->GetNbins();
+
+        Int_t maxx=nx/2;
+        Int_t maxy=ny/2;
+
+        Int_t max = h4->GetBin(maxx, maxy);
+
+        for (int ix=0; ix<nx; ix++)
+            for (int iy=0; iy<ny; iy++)
+            {
+                const Int_t n = h4->GetBin(ix+1, iy+1);
+
+                const Double_t s = h3->GetBinContent(n);
+                const Double_t b = h2->GetBinContent(n);
+
+                const Double_t k = b==0 ? 0 : s/b;
+                const Double_t f = s+k*k*b;
+
+                const Double_t sig = f==0 ? 0 : (s-b)/TMath::Sqrt(f);
+
+                //if (b>0 && s>0)
+                //    *fLog << dbg << b << " " << s << endl;
+
+                // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
+                h4->SetBinContent(n, sig);
+
+                if (sig>h4->GetBinContent(max) && sig!=0)
+                {
+                    max = n;
+                    maxx=ix;
+                    maxy=iy;
+                }
+            }
+
+        padsave->cd(1);
+        if ((h1 = (TH1D*)gPad->FindObject("Alpha")))
+        {
+            //h1->Reset();
+
+            const Double_t x = fHist.GetXaxis()->GetBinCenter(maxx);
+            const Double_t y = fHist.GetYaxis()->GetBinCenter(maxy);
+            const Double_t s = h4->GetBinContent(max);
+
+            TH1 *h = fHist.ProjectionZ("Alpha", maxx, maxx, maxy, maxy);
+            h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma=%.1f)", x, y, s));
+        }
+    }
+
+    gPad = padsave;
+}
+
+// --------------------------------------------------------------------------
+//
+// Draw the histogram
+//
+void MHFalseSource::Draw(Option_t *opt)
+{
+    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
+    pad->SetBorderMode(0);
+
+    AppendPad("");
+
+    pad->Divide(2, 2);
+
+    // draw the 2D histogram Sigmabar versus Theta
+    pad->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(4);
+    gPad->SetBorderMode(0);
+    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
+    TH1 *h2 = fHist.Project3D("xy_off");
+    h2->SetDirectory(NULL);
+    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
+    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
+    h2->Draw("colz");
+    h2->SetBit(kCanDelete);
+
+    pad->cd(3);
+    gPad->SetBorderMode(0);
+    fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
+    TH1 *h3 = fHist.Project3D("xy_on");
+    fHist.GetZaxis()->SetRange(0,9999);
+    h3->SetDirectory(NULL);
+    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
+    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
+    h3->Draw("colz");
+    h3->SetBit(kCanDelete);
+
+    pad->cd(2);
+    gPad->SetBorderMode(0);
+    fHist.GetZaxis()->SetRange(0,0);
+    TH1 *h4 = fHist.Project3D("xy_sig"); // Do this to get the correct binning....
+    fHist.GetZaxis()->SetRange(0,9999);
+    h4->SetTitle("Significance");
+    h4->SetDirectory(NULL);
+    h4->SetXTitle(fHist.GetXaxis()->GetTitle());
+    h4->SetYTitle(fHist.GetYaxis()->GetTitle());
+    h4->Draw("colz");
+    h4->SetBit(kCanDelete);
+
+}
+
+// --------------------------------------------------------------------------
+//
+// Everything which is in the main pad belongs to this class!
+//
+Int_t MHFalseSource::DistancetoPrimitive(Int_t px, Int_t py)
+{
+    return 0;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set all sub-pads to Modified()
+//
+void MHFalseSource::Modified()
+{
+    if (!gPad)
+        return;
+
+    TVirtualPad *padsave = gPad;
+    padsave->Modified();
+    padsave->cd(1);
+    gPad->Modified();
+    padsave->cd(2);
+    gPad->Modified();
+    padsave->cd(3);
+    gPad->Modified();
+    padsave->cd(4);
+    gPad->Modified();
+    gPad->cd();
+}
+
+Double_t fcn(Double_t *arg, Double_t *p)
+{
+    const Double_t x = arg[0];
+
+    const Double_t dx = (x-p[1])/p[2];
+
+    const Double_t f1 = p[0]*TMath::Exp(-0.5*dx*dx);
+    const Double_t f2 = p[3]*x*x + p[4];
+
+    return f1 + f2;
+}
+/*
+Double_t FcnIntegral(Double_t *arg, Double_t *p)
+{
+    static const Double_t sqrt2   = TMath::Sqrt(2);
+    static const Double_t sqrt2pi = TMath::Sqrt(TMath::TwoPi());
+
+    const Double_t x = arg[0];
+
+    const Double_t dx = (x-p[1])/p[2];
+
+    const Double_t f1 = p[0]*p[2]*sqrt2pi*TMath::Erf(dx/sqrt2)/2;
+    const Double_t f2 = p[3]*x*2;
+
+    return f1+f2;
+}
+*/
+/*
+void MHFalseSource::FitSignificance()
+{
+    TH1D h0("A",     "Parameter A",     50, 0, 10000);
+    TH1D h1("mu",    "Parameter mu",    50, -1, 1);
+    TH1D h2("sigma", "Parameter sigma", 50, 0, 20);
+    TH1D h3("b",     "Parameter b",     50, 0.001, 0.01);
+    h0.SetDirectory(NULL);
+    h1.SetDirectory(NULL);
+    h2.SetDirectory(NULL);
+    h3.SetDirectory(NULL);
+
+    TH1 *hist = fHist.Project3D("xy_fit");
+
+    //                      xmin, xmax, npar
+    TF1 func("MyFunc", fcn, 0,    90,   5);
+
+    TArrayD samx(50);
+    TArrayD samw(50);
+
+    // func.CalcGaussLegendreSamplingPoints(50, samx.GetArray(), samw.GetArray());
+    // func.IntegralFast(50, samx.GetArray(), samw.GetArray(), 0, 15);
+
+    func.SetParName(0, "A");
+    func.SetParName(1, "mu");
+    func.SetParName(2, "sigma");
+    func.SetParName(3, "a");
+    func.SetParName(4, "b");
+
+    const Int_t nx = fHist.GetXaxis()->GetNbins();
+    const Int_t ny = fHist.GetYaxis()->GetNbins();
+
+    TH1 *h=0;
+    for (int ix=0; ix<nx; ix++)
+        for (int iy=0; iy<ny; iy++)
+        {
+            h = fHist.ProjectionZ("AlphaFit", ix+1, ix+1, iy+1, iy+1);
+
+            // Check for the regios which is not filled...
+            if (h->GetBinContent(1)==0)
+                continue;
+
+            func.SetParLimits(3, 0, 1);
+
+            // First fit a polynom in the off region
+            func.FixParameter(0, 0);
+            func.FixParameter(1, 0);
+            func.FixParameter(2, 1);
+            func.ReleaseParameter(3);
+            func.ReleaseParameter(4);
+
+            h->Fit(&func, "N0Q", "", 35, 80);
+            *fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
+
+            // Now fit a gaus in the on region on top of the polynom
+            func.ReleaseParameter(0);
+            //func.ReleaseParameter(1);
+            func.ReleaseParameter(2);
+
+            func.SetParameter(0, h->GetBinContent(1));
+            func.SetParameter(2, 10);
+
+            func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
+            func.SetParLimits(2, 0, 90);
+
+            func.FixParameter(3, func.GetParameter(3));
+            func.FixParameter(4, func.GetParameter(4));
+
+            h->Fit(&func, "N0Q", "", 0, 20);
+            *fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
+
+            // Fill results into some histograms
+            h0.Fill(func.GetParameter(0));
+            h1.Fill(func.GetParameter(1));
+            h2.Fill(func.GetParameter(2));
+            h3.Fill(func.GetParameter(3));
+
+            *fLog << dbg << "     " << func.GetChisquare() << " " << func.GetNDF() << " " << func.GetNpx() << " " << func.GetNumberFreeParameters() << " " << func.GetNumberFitPoints() << endl;
+
+            const Int_t n = hist->GetBin(ix+1, iy+1);
+            hist->SetBinContent(n, func.GetParameter(0));
+        }
+
+    hist->SetMinimum(0);
+    hist->SetMaximum(10000);
+
+    TCanvas *c=new TCanvas;
+
+    c->Divide(2,2);
+    c->cd(1);
+    h0.DrawCopy();
+    c->cd(2);
+    hist->Draw("colz");
+    hist->SetDirectory(NULL);
+    hist->SetBit(kCanDelete);
+    c->cd(3);
+    h2.DrawCopy();
+    c->cd(4);
+    h3.DrawCopy();
+}
+*/
Index: /trunk/MagicSoft/Mars/mhist/MHFalseSource.h
===================================================================
--- /trunk/MagicSoft/Mars/mhist/MHFalseSource.h	(revision 3548)
+++ /trunk/MagicSoft/Mars/mhist/MHFalseSource.h	(revision 3548)
@@ -0,0 +1,70 @@
+#ifndef MARS_MHFalseSource
+#define MARS_MHFalseSource
+
+#ifndef MARS_MH
+#include "MH.h"
+#endif
+
+#ifndef ROOT_TH3
+#include <TH3.h>
+#endif
+
+class TH2D;
+
+class MHillasSrc;
+class MEnergyEst;
+class MParList;
+class MTime;
+class MPointingPos;
+class MObservatory;
+
+class MHFalseSource : public MH
+{
+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
+
+    Float_t fMm2Deg;            // conversion factor for display in degrees
+    Bool_t  fUseMmScale;        // which scale to use?
+
+    Float_t fAlphaCut;          // Alpha cut
+
+    Float_t fBgMean;            // Background mean
+
+    TH3D    fHist;              // Alpha vs. x and y
+
+    Int_t DistancetoPrimitive(Int_t px, Int_t py);
+    void Modified();
+
+    void ProjectOff(TH2D *h);
+    void ProjectOn(TH2D *h);
+
+public:
+    MHFalseSource(const char *name=NULL, const char *title=NULL);
+
+    Bool_t SetupFill(const MParList *pList);
+    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 SetAlphaCut(Float_t alpha); //*MENU*
+    void SetAlphaPlus5()  { SetAlphaCut(fAlphaCut+5); } //*MENU*
+    void SetAlphaMinus5() { SetAlphaCut(fAlphaCut-5); } //*MENU*
+
+    void SetBgMean(Float_t alpha); //*MENU*
+    void SetBgMeanPlus5()  { SetBgMean(fBgMean+5); } //*MENU*
+    void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU*
+
+    //void FitSignificance(); //*MENU*
+
+    void Paint(Option_t *opt="");
+    void Draw(Option_t *option="");
+
+    ClassDef(MHFalseSource, 1) //3D-histogram in alpha, x and y
+};
+
+#endif