source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/radiativetransfer/src/PureMCRadiativeTransfer.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: 13.8 KB
Line 
1// $Id: PureMCRadiativeTransfer.cc 2767 2006-11-15 15:50:14Z moreggia $
2// Author: Sylvain Moreggia   2006/03/09
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: PureMCRadiativeTransfer                                                           *
8 *  Package: <packagename>                                                   *
9 *  Coordinator: <coordinator>                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// PureMCRadiativeTransfer
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 "PureMCRadiativeTransfer.hh"
27#include "Atmosphere.hh" 
28#include "RadiativeFactory.hh"
29#include "ListPhotonsInAtmosphere.hh"
30#include "ListPhotonsOnPupil.hh"
31#include "EsafRandom.hh"
32#include <TROOT.h>
33#include <TH1D.h>
34#include <TH2D.h>
35#include <TProfile.h>
36#include <math.h> 
37#include "ParentPhoton.hh"
38#include "SinglePhoton.hh"
39#include "Config.hh"
40#include "LowtranRadiativeProcessesCalculator.hh"
41#include "FastRadiativeProcessesCalculator.hh"
42#include "EConst.hh"
43#include "BunchOfPhotons.hh"
44#include "Ground.hh"
45
46#include "EEvent.hh"
47#include "EAtmosphere.hh"
48#include "EAtmosphereBunchAdder.hh"
49#include "EAtmosphereSingleAdder.hh"
50
51#include "TBenchmark.h"
52
53using namespace TMath;
54using namespace sou;
55using namespace EConst;
56
57
58ClassImp(PureMCRadiativeTransfer)
59
60//_____________________________________________________________________________
61PureMCRadiativeTransfer::PureMCRadiativeTransfer() : RadiativeTransfer(), fTotalList(0), fPropagator(0),
62                                             fNbBunch(0), fNbTot(0), fNbSingles(0) {
63    //
64    // Constructor
65    //
66
67    Msg(EsafMsg::Info) << "Enabled" << MsgDispatch;
68    fGround = RadiativeFactory::Get()->GetGround();
69    if(!fGround) Msg(EsafMsg::Panic) << "Pb of memory allocation for fGround" << MsgDispatch;
70   
71    Double_t maxtof = Conf()->GetNum("PureMCRadiativeTransfer.fTofCut")*microsecond;
72    Double_t TOA = Conf()->GetNum("PureMCRadiativeTransfer.fTOA")*km;
73    fPropagator = new SinglePhotonPropagator(fGround,maxtof,TOA);
74   
75    fMaxPhNb = Conf()->GetNum("PureMCRadiativeTransfer.fMaxPhNb");   
76}
77
78//_____________________________________________________________________________
79PureMCRadiativeTransfer::~PureMCRadiativeTransfer() {
80    //
81    // Destructor
82    //
83    SafeDelete(fPropagator);
84}
85
86//______________________________________________________________________________
87Bool_t PureMCRadiativeTransfer::ClearMemory() {
88    //
89    // Release all the mem hold in the buffers
90    //
91
92    Reset();
93    vector<SinglePhoton*> emptySingle;
94    emptySingle.swap(fTempList);
95   
96    return kTRUE;
97}
98
99
100//_______________________________________________________________________________________________________________________
101PhotonsOnPupil* PureMCRadiativeTransfer::Get(PhotonsInAtmosphere* photons, const DetectorGeometry* dg) {
102    //
103    // Transport photons in atmosphere : PhotonsOnPupil generated MUST be considered as fake photons
104    //
105
106
107
108    ////////////////////////////////////////////////////////////////////////////////////////////
109    ///////////////////////////////// Initializations /////////////////////////////////////////
110    ////////////////////////////////////////////////////////////////////////////////////////////
111
112
113// 1. Preprocesses
114    // copy DetectorGeometry
115    CopyDetectorGeometry(dg);
116       
117    // if given list of photons is empty 
118    if(!photons) return NULL;
119
120    if(photons->GetType() == "empty") Msg(EsafMsg::Panic) <<"When NoLightSource used, NoRadiativeTransfer should be used"<< MsgDispatch;
121    if(photons->GetType() != "list") Msg(EsafMsg::Panic) <<"Wrong PhotonsInAtmosphere format. ListPhotonsInAtmosphere expected."<< MsgDispatch;
122    fTotalList = (ListPhotonsInAtmosphere*) photons;
123
124
125
126// 2. Get some infos about created light
127    // look at number of bunches and photons
128    fNbBunch = fTotalList->GetListOfBunch().size();
129    if(fNbBunch > 1) Msg(EsafMsg::Panic) << "PureMCRadiativeTransfer code CANNOT handle more than one bunch of photons"<< MsgDispatch;
130    fNbTot = 0;
131    fNbSingles = 0;
132    Double_t nf(0), nc(0);
133    for(size_t i=0; i<fNbBunch; i++) {
134        fNbTot += (fTotalList->GetListOfBunch()[i])->GetWeight();
135        if(fTotalList->GetListOfBunch()[i]->GetType() == Fluo)
136             nf += fTotalList->GetListOfBunch()[i]->GetWeight();
137        else nc += fTotalList->GetListOfBunch()[i]->GetWeight();
138    }
139
140// some prints
141#ifdef DEBUG
142    Msg(EsafMsg::Debug) << "SIZE OF LIST OF BUNCHES = "<<fNbBunch << MsgDispatch;
143    Msg(EsafMsg::Debug) << "Nb fluo = " <<nf << MsgDispatch;
144    Msg(EsafMsg::Debug) << "Nb ckov = " <<nc << MsgDispatch;
145
146    if ( fNbBunch >=1 ) {
147        EarthVector track = fTotalList->GetListOfBunch()[fNbBunch - 1]->GetPos() 
148                          - fTotalList->GetListOfBunch()[0]->GetPos();
149        Msg(EsafMsg::Debug) << "size of the photon track = "<<track.Mag()/km <<" km"<< MsgDispatch;
150    }
151#endif
152    Msg(EsafMsg::Info) << "TOTAL number of photons ="<<fNbTot << MsgDispatch;
153   
154   
155   
156
157// 3. root output init
158    EEvent* ev = EEvent::GetCurrent();
159    if(ev->GetAtmosphere()) ev->GetAtmosphere()->SetMaxScatOrder(0);
160
161    ////////////////////////////////////////////////////////////////////////////////////////////
162    //////////////// ALGO STARTS HERE FOR SINGLES created by LightSource ///////////////////////
163    ////////////////////////////////////////////////////////////////////////////////////////////
164    // add photons into fTemplist and remove them from fTotalList
165    for(size_t j=0; j < fTotalList->GetSingleEntries(); j++) {
166        fTempList.push_back(fTotalList->GetSingle());
167        fNbSingles++;
168    }
169    fTotalList->Clear(2);  // bunches list not cleared, but singles list is cleared (but pointers not destroyed)
170   
171    // scattering simu for single photons in fTempList
172    fPropagator->GoForPureMC(fTempList);
173
174    // Scan fTempList before filling fTotalList, applying relevant cuts :
175    //   - if photon status > 4, it is lost
176    //   - if upward-going photons : look if absorbed by ozone
177    ScanTempList();
178
179    // Move photons from fTempList to fTotalList : fTempList does not contain photons anymore
180    // Only photons which have successfully scattered are written in fTotalList, others are lost
181    fTotalList->Add(fTempList);
182
183    // save SinglePhotons after scattering simulation
184    if ( ev ) {
185        const vector<SinglePhoton*>& list_single_2 = fTotalList->GetListOfSingle();
186        size_t nbnotsaved = fTotalList->GetNbNotSaved();
187        for (size_t i = list_single_2.size() - nbnotsaved; i<list_single_2.size(); i++ ) {
188            EAtmosphereSingleAdder sa( list_single_2[i],true,false,0 );
189            ev->Fill(sa);
190            fTotalList->OneSingleSaved();
191        }
192    }
193
194
195
196    ////////////////////////////////////////////////////////////////////////////////////////////
197    ///////////////////////////////// ALGO STARTS HERE FOR BUNCHES /////////////////////////////
198    ////////////////////////////////////////////////////////////////////////////////////////////
199    BunchOfPhotons* bunch = fTotalList->GetBunch(0);
200    if(bunch) {
201        Double_t remainingPhotons = -HUGE;
202        Int_t split_level = Int_t(ceil(fNbTot / fMaxPhNb));
203        Msg(EsafMsg::Warning) << "NB of photons to propagate after reduction "<<fNbTot  << MsgDispatch;
204        Msg(EsafMsg::Warning) << "Each scattering order simu will be split in "<<split_level<<" parts" << MsgDispatch;
205
206        // saves bunch parameters at creation
207        if ( ev ) {
208            EAtmosphereBunchAdder ba( bunch, true );
209            ev->Fill(ba);
210        }
211
212
213        // 1. SCATTERING SIMULATION
214    #ifdef DEBUG
215        TString jobMCRT = "Time spent for Scattering simu:";
216        TBenchmark gB;
217        gB.Start(jobMCRT);
218    #endif
219        fTotalList->ResetBunchCounter();  // for below extraction from list
220        for(Int_t level=0; level < split_level; level++) {
221            Msg(EsafMsg::Info) << "************ Create single photons, Level : "<<level+1<<" ************"<< MsgDispatch;
222
223            // 1.1 generate single photons
224            remainingPhotons = ConvertIntoSingles(*bunch,remainingPhotons);
225
226            Msg(EsafMsg::Info) << "NB of SinglePhoton to PROPAGATE = "<<fNbSingles<< MsgDispatch;
227            fNbSingles = 0;
228
229
230            // 1.2 scattering simu for single photons in fTempList
231            fPropagator->GoForPureMC(fTempList);
232
233            // 1.3 Scan fTempList before filling fTotalList, applying relevant cuts :
234            //     - if photon status > 4, it is lost
235            ScanTempList();
236
237            // 2.4 Move photons from fTempList to fTotalList : fTempList does not contain photons anymore
238            //     Only photons which have successfully scattered are written in fTotalList, others are lost
239            fTotalList->Add(fTempList);
240
241    #ifdef DEBUG
242            gB.Stop(jobMCRT);
243            MsgForm(EsafMsg::Info,"Time spent for this level :  REAL=%6.2f s   CPU=%6.2f s",gB.GetRealTime(jobMCRT),gB.GetCpuTime(jobMCRT));
244    #endif
245        }
246
247
248        // 2. saves single photon parameters after scattering simulation
249        if ( ev ) {
250            const vector<SinglePhoton*>& list_single = fTotalList->GetListOfSingle();
251            size_t nbnotsaved = fTotalList->GetNbNotSaved();
252            for (size_t i = list_single.size() - nbnotsaved; i<list_single.size(); i++ ) {
253                EAtmosphereSingleAdder sa( list_single[i],true,true );
254                ev->Fill(sa);
255                fTotalList->OneSingleSaved();
256            }
257        }
258    }
259
260#ifdef DEBUG
261    Msg(EsafMsg::Debug) <<"final list of single, size = " <<fTotalList->GetListOfSingle().size() <<MsgDispatch;
262#endif
263    if(fTotalList->GetNbNotSaved()) Msg(EsafMsg::Warning) <<fTotalList->GetNbNotSaved()<<" SinglePhoton not saved into root" <<MsgDispatch;
264
265
266   
267    // DUMMY (saves single photons after propagation until detector)
268    if ( ev ) {
269        const vector<SinglePhoton*>& list_single_3 = fTotalList->GetListOfSingle();
270        for ( size_t i=0; i<list_single_3.size(); i++ ) {
271            EAtmosphereSingleAdder sa( list_single_3[i],false,false,i );
272            ev->Fill(sa);
273        }
274    }
275   
276
277    //FIXME: small mem leak here.
278    return new ListPhotonsOnPupil;
279}
280
281//_______________________________________________________________________________________________________________________
282Double_t PureMCRadiativeTransfer::ConvertIntoSingles(const BunchOfPhotons& b, Double_t remainingPhotons) {
283    //
284    // Photon features are set here (lambda, position, direction)
285    // SinglePhoton are stored into fTempList for scattering simu
286    //
287    // if remainingPhotons > -HUGE, means the nb of photons to create is specified
288    //
289
290    Int_t nb = 0;
291    Double_t rtn = -HUGE;
292    Double_t wl, date, tof;
293    EarthVector showerpos, dir, diff;
294    SinglePhoton* s = 0;
295    UInt_t bid = b.GetId();
296    PhotonType type = b.GetType();
297    tof = 0;
298   
299    if(remainingPhotons > -HUGE) nb = Int_t(remainingPhotons);
300    else nb = Int_t(b.GetWeight());
301   
302    for(Int_t i=0; i<nb; i++) {
303        // get features at creation
304        showerpos = b.RandomPosInShower();
305        if(fGround->IsUnderGround(showerpos)) continue;
306        diff = showerpos - b.GetShowerPos();
307        date = b.GetDate() + diff.Dot(b.GetDir().Unit())/Clight();
308        dir  = b.RandomDirection();
309        wl   = b.GetWlSpectrum().GetLambda();
310        s    = new SinglePhoton(type,date,tof,wl,showerpos,showerpos,dir,None,bid,b.GetShowerAge());
311        fNbSingles++;
312        fTempList.push_back(s);
313        if(fNbSingles >= fMaxPhNb) {
314            rtn = nb - 1 - i;
315            break;
316        }
317    }
318   
319    return rtn;
320}
321
322//_______________________________________________________________________________________________________________________
323void PureMCRadiativeTransfer::Reset() {
324    //
325    // get ready for next event
326    //
327    if(fGround) fGround->Reset();
328    if(fPropagator) fPropagator->Reset();
329    for(size_t i=0; i<fTempList.size(); i++) {
330        SafeDelete(fTempList[i]);
331        Msg(EsafMsg::Panic) <<"fTempList SHOULD BE EMPTY, remains : "<<fTempList.size()<<" photons inside"<<MsgDispatch;
332       
333    }
334    fTempList.clear();
335    fNbBunch = 0;
336    fNbTot = 0;
337    fNbSingles = 0;
338}
339
340//_______________________________________________________________________________________________________________________
341void PureMCRadiativeTransfer::ScanTempList() {
342    //
343    // remove unrelevant photons (status == LostByCut)
344    // and use Lowtran to handle Ozone absorption : calculated only for atmo outgoing photons
345    //
346   
347    // for future info dumping
348    Bool_t next(0),subnext(0);
349    Double_t nbTracked(0), fraction(0), left(0),right(10),subleft(0),subright(2.5);
350    size_t nbPh = fTempList.size();
351    Msg(EsafMsg::Info) << "Scan Temp List (remove unrelevant photons) : 0%" << MsgFlush;
352   
353    // init
354    SinglePhoton* p = 0;
355   
356    // scan temp list
357    for(UInt_t m=0; m < fTempList.size(); m++) {
358        p = fTempList[m];
359        if(p->Status() == LostByCut) SafeDelete(fTempList[m]);
360        nbTracked++;
361
362        // dump info
363        fraction = 100*nbTracked/Double_t(nbPh);
364        if (left<fraction && fraction <=right) 
365            next = kFALSE;
366        else if (fraction>right) {
367            next = kTRUE;
368            left = right;
369            right += 10;
370        }
371        if (subleft<fraction  && fraction <=subright) {
372            subnext = kFALSE;
373        }
374        if (fraction > subright) {
375            subnext = kTRUE;
376            subleft = subright;
377            subright +=2;
378        }
379        if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
380        if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;       
381    }
382    Msg(EsafMsg::Info) << " -done" << MsgDispatch;
383}
384
385
Note: See TracBrowser for help on using the repository browser.