Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

SinglePhotonPropagator - source file

// $Id: SinglePhotonPropagator.cc,v 1.15 2005/11/03 12:01:33 moreggia Exp $
// Author: S. Moreggia    2005/08/16

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: SinglePhotonPropagator                                                           *
 *  Package: <packagename>                                                   *
 *  Coordinator: <coordinator>                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// SinglePhotonPropagator
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "SinglePhotonPropagator.hh"
#include "SinglePhoton.hh"
#include "FastRadiativeProcessesCalculator.hh"
#include "Clouds.hh"
#include "Atmosphere.hh"
#include "EsafRandom.hh"
#include "Ground.hh"
#include "EConst.hh"

using namespace sou;
using namespace TMath;
using namespace EConst;

ClassImp(SinglePhotonPropagator)

//_____________________________________________________________________________
 SinglePhotonPropagator::SinglePhotonPropagator() : EsafMsgSource() {
    //
    // Constructor
    //
    Msg(EsafMsg::Panic) <<"This constructor SHOULD not be used : ground must be given in argument" << MsgDispatch;

}
//_____________________________________________________________________________
 SinglePhotonPropagator::SinglePhotonPropagator(const Ground* g, Double_t maxphfunc, Double_t maxtof) {
    //
    // Constructor which copy RT ground description
    //
    CopyGround(g);
    fCalcul = new FastRadiativeProcessesCalculator();
    fDetGeom = 0;
    fMaxPhaseFunction = maxphfunc;
    fTofCut = maxtof;
}

//_____________________________________________________________________________
 SinglePhotonPropagator::~SinglePhotonPropagator() {
    //
    // Destructor
    //
    SafeDelete(fCalcul);
}

//_____________________________________________________________________________
 void SinglePhotonPropagator::Go(const vector<SinglePhoton*>& list,Int_t scatorder) const {
    //
    // MC code for photon scattering simulation
    //

    size_t nbPh = list.size();
    Msg(EsafMsg::Info) <<"Nb of SinglePhoton to propagate = " << nbPh << MsgDispatch;
    // for future info dumping
    Bool_t next(0),subnext(0);
    Double_t nbTracked(0), fraction(0), left(0),right(10),subleft(0),subright(2.5);
    Msg(EsafMsg::Info) << "Scattering Simulation: 0%" << MsgFlush;
    
    // loop over the list
    SinglePhoton* s = 0;
    for(size_t i=0; i<nbPh; i++) {
        s = list[i];
	
	DoOnePhoton(*s,scatorder);
	nbTracked++;

	// dump info
        fraction = 100*nbTracked/Double_t(nbPh);
        if (left<fraction && fraction <=right) 
            next = kFALSE;
        else if (fraction>right) {
            next = kTRUE;
            left = right;
            right += 10;
        }
        if (subleft<fraction  && fraction <=subright) {
            subnext = kFALSE;
        }
        if (fraction > subright) {
            subnext = kTRUE;
            subleft = subright;
            subright +=2;
        }
        if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
        if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;	 
    }
    Msg(EsafMsg::Info) << " -done" << MsgDispatch;
}

//_____________________________________________________________________________
 void SinglePhotonPropagator::DoOnePhoton(SinglePhoton& s,Int_t scatorder) const {
    //
    // MC propagation : handle rayleigh, clouds and ground for fluo and cerenkov
    // Each scattering order simulation is independent from another
    //

    const Clouds* clouds = Atmosphere::Get()->GetClouds();
    const Atmosphere* atmo = Atmosphere::Get();
    Int_t nbInter = 0;
    Double_t TOF = 0.;
    Int_t invert_status(0);   // for Atmosphere::InvertGrammage() method
    EarthVector dir(1.);
    Int_t count_debug(0); // to stop simu if infinite loop
    Bool_t IsCloudInter = false;
    
    // init
    TRandom* rndm = EsafRandom::Get();
    EarthVector nextimpact(1), scatpos_rayl(0.,0,HUGE), scatpos_clouds(0,0,HUGE);
    Medium nextmedium = CLEARSKY;

    // propagation stops when :
    // - photon is absorbed
    // - photon is out of atmosphere
    // - relevant scattering order is reached
    while(true) {
        
	// if weird situations appear
	if(nbInter == scatorder) Msg(EsafMsg::Warning) <<"<DoOnePhoton> Loop should have stopped before"<<MsgDispatch;
        if(count_debug++ > 100) {
            Msg(EsafMsg::Warning) <<"<DoOnePhoton> infinite loop broken by hand "<<MsgDispatch;
            break;
        }
	
	// cut in time of flight defined here
	if(TOF > fTofCut) {
	    s.SetStatus(LostByCut);
#ifdef DEBUG
            //Msg(EsafMsg::Debug) <<"<DoOnePhoton> TOF cut reached : photon is lost"<<MsgDispatch;
#endif
	    break;
	}
		
	// if within clouds
	//     - single scattering albedo ~ 1  ->  no absorption by water droplets or ice crystals (cf. Liou)
	//     - scattering cross-section is ~ lambda independent in UV band (cf. Liou)
	// so far works if clouds is a SINGLE LAYER medium
	if(clouds->IsInClouds(s.Pos()) || (nextmedium == CLOUDY) ) {
	    // 1. get ground impact
	    nextmedium = GetNextImpact(s.Pos(),s.Dir(),nextimpact,true);
	    
	    // 2. get random clouds scattering position,
	    // check it is within clouds.        If not  ->  scatpos_clouds = (0,0,HUGE)
	    // check it is before ground impact. If not  ->  scatpos_clouds = (0,0,HUGE)
	    scatpos_clouds = clouds->RandomScatPos(s.Pos(),s.Dir());
	    if(!clouds->IsInClouds(scatpos_clouds))    scatpos_clouds.SetXYZ(0,0,HUGE);
	    if(fGround->IsUnderGround(scatpos_clouds)) scatpos_clouds.SetXYZ(0,0,HUGE);
	    
	    // 3. get random rayleigh scattering position
	    // check it is before ground impact. If not  ->  grammage_rayl = -1
	    Double_t grammage_rayl = fCalcul->RandomGrammage_RaylScat(s.Wl());
	    if(atmo->Grammage(s.Pos(),nextimpact) < grammage_rayl) grammage_rayl = -1;
	    
	    // 4. if ground reached before any scattering,
	    // position sets to impact and goes to next step
	    if((scatpos_clouds.Z() == HUGE) && (grammage_rayl == -1) ) {
	        TOF += (nextimpact - s.Pos()).Mag() / Clight();
		s.AddToPosTof(nextimpact - s.Pos());
		continue;
	    }
	    
	    // else photon interacts with clouds or air
	    s.AddInteraction();
	    nbInter++;
	    
	    // 5. compare the two scattering positions and choose the nearest
	    IsCloudInter = false;
	    if(scatpos_clouds.Z() == HUGE) IsCloudInter = false;
	    else if((atmo->Grammage(s.Pos(),scatpos_clouds) < grammage_rayl) || (grammage_rayl == -1)) IsCloudInter = true;
	    
	    // then the corresponding type of scattering is applied
	    
	    // 5.1 clouds scattering	    
	    if(IsCloudInter) {
	        TOF += (scatpos_clouds - s.Pos()).Mag() / Clight();
	        s.AddToPosTof(scatpos_clouds - s.Pos());
		s.AddHistory(CloudScat);
	        // 5.10 if photon has come to its order of scattering
	        // look if it points toward detector and stop its simu
		if(scatorder == nbInter) {
	            if(rndm->Rndm() < (clouds->PhaseFunction(s.Dir(),(EUSO() - s.Pos()).Unit()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
	        	s.SetStatus(CloudScat);
			s.SetDir(EUSO() - s.Pos());
		    }
		    else s.SetStatus(LostByCut);
		    break;
		}
		// 5.11 else chose a new direction, and goes to next step
		else {
		    dir = s.Dir();
		    clouds->RandomDir(dir);
		    s.SetDir(dir);
		    s.SetStatus(LostByCut);
		    nextmedium = CLOUDY;
		}
	    }
	    
	    // 5.2 or rayleigh scattering
	    else {
	        if(grammage_rayl > -1) invert_status = atmo->InvertGrammage(s.Pos(),s.Dir(),grammage_rayl,scatpos_rayl,fTofCut - TOF);
		else Msg(EsafMsg::Warning) <<"<DoOnePhoton> This situation SHOULD NOT occur"<<MsgDispatch;
	        if(invert_status > 2) Msg(EsafMsg::Warning) <<"<DoOnePhoton> Atmo::InvertGrammage() returns status = " <<invert_status<<MsgDispatch;
                if(invert_status == 1) {
	            s.SetStatus(LostByCut);
		    break;
	        }
	        TOF += (scatpos_rayl - s.Pos()).Mag() / Clight();
	        s.AddToPosTof(scatpos_rayl - s.Pos());
		s.AddHistory(RaylScat);
	        // 5.20 if photon has come to its order of scattering
	        // look if it points toward detector and stop its simu
		if(scatorder == nbInter) {
	            if(rndm->Rndm() < (fCalcul->RayleighPhaseFunction(s.Dir(),(EUSO() - s.Pos()).Unit()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
	        	s.SetStatus(RaylScat);
			s.SetDir(EUSO() - s.Pos());
		    }
		    else s.SetStatus(LostByCut);
		    break;
		}
		// 5.21 else chose a new direction, and goes to next step
		else {
		    dir = s.Dir();
		    fCalcul->RandomDir("rayleigh",dir);
		    s.SetDir(dir);
		    s.SetStatus(LostByCut);
		    nextmedium = CLEARSKY;
		}
	    }
	}
	
	
	
	
	// if clear sky
	else if(nextmedium == CLEARSKY) {
	    // 1. get next medium type, and impact position upon it
	    nextmedium = GetNextImpact(s.Pos(),s.Dir(),nextimpact);
	    
	    // 2. get random rayleigh scattering position, expressed in grammage
	    Double_t grammage_rayl = fCalcul->RandomGrammage_RaylScat(s.Wl());
	    
	    // 3. if photon does not rayleigh-interact before next impact
	    // position sets to impact and goes to next step
	    if(atmo->Grammage(s.Pos(),nextimpact) < grammage_rayl) {
	        TOF += (nextimpact - s.Pos()).Mag() / Clight();
	        s.AddToPosTof(nextimpact - s.Pos());
		continue;
	    }
	    
	    // 4. else photon is rayleigh scattered at corresponding position
	    s.AddInteraction();
	    nbInter++;
	    invert_status = atmo->InvertGrammage(s.Pos(),s.Dir(),grammage_rayl,scatpos_rayl,fTofCut - TOF);
	    if(invert_status > 1) Msg(EsafMsg::Warning) <<"<DoOnePhoton> Atmo::InvertGrammage() returns status = " <<invert_status<<MsgDispatch;
            if(invert_status == 1) {
	        s.SetStatus(LostByCut);
		break;
	    }
	    TOF += (scatpos_rayl - s.Pos()).Mag() / Clight();
	    s.AddToPosTof(scatpos_rayl - s.Pos());
	    s.AddHistory(RaylScat);
	    
	    // 5. if photon has come to its order of scattering
	    // look if it points toward detector and stop its simu
	    if(scatorder == nbInter) {
	        if(rndm->Rndm() < (fCalcul->RayleighPhaseFunction(s.Dir(),(EUSO() - s.Pos()).Unit()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
	            s.SetStatus(RaylScat);
		    s.SetDir(EUSO() - s.Pos());
		}
		else s.SetStatus(LostByCut);
		break;
	    }
	    
	    // 6. else chose a new direction, and goes to next step
	    else {
		dir = s.Dir();
		fCalcul->RandomDir("rayleigh",dir);
		s.SetDir(dir);
		s.SetStatus(LostByCut);
		nextmedium = CLEARSKY;
	    }
		
	}
	
	
	
	
	
	
	// if on ground
	else if(nextmedium == GROUND) {
	    s.AddInteraction();
	    s.AddHistory(Reflected);
	    nbInter++;
	    // 1. if not absorbed at ground
	    if(rndm->Rndm() < fGround->Albedo(s.Pos())) {
		
	        // 1.1. if photon has come to its order of scattering
	        // look if it points toward detector and stop its simu
	        if(scatorder == nbInter) {
		    if(rndm->Rndm() < (fGround->Outgoing_phase_function(s.Pos(),EUSO()) * EusoOmega(s.Pos()) / (fMaxPhaseFunction*fOmegaMax)) ) {
		        s.SetStatus(Reflected);
		        s.SetDir(EUSO() - s.Pos());
		    }
		    else s.SetStatus(LostByCut);
		    break;
		}
		
		// 1.2. else get a new direction and go to next step
		else {
		    fGround->RandomDir(s.Pos(),dir);
		    s.SetDir(dir);
		    s.SetStatus(LostByCut);
		    nextmedium = CLEARSKY;
		}
	    }
	    
	    // 2. else photon is absorbed at ground
	    else {
	        s.SetAbsorbed();
		s.SetStatus(ScatAbsorb);
		break;
	    }
	}
	
	
	else if(nextmedium == TOA) {
	    s.AddHistory(OutAtmo);
	    s.SetStatus(OutAtmo);
	    break;
	}
	
	else if(nextmedium == NONE) Msg(EsafMsg::Panic) <<"<DoOnePhoton> next medium is NONE -> SHOULD NOT" <<MsgDispatch;
    }
}

//_____________________________________________________________________________
 Medium SinglePhotonPropagator::GetNextImpact(const EarthVector& pos,const EarthVector& dir,EarthVector& res,Bool_t without_clouds) const {
    //
    // Calculate the impact upon the next medium according to the given track
    // Return type of the next encountered medium
    //
    // - If pos is underground, NONE medium is returned and impact set to (0,0,-HUGE)
    // OUT OF TOA is handled
    //
    
    Medium rtn = GROUND;
    EarthVector gd_imp, cloud_imp;
    
    // get impact on ground
    gd_imp = fGround->GetImpact(pos,dir);
    res = gd_imp;
    // if pos is underground
    if(gd_imp.Z() == -HUGE) return NONE;
  
    // get impact on clouds, then make appropriate changes
    cloud_imp = Atmosphere::Get()->GetClouds()->GetCloudImpact(pos,dir);
    if( ((cloud_imp - pos).Mag() < (gd_imp - pos).Mag()) && !without_clouds ) {
	res = cloud_imp;
	rtn = CLOUDY;
    }
    
    // if no impact, means top of atmosphere reached
    if(res.Z() == HUGE) {
        rtn = TOA;
	res = Atmosphere::Get()->ImpactAtTOA(pos,dir,50*km);  // (should be same as in atmosphere::InvertGrammage() method) TOA at 50km to avoid pbic cases at high altitudes
    }
      
    return rtn;
}

//_____________________________________________________________________________
 void SinglePhotonPropagator::Reset() {
    //
    // ready for next event
    //
    
    fCalcul->Reset();
}



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