TimeModel.h

Go to the documentation of this file.

#ifndef h_TimeModel
#define h_TimeModel

using namespace std;

#include "TMath.h"
#include "TF1.h"
#include "TObject.h"

#include <evt/ShowerSimData.h>

class TProfile;

namespace utl {
    class TabulatedFunction;
};


class TimeModel {

 public:

    TimeModel (double angle_in, // angle in deg
               int primary,
               int Da_flag_in=1);

    TimeModel (double angle_in, // angle in deg
               const utl::TabulatedFunction& MuonProfile,
               double Xobs, 
               evt::ShowerSimData::ProfileType profType,
               int Da_flag_in=1);

    ~TimeModel();



    void Info(void);


    
    TProfile *Get_dNdz (float zMin, float zMax, int nBin);
    TProfile *Get_dNdX (float zMin, float zMax, int nBin,
                        double cosTheta, double HeightObs);
    TProfile *Get_dNdz_FromProfile();
    TProfile *Get_dNdX_FromProfile (double cosTheta, double HeightObs);
    
private:

    bool fProductionHeightFromProfile;
    utl::TabulatedFunction *fMuonProductionHeightDistribution;
    double fZmin;
    double fZmax;

    double m;
    double m2;
    double pk;
    double kappa;
    double lambda;
    double ptc;
    double gam;
    double c;
    double logt_e_up;

    double logMean;
    double logSigma;
    double logLambda;
    double angle;
    double zMean;
  
    double logt_g_up;
    double logt_g_low;

    double Delta;
        
    double MomentumNumber;
    double int_e;
    double int_g;

    int Da_flag;


    void MakeProductionHeightParameters (double thetadeg, 
                                         int primary /*,double logE*/);

    
    void MakeProductionHeightDistribution (double thetadeg, 
                                 const utl::TabulatedFunction& MuonProfile,
                                           double XobsVertical);
   
    void ConvertProductionHeightDistribution (double thetadeg, 
                     const utl::TabulatedFunction& MuonProdProfile,
                                              double Xobs); // vert. depth

    double L_t(double z,double r);
    double cosaDa(double z,double r,double Delta);
    double z_t(double t,double r);
    double dzdt(double t,double r);
    double dNdlogz(double logz);
    double dNdz(double z);
    double g_t(double t,double r,double Delta,double n);
    double fg_t(double *t,double *r);
   
    
    double E(double t,double x);
    double dEdt(double t,double x);
    double dNdE(double E,double x,double r);
    double dNdlogE(double logE,double x,double r);
    double fdNdlogE(double *logE,double *p);
    double e_t(double t,double r,double z,double n);
    double fe_t(double *t,double *r);
     
//    double TotaldNdt(double t,double r,double z,double Delta, double n);

 public:
    TF1 *Fe_logt;
    TF1 *FtMe_logt;
    TF1 *Fg_logt;
    TF1 *FtMg_logt;
    TF1 *FTotaldNdlogt;
   

 double fe_logt(double *logt,double *r);
 double ftMe_logt(double *logt,double *r);
 double fg_logt(double *logt,double *r);
 double ftMg_logt(double *logt,double *r);
 double fTotaldNdlogt(double *logt,double *r);


 
    double SetCoordinates(double r,double psi);//psi rad
    void SetMomentumNumber(double MomentumNumber_in);
    double GetMomentumNumber();
    
    double GetFirstTime(int N=1);
    double GetTimes(int N,double *at);
    double GetMeanTime(int N=1);
    void GetFirstAndMeanTime(double& t_first,double& t_mean,int N=1);
    double GetLastTime(int N=1);
    void GetMeanAndRMSOfFirstTime(double& mean_t1,double& RMS_t1,int N=1,int stats=1000);
    double GetDeltaTime();
    double Eval_e_t(double t);

    double Get_tM_g_t(double MomentumNumber_in);
    double Get_tM_e_t(double MomentumNumber_in);
    double Get_tM_e_t(double MomentumNumber_in,double z);

    double TotaldNdt(double t,double r,double Delta, double n);
 

    double function(double logt,double level);
    double get_function_zero(double level);
    double GetRiseTime(double down=0.1,double up=0.5);
    
};
    
#endif