Changeset 8883

Timestamp:

03/18/08 17:47:37 (17 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhflux/MMcSpectrumWeight.cc (modified) (2 diffs)
• mhflux/MMcSpectrumWeight.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -45 +45 @@
    * mpointing/MPointingDevCalc.cc:
      - added some more comments
+
+   * mhflux/MMcSpectrumWeight.[h,cc]:
+     - added the possibility to extend a spectrum by "fake" statistics
 
 

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

@@ -428 +428 @@
 // and fEnergyMax correctly describe the range.
 //
+// If scale is given the histogram statistics is further extended by the
+// new spectrum according to the scale factor (eg. 1.2: by 20%)
+//
 // In the 1D case it is assumed that the x-axis is a zenith angle binning.
 // In the 2D case the x-axis is assumed to be zenith angle, the y-axis
 // to be energy.
 //
-void MMcSpectrumWeight::CompleteEnergySpectrum(TH1 &h, Double_t emin) const
+void MMcSpectrumWeight::CompleteEnergySpectrum(TH1 &h, Double_t emin, Double_t scale) const
 {
     if (h.InheritsFrom(TH3::Class()))
     {
-        *fLog << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum doesn't support TH3." << endl;
         return;
     }
 
-    if (fEnergyMin <= emin)
+    if (fEnergyMin < emin)
+    {
+        *fLog << err << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum: fEnergyMin (";
+        *fLog << fEnergyMin << ") smaller than emin (" << emin << ")." << endl;
         return;
-
+    }
+
+    // Total number of events for the new spectrum in the same
+    // energy range as the current histogram is filled
     const Double_t norm = GetSpecNewIntegral();
 
+    // Check if it is only a histogram in ZA
     if (!h.InheritsFrom(TH2::Class()))
     {
-        h.Scale(GetSpecNewIntegral(emin, fEnergyMax)/norm);
+        // Warning: Simply scaling the zenith angle distribution might
+        // increase fluctuations for low statistics.
+        const Double_t f = GetSpecNewIntegral(emin, fEnergyMax)/norm;
+        h.Scale(f*scale);
         return;
     }
@@ -453 +465 @@
     const TAxis &axey = *h.GetYaxis();
 
+    // Find energy range between the minimum energy to be filles (emin)
+    // and the minimum energy corresponding to the data filled into
+    // this histogram (fEnergyMin)
     const Int_t first = axey.FindFixBin(emin);
     const Int_t last  = axey.FindFixBin(fEnergyMin);
+    const Int_t max   = axey.FindFixBin(fEnergyMax);
 
     for (int x=1; x<=h.GetNbinsX(); x++)
     {
-         const Double_t f = h.Integral(x, x, -1, 9999)/norm;
-
-        for (int y=first; y<=last; y++)
-        {
-            const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
-            const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMin ? fEnergyMin : axey.GetBinLowEdge(y+1);
-
-            h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi));
-        }
+        // Ratio between the number of events in the zenith angle
+        // bin corresponding to x and the new spectrum.
+        const Double_t f = h.Integral(x, x, -1, 9999)/norm;
+
+        // Fill histogram with the "new spectrum" between
+        // emin and fEnergyMin.
+        if (emin<fEnergyMin)
+            for (int y=first; y<=last; y++)
+            {
+                // Check if the bin is only partly filled by the energy range
+                const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
+                const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMin ? fEnergyMin : axey.GetBinLowEdge(y+1);
+
+                // Add the new spectrum extending the existing spectrum
+                h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi));
+            }
+
+        // If scale is >1 we also have to increse the statistics f the
+        // histogram according to scale.
+        if (scale>1)
+            for (int y=first; y<=max; y++)
+            {
+                // Check if the bin is only partly filled by the energy range
+                const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
+                const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMax ? fEnergyMax : axey.GetBinLowEdge(y+1);
+
+                // Use the analytical solution to scale the histogram
+                h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi)*(scale-1));
+            }
     }
 }

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h

@@ -101 +101 @@
 
     // Functions
-    void CompleteEnergySpectrum(TH1 &h, Double_t emin) const;
+    void CompleteEnergySpectrum(TH1 &h, Double_t emin, Double_t scale=0) const;
 
     // TObject