source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/radiativetransfer/src/BunchRadiativeTransfer.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: 24.2 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: BunchRadiativeTransfer.cc 3023 2012-04-04 14:52:02Z mabl $
3// Sylvain Moreggia created Jan, 15 2004
4
5#include "BunchRadiativeTransfer.hh"
6#include "Atmosphere.hh" 
7#include "RadiativeFactory.hh"
8#include "ListPhotonsInAtmosphere.hh"
9#include "ListPhotonsOnPupil.hh"
10#include "EsafRandom.hh"
11#include <TROOT.h>
12#include <TH1D.h>
13#include <TH2D.h>
14#include <TProfile.h>
15#include <math.h> 
16#include "ParentPhoton.hh"
17#include "SinglePhoton.hh"
18#include "Config.hh"
19#include "LowtranRadiativeProcessesCalculator.hh"
20#include "EConst.hh"
21#include "BunchOfPhotons.hh"
22#include "Ground.hh"
23
24#include "EEvent.hh"
25#include "EAtmosphere.hh"
26#include "EAtmosphereBunchAdder.hh"
27#include "EAtmosphereSingleAdder.hh"
28
29#include "TBenchmark.h"
30
31ClassImp(BunchRadiativeTransfer)
32
33using namespace TMath;
34using namespace sou;
35using namespace EConst;
36
37//#define DEBUG
38
39//_______________________________________________________________________________________________________________________
40BunchRadiativeTransfer::BunchRadiativeTransfer() : RadiativeTransfer(), fPhotons(0), fTotalList(0),
41                                                   fCSpropag(0), fICpropag(0), fNbBunch(0), fNbTot(0) {
42    //
43    // ctor
44    //
45    Msg(EsafMsg::Info) << "Enabled" << MsgDispatch;
46    fGround = RadiativeFactory::Get()->GetGround();
47    if(!fGround) Msg(EsafMsg::Panic) << "Pb of memory allocation for fGround" << MsgDispatch;
48   
49    fCSpropag = RadiativeFactory::Get()->GetClearSkyPropagator(fGround);
50    fICpropag = RadiativeFactory::Get()->GetInCloudsPropagator(fGround);
51    if(!fCSpropag) Msg(EsafMsg::Panic) << "Problem of memory allocation for CSpropag" << MsgDispatch;
52    if(!fICpropag) Msg(EsafMsg::Panic) << "Problem of memory allocation for ICpropag" << MsgDispatch;
53    fTransToDetec = RadiativeFactory::Get()->GetRadiativeProcessesCalculator();
54   
55    // detector geometry : decoupled or optimized
56    // decoupled : detector seems to be a sphere. Allows to study several detector geometries with the same RadiativeTransfer simulation
57    // optimized : true detector geometry is used in RadiativeTransfer simulation
58    string name = Conf()->GetStr("BunchRadiativeTransfer.fDecoupled");
59    if(name == "decoupled") fDecoupled = true;
60    else if(name == "optimized") fDecoupled = false;
61    else Msg(EsafMsg::Panic) << "Wrong option for BunchRadiativeTransfer.fDecoupled" << MsgDispatch;
62   
63    // cloud treatment : if "yes" --> cloud scattering simu, if "no" --> just transmission effect
64    fCloudStatus = Conf()->GetBool("BunchRadiativeTransfer.fCloudStatus");
65    if(!fCloudStatus && fCSpropag->GetOrder() > 1)
66         Msg(EsafMsg::Panic) << "<ctor> No cloud scattering simulation possible only with O1_CskyPropagator" << MsgDispatch;
67   
68}
69
70//_______________________________________________________________________________________________________________________
71BunchRadiativeTransfer::~BunchRadiativeTransfer() {
72    //
73    // dtor
74    //   
75    SafeDelete(fPhotons);
76    SafeDelete(fCSpropag);
77    SafeDelete(fICpropag);
78    SafeDelete(fTransToDetec);
79}
80
81//_______________________________________________________________________________________________________________________
82PhotonsOnPupil* BunchRadiativeTransfer::Get(PhotonsInAtmosphere* photons, const DetectorGeometry* dg) {
83    //
84    // Transport photons from source to Euso pupil
85    //
86
87
88    // set DetectorGeometry for all propagators
89    CopyDetectorGeometry(dg);
90   
91   
92    // in case ground detector simulation    //GRNDetec
93    if(EUSO().Zv() == 0) fDetAtGrnd = true;
94    else fDetAtGrnd = false;
95    if(fDetAtGrnd && !fDecoupled) Msg(EsafMsg::Panic) << "In Ground detector mode, decoupled mode must be switched on" << MsgDispatch;
96
97
98    // if given list of photons is empty 
99    if(!photons) return NULL;
100   
101
102#ifdef DEBUG
103    TString jobRT = "Time spent for transfer process:";
104    TBenchmark gB;
105    gB.Start(jobRT);
106#endif
107
108    if(photons->GetType() == "empty") Msg(EsafMsg::Panic) <<"When NoLightSource used, NoRadiativeTransfer must be used"<< MsgDispatch;
109    if(photons->GetType() != "list") Msg(EsafMsg::Panic) <<"Wrong PhotonsInAtmosphere format. ListPhotonsInAtmosphere expected."<< MsgDispatch;
110    fTotalList = (ListPhotonsInAtmosphere*) photons;
111    fCSpropag->SetEndFoV(fTotalList->GetTrackStep(fTotalList->GetNbTrackSteps()-1)); //TOFIX when general FoV handling
112    fICpropag->SetEndFoV(fTotalList->GetTrackStep(fTotalList->GetNbTrackSteps()-1)); //TOFIX when general FoV handling
113
114
115    // if fTransToDetec is LOWTRAN7
116    if(fTransToDetec->GetName() == "lowtran") {
117        LowtranRadiativeProcessesCalculator* TransCalcul = (LowtranRadiativeProcessesCalculator*) fTransToDetec;
118        // Build a datacard of VERTICAL Lowtran transmission for PropagateToDetector() method
119        if(!fDecoupled) {
120            Msg(EsafMsg::Info) <<"Building LOWTRAN vertical transmission datacard"<< MsgDispatch;
121            TransCalcul->MakeVerticalDatacard(fGround,EUSO().Zv());
122            Msg(EsafMsg::Info) <<"LOWTRAN datacard built"<< MsgDispatch;
123        }
124        // Build a datacard of ALONG TRACK Lowtran transmission for PropagateToDetector() method
125        else {
126            Msg(EsafMsg::Info) <<"Building LOWTRAN datacard ALONG TRACK (not vertical)"<< MsgDispatch;
127            TransCalcul->MakeTrackToDetectorDatacard(fGround,EUSO(),*fTotalList);
128            Msg(EsafMsg::Info) <<"LOWTRAN datacard built"<< MsgDispatch;
129        }
130    }
131
132
133    // number of bunches
134    fNbBunch = fTotalList->GetListOfBunch().size();
135    fNbTot = 0;
136    Double_t nf = 0;
137    Double_t nc = 0;
138    for(size_t i=0; i<fNbBunch; i++) {
139        fNbTot += (fTotalList->GetListOfBunch()[i])->GetWeight();
140        if(fTotalList->GetListOfBunch()[i]->GetType() == Fluo) nf += fTotalList->GetListOfBunch()[i]->GetWeight();
141        else nc += fTotalList->GetListOfBunch()[i]->GetWeight();
142    }
143#ifdef DEBUG
144    Msg(EsafMsg::Debug) << "SIZE OF LIST OF BUNCHES = "<<fNbBunch << MsgDispatch;
145    Msg(EsafMsg::Debug) << "Nb fluo = " <<nf << MsgDispatch;
146    Msg(EsafMsg::Debug) << "Nb ckov = " <<nc << MsgDispatch;
147#endif
148    Msg(EsafMsg::Info) << "TOTAL number of photons ="<<fNbTot << MsgDispatch;
149
150    if ( fNbBunch >=1 ) {
151        EarthVector track = fTotalList->GetListOfBunch()[fNbBunch - 1]->GetPos() 
152            - fTotalList->GetListOfBunch()[0]->GetPos();
153#ifdef DEBUG
154        Msg(EsafMsg::Debug) << "size of the photon track = "<<track.Mag()/km <<" km"<< MsgDispatch;
155#endif
156    }
157
158
159    EEvent* ev = EEvent::GetCurrent();
160    if(ev && ev->GetAtmosphere()) ev->GetAtmosphere()->SetMaxScatOrder(1);
161
162
163    // if some SinglePhotons created directly by LightSource module, they are saved into root
164    if ( ev ) {
165        const vector<SinglePhoton*>& list_single_2 = fTotalList->GetListOfSingle();
166        size_t nbnotsaved = fTotalList->GetNbNotSaved();
167        for (size_t i = list_single_2.size() - nbnotsaved; i<list_single_2.size(); i++ ) {
168            EAtmosphereSingleAdder sa( list_single_2[i],true,false,0 );
169            ev->Fill(sa);
170            fTotalList->OneSingleSaved();
171        }
172    }
173
174
175    // loop to transport all the bunches
176    BunchOfPhotons* bunch = 0;
177    Double_t nbTracked(0);
178    Int_t progress = 0;
179
180    while(true) {
181        bunch = fTotalList->GetBunch();
182        if(!bunch) break;
183        nbTracked++;
184
185        // saves bunch parameters at creation
186        if ( ev ) {
187            EAtmosphereBunchAdder ba( bunch, true );
188            ev->Fill(ba);
189        }
190        // propagation of bunches and creation of single photons
191        BunchPropagation(*bunch);
192
193#ifdef DEBUG
194        Msg(EsafMsg::Debug) <<"after the process of this bunch, size of list of single = " <<fTotalList->GetListOfSingle().size() <<MsgDispatch;
195#endif
196
197        // saves single photon parameters at creation
198        if ( ev ) {
199            const vector<SinglePhoton*>& list_single = fTotalList->GetListOfSingle();
200            size_t nbnotsaved = fTotalList->GetNbNotSaved();
201            for (size_t i = list_single.size() - nbnotsaved; i<list_single.size(); i++ ) {
202                EAtmosphereSingleAdder sa( list_single[i],true,true );
203                ev->Fill(sa);
204                fTotalList->OneSingleSaved();
205            }
206        }
207
208        if (100.*nbTracked/fNbBunch >= progress) {
209            Msg(EsafMsg::Info).SetProgress(progress);
210            Msg(EsafMsg::Info) << "Bunches Processing:" << MsgCount;
211            progress+=10;
212        }
213    }
214
215#ifdef DEBUG
216    Msg(EsafMsg::Debug) <<"final list of single, size = " <<fTotalList->GetListOfSingle().size() <<MsgDispatch;
217#endif
218    if(fTotalList->GetNbNotSaved()) Msg(EsafMsg::Warning) <<fTotalList->GetNbNotSaved()<<" SinglePhoton not saved into root" <<MsgDispatch;
219
220    // propagates single photons
221    PropagationOfSingles();
222
223    // saves single photons after propagation in atmosphere
224    if ( ev ) {
225        const vector<SinglePhoton*>& list_single_3 = fTotalList->GetListOfSingle();
226        for ( size_t i=0; i<list_single_3.size(); i++ ) {
227            EAtmosphereSingleAdder sa( list_single_3[i],false,false,i );
228            ev->Fill(sa);
229        }
230    }
231
232
233#ifdef DEBUG
234    gB.Stop(jobRT);
235    MsgForm(EsafMsg::Debug,"Time spent for transfer process:  REAL=%6.2f s   CPU=%6.2f s",gB.GetRealTime(jobRT),gB.GetCpuTime(jobRT));
236#endif
237
238    return fPhotons;
239}
240
241//______________________________________________________________________________
242Bool_t BunchRadiativeTransfer::ClearMemory() {
243    //
244    // Release all the mem hold in the buffers
245    //
246
247    Reset();
248    return fPhotons ? fPhotons->ClearMemory() : kTRUE;
249}
250
251//______________________________________________________________________________
252void BunchRadiativeTransfer::Reset() {
253    //
254    // get ready for next event
255    // NB : if fTransToDetec is LOWTRAN7 --> not really reset, its datacard is used for all the events of a run
256    //
257    if(fGround)     fGround->Reset();
258    if(fPhotons)    fPhotons->Clear();
259    if(fCSpropag)   fCSpropag->Reset();
260    if(fICpropag)   fICpropag->Reset();
261    if(fTransToDetec) fTransToDetec->Reset();
262    fNbBunch = 0;
263    fNbTot = 0;
264}
265
266//_______________________________________________________________________________________________________________________
267void BunchRadiativeTransfer::DirectToEuso(const BunchOfPhotons& bigbunch) const {
268    //
269    // Creates a list of SinglePhoton considering the EUSO solid angle
270    // handles fluo and cerenkov (angular distrib used for the later)
271    //
272   
273    Int_t nb = 0;
274    EarthVector towardEUSO = (EUSO() - bigbunch.GetPos()).Unit();
275    Double_t theta = fabs(towardEUSO.Angle(bigbunch.GetDir())); // to keep theta within 0-Pi()
276    if(theta > Pi()) Msg(EsafMsg::Warning) << "<DirectToEuso> Pb with angle definition" << MsgDispatch;
277    Double_t angular_distrib_value = bigbunch.AngularDist_OverTwoPi(theta);
278    nb = EsafRandom::Get()->Poisson(bigbunch.GetWeight() * EusoOmega(bigbunch.GetPos()) * angular_distrib_value);
279    GenerateDirectSingles(bigbunch,nb);
280
281#ifdef DEBUG
282    Msg(EsafMsg::Debug)  <<"  Direct gives -> " <<nb <<" photons" << MsgDispatch;
283#endif
284}
285
286//_______________________________________________________________________________________________________________________
287void BunchRadiativeTransfer::BunchPropagation( BunchOfPhotons& bigbunch ) {
288    //
289    // Propagation of a bunch through the atmosphere
290    //
291    // temp remark : bunch splitting not used so far, thus don't handle for root filling //TOFIX
292    //
293
294#ifdef DEBUG
295    Msg(EsafMsg::Debug)  << "\nPropagation of the BUNCH nb " << bigbunch.GetId() <<MsgDispatch;
296    Msg(EsafMsg::Debug)  << "Weight = " << bigbunch.GetWeight() <<MsgDispatch;
297    const ParentBunch* parentb = bigbunch.GetParent();
298    Double_t length = (parentb->GetShowerPosi() - parentb->GetShowerPosf()).Mag()/km;
299    Msg(EsafMsg::Debug)  << "Bunch longitudinal extension = "<<length <<" km" << MsgDispatch;
300#endif   
301
302    // if bunch is underground, its simulation stops here
303    if(fGround->IsUnderGround(bigbunch.GetParent()->GetShowerPosf()))   //TOFIX : a part can be above ground (also see related point in ListPinAtmo::fTrack, O2_CSPropag)
304        bigbunch.SetFate(3);
305    else {
306        // photons directly emitted within the Euso solid angle
307        DirectToEuso(bigbunch);
308        Medium state = CLEARSKY;
309
310        // if bunch weight < 1% "mean bunch weight" -> not propagated, because it doesn't contribute
311        if( (bigbunch.GetWeight() < 0.001*fNbTot/fNbBunch)) {
312            bigbunch.SetFate(2);
313            state = NONE;
314        }
315
316        // relevant propagator called
317        while(state > NONE) {
318            switch(state) {
319
320                case CLEARSKY :       state = fCSpropag->Go(bigbunch,*fTotalList);
321                                      break;
322
323                case CLOUDY   :       if(!fCloudStatus) Msg(EsafMsg::Panic) << "<BunchPropagation> SHOULD NOT OCCUR, no cloud treatment" << MsgDispatch;
324                                      state = fICpropag->Go(bigbunch,*fTotalList);
325                                      break;
326
327                case GROUND   :       GroundReflection(bigbunch);
328                                      state = NONE;
329                                      bigbunch.SetFate(4);
330                                      break;
331
332                case AEROSOLS :       Msg(EsafMsg::Panic) << "AerosolsPropagator (for MS) not implemented" << MsgDispatch;
333                                      break;
334
335                default       :       Msg(EsafMsg::Panic) << "Invalid type of medium in bunch propagation = "<<state << MsgDispatch;
336            }
337        }
338    }
339#ifdef DEBUG
340    switch(bigbunch.GetFate()) {
341        case 1    : Msg(EsafMsg::Debug)  << "NOT PROPAGATED -> FLUO" << MsgDispatch; break;
342        case 2    : Msg(EsafMsg::Debug)  << "NOT PROPAGATED -> TOO SMALL WEIGHT" << MsgDispatch; break;
343        case 3    : Msg(EsafMsg::Debug)  << "NOT PROPAGATED -> CREATED UNDERGROUND" << MsgDispatch; break;
344        case 4    : Msg(EsafMsg::Debug)  << "HAS REACHED GROUND" << MsgDispatch; break;
345        case 5    : Msg(EsafMsg::Debug)  << "GONE OUT Euso FoV" << MsgDispatch; break;
346        default   : Msg(EsafMsg::Debug)  << "PB with STATE, fate = "<< bigbunch.GetFate()<< MsgDispatch;
347    }
348#endif
349
350
351    // saves propagated bunch parameters
352    EEvent* ev = EEvent::GetCurrent();
353    if ( ev ) {
354        EAtmosphereBunchAdder ba(&bigbunch,false);
355        ev->Fill(ba);
356    }
357}
358
359//_____________________________________________________________________________________________________
360void BunchRadiativeTransfer::GroundReflection(const BunchOfPhotons& b) const {
361    //
362    // Process of a bunch reaching ground
363    //
364   
365    // initializations
366#ifdef DEBUG
367    Msg(EsafMsg::Debug) << "GroundReflection" << MsgDispatch; 
368#endif
369    if(fabs(b.GetPos().Zv() - fGround->Altitude(b.GetPos())) > 1)
370        Msg(EsafMsg::Warning) << "GroundReflection() must be called if bunch has reached the ground" << MsgDispatch;
371   
372    // determination of number of singlephotons produced
373    Double_t phfunc = fGround->Outgoing_phase_function(b.GetPos(),b.GetDir(),(EUSO() - b.GetPos()).Unit());
374    Double_t weight = b.GetWeight();
375    if(!fCloudStatus) weight = b.GetWeightNocloud();
376    Int_t nb = EsafRandom::Get()->Poisson(weight * fGround->Albedo(b.GetPos()) * EusoOmega(b.GetPos()) * phfunc);
377   
378    // creation of singlephotons
379    GenerateReflectedSingles(nb,b);
380}
381   
382//_______________________________________________________________________________________________________________________
383void BunchRadiativeTransfer::PropagationOfSingles() {
384    //
385    // Final phase of transfer. Propagate all the SinglePhoton til EUSO pupil
386    //
387
388    // initializations
389    if(!fPhotons) fPhotons = new ListPhotonsOnPupil((vector<ParentPhoton*>*)NULL);
390    if(!fPhotons) Msg(EsafMsg::Panic) << "BunchRadiativeTransfer::PropagationOfSingles, NULL fPhotons : Memory pb" << MsgDispatch;
391
392    // build PhotonsOnPupil's frame
393    BuildPupilFrame(fPhotons);
394   
395    EarthVector dirtest(1);
396    Double_t Trans[4];
397    Double_t TotTrans(0.);
398    TRandom* rndm = EsafRandom::Get();
399    SinglePhoton* p = 0;
400    size_t nTotal(0),nTracked(0);
401    Int_t progress = 0;
402    TVector3 local_dir, pos;
403
404    nTotal = fTotalList->GetSingleEntries();
405
406    // loop over the list of SinglePhoton
407    while(true) {
408        p = fTotalList->GetSingle();
409        if(!p) break;
410        nTracked++;
411        dirtest = (EUSO() - p->Pos()).Unit();
412        if((dirtest - p->Dir()).Mag() > TOLERANCE) Msg(EsafMsg::Warning) << "PropagationOfSingles : SinglePhoton must be directed toward EUSO" << MsgDispatch;
413
414        // calculate transmission from photon position to EUSO
415        TotTrans = fTransToDetec->Trans(*p,EUSO(),Trans);
416        p->SetLastTrans(TotTrans,"tot");
417        p->SetLastTrans(Trans[1],"rayl");
418        p->SetLastTrans(Trans[2],"ozone");
419        p->SetLastTrans(Trans[3],"aero");
420        p->SetLastTrans(TotTrans/(Trans[1]*Trans[2]*Trans[3]),"cloud");
421       
422        // check if photon has been absorbed previously during bunch propagation
423        // only for the case of fCloudStatus = false (special mode !!)
424        Double_t RNDM = rndm->Rndm();
425        if(!fCloudStatus && p->IsAbsorbed() && p->IsCloudAbsorbed()) p->SetAbsorbed();
426
427        // photon is transmitted or absorbed
428        else if(TotTrans < RNDM)  {
429            p->SetAbsorbed();
430            // precise if absorbed by clouds, only for direct fluo and reflected cerenkov (single scattering)
431            if(!fCloudStatus && ( Trans[1]*Trans[2]*Trans[3] > RNDM)) {
432                if(p->Status() == 0 || (p->Status() == 1 && p->Type() == 1 && p->NbOfInteractions() == 1) )
433                    p->SetCloudAbsorbed();
434            }
435        }
436       
437        // Sample a random photon position on pupil
438        pos =  p->Pos();
439        local_dir = p->Dir();
440        p->SetPosInAtmo(p->Pos());
441        p->AddToPosTof(EUSO() - p->Pos());
442        RamdomPosOnPupil(fPhotons,pos,local_dir);
443       
444        // check if photon is within the FoV
445        if(!GetDetGeometry()->IsInFoV(local_dir)) p->SetOutFoV();
446        else p->SetOutFoV(false);
447       
448        // if photon transmitted, becomes a photon on pupil
449        // FoV 'status' not relevant here (too simply treated, ONLY a flag in RT part)
450        // -->  will be considered in detector part
451        if(!p->IsAbsorbed() || p->IsCloudAbsorbed()) fPhotons->Add(*p,pos,local_dir);
452
453        //GRNDetec
454        // for ground detector simulation, cos(theta) effect is taken into account here
455        if(fDetAtGrnd) {
456            EarthVector enterdir = -dirtest;
457            if((rndm->Rndm() > cos(enterdir.Theta())) || (enterdir.Theta() > PiOver2())) p->SetAbsorbed();
458        }
459
460        // Dump CPU commentaries
461        if (100.*nTracked/nTotal >= progress) {
462            Msg(EsafMsg::Info).SetProgress(progress);
463            Msg(EsafMsg::Info) << "Singles Processing:" << MsgCount;
464            progress+=10;
465        }
466    }
467
468#ifdef DEBUG
469    Msg(EsafMsg::Debug) << "Number of ParentPhoton = " << fPhotons->GetNphotons() << MsgDispatch;
470#endif
471}
472
473//_______________________________________________________________________________________________________________________
474void BunchRadiativeTransfer::GenerateDirectSingles(const BunchOfPhotons& b,Int_t nb) const {
475    //
476    // SinglePhoton generation (BunchOfPhotons portion which is created within EUSO solid angle)
477    // SinglePhotons created added to fTotalList
478    //
479    Double_t wl, date, tof;
480    EarthVector showerpos, dir, diff;
481    SinglePhoton* s = 0;
482   
483    UInt_t bid = b.GetId();
484    PhotonType type = b.GetType();
485    tof = 0;
486   
487    for(Int_t i=0; i<nb; i++) {
488        // corrections for date and showerpos from BunchOfPhotons mean values and longitudinal dispersion
489        //showerpos = b.RandomPosInShower();
490        showerpos = DistributeToDensity(b.GetShowerPosi(),b.GetShowerPosf());//random position in shower ACCORDING TO DENSITY
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        fTotalList->Add(s);
498    }
499}
500
501//_____________________________________________________________________________________________________
502void BunchRadiativeTransfer::GenerateReflectedSingles(Int_t nb, const BunchOfPhotons& b) const {
503    //
504    // Creates SinglePhoton objects coming from a BunchOfPhotons reflection on ground
505    //   
506#ifdef DEBUG
507    Msg(EsafMsg::Debug)<<"nb of single produced in reflexion process = "<<nb <<MsgDispatch;
508#endif
509    TRandom* rndm = EsafRandom::Get();
510    Double_t wl, date, tof;
511    EarthVector showerpos, pos, dir, diff;
512    SinglePhoton* s = 0;
513   
514    UInt_t bid = b.GetId();
515    PhotonType type = b.GetType();
516    date = b.GetDate();
517   
518    for(Int_t i=0; i<nb; i++) {
519        // corrections for date and showerpos, from BunchOfPhotons mean values and longitudinal+lateral distributions at creation
520        showerpos = b.RandomPosInShower();
521       
522        // if showerpos is underground
523        if(fGround->IsUnderGround(showerpos)) continue;
524       
525        diff = showerpos - b.GetShowerPos();
526        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
527        // tof and position corrections, due to angular distributions at creation
528        pos = b.PosRandomAngCorrec(showerpos,b.GetPos());   // here is a fake correction in position, used to get direction for new impact calculation
529        pos = fGround->GetImpact(showerpos,(pos - showerpos).Unit());
530        // if no impact
531        if(pos.Z() == HUGE) continue;
532        tof = b.GetTof() * (showerpos - pos).Mag()/(b.GetPos() - b.GetShowerPos()).Mag();
533        dir = EUSO() - pos;
534        wl = b.GetWlSpectrum().GetLambda();
535        s = new SinglePhoton(type,date,tof,wl,showerpos,pos,dir,Reflected,bid,b.GetShowerAge());
536        s->AddInteraction();
537        s->AddHistory(Reflected);
538        if(!fCloudStatus) {
539            Double_t transcloud = b.GetWeight() / b.GetWeightNocloud();
540            if(rndm->Rndm() > transcloud) {
541                s->SetAbsorbed();
542                s->SetCloudAbsorbed();
543            }
544        }
545        fTotalList->Add(s);
546    } 
547}
548
549EarthVector BunchRadiativeTransfer::DistributeToDensity(EarthVector initial,EarthVector final) const {
550  //
551  // This function distributes the photons according to the slant depth inside a
552  // a step. This makes the photon production process more physical.
553  // The light curves looks like smoother.
554  //
555
556  // FIXME: This assumes a monotone dependence of the Air Density between
557  //         initial and final position!
558
559    const int numberOfSteps = 10000;
560
561
562
563    const Atmosphere* atm= Atmosphere::Get();
564    const Double_t density1=atm->Air_Density(initial);
565    const Double_t density2=atm->Air_Density(final);
566    const Bool_t invertedDensityDirection = density2 < density1;
567
568    const EarthVector v1 = invertedDensityDirection ? final : initial;
569    const EarthVector v2 = invertedDensityDirection ? initial : final;
570    const EarthVector stepSize = (v2 - v1) * (1.0 / numberOfSteps);
571
572    TRandom* rndm = EsafRandom::Get();
573    const Double_t rndDepth = rndm->Rndm();
574
575
576    const Double_t density_FINAL=density1+rndDepth*(density2-density1);//random sampled density
577
578    EarthVector pos = v1;
579
580    int firstStep = 0;
581    int lastStep = numberOfSteps;
582
583    while (firstStep <= lastStep) {
584
585       const int midStep = (firstStep + lastStep) / 2;
586       pos = v1 + midStep * stepSize;
587       const Double_t den_INTER=atm->Air_Density(pos);
588
589       if (density_FINAL > den_INTER) {
590
591         // repeat search in top half.
592         firstStep = midStep + 1;
593
594       } else if (density_FINAL < den_INTER) {
595         // repeat search in bottom half.
596         lastStep = midStep - 1;
597
598       } else {
599
600         // Found exact value, or return one step too low
601         break;
602       }
603     }
604
605#ifdef DEBUG
606
607    Msg(EsafMsg::Debug) << MsgDispatch;
608
609    Msg(EsafMsg::Debug) << "Distributing bunch position over depth:" << MsgDispatch;
610    Msg(EsafMsg::Debug) << "Initial: (" << initial.X() << ", "
611                                     << initial.Y() << ", "
612                                     << initial.Z() << ")"
613                                     << "\t Density: " << density1 << MsgDispatch;
614
615    Msg(EsafMsg::Debug) << "Final: (" << final.X() << ", "
616                                     << final.Y() << ", "
617                                     << final.Z() << ")"
618                                     << "\t Density: " << density2 << MsgDispatch;
619
620    Msg(EsafMsg::Debug) << "Random density: " << density_FINAL << MsgDispatch;
621    Msg(EsafMsg::Debug) << "Position of rnd density: ("
622                        << pos.X() << ", "
623                        << pos.Y() << ", "
624                        << pos.Z() << ")" << MsgDispatch;
625
626    // The desinity at the location should be equal to random density
627    Msg(EsafMsg::Debug) << "Achived density: " <<  atm->Air_Density(pos) << MsgDispatch;
628
629    // Newline for formating
630
631    Msg(EsafMsg::Debug) << MsgDispatch;
632
633
634#endif
635
636    return pos;
637}
Note: See TracBrowser for help on using the repository browser.