Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestClouds - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TestClouds.cc,v 1.16 2005/10/19 12:10:26 moreggia Exp $
// Sylvain Moreggia created Apr,  8 2004

#include "TestClouds.hh"
#include "EarthVector.hh"
#include "Config.hh"
#include "EConst.hh"
#include "EsafRandom.hh"

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

// ctor
//______________________________________________________________________________
 TestClouds::TestClouds() : Clouds() {
    //
    // ctor
    //
    Configure();
    Build();
}

//______________________________________________________________________________
 TestClouds::~TestClouds() {
    //
    // dtor
    //
}

//______________________________________________________________________________
 void TestClouds::Configure() {
    //
    // read config files
    //
    
    // clouds features
    fType = Conf()->GetStr("TestClouds.fType");
    ConfigFileParser* pConfig3 = Config::Get()->GetCF("Atmosphere","AtmosphereFactory");
    fFeat = pConfig3->GetStr("Clouds.fFeat");    
}

//______________________________________________________________________________
 void TestClouds::Build() {
    //
    // build the object according to config files
    //
    
    TRandom* rndm = EsafRandom::Get();
    Double_t ran = 0;
    
    fG = Conf()->GetNum("TestClouds.fG");  //TOFIX : should depend of clouds altitude

    // set by hand
    if(fFeat == "perso") {
        fOD          = Conf()->GetNum("TestClouds.fOD");
        fTopAltitude = Conf()->GetNum("TestClouds.fTopAltitude")*km;
        fThickness   = Conf()->GetNum("TestClouds.fThickness")*km;
    }


    // set according redbook
    Bool_t CS = false;
    if(fFeat == "redbook" || fFeat == "redbook_noCS") {
        // if clearsky allowed, CS set to true
	if(fFeat == "redbook") CS = true;
	do {
	    ran = rndm->Rndm();
	    if(ran <= 0.041) {
	        CS = true;
        	fOD = ran/0.041 * 0.1;
		ran = rndm->Rndm() * 0.041;
		if(ran <= 0.013)       fTopAltitude = ran/0.013 * 3*km;
		else if(ran <= 0.02)   fTopAltitude = (ran - 0.013)/(0.02 - 0.013) * (7 - 3)*km + 3*km;
		else if(ran <= 0.025)  fTopAltitude = (ran - 0.02)/(0.025 - 0.02) * (10 - 7)*km + 7*km;
		else                   fTopAltitude = (ran - 0.025)/(0.041 - 0.025) * (12 - 10)*km + 10*km;   // 12 (Liou p.180)
	    }
	    else if (ran <= 0.334) {
	        CS = true;
		fOD = (ran - 0.041)/(0.334 - 0.041) * (1 - 0.1) + 0.1;
		ran = rndm->Rndm() * 0.293;
		if(ran <= 0.056)       fTopAltitude = ran/0.056 * 3*km;
		else if(ran <= 0.078)  fTopAltitude = (ran - 0.056)/(0.078 - 0.056) * (7 - 3)*km + 3*km;
		else if(ran <= 0.114)  fTopAltitude = (ran - 0.078)/(0.114 - 0.078) * (10 - 7)*km + 7*km;
		else                   fTopAltitude = (ran - 0.114)/(0.293 - 0.114) * (12 - 10)*km + 10*km;   // 12 (Liou p.180)
	    }
	    else if(ran <= 0.479) {
	        CS = true;
		fOD = (ran - 0.334)/(0.479 - 0.334) * (2 - 1) + 1;
		ran = rndm->Rndm() * 0.145;
		if(ran <= 0.051)       fTopAltitude = ran/0.051 * 3*km;
		else if(ran <= 0.08)   fTopAltitude = (ran - 0.051)/(0.08 - 0.051) * (7 - 3)*km + 3*km;
		else if(ran <= 0.115)  fTopAltitude = (ran - 0.08)/(0.115 - 0.08) * (10 - 7)*km + 7*km;
		else                   fTopAltitude = (ran - 0.115)/(0.145 - 0.115) * (12 - 10)*km + 10*km;   // 12 (Liou p.180)

	    }
	    else if(ran <= 0.822) {
	        CS = true;
		fOD = (ran - 0.479)/(0.822 - 0.479) * (3 - 2) + 2;   // 3 is arbitrary <=> 95% absorption
		ran = rndm->Rndm() * 0.343;
		if(ran <= 0.161)       fTopAltitude = ran/0.161 * 3*km;
		else if(ran <= 0.219)  fTopAltitude = (ran - 0.161)/(0.219 - 0.161) * (7 - 3)*km + 3*km;
		else if(ran <= 0.282)  fTopAltitude = (ran - 0.219)/(0.282 - 0.219) * (10 - 7)*km + 7*km;
		else                   fTopAltitude = (ran - 0.282)/(0.343 - 0.282) * (12 - 10)*km + 10*km;   // 12 (Liou p.180)
	    }
	    else {
		fOD = 0; // means no clouds
		fTopAltitude = 0;
	    }

	    // thickness according to Liou (rough estimation)
	    ran = rndm->Rndm();
	    if(fTopAltitude == 0) fThickness = 0;
	    else if(fTopAltitude <= 3*km)  fThickness = ran * 2*km;
	    else if(fTopAltitude <= 7*km)  fThickness = ran * 4*km;
	    else                           fThickness = 2*km;
	} while(!CS);

    }
    
    // checks
    ConfigFileParser* pConfig1 = Config::Get()->GetCF("RadiativeTransfer","TestGround");
    ConfigFileParser* pConfig2 = Config::Get()->GetCF("RadiativeTransfer","RadiativeFactory");
    if(pConfig2->GetStr("RadiativeFactory.Ground") != "test") 
        Msg(EsafMsg::Warning) << "if TESTGROUND not used, be carefull of ground altitude w.r.t clouds" << MsgDispatch;
    Double_t altg = pConfig1->GetNum("TestGround.fAltitude")*km;
    if((fTopAltitude - fThickness) < altg) fThickness = fTopAltitude - altg - 10*m; // 10m to avoid pb at interface
    if(fThickness < 0) {
        // means no clouds
        Msg(EsafMsg::Info) << "Clouds and Ground are too close -> No clouds" << MsgDispatch;
	fOD = 0;
	fTopAltitude = altg - TOLERANCE;
	fThickness = 0;
    }
    
    // red book formula for clouds albedo as a function of OD
    fAlbedo = 0.11*fOD / (1 + 0.11*fOD);
    
    Msg(EsafMsg::Info) << "*** TestClouds built ***   OD = " << fOD <<" --  TopAlt = " << fTopAltitude <<" --  Thick = "<< fThickness << MsgDispatch;
}

//______________________________________________________________________________
 Bool_t TestClouds::IsInClouds(const EarthVector& pos) const {
    //
    // true if position is within clouds, false otherwise
    //
    Bool_t rtn = false;
    if(!fOD) return rtn;
    if(pos.IsUnderSeaLevel()) return false;
    if( (pos.Zv() <= (fTopAltitude + kTolerance)) && (pos.Zv() >= (fTopAltitude - fThickness - kTolerance)) ) rtn = true;
    
    return rtn;
}

//______________________________________________________________________________
 Double_t TestClouds::TotalOD(const EarthVector& pos,const EarthVector& dir) const {
    //
    // Returns the total OD along a track defined by pos and dir
    // Case with pos within clouds handled
    // NB : vertical OD is homogeneously distributed
    //
    if(!fOD) return 0;
    EarthVector exit, enter, direc;
    Double_t rtn(0), mag(0);
    direc = dir.Unit();
    enter = GetCloudImpact(pos,direc);
    
    // if no clouds along the track
    if(enter.Z() == HUGE) return 0;
    
    exit = GetCloudImpact(pos,direc,"outgoing");
    if(exit.Z() == HUGE) Msg(EsafMsg::Warning) << "if enter exists, exit should as well" << MsgDispatch;
    mag = (exit - enter).Mag();
    rtn = fOD * mag/fThickness;
     
    return rtn;
}

//______________________________________________________________________________
 EarthVector TestClouds::RandomScatPos(const EarthVector& pos,const EarthVector& dir) const {
    //
    // returns a scattering position along a given track
    // cloud is uniform so scattering length is independent of direction : beta is constant
    // Trans = exp(-OD) = exp(-beta*l)
    //
    if(!fOD) return EarthVector(0,0,HUGE);
    EarthVector rtn;
    EarthVector direc = dir.Unit();

    Double_t beta = fOD/fThickness;
    rtn = pos - (log(1 - EsafRandom::Get()->Rndm()) / beta)*direc;
    
    return rtn;
}

//______________________________________________________________________________
 Double_t TestClouds::PhaseFunction(const EarthVector& incom_dir,const EarthVector& outgo_dir) const {
    //
    // Mie scattering phase function - Henyey-Greenstein parametrization
    // Normalized when integrated over full solid angle
    //
    
    if(!fOD) return 0.;
    Double_t rtn;
    Double_t theta = outgo_dir.Angle(incom_dir);
    rtn = 1./(4*Pi()) * (1 - fG*fG) / pow(1 - 2*fG*cos(theta) + fG*fG,1.5);

    return rtn;
}

//______________________________________________________________________________
 void TestClouds::RandomDir(EarthVector& dir) const {
    //
    // 'dir' is the incoming direction -> this method set dir to the scattering outgoing direction
    // get a random direction, sampling clouds scattering phase function
    //
    
    if(!fOD) return;
    
    Double_t theta(0);
    Double_t u1, u2;
    TRandom* rndm = EsafRandom::Get();
    
    // get random theta, assuming isotropic distribution in phi
    // theta is the angle between incoming and outgoing directions
    // von neumann sampling -> much faster than using TF1
    do {
        u1 = -1 + rndm->Rndm()*2;
        u2 = rndm->Rndm() * 3./4.;
    } while(u2 >= (1./(4*Pi()) * (1 - fG*fG) / pow(1 - 2*fG*u1 + fG*fG,1.5)));
    
    theta = acos(u1);
    
    // then get a random direction
    EarthVector v = dir.Orthogonal();
    v.Rotate(rndm->Rndm()*2*Pi(),dir);
    EarthVector axis = v.Cross(dir);
    dir.Rotate(theta,axis);
}

//_____________________________________________________________________________________________________
 EarthVector TestClouds::GetCloudImpact(const EarthVector& pos, const EarthVector& dir, string opt) const {
    //
    // Impact on clouds top surface form a given (position,direction) pair in atmosphere
    // if no impact, (0,0,HUGE) is returned
    // upward bunch is handled, but no scattering to the detector, only absorption
    //
    // if opt == "default"         -> incoming impact returned
    //                                if pos is within clouds, pos is returned
    // if opt == "outgoing"        -> outgoing impact returned
    //

    // init
    EarthVector rtn(0,0,HUGE);
    if(!fOD) return rtn;
    Bool_t isinclouds, downward;
    Double_t Altimpact(0);
    EarthVector direc = dir.Unit();
        
    // says if dir is going up/downward
    downward = GoingDownward(pos,direc);
    
    // says if pos is within the layer (altitude round-off effects < 1mm are handled)
    isinclouds = IsInClouds(pos);
    
    // set the impact surface altitude according to bunch present situation (up/downward, above/below, enter/exit clouds)
    if(isinclouds) {
        if(opt == "default") return pos;
	else if(opt == "outgoing") {
	    if(!downward) Altimpact = fTopAltitude;
	    else Altimpact = (fTopAltitude - fThickness);
	}
	else Msg(EsafMsg::Panic) << "<GetCloudImpact()> wrong argument (internal pos)" << MsgDispatch;
    }
    else {
        // if not going toward clouds
	if(pos.Zv() < (fTopAltitude - fThickness) && downward) return rtn; //TOFIX : because of earth-sphericity, impact is possible in MScatt
	if(pos.Zv() > fTopAltitude && !downward) return rtn;
	
	if(opt == "outgoing" && !downward || opt == "default" && downward) Altimpact = fTopAltitude;
	else if(opt == "outgoing" && downward || opt == "default" && !downward) Altimpact = (fTopAltitude - fThickness);
	else Msg(EsafMsg::Panic) << "<GetCloudImpact()> wrong argument (external pos)" << MsgDispatch;
    }

    // now calculate impact
    Double_t mag(0);
    Double_t b = pos*direc + direc.Z()*EarthRadius();
    Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - Altimpact) - pow(Altimpact,2);
    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
    if(p.first == 0) {
#ifdef DEBUG
    //Msg(EsafMsg::Debug) <<"NO IMPACT ON Clouds" << MsgDispatch;
#endif
	return rtn;
    }
    if(p.first == 1) mag = p.second[0];
    else if(p.first ==2) {
	if((p.second[0]+kTolerance) * (p.second[1]+kTolerance) < 0) mag = max(p.second[0],p.second[1]);
	else mag = min(p.second[0],p.second[1]);
    }
    if(mag < -kTolerance) Msg(EsafMsg::Warning) << "PB in impact calculation GetCloudImpact(), [min,max] = ["<<mag<<", "<< max(p.second[0],p.second[1])<<"]" << MsgDispatch;
    rtn = pos + mag*direc;    

    return rtn;
}

//_____________________________________________________________________________
 Bool_t TestClouds::GoingDownward(const EarthVector& pos,const EarthVector& direc) const {
    //
    // to know if locally track is going up/downward
    // locally means at pos nadir
    // if perfectly horizontal (i.e. local angle of halfpi) returns false (i.e. upward)
    //
    if(!fOD) return false;
    Bool_t rtn;
    EarthVector dir = direc.Unit();
    // find local zenith angle
    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
    // returns accordingly
    if(angle < (TMath::PiOver2() - kTolerance)) rtn = false;
    else  rtn = true;
    
    return rtn;
}

//_____________________________________________________________________________
 void TestClouds::Reset() {
    //
    // reset clodus features
    //
    
    fOD = 0;
    fThickness = 0;
    fTopAltitude = 0;
    fAlbedo = 0;
    
    Build();
}


ClassImp(TestClouds)
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