Changeset 5901 for trunk

Timestamp:

01/19/05 13:17:09 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (3 diffs)
• mbase/MMath.cc (modified) (2 diffs)
• mhflux/MAlphaFitter.cc (modified) (10 diffs)
• mhflux/MAlphaFitter.h (modified) (5 diffs)
• mhflux/MHAlpha.cc (modified) (6 diffs)
• mhflux/MHFalseSource.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -21 +21 @@
                                                  -*-*- END OF LINE -*-*-
 
- 2205/01/19 Abelardo Moralejo
+ 2005/01/19 Thomas Bretz
+ 
+   * mbase/MMath.cc:
+     - added a comment to SignificanceLiMa, made by Daniel Mazin
+     - also check for b==0
+
+   * mhflux/MAlphaFitter.[h,cc]:
+     - fixed significance calculation in case of on-off data
+     - added fScaleFactor
+
+   * mhflux/MHAlpha.[h,cc], mhflux/MHFalseSource.cc:
+     - handle scale factor in case of on-off observations
+
+
+
+ 2005/01/19 Abelardo Moralejo
 
    * mtrigger/MFTriggerPattern.[cc,h]
@@ -27 +42 @@
        defined by N. Galante (see below)
 
- 2205/01/19 Nicola Galante
+
+
+ 2005/01/19 Nicola Galante
 
    * mtrigger/MFTriggerPattern.[cc,h]
@@ -35 +52 @@
      - fixed a bug in Process, both fMaskRequiredUnprescaled and
        fMaskRequiredPrescaled are checked simultaneously.
+
 
 

trunk/MagicSoft/Mars/mbase/MMath.cc

@@ -70 +70 @@
 //
 //  Returns -1 if sum<0 or alpha<0 or the argument of sqrt<0
-//  Returns  0 if s+b==0 or s==0
+//  Returns  0 if s+b==0, s==0 or b==0
+//
+// Here is some eMail written by Daniel Mazin about the meaning of the arguments:
+//
+//  > Ok. Here is my understanding:
+//  > According to Li&Ma paper (correctly cited in MMath.cc) alpha is the
+//  > scaling factor. The mathematics behind the formula 17 (and/or 9) implies
+//  > exactly this. If you scale OFF to ON first (using time or using any other
+//  > method), then you cannot use formula 17 (9) anymore. You can just try
+//  > the formula before scaling (alpha!=1) and after scaling (alpha=1), you
+//  > will see the result will be different.
+//
+//  > Here are less mathematical arguments:
+//
+//  >  1) the better background determination you have (smaller alpha) the more
+//  > significant is your excess, thus your analysis is more sensitive. If you
+//  > normalize OFF to ON first, you loose this sensitivity.
+//
+//  >  2) the normalization OFF to ON has an error, which naturally depends on
+//  > the OFF and ON. This error is propagating to the significance of your
+//  > excess if you use the Li&Ma formula 17 correctly. But if you normalize
+//  > first and use then alpha=1, the error gets lost completely, you loose
+//  > somehow the criteria of goodness of the normalization.
 //
 Double_t MMath::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha)
@@ -76 +98 @@
     const Double_t sum = s+b;
 
-    if (s==0 || sum==0)
+    if (s==0 || b==0 || sum==0)
         return 0;
 

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

@@ -65 +65 @@
     fChiSqBg=0;
     fIntegralMax=0;
+    fScaleFactor=1;
 
     fCoefficients.Reset();
@@ -210 +211 @@
 }
 
-Bool_t MAlphaFitter::Fit(TH1D &hon, TH1D &hof, Bool_t paint)
+Bool_t MAlphaFitter::Fit(TH1D &hon, TH1D &hof, Double_t alpha, Bool_t paint)
 {
     /*
@@ -224 +225 @@
     fit.SetPolynomOrder(1);
 
-    if (!fit.Fit(h, paint))
+    if (alpha<=0 || !fit.Fit(h, paint))
         return kFALSE;
 
@@ -237 +238 @@
     fEventsSignal     = hon.Integral(0, bin);
     fEventsExcess     = fEventsSignal-fEventsBackground;
-    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
+    fScaleFactor      = alpha;
+    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground/alpha, alpha);
 
     if (fEventsExcess<0)
@@ -399 +401 @@
     h0->SetDirectory(0);
 
-    Scale(*h0, *h1);
-
-    const Bool_t rc = Fit(*h1, *h0, paint);
+    const Bool_t rc = ScaleAndFit(*h1, h0, paint);
 
     delete h0;
@@ -419 +419 @@
     h0->SetDirectory(0);
 
-    Scale(*h0, *h1);
-
-    const Bool_t rc = Fit(*h1, *h0, paint);
+    const Bool_t rc = ScaleAndFit(*h1, h0, paint);
 
     delete h0;
@@ -439 +437 @@
     h0->SetDirectory(0);
 
-    Scale(*h0, *h1);
-
-    const Bool_t rc = Fit(*h1, *h0, paint);
+    const Bool_t rc = ScaleAndFit(*h1, h0, paint);
 
     delete h0;
@@ -449 +445 @@
 }
 
-void MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
+Double_t MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
 {
     Float_t scaleon = 1;
@@ -456 +452 @@
     {
     case kNone:
-        return;
+        return 1;
 
     case kEntries:
@@ -481 +477 @@
         break;
 
+        // This is just to make some compiler happy
     default:
-        return;
+        return 1;
     }
 
     if (scaleof!=0)
+    {
         of.Scale(scaleon/scaleof);
+        return scaleon/scaleof;
+    }
     else
+    {
         of.Reset();
-}
+        return 0;
+    }
+}

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h

@@ -40 +40 @@
     Double_t fEventsExcess;
     Double_t fEventsSignal;
-    Double_t fEventsBackground;
+    Double_t fEventsBackground; // fScaleFactor already applied
 
     Double_t fChiSqSignal;
     Double_t fChiSqBg;
     Double_t fIntegralMax;
+    Double_t fScaleFactor;  // Scale factor for off-data
 
     TArrayD fCoefficients;
@@ -62 +63 @@
         gROOT->GetListOfFunctions()->Remove(fFunc);
         fFunc->SetName("Dummy");
+
+        Clear();
     }
 
@@ -103 +106 @@
     Double_t GetChiSqSignal() const        { return fChiSqSignal; }
     Double_t GetChiSqBg() const            { return fChiSqBg; }
+    Double_t GetScaleFactor() const        { return fScaleFactor; }
 
     Double_t GetGausSigma() const          { return fCoefficients[2]; }
@@ -113 +117 @@
 
     Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
-    Bool_t Fit(TH1D &on, TH1D &off, Bool_t paint=kFALSE);
+    Bool_t Fit(TH1D &on, TH1D &off, Double_t alpha, Bool_t paint=kFALSE);
+    Bool_t Fit(TH1D &on, TH1D *off, Double_t alpha, Bool_t paint=kFALSE)
+    {
+        return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
+    }
     Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
     {
-        return off ? Fit(on, *off, paint) : Fit(on, paint);
+        return off ? Fit(on, *off, 1, paint) : Fit(on, paint);
+    }
+    Bool_t ScaleAndFit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
+    {
+        const Double_t alpha = off ? Scale(*off, on) : 1;
+        return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
     }
 
@@ -140 +153 @@
     }
 
-    void Scale(TH1D &off, const TH1D &on) const;
+    Double_t Scale(TH1D &off, const TH1D &on) const;
 
     ClassDef(MAlphaFitter, 1)

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -392 +392 @@
 
     TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
-    const Bool_t rc = fit.Fit(fHAlphaTime, h);
+    const Bool_t rc = fit.ScaleAndFit(fHAlphaTime, h);
     if (h)
         delete h;
@@ -544 +544 @@
         {
             TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
-            fFit.Scale(*hoff, *hon);
+            const Double_t alpha = fFit.Scale(*hoff, *hon);
 
             hon->SetMaximum();
             hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
+
+            // BE CARFEULL: This is a really weird workaround!
+            hoff->SetMaximum(alpha);
 
             if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
@@ -564 +567 @@
         if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
         {
-            // Do not store the 'final' result in fFit
-            MAlphaFitter fit(fFit);
+            // Check whether we are dealing with on-off analysis
             TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff");
-            fit.Fit(*h0, hoff, kTRUE);
-            fit.PaintResult();
+            if (hoff)
+            {
+                // Do not store the 'final' result in fFit
+                MAlphaFitter fit(fFit);
+
+                // BE CARFEULL: This is a really weird workaround!
+                const Double_t alpha = hoff->GetMaximum();
+                fit.Fit(*h0, hoff, alpha<0 ? 1: alpha, kTRUE);
+
+                fit.PaintResult();
+            }
         }
 
@@ -757 +768 @@
 
         TH1D *hof = 0;
+        Double_t alpha = 1;
 
         if (fOffData)
@@ -772 +784 @@
             hof->Draw("same");
 
-            fit.Scale(*hof, *hon);
+            alpha = fit.Scale(*hof, *hon);
 
             hon->SetMaximum();
@@ -796 +808 @@
         }
 
-        if (hof ? fit.Fit(*hon, *hof) : fit.Fit(*hon))
-                *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
+        if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
+            *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
         /*
         if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))

trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc

@@ -954 +954 @@
 
     TH1D *h=0, *hoff=0;
+    Double_t scale = 1;
     for (int ix=1; ix<=nx; ix++)
         for (int iy=1; iy<=ny; iy++)
@@ -969 +970 @@
             {
                 hoff = fHistOff->ProjectionZ("AlphaFitOff", ix, ix, iy, iy);
-                hoff->Scale(entries/hoff->GetEntries());
+                hoff->Scale(entries/h->GetEntries());
+                scale = fit.Scale(*hoff, *h);
             }
 
-            if (!fit.Fit(*h, hoff))
+            if (!fit.Fit(*h, hoff, scale))
                 continue;