Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestGround - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TestGround.cc,v 1.35 2005/10/06 08:09:55 moreggia Exp $
// S. Moreggia created 27 October 2003

#include "euso.hh"
#include "TestGround.hh"
#include "utils.hh"
#include <stdlib.h>
#include "EConst.hh"
#include "EsafRandom.hh"
#include "SinglePhoton.hh"
#include "ListPhotonsInAtmosphere.hh"
#include <TF1.h>

ClassImp(TestGround)

using EConst::EarthRadius;
using namespace TMath;
using namespace sou;


//_____________________________________________________________________________________________________
 TestGround::TestGround() {
    //
    // ctor
    //
    Configure();
    Msg(EsafMsg::Info) << "*** TestGround built ***" << MsgDispatch;
}

//_____________________________________________________________________________________________________
 TestGround::~TestGround(){
    //
    // dtor
    //
}

//_____________________________________________________________________________________________________
 void TestGround::Configure() {
    //
    // Configure the object
    //
    fType = Conf()->GetStr("TestGround.fType");
    fAlbedo  = Conf()->GetNum("TestGround.fAlbedo");
    fAltitude = Conf()->GetNum("TestGround.fAltitude")*km;
}

//_____________________________________________________________________________________________________
 void TestGround::RandomDir(const EarthVector& pos, EarthVector& incom_dir) const {
    //
    // set incom_dir to the outgoing direction of scattering, sampling Outgoing phase function
    //
    
    // get random LOCAL theta (not in the MES)
    TRandom* rndm = EsafRandom::Get();
    Double_t theta = acos(sqrt(rndm->Rndm())); // 2pi * phase_func(theta) * sin(theta)   --- cos(theta) = x ----->  dProba/dx = 2x
    
    // then get a random direction in the MES
    EarthVector local_Uz(pos);
    local_Uz(2) += EarthRadius();           // now local_Uz is the local vertical direction
    local_Uz.Unit();
    EarthVector v = local_Uz.Orthogonal();
    v.Rotate(rndm->Rndm()*2*Pi(),local_Uz);
    EarthVector axis = v.Cross(local_Uz);
    local_Uz.Rotate(theta,axis);
    incom_dir = local_Uz; // outgoing direction
    
#ifdef DEBUG
    if(IsUnderGround(pos + (100*m)*incom_dir)) Msg(EsafMsg::Warning) << "<RandomDir> Outgoing direction is not UPWARD, it should"<<MsgDispatch;
#endif
}

//_____________________________________________________________________________________________________
 Double_t TestGround::Outgoing_phase_function(const EarthVector& pos1, const EarthVector& pos2) const {
    //
    // Phase function of the outgoing light, emitted from pos1 toward pos2
    // earth sphericity taken into account
    // Phase function is normalized when integrated over an hemisphere, i.e. theta within [0,Pi()/2]
    //
    Double_t rtn(0);
    Double_t angle;
    EarthVector temp(pos1);
    temp(2) += EarthRadius();            // pos1 expressed in earth-centered frame -> gives local vertical in MES frame
    angle = (pos2 - pos1).Angle(temp);   // angle between outgoing direction and local vertical

    if(fType == "isotropic")
        rtn = 1./TwoPi();
    else if(fType == "lambertian")
        rtn = cos(angle)/Pi();
    else invalid_argument("unvalid ground type for brdf");

    return rtn;
}

//_____________________________________________________________________________________________________
 EarthVector TestGround::GetImpact(const EarthVector& pos, const EarthVector& dir) const {
    //
    // Impact on ground form a given (position,direction) pair in atmosphere
    //

    // if pos is underground
    EarthVector rtn;
    if((pos.Zv() + 1*mm) < fAltitude) {
        Msg(EsafMsg::Debug) << "Starting position is underground, so no ground impact calculated" << MsgDispatch;
        rtn.SetXYZ(0,0,-HUGE);
        return rtn;
    }

#ifdef DEBUG
   // Msg(EsafMsg::Debug) <<"Impact on ground = " << MsgDispatch;
#endif

    // to know if dir is locally going upward
    Double_t angle;
    EarthVector temp(pos);
    temp(2) += EarthRadius();  // pos expressed in earth-centered frame -> gives local vertical direction expressed in MES frame
    angle = dir.Angle(temp);   // angle between dir and local vertical
    if(angle < (PiOver2() - kTolerance)) {
      //Msg(EsafMsg::Debug) << MsgDispatch;
        //Msg(EsafMsg::Debug) <<"NO IMPACT ON GROUND, track is locally going upward" << MsgDispatch;
        rtn.SetXYZ(0,0,HUGE);
        return rtn;
    }

    // now impact calculation
    Double_t mag(0);
    EarthVector direc = dir.Unit();

    // spherical earth
    Double_t b = pos*direc + direc.Z()*EarthRadius();
    Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - fAltitude) - fAltitude*fAltitude;
    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
    if(p.first == 0) {
#ifdef DEBUG
        //Msg(EsafMsg::Debug) <<"NO IMPACT ON GROUND" << MsgDispatch;
#endif
        rtn.SetXYZ(0,0,HUGE);
        return rtn;
    }
    if(p.first == 1) mag = p.second[0];
    else if(p.first ==2) mag = min(p.second[0],p.second[1]);
    rtn = pos + mag*direc;

#ifdef DEBUG
    //Msg(EsafMsg::Debug) << rtn << MsgDispatch;
#endif

    return rtn;
}


About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005