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> |
---|
14 | using namespace sou; |
---|
15 | |
---|
16 | ClassImp(SlastLightToEuso) |
---|
17 | |
---|
18 | #ifdef ___USING_SLAST77___ |
---|
19 | // interface to fortran functions |
---|
20 | extern "C" int slastinit_(); |
---|
21 | extern "C" int slastsimulation_(); |
---|
22 | extern "C" int slastdetectorgeometry_(Float_t*,Float_t*,Float_t*,Float_t*); |
---|
23 | extern "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 | //_________________________________________________________________________________________ |
---|
85 | SlastLightToEuso::SlastLightToEuso() : LightToEuso("SLAST"),fPhotons(0),fTruth(0) { |
---|
86 | // |
---|
87 | // constructor |
---|
88 | // |
---|
89 | Init(); |
---|
90 | } |
---|
91 | |
---|
92 | //_________________________________________________________________________________________ |
---|
93 | void 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 | //_________________________________________________________________________________________ |
---|
109 | SlastLightToEuso::~SlastLightToEuso() { |
---|
110 | // |
---|
111 | // destructor |
---|
112 | // |
---|
113 | |
---|
114 | SafeDelete(fTruth); |
---|
115 | SafeDelete(fPhotons); |
---|
116 | } |
---|
117 | |
---|
118 | //______________________________________________________________________________ |
---|
119 | Bool_t SlastLightToEuso::ClearMemory() { |
---|
120 | // |
---|
121 | // |
---|
122 | |
---|
123 | return fPhotons ? fPhotons->ClearMemory() : kTRUE; |
---|
124 | } |
---|
125 | |
---|
126 | //_________________________________________________________________________________________ |
---|
127 | PhotonsOnPupil *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 | //_________________________________________________________________________________________ |
---|
238 | void SlastLightToEuso::Configure() { |
---|
239 | // |
---|
240 | // re-configure SLAST |
---|
241 | // |
---|
242 | } |
---|
243 | |
---|
244 | //_________________________________________________________________________________________ |
---|
245 | MCTruth* 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 | //_________________________________________________________________________________________ |
---|
276 | PhysicsData* SlastLightToEuso::GetPhysics() { |
---|
277 | // |
---|
278 | // returns simulated shower additional infos for this event |
---|
279 | // |
---|
280 | return NULL; |
---|
281 | } |
---|
282 | |
---|
283 | //_________________________________________________________________________________________ |
---|
284 | void 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 | |
---|