// ESAF : Euso Simulation and Analysis Framework
// $Id: TestCloudsPropagator.cc,v 1.8 2005/04/27 11:32:42 moreggia Exp $
// Author: Sylvain Moreggia Aug, 19 2004
//_____________________________________________________________________________
//
// TestCloudsPropagator
//
// Associated with TestClouds objects
// Generates clouds backscattered photons according to TestClouds albedo
// Lambertian angular response - depth reachec inside clouds taken into account
//
// Approx : - BunchOfPhotons are not attenuated when going through clouds (scattering strongly peaked forward)
// - SinglePhoton going toward EUSO is attenuated
// - if a bunch created inside clouds -> no clouds response to this bunch
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "TestCloudsPropagator.hh"
#include "TestClouds.hh"
#include "Atmosphere.hh"
#include "EsafRandom.hh"
#include "BunchRadiativeTransfer.hh"
#include "EConst.hh"
#include "SinglePhoton.hh"
#include "ListPhotonsInAtmosphere.hh"
#include "TestGround.hh"
using EConst::C;
ClassImp(TestCloudsPropagator)
//_____________________________________________________________________________
TestCloudsPropagator::TestCloudsPropagator() : InCloudsPropagator() {
//
// Constructor
//
}
//_____________________________________________________________________________
TestCloudsPropagator::TestCloudsPropagator(const Ground* g) : InCloudsPropagator(g) {
//
// Constructor
//
}
//_____________________________________________________________________________
TestCloudsPropagator::~TestCloudsPropagator() {
//
// Destructor
//
}
//_____________________________________________________________________________
Medium TestCloudsPropagator::Go(BunchOfPhotons& b,ListPhotonsInAtmosphere& list) const {
//
// Clouds response for a BunchOfPhotons, generating SinglePhotons (ListPhotonsInAtmosphere filled)
// TestClouds is considered as a lambertian object, with penetration in clouds taken into account
// if a bunch (or a singlephoton) is created below clouds top surface -> no clouds backscattering to EUSO
//
// bunch NOT transported here. Just clouds backscattering, and bunch weight attenuation
// transportation is done afterwards, with ClearSkyPropagator
//
// A bunch reaching this method has its position on clouds top surface
//
#ifdef DEBUG
Msg(EsafMsg::Debug) << "TestClouds propagation" << MsgDispatch;
#endif
// initializations
TRandom* rndm = EsafRandom::Get();
Medium rtn = CLEARSKY; // after clouds, bunch treatment should come back to clear sky (at least the part within clouds)
const TestClouds* clouds = (TestClouds*) Atmosphere::Get()->GetClouds();
const BunchOfPhotons& cst_b = b;
Double_t OD(0);
Int_t nb(0);
EarthVector exit;
// outgoing impact calculated. NextImpact set accordingly
exit = clouds->GetCloudImpact(cst_b.GetPos(),cst_b.GetDir(),"outgoing");
if(exit.Z() == HUGE) Msg(EsafMsg::Warning) << "outgoing impact should exist in this case" <<MsgDispatch;
b.SetNextImpact(exit);
OD = clouds->TotalOD(cst_b.GetPos(),cst_b.GetDir());
// if bunch created within clouds, no clouds backscattering, only absorption
// if created externally, here below is creation of backscattered photons
if(!clouds->IsInClouds(cst_b.GetShowerPos())) {
nb = EsafRandom::Get()->Poisson(cst_b.GetWeight() * clouds->GetAlbedo() * EusoOmega(cst_b.GetPos()) * Brdf(cst_b.GetPos()));
}
else nb = 0;
for(Int_t i=0; i<nb; i++) {
// init
Double_t wl, date, tof;
EarthVector showerpos, pos, dir, diff, step, b_impact, scat_repart;
SinglePhoton* s = 0;
UInt_t bid = cst_b.GetId();
PhotonType type = cst_b.GetType();
date = cst_b.GetDate();
// corrections for date and showerpos, from BunchOfPhotons mean values and longitudinal (not lateral) distribution at creation
showerpos = cst_b.RandomPosInShower();
// if showerpos is under clouds top surface, photons not created
if(fGround->IsUnderGround(showerpos)) continue;
if(showerpos.Zv() < clouds->GetTopAltitude()) continue;
diff = showerpos - cst_b.GetShowerPos();
date = cst_b.GetDate() + diff.Dot(cst_b.GetDir().Unit())/C();
// scattering repartition within clouds
b_impact = clouds->GetCloudImpact(cst_b.GetPos(),cst_b.GetDir(),"outgoing");
// if no impact
if(b_impact.Z() == HUGE) {
Msg(EsafMsg::Warning) << "Bunch cannot go out from clouds -> its propagation stops" << MsgDispatch;
return NONE;
}
step = b_impact - cst_b.GetPos(); // bunch travel within clouds
scat_repart = -(log(1 - rndm->Rndm()*(1 - exp(-OD))) / OD) * step; // repartition within clouds (single scattering exponential repartition)
// tof and position corrections, due to angular distributions at creation
pos = cst_b.PosRandomAngCorrec(showerpos,cst_b.GetPos()); // fake pos, used only to get direction, for new impact calculation
// impact of the photon on clouds surface - this is a trick in order to use the impact calculation method
pos = clouds->GetCloudImpact(showerpos - scat_repart,(pos - showerpos).Unit()) + scat_repart;
// if no impact
if(pos.Z() == HUGE) continue;
if(fGround->IsUnderGround(pos)) {
Msg(EsafMsg::Warning) << "Clouds Backscat photons produced under ground - Occurs if clouds bottom near ground" << MsgDispatch;
continue;
}
tof = cst_b.GetTof() * (showerpos - pos).Mag()/(cst_b.GetPos() - cst_b.GetShowerPos()).Mag();
dir = EUSO() - pos;
wl = cst_b.GetWlSpectrum().GetLambda();
s = new SinglePhoton(type,date,tof,wl,showerpos,pos,dir,CloudScat,bid);
list.Add(s);
}
b.Convolute(exp(-OD));
#ifdef DEBUG
Msg(EsafMsg::Debug)<<"nb of single produced in clouds process = "<<nb <<MsgDispatch;
#endif
return rtn;
}
//_____________________________________________________________________________________________________
Double_t TestCloudsPropagator::Brdf(const EarthVector& pos) const {
//
// Bidirectionnal reflectance distribution function
// Returns the value for light going toward EUSO detector from the position given in parameter
//
const TestClouds* clouds = (TestClouds*) Atmosphere::Get()->GetClouds();
Double_t rtn(0);
if(clouds->GetType() == "isotropic")
rtn = 1./twopi;
else if(clouds->GetType() == "lambertian")
rtn = cos((EUSO() - pos).Theta())/pi;
else Msg(EsafMsg::Panic) << "unvalid clouds type for brdf" <<MsgDispatch;
return rtn;
}