Changeset 8680 for trunk/MagicSoft/Mars/mhflux

Timestamp:

08/20/07 11:04:58 (17 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• MMcSpectrumWeight.cc (modified) (18 diffs)
• MMcSpectrumWeight.h (modified) (3 diffs)

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

@@ -19 +19 @@
 !   Author(s): Marcos Lopez 10/2003 <mailto:marcos@gae.ucm.es>
 !
-!   Copyright: MAGIC Software Development, 2000-2006
+!   Copyright: MAGIC Software Development, 2000-2007
 !
 !
@@ -40 +40 @@
 //  Method:
 //  -------
-//                                 
+//
 //   - Corsika spectrun: dN/dE = A * E^(a)
 //     with a = fOldSlope, and A = N/integral{E*de} from ELowLim to EUppLim
@@ -72 +72 @@
 //    MMcEvt
 //    MMcCorsikaRunHeader
+//    [MPointingPos]
+//    [MHillas]
 //
 //  Output Container:
@@ -81 +83 @@
 #include <TF1.h>
 #include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
 #include <TSpline.h>
 
@@ -109 +113 @@
     fNameMcEvt   = "MMcEvt";
 
-    fNewSlope    = -9;
-    fOldSlope    = -9;
+    fNewSlope    = -99;
+    fOldSlope    = -99;
 
     fEnergyMin   = -1;
@@ -210 +214 @@
 // The slope is returned as %.3f
 //
-TString MMcSpectrumWeight::GetFormulaSpecOld() const
-{
-    return Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fOldSlope);
+TString MMcSpectrumWeight::GetFormulaSpecOld(const char *name) const
+{
+    return Form("pow(%s.fEnergy, %.3f)", name, fOldSlope);
 }
 
@@ -225 +229 @@
 // unchanged.
 //
-TString MMcSpectrumWeight::GetFormulaSpecNew() const
-{
-    TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula.Data();
+TString MMcSpectrumWeight::GetFormulaSpecNew(const char *name) const
+{
+    TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", name, fNewSlope) : fFormula.Data();
     if (!fFormula.IsNull())
-        str.ReplaceAll("X", Form("(%s.fEnergy)", fNameMcEvt.Data()));
+        str.ReplaceAll("X", Form("(%s.fEnergy)", name));
 
     return str;
@@ -258 +262 @@
 // Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600))
 //
-TString MMcSpectrumWeight::GetFormulaWeights() const
+TString MMcSpectrumWeight::GetFormulaWeights(const char *name) const
 {
     if (GetFormulaSpecOld()==GetFormulaSpecNew())
@@ -268 +272 @@
     const Double_t norm = fNorm*iold/inew;
 
-    return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data());
+    return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew(name).Data(), GetFormulaSpecOld(name).Data());
 }
 
@@ -276 +280 @@
 // GetFormulaSpecNewX() describing the destination spectrum
 //
-Double_t MMcSpectrumWeight::GetSpecNewIntegral() const
+Double_t MMcSpectrumWeight::GetSpecNewIntegral(Double_t emin, Double_t emax) const
 {
     TF1 funcnew("Dummy", GetFormulaSpecNewX());
-    return funcnew.Integral(fEnergyMin, fEnergyMax);
+    return funcnew.Integral(emin, emax);
 }
 
@@ -287 +291 @@
 // GetFormulaSpecOldX() describing the simulated spectrum
 //
-Double_t MMcSpectrumWeight::GetSpecOldIntegral() const
+Double_t MMcSpectrumWeight::GetSpecOldIntegral(Double_t emin, Double_t emax) const
 {
     TF1 funcold("Dummy", GetFormulaSpecOldX());
-    return funcold.Integral(fEnergyMin, fEnergyMax);
+    return funcold.Integral(emin, emax);
 }
 
@@ -357 +361 @@
     if (fEnergyMax>fEnergyMin && !fAllowChange)
     {
+        if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10)
+        {
+            *fLog << err;
+            *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change (";
+            *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl;
+            return kFALSE;
+        }
+
         if (TMath::Abs(fOldSlope-rh.GetSlopeSpec())>1e-10)
         {
@@ -364 +376 @@
             return kFALSE;
         }
-        if (TMath::Abs(fEnergyMin-rh.GetELowLim())>1e-10)
-        {
-            *fLog << err;
-            *fLog << "ERROR - The minimum simulated Monte Carlo energy is not allowed to change ";
-            *fLog << "(" << fEnergyMin << " --> " << rh.GetELowLim() << ")... abort." << endl;
-            return kFALSE;
-        }
-        if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10)
-        {
-            *fLog << err;
-            *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change (";
-            *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl;
-            return kFALSE;
-        }
-        return kTRUE;
+
+        // No change happened
+        if (TMath::Abs(fEnergyMin-rh.GetELowLim())<=1e-10)
+            return kTRUE;
+
+        // The lower energy limit has changed
+        *fLog << warn;
+        *fLog << "The minimum simulated Monte Carlo energy has changed from ";
+        *fLog << fEnergyMin << "GeV to " << rh.GetELowLim() << "GeV." << endl;
+        fEnergyMin = rh.GetELowLim();
     }
 
@@ -385 +392 @@
     fEnergyMax = rh.GetEUppLim();
 
-    if (fNewSlope==-9)
-    {
-        *fLog << inf << "The new slope of the power law is undefined (-9)... no weighting applied." << endl;
+    if (fNewSlope==-99)
+    {
+        *fLog << inf << "A new slope for the power law has not been defined... no weighting applied." << endl;
         fNewSlope = fOldSlope;
     }
@@ -403 +410 @@
 // ---------------------------------------------------------------------------
 //
+// completes a simulated spectrum starting at an energy fEnergyMin down to
+// an energy emin.
+//
+// It is assumed that the contents of MMcSpectrumWeight for the new spectrum
+// correctly describe the spectrum within the histogram, and fEnergyMin
+// and fEnergyMax correctly describe the range.
+//
+// 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
+{
+    if (h.InheritsFrom(TH3::Class()))
+    {
+        *fLog << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum doesn't support TH3." << endl;
+        return;
+    }
+
+    if (fEnergyMin <= emin)
+        return;
+
+    const Double_t norm = GetSpecNewIntegral();
+
+    if (!h.InheritsFrom(TH2::Class()))
+    {
+        h.Scale(GetSpecNewIntegral(emin, fEnergyMax)/norm);
+        return;
+    }
+
+    const TAxis &axey = *h.GetYaxis();
+
+    const Int_t first = axey.FindFixBin(emin);
+    const Int_t last  = axey.FindFixBin(fEnergyMin);
+
+    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));
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+//
 // The current contants are printed
 //
 void MMcSpectrumWeight::Print(Option_t *o) const
 {
+    const TString opt(o);
+
+    const Bool_t hasnew = opt.Contains("new") || opt.IsNull();
+    const Bool_t hasold = opt.Contains("old") || opt.IsNull();
+
     *fLog << all << GetDescriptor() << endl;
-    *fLog << " Simulated energy range:   " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
-    *fLog << " Simulated spectral slope: " << fOldSlope << endl;
-    *fLog << " New spectral slope:       " << fNewSlope << endl;
+
+    if (hasold)
+    {
+        *fLog << " Simulated energy range:   " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
+        *fLog << " Simulated spectral slope: ";
+        if (fOldSlope==-99)
+            *fLog << "undefined" << endl;
+        else
+            *fLog << fOldSlope << endl;
+    }
+    if (hasnew)
+    {
+        *fLog << " New spectral slope:       ";
+        if (fNewSlope==-99)
+            *fLog << "undefined" << endl;
+        else
+            *fLog << fNewSlope << endl;
+    }
     *fLog << " Additional user norm.:    " << fNorm << endl;
     *fLog << " Spectra are normalized:   " << (fNormEnergy<0?"by integral":Form("at %.1fGeV", fNormEnergy)) << endl;
-    *fLog << " Old Spectrum:     " << GetFormulaSpecOldX() << "   (I=" << GetSpecOldIntegral() << ")" << endl;
-    *fLog << " New Spectrum:     " << GetFormulaSpecNewX() << "   (I=" << GetSpecNewIntegral() << ")" << endl;
+    if (hasold)
+    {
+        *fLog << " Old Spectrum: " << GetFormulaSpecOldX();
+        if (fEnergyMin>=0 && fEnergyMax>0)
+            *fLog << "   (I=" << GetSpecOldIntegral() << ")";
+        *fLog << endl;
+    }
+    if (hasnew)
+    {
+        *fLog << " New Spectrum: " << GetFormulaSpecNewX();
+        if (fEnergyMin>=0 && fEnergyMax>0)
+            *fLog << "   (I=" << GetSpecNewIntegral() << ")";
+        *fLog << endl;
+    }
     if (fFunc)
-        *fLog << " Weight function:  " << fFunc->GetTitle()   << endl;
+        *fLog << " Weight func:  " << fFunc->GetTitle()   << endl;
 }
 
@@ -436 +526 @@
 }
 
+/*
+ * This could be used to improve the Zd-weighting within a bin.
+ * Another option is to use more bins, or both.
+ * Note that it seems unnecessary, because the shape within the
+ * theta-bins should be similar in data and Monte Carlo... hopefully.
+ *
+void MMcSpectrumWeight::InitZdWeights()
+{
+    TH2D w(*fWeightsZd);
+
+    for (int i=1; i<=w.GetNbinsX(); i++)
+    {
+         const Double_t tmin = w.GetBinLowEdge(i)  *TMath::DegToRad();
+         const Double_t tmax = w.GetBinLowEdge(i+1)*TMath::DegToRad();
+
+         const Double_t wdth  = tmax-tmin;
+         const Double_t integ = cos(tmin)-cos(tmax);
+
+         w.SetBinContent(i, w.GetBinContent(i)*wdth/integ);
+    }
+
+    //  const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
+    //  const Double_t theta = fPointing->GetZd()*TMath::DegToRad();
+    //  w = sin(theta)*w.GetBinContent(i);
+}
+*/
+
 // ----------------------------------------------------------------------------
 //
@@ -447 +564 @@
     if (fWeightsZd)
     {
-        /*
-         const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
-
-         Double_t tmin = fWeightsZd->GetXaxis()->GetBinLowEdge(i)  *TMath::DegToRad();
-         Double_t tmax = fWeightsZd->GetXaxis()->GetBinLowEdge(i+1)*TMath::DegToRad();
-         Double_t wdth = fWeightsZd->GetXaxis()->GetBinWidth(i)    *TMath::DegToRad();
-
-         Double_t cont = fWeightsZd->GetBinContent(i);
-
-         Double_t integ = cos(tmin)-cos(tmax);
-
-         w = sin(tmax)*cont/integ*wdth;
-         */
-
         const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
         w = fWeightsZd->GetBinContent(i);
@@ -472 +575 @@
 
     const Double_t e = fMcEvt->GetEnergy();
+
     fWeight->SetVal(fFunc->Eval(e)*w);
 

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h

@@ -77 +77 @@
 
     // Getter
-    TString GetFormulaSpecOld() const;
-    TString GetFormulaSpecNew() const;
-    TString GetFormulaWeights() const;
+    TString GetFormulaSpecOld(const char *name) const;
+    TString GetFormulaSpecNew(const char *name) const;
+    TString GetFormulaWeights(const char *name) const;
+    TString GetFormulaSpecOld() const { return GetFormulaSpecOld(fNameMcEvt); }
+    TString GetFormulaSpecNew() const { return GetFormulaSpecNew(fNameMcEvt); }
+    TString GetFormulaWeights() const { return GetFormulaWeights(fNameMcEvt); }
 
     TString GetFormulaSpecOldX() const { return ReplaceX(GetFormulaSpecOld()); }
@@ -85 +88 @@
     TString GetFormulaWeightsX() const { return ReplaceX(GetFormulaWeights()); }
 
-    Double_t GetSpecNewIntegral() const;
-    Double_t GetSpecOldIntegral() const;
+    Double_t GetSpecNewIntegral(Double_t emin, Double_t emax) const;
+    Double_t GetSpecOldIntegral(Double_t emin, Double_t emax) const;
+
+    Double_t GetSpecNewIntegral() const { return GetSpecNewIntegral(fEnergyMin, fEnergyMax); }
+    Double_t GetSpecOldIntegral() const { return GetSpecOldIntegral(fEnergyMin, fEnergyMax); }
 
     Double_t CalcSpecNew(Double_t e) const;
@@ -93 +99 @@
     Double_t GetEnergyMin() const { return fEnergyMin; }
     Double_t GetEnergyMax() const { return fEnergyMax; }
+
+    // Functions
+    void CompleteEnergySpectrum(TH1 &h, Double_t emin) const;
 
     // TObject