Changeset 7217 for trunk/MagicSoft/Mars/mhflux

Timestamp:

07/25/05 10:35:31 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• MHDisp.cc (modified) (17 diffs)
• MHDisp.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mhflux/MHDisp.cc

@@ -51 +51 @@
 #include "MParameters.h"
 
-#include "MHillasExt.h"
-#include "MHillasSrc.h"
+#include "MHillas.h"
 #include "MSrcPosCam.h"
 #include "MPointingPos.h"
+#include "MPointingDev.h"
 
 ClassImp(MHDisp);
@@ -65 +65 @@
 //
 MHDisp::MHDisp(const char *name, const char *title)
-    : fHilExt(0), fDisp(0), fSmearing(-1), fWobble(kFALSE)
+    : fDisp(0), fDeviation(0), fSmearing(-1), fWobble(kFALSE)
 {
     //
@@ -106 +106 @@
     }
 
-    fHilExt = (MHillasExt*)plist->FindObject("MHillasExt");
-    if (!fHilExt)
-    {
-        *fLog << err << "MHillasExt not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fHilSrc = (MHillasSrc*)plist->FindObject("MHillasSrc");
-    if (!fHilSrc)
-    {
-        *fLog << err << "MHillasSrc not found... aborting." << endl;
-        return kFALSE;
-    }
+    fDeviation = (MPointingDev*)plist->FindObject("MPointingDev");
+    if (!fDeviation)
+        *fLog << warn << "MPointingDev not found... ignored." << endl;
 
     MBinning binsx, binsy;
@@ -160 +150 @@
     Double_t rho = 0;
     if (fTime && fObservatory && fPointPos)
-        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
+        rho = fPointPos->RotationAngle(*fObservatory, *fTime, fDeviation);
 
     // Vector of length disp in direction of shower
@@ -347 +337 @@
 //
 // Return the mean signal in h around (0,0) in a distance between
-// 0.33 and 0.47deg
+// 0.325 and 0.475deg
 //
 Double_t MHDisp::GetOffSignal(TH1 &h) const
@@ -360 +350 @@
         {
             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)
+            if (d>0.325 && d<0.475)
             {
                 sum += h.GetBinContent(x+1,y+1);
@@ -367 +356 @@
             }
         }
-
     return sum/cnt;
 }
@@ -415 +403 @@
             if (rmax<TMath::Hypot(px, py))
                 h2.SetBinContent(bin, 0);
+            else
+                if (h2.GetBinContent(bin)==0)
+                    h2.SetBinContent(bin, 1e-10);
         }
 }
@@ -427 +418 @@
     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));
-
+    const Int_t n = TMath::Nint(0.2/axex.GetBinWidth(1));
+
+    for (int x=1+n; x<=axex.GetNbins()-n; x++)
+        for (int y=1+n; y<=axey.GetNbins()-n; y++)
+        {
+            Double_t sig=0;
+            Double_t bg =0;
+
+            for (int dx=-n; dx<=n; dx++)
+                for (int dy=-n; dy<=n; dy++)
+                {
+                    if (TMath::Hypot((float)dx, (float)dy)>n)
+                        continue;
+
+                    const Int_t  bin = s.GetBin(x+dx,y+dy);
+
+                    sig += h1.GetBinContent(bin);
+                    bg  += h2.GetBinContent(bin);
+                }
+
+            const Double_t S = sig>0 ? MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale)) : 0;
+
+            const Int_t bin = s.GetBin(x,y);
             s.SetBinContent(bin, S);
         }
 }
 
+void MHDisp::Profile1D(const char *name, const TH2 &h) const
+{
+    if (!gPad)
+        return;
+
+    TH1D *hrc = dynamic_cast<TH1D*>(gPad->FindObject(name));
+    if (!hrc)
+        return;
+
+    hrc->Reset();
+
+    const Double_t max = h.GetMaximum();
+
+    MBinning(50, -max*1.1, max*1.1).Apply(*hrc);
+
+    for (int x=1; x<=h.GetXaxis()->GetNbins(); x++)
+        for (int y=1; y<=h.GetXaxis()->GetNbins(); y++)
+        {
+            const Int_t   bin  = h.GetBin(x,y);
+            const Double_t sig = h.GetBinContent(bin);
+            if (sig!=0)
+                hrc->Fill(sig);
+        }
+
+    gPad->SetLogy(hrc->GetMaximum()>0);
+
+    if (!fHistOff || hrc->GetRMS()<0.1)
+        return;
+
+    // ---------- Fix mean ----------
+    // TF1 *g = (TF1*)gROOT->GetFunction("gaus");
+    // g->FixParameter(1, 0);
+    // hrc->Fit("gaus", "BQI");
+
+    hrc->Fit("gaus", "QI");
+
+    TF1 *f = hrc->GetFunction("gaus");
+    if (f)
+    {
+        f->SetLineWidth(1);
+        f->SetLineColor(kBlue);
+    }
+}
 
 void MHDisp::Paint(Option_t *o)
@@ -460 +507 @@
     h1->SetContour(99);
 
+    TH2 *hx=0;
+
     Double_t scale = 1;
     if (fHistOff)
@@ -467 +516 @@
         TH2 *h=(TH2*)fHistOff->Project3D("yx_off");
 
-        scale = -1;
-
         const Double_t h1off = GetOffSignal(*h1);
         const Double_t hoff  = GetOffSignal(fHistBg);
@@ -474 +521 @@
         scale = hoff==0 ? -1 : -h1off/hoff;
 
-        const Bool_t sig = kFALSE;
-
-        if (!sig)
-            h1->Add(&fHistBg, scale);
-        else
-            MakeSignificance(*h1, *h1, fHistBg, scale);
+        hx = (TH2*)pad->GetPad(4)->FindObject("Alpha_yx_sig");
+        if (hx)
+        {
+            hx->SetContour(99);
+            MakeSignificance(*hx, *h1, fHistBg, scale);
+            MakeDot(*hx);
+            MakeSymmetric(hx);
+        }
+
+        h1->Add(&fHistBg, scale);
 
         MakeDot(*h1);
+        MakeSymmetric(h1);
 
         delete h;
-
-        MakeSymmetric(h1);
     }
 
@@ -515 +565 @@
         func.SetParameter(0, h2->GetBinContent(1));
         func.FixParameter(1, 0);
-        func.SetParameter(2, 0.15);
+        func.SetParameter(2, 0.12);
         if (fHistOff)
             func.FixParameter(3, 0);
-        func.SetParameter(4, h2->GetBinContent(15));
+        func.SetParameter(4, 0);//h2->GetBinContent(20));
         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)));
     }
+
+    pad->cd(5);
+    if (h1)
+        Profile1D("Exc", *h1);
+
+    pad->cd(6);
+    if (hx)
+        Profile1D("Sig", *hx);
 }
 
@@ -533 +591 @@
     AppendPad(o);
 
-    TString name = Form("%s_1", pad->GetName());
-    TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
+    // ----- Pad number 4 -----
+    TString name = Form("%s_4", pad->GetName());
+    TPad *p = new TPad(name,name, 0.5025, 0.3355, 0.995, 0.995, col, 0, 0);
+    p->SetNumber(4);
+    p->Draw();
+    p->cd();
+
+    TH1 *h3 = fHist.Project3D("yx_sig");
+    h3->SetTitle("Significance Map");
+    h3->SetDirectory(NULL);
+    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
+    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
+    h3->SetMinimum(0);
+    h3->Draw("colz");
+    h3->SetBit(kCanDelete);
+
+    if (fHistOff)
+        GetCatalog()->Draw("white mirror same *");
+
+    // ----- Pad number 1 -----
+    pad->cd();
+    name = Form("%s_1", pad->GetName());
+    p = new TPad(name,name, 0.005, 0.3355, 0.4975, 0.995, col, 0, 0);
     p->SetNumber(1);
     p->Draw();
     p->cd();
 
-    TH1 *h3 = fHist.Project3D("yx_on");
+    h3 = fHist.Project3D("yx_on");
     h3->SetTitle("Distribution of equivalent events");
     h3->SetDirectory(NULL);
@@ -551 +630 @@
         GetCatalog()->Draw("white mirror same *");
 
+    // ----- Pad number 2 -----
     pad->cd();
     name = Form("%s_2", pad->GetName());
-    p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
+    p = new TPad(name,name, 0.005, 0.005, 0.2485, 0.3315, col, 0, 0);
     p->SetNumber(2);
     p->Draw();
@@ -559 +639 @@
     h3->Draw("surf3");
 
+    // ----- Pad number 3 -----
     pad->cd();
     name = Form("%s_3", pad->GetName());
-    p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
+    p = new TPad(name,name, 0.2525, 0.005, 0.4985, 0.3315, col, 0, 0);
     p->SetNumber(3);
     p->SetGrid();
@@ -577 +658 @@
     h->SetBit(kCanDelete);
     h->Draw();
+
+    // ----- Pad number 5 -----
+    pad->cd();
+    name = Form("%s_5", pad->GetName());
+    p = new TPad(name,name, 0.5025, 0.005, 0.7485, 0.3315, col, 0, 0);
+    p->SetNumber(5);
+    p->SetGrid();
+    p->Draw();
+    p->cd();
+
+    h3 = new TH1D("Exc", "Distribution of excess events", 10, -1, 1);
+    h3->SetDirectory(NULL);
+    h3->SetXTitle("N");
+    h3->SetYTitle("Counts");
+    h3->Draw();
+    h3->SetBit(kCanDelete);
+
+    // ----- Pad number 6 -----
+    pad->cd();
+    name = Form("%s_6", pad->GetName());
+    p = new TPad(name,name, 0.7525, 0.005, 0.9995, 0.3315, col, 0, 0);
+    p->SetNumber(6);
+    p->SetGrid();
+    p->Draw();
+    p->cd();
+
+    h3 = new TH1D("Sig", "Distribution of significances", 10, -1, 1);
+    h3->SetDirectory(NULL);
+    h3->SetXTitle("N");
+    h3->SetYTitle("Counts");
+    h3->Draw();
+    h3->SetBit(kCanDelete);
 }
 

trunk/MagicSoft/Mars/mhflux/MHDisp.h

@@ -16 +16 @@
 class MHillasSrc;
 class MParameterD;
+class MPointingDev;
 
 class MHDisp : public MHFalseSource
 {
 private:
-    MHillasExt   *fHilExt;  //!
-    MHillasSrc   *fHilSrc;  //!
     MParameterD  *fDisp;    //!
+    MPointingDev *fDeviation; //!
 
     TH2D         fHistBg;
-
     TH2D         fHistBg1;
     TH2D         fHistBg2;
@@ -39 +38 @@
     Double_t DeltaPhiSrc(const TVector2 &v) const;
 
+    void Profile1D(const char *name, const TH2 &h) 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;