Changeset 7688 for trunk/MagicSoft/Mars/mhflux

Timestamp:

05/04/06 08:26:21 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• MMcSpectrumWeight.cc (modified) (9 diffs)
• MMcSpectrumWeight.h (modified) (4 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-2005
+!   Copyright: MAGIC Software Development, 2000-2006
 !
 !
@@ -26 +26 @@
 //////////////////////////////////////////////////////////////////////////////
 //
-//  MMcWeightEnergySlopeCalc
+//  MMcSpectrumWeight
 //               
 //  Change the spectrum of the MC showers simulated with Corsika (a power law)
@@ -104 +104 @@
     AddToBranchList("MMcEvt.fEnergy");
 
-    fNameWeight = "MWeight";
-    fNameMcEvt  = "MMcEvt";
-
-    fNewSlope  = -1;
-    fOldSlope  = -1;
-
-    fEnergyMin = -1;
-    fEnergyMax = -2;
-
-    fNorm      =  1;
+    fNameWeight  = "MWeight";
+    fNameMcEvt   = "MMcEvt";
+
+    fNewSlope    = -1;
+    fOldSlope    = -1;
+
+    fEnergyMin   = -1;
+    fEnergyMax   = -2;
+
+    fNorm        =  1;
+    fNormEnergy  =  1;
 
     fAllowChange = kFALSE;
@@ -235 +236 @@
 // Return the formula to calculate weights.
 // Is is compiled by
-//   o = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX());
-//   n = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX());
-//
-//   result: fNorm*o/n*GetFormulaNewSpec()/GetFormulaOldSpec()
+//   o1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX());
+//   n1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX());
+//   o2 = CalcSpecOld(fNormEnergy);
+//   n2 = CalcSpecNew(fNormEnergy);
+//
+//   result (fNormEnergy<0):
+//      fNorm*o1/n1*GetFormulaNewSpec()/GetFormulaOldSpec()
+//
+//   result (fNormEnergy>0):
+//      fNorm*o2/n2*GetFormulaNewSpec()/GetFormulaOldSpec()
 //
 // fNorm is 1 by default but can be overwritten using SetNorm()
@@ -254 +261 @@
         return Form("%.16f", fNorm);
 
-    const Double_t iold = GetSpecOldIntegral();
-    const Double_t inew = GetSpecNewIntegral();
+    const Double_t iold = fNormEnergy<0 ? GetSpecOldIntegral() : CalcSpecOld(fNormEnergy);
+    const Double_t inew = fNormEnergy<0 ? GetSpecNewIntegral() : CalcSpecNew(fNormEnergy);
 
     const Double_t norm = fNorm*iold/inew;
@@ -282 +289 @@
     TF1 funcold("Dummy", GetFormulaSpecOldX());
     return funcold.Integral(fEnergyMin, fEnergyMax);
+}
+
+// ---------------------------------------------------------------------------
+//
+// Returns the value of GetFormulaSpecNewX() at the energy e describing
+// the destination spectrum
+//
+Double_t MMcSpectrumWeight::CalcSpecNew(Double_t e) const
+{
+    TF1 funcnew("Dummy", GetFormulaSpecNewX());
+    return funcnew.Eval(e);
+}
+
+// ---------------------------------------------------------------------------
+//
+// Returns the value of GetFormulaSpecOldX() at the energy e describing
+// the simulated spectrum
+//
+Double_t MMcSpectrumWeight::CalcSpecOld(Double_t e) const
+{
+    TF1 funcnew("Dummy", GetFormulaSpecOldX());
+    return funcnew.Eval(e);
 }
 
@@ -392 +421 @@
     *fLog << " New spectral slope:       " << fNewSlope << endl;
     *fLog << " User normalization:       " << 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;
@@ -446 +476 @@
 //
 // Read the setup from a TEnv, eg:
-//   MMcSpectrumWeight.NewSlope: -2.6
-//   MMcSpectrumWeight.Norm:      1.0
-//   MMcSpectrumWeight.Formula:  pow(X, -2.6)
+//   MMcSpectrumWeight.NewSlope:   -2.6
+//   MMcSpectrumWeight.Norm:        1.0
+//   MMcSpectrumWeight.NormEnergy:  200
+//   MMcSpectrumWeight.Formula:     pow(X, -2.6)
 //
 Int_t MMcSpectrumWeight::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
@@ -463 +494 @@
         SetNorm(GetEnvValue(env, prefix, "Norm", fNorm));
     }
+    if (IsEnvDefined(env, prefix, "NormEnergy", print))
+    {
+        rc = kTRUE;
+        SetNormEnergy(GetEnvValue(env, prefix, "NormEnergy", fNormEnergy));
+    }
     if (IsEnvDefined(env, prefix, "Formula", print))
     {

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h

@@ -40 +40 @@
 
     Double_t fNorm;         // Normalization constant (additional normalization constant)
+    Double_t fNormEnergy;   // Energy at which the spectra are normalized (default -1 means the integral is used)
 
     TString fFormula;       // Text Formula for new spectrum: eg. "pow(MMcEvt.fEnergy, -2.0)"
@@ -66 +67 @@
     void SetNewSlope(Double_t s=-1) { fNewSlope = s; }
     void SetNorm(Double_t s=1) { fNorm = s; }
+    void SetNormEnergy(Double_t s=1) { fNormEnergy = s; }
     void SetFormula(const char *f="") { fFormula = f; }
     void SetEnergyRange(Double_t min=-2, Double_t max=-1) { fEnergyMin=min; fEnergyMax=max; }
@@ -71 +73 @@
     void SetWeightsZd(TH1 *h=0) { fWeightsZd = h; }
     void SetWeightsSize(TH1D *h=0);
+
     Bool_t Set(const MMcCorsikaRunHeader &h);
 
@@ -85 +88 @@
     Double_t GetSpecOldIntegral() const;
 
+    Double_t CalcSpecNew(Double_t e) const;
+    Double_t CalcSpecOld(Double_t e) const;
+
     Double_t GetEnergyMin() const { return fEnergyMin; }
     Double_t GetEnergyMax() const { return fEnergyMax; }