source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/lighttoeuso/src/SlastLightToEuso.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.2 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: SlastLightToEuso.cc 2807 2008-10-12 17:13:07Z naumov $
3//  created Oct, 28 2002
4#include "SlastLightToEuso.hh"
5#include "EEventTruthAdder.hh"
6#include "Config.hh"
7#include "ShowerTrack.hh"
8#include "EShower.hh"
9#include "EShowerFiller.hh"
10#include "Atmosphere.hh"
11#include "EConst.hh"
12#include "DetectorGeometry.hh"
13#include <TVector3.h>
14using namespace sou;
15
16ClassImp(SlastLightToEuso)
17
18#ifdef ___USING_SLAST77___
19// interface to fortran functions
20extern "C" int slastinit_();
21extern "C" int slastsimulation_();
22extern "C" int slastdetectorgeometry_(Float_t*,Float_t*,Float_t*,Float_t*);
23extern "C" {
24    struct slast2esaf_header {
25        Float_t Energy;
26        Float_t Theta;
27        Float_t Phi;
28        Float_t X1;
29        Float_t Rint[3]; 
30        Float_t AgeEarth;
31        Float_t Rimpact[3];
32        Float_t Rmax[3];
33        Float_t Xmax;
34        Float_t ShowerHitGround;
35    };
36    struct slast2esaf_shower {
37        int sn;
38        Float_t ShowerXi[10000];
39        Float_t ShowerXf[10000];
40        Float_t ShowerRxi[10000];
41        Float_t ShowerRyi[10000];
42        Float_t ShowerRzi[10000];
43        Float_t ShowerRxf[10000];
44        Float_t ShowerRyf[10000];
45        Float_t ShowerRzf[10000];
46        Float_t ShowerTimei[10000];
47        Float_t ShowerTimef[10000];
48        Float_t ShowerAgei[10000];
49        Float_t ShowerAgef[10000];
50        Float_t ShowerNe[10000];
51    }; 
52    struct slast2esaf_fluoresence {
53        int fn;
54        Float_t fxshower[500000];
55        Float_t fyshower[500000];
56        Float_t fzshower[500000];
57        Float_t fxfocal[500000];
58        Float_t fyfocal[500000];
59        Float_t fzfocal[500000];
60        Float_t ftheta1[500000];
61        Float_t fphi1[500000];
62        Float_t flambda[500000];
63        Float_t ftime[500000];
64    };
65    struct slast2esaf_cherenkov {
66        int cn;
67        Float_t cxshower[500000];
68        Float_t cyshower[500000];
69        Float_t czshower[500000];
70        Float_t cxfocal[500000];
71        Float_t cyfocal[500000];
72        Float_t czfocal[500000];
73        Float_t ctheta1[500000];
74        Float_t cphi1[500000];
75        Float_t clambda[500000];
76        Float_t ctime[500000];
77    };
78    extern struct slast2esaf_header      esafheader_;
79    extern struct slast2esaf_shower      esafshower_;
80    extern struct slast2esaf_fluoresence esaffluor_;
81    extern struct slast2esaf_cherenkov   esafcherenkov_;
82}
83#endif
84//_________________________________________________________________________________________
85SlastLightToEuso::SlastLightToEuso() : LightToEuso("SLAST"),fPhotons(0),fTruth(0) {
86    //
87    // constructor
88    //
89    Init();
90}
91
92//_________________________________________________________________________________________
93void SlastLightToEuso::Init() {
94    //
95    // Enabling the atmosphere
96    //
97
98    Msg(EsafMsg::Info) << "Atmosphere enabled " << Atmosphere::Get()->GetType() << MsgDispatch;
99    PrepareDataCards();
100#ifdef ___USING_SLAST77___
101    slastinit_();
102#endif
103    fPhotons = new ListPhotonsOnPupil;
104
105}
106
107
108//_________________________________________________________________________________________
109SlastLightToEuso::~SlastLightToEuso() {
110    //
111    // destructor
112    //
113
114    SafeDelete(fTruth);
115    SafeDelete(fPhotons);
116}
117
118//______________________________________________________________________________
119Bool_t SlastLightToEuso::ClearMemory() {
120    //
121    //
122
123    return fPhotons ? fPhotons->ClearMemory() : kTRUE;
124}
125
126//_________________________________________________________________________________________
127PhotonsOnPupil *SlastLightToEuso::Get(const DetectorGeometry* dg) {
128    //
129    // returns the appropriate PhotonsOnPupil object for one event
130    //
131#ifdef ___USING_SLAST77___
132
133    // Slast banner
134    Msg(EsafMsg::Info) << "SLAST.f77 called" << MsgDispatch;
135    // call SLAST fortran code
136    Float_t ISSX = (dg->GetPos()).X()/km; 
137    Float_t ISSY = (dg->GetPos()).Y()/km; 
138    Float_t ISSZ = (dg->GetPos()).Z()/km; 
139    Float_t radius = dg->GetRadius()*mm/m;
140    slastdetectorgeometry_(&radius,&ISSX,&ISSY,&ISSZ);
141    slastsimulation_();
142    Msg(EsafMsg::Info) << "SLAST.f77 completed succesfully" << MsgDispatch;
143    // convert data to ListPhotonsOnPupil and return it
144    fPhotons->Clear();
145    // loop on photons coming from slast
146    double x[3];
147    double y[3];
148    double th=0.,ph=0.,w=0.,t=0.;
149    int fn = esaffluor_.fn;
150    if( esaffluor_.fn == 0 && esafcherenkov_.cn == 0) return fPhotons;
151    if (fn > 500000) fn = 500000;
152    // Filling fluorescence information from SLAST
153    for(int i=0; i<fn; i++) 
154    {
155        x[0] = esaffluor_.fxshower[i]*m; // in m
156        x[1] = esaffluor_.fyshower[i]*m; // in m
157        x[2] = esaffluor_.fzshower[i]*m; // in m
158        y[0] = esaffluor_.fxfocal[i]*m;
159        y[1] = esaffluor_.fyfocal[i]*m;
160        y[2] = esaffluor_.fzfocal[i]*m;
161        th   = esaffluor_.ftheta1[i]*deg;
162        ph   = esaffluor_.fphi1[i]*deg;
163        w    = esaffluor_.flambda[i]*nm ;
164        t    = esaffluor_.ftime[i]*microsecond;
165        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,w,t,(ParentPhotonType)1);
166        fPhotons->Add( p );
167    }
168    int cn = esafcherenkov_.cn;
169    if (cn > 500000) cn = 500000;
170    // Filling cherenkov information from SLAST
171    for(int i=0; i<cn; i++) 
172    {
173        x[0] = esafcherenkov_.cxshower[i]*m;
174        x[1] = esafcherenkov_.cyshower[i]*m;
175        x[2] = esafcherenkov_.czshower[i]*m;
176        y[0] = esafcherenkov_.cxfocal[i]*m;
177        y[1] = esafcherenkov_.cyfocal[i]*m;
178        y[2] = esafcherenkov_.czfocal[i]*m;
179        th   = esafcherenkov_.ctheta1[i]*deg;
180        ph   = esafcherenkov_.cphi1[i]*deg;
181        w    = esafcherenkov_.clambda[i]*nm ;
182        t    = esafcherenkov_.ctime[i]*microsecond;
183        ParentPhoton *p = new ParentPhoton(i+fn,x,y,th,ph,w,t,(ParentPhotonType)2);
184        fPhotons->Add( p );
185    }
186
187    MsgForm(EsafMsg::Info,"Fluorescent %d Cherenkov %d Total %d photons generated",fn,cn,fn+cn);
188    // set smontecarlo truth for the shower (clearing previous event first)
189    SafeDelete(fTruth);
190
191    EEvent* event = EEvent::GetCurrent();
192
193    if ( event ) {
194        EEventTruthAdder a( GetTruth());
195        event->Fill( a );
196        // Fill the Shower Track Object
197
198        ShowerTrack *track = new ShowerTrack();
199        ShowerStep s;
200        int sn = esafshower_.sn;
201        if (sn > 10000) sn = 10000;
202        for(int i=0;i<sn;i++) {
203            s.fXi = esafshower_.ShowerXi[i]*g/cm2;
204            s.fXf = esafshower_.ShowerXf[i]*g/cm2;
205            s.fXYZi.SetXYZ(esafshower_.ShowerRxi[i]*m,esafshower_.ShowerRyi[i]*m,esafshower_.ShowerRzi[i]*m);
206            s.fXYZf.SetXYZ(esafshower_.ShowerRxf[i]*m,esafshower_.ShowerRyf[i]*m,esafshower_.ShowerRzf[i]*m);
207            s.fTimei = esafshower_.ShowerTimei[i]*microsecond;
208            s.fTimef = esafshower_.ShowerTimef[i]*microsecond;
209            s.fAgei  = esafshower_.ShowerAgei[i];
210            s.fAgef  = esafshower_.ShowerAgef[i];
211            s.fNelectrons = esafshower_.ShowerNe[i];
212            track->Add(s);
213        }
214        track->fEnergy    = esafheader_.Energy;
215        track->fTheta     = esafheader_.Theta;
216        track->fPhi       = esafheader_.Phi;
217        track->fX1        = esafheader_.X1*g/cm2;
218        track->fEthrEl    = 0;
219        track->fDirVers.SetXYZ(TMath::Sin(esafheader_.Theta)*TMath::Cos(esafheader_.Phi),
220                TMath::Sin(esafheader_.Theta)*TMath::Sin(esafheader_.Phi),
221                TMath::Cos(esafheader_.Theta));
222        track->fInitPos.SetXYZ(esafheader_.Rint[0]*m,esafheader_.Rint[1]*m,esafheader_.Rint[2]*m);
223        track->fHitGround = (esafheader_.ShowerHitGround == 1 ? kTRUE: kFALSE);
224
225        EShowerFiller f( track );
226        event->Fill( f );
227        delete track;
228    }
229
230    system("rm -rf TAPE*");
231    return fPhotons;
232#else
233    return NULL;
234#endif
235}
236
237//_________________________________________________________________________________________
238void SlastLightToEuso::Configure() {
239    //
240    // re-configure SLAST
241    //
242}
243
244//_________________________________________________________________________________________
245MCTruth* SlastLightToEuso::GetTruth() {
246    //
247    // returns Montecarlo shower truth informations for this event
248    //
249#ifdef ___USING_SLAST77___
250    if (!fTruth) {
251        fTruth = new MCTruth();
252        fTruth->SetLightToEuso( this );
253        fTruth->SetEnergy( esafheader_.Energy*eV);
254        fTruth->SetThetaPhi(esafheader_.Theta,esafheader_.Phi);
255        EVector v;
256        v[0] = esafheader_.Rint[0]*m;
257        v[1] = esafheader_.Rint[1]*m;
258        v[2] = esafheader_.Rint[2]*m;
259        fTruth->SetFirstInt(v,esafheader_.X1*g/cm2);
260        v[0] = esafheader_.Rimpact[0]*m;
261        v[1] = esafheader_.Rimpact[1]*m;
262        v[2] = esafheader_.Rimpact[2]*m;
263        fTruth->SetEarthImpact(v,esafheader_.AgeEarth);
264        v[0] = esafheader_.Rmax[0]*m;
265        v[1] = esafheader_.Rmax[1]*m;
266        v[2] = esafheader_.Rmax[2]*m;
267        fTruth->SetShowerMax(v,esafheader_.Xmax*g/cm2);
268    }
269    return fTruth;
270#else
271    return NULL;
272#endif
273}
274
275//_________________________________________________________________________________________
276PhysicsData* SlastLightToEuso::GetPhysics() {
277    //
278    // returns simulated shower additional infos for this event
279    //
280    return NULL;
281}
282
283//_________________________________________________________________________________________
284void SlastLightToEuso::PrepareDataCards() {
285    //
286    // interface to fortran code
287    //
288    //
289    // Prepare data card to read by slast77 engine
290    //
291    ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
292    Float_t  EarthRadius    = EConst::EarthRadius()/km;
293   
294    // reading common to slast77 and slast++ options from GeneratorLightToEuso file   
295    pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
296    Float_t  FOV                = pConfig->GetNum("GeneratorLightToEuso.FoV"); 
297    Float_t  InteractionVectorX = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorX");
298    Float_t  InteractionVectorY = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorY");
299    Float_t  InteractionVectorZ = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorZ");
300    Float_t  ThetaRangeMin  = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMin");
301    Float_t  ThetaRangeMax  = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMax");
302    Float_t  PhiRangeMin    = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMin");
303    Float_t  PhiRangeMax    = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMax");
304    Float_t  EnergyRangeMin = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMin");
305    Float_t  EnergyRangeMax = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMax");
306    Float_t  EnergySlope    = pConfig->GetNum("GeneratorLightToEuso.EnergySlope");
307    string fType            = pConfig->GetStr("GeneratorLightToEuso.UhecrType");
308    Float_t UhecrType = EConst::AtomicMass(fType);
309   
310    // reading specific to slast77 options from its own config
311    Float_t  WaveRangeMin   =  Conf()->GetNum("SlastLightToEuso.WaveRangeMin");
312    Float_t  WaveRangeMax   =  Conf()->GetNum("SlastLightToEuso.WaveRangeMax");
313    string DoCherenkov      = Conf()->GetStr("SlastLightToEuso.DoCherenkov");
314    string DoFluorescence   = Conf()->GetStr("SlastLightToEuso.DoFluorescence");
315    string AtmType          = Conf()->GetStr("SlastLightToEuso.AtmosphericType");
316    Float_t  AtmTemperature = Conf()->GetNum("SlastLightToEuso.AtmTemperature");
317    Float_t  Albedo         = Conf()->GetNum("SlastLightToEuso.Albedo");
318    Float_t  GTU            = Conf()->GetNum("SlastLightToEuso.GTU");
319    string AtmCurvature     = Conf()->GetStr("SlastLightToEuso.AtmCurvature");
320    string EnergyDistributionParametrization = Conf()->GetStr("SlastLightToEuso.EnergyDistributionParametrization");
321    string ShowerParametrization = Conf()->GetStr("SlastLightToEuso.ShowerParametrization");
322   
323    FILE *datacard = fopen("auxilar/global.datacard","w");
324    int doch=0,dofl=0,atm_type=2,curvature=2,energydistribution=1,showerparametrization=1;
325    fprintf(datacard,"LIST\n");
326    fprintf(datacard,"%s%f\n","ERAD ",EarthRadius);
327    fprintf(datacard,"%s%f\n","FOV  ",FOV);
328    fprintf(datacard,"%s%f %f\n","WAVE ",WaveRangeMin,WaveRangeMax);
329    if(DoCherenkov == "yes")    doch=1;
330    if(DoFluorescence == "yes") dofl=1;
331    fprintf(datacard,"%s%d\n","DOCH ",doch);
332    fprintf(datacard,"%s%d\n","DOFL ",dofl);
333    if(AtmType == "Isothermic") atm_type=1;
334    if(AtmType == "USStandard") atm_type=2;
335    fprintf(datacard,"%s%d\n","ATMO ",atm_type);
336    fprintf(datacard,"%s%f\n","TEMP ",AtmTemperature);
337    fprintf(datacard,"%s%f\n","ALBE ",Albedo);
338    fprintf(datacard,"%s%f\n","GTU  ",GTU);
339    if(AtmCurvature == "Curved") curvature = 1;
340    if(AtmCurvature == "Planar") curvature = 2;
341    fprintf(datacard,"%s%d\n","CURV ",curvature);
342    fprintf(datacard,"%s%f %f %f\n","RINT ",InteractionVectorX,InteractionVectorY,InteractionVectorZ);
343    fprintf(datacard,"%s%f %f\n","THET ",ThetaRangeMin,ThetaRangeMax);
344    fprintf(datacard,"%s%f %f\n","PHI  ",PhiRangeMin,PhiRangeMax);
345    fprintf(datacard,"%s%e %e\n","EINT ",EnergyRangeMin,EnergyRangeMax);
346    fprintf(datacard,"%s%f\n","ERAN ",EnergySlope);
347    fprintf(datacard,"%s%f\n","TYPE ",UhecrType);
348    if(EnergyDistributionParametrization == "Hillas") energydistribution=1;
349    fprintf(datacard,"%s%d\n","DIST ",energydistribution);
350    if(ShowerParametrization == "GIL") showerparametrization=1;
351    fprintf(datacard,"%s%d\n","SHOW ",showerparametrization);
352    fclose(datacard);
353}
354
Note: See TracBrowser for help on using the repository browser.