USCMuonProfile.h

Go to the documentation of this file.

#ifndef h_USCMuonProfile
#define h_USCMuonProfile

#include <utl/AugerUnits.h>
#include <tls/VMuonProfile.h>
#include <tls/MuonProfileUtilities.h>

#include <cmath>
#include <vector>

class TFile;

class TH2D;

namespace utl { class Branch; }

namespace USCMuonProfileNS {

 class USCMuonProfile : public tls::VMuonProfile {
  public:
    USCMuonProfile(utl::Branch branch);
    ~USCMuonProfile() {}

    inline
    void SetProfile(double theta, double phi, double energy)
    {
      if (theta < fThetaRange[0]) theta = GetThetaMin();
      else if (theta > fThetaRange[1]) theta = GetThetaMax();

      if (theta == fTheta && phi == fPhi && energy == fEnergy) return;

      fTheta  = theta;
      fPhi    = phi;
      fEnergy = energy;
      fBPsi   = tls::GetBpsi(theta, phi);

      fNMuonScaleEnergy = fNMuonScale*std::pow(energy/(1e19*utl::eV),fExponent);

      LoadHistograms(fTheta, fPhi);
    }

    double NMuon(double xpos, double ypos, double theta, double phi, double energy);
    double ThetaMuon(double xpos, double ypos, double theta, double phi, double energy);
    double EnergyMuon(double xpos, double ypos, double theta, double phi, double energy);

    double GetThetaMax();
    double GetThetaMin();

  private:
    TFile* fDensity;

    std::vector<double> fThetaRange;

    double fTheta, fPhi, fEnergy;
    double fBPsi;

    double fExponent;

    double fNMuonScale;
    double fNMuonScaleEnergy;

    unsigned int fIndexThetaLow, fIndexThetaUp;
    unsigned int fIndexPhiLow, fIndexPhiUp;

    TH2D* fHistDen1;
    TH2D* fHistDen2;
    TH2D* fHistDen3;
    TH2D* fHistDen4;

    void LoadHistograms(double theta, double phi);

    unsigned int Theta2IndexLower(double theta);
    unsigned int Theta2IndexUpper(double theta);
    unsigned int Phi2IndexLower(double theta);
    unsigned int Phi2IndexUpper(double theta);

    double Index2Theta(unsigned int itheta);
    double Index2Phi(unsigned int iphi);

    double PhiMap(double phi);

    // rounding function to handle the problems with floating point precission and the indexing system.
    inline
    double approx(double num, double prec = 1e-3)
    {
      return floor(num/prec+0.5)*prec;
    }
  };
}

#endif