Changeset 7173

Timestamp:

06/28/05 10:22:44 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (3 diffs)
• mfilter/MFMagicCuts.cc (modified) (2 diffs)
• mfilter/MFMagicCuts.h (modified) (1 diff)
• mhbase/MH.cc (modified) (5 diffs)
• mhbase/MH.h (modified) (1 diff)
• mhflux/MHDisp.cc (modified) (9 diffs)
• mhflux/MHDisp.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -21 +21 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2005/06/28 Thomas Bretz
+
+   * mfilter/MFMagicCuts.[h,cc]:
+     - first full implemtation... updates to come.
+
+   * mhbase/MH.[h,cc]:
+     - added new member function to set several palettes
+
+   * mhflux/MHDisp.[h,cc]:
+     - fixed z-axis dscription
+     - rotate filling off data by 180 deg.
+     - implemented subtracting off-data from on-data
+     - set different palettes
+
+
 
  2005/06/27 Thomas Bretz

trunk/MagicSoft/Mars/NEWS

@@ -12 +12 @@
 
    - general: MHCamera now displays the profiles in deg instead of mm
+
+   - general: MH::SetPalette offers a lot of new palettes
 
    - callisto: MCalibrationHiLoCam can now be printed from its context
@@ -27 +29 @@
      argument
 
+   - ganymed: The first version of MFMagicCuts have been released
+
    - ganymed: the Conc1 plot was incorrectly scaled in MHVsSize
 
@@ -32 +36 @@
      of bins for the signal region (NumBinsSignal) and the number of
      total bins (NumBInsTital) in the MHThetaSq histogram
+
+   - ganymed: optimized palettes for MHDisp
 
    - sponde: the zenith angle distribution is now weighted instead of

trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc

@@ -30 +30 @@
 #include "MFMagicCuts.h"
 
+#include <fstream>
+
+#include "MHMatrix.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MParList.h"
+
+#include "MGeomCam.h"
+#include "MHillasSrc.h"
+#include "MHillasExt.h"
+#include "MParameters.h"
+#include "MPointingPos.h"
+
 ClassImp(MFMagicCuts);
 
 using namespace std;
 
-
 // --------------------------------------------------------------------------
 //
@@ -40 +54 @@
 //
 MFMagicCuts::MFMagicCuts(const char *name, const char *title)
+    : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fThetaSq(0),
+    fDisp(0), fPointing(0), fMatrix(0), fVariables(30), fThetaCut(kNone),
+    fHadronnessCut(kAll)
 {
     fName  = name  ? name  : "MFMagicCuts";
     fTitle = title ? title : "Class to evaluate the MagicCuts";
-}
+
+    fMatrix = NULL;
+
+    AddToBranchList("MHillas.fWidth");
+    AddToBranchList("MHillas.fLength");
+    AddToBranchList("MHillas.fSize");
+    AddToBranchList("MHillasSrc.fAlpha");
+    AddToBranchList("MHillasSrcAnti.fAlpha");
+    AddToBranchList("MPointingPos.fZd");
+    AddToBranchList("MHillasExt.fMaxDist");
+    AddToBranchList("MHillasExt.fM3Long");
+
+    fVariables[0] =  1.42547;      // Xi
+    fVariables[1] =  0.233773;     // Theta^2
+    fVariables[2] =  0.2539460;    // Area[0]
+    fVariables[3] =  5.2149800;    // Area[1]
+    fVariables[4] =  0.1139130;    // Area[2]
+    fVariables[5] = -0.0889;       // M3long
+    fVariables[6] =  0.18;         // Xi(theta)
+}
+
+// --------------------------------------------------------------------------
+//
+// The theta cut value GetThetaSqCut() is propageted to the parameter list
+// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else.
+//
+Int_t MFMagicCuts::PreProcess(MParList *pList)
+{
+    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
+    if (!cam)
+    {
+        *fLog << err << "MGeomCam not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fMm2Deg = cam->GetConvMm2Deg();
+
+    fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
+    if (!fThetaSq)
+        return kFALSE;
+    fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
+    if (!fDisp)
+        return kFALSE;
+
+    // propagate Theta cut to the parameter list
+    // so that it can be used somewhere else
+    MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
+    if (!thetacut)
+        return kFALSE;
+    thetacut->SetVal(GetThetaSqCut());
+    //thetacut->SetReadyToSave();
+
+    MParameterD *m3lcut = (MParameterD*)pList->FindCreateObj("MParameterD", "M3lCut");
+    if (!m3lcut)
+        return kFALSE;
+    m3lcut->SetVal(fVariables[5]);
+
+    MParameterD *dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXi");
+    if (!dispxi)
+        return kFALSE;
+    dispxi->SetVal(fVariables[0]);
+
+    dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXiTheta");
+    if (!dispxi)
+        return kFALSE;
+    dispxi->SetVal(fVariables[6]);
+
+    Print();
+
+    if (fMatrix)
+        return kTRUE;
+
+    //-----------------------------------------------------------
+
+    fHil = (MHillas*)pList->FindObject("MHillas");
+    if (!fHil)
+    {
+        *fLog << err << "MHillas not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
+    if (!fHilExt)
+    {
+        *fLog << err << "MHillasExt not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
+    if (!fPointing)
+    {
+        *fLog << err << "MPointingPos not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
+    if (!fHilSrc)
+    {
+        *fLog << err << "MHillasSrc not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    if (fThetaCut&kOff)
+    {
+        fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
+        if (!fHilAnti)
+        {
+            *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
+            return kFALSE;
+        }
+    }
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the mapped value from the Matrix
+//
+Double_t MFMagicCuts::GetVal(Int_t i) const
+{
+    return (*fMatrix)[fMap[i]];
+}
+
+// --------------------------------------------------------------------------
+//
+// You can use this function if you want to use a MHMatrix instead of the
+// given containers. This function adds all necessary columns to the
+// given matrix. Afterward you should fill the matrix with the corresponding
+// data (eg from a file by using MHMatrix::Fill). If you now loop
+// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
+// will take the values from the matrix instead of the containers.
+//
+void MFMagicCuts::InitMapping(MHMatrix *mat)
+{
+    if (fMatrix)
+      return;
+
+    fMatrix = mat;
+
+    fMap[kESize]   = fMatrix->AddColumn("log10(MHillas.fSize)");
+    fMap[kEArea]   = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
+    fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
+    fMap[kEWdivL]  = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
+    fMap[kEZd]     = fMatrix->AddColumn("MPointingPos.GetZdRad");
+
+    fMap[kEAlpha]  = fMatrix->AddColumn("MHillasSrc.fAlpha");
+    fMap[kEDist]   = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
+    //fMap[kDCA]    = fMatrix->AddColumn("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg");
+    if (fThetaCut&kOff)
+    {
+        fMap[kEAlphaAnti]  = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
+        fMap[kEDistAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
+        fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
+        //fMap[kDCAAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDCA*MGeomCam.fConvMm2Deg");
+    }
+
+    //fMap[kLength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
+    //fMap[kWidth]  = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns theta squared based on the formula:
+//    p := c*(1-width/length)
+//    return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
+//
+// If par!=NULL p is stored in this MParameterD
+//
+Double_t MFMagicCuts::GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par) const
+{
+    // wl = width/l [1]
+    // d  = dist    [deg]
+    // a  = alpha   [deg]
+    const Double_t p = c*(1-wl);
+    if (par)
+        par->SetVal(p);
+    return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns squared cut value in theta: fVariables[1]^2
+//
+Double_t MFMagicCuts::GetThetaSqCut() const
+{
+    return fVariables[1]*fVariables[1];
+}
+
+// ---------------------------------------------------------------------------
+//
+// Evaluate dynamical magic-cuts
+//
+Int_t MFMagicCuts::Process()
+{
+    // Get some variables
+    const Double_t wdivl  = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
+    const Double_t lgsize = fMatrix ? GetVal(kESize)  : TMath::Log10(fHil->GetSize());
+    const Double_t zd     = fMatrix ? GetVal(kEZd)    : fPointing->GetZdRad();
+
+    // For simplicity
+    const Double_t *c = fVariables.GetArray();
+
+    // Value for Theta cut
+    const Double_t cut = GetThetaSqCut(); //c[1]*c[1];
+    const Double_t xi  = (c[0]+c[8]*(lgsize-c[7])*(lgsize-c[7])) - c[6]*zd*zd;
+
+    // Default if we return before the end
+    fResult = kFALSE;
+
+    // ------------------- Most efficient cut -----------------------
+    // ---------------------- Theta cut ON --------------------------
+    const Double_t dist    = fMatrix ? GetVal(kEDist)  : fHilSrc->GetDist()*fMm2Deg;
+    const Double_t alpha   = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha();
+
+    const Double_t thetasq = GetThetaSq(xi, wdivl, dist, alpha, fDisp);
+
+    fThetaSq->SetVal(thetasq);
+
+    if (fThetaCut&kOn && thetasq >= cut)
+        return kTRUE;
+
+    // --------------------- Efficient  cut -------------------------
+    // ---------------------- Area/M3l cut --------------------------
+    if (fHadronnessCut&kArea)
+    {
+        const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
+        const Double_t A    = lgsize>c[3] ? c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
+        if (area>=A)
+            return kTRUE;
+    }
+    if (fHadronnessCut&kM3Long)
+    {
+        const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
+        if (m3long<=c[5])
+            return kTRUE;
+    }
+
+    // ------------------- Least efficient cut -----------------------
+    // ---------------------- Theta cut OFF --------------------------
+    if (fThetaCut&kOff)
+    {
+        const Double_t distanti    = fMatrix ? GetVal(kEDistAnti)   : fHilAnti->GetDist()*fMm2Deg;
+        const Double_t alphaanti   = fMatrix ? GetVal(kEAlphaAnti)  : fHilAnti->GetAlpha();
+        const Double_t thetasqanti = GetThetaSq(xi, wdivl, distanti, alphaanti);
+        const Double_t m3lanti     = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
+
+        if (thetasqanti<cut && (fHadronnessCut&kM3Long || m3lanti>c[5]))
+            return kTRUE;
+    }
+
+    fResult = kTRUE;
+
+    return kTRUE;
+}
+
+void MFMagicCuts::SetVariables(const TArrayD &arr)
+{
+    fVariables = arr;
+}
+
+TString MFMagicCuts::GetParam(Int_t i) const
+{
+    if (i>=fVariables.GetSize())
+        return "";
+
+    return Form("Param[%2d]=%+f", i, fVariables[i]);
+}
+
+void MFMagicCuts::Print(Option_t *o) const
+{
+    *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;
+
+    *fLog << "Using Theta cut: ";
+    switch (fThetaCut)
+    {
+    case kNone:
+        *fLog << "none" << endl;
+        break;
+    case kOn:
+        *fLog << "on" << endl;
+        break;
+    case kOff:
+        *fLog << "off" << endl;
+        break;
+    case kWobble:
+        *fLog << "on and off (wobble)" << endl;
+        break;
+    }
+    *fLog << "Using hadronness cut: ";
+    switch (fHadronnessCut)
+    {
+    case kNoCut:
+        *fLog << "none" << endl;
+        break;
+    case kArea:
+        *fLog << "area" << endl;
+        break;
+    case kM3Long:
+        *fLog << "m3long" << endl;
+        break;
+    case kAll:
+        *fLog << "all" << endl;
+        break;
+    }
+    *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
+
+    const Int_t n = fVariables.GetSize();
+    for (int i=0; i<n; i+=3)
+    {
+        *fLog << setw(25) << GetParam(i);
+        *fLog << setw(25) << GetParam(i+1);
+        *fLog << setw(25) << GetParam(i+2);
+        *fLog << endl;
+    }
+    *fLog << setiosflags(ios::right);
+}
+
+Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
+{
+    ifstream fin(fname);
+    if (!fin)
+    {
+        *fLog << err << "Cannot open file " << fname << ": ";
+        *fLog << strerror(errno) << endl;
+        return kFALSE;
+    }
+
+    for (int i=0; i<fVariables.GetSize(); i++)
+    {
+        fin >> fVariables[i];
+        if (!fin)
+        {
+            *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
+            return kERROR;
+        }
+    }
+    return kTRUE;
+}
+
+Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
+{
+    ofstream fout(fname);
+    if (!fout)
+    {
+        *fLog << err << "Cannot open file " << fname << ": ";
+        *fLog << strerror(errno) << endl;
+        return kFALSE;
+    }
+
+    for (int i=0; i<fVariables.GetSize(); i++)
+        fout << setw(11) << fVariables[i] << endl;
+
+    return kTRUE;
+}
+
+Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    if (MFilter::ReadEnv(env, prefix, print)==kERROR)
+        return kERROR;
+
+
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "ThetaCut", print))
+    {
+        TString str = GetEnvValue(env, prefix, "ThetaCut", "");
+        str.ToLower();
+        str=str.Strip(TString::kBoth);
+
+        if (str==(TString)"nocut" || str==(TString)"none")
+            fThetaCut = kNone;
+        if (str==(TString)"on")
+            fThetaCut = kOn;
+        if (str==(TString)"off")
+            fThetaCut = kOff;
+        if (str==(TString)"wobble")
+            fThetaCut = kWobble;
+        rc = kTRUE;
+    }
+    if (IsEnvDefined(env, prefix, "HadronnessCut", print))
+    {
+        TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
+        str.ToLower();
+        str=str.Strip(TString::kBoth);
+
+        if (str==(TString)"nocut" || str==(TString)"none")
+            fHadronnessCut = kNoCut;
+        if (str==(TString)"area")
+            fHadronnessCut = kArea;
+        if (str==(TString)"m3long")
+            fHadronnessCut = kM3Long;
+        if (str==(TString)"all")
+            fHadronnessCut = kAll;
+
+        rc = kTRUE;
+    }
+
+    if (IsEnvDefined(env, prefix, "File", print))
+    {
+        const TString fname = GetEnvValue(env, prefix, "File", "");
+        if (!CoefficentsRead(fname))
+            return kERROR;
+
+        return kTRUE;
+    }
+
+    for (int i=0; i<fVariables.GetSize(); i++)
+    {
+        if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
+        {
+            // It is important that the last argument is a floating point
+            fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
+            rc = kTRUE;
+        }
+    }
+    return rc;
+    //return kTRUE; // means: can use default values
+    //return rc;  // means: require something in resource file
+}

trunk/MagicSoft/Mars/mfilter/MFMagicCuts.h

@@ -6 +6 @@
 #endif
 
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+class MParList;
+
+class MHillas;
+class MHillasSrc;
+class MHillasExt;
+class MParameterD;
+class MPointingPos;
 class MHMatrix;
 
 class MFMagicCuts : public MFilter
 {
+public:
+    // Possible kind of theta cuts
+    enum ThetaCut_t {
+        kNone  =BIT(0),
+        kOn    =BIT(1),
+        kOff   =BIT(2),
+        kWobble=kOn|kOff
+    };
+    // Possible kind of hadronness cuts
+    enum HadronnessCut_t {
+        kNoCut =BIT(0),
+        kArea  =BIT(1),
+        kM3Long=BIT(2),
+        kAll   =kArea|kM3Long
+    };
+
 private:
+    // Elements for mapping. kLastElement must not be used and must be
+    // the last on in the list. It is used as counter for fMap.
+    enum {
+        kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist,
+        kEM3Long, kEM3LongAnti, kEDistAnti, kEWdivL, kEZd,
+        kLastElement
+    };
 
-    virtual Bool_t IsExpressionTrue() const { return kTRUE; }
+    MHillas        *fHil;              //! Pointer to MHillas container
+    MHillasSrc     *fHilSrc;           //! Pointer to MHillasSrc container
+    MHillasSrc     *fHilAnti;          //! Pointer to MHillasSrc container called MHillasSrcAnti
+    MHillasExt     *fHilExt;           //! Pointer to MHillasExt container
+    MParameterD    *fThetaSq;          //! Pointer to MParameterD container called ThetaSq
+    MParameterD    *fDisp;             //! Pointer to MParameterD container called Disp
+    MPointingPos   *fPointing;         //! Pointer to MPointingPos container
+
+    Float_t         fMm2Deg;            //! Conversion factor from mm to deg, from MGeomCam
+    Bool_t          fResult;            //! Result of the filter evaluation
+
+    Int_t           fMap[kLastElement]; //! Mapping table for fast optimization
+    MHMatrix       *fMatrix;            //! Matrix thorugh which the mapped elements are accessed
+
+    TArrayD         fVariables;         // Coefficients of cuts
+
+    ThetaCut_t      fThetaCut;          // Which kind of theta cut should be evaluated
+    HadronnessCut_t fHadronnessCut;     // Which kind of hadronness cut should be evaluated
+
+    // MTask
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
+
+    // MFilter
+    Bool_t IsExpressionTrue() const { return fResult; }
+
+    // MFMagicCuts
+    Double_t GetVal(Int_t i) const;
+    TString  GetParam(Int_t i) const;
+    Double_t GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par=0) const;
 
 public:
     MFMagicCuts(const char *name=NULL, const char *title=NULL);
 
-    void InitMapping(MHMatrix *mat) { }
-    void StopMapping() { InitMapping(NULL); }
+    // Getter
+    Double_t GetThetaSqCut() const;
 
-    ClassDef(MFMagicCuts, 0) // A filter to evaluate the MagicCuts
+    // Setter
+    void   SetThetaCut(ThetaCut_t c) { fThetaCut=c; }
+    void   SetHadronnessCut(HadronnessCut_t c) { fHadronnessCut=c; }
+
+    // MFMagicCuts
+    void   InitMapping(MHMatrix *mat);
+    void   StopMapping() { InitMapping(NULL); }
+
+    Bool_t CoefficentsRead(const char *fname);
+    Bool_t CoefficentsWrite(const char *fname) const;
+
+    // MParContainer
+    void   SetVariables(const TArrayD &arr);
+
+    Int_t  ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
+
+    // TObject
+    void   Print(Option_t *o="") const;
+
+    ClassDef(MFMagicCuts, 1) // A filter to evaluate the MagicCuts
 };
 

trunk/MagicSoft/Mars/mhbase/MH.cc

@@ -1204 +1204 @@
 // --------------------------------------------------------------------------
 //
-// M.Gaug added this withouz Documentation
+// 
 //
 TH1I* MH::ProjectArray(const TArrayF &array, Int_t nbins, 
@@ -1251 +1251 @@
 // --------------------------------------------------------------------------
 //
-// M.Gaug added this withouz Documentation
+// 
 //
 TH1I* MH::ProjectArray(const TArrayD &array, Int_t nbins, const char* name, const char* title)
@@ -1297 +1297 @@
 // --------------------------------------------------------------------------
 //
-// M.Gaug added this withouz Documentation
+// 
 //
 TH1I* MH::ProjectArray(const MArrayF &array, Int_t nbins, 
@@ -1307 +1307 @@
 // --------------------------------------------------------------------------
 //
-// M.Gaug added this withouz Documentation
+// 
 //
 TH1I* MH::ProjectArray(const MArrayD &array, Int_t nbins, const char* name, const char* title)
@@ -1324 +1324 @@
     *fLog << "%) Evts skipped: " << str << endl;
 }
+
+// --------------------------------------------------------------------------
+//
+// Calls gStyle->SetPalette. Allowed palettes are:
+//  pretty
+//  deepblue:  darkblue -> lightblue
+//  lightblue: black -> blue -> white
+//  greyscale: black -> white
+//  glow1:     black -> darkred -> orange -> yellow -> white
+//  glow2:
+//  glowsym:   lightblue -> blue -> black -> darkred -> orange -> yellow -> white
+//  redish:    darkred -> lightred
+//  bluish:    darkblue -> lightblue
+//  small1:
+//
+// If the palette name contains 'inv' the order of the colors is inverted.
+//
+// The second argument determines the number of colors for the palette.
+// The default is 50. 'pretty' always has 50 colors.
+//
+// (Remark: Drawing 3D object like TH2D with surf3 allows a maximum
+//          of 99 colors)
+//
+void MH::SetPalette(TString paletteName, Int_t ncol)
+{
+    Bool_t found=kFALSE;
+
+    paletteName.ToLower();
+
+    const Bool_t inverse = paletteName.Contains("inv");
+
+    if (paletteName.Contains("pretty"))
+    {
+        gStyle->SetPalette(1, 0);
+        ncol=50;
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("deepblue"))
+    {
+        Double_t s[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+        Double_t r[5] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
+        Double_t g[5] = { 0.01, 0.02, 0.39, 0.68, 0.97 };
+        Double_t b[5] = { 0.17, 0.39, 0.62, 0.79, 0.97 };
+        gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("lightblue"))
+    {
+        Double_t s[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+        Double_t r[5] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
+        Double_t g[5] = { 0.00, 0.02, 0.40, 0.70, 1.00 };
+        Double_t b[5] = { 0.00, 0.27, 0.51, 0.81, 1.00 };
+        gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("greyscale"))
+    {
+        double s[2] = {0.00, 1.00};
+        double r[2] = {0.00, 1.00};
+        double g[2] = {0.00, 1.00};
+        double b[2] = {0.00, 1.00};
+        gStyle->CreateGradientColorTable(2, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("glow1"))
+    {
+        double s[5] = {0., 0.10, 0.45, 0.75, 1.00};
+        double r[5] = {0., 0.35, 0.85, 1.00, 1.00};
+        double g[5] = {0., 0.10, 0.20, 0.73, 1.00};
+        double b[5] = {0., 0.03, 0.06, 0.00, 1.00};
+        gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("glow2"))
+    {
+        double s[4] = {0.00, 0.50, 0.75, 1.00};
+        double r[4] = {0.24, 0.67, 1.00, 1.00};
+        double g[4] = {0.03, 0.04, 0.80, 1.00};
+        double b[4] = {0.03, 0.04, 0.00, 1.00};
+        gStyle->CreateGradientColorTable(4, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("glowsym"))
+    {
+        double s[8] = {0.00, 0.17, 0.39, 0.50, 0.55, 0.72, 0.88, 1.00};
+        double r[8] = {0.09, 0.18, 0.09, 0.00, 0.35, 0.85, 1.00, 1.00};
+        double g[8] = {0.70, 0.40, 0.02, 0.00, 0.10, 0.20, 0.73, 1.00};
+        double b[8] = {0.81, 0.51, 0.27, 0.00, 0.03, 0.06, 0.00, 1.00};
+        gStyle->CreateGradientColorTable(8, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("redish"))
+    {
+        double s[3] = {0., 0.5, 1.};
+        double r[3] = {0., 1.0, 1.};
+        double g[3] = {0., 0.0, 1.};
+        double b[3] = {0., 0.0, 1.};
+        gStyle->CreateGradientColorTable(3, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("bluish"))
+    {
+        double s[3] = {0., 0.5, 1.};
+        double r[3] = {0., 0.0, 1.};
+        double g[3] = {0., 0.0, 1.};
+        double b[3] = {0., 1.0, 1.};
+        gStyle->CreateGradientColorTable(3, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (paletteName.Contains("small1"))
+    {
+        double s[4] = {0.00, 0.50, 0.95, 1.};
+        double r[4] = {0.04, 0.28, 0.98, 1.};
+        double g[4] = {0.28, 0.93, 0.03, 1.};
+        double b[4] = {0.79, 0.11, 0.03, 1.};
+        gStyle->CreateGradientColorTable(4, s, r, g, b, ncol);
+        found=kTRUE;
+    }
+
+    if (inverse)
+    {
+        TArrayI c(ncol);
+        for (int i=0; i<ncol; i++)
+            c[ncol-i-1] = gStyle->GetColorPalette(i);
+        gStyle->SetPalette(ncol, c.GetArray());
+    }
+
+    if (!found)
+        gLog << warn << "MH::SetPalette: Palette " << paletteName << " unknown... ignored." << endl;
+}

trunk/MagicSoft/Mars/mhbase/MH.h

@@ -117 +117 @@
     static TObject *FindObjectInPad(const char *name, TVirtualPad *pad=NULL);
 
+    static void SetPalette(TString paletteName, Int_t ncol=50);
+
     ClassDef(MH, 2) //A base class for Mars histograms
 };

trunk/MagicSoft/Mars/mhflux/MHDisp.cc

@@ -64 +64 @@
 //
 MHDisp::MHDisp(const char *name, const char *title)
-    : fHilExt(0), fDisp(0)
+    : fHilExt(0), fDisp(0)//, fIsWobble(kFALSE)
 {
     //
@@ -78 +78 @@
     fHist.SetXTitle("x [\\circ]");
     fHist.SetYTitle("y [\\circ]");
-    fHist.SetZTitle("\\vartheta [deg^{2}]");
+    fHist.SetZTitle("Eq.cts");
 }
 
@@ -130 +130 @@
         return kFALSE;
     }
+
+    //const Double_t xmax = fHist.GetXaxis()->GetXmax();
 
     // Initialize all bins with a small (=0) value otherwise
     // these bins are not displayed
     for (int x=0; x<fHist.GetNbinsX(); x++)
-        for (int y=0; y<fHist.GetNbinsX(); y++)
-            fHist.Fill(fHist.GetXaxis()->GetBinCenter(x+1),
-                       fHist.GetYaxis()->GetBinCenter(y+1),
-                       0.0, 1e-10);
+        for (int y=0; y<fHist.GetNbinsY(); y++)
+        {
+            const Double_t cx = fHist.GetXaxis()->GetBinCenter(x+1);
+            const Double_t cy = fHist.GetYaxis()->GetBinCenter(y+1);
+            //if (TMath::Hypot(cx, cy)<xmax)
+                fHist.Fill(cx, cy, 0.0, 1e-10);
+        }
 
     return kTRUE;
@@ -160 +165 @@
         rho = fPointPos->RotationAngle(*fObservatory, *fTime);
 
+    // FIXME: Do wobble-rotation when subtracting?
+    if (!fHistOff/* && fIsWobble*/)
+        rho += TMath::Pi();
+
     // Get Disp from Parlist
     const Double_t disp = fDisp->GetVal();
 
-    // Calculate where disp is pointing
+    // Calculate both postiions where disp could point
     TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta());
     TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta());
-    pos1 *=  disp;
-    pos2 *= -disp;
-
+    pos1 *=  disp; // Vector of length  disp in direction of shower
+    pos2 *= -disp; // Vector of length -disp in direction of shower
+
+    // Move origin of vector to center-of-gravity of shower
     pos1 += hil->GetMean()*fMm2Deg;
     pos2 += hil->GetMean()*fMm2Deg;
@@ -207 +217 @@
         srcp = fSrcPos->GetXY();
 
+    // Derotate all position around camera center by -rho
     if (rho!=0)
     {
@@ -214 +225 @@
     }
 
+    // Shift the source position to 0/0
     pos1 -= srcp*fMm2Deg;
     pos2 -= srcp*fMm2Deg;
 
-    fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight);
-    fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight));
+    //const Double_t xmax = fHist.GetXaxis()->GetXmax();
+
+    // Fill histograms
+    //if (pos1.Mod()<xmax)
+        fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight);
+    //if (pos2.Mod()<xmax)
+        fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight));
 
     return kTRUE;
 }
-/*
-static Double_t FcnGauss2d(Double_t *x, Double_t *par)
-{
-    TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad());
-
-    const Double_t g0 = TMath::Gaus(v.X(), par[1], par[2]);
-    const Double_t g1 = TMath::Gaus(v.Y(), par[3], par[4]);
-
-    //Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
-    return par[0]*g0*g1  + par[6]*(v.X()*v.X() + v.Y()*v.Y()) +  par[7];
-}*/
+
+Double_t MHDisp::GetOffSignal(TH1 &h) const
+{
+    const TAxis &axex = *h.GetXaxis();
+    const TAxis &axey = *h.GetYaxis();
+
+    Double_t sum = 0;
+    for (int x=0; x<h.GetNbinsX(); x++)
+        for (int y=0; y<h.GetNbinsY(); y++)
+        {
+            if (TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1))>0.35)
+                sum += h.GetBinContent(x+1,y+1);
+        }
+
+    return sum;
+}
 
 void MHDisp::Paint(Option_t *o)
@@ -240 +262 @@
     pad->cd(1);
 
+    // Project on data onto yx-plane
     fHist.GetZaxis()->SetRange(0,9999);
     TH1 *h1=fHist.Project3D("yx_on");
-    gStyle->SetPalette(1, 0);
+
+    // Set Glow-palette (PaintSurface doesn't allow more than 99 colors)
+    MH::SetPalette(fHistOff?"glowsym":"glow1", 99);
+    h1->SetContour(99);
+
+    Double_t scale = 1;
+    if (fHistOff)
+    {
+        // Project off data onto yx-plane and subtract it from on-data
+        fHistOff->GetZaxis()->SetRange(0,9999);
+        TH1 *h=fHistOff->Project3D("yx_off");
+
+        scale = -1;
+
+        //if (!fIsWobble)
+        {
+            const Double_t h1off = GetOffSignal(*h1);
+            const Double_t hoff  = GetOffSignal(*h);
+            scale = hoff==0?-1:-h1off/hoff;
+        }
+
+        h1->Add(h, scale);
+        delete h;
+
+        // Force calculation of minimum, maximum
+        h1->SetMinimum();
+        h1->SetMaximum();
+
+        // Find maximum
+        const Double_t max = TMath::Max(TMath::Abs(h1->GetMinimum()),
+                                        TMath::Abs(h1->GetMaximum()));
+
+        // Set new minimum, maximum
+        h1->SetMinimum(-max);
+        h1->SetMaximum( max);
+    }
 
     Int_t ix, iy, iz;
@@ -304 +362 @@
         const Double_t s  = MMath::SignificanceLiMa(e, b);
 
-        h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f", x0, y0, func.GetParameter(2), s, e-b, b));
+        h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f f=%.2f", x0, y0, func.GetParameter(2), s, e-b, b, scale));
     }
    /*
@@ -386 +444 @@
     h->Draw();
 }
+
+// --------------------------------------------------------------------------
+//
+// The following resources are available:
+//
+//    MHDisp.Wobble: on/off
+//
+/*
+Int_t MHFalseSource::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "DistMin", print))
+    {
+        rc = kTRUE;
+        SetMinDist(GetEnvValue(env, prefix, "DistMin", fMinDist));
+    }
+    if (IsEnvDefined(env, prefix, "DistMax", print))
+    {
+        rc = kTRUE;
+        SetMaxDist(GetEnvValue(env, prefix, "DistMax", fMaxDist));
+    }
+
+    if (IsEnvDefined(env, prefix, "DWMin", print))
+    {
+        rc = kTRUE;
+        SetMinDW(GetEnvValue(env, prefix, "DWMin", fMinDist));
+    }
+    if (IsEnvDefined(env, prefix, "DWMax", print))
+    {
+        rc = kTRUE;
+        SetMaxDW(GetEnvValue(env, prefix, "DWMax", fMaxDist));
+    }
+
+    return rc;
+}
+*/

trunk/MagicSoft/Mars/mhflux/MHDisp.h

@@ -13 +13 @@
 {
 private:
-    MHillasExt  *fHilExt; //!
-    MParameterD *fDisp;   //!
+    MHillasExt  *fHilExt;  //!
+    MParameterD *fDisp;    //!
 
-    Double_t fM3lCut;     //!
-    Double_t fXi;         //!
-    Double_t fXiTheta;    //!
+    Double_t     fM3lCut;  //!
+    Double_t     fXi;      //!
+    Double_t     fXiTheta; //!
+
+    //Bool_t       fIsWobble;
+
+    // MHDisp
+    Double_t GetOffSignal(TH1 &h) const;
 
 public:
@@ -27 +32 @@
     Bool_t Fill(const MParContainer *par, const Stat_t w=1);
 
+    // TObject
     void Paint(Option_t *o="");
     void Draw(Option_t *o="");
+
+    // MParContainer
+    //Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
 
     ClassDef(MHDisp, 1) //3D-histogram in alpha, x and y