Ignore:
Timestamp:
03/18/08 17:47:37 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhflux
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

    r8746 r8883  
    428428// and fEnergyMax correctly describe the range.
    429429//
     430// If scale is given the histogram statistics is further extended by the
     431// new spectrum according to the scale factor (eg. 1.2: by 20%)
     432//
    430433// In the 1D case it is assumed that the x-axis is a zenith angle binning.
    431434// In the 2D case the x-axis is assumed to be zenith angle, the y-axis
    432435// to be energy.
    433436//
    434 void MMcSpectrumWeight::CompleteEnergySpectrum(TH1 &h, Double_t emin) const
     437void MMcSpectrumWeight::CompleteEnergySpectrum(TH1 &h, Double_t emin, Double_t scale) const
    435438{
    436439    if (h.InheritsFrom(TH3::Class()))
    437440    {
    438         *fLog << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum doesn't support TH3." << endl;
    439441        return;
    440442    }
    441443
    442     if (fEnergyMin <= emin)
     444    if (fEnergyMin < emin)
     445    {
     446        *fLog << err << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum: fEnergyMin (";
     447        *fLog << fEnergyMin << ") smaller than emin (" << emin << ")." << endl;
    443448        return;
    444 
     449    }
     450
     451    // Total number of events for the new spectrum in the same
     452    // energy range as the current histogram is filled
    445453    const Double_t norm = GetSpecNewIntegral();
    446454
     455    // Check if it is only a histogram in ZA
    447456    if (!h.InheritsFrom(TH2::Class()))
    448457    {
    449         h.Scale(GetSpecNewIntegral(emin, fEnergyMax)/norm);
     458        // Warning: Simply scaling the zenith angle distribution might
     459        // increase fluctuations for low statistics.
     460        const Double_t f = GetSpecNewIntegral(emin, fEnergyMax)/norm;
     461        h.Scale(f*scale);
    450462        return;
    451463    }
     
    453465    const TAxis &axey = *h.GetYaxis();
    454466
     467    // Find energy range between the minimum energy to be filles (emin)
     468    // and the minimum energy corresponding to the data filled into
     469    // this histogram (fEnergyMin)
    455470    const Int_t first = axey.FindFixBin(emin);
    456471    const Int_t last  = axey.FindFixBin(fEnergyMin);
     472    const Int_t max   = axey.FindFixBin(fEnergyMax);
    457473
    458474    for (int x=1; x<=h.GetNbinsX(); x++)
    459475    {
    460          const Double_t f = h.Integral(x, x, -1, 9999)/norm;
    461 
    462         for (int y=first; y<=last; y++)
    463         {
    464             const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
    465             const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMin ? fEnergyMin : axey.GetBinLowEdge(y+1);
    466 
    467             h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi));
    468         }
     476        // Ratio between the number of events in the zenith angle
     477        // bin corresponding to x and the new spectrum.
     478        const Double_t f = h.Integral(x, x, -1, 9999)/norm;
     479
     480        // Fill histogram with the "new spectrum" between
     481        // emin and fEnergyMin.
     482        if (emin<fEnergyMin)
     483            for (int y=first; y<=last; y++)
     484            {
     485                // Check if the bin is only partly filled by the energy range
     486                const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
     487                const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMin ? fEnergyMin : axey.GetBinLowEdge(y+1);
     488
     489                // Add the new spectrum extending the existing spectrum
     490                h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi));
     491            }
     492
     493        // If scale is >1 we also have to increse the statistics f the
     494        // histogram according to scale.
     495        if (scale>1)
     496            for (int y=first; y<=max; y++)
     497            {
     498                // Check if the bin is only partly filled by the energy range
     499                const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
     500                const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMax ? fEnergyMax : axey.GetBinLowEdge(y+1);
     501
     502                // Use the analytical solution to scale the histogram
     503                h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi)*(scale-1));
     504            }
    469505    }
    470506}
  • trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h

    r8680 r8883  
    101101
    102102    // Functions
    103     void CompleteEnergySpectrum(TH1 &h, Double_t emin) const;
     103    void CompleteEnergySpectrum(TH1 &h, Double_t emin, Double_t scale=0) const;
    104104
    105105    // TObject
Note: See TracChangeset for help on using the changeset viewer.