source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/radiativetransfer/src/AlongTrack_CSPropagator.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 12.4 KB
Line 
1// $Id: AlongTrack_CSPropagator.cc 2751 2006-07-13 12:43:10Z moreggia $
2// Author: Sylvain Moreggia   2004/09/29
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: AlongTrack_CSPropagator                                              *
8 *  Package: RadiativeTransfer                                               *
9 *  Coordinator: Sylvain Moreggia                                            *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// AlongTrack_CSPropagator
16//
17// <extensive class description>
18//
19//   Config file parameters
20//   ======================
21//
22//   <parameter name>: <parameter description>
23//   -Valid options: <available options>
24//
25
26#include "AlongTrack_CSPropagator.hh"
27#include "TBenchmark.h"
28#include "SinglePhoton.hh"
29#include "BunchOfPhotons.hh"
30#include <math.h>
31#include "EsafRandom.hh"
32#include "ListPhotonsInAtmosphere.hh"
33#include "Config.hh"
34#include "EsafSpectrum.hh"
35#include "RadiativeFactory.hh"
36#include "EConst.hh"
37#include "Atmosphere.hh"
38#include "Clouds.hh"
39#include "LowtranRadiativeProcessesCalculator.hh"
40
41ClassImp(AlongTrack_CSPropagator)
42
43//_____________________________________________________________________________
44AlongTrack_CSPropagator::AlongTrack_CSPropagator() : O2_ClearSkyPropagator() {
45    //
46    // Constructor (should not be used)
47    //
48    fOrder = 2;
49    fFlag = false;
50    fFinal_medium = CLEARSKY;
51    fFirstBunch = 0;
52}
53
54//_____________________________________________________________________________
55AlongTrack_CSPropagator::AlongTrack_CSPropagator(const Ground* g) : O2_ClearSkyPropagator(g) {
56    //
57    // Constructor, copy RadiatvieTransfer ground description
58    //
59    fOrder = 2;
60    fFlag = false;
61    fFinal_medium = CLEARSKY;
62    fFirstBunch = 0;
63}
64
65//_____________________________________________________________________________
66AlongTrack_CSPropagator::~AlongTrack_CSPropagator() {
67    //
68    // Destructor
69    //
70    SafeDelete(fFirstBunch);
71}
72
73//_____________________________________________________________________________
74Medium AlongTrack_CSPropagator::Go(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
75    //
76    // call preprocessing, then propagation method
77    //
78   
79    // to assess backscattering algorithm time consumption
80    #ifdef DEBUG
81    char prt[100];
82    sprintf(prt,"Bunch ALONG TRACK propagation");
83    TBenchmark bench;
84    bench.Start(prt);
85    #endif
86   
87    // initializations
88    static UInt_t bunch_id = 0; // to know if bunch has already come here (for clouds etc..)
89    const BunchOfPhotons& b = bunch;
90    EarthVector impact;
91    Medium rtn = NONE;
92   
93    // if fluo, not propagated  //TOFIX: fluo propagation should be handled in the future
94    if(b.GetType() == Fluo) {
95        bunch.SetFate(1);
96        return NONE;
97    }
98
99    // run preprocess -> generation of lowtran table along the track
100    if(!fFlag) PreProcess(b);
101    // if PreProcess failed because track beginning is under ground
102    if(!fFlag) {
103        bunch.SetFate(3);
104        return NONE;
105    }
106   
107    // if first treatment for this bunch,
108    // get next impact and next medium encountered
109    //TOFIX : valid algorithm if a SINGLE intermediate medium exists
110    if(bunch_id != b.GetId()) {
111        bunch_id = b.GetId();
112        bunch.SetFinalImpact(fFirstBunch->GetFinalImpact());
113        rtn = GetNextImpact(b,impact);
114        bunch.SetNextImpact(impact);
115        if(impact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
116        if(b.GetFinalImpact().Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
117        // if no intermediate medium, go to final impact
118        if((b.GetNextImpact() - b.GetFinalImpact()).Mag() == 0) rtn = fFinal_medium;
119    }
120    // otherwise go to final impact
121    else {
122        bunch.SetNextImpact(b.GetFinalImpact());
123        rtn = fFinal_medium;
124    }
125   
126    // test if propagator is really dealing with a track
127    if(!IsAlongTrack(b)) Msg(EsafMsg::Warning) << "Go() : AlongTrack used, but it does not seem to be a track" << MsgDispatch;
128
129    // if bunch is under ground, propagation is stopped
130    if((b.GetFinalImpact() - fFirstBunch->GetPos()).Mag() < (b.GetPos() - fFirstBunch->GetPos()).Mag()) {
131        bunch.SetFate(3);
132        return NONE;
133    }
134   
135    // if bunch showerpos is within clouds, switch off the present clear sky propagation
136    if(rtn == CLOUDY && (b.GetNextImpact() - b.GetPos()).Mag() < TOLERANCE) {
137        return rtn;
138    }
139
140    // run propagation
141    Propagate(bunch,list);
142   
143    #ifdef DEBUG
144    bench.Stop(prt);
145    MsgForm(EsafMsg::Debug,"REAL=%6.2f s   CPU=%6.2f s",bench.GetRealTime(prt),bench.GetCpuTime(prt)); 
146    #endif
147
148    if(rtn == GROUND) bunch.SetFate(4);
149    // in the case the impact is out of FoV, NONE is returned to stop the propagation (no ground reflexion)
150    if(rtn == CLEARSKY) {
151        rtn = NONE;
152        bunch.SetFate(5);
153    }
154
155    return rtn;
156}
157 
158//_____________________________________________________________________________
159void AlongTrack_CSPropagator::Propagate(BunchOfPhotons& bunch,ListPhotonsInAtmosphere& list) const {
160    //
161    // Propagate along a track a BunchOfPhotons in clear sky conditions,
162    // generating SinglePhotons all along the travel -> ListPhotonsInAtmosphere filled
163    // NB: bunch can be modified, its copy b cannot
164    //
165   
166#ifdef DEBUG
167    Msg(EsafMsg::Debug) <<"Propagate()" << MsgDispatch;
168#endif
169
170    // initializations
171    const BunchOfPhotons& b = bunch;
172    const EarthVector& impact = b.GetNextImpact();
173    Double_t conv = 1.;
174    Int_t nb = b.GetWlSpectrum().GetAxis().GetNbins();
175    vector<Double_t*> coeff;
176    Double_t* TotalTrans = new Double_t[nb];
177    Double_t* RayleighTrans = new Double_t[nb];
178    Double_t* OzoneTrans = new Double_t[nb];
179    Double_t* AerosolTrans = new Double_t[nb];
180    if(!TotalTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for TotalTrans" << MsgDispatch;
181    if(!RayleighTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for RayleighTrans" << MsgDispatch;
182    if(!OzoneTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for OzoneTrans" << MsgDispatch;
183    if(!AerosolTrans) Msg(EsafMsg::Panic) << "Pb of memory allocation for AerosolTrans" << MsgDispatch;
184    coeff.push_back(TotalTrans);
185    coeff.push_back(RayleighTrans);
186    coeff.push_back(OzoneTrans);
187    coeff.push_back(AerosolTrans);
188    Double_t raylScat(0),aerScat(0),phfunc_rayl(0),phfunc_aer(0);
189    Int_t NbRaylBackscat = 0;
190    Int_t Nb_Aer_Backscat = 0;
191    EarthVector pos, step, entry, exit;
192    if(list.GetNbTrackSteps() == 0) Msg(EsafMsg::Panic) << "this propagator needs ListPhotonsInAtmosphere::fTrack" << MsgDispatch;
193   
194   
195    // find the matching track step
196    UInt_t index = list.TrackStepIndex(b.GetPos());
197    entry = b.GetPos();
198    if(index < (list.GetNbTrackSteps()-1)) exit = list.GetTrackStep(++index);
199    else {
200         Msg(EsafMsg::Warning) << "pb with array index : a bunch position is OUTSIDE light track"<<MsgDispatch;
201         exit = impact;
202#ifdef DEBUG
203Msg(EsafMsg::Debug) << "POSITION OF BUNCH WHEN PROBLEM = " <<b.GetPos() <<MsgDispatch;
204Msg(EsafMsg::Debug) << "POSITION OF BUNCH WHEN PROBLEM = " <<(b.GetPos() - list.GetTrackStep(0)).Mag()*1e-6 <<MsgDispatch;
205Msg(EsafMsg::Debug) << "ALTITUDE OF BUNCH WHEN PROBLEM = " <<b.GetPos().Zv() <<MsgDispatch;
206#endif
207    }
208   
209#ifdef DEBUG
210    // to get nb of iteration for propagation of this bunch
211    Float_t count = 0;
212#endif
213
214    // loop of step by step propagation for backscattering simulation
215    while((impact - b.GetPos()).Dot(impact - b.GetShowerPos()) > TOLERANCE) {
216#ifdef DEBUG
217        count++;
218#endif
219        step = exit - entry;
220        if((impact - (b.GetPos()+step)).Dot(impact - b.GetPos()) < TOLERANCE) step = impact - b.GetPos();
221        fCalcul->Trans(b,b.GetPos()+step,coeff);
222        phfunc_rayl = fCalcul->RayleighPhaseFunction(b.GetDir(),EUSO() - b.GetPos() - 0.5*step);
223        phfunc_aer = fCalcul->Mie_HG_PhaseFunction(b.GetDir(),EUSO() - b.GetPos() - 0.5*step,0.7);
224       
225        // position and time of flight incrementation (bunch position is at the step end)
226        bunch.AddToPosTof(step);
227        bunch.SetPos(exit); // often useless, but SAFER : handle misaligned bunches
228       
229        // cerenkov backscattering (rayleigh - aerosol mie)
230        for(Int_t i=0; i<nb; i++) {
231            raylScat = (1 - coeff[1][i]);
232            aerScat  = (1 - pow(coeff[3][i],fAerosol_correction_factor));
233            NbRaylBackscat  = EsafRandom::Get()->Poisson(b.GetWeight(i) * EusoOmega(b.GetPos() - 0.5*step) * raylScat*phfunc_rayl);
234            Nb_Aer_Backscat = EsafRandom::Get()->Poisson(b.GetWeight(i) * EusoOmega(b.GetPos() - 0.5*step) * aerScat*phfunc_aer);
235            // Singles generation (bunch pos at step end) -- aerosols-ozone absorption also taken into account there
236            if(NbRaylBackscat) GenerateRayleighSingles(i,NbRaylBackscat,raylScat,step,b,list,coeff[3][i],coeff[2][i]);
237            if(Nb_Aer_Backscat) GenerateAerosolSingles(i,Nb_Aer_Backscat,aerScat,step,b,list,coeff[1][i],coeff[2][i]);
238        }
239        // convolutions with transmission coefficients
240        conv = bunch.ConvoluteWl(coeff[0]);
241        bunch.Convolute(conv);
242        entry = exit;
243        if(index != list.GetNbTrackSteps()-1) {
244            index++;
245            exit = list.GetTrackStep(index);
246        }
247        else exit = impact;
248    }
249    // to have position precisely, for later processes
250    bunch.SetPos(impact);
251   
252    #ifdef DEBUG
253    Msg(EsafMsg::Debug) <<"Nb iteration = " <<count << MsgDispatch;
254    Msg(EsafMsg::Debug) <<"after backscat, SinglePhoton list size = "<<list.GetListOfSingle().size() << MsgDispatch;
255    Msg(EsafMsg::Debug) <<"Nb of Singles CORRECTED = "<<fNbSinglesCorrected << MsgDispatch;
256    #endif
257
258    // cleaning of coeff
259    for(UInt_t m=0; m<coeff.size(); m++) {
260        if(coeff[m]) delete[] coeff[m];
261        coeff[m] = 0;
262    }
263    coeff.clear();
264}
265
266//_____________________________________________________________________________
267void AlongTrack_CSPropagator::PreProcess(const BunchOfPhotons& b) const {
268    //
269    // Set fFlag and fImpact, Call RadiativeProcessesCalculator::PreProcess()
270    // Called once (for the first bunch) and used for propagation of all the bunches
271    //
272   
273#ifdef DEBUG
274    Msg(EsafMsg::Debug) <<"PreProcess()" << MsgDispatch;
275#endif
276
277    // get final impact
278    EarthVector finalImpact;
279    fFinal_medium = GetFinalImpact(b,finalImpact);
280    Msg(EsafMsg::Info) <<"PreProcess() : FINAL MEDIUM FOR track IMPACT = "<< fFinal_medium << MsgDispatch;
281    Msg(EsafMsg::Info) <<"PreProcess() : FINAL track IMPACT = "<< finalImpact << MsgDispatch;
282    fFirstBunch = new BunchOfPhotons(b);
283    fFirstBunch->SetFinalImpact(finalImpact);
284    if(finalImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "Bunches should always have an \"impact\"" << MsgDispatch;
285   
286    if(fFinal_medium != NONE) {
287        // create lowtran trans table along the whole track (until final track impact)
288        if(fCalcul->GetName() == "lowtran") {
289            LowtranRadiativeProcessesCalculator* calcul = (LowtranRadiativeProcessesCalculator*) fCalcul;
290            calcul->PreProcess(*fFirstBunch);
291        }
292        fFlag = true;
293    }
294    // if first bunch is under ground, it is not used to initialize the propagator. The next bunch will be
295    else {
296        Msg(EsafMsg::Info) <<"PreProcess() : track starts underground" << MsgDispatch;
297        SafeDelete(fFirstBunch);
298    }
299}
300
301//_____________________________________________________________________________
302Bool_t AlongTrack_CSPropagator::IsAlongTrack(const BunchOfPhotons& b) const {
303    //
304    // Check if bunches are along a track, with first at the track beginning
305    //
306    if(b.GetId() == fFirstBunch->GetId()) return true;
307    Bool_t rtn = false;
308    EarthVector diff = (b.GetPos() - fFirstBunch->GetPos()).Unit();
309   
310    // check if bunch is on the track (defined by firstbunch pos and dir)
311    // and check if firstbunch is really the first one
312    if((diff - fFirstBunch->GetDir()).Mag() < 1.e5*TOLERANCE) rtn = true;
313    else {
314#ifdef DEBUG
315        Msg(EsafMsg::Debug) << "Bunch nb "<<b.GetId()<<" not perfectly along track -> "<<(diff - fFirstBunch->GetDir()).Mag() <<MsgDispatch;
316#endif
317    }
318   
319    return rtn;
320}
321
322//_____________________________________________________________________________
323void AlongTrack_CSPropagator::Reset() {
324    //
325    // Reset method
326    //
327    fCalcul->Reset();
328    fFlag = false;
329    fFinal_medium = CLEARSKY;
330    SafeDelete(fFirstBunch);
331    fNbSinglesCorrected = 0;
332}
333
334
Note: See TracBrowser for help on using the repository browser.