TimeModel.cc

Go to the documentation of this file.

#include "TimeModel.h"
#include "TMath.h"
#include <TProfile.h>

#include <utl/TabulatedFunction.h>
#include <utl/UTMPoint.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
using namespace utl;

#include <det/Detector.h>

#include <evt/ShowerSimData.h>
using namespace evt;

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
using namespace atm;

#include <vector>
using namespace std;



TimeModel *gThisTimeModel;


double LinkTofe_logt(double *x,double *p) {
    return gThisTimeModel->fe_logt(x,p);
}

double LinkToftMe_logt(double *x,double *p) {
    return gThisTimeModel->ftMe_logt(x,p);
}

double LinkTofg_logt(double *x,double *p) {
    return gThisTimeModel->fg_logt(x,p);
}

double LinkToftMg_logt(double *x,double *p) {
    return gThisTimeModel->ftMg_logt(x,p);
}

double LinkTofTotaldNdlogt(double *x,double *p) {
    return gThisTimeModel->fTotaldNdlogt(x,p);
}



TimeModel::TimeModel (double angle_in,        // in deg
                      int primary,            // 0:photon, 1:proton, 2: iron
                      int Da_flag_in) :       // volume detector: 1 flat: 0
    fMuonProductionHeightDistribution (0) {

            gThisTimeModel=this;

            Da_flag=Da_flag_in;

            /* Model Constants */
             m=0.105;
             m2=0.011;
             pk=0.0002;
             kappa=0.8;
             lambda=1;
             ptc=0.17;
             gam=2.6;
             c=0.2997924580;

             logt_e_up=3.23;
             /*********************/
             MomentumNumber=1;
            
            
            angle = angle_in;   // in degrees
            fProductionHeightFromProfile = false; // use parameters
            MakeProductionHeightParameters (angle, primary /*,logE*/);

            zMean = pow (10,logMean);


            char nome[100];
            sprintf(nome,"e_logt_%.1f",angle);
            Fe_logt=new TF1(nome,LinkTofe_logt,-6.,5.,2);
         
            sprintf(nome,"tMe_logt_%.1f",angle);
            FtMe_logt=new TF1(nome,LinkToftMe_logt,-6.,5.,2);
          
            sprintf(nome,"g_logt_%.1f",angle);
            Fg_logt=new TF1(nome,LinkTofg_logt,-4.,6.,2);
         
            sprintf(nome,"tMg_logt_%.1f",angle);
            FtMg_logt=new TF1(nome,LinkToftMg_logt,-4.,6.,2);

            sprintf(nome,"TotaldNdlogt_%.1f",angle);
            FTotaldNdlogt=new TF1(nome,LinkTofTotaldNdlogt,-4.,6.,3);

}


TimeModel::TimeModel (double angle_in,        // in deg
                      const utl::TabulatedFunction& MuonProfile,  
                      double XobsVertical,
                      evt::ShowerSimData::ProfileType profType,
                      int Da_flag_in) : 
    fMuonProductionHeightDistribution (0) {

    
            gThisTimeModel=this;
            
            fProductionHeightFromProfile = true;
            if (profType == evt::ShowerSimData::eMuon) { 

                MakeProductionHeightDistribution (angle_in,
                                                  MuonProfile,
                                                  XobsVertical);

            } else {

                ConvertProductionHeightDistribution (angle_in,
                                                     MuonProfile,
                                                     XobsVertical);

            }

            Da_flag=Da_flag_in;

            /* Model Constants */
             m=0.105;
             m2=0.011;
             pk=0.0002;
             kappa=0.8;
             lambda=1;
             ptc=0.17;
             gam=2.6;
             c=0.2997924580;

             logt_e_up=3.23;
             /*********************/
             MomentumNumber=1;
            
            
            angle = angle_in;   // in degrees

            zMean = pow(10,logMean);



            char nome[100];
            sprintf(nome,"e_logt_%.1f",angle);
            Fe_logt=new TF1(nome,LinkTofe_logt,-6.,5.,2);
         
            sprintf(nome,"tMe_logt_%.1f",angle);
            FtMe_logt=new TF1(nome,LinkToftMe_logt,-6.,5.,2);
          
            sprintf(nome,"g_logt_%.1f",angle);
            Fg_logt=new TF1(nome,LinkTofg_logt,-4.,6.,2);
         
            sprintf(nome,"tMg_logt_%.1f",angle);
            FtMg_logt=new TF1(nome,LinkToftMg_logt,-4.,6.,2);

            sprintf(nome,"TotaldNdlogt_%.1f",angle);
            FTotaldNdlogt=new TF1(nome,LinkTofTotaldNdlogt,-4.,6.,3);

}



TimeModel::~TimeModel() {

            gThisTimeModel=this;

            delete Fe_logt;
            delete FtMe_logt;
            delete Fg_logt;
            delete FtMg_logt;
            
            delete FTotaldNdlogt;
            
            delete fMuonProductionHeightDistribution;
}




void TimeModel::Info(void) {
    
    gThisTimeModel=this;

            printf("\n");
            printf("*******************************\n");
            printf("* MODEL CONSTANTS             *\n");
            printf("*******************************\n");
            printf("* pk             %5.2f GeV/km *\n",pk*1000);
            printf("* kappa         %5.1f         *\n",kappa);
            printf("* lambda        %5.1f         *\n",lambda);
            printf("* gamma         %5.1f         *\n",gam);
            printf("* ptc or Q       %5.2f GeV    *\n",ptc);
            printf("*******************************\n");
            printf("*******************************\n");
            printf("* SHOWER CONSTANTS            *\n");
            printf("*******************************\n");
            printf("* angle         %5.1f  deg    *\n",angle);
            printf("* logMean /m      %5.3f       *\n",logMean);
            printf("* logSigma /m     %5.3f       *\n",logSigma);
            printf("* logLambda /m    %5.3f       *\n",logLambda);
            printf("*******************************\n");
            if(Da_flag==1)
            printf("* D(a)=1                      *\n"); else
            printf("* D(a)=cos(a)+Delta/L         *\n");
            printf("*******************************\n");        
            /*
            printf("*******************************\n");
            printf("* COORDINATES                 *\n");
            printf("*******************************\n");
            printf("* r        %10.1f  m    *\n",r);
            printf("* psi           %5.1f  deg    *\n",psi);
            printf("*******************************\n");
            */

}







void TimeModel::MakeProductionHeightParameters (double thetadeg, 
                                                int primary /*, 
                                                double logE*/) {

    double c = cos (thetadeg*TMath::DegToRad());

    switch (primary) {
        
        case 0:
        
            logMean=
                +5.712
                -6.815 *c
                +13.82 *c*c
                -15.28 *c*c*c
                +6.124 *c*c*c*c;
            
            logSigma=
                +0.01546 
                +0.594   *c
                -1.499   *c*c
                +1.000   *c*c*c
                -1.571   *c*c*c*c
                -1.496   *c*c*c*c*c;
            /*
            Norm=
                -6.79e5
                +1.81e7 *c
                -9.30e7  *c*c
                +2.18e8 *c*c*c
                -1.33e8 *c*c*c*c;
            */
            logLambda=
                -0.198
                +2.82  *c
                -19.25 *c*c
                +54.26 *c*c*c
                -57.43 *c*c*c*c
                -20.61 *c*c*c*c*c;
            break;



        case 1:
        
            logMean=
                +5.682 
                -5.524 *c
                +9.180 *c*c
                -8.277 *c*c*c
                +2.797 *c*c*c*c;
            
            logSigma=
                +0.03266 
                +0.456   *c
                -1.084   *c*c
                +1.199   *c*c*c
                -0.4366  *c*c*c*c;
            /*
            Norm=
                -3.673E5 
                +4.467E7 *c
                -9.47E7  *c*c
                +2.917E8 *c*c*c
                -1.781E8 *c*c*c*c;
            */
            logLambda=
                -0.01111 
                +0.2581 *c
                -0.7385 *c*c
                +3.362  *c*c*c
                -2.261  *c*c*c*c;
            break;


        case 2:
        
            logMean=
                +5.692 
                -5.490 *c
                +9.147 *c*c
                -8.228 *c*c*c
                +2.782 *c*c*c*c;
            
            logSigma=
                +0.03432 
                +0.457   *c
                -1.028   *c*c
                +1.047   *c*c*c
                -0.3476  *c*c*c*c;
            /*
            Norm=
                -6.72e4 
                +5.23e7  *c
                -8.03e7  *c*c
                +2.82e8  *c*c*c
                -1.71e8  *c*c*c*c;
            */
            logLambda=
                -0.00579 
                +0.154 *c
                -0.239 *c*c
                +2.292  *c*c*c
                -1.626  *c*c*c*c;
            break;
    }
    
    
    if (logLambda<0.03) 
        logLambda = 0.03;

}





void TimeModel::MakeProductionHeightDistribution (double thetadeg, 
                          const utl::TabulatedFunction& MuonProfile,
                                                  double Xobs) { // vert. depth


    double cosTheta = cos (thetadeg*deg);

    // convert X (slant depth) to z / m
    // taking care of earth curvature
    det::Detector& Det = det::Detector::GetInstance();
    
    const ProfileResult heightProfile 
        = Det.GetAtmosphere().EvaluateHeightVsDepth();
    const ProfileResult depthProfile 
        = Det.GetAtmosphere().EvaluateDepthVsHeight();

    double HeightXobs = heightProfile.Y (Xobs);

    vector <double> ValueZ;
    vector <double> ValuedNdZ;

    double AtmDepthMin = depthProfile.Y (depthProfile.MaxX());
    double AtmDepthMax = depthProfile.Y (depthProfile.MinX());

    logMean = 0;
    double Norm = 0;

    bool First = true;
    double z_old = 0, HeightX_old = 0, Nmu_old = 0; 
    for (TabulatedFunction::ConstIterator iProfile = MuonProfile.Begin();
         iProfile != MuonProfile.End();
         ++iProfile) {

        double Nmu = iProfile->Y();
        double Xslant = iProfile->X();
        double Xvert = Xslant*cosTheta;
        
        if (Xvert<AtmDepthMin || Xvert>AtmDepthMax)
            continue;

        double HeightX = heightProfile.Y (Xvert);
        double z = (HeightX-HeightXobs)/cosTheta/meter;
        
        // only look at shower not before
        if (Nmu==0) {
            continue; 
        }

        /*
         // only look above the ground
         if (z<-1000) {
         continue; 
         }
        */
        
        // first bin is for init only
        if (First) {

            First = false;

            z_old = z;
            HeightX_old = HeightX;
            Nmu_old = Nmu;
            
            fZmax = 0;
            fZmin = z;
            continue;
        }
        
        double dz = z - z_old;
        double dN = Nmu_old - Nmu;
        double zMean = (z_old + z) / 2.;

        if (zMean>0) {
            logMean += log10 (zMean) * (Nmu+Nmu_old)/2.;
            Norm += (Nmu+Nmu_old)/2.;
        }

        ValueZ.insert (ValueZ.begin(), zMean);
        ValuedNdZ.insert (ValuedNdZ.begin(), dN/dz);

        if (zMean<fZmin) fZmin = zMean;
        if (zMean>fZmax) fZmax = zMean;
        
        z_old = z;
        HeightX_old = HeightX;
        Nmu_old = Nmu;
    }

    fMuonProductionHeightDistribution 
        = new TabulatedFunction (ValueZ, ValuedNdZ);
    
    
    fZmin = *(fMuonProductionHeightDistribution->XBegin());
    fZmax = *(fMuonProductionHeightDistribution->XEnd()-1);
    

    /*
    int n=ValueZ.size();    
      for (int i=0; i<n; i++) {
      cout << " B> " << ValueZ [i] << " " << ValuedNdZ[i] << endl;
      }
    */

    ValueZ.clear();
    ValuedNdZ.clear();

    logMean /= Norm;
}





void TimeModel::ConvertProductionHeightDistribution (double thetadeg, 
                          const utl::TabulatedFunction& MuonProdProfile,
                                                  double Xobs) { // vert. depth


    double cosTheta = cos (thetadeg*deg);

    // convert X (slant depth) to z / m
    // taking care of earth curvature
    det::Detector& Det = det::Detector::GetInstance();
    
    const ProfileResult heightProfile 
        = Det.GetAtmosphere().EvaluateHeightVsDepth();
    const ProfileResult depthProfile 
        = Det.GetAtmosphere().EvaluateDepthVsHeight();

    double HeightXobs = heightProfile.Y (Xobs);

    vector <double> ValueZ;
    vector <double> ValuedNdZ;

    double AtmDepthMin = depthProfile.Y (depthProfile.MaxX());
    double AtmDepthMax = depthProfile.Y (depthProfile.MinX());

    logMean = 0;
    double Norm = 0;

    double dXslant = abs (MuonProdProfile[1].X() - MuonProdProfile[0].X());
    double dXvert = dXslant * cosTheta;

    for (TabulatedFunction::ConstIterator iProfile = MuonProdProfile.Begin();
         iProfile != MuonProdProfile.End();
         ++iProfile) {

        double nDmu = iProfile->Y();
        double Xslant = iProfile->X();
        double Xvert = Xslant*cosTheta;

        
        if (Xvert-dXvert/2<AtmDepthMin || Xvert+dXvert/2>AtmDepthMax)
            continue;

        double Height = heightProfile.Y (Xvert);
        double z = (Height-HeightXobs)/cosTheta/meter;

        double Height1 = heightProfile.Y (Xvert-dXvert/2);
        double z1 = (Height1-HeightXobs)/cosTheta/meter;
        
        double Height2 = heightProfile.Y (Xvert+dXvert/2);
        double z2 = (Height2-HeightXobs)/cosTheta/meter;
        
        double dz = abs (z1-z2);

        // only look at shower not before
        if (nDmu==0) {
            continue; 
        }
        
        if (z>0) {
            logMean += log10 (z) * (nDmu);
            Norm += (nDmu);
        }

        ValueZ.insert (ValueZ.begin(), z);
        ValuedNdZ.insert (ValuedNdZ.begin(), nDmu*dXslant/dz);
        
        if (z<fZmin) fZmin = z;
        if (z>fZmax) fZmax = z;
        
    }

    fMuonProductionHeightDistribution 
        = new TabulatedFunction (ValueZ, ValuedNdZ);
    
    fZmin = *(fMuonProductionHeightDistribution->XBegin());
    fZmax = *(fMuonProductionHeightDistribution->XEnd()-1);

    /*
    int n=ValueZ.size();    
      for (int i=0; i<n; i++) {
      cout << " B> " << ValueZ [i] << " " << ValuedNdZ[i] << endl;
      }
    */

    ValueZ.clear();
    ValuedNdZ.clear();

    logMean /= Norm;
}




 double TimeModel::L_t(double z,double r)
{
    return sqrt(z*z+r*r);
}


 double TimeModel::cosaDa(double Z,double r,double Delta)
{ 
      if (Z<0) return 0; 
    // isto permite muons back!?
    double L=L_t(Z,r);
    double cosa2=1-(r*r)/(L*L);
    
     if (Da_flag==1) return sqrt(cosa2);

    double ret=cosa2+Delta/L*sqrt(cosa2);
  
    return ret;

}

double TimeModel::z_t(double t,double r)
{
    return 0.5*(r*r/(c*t)-c*t);
}

double TimeModel::dzdt(double t,double r)
{
    return 0.5*(r*r/(c*t*t)+c);
}


// RU Wed May 18 12:15:06 CEST 2005
// to get some plots
TProfile *TimeModel::Get_dNdz (float zMin, float zMax, int nBin) {
    /*
      LOG scale

    if (zMin<=0) zMin = 1e0;
    double logZmin = log10 (zMin);
    double logZmax = log10 (zMax);

    TProfile *hist = new TProfile ("lorenzo_dNdz", 
                                   "lorenzo model mu prod. height",
                                   nBin, logZmin, logZmax);

    double dz = (logZmax-logZmin)/nBin;
    for (double z = logZmin; z<logZmax; z+=dz) {
        
        hist->Fill (z+dz/2, dNdlogz (z));
    }
    */

    TProfile *hist = new TProfile ("lorenzo_dNdz", 
                                   "lorenzo model mu prod. height",
                                   nBin, zMin, zMax);

    double dz = (zMax-zMin)/nBin;
    for (double z = zMin; z<zMax; z+=dz) {
        
        hist->Fill (z+dz/2, dNdz (z));
    }
    
    return hist;
}


TProfile *TimeModel::Get_dNdX (float zMin, float zMax, int nBin,
                               double cosTheta, double HeightObs) {
    
    det::Detector& Det = det::Detector::GetInstance();

    const ProfileResult depthProfile = 
        Det.GetAtmosphere().EvaluateDepthVsHeight();

    double HeightX, Xvert;
    HeightX = zMin*meter*cosTheta + HeightObs;
    Xvert = depthProfile.Y (HeightX);
    double xMax = Xvert/cosTheta;
    HeightX = zMax*meter*cosTheta + HeightObs;
    Xvert = depthProfile.Y (HeightX);
    double xMin = Xvert/cosTheta;

    TProfile *hist = new TProfile ("lorenzo_dNdX", 
                                   "lorenzo model mu prod. height",
                                   nBin, xMin/g*cm*cm, xMax/g*cm*cm);
    double dz = (zMax-zMin)/nBin;
    for (double z = zMin; z<zMax; z+=dz) {
        
        HeightX = z*meter*cosTheta + HeightObs;
        Xvert = depthProfile.Y (HeightX);
        double X = Xvert/cosTheta;
        
        hist->Fill (X/g*cm*cm, dNdz (z));
    }
                           
    return hist;
}

TProfile *TimeModel::Get_dNdz_FromProfile() {

    int nBin = fMuonProductionHeightDistribution->GetNPoints();
    double zMin = (*fMuonProductionHeightDistribution) [0].X();
    double zMax = (*fMuonProductionHeightDistribution) [nBin-1].X();
    TProfile *hist = new TProfile ("dNdz_prof", "muon production height", 
                                   nBin, zMin, zMax);

    for (utl::TabulatedFunction::Iterator x = 
             fMuonProductionHeightDistribution->Begin();
         x!=fMuonProductionHeightDistribution->End();
         ++x) {
        
        hist->Fill (x->X()/meter, x->Y());
    }

    return hist;
}

TProfile *TimeModel::Get_dNdX_FromProfile (double cosTheta, double HeightObs) {

    det::Detector& Det = det::Detector::GetInstance();
    const ProfileResult depthProfile =
        Det.GetAtmosphere().EvaluateDepthVsHeight();

    int nBin = fMuonProductionHeightDistribution->GetNPoints();
    double zMin = (*fMuonProductionHeightDistribution) [0].X();
    double zMax = (*fMuonProductionHeightDistribution) [nBin-1].X();

    double HeightX, Xvert;
    HeightX = zMin*meter*cosTheta + HeightObs;
    Xvert = depthProfile.Y (HeightX);
    double xMax = Xvert/cosTheta;
    HeightX = zMax*meter*cosTheta + HeightObs;
    Xvert = depthProfile.Y (HeightX);
    double xMin = Xvert/cosTheta;

    TProfile *hist = new TProfile ("dNdX_prof", "charged profile", 
                                   nBin, xMin/g*cm*cm, xMax/g*cm*cm);

    for (utl::TabulatedFunction::Iterator x = 
             fMuonProductionHeightDistribution->Begin();
         x != fMuonProductionHeightDistribution->End();
         ++x) {

        HeightX = x->X()*meter*cosTheta + HeightObs;
        Xvert = depthProfile.Y (HeightX);
        double X = Xvert/cosTheta;

        hist->Fill (X/g*cm*cm, x->Y());
    }

    return hist;
}
// RU --------------------


double TimeModel::dNdlogz(double logz)
{
    if(logLambda<=0)
    {
        double e=(logz-logMean)/logSigma;
        return (exp(-0.5*e*e))/(logSigma*2.50662827463100024);
    }
    else
    { // a nova parametrizacion incluida na tese
      // estan definidos desde fora
        double norm=1;
        double s2=sqrt(2.);
        
        double f1=(logz-logMean)/logLambda;
        double f2=.5*(logSigma*logSigma)/(logLambda*logLambda);
        double f3=(logz-logMean)/(s2*logSigma);
        double f4=logSigma/(s2*logLambda);
        
        return norm*exp(f1+f2)*TMath::Erfc(f3+f4)/(2*logLambda);
        
    }
}

double TimeModel::dNdz (double z) {

    if(z<0) 
        return 0;
    
    if (fProductionHeightFromProfile) { // from dN/dLog10z profile
        
        if (z<fZmin || z>fZmax)
            return 0.;

        //cout << " z " << z << " " << fZmin << " " << fZmax << endl;

        double value = fMuonProductionHeightDistribution->Y (z);

        // don't look at decaying muon (at least if they are outnumbering 
        // the muon-production (!) and dN/dz<0
        if (value<0)
            value = 0;
        
        return value;

    } else { // from dN/dlog10z parametrisation
        
        double logz = log(z)/log(10.);
        return 1/(z*log (10.))*dNdlogz (logz);
    }
}


double TimeModel::g_t(double t,double r,double Delta,double n)
{  
    double Z=z_t(t,r);
    if(t<=0) 
        return 0;   /* || z<=0)  antiguo*/
    else 
        return n*cosaDa(Z,r,Delta)*dNdz(Z+Delta)*dzdt(t,r);
    //dzdt=dZdt
}

double TimeModel::fg_t(double *t,double *r)
{
    return g_t(t[0],r[0],r[1],1);
}

double TimeModel::fg_logt(double *logt,double *r)
{
    double t=pow(10,logt[0]);
    return t*log(10.)*g_t(t,r[0],r[1],1);
    
}
double TimeModel::ftMg_logt(double *logt,double *r)
{
    double t=pow(10,logt[0]);  
    return pow(t,MomentumNumber)*t*log(10.)*g_t(t,r[0],r[1],1);
}


double TimeModel::E(double t,double x)
{   
    return 0.5*pk*x*(1+sqrt(1+2*m2/(pk*pk*x*c*t)));
}
double TimeModel::dEdt(double t,double x)
{
    return 0.5*pk*x*1/sqrt(1+2*m2/(pk*pk*x*c*t))*m2/(pk*pk*x*c*t*t);
}

double TimeModel::dNdE(double E,double x,double r)
{
    if (E<pk*x || x<=0 || r>=x) return 0; // os dous ultimos xa non poden ser por construccion
    return (1-r*r/(x*x))*pow(r/x,lambda)*pow(E,-gam-kappa+lambda+1)*pow(E-pk*x,kappa)*exp(-E*r/(x*ptc));
}

double TimeModel::dNdlogE(double logE,double x,double r)
{
    double E=pow(10.,logE);
    return E*log(10.)*dNdE(E,x,r);
}

double TimeModel::fdNdlogE(double *logE,double *p)
{
    return dNdlogE(logE[0],p[0],p[1]);
}

double TimeModel::e_t(double t,double r,double z,double n)
{
    
    if(t<=0) return 0;
    else 
    {
        double x=sqrt(z*z+r*r);
        return n*dNdE(E(t,x),x,r)*dEdt(t,x);
    }
    
}

double TimeModel::fe_t(double *t,double *r)
{
    return e_t(t[0],r[0],r[1],r[2]);
}


 double TimeModel::fe_logt(double *logt,double *r)
{
    double t=pow(10.,logt[0]);
    return t*log(10.)*e_t(t,r[0],r[1],1);
}
  double TimeModel::ftMe_logt(double *logt,double *r)
{
    double t=pow(10,logt[0]);
    return pow(t,MomentumNumber)*t*log(10.)*e_t(t,r[0],r[1],1);
    
}
   




double TimeModel::SetCoordinates(double r,double psi)//psi rad
{ 
    gThisTimeModel=this;
            
            Delta=r*cos(psi)*tan(TMath::DegToRad()*angle);
            
            Fe_logt->SetParameter(0,r);
            Fe_logt->SetParameter(1,zMean-Delta);
            FtMe_logt->SetParameter(0,r);
            FtMe_logt->SetParameter(1,zMean-Delta);
           
            Fg_logt->SetParameter(0,r);
            Fg_logt->SetParameter(1,Delta);
            FtMg_logt->SetParameter(0,r);
            FtMg_logt->SetParameter(1,Delta);
            
            FTotaldNdlogt->SetParameter(0,r);
            FTotaldNdlogt->SetParameter(1,Delta);
            FTotaldNdlogt->SetParameter(2,1);


            double Integral;

            double LZ=(pow(10,logMean)-Delta);
          

            logt_g_up=log10(8000*(sqrt(r*r+LZ*LZ)-LZ)); /* limite superior, para estar seguro*/     
            logt_g_low=log10(0.2*(sqrt(r*r+LZ*LZ)-LZ)); /* limite superior, para estar seguro*/

            {

                Integral=Fg_logt->Integral(logt_g_low,logt_g_up);
                                
                //      if (Integral>.9) break;
                
        
            }
            
            
            Fg_logt->SetRange(logt_g_low,logt_g_up);
            FTotaldNdlogt->SetRange(logt_g_low,logt_g_up);

            //   printf("Draw()!!\n");
            //  Fg_logt->Draw();
            
              

            int_e=Fe_logt->Integral(-6,logt_e_up,(double *)NULL,0.0000001);
            int_g=Fg_logt->Integral(logt_g_low,logt_g_up,(double *)NULL,0.0000001);


            return Integral;
            
        }

void TimeModel::SetMomentumNumber(double MomentumNumber_in)
{
    MomentumNumber=MomentumNumber_in;
}

double TimeModel::GetMomentumNumber()
{
    return MomentumNumber;
}

double TimeModel::GetFirstTime(int N/*=1*/)
        {
            gThisTimeModel=this;
            double t_first=1e25;
            double t_i;
            for(int i=0;i<N;i++)
            {
                 t_i=pow(10,Fg_logt->GetRandom())+pow(10,Fe_logt->GetRandom());
                 if(t_i<t_first) t_first=t_i;
            }
            
            return t_first-Delta/c;
        }

  double TimeModel::GetTimes(int N,double *at)
        {
            gThisTimeModel=this;
            double t_first=1e25;
            double t_i;
            for(int i=0;i<N;i++)
            {
                 t_i=pow(10,Fg_logt->GetRandom())+pow(10,Fe_logt->GetRandom());
                 at[i]=t_i-Delta/c;
                 if(t_i<t_first) t_first=t_i;
            }
            
            return t_first-Delta/c;
        }


double TimeModel::GetMeanTime(int N/*=1*/)
        {
            gThisTimeModel=this;
            double t_mean=0;
            for(int i=0;i<N;i++)
            {
                 t_mean+=pow(10,Fg_logt->GetRandom())+pow(10,Fe_logt->GetRandom());
                 
            }
            return t_mean/N-Delta/c;
        }
 
void TimeModel::GetFirstAndMeanTime(double& t_first,double& t_mean,int N/*=1*/)
        {
            gThisTimeModel=this;        
            double t_i;
            t_mean=0; 
            t_first=1e25;
            for(int i=0;i<N;i++)
            {
                t_i=pow(10,Fg_logt->GetRandom())+pow(10,Fe_logt->GetRandom());
                if(t_i<t_first) t_first=t_i;
                t_mean+=t_i;
            }
            t_first-=Delta/c;
            t_mean=t_mean/N-Delta/c;
            
        }
double TimeModel::GetLastTime(int N/*=1*/)
           {
               gThisTimeModel=this;
               double t_last=0;
               double t_i;
               for(int i=0;i<N;i++)
               {
                   t_i=pow(10,Fg_logt->GetRandom())+pow(10,Fe_logt->GetRandom());
                   if(t_i>t_last) t_last=t_i;
               }
               return t_last-Delta/c;
           }

void TimeModel::GetMeanAndRMSOfFirstTime(double& mean_t1,double& RMS_t1,int N/*=1*/,int stats/*=1000*/)
        {
            gThisTimeModel=this;
            double t1sum=0,t1sum2=0,t1;
            
            for(int i=0;i<stats;i++)
            {  
                t1=GetFirstTime(N);
                t1sum+=t1;
                t1sum2+=t1*t1;
            }
            
            mean_t1= t1sum/stats;
            RMS_t1=sqrt(t1sum2/stats-mean_t1*mean_t1);
        
        }
double TimeModel::GetDeltaTime()
{
    gThisTimeModel=this;
    return Delta/c;
}
double TimeModel::Eval_e_t(double t)
{
    gThisTimeModel=this;
    return e_t(t,Fe_logt->GetParameter(0),zMean-Delta,1);
}

double TimeModel::Get_tM_g_t(double MomentumNumber_in)
{
    gThisTimeModel=this;
    MomentumNumber=MomentumNumber_in;
    return  FtMg_logt->Integral(logt_g_low,logt_g_up,(double *)NULL,0.0000001)/int_g;

}

double TimeModel::Get_tM_e_t(double MomentumNumber_in)
{
    gThisTimeModel=this;
    MomentumNumber=MomentumNumber_in;
    return  FtMe_logt->Integral(-6,logt_e_up,(double *)NULL,0.0000001)/int_e;
}


double TimeModel::Get_tM_e_t(double MomentumNumber_in,double z)
{
    gThisTimeModel=this;
    MomentumNumber=MomentumNumber_in;
    double auxz=Fe_logt->GetParameter(1);
    Fe_logt->SetParameter(1,z-Delta);
    double ret=FtMe_logt->Integral(-6,logt_e_up,(double *)NULL,0.0000001)/int_e;
    Fe_logt->SetParameter(1,auxz);
    return ret;
}

double TimeModel::TotaldNdt(double t,double r,double Delta, double n)
{
double steps=850;
double tmin=0;
 double tmax=1700;
        
        
        if (t<tmax) tmax=t;
        
        double h=(tmax-tmin)/steps;
        int N=int(steps);
        double sum=0,t_i=tmin+h/2;
        
        for(int i=0;i<N;i++)
        {
            sum+=g_t(t-t_i,r,Delta,1)*e_t(t_i,r,zMean-Delta,1);
            t_i+=h;
        }
        
    return n*sum*h;
 

}

double TimeModel::fTotaldNdlogt(double *logt,double *r)
{

     double t=pow(10.,logt[0]);
    return t*log(10.)*TotaldNdt(t,r[0],r[1],1);
}

double TimeModel::function(double logt,double level)
{
    

    
/*
    double integral=Fg_logt->Integral(logtd,logt);
    return integral/(int_g) - level;
*/
    
    double integral=FTotaldNdlogt->Integral(logt_g_low,logt);
    return integral/(int_g*int_e) - level;
  
}

double TimeModel::get_function_zero(double level)
{
    double xmin=logt_g_low;
    double xmax=logt_g_up;
    double precision=1.e-4;
    double f1=function(xmin,level);
    double f2=function(xmax,level);
    double x1=xmin;
    double x2=xmax;

    if (f1*f2<0) {
        double x0=(x1+x2)/2;
        double f0=function(x0,level);
        double epsilon= (x0-x1)/2;
        while (epsilon>=precision) {
            if (f1*f0<0) {
                x2=x0;
                f2=f0;
            }
            else
            {
                x1=x0;
                f1=f0;
            }
            epsilon= (x2-x1)/2;
            x0=(x1+x2)/2;
            f0=function(x0,level);
        }

    }
    else {
       printf("No hay zeros en el rango \n");
        return pow(10,xmin);
    }
//    printf("level %.1f time=%.1f\n",level,pow(10,x1));
    return pow(10,x1);
}


double TimeModel::GetRiseTime(double down,double up)
{    
    return (get_function_zero(up)-get_function_zero(down));
}