Changeset 7594 for trunk/MagicSoft/Mars/mhflux/MHPhi.cc

Timestamp:

03/13/06 15:51:11 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhflux/MHPhi.cc (modified) (10 diffs)

trunk/MagicSoft/Mars/mhflux/MHPhi.cc

@@ -33 +33 @@
 // http://www.astro.uni-wuerzburg.de/results/ringmethod/
 //
+// Class Version 2:
+//   + TH1D fHPhi;             // Phi plot of the signal w.r.t. source
+//   + TH1D fHPhiOff;          // Phi plot of the signal w.r.t. source+180deg
+//   + TH1D fHTemplate;        // Template how the background should look in the ideal case
+// 
+//   + TH1D fHInhom;           // Phi plot off the signal w.r.t. source (out of the ring with the signal)
+//   + TH1D fHInhomOff;        // Phi plot off the signal w.r.t. source+180deg (out of the ring with the signal)
+// 
+//   + Int_t   fNumBinsSignal; // Number of bins for signal region
+//   + Float_t fThetaCut;      // Theta cut
+//   + Float_t fDistSrc;       // Nominal distance of source from center in wobble
+// 
+//   + Bool_t  fUseAntiPhiCut; // Whether to use a anti-phi cut or not
+//
 ////////////////////////////////////////////////////////////////////////////
 #include "MHPhi.h"
@@ -63 +77 @@
 //
 MHPhi::MHPhi(const char *name, const char *title)
-: fHillas(0), fSrcPos(0), fDisp(0)//, fOnOffMode(kTRUE), fIsOffLoop(kFALSE)
+    : fHillas(0), fSrcPos(0), fDisp(0),
+    fNumBinsSignal(2), fThetaCut(0.23), fDistSrc(0.4),
+    fUseAntiPhiCut(kTRUE)
 {
     fName  = name  ? name  : "MHPhi";
@@ -70 +86 @@
     // Init Graphs
     fHPhi.SetNameTitle("Phi", "\\Delta\\Phi-Distribution");
-
-    fHPhi.SetXTitle("\\Delta\\Phi [\\circ]");
+    fHPhi.SetXTitle("\\Delta\\Phi' [\\circ]");
     fHPhi.SetYTitle("Counts");
 
@@ -79 +94 @@
     fHPhi.SetBit(TH1::kNoStats);
     fHPhi.SetMarkerStyle(kFullDotMedium);
-
-    fHPhi.GetYaxis()->SetTitleOffset(1.2);
-
-    /*
-    fNameParameter = "Disp";
-
-    fHist.SetNameTitle("Phi", "\\Delta\\Phi-Distribution");
-    fHist.SetZTitle("\\Delta\\Phi [\\circ]");
-    fHist.SetDirectory(NULL);
-
-    // Main histogram
-    fHistTime.SetName("Phi");
-    fHistTime.SetXTitle("\\Delta\\Phi [\\circ]");
-    fHistTime.SetDirectory(NULL);
-
-    MBinning binsa, binse, binst;
-    //binsa.SetEdges(75, 0, 1.5);
-    //binsa.SetEdges(arr);
-    binse.SetEdgesLog(15, 10, 100000);
-    binst.SetEdgesASin(67, -0.005, 0.665);
-    //binsa.Apply(fHistTime);
-
-    MH::SetBinning(&fHist, &binst, &binse, &binsa);
-     */
-}
-
-/*
-Double_t MHPhi::GetVal() const
-{
-    const Dopuble_t disp = static_cast<const MParameterD*>(fParameter)->GetVal();
-
-    const TVector2 pos = fHillas->GetMean()*fConvMm2Deg + fHillas->GetNormAxis()*disp;
-    const TVector2 src = fSrcPos->GetXY()*fConvMm2Deg;
-
-    // Calculate radial distance.
-    const Double_t d = pos.Mod() - src.Mod();
-
-    if (d<-fThetaCut*0.913 || d>fThetaCut)
-        return kTRUE;
-
-    const Double_t delta = src.DeltaPhi(pos)*TMath::RadToDeg();
-    const Double_t absd  = TMath::Abs(delta)
-
-    return fHistOff ? absd : 180-absd;
-}
-*/
+    fHPhi.SetLineColor(kBlue);
+    fHPhi.SetMarkerColor(kBlue);
+    fHPhi.GetYaxis()->SetTitleOffset(1.3);
+
+    fHPhiOff.SetMinimum(0);
+    fHPhiOff.SetDirectory(0);
+    fHPhiOff.Sumw2();
+    fHPhiOff.SetBit(TH1::kNoStats);
+    fHPhiOff.SetLineColor(kRed);
+    fHPhiOff.SetMarkerColor(kRed);
+
+    fHTemplate.SetMinimum(0);
+    fHTemplate.SetDirectory(0);
+    fHTemplate.SetBit(TH1::kNoStats);
+    fHTemplate.SetLineColor(kGreen);
+
+    fHInhom.SetNameTitle("Inhomogeneity", "\\Delta\\Phi-Distribution");
+    fHInhom.SetXTitle("\\Delta\\Phi' [\\circ]");
+    fHInhom.SetYTitle("Counts");
+    fHInhom.Sumw2();
+    fHInhom.SetMinimum(0);
+    fHInhom.SetDirectory(0);
+    fHInhom.SetBit(TH1::kNoStats);
+    fHInhom.GetYaxis()->SetTitleOffset(1.3);
+
+    fHInhomOff.Sumw2();
+    fHInhomOff.SetMinimum(0);
+    fHInhomOff.SetDirectory(0);
+    fHInhomOff.SetBit(TH1::kNoStats);
+    fHInhomOff.SetLineColor(kRed);
+    fHInhomOff.SetMarkerColor(kRed);
+}
 
 // --------------------------------------------------------------------------
@@ -161 +162 @@
     fConvMm2Deg = geom->GetConvMm2Deg();
 
-    fNumBinsSignal = 2;
-    fThetaCut      = 0.21/1.2;
-    fDistSrc       = 0.4;
-    //fIsOffLoop = !fIsOffLoop;
+    MParameterD *cut = (MParameterD*)plist->FindObject("ThetaSquaredCut", "MParameterD");
+    if (!cut)
+        *fLog << warn << "ThetaSquareCut [MParameterD] not found... using default theta<" << fThetaCut << "." << endl;
+    else
+        fThetaCut = TMath::Sqrt(cut->GetVal());
 
     const Double_t w  = TMath::ATan(fThetaCut/fDistSrc);
@@ -170 +172 @@
     const Int_t    n  = TMath::Nint(TMath::Ceil(180/sz));
 
-    MBinning(n, 0, n*sz).Apply(fHPhi);
-    /*
-
-    // Get Histogram binnings
-    MBinning binst, binse;
-    binst.SetEdges(fHist, 'x');
-    binse.SetEdges(fHist, 'y');
-
-    MBinning binsa(n, 0, n*sz);
-
-    // Apply binning
-    binsa.Apply(fHistTime);
-    MH::SetBinning(&fHist, &binst, &binse, &binsa);
-
-    // Remark: Binnings might be overwritten in MHAlpha::SetupFill
-    return MHAlpha::SetupFill(pl);
-     */
+    MBinning(n+3, 0, (n+3)*sz).Apply(fHPhi);
+    MBinning(n+3, 0, (n+3)*sz).Apply(fHPhiOff);
+    MBinning(n+3, 0, (n+3)*sz).Apply(fHTemplate);
+    MBinning(n+3, 0, (n+3)*sz).Apply(fHInhom);
+    MBinning(n+3, 0, (n+3)*sz).Apply(fHInhomOff);
+
     return kTRUE;
 }
@@ -197 +188 @@
 Bool_t MHPhi::Fill(const MParContainer *par, const Stat_t weight)
 {
+    // Here we calculate an upper phi cut to take a
+    // possible anti-theta cut into account
+    const Double_t ulim = fUseAntiPhiCut ? 180-fHPhi.GetBinLowEdge(fNumBinsSignal+1)*1.1 : 180;
+
+    // Calculate the shower origin and the source position in units of deg
     const TVector2 pos = fHillas->GetMean()*fConvMm2Deg + fHillas->GetNormAxis()*fDisp->GetVal();
     const TVector2 src = fSrcPos->GetXY()*fConvMm2Deg;
 
-    // Calculate radial distance.
+    // Calculate radial distance between shower origin and source
     const Double_t d = pos.Mod() - src.Mod();
 
-    if (d<-fThetaCut*0.913 || d>fThetaCut)
+    // define an upper and lower cut for the radial distance between both
+    const Double_t dR = fThetaCut;
+    const Double_t dr = fThetaCut*0.913;
+
+    // calculate the phi-angle of the shower origin w.r.t. the source position
+    const Double_t delta = src.DeltaPhi(pos)*TMath::RadToDeg();
+
+    // Define the radii of the upper and lower ring border
+    const Double_t R = src.Mod()+dR;
+    const Double_t r = src.Mod()-dr;
+
+    // Calculate a scale to scale all source positions to the
+    // nominal distance to center
+    const Double_t scale = src.Mod()/fDistSrc;
+
+    // Fill a phi-histograms with all events outside the ring
+    // Take the upper phi cut into account
+    if ((d<-dr || d>dR)/*TMath::Abs(d)>fThetaCut*1.2*/ && TMath::Abs(delta)<ulim)
+    {
+        fHInhom.Fill(TMath::Abs(delta)*scale,  weight);
+        fHInhomOff.Fill((ulim-TMath::Abs(delta))*scale,  weight);
+    }
+
+    // Now forget about all events which are not inside the ring
+    if (d<-dr || d>dR)
         return kTRUE;
 
-    const Double_t delta = src.DeltaPhi(pos)*TMath::RadToDeg();
-
-    fHPhi.Fill(TMath::Abs(delta), weight);
-
-    // const Double_t absd = TMath::Abs(delta)
-    // fHPhi.Fill(fHistOff ? absd : 180-absd, weight);
+    // Fill the histograms for on and off with the scaled phi
+    // only if we are below the upper phi cut
+    if (TMath::Abs(delta)<ulim)
+    {
+        fHPhi.Fill(         TMath::Abs(delta)*scale,  weight);
+        fHPhiOff.Fill((ulim-TMath::Abs(delta))*scale, weight);
+    }
+
+    // Calculate the maximum scaled phi taking the upper phi cut into account
+    const Double_t max = scale*ulim;
+
+    // Fill a template, this is how the phi-plot would look like
+    // without a signal and for an ideal camera.
+    const Int_t n = fHTemplate.GetNbinsX();
+    TArrayD arr(n);
+    for (int i=1; i<=n; i++)
+    {
+        const Double_t hi = fHTemplate.GetBinLowEdge(i+1);
+        const Double_t lo = fHTemplate.GetBinLowEdge(i);
+
+        // Decide whether the bin is fully contained in the upper phi-cut or
+        // the maximum possible phi is inside the bin
+        if (hi<max)
+            arr[i-1] = 1;
+        else
+            if (lo<max) // if its inside calculate the ratio
+                arr[i-1] = (max-lo)/fHTemplate.GetBinWidth(i+1);
+            else
+                break;
+    }
+
+    // The template is scaled with the current ring width. The upper phi-
+    // cut must not be taken into account because it is just a constant
+    // for all events.
+    const Double_t sum = arr.GetSum();
+    for (int i=1; i<=n; i++)
+        fHTemplate.AddBinContent(i, (R*R-r*r)*arr[i-1]/sum);
 
     return kTRUE;
@@ -227 +278 @@
     pad->SetBorderMode(0);
 
-    AppendPad("combine");
+    AppendPad("update");
+
+    pad->Divide(2,2);
+
+    // --------------------------
+
+    pad->cd(1);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
 
     fHPhi.Draw();
-
-    AppendPad(opt);
-}
-
-void MHPhi::Paint(Option_t *o)
-{
-    //TString opt(o);
-    //opt.ToLower();
-
-    // if (opt=="combine" && fHistOff)
-    // {
-    //    fHPhi.Add(fHist, fHistOff);
-    //    return;
-    // }
-
-    const Bool_t wobble = TString(o).Contains("anticut", TString::kIgnoreCase);
-
-    const Double_t cut  = fHPhi.GetBinLowEdge(fNumBinsSignal+1);
-
-    const Int_t maxbin  = wobble ? fHPhi.GetXaxis()->FindFixBin(180-cut)-1 : fHPhi.GetNbinsX();
-    const Double_t cut2 = wobble ? fHPhi.GetBinLowEdge(maxbin+1) : 180;
-
-    const Double_t sig  = fHPhi.Integral(1, fNumBinsSignal);
-    const Double_t bg   = fHPhi.Integral(1+fNumBinsSignal, maxbin);
-
-    const Double_t f    = cut/(cut2-cut);
-
-    const Double_t S0   = MMath::SignificanceLiMaSigned(sig, bg*f);
-    const Double_t S    = MMath::SignificanceLiMaSigned(sig, bg, f);
-
-    const TString  fmt  = Form("\\sigma_{L/M}=%.1f (\\sigma_{0}=%.1f)  \\Delta\\Phi_{on}<%.1f\\circ  \\Delta\\Phi_{off}<%.1f\\circ  E=%.0f  B=%.0f  f=%.2f",
-                               S, S0, cut, cut2, sig-bg*f, bg*f, f);
+    fHPhiOff.Draw("same");
+
+    TH1D *h1 = new TH1D(fHTemplate);
+    h1->SetName("Template");
+    h1->SetBit(kCanDelete);
+    h1->SetDirectory(0);
+    h1->Draw("same");
+
+    // --------------------------
+
+    pad->cd(2);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
+
+    fHInhom.Draw();
+    fHInhomOff.Draw("same");
+
+    TH1D *h2 = new TH1D(fHTemplate);
+    h2->SetName("Template");
+    h2->SetBit(kCanDelete);
+    h2->SetDirectory(0);
+    h2->Draw("same");
+
+    // --------------------------
+
+    pad->cd(3);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
+
+    fHPhi.Draw();
+    TH1D *h4 = new TH1D(fHInhom);
+    h4->SetName("Inhomogeneity");
+    h4->SetBit(kCanDelete);
+    h4->SetDirectory(0);
+    h4->Draw("same");
+
+    h1->Draw("same");
+
+    // --------------------------
+
+    pad->cd(4);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
+
+    TH1D *h3 = new TH1D(fHPhi);
+    h3->SetName("Result");
+    h3->SetBit(kCanDelete);
+    h3->SetDirectory(0);
+    h3->Draw();
+
+    h1->Draw("same");
+
+    // --------------------------
+
+    pad->cd();
+    AppendPad("result");
+}
+
+void MHPhi::PaintUpdate() const
+{
+    TVirtualPad *pad1 = gPad->GetPad(1);
+    TVirtualPad *pad2 = gPad->GetPad(2);
+    TVirtualPad *pad3 = gPad->GetPad(3);
+    TVirtualPad *pad4 = gPad->GetPad(4);
+
+    Double_t sig2 = 0;
+    Double_t bg2  = 0;
+    Double_t f2   = 1;
+    TH1D *htemp = pad1 ? dynamic_cast<TH1D*>(pad1->FindObject("Template")) : NULL;
+    if (htemp)
+    {
+        fHTemplate.Copy(*htemp);
+        htemp->SetName("Template");
+        htemp->SetDirectory(0);
+
+        Double_t sc1 = 1;
+        Double_t sc2 = 1;
+
+        TH1D *res = pad4 ? dynamic_cast<TH1D*>(pad4->FindObject("Result")) : NULL;
+        if (res)
+        {
+            fHPhi.Copy(*res);
+
+            // Scale inhomogeneity to match the phi-plot in the off-region
+            sc1 = res->Integral(fNumBinsSignal+1, 9999)/fHInhom.Integral(fNumBinsSignal+1, 9999);
+            // Scale inhomogeneity to match the phi-plot in the off-region
+            sc2 = fHInhom.Integral()/htemp->Integral();
+
+            res->Add(&fHInhom, -sc1);
+            res->Add(htemp,     sc1*sc2);
+            res->SetName("Result");
+            res->SetDirectory(0);
+
+            htemp->Scale(res->Integral(fNumBinsSignal+1, 9999)/htemp->Integral(fNumBinsSignal+1, 9999));
+
+            sig2 = res->Integral(1, fNumBinsSignal);
+            bg2  = fHPhi.Integral(fNumBinsSignal+1, 9999);
+            f2   = htemp->Integral(1, fNumBinsSignal)/htemp->Integral(fNumBinsSignal+1, 9999);
+        }
+
+        TH1D *hinhom = pad3 ? dynamic_cast<TH1D*>(pad3->FindObject("Inhomogeneity")) : NULL;
+        if (hinhom)
+        {
+            fHInhom.Copy(*hinhom);
+            hinhom->SetName("Inhomogeneity");
+            hinhom->SetDirectory(0);
+            hinhom->Scale(sc1);
+        }
+    }
+
+    htemp = pad2 ? dynamic_cast<TH1D*>(pad2->FindObject("Template")) : NULL;
+    if (htemp)
+    {
+        fHTemplate.Copy(*htemp);
+        htemp->Scale(fHInhom.Integral()/htemp->Integral());
+        htemp->SetName("Template");
+        htemp->SetDirectory(0);
+    }
+}
+
+void MHPhi::PaintText(const TH1D &res) const
+{
+    const Double_t cut  = res.GetBinLowEdge(fNumBinsSignal+1);
+
+    const Double_t sig  = res.Integral(1, fNumBinsSignal);
+    const Double_t bg   = res.Integral(fNumBinsSignal+1, 9999);
+
+    const Double_t f    = fHTemplate.Integral(1, fNumBinsSignal)/fHTemplate.Integral(fNumBinsSignal+1, 9999);
+
+    const Double_t S0   = MMath::SignificanceLiMaSigned(sig,  bg*f);
+    const Double_t S    = MMath::SignificanceLiMaSigned(sig,  bg, f);
+
+    const TString fmt = Form("\\sigma_{L/M}=%.1f (\\sigma_{0}=%.1f)  \\Delta\\Phi_{on}<%.1f\\circ  E=%.0f  B=%.0f  f=%.2f",
+                             S, S0, cut, sig-bg*f, bg*f, f);
 
     const Double_t b = bg             *f/fNumBinsSignal;
     const Double_t e = TMath::Sqrt(bg)*f/fNumBinsSignal;
 
-    TLatex text(0.27, 0.94, fmt);
+    TLatex text(0.275, 0.95, fmt);
     text.SetBit(TLatex::kTextNDC);
-    text.SetTextSize(0.035);
+    text.SetTextSize(0.042);
     text.Paint();
 
@@ -274 +435 @@
     line.SetLineColor(14);
     line.PaintLine(cut, gPad->GetUymin(), cut, gPad->GetUymax());
-    if (maxbin<fHPhi.GetNbinsX())
-        line.PaintLine(cut2, gPad->GetUymin(), cut2, gPad->GetUymax());
-    line.SetLineColor(kBlue);
+
+    // Here we calculate an upper phi cut to take a
+    // possible anti-theta cut into account
+    const Double_t ulim = fUseAntiPhiCut ? 180-fHPhi.GetBinLowEdge(fNumBinsSignal+1)*1.1 : 180;
+    line.SetLineStyle(kDotted);
+    line.PaintLine(ulim, gPad->GetUymin(), ulim, gPad->GetUymax());
+
+    line.SetLineStyle(kSolid);
+    line.SetLineColor(kBlack);
     line.PaintLine(0, b, cut, b);
     line.PaintLine(cut/2, b-e, cut/2, b+e);
@@ -283 +450 @@
 
     TMarker m;
-    m.SetMarkerColor(kBlue);
     m.SetMarkerStyle(kFullDotMedium);
     m.PaintMarker(cut/2, b);
 }
+
+void MHPhi::PaintResult() const
+{
+    TVirtualPad *pad = gPad;
+
+    pad->cd(1);
+    PaintText(fHPhi);
+
+    pad->cd(4);
+    TH1D *res = gPad ? dynamic_cast<TH1D*>(gPad->FindObject("Result")) : NULL;
+    if (res)
+        PaintText(*res);
+}
+
+void MHPhi::Paint(Option_t *o)
+{
+    TString opt(o);
+    if (opt=="update")
+        PaintUpdate();
+    if (opt=="result")
+        PaintResult();
+}
+
+Int_t MHPhi::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "NumBinsSignal", print))
+    {
+        SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal));
+        rc = kTRUE;
+    }
+    if (IsEnvDefined(env, prefix, "UseAntiPhiCut", print))
+    {
+        SetUseAntiPhiCut(GetEnvValue(env, prefix, "UseAntiPhiCut", (Int_t)fUseAntiPhiCut));
+        rc = kTRUE;
+    }
+
+    return rc;
+}