// $Id: AlongTrack_CSPropagator.cc,v 1.25 2006/03/10 13:53:27 moreggia Exp $
// Author: Sylvain Moreggia   2004/09/29

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: AlongTrack_CSPropagator                                              *
 *  Package: RadiativeTransfer                                               *
 *  Coordinator: Sylvain Moreggia                                            *
 *                                                                           *
 *****************************************************************************/

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

#include "AlongTrack_CSPropagator.hh"
#include "TBenchmark.h"
#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"

ClassImp(AlongTrack_CSPropagator)

//_____________________________________________________________________________
 AlongTrack_CSPropagator::AlongTrack_CSPropagator() : O2_ClearSkyPropagator() {
    //
    // Constructor (should not be used)
    //
    fOrder = 2;
    fFlag = false;
    fFinal_medium = CLEARSKY;
    fFirstBunch = 0;
}

//_____________________________________________________________________________
 AlongTrack_CSPropagator::AlongTrack_CSPropagator(const Ground* g) : O2_ClearSkyPropagator(g) {
    //
    // Constructor, copy RadiatvieTransfer ground description
    //
    fOrder = 2;
    fFlag = false;
    fFinal_medium = CLEARSKY;
    fFirstBunch = 0;
}

//_____________________________________________________________________________
 AlongTrack_CSPropagator::~AlongTrack_CSPropagator() {
    //
    // Destructor
    //
    SafeDelete(fFirstBunch);
}

//_____________________________________________________________________________
 Medium AlongTrack_CSPropagator::Go(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
    //
    // call preprocessing, then propagation method
    //
    
    // to assess backscattering algorithm time consumption
    #ifdef DEBUG
    char prt[100];
    sprintf(prt,"Bunch ALONG TRACK 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..)
    const BunchOfPhotons& b = bunch;
    EarthVector impact;
    Medium rtn = NONE;
    
    // if fluo, not propagated  //TOFIX: fluo propagation should be handled in the future
    if(b.GetType() == Fluo) {
	bunch.SetFate(1);
        return NONE;
    }

    // run preprocess -> generation of lowtran table along the track
    if(!fFlag) PreProcess(b);
    // if PreProcess failed because track beginning is under ground
    if(!fFlag) {
        bunch.SetFate(3);
	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_id = b.GetId();
	bunch.SetFinalImpact(fFirstBunch->GetFinalImpact());
        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;
    }
    
    // test if propagator is really dealing with a track
    if(!IsAlongTrack(b)) Msg(EsafMsg::Warning) << "Go() : AlongTrack used, but it does not seem to be a track" << MsgDispatch;

    // if bunch is under ground, propagation is stopped
    if((b.GetFinalImpact() - fFirstBunch->GetPos()).Mag() < (b.GetPos() - fFirstBunch->GetPos()).Mag()) {
        bunch.SetFate(3);
	return NONE;
    }
    
    // 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
    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 the impact is out of FoV, NONE is returned to stop the propagation (no ground reflexion)
    if(rtn == CLEARSKY) {
        rtn = NONE;
	bunch.SetFate(5);
    }

    return rtn;
}
 
//_____________________________________________________________________________
 void AlongTrack_CSPropagator::Propagate(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
    //
    // Propagate along a track a BunchOfPhotons in clear sky conditions,
    // generating SinglePhotons all along the travel -> ListPhotonsInAtmosphere filled
    // NB: bunch can be modified, its copy b cannot
    //
    
#ifdef DEBUG
    Msg(EsafMsg::Debug) <<"Propagate()" << MsgDispatch;
#endif

    // initializations
    const BunchOfPhotons& b = bunch;
    const EarthVector& impact = b.GetNextImpact();
    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),aerScat(0),phfunc_rayl(0),phfunc_aer(0);
    Int_t NbRaylBackscat = 0;
    Int_t Nb_Aer_Backscat = 0;
    EarthVector pos, step, entry, exit;
    if(list.GetNbTrackSteps() == 0) Msg(EsafMsg::Panic) << "this propagator needs ListPhotonsInAtmosphere::fTrack" << MsgDispatch;
    
    // find the matching track step
    UInt_t index = list.TrackStepIndex(b.GetPos());
    entry = b.GetPos();
    if(index < (list.GetNbTrackSteps()-1)) exit = list.GetTrackStep(++index);
    else {
         Msg(EsafMsg::Warning) << "pb with array index : a bunch position is OUTSIDE light track"<<MsgDispatch;
	 exit = impact;
#ifdef DEBUG
Msg(EsafMsg::Debug) << "POSITION OF BUNCH WHEN PROBLEM = " <<b.GetPos() <<MsgDispatch;
Msg(EsafMsg::Debug) << "POSITION OF BUNCH WHEN PROBLEM = " <<(b.GetPos() - list.GetTrackStep(0)).Mag()*1e-6 <<MsgDispatch;
Msg(EsafMsg::Debug) << "ALTITUDE OF BUNCH WHEN PROBLEM = " <<b.GetPos().Zv() <<MsgDispatch;
#endif
    }
    
#ifdef DEBUG
    // to get nb of iteration for propagation of this bunch
    Float_t count = 0;
#endif

    // loop of step by step propagation for backscattering simulation
    while((impact - b.GetPos()).Dot(impact - b.GetShowerPos()) > TOLERANCE) {
#ifdef DEBUG
        count++;
#endif
	step = exit - entry;
	if((impact - (b.GetPos()+step)).Dot(impact - b.GetPos()) < TOLERANCE) step = impact - b.GetPos();
	fCalcul->Trans(b,b.GetPos()+step,coeff);
	phfunc_rayl = fCalcul->RayleighPhaseFunction(b.GetDir(),EUSO() - b.GetPos() - 0.5*step);
	phfunc_aer = fCalcul->Mie_HG_PhaseFunction(b.GetDir(),EUSO() - b.GetPos() - 0.5*step,0.7);
	
	// position and time of flight incrementation (bunch position is at the step end)
	bunch.AddToPosTof(step);
	bunch.SetPos(exit); // often useless, but SAFER : handle misaligned bunches
	
	// cerenkov backscattering (rayleigh-mie)
	for(Int_t i=0; i<nb; i++) {
	    raylScat = (1 - coeff[1][i]);
	    aerScat  = (1 - coeff[3][i]);
	    NbRaylBackscat  = EsafRandom::Get()->Poisson(b.GetWeight(i) * EusoOmega(b.GetPos() - 0.5*step) * raylScat*phfunc_rayl);
	    Nb_Aer_Backscat = EsafRandom::Get()->Poisson(b.GetWeight(i) * EusoOmega(b.GetPos() - 0.5*step) * aerScat*phfunc_aer);
	    // 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]);
	    if(Nb_Aer_Backscat) GenerateAerosolSingles(i,Nb_Aer_Backscat,aerScat,step,b,list,coeff[1][i],coeff[2][i]);
	}
	// convolutions with transmission coefficients
	conv = bunch.ConvoluteWl(coeff[0]);
	bunch.Convolute(conv);
	entry = exit;
	if(index != list.GetNbTrackSteps()-1) {
	    index++;
	    exit = list.GetTrackStep(index);
        }
	else exit = impact;
    }
    // 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;
    Msg(EsafMsg::Debug) <<"Nb of Singles CORRECTED = "<<fNbSinglesCorrected << 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 AlongTrack_CSPropagator::PreProcess(const BunchOfPhotons& b) const {
    //
    // Set fFlag and fImpact, Call RadiativeProcessesCalculator::PreProcess()
    // Called once (for the first bunch) and used for propagation of all the bunches
    //
    
#ifdef DEBUG
    Msg(EsafMsg::Debug) <<"PreProcess()" << MsgDispatch;
#endif

    // get final impact
    EarthVector finalImpact;
    fFinal_medium = GetFinalImpact(b,finalImpact);
    Msg(EsafMsg::Info) <<"PreProcess() : FINAL MEDIUM FOR track IMPACT = "<< fFinal_medium << MsgDispatch;
    Msg(EsafMsg::Info) <<"PreProcess() : FINAL track IMPACT = "<< finalImpact << MsgDispatch;
    fFirstBunch = new BunchOfPhotons(b);
    fFirstBunch->SetFinalImpact(finalImpact);
    if(finalImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an "impact"" << MsgDispatch;
    
    if(fFinal_medium != NONE) {
        // create lowtran trans table along the whole track (until final track impact)
        if(fCalcul->GetName() == "lowtran") {
	    LowtranRadiativeProcessesCalculator* calcul = (LowtranRadiativeProcessesCalculator*) fCalcul;
	    calcul->PreProcess(*fFirstBunch);
	}
        fFlag = true;
    }
    // if first bunch is under ground, it is not used to initialize the propagator. The next bunch will be
    else {
        Msg(EsafMsg::Info) <<"PreProcess() : track starts underground" << MsgDispatch;
        SafeDelete(fFirstBunch);
    }
}

//_____________________________________________________________________________
 Bool_t AlongTrack_CSPropagator::IsAlongTrack(const BunchOfPhotons& b) const {
    //
    // Check if bunches are along a track, with first at the track beginning
    //
    if(b.GetId() == fFirstBunch->GetId()) return true;
    Bool_t rtn = false;
    EarthVector diff = (b.GetPos() - fFirstBunch->GetPos()).Unit();
    
    // check if bunch is on the track (defined by firstbunch pos and dir)
    // and check if firstbunch is really the first one
    if((diff - fFirstBunch->GetDir()).Mag() < 1.e5*TOLERANCE) rtn = true;
    else {
#ifdef DEBUG
        Msg(EsafMsg::Debug) << "Bunch nb "<<b.GetId()<<" not perfectly along track -> "<<(diff - fFirstBunch->GetDir()).Mag() <<MsgDispatch;
#endif
    }
   
    return rtn;
}

//_____________________________________________________________________________
 void AlongTrack_CSPropagator::Reset() {
    //
    // Reset method
    //
    fCalcul->Reset();
    fFlag = false;
    fFinal_medium = CLEARSKY;
    SafeDelete(fFirstBunch);
    fNbSinglesCorrected = 0;
}