Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestCloudsPropagator - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TestCloudsPropagator.cc,v 1.13 2005/10/21 13:27:38 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 reached inside clouds taken into account
//
// Approx : - BunchOfPhotons are attenuated after having gone through clouds (not during they're inside)
//          - 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 namespace TMath;
using namespace EConst;

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())/Clight();
			
	// 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);
    Double_t angle;
    EarthVector temp(pos);
    temp(2) += EarthRadius();            // pos expressed in earth-centered frame -> gives local vertical in MES frame
    angle = (EUSO() - pos).Angle(temp);  // angle between outgoing direction and local vertical

    if(clouds->GetType() == "isotropic")
        rtn = 1./TwoPi();
    else if(clouds->GetType() == "lambertian")
        rtn = cos(angle)/Pi();
    else Msg(EsafMsg::Panic) << "unvalid clouds type for brdf" <<MsgDispatch;

    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