Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc	(revision 4830)
+++ /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc	(revision 4830)
@@ -0,0 +1,1216 @@
+/* ======================================================================== *\
+!
+! *
+! * 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>
+!   Author(s): Daniel Mazin, 8/2004 <mailto:mazin@mppmu.mpg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// MSkyPlot
+//
+// 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.
+//
+// MSkyPlot 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'
+//
+// 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 border of the camera. This, however, produces more noise
+// at the border.
+//
+// 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=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("MSkyPlot", "MHillas");
+// or
+//    MSkyPlot 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
+//  - more drawing options...
+//  - Move Significance() to a more 'general' place and implement
+//    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
+//  - ERRORS on each alpha bin to estimate the chi^2 correctly!
+//    (sqrt(N)/binwidth) needed for WOlfgangs proposed caluclation
+//    of alpha(Li/Ma)
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MSkyPlot.h"
+
+#include <TF1.h>
+#include <TH2.h>
+#include <TH1.h>
+#include <TGraph.h>
+#include <TStyle.h>
+#include <TLatex.h>
+#include <TCanvas.h>
+#include <TPaveText.h>
+#include <TStopwatch.h>
+#include <TFile.h>
+
+#include "MGeomCam.h"
+#include "MSrcPosCam.h"
+#include "MReportDrive.h"
+#include "MHillasSrc.h"
+#include "MHillas.h"
+#include "MHillasExt.h"
+#include "MNewImagePar.h"
+#include "MTime.h"
+#include "MObservatory.h"
+#include "MPointingPos.h"
+#include "MAstroCatalog.h"
+#include "MAstroSky2Local.h"
+#include "MStarCamTrans.h"
+#include "MStatusDisplay.h"
+#include "MMath.h"
+#include "MSupercuts.h"
+
+#include "MBinning.h"
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MSkyPlot);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default Constructor
+//
+MSkyPlot::MSkyPlot(const char *name, const char *title)
+    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1) 
+//    fAlphaCut(12.5), BgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
+{
+
+*fLog << warn << "entering default constructor in MSkyPlot" << endl; 
+    //
+    //   set the name and title of this object
+    //
+    fName  = name  ? name  : "MSkyPlot";
+    fTitle = title ? title : "sky plot vs x, y";
+
+    fSetCenter=kTRUE;
+
+    fHistSignif.SetDirectory(NULL);
+    fHistNexcess.SetDirectory(NULL);
+    fHistOn.SetDirectory(NULL);
+
+    fHistSignif.SetName("SkyPlotSignif");
+    fHistSignif.SetTitle("Sky Plot of significance vs x, y");
+    fHistSignif.SetXTitle("x [\\circ]");
+    fHistSignif.SetYTitle("y [\\circ]");
+
+    fHistNexcess.SetName("SkyPlotNexcess");
+    fHistNexcess.SetTitle("Sky Plot of number of excess vs x, y");
+    fHistNexcess.SetXTitle("x [\\circ]");
+    fHistNexcess.SetYTitle("y [\\circ]");
+
+    fHistOn.SetName("SkyPlotOn");
+    fHistOn.SetTitle("Sky Plot of number of On events vs x, y");
+    fHistOn.SetXTitle("x [\\circ]");
+    fHistOn.SetYTitle("y [\\circ]");
+
+    // set some values for the sky plot geometry:
+    fMinXGrid    =-1.;  		// [deg] 
+    fMaxXGrid    = 1.;             //  [deg] , right edge of the skyplot
+    fMinYGrid    =-1.;            //  [deg] , upper edge of the skyplot
+    fMaxYGrid    = 1.;             //  [deg] , lower edge of the skyplot
+    fBinStepGrid = 0.05;         //  [deg], 
+    fAlphaONMax  = 5.;           //  [deg] , upper cut for alpha ON region in the alpha plot
+                      // [deg], ON region in the alpha plot, maybe 5 deg is better
+                      // NOTE: up to now only values of 5, 10, 15, 20 degrees are possible
+// ra,dec lines from wolfgang:
+    fGridBinning = 0.50;   // degrees
+    fGridFineBin = 0.01;   // degrees
+
+
+    // some filter cuts:
+    fSizeMin    = 2000.;             // min size in photons
+    fSizeMax    = 100000.;             // max size in photons
+    fLeakMax    = 0.25;             // leakmax  in per cent
+    fMaxDist = 1.4;          // dist max cut (ever possible)
+    fMinDist = 0.1;          // dist max cut (ever possible)
+
+    fNumBinsAlpha   = 36; 	 // number of bins for alpha histograms
+    fAlphaLeftEdge  = 0.;        // [deg] left edge
+    fAlphaRightEdge = 90.;       // [deg] left edge
+
+    fAlphaBgLow     = 30.;
+    fAlphaBgUp      = 90.;
+
+    {
+// SET DEFAULT VALUES HERE	
+    fLengthUp[0] =  0.2;
+    fLengthUp[1] =  0.0;
+    fLengthUp[2] =  0.0;
+    fLengthUp[3] =  0.0;
+    fLengthUp[4] =  0.0;
+    fLengthUp[5] =  0.0;
+    fLengthUp[6] =  0.0;
+    fLengthUp[7] =  0.0;
+
+    fLengthLo[0] =  0.;
+    fLengthLo[1] =  0.;
+    fLengthLo[2] =  0.;
+    fLengthLo[3] =  0.;
+    fLengthLo[4] =  0.;
+    fLengthLo[5] =  0.;
+    fLengthLo[6] =  0.;
+    fLengthLo[7] =  0.;
+
+    fWidthUp[0] =  0.1;
+    fWidthUp[1] =  0.0;
+    fWidthUp[2] =  0.0;
+    fWidthUp[3] =  0.0;
+    fWidthUp[4] =  0.0;
+    fWidthUp[5] =  0.0;
+    fWidthUp[6] =  0.0;
+    fWidthUp[7] =  0.0;
+
+    fWidthLo[0] =  0.;
+    fWidthLo[1] =  0.;
+    fWidthLo[2] =  0.;
+    fWidthLo[3] =  0.;
+    fWidthLo[4] =  0.;
+    fWidthLo[5] =  0.;
+    fWidthLo[6] =  0.;
+    fWidthLo[7] =  0.;
+
+    fDistUp[0] =  1.e10;
+    fDistUp[1] =  0.0;
+    fDistUp[2] =  0.0;
+    fDistUp[3] =  0.0;
+    fDistUp[4] =  0.0;
+    fDistUp[5] =  0.0;
+    fDistUp[6] =  0.0;
+    fDistUp[7] =  0.0;
+
+    fDistLo[0] =  0.0;
+    fDistLo[1] =  0.0;
+    fDistLo[2] =  0.0;
+    fDistLo[3] =  0.0;
+    fDistLo[4] =  0.0;
+    fDistLo[5] =  0.0;
+    fDistLo[6] =  0.0;
+    fDistLo[7] =  0.0;
+    }
+}
+
+void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root")
+{
+    //--------------------------------
+    // create container for the supercut parameters
+    // and set them to their initial values
+    MSupercuts super;
+
+  // read the cuts coefficients
+    TFile inparam(parSCinit);
+    super.Read("MSupercuts");
+    inparam.Close();
+    *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file "
+          << parSCinit << endl;
+
+    TArrayD params = super.GetParameters();
+    TArrayD steps  = super.GetStepsizes();
+ // TMP2
+    if (params.GetSize() == steps.GetSize())
+    {
+            *fLog << "SC parameters and Setps are: " << endl;
+            for (Int_t z = 0; z < params.GetSize(); z++)
+            {
+                cout << params[z] << setw(20) << steps[z] << endl;
+            }
+    }
+ // ENDTMP2
+    for (Int_t i=0; i<8; i++)
+    {
+	fLengthUp[i]  = params[i];
+        fLengthLo[i] = params[i+8];
+	fWidthUp[i]   = params[i+16];
+	fWidthLo[i]  = params[i+24];
+        fDistUp[i]    = params[i+32];
+        fDistLo[i]   = params[i+40];
+    }
+}
+
+void MSkyPlot::SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step)
+{
+    Float_t temp;
+
+    if (xmax<xmin)
+    {
+        *fLog << warn << "SetSkyPlot: xmax is smaller xmin ... exchanging them." << endl;
+	temp = xmax;
+	xmax = xmin;
+	xmin = temp;
+    }	
+
+    if (ymax<ymin)
+    {
+        *fLog << warn << "SetSkyPlot: ymax is smaller ymin ... exchanging them." << endl;
+	temp = ymax;
+	ymax = ymin;
+	ymin = temp;
+    }	
+    
+    if (step<0)
+        *fLog << warn << "SetSkyPlot: step<0... taking absolute value." << endl;
+
+    fBinStepGrid = TMath::Abs(step);
+    fMinXGrid    = xmin;
+    fMaxXGrid    = xmax;
+    fMinYGrid    = ymin;
+    fMaxYGrid    = ymax;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set the alpha cut (|alpha|<fAlphaCut) which is use to estimate the
+// number of excess events
+//
+void MSkyPlot::SetAlphaCut(Float_t alpha)
+{
+    if (alpha<0)
+        *fLog << warn << "Alpha<0... taking absolute value." << endl;
+
+    fAlphaONMax = TMath::Abs(alpha);
+}
+
+// calculate the values for the cuts:
+Double_t MSkyPlot::CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2)
+{
+
+    Double_t  limit = a[0] + a[1] * dd2 +
+                      ls  * (a[3] + a[4] * dd2) +
+                      ls2 * (a[6] + a[7] * dd2);
+
+     return limit;
+}
+
+
+
+// --------------------------------------------------------------------------
+//
+// Set binnings (takes BinningFalseSource) and prepare filling of the
+// histogram.
+//
+// Also search for MTime, MObservatory and MPointingPos
+// 
+Int_t MSkyPlot::PreProcess(MParList *plist)
+{
+*fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl; 
+    fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
+    if (!fGeomCam)
+    {
+        *fLog << err << "MGeomCam not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fMm2Deg = fGeomCam->GetConvMm2Deg();
+
+    fRepDrive = (MReportDrive*)plist->FindObject(AddSerialNumber("MReportDrive"));
+    if (!fRepDrive)
+        *fLog << warn << "MReportDrive not found... no sky map." << endl;
+
+
+    fSrcPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
+    if (!fSrcPosCam)
+        *fLog << warn << "MSrcPosCam not found... no sky map." << endl;
+
+    fPntPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MPntPosCam"));
+    if (!fPntPosCam)
+        *fLog << warn << "MPntPosCam not found... no sky map." << endl;
+
+    fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
+    if (!fPointPos)
+        *fLog << warn << "MPointingPos not found... no sky map." << endl;
+
+    fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
+    if (!fTime)
+        *fLog << warn << "MTime not found... no sky map." << endl;
+
+    fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
+    if (!fObservatory)
+        *fLog << warn << "MObservatory not found... no sky map." << endl;
+
+
+    fHillas = (MHillas*)plist->FindObject(AddSerialNumber("MHillas"));
+    if (!fHillas)
+        *fLog << err << "MHillas not found... no sky map." << endl;
+
+    fHillasExt = (MHillasExt*)plist->FindObject(AddSerialNumber("MHillasExt"));
+    if (!fHillasExt)
+        *fLog << err << "MHillasExt not found... no sky map." << endl;
+
+    fHillasSrc = (MHillasSrc*)plist->FindObject(AddSerialNumber("MHillasSrc"));
+    if (!fHillasSrc)
+        *fLog << err << "MHillasSrc not found... no sky map." << endl;
+
+    fNewImagePar = (MNewImagePar*)plist->FindObject(AddSerialNumber("MNewImagePar"));
+    if (!fNewImagePar)
+        *fLog << err << "MNewImagePar not found... no sky map." << endl;
+
+
+
+    // FIXME: Because the pointing position could change we must check
+    // for the current pointing position and add a offset in the
+    // Fill function!
+    kSaveAlphaPlots=kTRUE;
+    kSaveSkyPlots=kTRUE;
+
+   // prepare skyplot
+    fNumStepsX     =  (int) ((fMaxXGrid - fMinXGrid) / fBinStepGrid + 1.5);
+    fNumStepsY     =  (int) ((fMaxYGrid - fMinYGrid) / fBinStepGrid + 1.5);
+    fNumalphahist  =  (int) (fNumStepsX * fNumStepsY  + 0.5);
+
+    fHistSignif.SetName("SPSignif2ndOrder");
+    fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)");
+    fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
+			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
+    fHistSignif.SetDirectory(NULL);
+    fHistSignif.SetFillStyle(4000);
+    fHistSignif.UseCurrentStyle();
+
+    fHistNexcess.SetName("SPNex2ndOrder");
+    fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)");
+    fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
+			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
+    fHistNexcess.SetDirectory(NULL);
+    fHistNexcess.SetFillStyle(4000);
+    fHistNexcess.UseCurrentStyle();
+
+    fHistOn.SetName("SPOnCounted");
+    fHistOn.SetTitle("Sky Plot of ON events (counted)");
+    fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
+			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
+    fHistOn.SetDirectory(NULL);
+    fHistOn.SetFillStyle(4000);
+    fHistOn.UseCurrentStyle();
+
+
+    // prepare alpha plots
+    	// temp histogram for an alpha plot
+    TH1I *temp = new TH1I("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
+    temp->SetDirectory(NULL);
+    fHistAlphaBinWidth = temp->GetBinWidth(1);
+    // create alpha histograms
+    for (Int_t i=0; i< fNumalphahist; i++) {
+    	fHistAlpha.push_back(*temp);
+	fHistAlpha[i].SetDirectory(NULL);
+    }
+
+    delete temp;
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Fill the histogram. For details see the code or the class description
+// 
+Int_t MSkyPlot::Process()
+{
+// check whether MPointingPos comtains something:
+  if (TMath::Abs(fPointPos->GetRa())<1e-3 && TMath::Abs(fPointPos->GetDec())<1e-3) 
+  {
+     *fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
+     return kTRUE;
+  }
+
+    // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
+  if (fSetCenter==kTRUE)
+  {
+     fRa0  = fRepDrive->GetRa();  // [h]
+     fDec0 = fRepDrive->GetDec(); // [deg]
+     fSetCenter=kFALSE;
+  }
+
+
+// some filter cuts:
+  if ( fHillas->GetSize() > fSizeMin && fHillas->GetSize() < fSizeMax && fNewImagePar->GetLeakage1() < fLeakMax)
+  {
+
+    Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
+    Double_t dx, dy, beta, tanbeta, alphapar, distpar;
+    Double_t logsize, lgsize, lgsize2, dist2; 
+    const Double_t log3000 = log(3000.);
+    const Float_t fDistOffset = 0.9;
+
+    // Get camera rotation angle
+    Double_t rho = 0;
+    if (fTime && fObservatory && fPointPos)
+{
+        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
+*fLog << inf << " rho = " << rho << ", Zd = " << fPointPos->GetZd() << 
+                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
+}
+    //if (fPointPos)
+    //    rho = fPointPos->RotationAngle(*fObservatory);
+/*
+//TEMP
+//        theta = mcevt->GetTelescopeTheta();
+    Double_t theta, phi, sing, cosg;
+    theta = fPointPos->GetZd()*TMath::Pi()/180.;
+    phi = fPointPos->GetAz()*TMath::Pi()/180.;
+//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
+    fObservatory->RotationAngle(theta, phi, sing, cosg);
+      
+// conclusion: diffference in rho = 7 deg
+  *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ",  "
+        << phi << ",  " << cosg << ",  " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;
+
+    Double_t a1 =  0.876627;
+    Double_t a3 = -0.481171;
+    Double_t thetaTel=theta, phiTel=phi;
+
+    Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
+                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
+                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
+    Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
+    Double_t sinal =    a1 * sin(phiTel) * denom;
+
+  *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ",  "
+        << phiTel << ",  " << cosal << ",  " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;
+
+
+
+// END TEMP
+*/
+
+    // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 
+
+// Get Rot. Angle:
+	cosgam = TMath::Cos(rho);
+	singam = TMath::Sin(rho);
+// Get x_0, y_0 for RA_0,DEC_0 of the current event
+        x_0 = fPntPosCam->GetX()*fMm2Deg; 
+        y_0 = fPntPosCam->GetY()*fMm2Deg; 
+
+    Int_t index = 0;  // index for alpha histograms
+    // loop over xtilde
+    for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)   
+    {
+	ysource = fMinYGrid + gridy * fBinStepGrid;
+        // loop over ytilde
+	for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
+	{
+	   xsource = fMinXGrid + gridx * fBinStepGrid;
+ /*     derotation    : rotate  into camera coordinates */
+ /*     formel: (x_cam)    (cos(gam)  -sin(gam))   (xtilde)   (x_0)
+                (     )  = (                   ) * (      ) + (   ) 
+	        (y_cam)    (sin(gam)   sin(gam))   (ytilde)   (y_0)
+*/
+           xsourcam = cosgam * xsource - singam * ysource + x_0;
+           ysourcam = singam * xsource + cosgam * ysource + y_0;
+ /*    end derotatiom    */
+
+       /* calculate ALPHA und DIST according to the source position */
+           dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
+           dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
+           tanbeta = dy / dx ;
+           beta = TMath::ATan(tanbeta);
+           alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
+           distpar = sqrt( dy*dy + dx*dx );
+           if(alphapar >  90.) alphapar -= 180.;
+           if(alphapar < -90.) alphapar += 180.;
+           alphapar = TMath::Abs(alphapar);
+
+// apply cuts
+           logsize = log(fHillas->GetSize());
+           lgsize = logsize-log3000;
+           lgsize2 = lgsize*lgsize;
+           dist2 = (distpar-fDistOffset)*(distpar-fDistOffset);
+
+           if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
+                (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
+           if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
+                (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
+           if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
+                      distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
+                      distpar < fMaxDist && distpar > fMinDist)
+	   {
+// gamma candidates!
+		fHistAlpha[index].Fill(alphapar);
+		index++;
+	   }
+	}
+     }
+  }
+
+    return kTRUE;
+}
+
+// calculate number of events below alphacut, number of excess events, significance
+
+Int_t MSkyPlot::PostProcess()
+{
+
+    Int_t nrow, ncolumn;
+    Double_t Non, chisquarefit, numdegfreed, Noff_fit, Nex, Sign; 
+    const Int_t onbinsalpha = TMath::Nint(fAlphaONMax/fHistAlphaBinWidth);
+
+// fit function for the background
+    TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", fAlphaBgLow, fAlphaBgUp);
+    fitbgpar->SetLineColor(2);
+
+// number degrees of freedom:
+    numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.; 
+
+    vector <TH1I>::iterator alpha_iterator;
+    int index2 = 0;  // index of the TH2F histograms
+
+    alpha_iterator = fHistAlpha.begin();
+
+    while( alpha_iterator != fHistAlpha.end() ) {
+	 // find the bin in the 2dim histogram
+         nrow    = index2/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
+         ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
+
+         // number of ON events below fAlphaONMax
+         Non = 0.;
+         for(Int_t i=1; i<=onbinsalpha;i++) Non += (*alpha_iterator).GetBinContent(i);
+
+         fHistOn.SetBinContent(ncolumn, nrow, Non);		// fill in the fHistOn
+
+         // fit parabola to a background region
+         (*alpha_iterator).Fit(fitbgpar,"NWR");  // NWR OK?????????????????????????
+         // get Chi2
+         chisquarefit = fitbgpar->GetChisquare();
+         if (chisquarefit/numdegfreed<2. && chisquarefit/numdegfreed>0.01);
+         else  *fLog << warn << "Fit is bad, chi2/ndf = " << chisquarefit/numdegfreed << endl;
+
+         // estimate Noff from the fit:
+       	 Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * pow(fAlphaONMax,3.) +
+       	            (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;
+
+	 Nex = Non - Noff_fit;
+
+         fHistNexcess.SetBinContent(ncolumn, nrow, Nex);	// fill in the fHistNexcess
+
+         if (Noff_fit<0.) Sign = 0.;
+//         else Sign = LiMa17(Non,Noff_fit,1.);
+         else Sign = MMath::SignificanceLiMaSigned(Non, Noff_fit, 1.);
+
+         fHistSignif.SetBinContent(ncolumn, nrow, Sign);	// fill in the fHistSignif
+
+         alpha_iterator++;
+         index2++;
+    }
+
+//temp
+   Double_t decl, hang;
+   MStarCamTrans mstarc(*fGeomCam, *fObservatory);
+   mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);
+
+*fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
+*fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
+//endtemp
+
+// save alpha plots:
+  if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots();
+
+// save sky plots:
+  if(kSaveSkyPlots==kTRUE) SaveSkyPlots();
+
+  fHistAlpha.clear();
+  return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// 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 *MSkyPlot::GetCatalog()
+{
+    const Double_t maxr = 0.98*TMath::Abs(fHistSignif.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;
+}
+
+
+void MSkyPlot::SaveSkyPlots(const TString skyplotfilename="skyplots.root")
+{
+   TFile rootfile(skyplotfilename, "RECREATE",
+                   "sky plots in some variations");
+     fHistSignif.Write();
+     fHistNexcess.Write();
+     fHistOn.Write();
+
+// from Wolfgang: 
+   //--------------------------------------------------------------------
+    // Dec-Hour lines
+    Double_t rho, sinrho, cosrho;
+//    if (fTime && fObservatory && fPointPos)
+        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
+    sinrho=TMath::Sin(rho);
+    cosrho=TMath::Cos(rho);
+
+    Double_t theta, phi, sing, cosg;
+    theta = fPointPos->GetZd()*TMath::Pi()/180.;
+    phi = fPointPos->GetAz()*TMath::Pi()/180.;
+//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
+    fObservatory->RotationAngle(theta, phi, sing, cosg);
+
+*fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
+*fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
+
+    Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
+    Double_t gridbinning = fGridBinning;
+    Double_t gridfinebin = fGridFineBin;
+    Int_t    numgridlines = (Int_t)(4.0/gridbinning);
+    Double_t aberr = 1.07;
+    Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
+
+    Double_t declin, hangle;  // [deg], [h]
+    MStarCamTrans mstarc(*fGeomCam, *fObservatory);
+    mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
+
+    TGraph *graph1;
+    TLatex *pix;
+
+    Char_t tit[100];
+    Double_t xtxt;
+    Double_t ytxt;
+
+    Double_t xlo = -534.0 * mmtodeg;
+    Double_t xup =  534.0 * mmtodeg;
+
+    Double_t ylo = -534.0 * mmtodeg;
+    Double_t yup =  534.0 * mmtodeg;
+
+    Double_t xx, yy;
+
+
+    // direction for the center of the camera
+    Double_t dec0 = declin;
+    Double_t h0   = hangle*360./24.;
+    //trans.LocToCel(theta0, phi0, dec0, h0);
+
+    TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
+    c1->cd(1);
+    fHistOn.Draw("colz");
+
+    //-----   lines for fixed dec   ------------------------------------
+
+    const Int_t Ndec = numgridlines;
+    Double_t dec[Ndec];
+    Double_t ddec = gridbinning;
+
+    Int_t nh = (Int_t)(4.0/gridfinebin);
+    const Int_t Nh   = nh+1;
+    Double_t h[Nh];
+    Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
+    if ( dh > 360.0/((Double_t)(Nh-1)) )
+      dh = 360.0/((Double_t)(Nh-1));
+
+
+    for (Int_t j=0; j<Ndec; j++)
+    {
+      dec[j] = dec0 + ((Double_t)j)*ddec
+                        - ((Double_t)(Ndec/2)-1.0)*ddec;
+    }
+
+    for (Int_t k=0; k<Nh; k++)
+    {
+      h[k] = h0 + ((Double_t)k)*dh - ((Double_t)(Nh/2)-1.0)*dh;
+    }
+
+    Double_t xh[Nh];
+    Double_t yh[Nh];
+
+
+
+    for (Int_t j=0; j<Ndec; j++)
+    {
+      if (fabs(dec[j]) > 90.0) continue;
+
+      for (Int_t k=0; k<Nh; k++)
+      {
+        Double_t hh0 = h0   *24.0/360.0;
+        Double_t hx = h[k]*24.0/360.0;
+        mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
+                           xx, yy);
+        xx = xx * mmtodeg * aberr;
+        yy = yy * mmtodeg * aberr;
+//        xh[k] = xx * mmtodeg * aberr;
+//        yh[k] = yy * mmtodeg * aberr;
+        xh[k] = cosg * xx + sing * yy;
+        yh[k] =-sing * xx + cosg * yy;
+//        xh[k] = cosrho * xx + sinrho * yy;
+//        yh[k] =-sinrho * xx + cosrho * yy;
+
+
+        //gLog << "dec0, h0 = " << dec0 << ",  " << h0
+        //     << " : dec, h, xh, yh = " << dec[j] << ",  "
+        //     << h[k] << ";   " << xh[k] << ",  " << yh[k] << endl;
+      }
+
+//      c1->cd(2);
+      graph1 = new TGraph(Nh,xh,yh);
+      graph1->SetLineColor(j+1);
+      graph1->SetLineStyle(1);
+      graph1->SetMarkerColor(j+1);
+      graph1->SetMarkerSize(.2);
+      graph1->SetMarkerStyle(20);
+      graph1->Draw("PL");
+      //delete graph1;
+
+      sprintf(tit,"Dec = %6.2f", dec[j]);
+      xtxt = xlo + (xup-xlo)*0.1;
+      ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
+      pix = new TLatex(xtxt, ytxt, tit);
+      pix->SetTextColor(j+1);
+      pix->SetTextSize(.03);
+      pix->Draw("");
+      //delete pix;
+
+    }
+
+    //-----   lines for fixed hour angle   ------------------------------------
+
+    Int_t ndec1 = (Int_t)(4.0/gridfinebin);
+    const Int_t Ndec1 = ndec1;
+    Double_t dec1[Ndec1];
+    Double_t ddec1 = gridfinebin;
+
+    const Int_t Nh1   = numgridlines;
+    Double_t h1[Nh1];
+    Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
+    if ( dh1 > 360.0/((Double_t)Nh1) )
+      dh1 = 360.0/((Double_t)Nh1);
+
+
+    for (Int_t j=0; j<Ndec1; j++)
+    {
+      dec1[j] = dec0 + ((Double_t)j)*ddec1
+                        - ((Double_t)(Ndec1/2)-1.0)*ddec1;
+    }
+
+    for (Int_t k=0; k<Nh1; k++)
+    {
+      h1[k] = h0 + ((Double_t)k)*dh1 - ((Double_t)(Nh1/2)-1.0)*dh1;
+    }
+
+
+    Double_t xd[Ndec1];
+    Double_t yd[Ndec1];
+
+    for (Int_t k=0; k<Nh1; k++)
+    {
+      Int_t count = 0;
+      for (Int_t j=0; j<Ndec1; j++)
+      {
+        if (fabs(dec1[j]) > 90.0) continue;
+
+        Double_t hh0 = h0   *24.0/360.0;
+        Double_t hhx = h1[k]*24.0/360.0;
+        mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
+                           xx, yy);
+//        xd[count] = xx * mmtodeg * aberr;
+//        yd[count] = yy * mmtodeg * aberr;
+
+        xx = xx * mmtodeg * aberr;
+        yy = yy * mmtodeg * aberr;
+        xd[count] = cosg * xx + sing * yy;
+        yd[count] =-sing * xx + cosg * yy;
+
+        //gLog << "dec0, h0 = " << dec0 << ",  " << h0
+        //     << " : dec1, h1, xd, yd = " << dec1[j] << ",  "
+        //     << h1[k] << ";   " << xd[count] << ",  " << yd[count] << endl;
+
+        count++;
+      }
+
+//      c1->cd(2);
+      graph1 = new TGraph(count,xd,yd);
+      graph1->SetLineColor(k+1);
+      graph1->SetLineStyle(2);
+      graph1->SetMarkerColor(k+1);
+      graph1->SetMarkerSize(.2);
+      graph1->SetMarkerStyle(20);
+      graph1->Draw("PL");
+      //delete graph1;
+
+      sprintf(tit,"H = %6.2f", h1[k]);
+      Double_t xtxt = xlo + (xup-xlo)*0.75;
+      Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
+      pix = new TLatex(xtxt, ytxt, tit);
+      pix->SetTextColor(k+1);
+      pix->SetTextSize(.03);
+      pix->Draw("");
+      //delete pix;
+
+    }
+
+//    c1->cd(2);
+    sprintf(tit,"Dec0 = %6.2f    H0 = %6.2f", dec0, h0);
+    xtxt = xlo + (xup-xlo)*0.05;
+    ytxt = ylo + (yup-ylo)*0.92;
+    pix = new TLatex(xtxt, ytxt, tit);
+    pix->SetTextColor(1);
+    pix->SetTextSize(.06);
+    pix->Draw("");
+    //delete pix;
+
+    sprintf(tit,"              [deg]             [deg]");
+    xtxt = xlo + (xup-xlo)*0.05;
+    ytxt = ylo + (yup-ylo)*0.86;
+    pix = new TLatex(xtxt, ytxt, tit);
+    pix->SetTextColor(1);
+    pix->SetTextSize(.06);
+    pix->Draw("");
+
+    c1->Write();
+
+    TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
+    c2->Divide(1,2);
+    c2->SetBorderMode(0);
+    c2->cd(1);
+    fHistSignif.Draw("colz");
+    c2->cd(2);
+    fHistNexcess.Draw("colz");
+    c2->Write();
+
+   rootfile.Close();
+   delete c1;
+   delete c2;
+   delete graph1;
+   delete pix;
+}
+
+void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root")
+{
+   TFile rootfile(alphaplotfilename, "RECREATE",
+                   "all the alpha plots");
+
+    vector <TH1I>::iterator alpha_iterator;
+    int index = 0;  // index of the TH2F histograms
+    Char_t strtitle[100];
+    Char_t strname[100];
+    Float_t xpos, ypos, signif, nex;
+    Int_t nrow, ncolumn, non;
+
+
+    alpha_iterator = fHistAlpha.begin();
+
+    while( alpha_iterator != fHistAlpha.end() ) {
+
+          nrow    = index/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
+          ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
+          xpos    = fMinXGrid + fBinStepGrid*(ncolumn-1);
+          ypos    = fMinYGrid + fBinStepGrid*(nrow-1);
+          non     = TMath::Nint(fHistOn.GetBinContent(ncolumn,nrow));
+          nex     = fHistNexcess.GetBinContent(ncolumn,nrow);
+          signif  = fHistSignif.GetBinContent(ncolumn,nrow);
+
+          sprintf(strname,"AlphaPlotX%.2fY%.2f", xpos, ypos); 
+          sprintf(strtitle,"for x= %.2f deg, y= %.2f deg: S= %.2f sigma, Nex= %.2f, Non= %d", 
+			xpos, ypos, signif, nex, non);
+           
+	  (*alpha_iterator).SetName(strname);
+	  (*alpha_iterator).SetTitle(strtitle);
+          (*alpha_iterator).Write();
+	  alpha_iterator++;
+          index++;
+    }
+
+   rootfile.Close();
+}
+
+// --------------------------------------------------------------------------
+//
+// Draw the histogram
+//
+void MSkyPlot::Draw(Option_t *opt)
+{
+/*
+    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
+    pad->SetBorderMode(0);
+
+    AppendPad("");
+
+    pad->Divide(1, 2, 0, 0.03);
+
+    TObject *catalog = GetCatalog();
+
+    // Initialize upper part
+    pad->cd(1);
+    gPad->SetBorderMode(0);
+    gPad->Divide(3, 1);
+*/
+/*
+    // PAD #1
+    pad->GetPad(1)->cd(1);
+    gPad->SetBorderMode(0);
+    fHist.GetZaxis()->SetRangeUser(0,fAlphaONMax);
+    TH1 *h3 = fHist.Project3D("yx_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);
+    catalog->Draw("mirror same");
+
+    // PAD #2
+    pad->GetPad(1)->cd(2);
+    gPad->SetBorderMode(0);
+    fHist.GetZaxis()->SetRange(0,0);
+    TH1 *h4 = fHist.Project3D("yx_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);
+    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");
+*/
+}
+
+// --------------------------------------------------------------------------
+//
+// Everything which is in the main pad belongs to this class!
+//
+Int_t MSkyPlot::DistancetoPrimitive(Int_t px, Int_t py)
+{
+    return 0;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set all sub-pads to Modified()
+//
+/*
+void MSkyPlot::Modified()
+{
+    if (!gPad)
+        return;
+
+    TVirtualPad *padsave = gPad;
+    padsave->Modified();
+    padsave->GetPad(1)->cd(1);
+    gPad->Modified();
+    padsave->GetPad(1)->cd(3);
+    gPad->Modified();
+    padsave->GetPad(2)->cd(1);
+    gPad->Modified();
+    padsave->GetPad(2)->cd(2);
+    gPad->Modified();
+    padsave->GetPad(2)->cd(3);
+    gPad->Modified();
+    gPad->cd();
+}
+*/
+// --------------------------------------------------------------------------
+//
+// This is a preliminary implementation of a alpha-fit procedure for
+// all possible source positions. It will be moved into its own
+// more powerfull class soon.
+//
+// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
+//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
+// or
+//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
+//
+// Parameter [1] is fixed to 0 while the alpha peak should be
+// symmetric around alpha=0.
+//
+// Parameter [4] is fixed to 0 because the first derivative at
+// alpha=0 should be 0, too.
+//
+// In a first step the background is fitted between bgmin and bgmax,
+// while the parameters [0]=0 and [2]=1 are fixed.
+//
+// In a second step the signal region (alpha<sigmax) is fittet using
+// the whole function with parameters [1], [3], [4] and [5] fixed.
+//
+// The number of excess and background events are calculated as
+//   s = int(0, sigint, gaus(0)+pol2(3))
+//   b = int(0, sigint,         pol2(3))
+//
+// The Significance is calculated using the Significance() member
+// function.
+//
+
+// I might need this
+/*
+    TCanvas *c=new TCanvas;
+
+    gStyle->SetPalette(1, 0);
+
+    c->Divide(3,2, 0, 0);
+    c->cd(1);
+    gPad->SetBorderMode(0);
+    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);
+    TVirtualPad *p=gPad;
+    p->SetBorderMode(0);
+    p->cd(1);
+    gPad->SetBorderMode(0);
+    h0b.DrawCopy();
+    h0a.DrawCopy("same");
+    p->cd(2);
+    gPad->SetBorderMode(0);
+    h3.DrawCopy();
+    p->cd(3);
+    gPad->SetBorderMode(0);
+    h2.DrawCopy();
+
+*/
Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h	(revision 4830)
+++ /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h	(revision 4830)
@@ -0,0 +1,135 @@
+#ifndef MARS_MSkyPlot
+#define MARS_MSkyPlot
+
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+#include <vector>
+
+#ifndef ROOT_TH2
+#include <TH2.h>
+#endif
+
+#ifndef ROOT_TH1
+#include <TH1.h>
+#endif
+
+class TH2D;
+class TH1I;
+
+class MParList;
+class MTime;
+class MPointingPos;
+class MSrcPosCam;
+class MReportDrive;
+class MObservatory;
+class MHillas;
+class MHillasExt;
+class MHillasSrc;
+class MNewImagePar;
+class MGeomCam;
+
+class MSkyPlot : public MTask
+{
+private:
+    MGeomCam      *fGeomCam;        //! container to take the event time from
+    MTime         *fTime;        //! container to take the event time from
+    MPointingPos  *fPointPos;    //! container to take pointing position from
+    MReportDrive  *fRepDrive;      
+    MSrcPosCam    *fSrcPosCam;    //! container with x and y of the source
+    MSrcPosCam    *fPntPosCam;    //! container with x and y of the position MReportDrive.GetRa, MReportDrive.GetDec
+    MObservatory  *fObservatory; //! container to take observatory location from
+    MHillas       *fHillas;      
+    MHillasExt    *fHillasExt;
+    MHillasSrc    *fHillasSrc;
+    MNewImagePar  *fNewImagePar;
+
+    Float_t        fMm2Deg;             // conversion factor for display in degrees
+    Double_t       fGridBinning;   // degrees
+    Double_t       fGridFineBin;   // degrees
+
+//    Float_t fAlphaCut;           // Alpha cut
+//    Float_t fBgMean;             // Background mean
+
+//    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
+
+    vector <TH1I> fHistAlpha;           // vector of histograms for alpha
+    Int_t fNumalphahist;		// number of histograms for alpha
+    Int_t fNumBinsAlpha;
+    Float_t fHistAlphaBinWidth;
+    Float_t fAlphaLeftEdge;
+    Float_t fAlphaRightEdge;
+    Float_t fAlphaONMax;			//  [deg] , upper cut for alpha ON region in the alpha plot 
+                      // [deg], ON region in the alpha plot, maybe 5 deg is better
+                      // NOTE: up to now only values of 5, 10, 15, 20 degrees are possible
+    Float_t fAlphaBgLow;	   // lower limit for bg region in the ON alpha plot 
+    Float_t fAlphaBgUp;		   // upper limit for bg region in the ON alpha plot
+    
+    TH2D     fHistSignif;                // sky plot of significance vs. x and y
+    TH2D     fHistNexcess;               // sky plot of number of excess events vs. x and y
+    TH2D     fHistOn;                    // sky plot of events below fAlphaONMax vs. x and y
+    Bool_t   fSetCenter;                 // used to set the center of these histograms once
+    Double_t fRa0;
+    Double_t fDec0;
+    Bool_t kSaveAlphaPlots;
+    Bool_t kSaveSkyPlots;
+
+    Float_t fMinXGrid;			//  [deg] , left edge of the skyplot
+    Float_t fMaxXGrid;			//  [deg] , right edge of the skyplot
+    Float_t fMinYGrid;			//  [deg] , upper edge of the skyplot
+    Float_t fMaxYGrid;			//  [deg] , lower edge of the skyplot
+    Float_t fBinStepGrid;			//  [deg] , grid size
+    Float_t fNumStepsX;			// number of bins in x 
+    Float_t fNumStepsY;			// number of bins in y
+
+
+    // some filter cuts:
+    Float_t fSizeMin;		// min size in photons
+    Float_t fSizeMax;		// max size in photons
+    Float_t fLeakMax;		// leakmax  in per cent
+    Float_t fMaxDist;          // dist max cut (ever possible)
+    Float_t fMinDist;          // dist min cut (ever possible)
+
+    // coefficients for the cuts:
+    Double_t fLengthUp[8], fLengthLo[8], fWidthUp[8], fWidthLo[8], fDistUp[8], fDistLo[8];
+
+    Double_t fRa;
+    Double_t fDec;
+
+    Int_t DistancetoPrimitive(Int_t px, Int_t py);
+
+    TObject *GetCatalog();
+
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
+    Int_t PostProcess();
+
+
+public:
+    MSkyPlot(const char *name=NULL, const char *title=NULL);
+
+//    TH1 *GetHistByName(const TString name) { return &fHist; }
+
+    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 ReadCuts(const TString parSCinit);
+    void SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step);
+    Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2);
+    void SaveSkyPlots(const TString skyplotfilename="skyplots.root");
+    void SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root");
+
+    void SetAlphaCut(Float_t alpha); 
+
+    void Draw(Option_t *option="");
+
+    ClassDef(MSkyPlot, 0) //2D-histogram in alpha, x and y
+};
+
+#endif
Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MpiLinkDef.h
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mmpi/MpiLinkDef.h	(revision 4830)
+++ /trunk/MagicSoft/Mars/mtemp/mmpi/MpiLinkDef.h	(revision 4830)
@@ -0,0 +1,9 @@
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class MSkyPlot+;
+
+#endif