// ESAF : Euso Simulation and Analysis Framework
// $Id: TestGround.cc,v 1.29 2005/04/14 15:57:47 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"
ClassImp(TestGround)
using EConst::EarthRadius;
using EConst::C;
//_____________________________________________________________________________________________________
TestGround::TestGround() {
//
// ctor
//
Msg(EsafMsg::Info) << "*** TestGround built ***" << MsgDispatch;
Configure();
}
//_____________________________________________________________________________________________________
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;
}
//_____________________________________________________________________________________________________
Double_t TestGround::Brdf(const EarthVector& pos1, const EarthVector& pos2) const {
//
// Bidirectionnal reflectance distribution function
// Returns the value for light going toward EUSO detector from the position given in parameter
//
Double_t rtn(0);
if(fType == "isotropic")
rtn = 1./twopi;
else if(fType == "lambertian")
rtn = cos((pos2 - pos1).Theta())/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() < 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 < halfpi) {
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;
}