Changeset 7193

Timestamp:

07/14/05 18:45:35 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (1 diff)
• mhflux/MHDisp.cc (modified) (15 diffs)
• mhflux/MHDisp.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -45 +45 @@
      - comment have been updated.
      - some old code had been removed
+
+   * mhflux/MHDisp.[h,cc]:
+     - updated with a more appropriate calculation of a background model
 
 

trunk/MagicSoft/Mars/NEWS

@@ -115 +115 @@
      more reliable in terms of gain-fluctuations and sudden VCSEL gain
      changes.
+
+   - ganymed: The false source plot (MHDisp) is now based on Disp
+     and a background model determied in the first loop is
+     subtracted
 
    - sponde: the zenith angle distribution is now weighted instead of

trunk/MagicSoft/Mars/mhflux/MHDisp.cc

@@ -46 +46 @@
 
 #include "MMath.h"
+#include "MBinning.h"
 
 #include "MParList.h"
@@ -64 +65 @@
 //
 MHDisp::MHDisp(const char *name, const char *title)
-    : fHilExt(0), fDisp(0)//, fIsWobble(kFALSE)
+    : fHilExt(0), fDisp(0), fSmearing(-1), fWobble(kFALSE)
 {
     //
@@ -72 +73 @@
     fTitle = title ? title : "3D-plot using Disp vs x, y";
 
-    fHist.SetDirectory(NULL);
-
     fHist.SetName("Alpha");
     fHist.SetTitle("3D-plot of ThetaSq vs x, y");
@@ -79 +78 @@
     fHist.SetYTitle("y [\\circ]");
     fHist.SetZTitle("Eq.cts");
+
+    fHist.SetDirectory(NULL);
+    fHistBg.SetDirectory(NULL);
+    fHistBg1.SetDirectory(NULL);
+    fHistBg2.SetDirectory(NULL);
+
+    fHist.SetBit(TH1::kNoStats);
 }
 
@@ -100 +106 @@
     }
 
-    MParameterD *p = (MParameterD*)plist->FindObject("M3lCut", "MParameterD");
-    if (!p)
-    {
-        *fLog << err << "M3lCut [MParameterD] not found... abort." << endl;
-        return kFALSE;
-    }
-    fM3lCut = TMath::Abs(p->GetVal());
-
-    p = (MParameterD*)plist->FindObject("DispXi", "MParameterD");
-    if (!p)
-    {
-        *fLog << err << "DispXi [MParameterD] not found... abort." << endl;
-        return kFALSE;
-    }
-    fXi = p->GetVal();
-
-    p = (MParameterD*)plist->FindObject("DispXiTheta", "MParameterD");
-    if (!p)
-    {
-        *fLog << err << "DispXiTheta [MParameterD] not found... abort." << endl;
-        return kFALSE;
-    }
-    fXiTheta = p->GetVal();
-
     fHilExt = (MHillasExt*)plist->FindObject("MHillasExt");
     if (!fHilExt)
@@ -131 +113 @@
     }
 
-    //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.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);
-        }
+    fHilSrc = (MHillasSrc*)plist->FindObject("MHillasSrc");
+    if (!fHilSrc)
+    {
+        *fLog << err << "MHillasSrc not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    MBinning binsx, binsy;
+    binsx.SetEdges(fHist, 'x');
+    binsy.SetEdges(fHist, 'y');
+    MH::SetBinning(&fHistBg, &binsx, &binsy);
+
+    if (!fHistOff)
+    {
+        MH::SetBinning(&fHistBg1, &binsx, &binsy);
+        MH::SetBinning(&fHistBg2, &binsx, &binsy);
+    }
 
     return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Calculate the delta angle between fSrcPos->GetXY() and v.
+// Return result in deg.
+//
+Double_t MHDisp::DeltaPhiSrc(const TVector2 &v) const
+{
+    return TMath::Abs(fSrcPos->GetXY().DeltaPhi(v))*TMath::RadToDeg();
 }
 
@@ -165 +162 @@
         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 both postiions where disp could point
-    TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta());
-    TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta());
-    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;
-
-    // gammaweight: If we couldn't decide which position makes the
-    // event a gamma, both position are assigned 'half of a gamma'
-    Double_t gweight = 0.5;
-
-    // Check whether our g/h-separation allows to asign definitly
-    // to one unique position. Therefor we requeire that the event
-    // would be considered a gamma for one, but not for the other
-    // position. This can only be the case if the third moment
-    // has a value higher than the absolute cut value.
-    if (TMath::Abs(fHilExt->GetM3Long()*fMm2Deg) > fM3lCut)
-    {
-        // Because at one position the event is considered a gamma
-        // we have to find out which position it is...
-        MSrcPosCam src;
-        MHillasSrc hsrc;
-        hsrc.SetSrcPos(&src);
-
-        // Calculate the sign for one of the desired source positions
-        // The other position must have the opposite sign
-        src.SetXY(pos1/fMm2Deg);
-        if (hsrc.Calc(*hil)>0)
-        {
-            *fLog << warn << "Calculation of MHillasSrc failed." << endl;
-            return kFALSE;
-        }
-        const Double_t m3l = fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, hsrc.GetCosDeltaAlpha());
-
-        gweight = m3l>fM3lCut ? 1 : 0;
-    }
-
-    // Now we can safly derotate both position...
-    TVector2 srcp;
-    if (fSrcPos)
-        srcp = fSrcPos->GetXY();
-
-    // Derotate all position around camera center by -rho
-    if (rho!=0)
-    {
+    // Vector of length disp in direction of shower
+    const TVector2 p(hil->GetCosDelta(), hil->GetSinDelta());
+
+    // Move origin of vector to center-of-gravity of shower and derotate
+    TVector2 pos1 = hil->GetMean()*fMm2Deg + p*fDisp->GetVal();
+
+    Double_t w0 = 1;
+    if (fWobble)
+    {
+        const Double_t delta = DeltaPhiSrc(pos1);
+
+        // Skip off-data not in the same half than the source (here: anti-source)
+        // Could be replaced by some kind of theta cut...
+        if (!fHistOff)
+        {
+            if (delta>165)
+                return kTRUE;
+
+            // Because afterwards the plots are normalized with the number
+            // of entries the  non-overlap  region corresponds to half the
+            // contents, the overlap region  to  full contents.  By taking
+            // the mean of both distributions  we get the same result than
+            // we would have gotten using full  off-events with a slightly
+            // increased uncertainty
+            // FIXME: The delta stuff could be replaced by a 2*antitheta cut...
+            w0 = delta>15 ? 1 : 2;
+        }
+
+        // Define by the source position which histogram to fill
+        if (DeltaPhiSrc(fFormerSrc)>90)
+            fHalf = !fHalf;
+        fFormerSrc = fSrcPos->GetXY();
+    }
+
+    // If on-data: Derotate source and P. Translate P to center.
+    TVector2 m;
+    if (fHistOff)
+    {
+        // Derotate all position around camera center by -rho
+        // Move origin of vector to center-of-gravity of shower and derotate
         pos1=pos1.Rotate(-rho);
-        pos2=pos2.Rotate(-rho);
-        srcp=srcp.Rotate(-rho);
-    }
-
-    // Shift the source position to 0/0
-    pos1 -= srcp*fMm2Deg;
-    pos2 -= srcp*fMm2Deg;
-
-    //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));
+
+        // Shift the source position to 0/0
+        if (fSrcPos)
+        {
+            // m: Position of the camera center in the FS plot
+            m = fSrcPos->GetXY().Rotate(-rho)*fMm2Deg;
+            pos1 -= m;
+        }
+    }
+
+    if (fSmearing<=0)
+    {
+        fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*w0);
+        return kTRUE;
+    }
+
+    // -------------------------------------------------
+    //  The following algorithm may look complicated...
+    //            It is optimized for speed
+    // -------------------------------------------------
+
+    // Define a vector used to calculate rotated x and y component
+    const TVector2 rot(TMath::Sin(rho), TMath::Cos(rho));
+
+    // Fold event position with the PSF and fill it into histogram
+    const TAxis &axex = *fHist.GetXaxis();
+    const TAxis &axey = *fHist.GetYaxis();
+
+    const Int_t nx = axex.GetNbins();
+    const Int_t ny = axey.GetNbins();
+
+    // normg: Normalization for Gauss [sqrt(2*pi)*sigma]^2
+    const Double_t weight = axex.GetBinWidth(1)*axey.GetBinWidth(1)*w*w0;
+    const Double_t psf    = 2*fSmearing*fSmearing;
+    const Double_t pi2    = fSmearing * TMath::Pi()/2;
+    const Double_t normg  = pi2*pi2 * TMath::Sqrt(TMath::TwoPi()) / weight;
+    const Int_t    bz     = fHist.GetZaxis()->FindFixBin(0);
+
+    TH1 &bg = fHalf ? fHistBg1 : fHistBg2;
+
+    // To calculate significance map smear with 2*theta-cut and
+    // do not normalize gauss, for event map use theta-cut/2 instead
+    for (int x=1; x<=nx; x++)
+    {
+        const Double_t cx = axex.GetBinCenter(x);
+        const Double_t px = cx-pos1.X();
+
+        for (int y=1; y<=ny; y++)
+        {
+            const Double_t cy = axey.GetBinCenter(y);
+            const Double_t sp = Sq(px, cy-pos1.Y());
+            const Double_t dp = sp/psf;
+
+            // Values below 1e-3 (>3.5sigma) are not filled into the histogram
+            if (dp<4)
+            {
+                const Double_t rc = TMath::Exp(-dp)/normg;
+                if (!fHistOff)
+                    bg.AddBinContent(bg.GetBin(x, y), rc);
+                else
+                    fHist.AddBinContent(fHist.GetBin(x, y, bz), rc);
+            }
+
+            if (!fHistOff)
+                continue;
+
+            // If we are filling the signal already (fHistOff!=NULL)
+            // we also fill the background by projecting the
+            // background in the camera into the sky plot.
+            const TVector2 v = TVector2(cx+m.X(), cy+m.Y());
+
+            // Speed up further: Xmax-fWobble
+            if (v.Mod()>axex.GetXmax()) // Gains 10% speed
+                continue;
+
+            const Int_t    bx = axex.FindFixBin(v^rot);
+            const Int_t    by = axey.FindFixBin(v*rot);
+            const Double_t bg = fHistOff->GetBinContent(bx, by, bz);
+
+            fHistBg.AddBinContent(fHistBg.GetBin(x, y), bg);
+        }
+    }
+    fHist.SetEntries(fHist.GetEntries()+1);
+
+    if (!fHistOff)
+        bg.SetEntries(fHistBg1.GetEntries()+1);
 
     return kTRUE;
 }
 
+// --------------------------------------------------------------------------
+//
+// Compile the background in camera coordinates from the two background
+// histograms
+//
+Bool_t MHDisp::Finalize()
+{
+    if (fHistOff)
+        return kTRUE;
+
+    const Int_t bz = fHist.GetZaxis()->FindFixBin(0);
+
+    const Double_t n1 = fHistBg1.GetEntries();
+    const Double_t n2 = fHistBg2.GetEntries();
+
+    for (int x=1; x<=fHist.GetNbinsX(); x++)
+        for (int y=1; y<=fHist.GetNbinsY(); y++)
+        {
+            const Int_t bin1 = fHistBg1.GetBin(x, y);
+
+            const Double_t rc =
+                (n1==0?0:fHistBg1.GetBinContent(bin1)/n1)+
+                (n2==0?0:fHistBg2.GetBinContent(bin1)/n2);
+
+            fHist.SetBinContent(x, y, bz, rc/2);
+        }
+
+    fHist.SetEntries(n1+n2);
+
+    // Result corresponds to one smeared background event
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Return the mean signal in h around (0,0) in a distance between
+// 0.33 and 0.47deg
+//
 Double_t MHDisp::GetOffSignal(TH1 &h) const
 {
@@ -245 +338 @@
     const TAxis &axey = *h.GetYaxis();
 
+    Int_t    cnt = 0;
     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)
+            const Double_t d = TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1));
+            if (d>0.33 && d<0.47)
+            //if (d>0.4 && d<0.8)
+            {
                 sum += h.GetBinContent(x+1,y+1);
-        }
-
-    return sum;
-}
+                cnt++;
+            }
+        }
+
+    return sum/cnt;
+}
+
+// --------------------------------------------------------------------------
+//
+// Fill h2 with the radial profile of h1 around (x0, y0)
+//
+void MHDisp::Profile(TH1 &h2, const TH2 &h1, Axis_t x0, Axis_t y0) const
+{
+    const TAxis &axex = *h1.GetXaxis();
+    const TAxis &axey = *h1.GetYaxis();
+
+    h2.Reset();
+
+    for (int x=1; x<=axex.GetNbins(); x++)
+        for (int y=1; y<=axey.GetNbins(); y++)
+        {
+            const Double_t dx = axex.GetBinCenter(x)-x0;
+            const Double_t dy = axey.GetBinCenter(y)-y0;
+
+            const Double_t r  = TMath::Hypot(dx, dy);
+
+            h2.Fill(r, h1.GetBinContent(x, y));
+        }
+}
+
+// --------------------------------------------------------------------------
+//
+// Remove contents of histogram around a circle.
+//
+void MHDisp::MakeDot(TH2 &h2) const
+{
+    const TAxis &axex = *h2.GetXaxis();
+    const TAxis &axey = *h2.GetYaxis();
+
+    const Double_t rmax = fWobble ? axex.GetXmax()-0.4 : axex.GetXmax();
+
+    for (int x=1; x<=axex.GetNbins(); x++)
+        for (int y=1; y<=axey.GetNbins(); y++)
+        {
+            const Int_t bin = h2.GetBin(x,y);
+
+            const Double_t px = h2.GetBinCenter(x);
+            const Double_t py = h2.GetBinCenter(y);
+
+            if (rmax<TMath::Hypot(px, py))
+                h2.SetBinContent(bin, 0);
+        }
+}
+
+// --------------------------------------------------------------------------
+//
+// Calculate from signal and background the significance map
+//
+void MHDisp::MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale) const
+{
+    const TAxis &axex = *s.GetXaxis();
+    const TAxis &axey = *s.GetYaxis();
+
+    for (int x=1; x<=axex.GetNbins(); x++)
+        for (int y=1; y<=axey.GetNbins(); y++)
+        {
+            const Int_t  bin = s.GetBin(x,y);
+
+            const Double_t sig = h1.GetBinContent(bin);
+            const Double_t bg  = h2.GetBinContent(bin);
+
+            const Double_t S = MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale));
+
+            s.SetBinContent(bin, S);
+        }
+}
+
 
 void MHDisp::Paint(Option_t *o)
 {
+    // Compile Background if necessary
+    Finalize();
+
+    // Paint result
     TVirtualPad *pad = gPad;
 
@@ -264 +438 @@
     // Project on data onto yx-plane
     fHist.GetZaxis()->SetRange(0,9999);
-    TH1 *h1=fHist.Project3D("yx_on");
+    TH2 *h1=(TH2*)fHist.Project3D("yx_on");
 
     // Set Glow-palette (PaintSurface doesn't allow more than 99 colors)
@@ -275 +449 @@
         // Project off data onto yx-plane and subtract it from on-data
         fHistOff->GetZaxis()->SetRange(0,9999);
-        TH1 *h=fHistOff->Project3D("yx_off");
+        TH2 *h=(TH2*)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);
+        const Double_t h1off = GetOffSignal(*h1);
+        const Double_t hoff  = GetOffSignal(fHistBg);
+
+        scale = hoff==0 ? -1 : -h1off/hoff;
+
+        const Bool_t sig = kFALSE;
+
+        if (!sig)
+            h1->Add(&fHistBg, scale);
+        else
+            MakeSignificance(*h1, *h1, fHistBg, scale);
+
+        MakeDot(*h1);
+
         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;
-    TF1 func("fcn", "gaus + [3]*x*x + [4]");
+        MakeSymmetric(h1);
+    }
 
     pad->cd(3);
     TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
-    /*
-    if (!h2)
-    {
-        const Double_t maxr = TMath::Hypot(h1->GetXaxis()->GetXmax(), h1->GetYaxis()->GetXmax());
-        const Int_t    nbin = (h1->GetNbinsX()+h1->GetNbinsY())/2;
-        TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
-        //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
-        h->Sumw2();
-        h->SetXTitle("\\vartheta [\\circ]");
-        h->SetYTitle("<cts>/\\Delta R");
-        h->SetBit(kCanDelete);
-        h->Draw();
-    }*/
+
+    TString opt(o);
+    opt.ToLower();
+
     if (h1 && h2)
     {
-        h2->Reset();
-
+        Int_t ix, iy, iz;
         h1->GetMaximumBin(ix, iy, iz);
 
         const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
         const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
-        const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);
-
-        for (int x=0; x<h1->GetNbinsX(); x++)
-            for (int y=0; y<h1->GetNbinsY(); y++)
-            {
-                const Double_t r = TMath::Hypot(h1->GetXaxis()->GetBinCenter(x+1)-x0,
-                                                h1->GetYaxis()->GetBinCenter(y+1)-y0);
-                h2->Fill(r, h1->GetBinContent(x+1, y+1));
-            }
+        //const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);
+
+        Profile(*h2, *h1, 0, 0);
+        //Profile(*h2, *h1, x0, y0);
 
         // Replace with MAlphaFitter?
+        TF1 func("fcn", "gaus + [3]*x*x + [4]");
         func.SetLineWidth(1);
         func.SetLineColor(kBlue);
@@ -347 +500 @@
         func.FixParameter(1, 0);
         func.SetParameter(2, 0.15);
-        func.SetParameter(4, h2->GetBinContent(10));
-        h2->Fit(&func, "IMQ", "", 0, 1.0);
-
-        // No wintegrate the function f(x) per Delta Area
-        // which is f(x)/(pi*delta r*(2*r+delta r))
-        TF1 func2("fcn2", Form("(gaus + [3]*x*x + [4])/(2*x+%.5f)", w0));
-        for (int i=0; i<5; i++)
-            func2.SetParameter(i, func.GetParameter(i));
-
-        const Double_t r0 = 2*func2.GetParameter(2);
-        const Double_t e  = func2.Integral(0, r0)/(w0*TMath::Pi());
-        func2.SetParameter(0, 0);
-        const Double_t b  = func2.Integral(0, r0)/(w0*TMath::Pi());
-        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 f=%.2f", x0, y0, func.GetParameter(2), s, e-b, b, scale));
-    }
-   /*
-    if (h1)
-    {
-        const Double_t maxr = 0.9*TMath::Abs(fHist.GetXaxis()->GetXmax());
-
-        TF2 f2d("Gaus2D", FcnGauss2d, -maxr, maxr, -maxr, maxr, 8);
-        f2d.SetLineWidth(1);
-
-        f2d.SetParameter(0, h1->GetMaximum()*5); // A
-        f2d.SetParLimits(1, h1->GetXaxis()->GetBinCenter(ix)-h1->GetXaxis()->GetBinWidth(ix)*5, h1->GetXaxis()->GetBinCenter(ix)+h1->GetXaxis()->GetBinWidth(ix));  // mu_1
-        f2d.SetParLimits(3, h1->GetYaxis()->GetBinCenter(iy)-h1->GetYaxis()->GetBinWidth(iy)*5, h1->GetYaxis()->GetBinCenter(iy)+h1->GetYaxis()->GetBinWidth(iy));  // mu_2
-        f2d.SetParLimits(2, 0, func.GetParameter(2)*5);          // sigma_1
-        f2d.SetParLimits(4, 0, func.GetParameter(2)*5);          // sigma_2
-        f2d.SetParLimits(5, 0, 45);            // phi
-        f2d.SetParLimits(6, 0, func.GetParameter(3)*5);
-        f2d.SetParLimits(7, 0, func.GetParameter(4)*5);
-
-        f2d.SetParameter(0, h1->GetMaximum()); // A
-        f2d.SetParameter(1, h1->GetXaxis()->GetBinCenter(ix)); // mu_1
-        f2d.SetParameter(2, func.GetParameter(2));             // sigma_1
-        f2d.SetParameter(3, h1->GetYaxis()->GetBinCenter(iy)); // mu_2
-        f2d.SetParameter(4, func.GetParameter(2));             // sigma_2
-        f2d.FixParameter(5, 0);                                // phi
-        f2d.SetParameter(6, func.GetParameter(3));
-        f2d.SetParameter(7, func.GetParameter(4));
-        h1->Fit(&f2d, "Q", "cont2");
-        //f2d.DrawCopy("cont2same");
-    }*/
+        if (fHistOff)
+            func.FixParameter(3, 0);
+        func.SetParameter(4, h2->GetBinContent(15));
+        h2->Fit(&func, "IQ", "", 0, 1.0);
+
+        h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ f=%.2f", x0, y0, func.GetParameter(2), TMath::Abs(scale)));
+    }
 }
 
@@ -400 +515 @@
     pad->SetBorderMode(0);
 
-    AppendPad("");
+    AppendPad(o);
 
     TString name = Form("%s_1", pad->GetName());
@@ -416 +531 @@
     h3->Draw("colz");
     h3->SetBit(kCanDelete);
-//    catalog->Draw("mirror same *");
+
+    if (fHistOff)
+        GetCatalog()->Draw("white mirror same *");
 
     pad->cd();
@@ -430 +547 @@
     p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
     p->SetNumber(3);
+    p->SetGrid();
     p->Draw();
     p->cd();
@@ -449 +567 @@
 // 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))
+//    MHDisp.Smearing: 0.1
+//    MHDisp.Wobble:   on/off
+//
+Int_t MHDisp::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Int_t rc = MHFalseSource::ReadEnv(env, prefix, print);
+    if (rc==kERROR)
+        return kERROR;
+
+    if (IsEnvDefined(env, prefix, "Smearing", print))
     {
         rc = kTRUE;
-        SetMinDist(GetEnvValue(env, prefix, "DistMin", fMinDist));
-    }
-    if (IsEnvDefined(env, prefix, "DistMax", print))
+        SetSmearing(GetEnvValue(env, prefix, "Smearing", fSmearing));
+    }
+    if (IsEnvDefined(env, prefix, "Wobble", 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));
+        SetWobble(GetEnvValue(env, prefix, "Wobble", fWobble));
     }
 
     return rc;
 }
-*/
+

trunk/MagicSoft/Mars/mhflux/MHDisp.h

@@ -5 +5 @@
 #include "MHFalseSource.h"
 #endif
+#ifndef ROOT_TH2
+#include <TH2.h>
+#endif
+#ifndef ROOT_TVector2
+#include <TVector2.h>
+#endif
 
 class MParList;
 class MHillasExt;
+class MHillasSrc;
 class MParameterD;
 
@@ -13 +20 @@
 {
 private:
-    MHillasExt  *fHilExt;  //!
-    MParameterD *fDisp;    //!
+    MHillasExt   *fHilExt;  //!
+    MHillasSrc   *fHilSrc;  //!
+    MParameterD  *fDisp;    //!
 
-    Double_t     fM3lCut;  //!
-    Double_t     fXi;      //!
-    Double_t     fXiTheta; //!
+    TH2D         fHistBg;
 
-    //Bool_t       fIsWobble;
+    TH2D         fHistBg1;
+    TH2D         fHistBg2;
+
+    TVector2     fFormerSrc;
+    Bool_t       fHalf;
+
+    Double_t     fSmearing;
+    Bool_t       fWobble;
 
     // MHDisp
     Double_t GetOffSignal(TH1 &h) const;
+    Double_t DeltaPhiSrc(const TVector2 &v) const;
+
+    void Profile(TH1 &h2, const TH2 &h1, Axis_t x0=0, Axis_t y0=0) const;
+    void MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale=1) const;
+    void MakeDot(TH2 &h2) const;
+
+    Double_t Sq(Double_t x, Double_t y) const { return x*x+y*y; }
 
 public:
     MHDisp(const char *name=NULL, const char *title=NULL);
 
+    // MHDisp
+    void SetSmearing(Double_t s=-1) { fSmearing=s; }
+    void SetWobble(Bool_t w=kTRUE)  { fWobble=w;   }
+
+    // MHFalseSource (to be moved!)
+    void SetOffData(const MHFalseSource &fs)
+    {
+        MHFalseSource::SetOffData(fs);
+        if (dynamic_cast<const MHDisp*>(&fs))
+        {
+            const MHDisp &h = dynamic_cast<const MHDisp&>(fs);
+            fWobble   = h.fWobble;
+            fSmearing = h.fSmearing;
+
+            h.fHistBg1.Copy(fHistBg1);
+            h.fHistBg2.Copy(fHistBg2);
+
+            fHistBg1.SetDirectory(0);
+            fHistBg2.SetDirectory(0);
+        }
+    }
+
     // MH
     Bool_t SetupFill(const MParList *pList);
     Bool_t Fill(const MParContainer *par, const Stat_t w=1);
+    Bool_t Finalize();
 
     // TObject
@@ -37 +80 @@
 
     // MParContainer
-    //Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
 
-    ClassDef(MHDisp, 1) //3D-histogram in alpha, x and y
+    ClassDef(MHDisp, 2) //Class to provide a Disp map
 };