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

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

ESAF version compilable on mac OS
File size: 32.4 KB
Line 
1// $Id: MCRadiativeTransfer.cc 2767 2006-11-15 15:50:14Z moreggia $
2// Author: Sylvain Moreggia   2005/08/16
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: MCRadiativeTransfer                                                  *
8 *  Package: radiativetransfer                                               *
9 *  Coordinator: Sylvain Moreggia                                            *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// MCRadiativeTransfer
16//
17// Detection from a ground detector is now possible in this mode
18//
19//
20//   Config file parameters
21//   ======================
22//
23//   <parameter name>: <parameter description>
24//
25//   - fMScattAnalysis >0 means Detector geometry is completely switched off
26//         - no direct photons
27//         - scattered SinglePhotons are located at their last scattering position
28//         - it is a pure Multiple Scattering Analysis : no propagation to detector (actually there is, but transmissions=1)
29//         - MaxPhi=1
30//         - MaxOmega=1e-20 for fPropagator and =1 for PhotonsInOmega()
31//         - if == 1 : simu checks if photon goes out freely of atmosphere. If it does, it is counted, if not it is lost
32//         - if == 2 : does not check -> allow to study scattering without being biased by last transfer
33//
34
35#include "MCRadiativeTransfer.hh"
36#include "Atmosphere.hh" 
37#include "RadiativeFactory.hh"
38#include "ListPhotonsInAtmosphere.hh"
39#include "ListPhotonsOnPupil.hh"
40#include "EsafRandom.hh"
41#include <TROOT.h>
42#include <TH1D.h>
43#include <TH2D.h>
44#include <TProfile.h>
45#include <math.h> 
46#include "ParentPhoton.hh"
47#include "SinglePhoton.hh"
48#include "Config.hh"
49#include "LowtranRadiativeProcessesCalculator.hh"
50#include "FastRadiativeProcessesCalculator.hh"
51#include "EConst.hh"
52#include "BunchOfPhotons.hh"
53#include "Ground.hh"
54
55#include "EEvent.hh"
56#include "EAtmosphere.hh"
57#include "EAtmosphereBunchAdder.hh"
58#include "EAtmosphereSingleAdder.hh"
59
60#include "TBenchmark.h"
61
62ClassImp(MCRadiativeTransfer)
63
64using namespace TMath;
65using namespace sou;
66using namespace EConst;
67
68//_____________________________________________________________________________
69MCRadiativeTransfer::MCRadiativeTransfer() : RadiativeTransfer() , fPhotons(0), fTotalList(0), fPropagator(0), fTransToDetec(0),
70                                             fNbBunch(0), fNbTot(0), fNbSingles(0), fOmegaMax(0), fDistance_cut(0) {
71    //
72    // ctor
73    //
74    Msg(EsafMsg::Info) << "Enabled" << MsgDispatch;
75    fGround = RadiativeFactory::Get()->GetGround();
76    if(!fGround) Msg(EsafMsg::Panic) << "Pb of memory allocation for fGround" << MsgDispatch;
77   
78    // init some DM
79    Double_t maxtof = Conf()->GetNum("MCRadiativeTransfer.fTofCut")*microsecond;
80    fMScattAnalysis = Int_t(Conf()->GetNum("MCRadiativeTransfer.fMScattAnalysis"));
81    fKeepAll = Conf()->GetBool("MCRadiativeTransfer.fKeepAll");
82    Double_t TOA = Conf()->GetNum("MCRadiativeTransfer.fTOA")*km;
83    fPropagator = new SinglePhotonPropagator(fGround,maxtof,TOA);
84   
85    // config for last trans to detector
86    fTransMode = Conf()->GetStr("MCRadiativeTransfer.fTransMode");
87    if(fTransMode == "fast")
88        fTransToDetec = new FastRadiativeProcessesCalculator();
89    else if(fTransMode == "lowtran")
90        fTransToDetec = new LowtranRadiativeProcessesCalculator();
91    else Msg(EsafMsg::Panic) << "Wrong option for MCRadiativeTransfer.fTransMode" << MsgDispatch;
92   
93   
94    fMaxPhNb = Conf()->GetNum("MCRadiativeTransfer.fMaxPhNb");
95   
96    fScatOrder = Int_t(Conf()->GetNum("MCRadiativeTransfer.fScatOrder"));
97    Msg(EsafMsg::Info) << "Max Scattering order to be simulated = " << fScatOrder << MsgDispatch;
98   
99    // detector geometry : decoupled or optimized
100    // decoupled : - detector seems to be a sphere. Allows to study several detector geometries with the same RadiativeTransfer simulation
101    //             - MUST be used in ground detector mode
102    // optimized : - true detector geometry is used in RadiativeTransfer simulation
103    //             - not available for ground detector
104    string name = Conf()->GetStr("MCRadiativeTransfer.fDecoupled");
105    if(name == "decoupled") fDecoupled = true;
106    else if(name == "optimized") fDecoupled = false;
107    else Msg(EsafMsg::Panic) << "Wrong option for MCRadiativeTransfer.fDecoupled" << MsgDispatch;
108   
109}
110
111//_____________________________________________________________________________
112MCRadiativeTransfer::~MCRadiativeTransfer() {
113    //
114    // Destructor
115    //
116    SafeDelete(fPhotons);
117    SafeDelete(fPropagator);
118    SafeDelete(fTransToDetec);
119}
120
121//______________________________________________________________________________
122Bool_t MCRadiativeTransfer::ClearMemory() {
123    //
124    // Release all the mem hold in the buffers
125    //
126
127    Reset();
128    vector<SinglePhoton*> emptySingle;
129    emptySingle.swap(fTempList);
130   
131    return fPhotons ? fPhotons->ClearMemory() : kTRUE;
132}
133
134
135//_______________________________________________________________________________________________________________________
136PhotonsOnPupil* MCRadiativeTransfer::Get(PhotonsInAtmosphere* photons, const DetectorGeometry* dg) {
137    //
138    // Transport photons from source to Euso pupil
139    //
140
141
142
143    ////////////////////////////////////////////////////////////////////////////////////////////
144    ///////////////////////////////// Initializations /////////////////////////////////////////
145    ////////////////////////////////////////////////////////////////////////////////////////////
146
147
148// 1. Preprocesses
149    // copy DetectorGeometry
150    CopyDetectorGeometry(dg);
151   
152    // in case ground detector simulation
153    if(EUSO().Zv() == 0) fDetAtGrnd = true;
154    else fDetAtGrnd = false;
155    if(fDetAtGrnd && !fDecoupled) Msg(EsafMsg::Panic) << "In Ground detector mode, decoupled mode must be switched on" << MsgDispatch;
156    if(fDetAtGrnd) {
157        Msg(EsafMsg::Warning) << "Ground detector mode enabled" << MsgDispatch;
158        if(fTransMode == "lowtran") Msg(EsafMsg::Panic) << "Ground detector mode requires option  :  MCRadiativeTransfer.fTransMode = fast" << MsgDispatch;
159    }
160   
161    // if given list of photons is empty 
162    if(!photons) return NULL;
163
164    if(photons->GetType() == "empty") Msg(EsafMsg::Panic) <<"When NoLightSource used, NoRadiativeTransfer should be used"<< MsgDispatch;
165    if(photons->GetType() != "list") Msg(EsafMsg::Panic) <<"Wrong PhotonsInAtmosphere format. ListPhotonsInAtmosphere expected."<< MsgDispatch;
166    fTotalList = (ListPhotonsInAtmosphere*) photons;
167
168
169    // if Lowtran FORTRAN code used : build a datacard of VERTICAL transmission for PropagateAllToDetector() and PropagateTempListToDetector(..) method
170    // NB : internally checks if datacard already exists. If so, does not rebuild it
171    if(!fDecoupled && !fMScattAnalysis && fTransMode == "lowtran") {
172        Msg(EsafMsg::Info) <<"Building LOWTRAN vertical transmission datacard"<< MsgDispatch;
173        LowtranRadiativeProcessesCalculator* LowtranToDetec = (LowtranRadiativeProcessesCalculator*)fTransToDetec;
174        LowtranToDetec->MakeVerticalDatacard(fGround,EUSO().Zv());
175        Msg(EsafMsg::Info) <<"LOWTRAN datacard built"<< MsgDispatch;
176    }
177    else if((fDecoupled || fMScattAnalysis) && fTransMode == "lowtran") {
178        Msg(EsafMsg::Warning) <<"Interface with LOWTRAN fortran code not well ADAPTED to decoupled mode (including ground mode) and MScattAnalysis mode"<< MsgDispatch;
179        Msg(EsafMsg::Warning) <<"Should better use configuration  :  MCRadiativeTransfer.fTransMode = fast"<< MsgDispatch;
180    }
181       
182
183
184
185
186// 2. Get some infos about created light
187    // look at number of bunches and photons
188    fNbBunch = fTotalList->GetListOfBunch().size();
189    fNbTot = 0;
190    fNbSingles = 0;
191    Double_t nf(0), nc(0);
192    for(size_t i=0; i<fNbBunch; i++) {
193        fNbTot += (fTotalList->GetListOfBunch()[i])->GetWeight();
194        if(fTotalList->GetListOfBunch()[i]->GetType() == Fluo)
195             nf += fTotalList->GetListOfBunch()[i]->GetWeight();
196        else nc += fTotalList->GetListOfBunch()[i]->GetWeight();
197    }
198
199// some prints
200#ifdef DEBUG
201    Msg(EsafMsg::Debug) << "SIZE OF LIST OF BUNCHES = "<<fNbBunch << MsgDispatch;
202    Msg(EsafMsg::Debug) << "Nb fluo = " <<nf << MsgDispatch;
203    Msg(EsafMsg::Debug) << "Nb ckov = " <<nc << MsgDispatch;
204
205    if ( fNbBunch >=1 ) {
206        EarthVector track = fTotalList->GetListOfBunch()[fNbBunch - 1]->GetPos() 
207                          - fTotalList->GetListOfBunch()[0]->GetPos();
208        Msg(EsafMsg::Debug) << "size of the photon track = "<<track.Mag()/km <<" km"<< MsgDispatch;
209    }
210#endif
211    Msg(EsafMsg::Info) << "TOTAL number of photons ="<<fNbTot << MsgDispatch;
212   
213   
214   
215   
216   
217// 3. Prepare Monte Carlo algorithm
218    // set omega max
219    if(fMScattAnalysis == 0) {
220        // fOmegaMax setting for space detector : last scattering layer is taken at 30km altitude
221        //TOFIX : if a scattering exists higher in altitude
222        fOmegaMax = EusoOmega(EarthVector(0,0,30*km));
223        if(fDetAtGrnd) {
224            fDistance_cut = Conf()->GetNum("MCRadiativeTransfer.fDistance_cut")*km;
225            fOmegaMax = EusoOmega(EarthVector(0,0,fDistance_cut));
226        }
227        fPropagator->SetOmegaMax(fOmegaMax);
228    }
229    else {
230        fOmegaMax = 1;
231        fPropagator->SetOmegaMax(1e-20);  // means that fPropagator will keep every photon that reaches the relevant scat order
232    }
233   
234   
235    // set phase function maximum value for the set of phase functions involved in the simulation :
236    // molecular, ground, clouds, (aerosols to come)
237    if(fMScattAnalysis == 0) {
238        Double_t molecular_value = fPropagator->GetRadiativeCalculator()->GetMaxPhaseFunction("rayleigh");
239        Double_t clouds_value = Atmosphere::Get()->GetClouds()->GetMaxPhaseFunction("truncated");
240        Double_t aerosol_value = Atmosphere::Get()->GetAerosol()->GetMaxPhaseFunction(GetWlRangeMin(),GetWlRangeMax());
241        Double_t ground_value = fGround->GetMaxOutgoing_phase_function();
242        fMaxPhaseFunction = Max(molecular_value,aerosol_value);
243        fMaxPhaseFunction = Max(fMaxPhaseFunction,clouds_value);
244        fMaxPhaseFunction = Max(fMaxPhaseFunction,ground_value);
245        Msg(EsafMsg::Info) << "MAXIMUM OF PHASE FUNCTION = "<<fMaxPhaseFunction << MsgDispatch;
246        Msg(EsafMsg::Info) << "molecular_value = "<<molecular_value << MsgDispatch;
247        Msg(EsafMsg::Info) << "clouds_value = "<<clouds_value << MsgDispatch;
248        Msg(EsafMsg::Info) << "aerosol_value = "<<aerosol_value << MsgDispatch;
249        Msg(EsafMsg::Info) << "ground_value = "<<ground_value << MsgDispatch;
250    }
251    else fMaxPhaseFunction = 1;
252    fPropagator->SetMaxPhaseFunction(fMaxPhaseFunction);
253
254
255
256
257// 4. root output init
258    EEvent* ev = EEvent::GetCurrent();
259    if(ev->GetAtmosphere()) ev->GetAtmosphere()->SetMaxScatOrder(fScatOrder);
260
261    // if some SinglePhotons created directly by LightSource module, they are saved into root
262    if ( ev ) {
263        const vector<SinglePhoton*>& list_single_2 = fTotalList->GetListOfSingle();
264        size_t nbnotsaved = fTotalList->GetNbNotSaved();
265        for (size_t i = list_single_2.size() - nbnotsaved; i<list_single_2.size(); i++ ) {
266            EAtmosphereSingleAdder sa( list_single_2[i],true,false,0 );
267            ev->Fill(sa);
268            fTotalList->OneSingleSaved();
269        }
270    }
271
272
273
274
275    ////////////////////////////////////////////////////////////////////////////////////////////
276    ///////////////////////////////// ALGO STARTS HERE /////////////////////////////////////////
277    ////////////////////////////////////////////////////////////////////////////////////////////
278    BunchOfPhotons* bunch = 0;
279    Bool_t next(0),subnext(0);
280    Double_t nbTracked(0), fraction(0), left(0),right(10),subleft(0),subright(2.5);
281    Int_t split_level = Int_t(ceil(fNbTot * fOmegaMax * fMaxPhaseFunction / fMaxPhNb)) + 1;  // +1 in case of round-off
282    Msg(EsafMsg::Warning) << "NB of photons to propagate after reduction "<<fNbTot * fOmegaMax * fMaxPhaseFunction<< MsgDispatch;
283    Msg(EsafMsg::Warning) << "Each scattering order simu will be split in "<<split_level<<" parts" << MsgDispatch;
284   
285    // 1. Direct photons simulation
286    while(true) {
287        bunch = fTotalList->GetBunch();
288        if(!bunch) break;
289
290        // saves bunch parameters at creation
291        if ( ev ) {
292            EAtmosphereBunchAdder ba( bunch, true );
293            ev->Fill(ba);
294        }
295
296        // photons going directly toward detector
297        if(!fMScattAnalysis) DirectToEuso(*bunch);  // no direct light in MScattAnalysis mode
298    }
299
300    // 2.3 Scan fTempList before filling fTotalList, applying relevant cuts :
301    //     - if photon status > 4, it is lost
302    //     - if GROUND detector, check if photon is a real candidate. If not, it is lost
303    //     - if fKeepAll == false, photons are propagated to detector and those absorbed are lost
304    Msg(EsafMsg::Info)  <<"Nb of direct photons before cuts = " <<fTempList.size() << MsgDispatch;
305    ScanTempList();
306    fTotalList->Add(fTempList);
307   
308
309
310    // 2. SCATTERING SIMULATION loop
311    for(Int_t ord=1; ord <= fScatOrder; ord++) {
312        Msg(EsafMsg::Info) << "\nSCATTERING ORDER = " << ord<< MsgDispatch;
313#ifdef DEBUG
314        TString jobMCRT = "Time spent for Scattering simu:";
315        TBenchmark gB;
316        gB.Start(jobMCRT);
317#endif
318        nbTracked = 0;
319        fraction = 0;
320        left = 0;
321        right = 10;
322        subleft = 0;
323        subright = 2.5;
324        fTotalList->ResetBunchCounter();  // for below extraction from list
325        Int_t loop_j = 0;
326        Int_t remainingPhotons = -1;
327        for(Int_t level=0; level < split_level; level++) {
328            Msg(EsafMsg::Info) << "************ Starting bunches transformation Level nb "<<level+1<<" ************"<< MsgDispatch;
329            if(level == 0) Msg(EsafMsg::Info) << " : 0%" << MsgFlush;
330            for(UInt_t i=loop_j; i <fTotalList->GetBunchEntries(); i++) {
331                bunch = fTotalList->GetBunch(loop_j);
332                loop_j++;
333                nbTracked++;
334
335                // 2.1 generate reduce nb of photons
336                remainingPhotons = PhotonsInOmega(*bunch,remainingPhotons);
337               
338                // if huge nb of photons (>1e6), present scattering order simu is splitted
339                if(remainingPhotons > -1) {
340                    Msg(EsafMsg::Warning) << "\nSplit level nb " <<level+1<<" reached" << MsgDispatch;
341                    loop_j--;
342                    nbTracked--;
343                    break;
344                }
345
346                // dump info
347                fraction = 100*nbTracked/fNbBunch;
348                if (left<fraction && fraction <=right) 
349                    next = kFALSE;
350                else if (fraction>right) {
351                    next = kTRUE;
352                    left = right;
353                    right += 10;
354                }
355                if (subleft<fraction  && fraction <=subright) {
356                    subnext = kFALSE;
357                }
358                if (fraction > subright) {
359                    subnext = kTRUE;
360                    subleft = subright;
361                    subright +=2;
362                }
363                if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
364                if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;       
365            }
366            Msg(EsafMsg::Info) << " -done" << MsgDispatch;
367            Msg(EsafMsg::Info) << "NB of SinglePhoton to PROPAGATE = "<<fNbSingles<< MsgDispatch;
368            fNbSingles = 0;  // reset
369
370
371            // 2.2 scattering simu for this bunch at the relevant scattering order
372            fPropagator->Go(fTempList,ord);
373
374            // 2.2' in MScattAnalysis MODE 1
375            // check if photons would scatter at next order
376            // if not they are KEPT, otherwise they are flagged as "LostByCut"
377            if(fMScattAnalysis == 1) fPropagator->CheckIfScatAgain(fTempList);
378
379
380            // 2.3 Scan fTempList before filling fTotalList, applying relevant cuts :
381            //     - if photon status > 4, it is lost
382            //     - if GROUND detector, check if photon is a real candidate. If not, it is lost
383            //     - if fKeepAll == false, photons are propagated to detector and those absorbed are lost
384            ScanTempList();
385
386            // 2.4 Move photons from fTempList to fTotalList : fTempList does not contain photons anymore
387            //     Only photons which have successfully scattered are written in fTotalList, others are lost
388            fTotalList->Add(fTempList);
389           
390#ifdef DEBUG
391            gB.Stop(jobMCRT);
392            MsgForm(EsafMsg::Info,"Time spent for this scattering order :  REAL=%6.2f s   CPU=%6.2f s",gB.GetRealTime(jobMCRT),gB.GetCpuTime(jobMCRT));
393#endif
394        }
395    }
396
397
398
399
400    // 3. saves single photon parameters after scattering simulation
401    if ( ev ) {
402        const vector<SinglePhoton*>& list_single = fTotalList->GetListOfSingle();
403        size_t nbnotsaved = fTotalList->GetNbNotSaved();
404        for (size_t i = list_single.size() - nbnotsaved; i<list_single.size(); i++ ) {
405            EAtmosphereSingleAdder sa( list_single[i],true,true );
406            ev->Fill(sa);
407            fTotalList->OneSingleSaved();
408        }
409    }
410
411#ifdef DEBUG
412    Msg(EsafMsg::Debug) <<"final list of single, size = " <<fTotalList->GetListOfSingle().size() <<MsgDispatch;
413#endif
414    if(fTotalList->GetNbNotSaved()) Msg(EsafMsg::Warning) <<fTotalList->GetNbNotSaved()<<" SinglePhoton not saved into root" <<MsgDispatch;
415
416
417
418
419    // 4. propagates all direct and scattered photons until detector (Lowtran is used, c.f. Ozone)
420    //    ALLOW to study transmission values along the last travel to detector
421    //    IS fake if fMScattAnalysis > 0
422    if(fKeepAll) PropagateAllToDetector();
423       
424   
425    // 5. Convert PinA in PonP
426    ConvertIntoPonP();
427
428    // saves single photons after propagation until detector
429    if ( ev ) {
430        const vector<SinglePhoton*>& list_single_3 = fTotalList->GetListOfSingle();
431        for ( size_t i=0; i<list_single_3.size(); i++ ) {
432            EAtmosphereSingleAdder sa( list_single_3[i],false,false,i );
433            ev->Fill(sa);
434        }
435    }
436   
437
438    return fPhotons;
439}
440
441//_______________________________________________________________________________________________________________________
442void MCRadiativeTransfer::Reset() {
443    //
444    // get ready for next event
445    // NB : when Lowtran used --> some fTransToDetec internal datacard not reset (used for all the events of a run)
446    //
447    if(fGround) fGround->Reset();
448    if(fPhotons) fPhotons->Clear();
449    if(fPropagator) fPropagator->Reset();
450    for(size_t i=0; i<fTempList.size(); i++) {
451        SafeDelete(fTempList[i]);
452        Msg(EsafMsg::Panic) <<"fTempList SHOULD BE EMPTY, remains : "<<fTempList.size()<<" photons inside"<<MsgDispatch;
453       
454    }
455    if(fTransToDetec) fTransToDetec->Reset();
456    fTempList.clear();
457    fNbBunch = 0;
458    fNbTot = 0;
459    fNbSingles = 0;
460    fOmegaMax = 0;
461}
462
463//_______________________________________________________________________________________________________________________
464void MCRadiativeTransfer::DirectToEuso(const BunchOfPhotons& b) {
465    //
466    // Creates a list of SinglePhoton considering the EUSO solid angle
467    // handles fluo and cerenkov (angular distrib used for the later)
468    // THESE SINGLEPHOTON WILL NOT UNDERGO SCATTERING SIMU
469    //
470   
471    Int_t nb = 0;
472    EarthVector towardEUSO = (EUSO() - b.GetPos()).Unit();
473    Double_t theta = fabs(towardEUSO.Angle(b.GetDir())); // to keep theta within 0-Pi()
474    if(theta > Pi()) Msg(EsafMsg::Warning) << "<DirectToEuso> Pb with angle definition" << MsgDispatch;
475   
476    Double_t angular_distrib_value = b.AngularDist_OverTwoPi(theta);
477    nb = EsafRandom::Get()->Poisson(b.GetWeight() * EusoOmega(b.GetPos()) * angular_distrib_value);
478
479    // SinglePhotons created and added to fTotalList (THEY WILL NOT UNDERGO SCATTERING SIMU)
480    Double_t wl, date, tof;
481    EarthVector showerpos, dir, diff;
482    SinglePhoton* s = 0;
483   
484    UInt_t bid = b.GetId();
485    PhotonType type = b.GetType();
486    tof = 0;
487   
488    for(Int_t i=0; i<nb; i++) {
489        // corrections for date and showerpos from BunchOfPhotons mean values and longitudinal dispersion
490        showerpos = b.RandomPosInShower();
491        if(fGround->IsUnderGround(showerpos)) continue;
492        diff = showerpos - b.GetShowerPos();
493        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
494        dir = EUSO() - showerpos;
495        wl = b.GetWlSpectrum().GetLambda();
496        s = new SinglePhoton(type,date,tof,wl,showerpos,showerpos,dir,Direct,bid,b.GetShowerAge());
497        fTempList.push_back(s);
498    }
499}
500
501//_______________________________________________________________________________________________________________________
502Int_t MCRadiativeTransfer::PhotonsInOmega(const BunchOfPhotons& b, Int_t remainingPhotons) {
503    //
504    // Detector solid angle is applied here -> the resulting photons will undergo scattering simu
505    // Photon features are set here (lambda, position, direction)
506    // SinglePhoton are stored into fTempList for scattering simu
507    //
508    // if remainingPhotons > -1, means the nb of photons to create is specified
509    //
510
511    Int_t nb = 0;
512    Int_t rtn = -1;
513    Double_t wl, date, tof;
514    EarthVector showerpos, dir, diff;
515    SinglePhoton* s = 0;
516    UInt_t bid = b.GetId();
517    PhotonType type = b.GetType();
518    tof = 0;
519   
520    if(fMScattAnalysis) nb = Int_t(b.GetWeight());
521    else if(remainingPhotons > -1) nb = remainingPhotons;
522    else nb = EsafRandom::Get()->Poisson(b.GetWeight() * fOmegaMax * fMaxPhaseFunction);
523   
524    for(Int_t i=0; i<nb; i++) {
525        // get features at creation
526        showerpos = b.RandomPosInShower();
527        if(fGround->IsUnderGround(showerpos)) continue;
528        diff = showerpos - b.GetShowerPos();
529        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
530        dir  = b.RandomDirection();
531        wl   = b.GetWlSpectrum().GetLambda();
532        s    = new SinglePhoton(type,date,tof,wl,showerpos,showerpos,dir,None,bid,b.GetShowerAge());
533        fNbSingles++;
534        fTempList.push_back(s);
535        if(fNbSingles >= fMaxPhNb) {
536            rtn = nb - 1 - i;
537            break;
538        }
539    }
540   
541#ifdef DEBUG
542    //Msg(EsafMsg::Debug)  <<"  InOmega gives -> " <<nb <<" photons" << MsgDispatch;
543#endif
544
545    return rtn;
546}
547
548//_______________________________________________________________________________________________________________________
549void MCRadiativeTransfer::PropagateAllToDetector() {
550    //
551    // Final phase of transfer. Propagate all the SinglePhoton til EUSO pupil
552    // USE LOWTRAN (because of Ozone transmission)
553    // DISABLED if fKeepAll == false, PropagateTempListToDetector() method used instead
554    // FAKE if fMScattAnalysis > 0
555    //
556
557    EarthVector dirtest(1);
558    Double_t Trans[4];
559    Double_t TotTrans(0.);
560    TRandom* rndm = EsafRandom::Get();
561    SinglePhoton* p = 0;
562    size_t nTotal(0),nTracked(0);
563    Float_t fraction(0); 
564    Bool_t next(0),subnext(0);
565    Float_t left(0),right(10),subleft(0),subright(2.5);
566    nTotal = fTotalList->GetSingleEntries();
567    Msg(EsafMsg::Info) << "Propagation All photons to Detector: 0%" << MsgFlush;
568
569
570    // loop over the list of SinglePhoton
571    while(true) {
572        p = fTotalList->GetSingle();
573        if(!p) break;
574        nTracked++;
575               
576        // photon direction is toward EUSO
577        p->SetDir((EUSO() - p->Pos()).Unit());
578
579        // calculate transmission from photon position to EUSO
580        if(!fMScattAnalysis) TotTrans = fTransToDetec->Trans(*p,EUSO(),Trans);
581        else {// MScatt analysis : does not propagate photons, detector does not exist
582            TotTrans = 1;
583            Trans[0] = 1; Trans[1] = 1; Trans[2] = 1; Trans[3] = 1; 
584        }
585        p->SetLastTrans(TotTrans,"tot");
586        p->SetLastTrans(Trans[1],"rayl");
587        p->SetLastTrans(Trans[2],"ozone");
588        p->SetLastTrans(Trans[3],"aero");
589        p->SetLastTrans(TotTrans/(Trans[1]*Trans[2]*Trans[3]),"cloud");
590
591        // tell if photon is absorbed
592        if(TotTrans < rndm->Rndm())  p->SetAbsorbed();
593       
594        // Dump CPU commentaries
595        fraction = 100*Float_t(nTracked)/Float_t(nTotal);
596        if (left<fraction && fraction <=right) 
597            next = kFALSE;
598        else if (fraction>right) {
599            next = kTRUE;
600            left = right;
601            right += 10;
602        }
603        if (subleft<fraction  && fraction <=subright) {
604            subnext = kFALSE;
605        }
606        if (fraction > subright) {
607            subnext = kTRUE;
608            subleft = subright;
609            subright +=2;
610        }
611        if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
612        if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;       
613    }
614    Msg(EsafMsg::Info) << " -done" << MsgDispatch;
615}
616
617//_______________________________________________________________________________________________________________________
618void MCRadiativeTransfer::ScanTempList() {
619    //
620    // remove unrelevant photons (status > 4)
621    // if fKeepAll == false : call PropagateTempListToDetector() --> photons absorbed are not kept
622    //
623    // if GROUND detector mode : check geometrical cuts
624    // check if photon is a candidate : if not, it is lost
625    //
626   
627    SinglePhoton* p = 0;
628    for(UInt_t m=0; m < fTempList.size(); m++) {
629        p = fTempList[m];
630        if(p->Status() > 4) {
631            SafeDelete(fTempList[m]);
632        }
633        else if(fDetAtGrnd && !IsSeenByGroundDetector(*p)) {
634            SafeDelete(fTempList[m]);
635        }
636    }
637   
638    if(!fKeepAll) PropagateTempListToDetector();
639}
640
641//_______________________________________________________________________________________________________________________
642void MCRadiativeTransfer::PropagateTempListToDetector() {
643    //
644    // same as PropagateToDetector but :
645    // - act on fTempList
646    // - disable fine transmission study
647    //
648
649    // init
650    EarthVector dirtest(1);
651    Double_t Trans[4];
652    Double_t TotTrans(0.);
653    TRandom* rndm = EsafRandom::Get();
654    SinglePhoton* p = 0;
655    size_t nTotal(0),nTracked(0);
656    Float_t fraction(0); 
657    Bool_t next(0),subnext(0);
658    Float_t left(0),right(10),subleft(0),subright(2.5);
659    nTotal = fTempList.size();
660    Msg(EsafMsg::Info) << "Propagation TempList to Detector: 0%" << MsgFlush;
661   
662
663    // loop over the list of SinglePhoton
664    for(UInt_t m=0; m < fTempList.size(); m++) {
665        p = fTempList[m];
666        nTracked++;
667        if(!p) continue;
668               
669        // photon direction is toward EUSO
670        p->SetDir((EUSO() - p->Pos()).Unit());
671
672        // calculate transmission from photon position to EUSO
673        if(!fMScattAnalysis) TotTrans = fTransToDetec->Trans(*p,EUSO(),Trans);
674        else {// MScatt analysis : does not propagate photons, detector does not exist
675            TotTrans = 1;
676            Trans[0] = 1; Trans[1] = 1; Trans[2] = 1; Trans[3] = 1; 
677        }
678        // fine transmission study is disabled
679        p->SetLastTrans(-1,"tot");
680        p->SetLastTrans(-1,"rayl");
681        p->SetLastTrans(-1,"ozone");
682        p->SetLastTrans(-1,"aero");
683        p->SetLastTrans(-1,"cloud");
684
685        // photon is transmitted or absorbed
686        if(TotTrans < rndm->Rndm())  p->SetAbsorbed();
687       
688        if(p->IsAbsorbed()) SafeDelete(fTempList[m]);
689       
690       
691        // Dump CPU commentaries
692        fraction = 100*Float_t(nTracked)/Float_t(nTotal);
693        if (left<fraction && fraction <=right) 
694            next = kFALSE;
695        else if (fraction>right) {
696            next = kTRUE;
697            left = right;
698            right += 10;
699        }
700        if (subleft<fraction  && fraction <=subright) {
701            subnext = kFALSE;
702        }
703        if (fraction > subright) {
704            subnext = kTRUE;
705            subleft = subright;
706            subright +=2;
707        }
708        if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
709        if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;       
710    }
711    Msg(EsafMsg::Info) << " -done" << MsgDispatch;
712}
713
714//_______________________________________________________________________________________________________________________
715Bool_t MCRadiativeTransfer::IsSeenByGroundDetector(const SinglePhoton& p) const {
716    //
717    // for ground detector only, assumed to be half a sphere
718    // if dist(photon,detec) < fDistance_cut, photon is lost
719    // distribute photons on a sphere : if photons are on the bottom hemisphere, they are lost
720    //
721
722    // photons with pos.Z < 0 won't be seen by ground detector
723    if(p.Pos().Z() < 0) return false;
724
725    Double_t dist_to_detec = (EUSO() - p.Pos()).Mag();
726
727    // check fDistance_cut
728    if(dist_to_detec < fDistance_cut) return false;
729   
730   
731   
732    // distribute photon on a local sphere (radius = detector radius)
733    // local frame is centered on sphere origin
734    TRandom* rndm = EsafRandom::Get();
735    Double_t Rsphere = EusoRadius();
736    EarthVector pos_on_sphere(1);
737    EarthVector pos_in_MESframe(1);
738    EarthVector pos_local(0,0,dist_to_detec);  // photon position in local frame
739    EarthVector dir_local(1);
740   
741      // position of photon on a horizontal disk cutting sphere in two hemisphere (<=> as sphere is seen from photon position)
742      // this disk contains the local frame origin
743    Double_t r(0.), phi(0.); // polar coordinates
744    r = Rsphere * sqrt(rndm->Rndm());
745    phi = TwoPi() * rndm->Rndm();
746    EarthVector pos_on_plane(r*cos(phi),r*sin(phi),0.);
747    dir_local = (pos_on_plane - pos_local).Unit();
748   
749      // then get photon position on the local sphere
750    Double_t mag(0);
751    Double_t b = pos_on_plane*dir_local;
752    Double_t c = pos_on_plane.Mag2() - Rsphere*Rsphere;
753    pair<Int_t,Double_t*>& solu = findRoots(1.,2*b,c);
754    if(solu.first == 0) {
755        Msg(EsafMsg::Warning) << "<IsSeenByGroundDetector> no impact on local sphere : SHOULD NOT OCCUR, photon is lost" << MsgDispatch;
756        return false;
757    }
758    if(solu.first == 1) mag = solu.second[0];
759    else if(solu.first ==2) mag = min(solu.second[0],solu.second[1]);
760    if((mag - kTolerance) > 0) {
761        Msg(EsafMsg::Warning) << "<IsSeenByGroundDetector> mag should always be < 0" << MsgDispatch;
762        Msg(EsafMsg::Warning) << "<IsSeenByGroundDetector> HERE mag = " <<mag <<"   and it SHOUD NOT OCCUR" << MsgDispatch;
763    }
764    pos_on_sphere = pos_on_plane + mag*dir_local;
765
766      // define rotation and apply it to pos_on_sphere in order to get photons position in MES
767    EarthVector Uz(0,0,1);  // main axis of local sphere
768    EarthVector axis = p.Pos().Cross(Uz);
769    pos_in_MESframe = pos_on_sphere;
770    pos_in_MESframe.Rotate(-p.Pos().Theta(),axis);
771   
772   
773    // check if pos_in_MESframe is in the upper hemisphere
774    if(pos_in_MESframe.Z() >= 0) return true;
775    else return false;
776}
777
778//_______________________________________________________________________________________________________________________
779void MCRadiativeTransfer::ConvertIntoPonP() {
780    //
781    // generate PhotonsOnPupil (fPhotons) from ListPhotonsInAtmosphere (fTotalList)
782    // SinglePhoton are destroyed after their conversion into ParentPhoton
783    //
784   
785    // initializations
786    if(!fPhotons) fPhotons = new ListPhotonsOnPupil((vector<ParentPhoton*>*)NULL);
787    if(!fPhotons) Msg(EsafMsg::Panic) << "MCRadiativeTransfer::PropagateAllToDetector, NULL fPhotons : Memory pb" << MsgDispatch;
788   
789    size_t nTotal(0),nTracked(0);
790    TVector3 local_dir, pos;
791    Float_t fraction(0); 
792    Bool_t next(0),subnext(0);
793    Float_t left(0),right(10),subleft(0),subright(2.5);
794    nTotal = fTotalList->GetSingleEntries();
795    SinglePhoton* p = 0;
796
797    // build PhotonsOnPupil's frame
798    BuildPupilFrame(fPhotons);
799    Msg(EsafMsg::Info) << "Convert PhotonsInAtmosphere into PhotonsOnPupil: 0%" << MsgFlush;
800   
801
802    // convert photons
803    fTotalList->ResetSingleCounter();   // for below extraction from list
804    while(true) {
805        p = fTotalList->GetSingle();
806        if(!p) break;
807        nTracked++;
808       
809        // Sample a random photon position on pupil
810        pos =  p->Pos();
811        local_dir = p->Dir();
812        p->SetPosInAtmo(p->Pos());
813        p->AddToPosTof(EUSO() - p->Pos());
814        RamdomPosOnPupil(fPhotons,pos,local_dir);
815
816
817        // check if photon is within the FoV
818        if(!GetDetGeometry()->IsInFoV(local_dir)) p->SetOutFoV();
819        else p->SetOutFoV(false);
820
821
822        // if photon transmitted, becomes a photon on pupil
823        // FoV 'status' not relevant here (too simply treated, ONLY a flag in RT part)
824        // -->  will be considered in detector part
825        if(!p->IsAbsorbed()) fPhotons->Add(*p,pos,local_dir);
826
827
828        // Dump CPU commentaries
829        fraction = 100*Float_t(nTracked)/Float_t(nTotal);
830        if (left<fraction && fraction <=right) 
831            next = kFALSE;
832        else if (fraction>right) {
833            next = kTRUE;
834            left = right;
835            right += 10;
836        }
837        if (subleft<fraction  && fraction <=subright) {
838            subnext = kFALSE;
839        }
840        if (fraction > subright) {
841            subnext = kTRUE;
842            subleft = subright;
843            subright +=2;
844        }
845        if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
846        if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;       
847    }
848    Msg(EsafMsg::Info) << " -done" << MsgDispatch;
849#ifdef DEBUG
850    Msg(EsafMsg::Debug) << "Number of ParentPhoton = " << fPhotons->GetNphotons() << MsgDispatch;
851#endif
852}
853
854
Note: See TracBrowser for help on using the repository browser.