source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/radiativetransfer/src/O2_ClearSkyPropagator.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: O2_ClearSkyPropagator.cc 2767 2006-11-15 15:50:14Z moreggia $
// Sylvain Moreggia created Jun,  2 2004

#include "O2_ClearSkyPropagator.hh"
#include "SinglePhoton.hh"
#include "BunchOfPhotons.hh"
#include <math.h>
#include "EsafRandom.hh"
#include "ListPhotonsInAtmosphere.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "RadiativeFactory.hh"
#include "EConst.hh"
#include "Atmosphere.hh"
#include "Clouds.hh"
#include "LowtranRadiativeProcessesCalculator.hh"
#include "TestGround.hh"
#include "TBenchmark.h"

ClassImp(O2_ClearSkyPropagator)

using namespace sou;
using namespace EConst;

//_____________________________________________________________________________
O2_ClearSkyPropagator::O2_ClearSkyPropagator() : ClearSkyPropagator() {
    //
    // ctor (should not be used)
    //
    fOrder = 2;
    fFinal_medium = NONE;
}

//_____________________________________________________________________________
O2_ClearSkyPropagator::O2_ClearSkyPropagator(const Ground* g) : ClearSkyPropagator(g) {
    //
    // ctor, copy RadiatvieTransfer ground description
    //
    fOrder = 2;
    fFinal_medium = NONE;
    fNbSinglesCorrected = 0;
    
    ConfigFileParser *pConfig = Config::Get()->GetCF("Atmosphere", "LowtranAtmosphere");
    fAerosolModel = (Int_t)pConfig->GetNum("LowtranAtmosphere.Ihaze");
    
    // if model is outside (0,1,2,4,5), it is not used for scattering simulation : ratio aborption / scattering is currently not known by me
    if(fAerosolModel != 0 && fAerosolModel != 1 && fAerosolModel != 2 && fAerosolModel != 4 && fAerosolModel != 5) {
        fAerosolModel = 0;
        Msg(EsafMsg::Warning) <<"Aerosol model is such that scattering simu cannot be done using BUNCH algorithm" << MsgDispatch;
        Msg(EsafMsg::Warning) <<"Aerosol model will be used only for last transmission to detector" << MsgDispatch;
    }
    // Aerosol case : absorption / extinction ratio treated in a simple way (cause detailled treatment not supported by the current interface with LOWTRAN)
    // rural    : 95 % scattering, ~independent of wavelength
    // urban    : 70 % scattering, ~independent of wavelength
    // maritime : 98 % scattering, ~independent of wavelength
    if(fAerosolModel == 0)      fAerosol_correction_factor = 0.;
    else if(fAerosolModel == 1) fAerosol_correction_factor = 0.95;
    else if(fAerosolModel == 2) fAerosol_correction_factor = 0.95;
    else if(fAerosolModel == 4) fAerosol_correction_factor = 0.98;
    else if(fAerosolModel == 5) fAerosol_correction_factor = 0.70;
}

//_____________________________________________________________________________
O2_ClearSkyPropagator::~O2_ClearSkyPropagator() {
    //
    // dtor
    //
}

//_____________________________________________________________________________
Medium O2_ClearSkyPropagator::Go(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
    //
    // call methods for propagation
    // NB: bunch can be modified, its copy b cannot
    //
    
    // to assess backscattering algorithm time consumption
    #ifdef DEBUG
    char prt[100];
    sprintf(prt,"Bunch order2 propagation");
    TBenchmark bench;
    bench.Start(prt);
    #endif
    
    // initializations
    static UInt_t bunch_id = 0; // to know if bunch has already come here (for clouds etc..)
    Medium rtn = NONE;
    const BunchOfPhotons& b = bunch;
    EarthVector impact;
    
    // if fluo, not propagated  //TOFIX: fluo propagation should be handled in the future
    if(b.GetType() == Fluo) {
        bunch.SetFate(1);
        return NONE;
    }
        
    // if first treatment for this bunch,
    // get next impact and next medium encountered
    //TOFIX : valid algorithm if a SINGLE intermediate medium exists
    if(bunch_id != b.GetId()) {
        // bunch preprocessing -> generate a lowtran table along the bunch track
        BunchPreProcess(bunch);
        if(fFinal_medium == NONE) return NONE;
        bunch_id = b.GetId();
        rtn = GetNextImpact(b,impact);
        bunch.SetNextImpact(impact);
        if(impact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
        if(b.GetFinalImpact().Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
        // if no intermediate medium, go to final impact
        if((b.GetNextImpact() - b.GetFinalImpact()).Mag() == 0) rtn = fFinal_medium;
    }
    // otherwise go to final impact
    else {
        bunch.SetNextImpact(b.GetFinalImpact());
        rtn = fFinal_medium;
    }
    
    // if bunch showerpos is within clouds, switch off the present clear sky propagation
    if(rtn == CLOUDY && (b.GetNextImpact() - b.GetPos()).Mag() < TOLERANCE) {
        return rtn;
    }

    // run propagation algorithms
    Propagate(bunch,list);

    #ifdef DEBUG
    bench.Stop(prt);
    MsgForm(EsafMsg::Debug,"REAL=%6.2f s   CPU=%6.2f s",bench.GetRealTime(prt),bench.GetCpuTime(prt));  
    #endif
    
    if(rtn == GROUND) bunch.SetFate(4);
    // in the case impact is out of FoV, rtn set to NONE to stop the propagation (no ground reflexion)
    if(rtn == CLEARSKY) {
        rtn = NONE;
        bunch.SetFate(5);
    }

    return rtn;
}

//_____________________________________________________________________________
void O2_ClearSkyPropagator::Propagate(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
    //
    // Propagate a BunchOfPhotons in clear sky conditions,
    // generating SinglePhotons all along the travel -> ListPhotonsInAtmosphere filled
    // NB: bunch can be modified, its copy b cannot
    //
    
    Int_t status = -100;  // for Atmosphere::InvertGrammage() method
    const BunchOfPhotons& b = bunch;
    const EarthVector& impact = b.GetNextImpact();
    const Atmosphere* atmo = Atmosphere::Get();
    // initializations
    Double_t conv = 1.;
    Int_t nb = b.GetWlSpectrum().GetAxis().GetNbins();
    vector<Double_t*> coeff;
    Double_t* TotalTrans = new Double_t[nb];
    Double_t* RayleighTrans = new Double_t[nb];
    Double_t* OzoneTrans = new Double_t[nb];
    Double_t* AerosolTrans = new Double_t[nb];
    if(!TotalTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for TotalTrans" << MsgDispatch;
    if(!RayleighTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for RayleighTrans" << MsgDispatch;
    if(!OzoneTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for OzoneTrans" << MsgDispatch;
    if(!AerosolTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for AerosolTrans" << MsgDispatch;
    coeff.push_back(TotalTrans);
    coeff.push_back(RayleighTrans);
    coeff.push_back(OzoneTrans);
    coeff.push_back(AerosolTrans);
    Double_t raylScat = 0;
    Double_t phfunc = 0;
    Int_t NbRaylBackscat = 0;
    EarthVector pos, step, entry, exit, dir;
    Double_t depthstep = Conf()->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
    
    entry = b.GetPos() - 0.5*(b.GetParent()->GetShowerPosf() - b.GetParent()->GetShowerPosi());
    dir = (impact - entry).Unit();
    status = atmo->InvertGrammage(entry,dir,depthstep,exit);
    
#ifdef DEBUG
    // to get nb of iteration for propagation of this bunch
    Float_t count = 0;
#endif
    
    while((impact - b.GetPos()).Mag() > TOLERANCE) {
#ifdef DEBUG
        count++;
#endif
        step = exit - entry;
        if((impact - b.GetPos()).Mag() < step.Mag()) step = impact - b.GetPos();
        fCalcul->Trans(b,b.GetPos()+step,coeff);
        phfunc = fCalcul->RayleighPhaseFunction(b.GetDir(),EUSO() - b.GetPos() - 0.5*step);
        
        // position and time of flight incrementation (bunch position is at the step end)
        bunch.AddToPosTof(step);
        
        // cerenkov backscattering (only rayleigh so far) //TOFIX
        for(Int_t i=0; i<nb; i++) {
            raylScat = (1 - coeff[1][i]);
            NbRaylBackscat  = EsafRandom::Get()->Poisson(b.GetWeight(i) * EusoOmega(b.GetPos() - 0.5*step) * raylScat*phfunc);
            // Singles generation (bunch pos at step end) -- aerosols-ozone absorption also taken into account there
            if(NbRaylBackscat) GenerateRayleighSingles(i,NbRaylBackscat,raylScat,step,b,list,coeff[3][i],coeff[2][i]);
        }
        // convolutions with transmission coefficients
        conv = bunch.ConvoluteWl(coeff[0]);
        bunch.Convolute(conv);
        entry = exit;
        status = atmo->InvertGrammage(entry,dir,depthstep,exit);
    }
    // to have position precisely, for later processes
    bunch.SetPos(impact);
    
#ifdef DEBUG
    Msg(EsafMsg::Debug) <<"Nb iteration = " <<count << MsgDispatch;
    Msg(EsafMsg::Debug) <<"after backscat, SinglePhoton list size = "<<list.GetListOfSingle().size() << MsgDispatch;
#endif

    // cleaning of coeff
    for(UInt_t m=0; m<coeff.size(); m++) {
        if(coeff[m]) delete[] coeff[m];
        coeff[m] = 0;
    }
    coeff.clear();
}

//_____________________________________________________________________________
void O2_ClearSkyPropagator::GenerateRayleighSingles(Int_t specbin, Int_t nb, Double_t scattrate, const EarthVector& step,
                                            const BunchOfPhotons& b, ListPhotonsInAtmosphere& list, Double_t Ta,Double_t To) const {
    //
    // Generate a SinglePhoton list, coming from Rayleigh scattering process applied to the given BunchOfPhotons
    // Tuning of scattering position done thanks to additional data : scattrate (rate of scattering),
    // step (step used for backscattering calculations)
    // NB : bunch position is at the end of the considered propagation step
    //

    Double_t wl, date, tof, scatlength;
    EarthVector showerpos, pos, dir, diff;
    SinglePhoton* s = 0;
    TRandom* rndm = EsafRandom::Get();
    
    UInt_t bid = b.GetId();
    PhotonType type = b.GetType();
    
    for(Int_t i=0; i<nb; i++) {
        // corrections for date and showerpos, from BunchOfPhotons mean values, longitudinal and lateral distributions at creation
        showerpos = b.RandomPosInShower();
        if(fGround->IsUnderGround(showerpos)) continue;
        diff = showerpos - b.GetShowerPos();
        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
        
        // calculation of the scattering position within the step
        pos = b.GetPos() - step;    // as if bunch was at the beginning of the step
        scatlength = (log(1 - rndm->Rndm()*scattrate)/log(1 - scattrate));
        pos += scatlength *step;
        
        // tell if this photon has been really transmitted until the rayleigh scat position (ozone - aerosols)
        // if not  ->  photon is killed
        if(rndm->Rndm() > exp(scatlength*log(Ta))) continue;
        if(rndm->Rndm() > exp(scatlength*log(To))) continue;
        
        // position correction, due to angular distributions at creation
        pos = b.PosRandomAngCorrec(showerpos,pos);
        
        // if underground, photon is kept for another try
        if(fGround->IsUnderGround(pos)) {
            EarthVector axis, impact;
            impact = b.GetNextImpact();
            if(impact.Z() == HUGE) {
                Msg(EsafMsg::Warning) <<"\nno bunch impact, but singles underground" <<MsgDispatch;
                continue;
            }
            else {
                EarthVector localvertical(impact);
                localvertical(2) += EarthRadius();
                axis = localvertical.Cross(b.GetDir());
                if(axis.Mag()) {    // if localvertical and bunch dir not colinear
                    EarthVector v1 = pos - impact;
                    v1.Rotate(TMath::Pi(),axis);
                    pos = v1 + impact;
                }
                else Msg(EsafMsg::Warning) <<"<GenerateSingles> OK if vertical track" <<MsgDispatch;
                if(fGround->IsUnderGround(pos)) {
                    Msg(EsafMsg::Warning) <<"GenerateSingles() -> pos is still underground" <<MsgDispatch;
                    continue;
                }
                fNbSinglesCorrected++;
            }
        }
        tof = b.GetTof() * (showerpos - pos).Mag()/(b.GetPos() - b.GetShowerPos()).Mag();
        dir = EUSO() - pos;
        wl = b.GetWlSpectrum().GetLambda(specbin);
        s = new SinglePhoton(type,date,tof,wl,showerpos,pos,dir,RaylScat,bid,b.GetShowerAge());
        s->AddInteraction();
        s->AddHistory(RaylScat);
        list.Add(s);
    } 
}

//_____________________________________________________________________________
void O2_ClearSkyPropagator::GenerateAerosolSingles(Int_t specbin, Int_t nb, Double_t scattrate, const EarthVector& step,
                                            const BunchOfPhotons& b, ListPhotonsInAtmosphere& list, Double_t Tr,Double_t To) const {
    //
    // Generate a SinglePhoton list, coming from Aerosols Mie scattering process applied to the given BunchOfPhotons
    // Tuning of scattering position done thanks to additional data : scattrate (rate of scattering),
    // NB : bunch position is at the end of the considered propagation step
    //

    Double_t wl, date, tof, scatlength;
    EarthVector showerpos, pos, dir, diff;
    SinglePhoton* s = 0;
    TRandom* rndm = EsafRandom::Get();
    
    UInt_t bid = b.GetId();
    PhotonType type = b.GetType();
    
    for(Int_t i=0; i<nb; i++) {
        // corrections for date and showerpos, from BunchOfPhotons mean values, longitudinal and lateral distributions at creation
        showerpos = b.RandomPosInShower();
        if(fGround->IsUnderGround(showerpos)) continue;
        diff = showerpos - b.GetShowerPos();
        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
        
        // calculation of the scattering position within the step
        pos = b.GetPos() - step;    // as if bunch was at the beginning of the step
        scatlength = (log(1 - rndm->Rndm()*scattrate)/log(1 - scattrate));
        pos += scatlength *step;
        
        // tell if this photon has been really transmitted until the mie scat position (ozone - molecular)
        // if not  ->  photon is killed
        if(rndm->Rndm() > exp(scatlength*log(Tr))) continue;
        if(rndm->Rndm() > exp(scatlength*log(To))) continue;
        
        // position correction, due to angular distributions at creation
        pos = b.PosRandomAngCorrec(showerpos,pos);
        
        // if underground, photon is kept for another try
        if(fGround->IsUnderGround(pos)) {
            EarthVector axis, impact;
            impact = b.GetNextImpact();
            if(impact.Z() == HUGE) {
                Msg(EsafMsg::Warning) <<"\nno bunch impact, but singles underground" <<MsgDispatch;
                continue;
            }
            else {
                EarthVector localvertical(impact);
                localvertical(2) += EarthRadius();
                axis = localvertical.Cross(b.GetDir());
                if(axis.Mag()) {    // if localvertical and bunch dir not colinear
                    EarthVector v1 = pos - impact;
                    v1.Rotate(TMath::Pi(),axis);
                    pos = v1 + impact;
                }
                if(fGround->IsUnderGround(pos)) {
                    Msg(EsafMsg::Warning) <<"GenerateSingles() -> pos is still underground" <<MsgDispatch;
                    continue;
                }
            }
        }
        tof = b.GetTof() * (showerpos - pos).Mag()/(b.GetPos() - b.GetShowerPos()).Mag();
        dir = EUSO() - pos;
        wl = b.GetWlSpectrum().GetLambda(specbin);
        s = new SinglePhoton(type,date,tof,wl,showerpos,pos,dir,AeroScat,bid,b.GetShowerAge());
        s->AddInteraction();
        s->AddHistory(AeroScat);
        list.Add(s);
    } 
}

//_____________________________________________________________________________
Medium O2_ClearSkyPropagator::BunchPreProcess(BunchOfPhotons& b) const {
    //
    // Preprocess run once for each bunch before propagation
    // Call RadiativeProcessesCalculator::BunchPreProcess
    // Get the Medium of final impact
    //
#ifdef DEBUG
    Msg(EsafMsg::Debug) <<"IN O2_ClearSkyPropagator::bunchPreProcess" << MsgDispatch;
#endif
    const BunchOfPhotons& bunch = b;
    EarthVector finalimpact;
    
    // assess the next impact
    fFinal_medium = GetFinalImpact(bunch,finalimpact);
    b.SetFinalImpact(finalimpact);
    if(finalimpact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
   
    // if bunch is underground, its simulation is stopped now
    if(fFinal_medium == NONE) {
        b.SetFate(3);
        return fFinal_medium;
    }
    
    // run Calculator BunchPreProcessing
    if(fCalcul->GetName() == "lowtran") {
        LowtranRadiativeProcessesCalculator* calcul = (LowtranRadiativeProcessesCalculator*) fCalcul;
        calcul->PreProcess(bunch);
    }
    
    return fFinal_medium;
}

//_____________________________________________________________________________
void O2_ClearSkyPropagator::Reset() {
    //
    // reset imbricated objects
    //
    fCalcul->Reset();
    fFinal_medium = NONE;
    fNbSinglesCorrected = 0;
}