]>
Commit | Line | Data |
---|---|---|
03ca248b | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* | |
17 | $Log$ | |
a7e55c80 | 18 | Revision 1.8 2002/12/06 12:28:44 morsch |
19 | Region to media mapping corrected and improved. | |
20 | ||
6d4d27f2 | 21 | Revision 1.7 2002/12/06 12:21:32 morsch |
22 | User stepping methods added (E. Futo) | |
23 | ||
fa3d1cc7 | 24 | Revision 1.6 2002/11/21 18:40:06 iglez2 |
25 | Media handling added | |
26 | ||
27b2f7fe | 27 | Revision 1.5 2002/11/07 17:59:10 iglez2 |
28 | Included the geometry through geant4_vmc/FLUGG | |
29 | ||
bf3aa28e | 30 | Revision 1.4 2002/11/04 16:00:46 iglez2 |
31 | The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent. | |
32 | ||
f9cb2fec | 33 | Revision 1.3 2002/10/22 15:12:14 alibrary |
34 | Introducing Riostream.h | |
35 | ||
eae0fe66 | 36 | Revision 1.2 2002/10/14 14:57:40 hristov |
37 | Merging the VirtualMC branch to the main development branch (HEAD) | |
38 | ||
b9d0a01d | 39 | Revision 1.1.2.8 2002/10/08 16:33:17 iglez2 |
40 | LSOUIT is set to true before the second call to flukam. | |
41 | ||
42 | Revision 1.1.2.7 2002/10/08 09:30:37 iglez2 | |
43 | Solved stupid missing ; | |
44 | ||
45 | Revision 1.1.2.6 2002/10/07 13:40:22 iglez2 | |
46 | First implementations of the PDG <--> Fluka Id conversion routines | |
47 | ||
48 | Revision 1.1.2.5 2002/09/26 16:26:03 iglez2 | |
49 | Added verbosity | |
50 | Call to gAlice->Generator()->Generate() | |
51 | ||
52 | Revision 1.1.2.4 2002/09/26 13:22:23 iglez2 | |
53 | Naive implementation of ProcessRun and ProcessEvent | |
54 | Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init() | |
55 | ||
56 | Revision 1.1.2.3 2002/09/20 15:35:51 iglez2 | |
57 | Modification of LFDRTR. Value is passed to FLUKA !!! | |
58 | ||
59 | Revision 1.1.2.2 2002/09/18 14:34:44 iglez2 | |
60 | Revised version with all pure virtual methods implemented | |
61 | ||
62 | Revision 1.1.2.1 2002/07/24 08:49:41 alibrary | |
63 | Adding TFluka to VirtualMC | |
64 | ||
65 | Revision 1.1 2002/07/05 13:10:07 morsch | |
66 | First commit of Fluka interface. | |
67 | ||
03ca248b | 68 | */ |
69 | ||
eae0fe66 | 70 | #include <Riostream.h> |
b9d0a01d | 71 | |
6d4d27f2 | 72 | #include "TClonesArray.h" |
03ca248b | 73 | #include "TFluka.h" |
b9d0a01d | 74 | #include "TCallf77.h" //For the fortran calls |
75 | #include "Fdblprc.h" //(DBLPRC) fluka common | |
b9d0a01d | 76 | #include "Fepisor.h" //(EPISOR) fluka common |
fa3d1cc7 | 77 | #include "Ffinuc.h" //(FINUC) fluka common |
78 | #include "Fiounit.h" //(IOUNIT) fluka common | |
79 | #include "Fpaprop.h" //(PAPROP) fluka common | |
f9cb2fec | 80 | #include "Fpart.h" //(PART) fluka common |
fa3d1cc7 | 81 | #include "Ftrackr.h" //(TRACKR) fluka common |
6d4d27f2 | 82 | #include "Fpaprop.h" //(PAPROP) fluka common |
fa3d1cc7 | 83 | #include "Ffheavy.h" //(FHEAVY) fluka common |
b9d0a01d | 84 | |
fa3d1cc7 | 85 | #include "TVirtualMC.h" |
bf3aa28e | 86 | #include "TG4GeometryManager.h" //For the geometry management |
87 | #include "TG4DetConstruction.h" //For the detector construction | |
88 | ||
89 | #include "FGeometryInit.hh" | |
fa3d1cc7 | 90 | #include "TLorentzVector.h" |
6d4d27f2 | 91 | #include "FlukaVolume.h" |
bf3aa28e | 92 | |
b9d0a01d | 93 | // Fluka methods that may be needed. |
94 | #ifndef WIN32 | |
95 | # define flukam flukam_ | |
96 | # define fluka_openinp fluka_openinp_ | |
97 | # define fluka_closeinp fluka_closeinp_ | |
f9cb2fec | 98 | # define mcihad mcihad_ |
99 | # define mpdgha mpdgha_ | |
b9d0a01d | 100 | #else |
101 | # define flukam FLUKAM | |
102 | # define fluka_openinp FLUKA_OPENINP | |
103 | # define fluka_closeinp FLUKA_CLOSEINP | |
f9cb2fec | 104 | # define mcihad MCIHAD |
105 | # define mpdgha MPDGHA | |
b9d0a01d | 106 | #endif |
107 | ||
108 | extern "C" | |
109 | { | |
110 | // | |
111 | // Prototypes for FLUKA functions | |
112 | // | |
113 | void type_of_call flukam(const int&); | |
114 | void type_of_call fluka_openinp(const int&, DEFCHARA); | |
115 | void type_of_call fluka_closeinp(const int&); | |
f9cb2fec | 116 | int type_of_call mcihad(const int&); |
117 | int type_of_call mpdgha(const int&); | |
b9d0a01d | 118 | } |
119 | ||
120 | // | |
121 | // Class implementation for ROOT | |
122 | // | |
03ca248b | 123 | ClassImp(TFluka) |
b9d0a01d | 124 | |
125 | // | |
bf3aa28e | 126 | //---------------------------------------------------------------------------- |
127 | // TFluka constructors and destructors. | |
b9d0a01d | 128 | //____________________________________________________________________________ |
129 | TFluka::TFluka() | |
130 | :TVirtualMC(), | |
131 | fVerbosityLevel(0), | |
bf3aa28e | 132 | fInputFileName(""), |
27b2f7fe | 133 | fDetector(0), |
134 | fCurrentFlukaRegion(-1) | |
b9d0a01d | 135 | { |
136 | // | |
137 | // Default constructor | |
138 | // | |
139 | } | |
140 | ||
b9d0a01d | 141 | TFluka::TFluka(const char *title, Int_t verbosity) |
142 | :TVirtualMC("TFluka",title), | |
143 | fVerbosityLevel(verbosity), | |
bf3aa28e | 144 | fInputFileName(""), |
27b2f7fe | 145 | fDetector(0), |
146 | fCurrentFlukaRegion(-1) | |
b9d0a01d | 147 | { |
148 | if (fVerbosityLevel >=3) | |
149 | cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl; | |
150 | ||
bf3aa28e | 151 | |
152 | // create geometry manager | |
153 | if (fVerbosityLevel >=2) | |
154 | cout << "\t* Creating G4 Geometry manager..." << endl; | |
155 | fGeometryManager = new TG4GeometryManager(); | |
156 | if (fVerbosityLevel >=2) | |
157 | cout << "\t* Creating G4 Detector..." << endl; | |
158 | fDetector = new TG4DetConstruction(); | |
159 | FGeometryInit* geominit = FGeometryInit::GetInstance(); | |
160 | if (geominit) | |
161 | geominit->setDetConstruction(fDetector); | |
162 | else { | |
163 | cerr << "ERROR: Could not create FGeometryInit!" << endl; | |
164 | cerr << " Exiting!!!" << endl; | |
165 | abort(); | |
166 | } | |
167 | ||
b9d0a01d | 168 | if (fVerbosityLevel >=3) |
169 | cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl; | |
6d4d27f2 | 170 | |
171 | fVolumeMediaMap = new TClonesArray("FlukaVolume",1000); | |
172 | fNVolumes = 0; | |
173 | fMediaByRegion = 0; | |
b9d0a01d | 174 | } |
175 | ||
bf3aa28e | 176 | TFluka::~TFluka() { |
177 | if (fVerbosityLevel >=3) | |
178 | cout << "==> TFluka::~TFluka() destructor called." << endl; | |
179 | ||
180 | delete fGeometryManager; | |
6d4d27f2 | 181 | fVolumeMediaMap->Delete(); |
182 | delete fVolumeMediaMap; | |
183 | ||
bf3aa28e | 184 | |
185 | if (fVerbosityLevel >=3) | |
186 | cout << "<== TFluka::~TFluka() destructor called." << endl; | |
187 | } | |
188 | ||
189 | // | |
190 | //_____________________________________________________________________________ | |
191 | // TFluka control methods | |
b9d0a01d | 192 | //____________________________________________________________________________ |
193 | void TFluka::Init() { | |
194 | if (fVerbosityLevel >=3) | |
195 | cout << "==> TFluka::Init() called." << endl; | |
196 | ||
197 | if (fVerbosityLevel >=2) | |
198 | cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F') | |
199 | << ") in fluka..." << endl; | |
200 | GLOBAL.lfdrtr = true; | |
201 | ||
202 | if (fVerbosityLevel >=2) | |
203 | cout << "\t* Opening file " << fInputFileName << endl; | |
204 | const char* fname = fInputFileName; | |
205 | fluka_openinp(lunin, PASSCHARA(fname)); | |
206 | ||
207 | if (fVerbosityLevel >=2) | |
208 | cout << "\t* Calling flukam..." << endl; | |
bf3aa28e | 209 | flukam(1); |
b9d0a01d | 210 | |
211 | if (fVerbosityLevel >=2) | |
212 | cout << "\t* Closing file " << fInputFileName << endl; | |
213 | fluka_closeinp(lunin); | |
214 | ||
215 | if (fVerbosityLevel >=3) | |
216 | cout << "<== TFluka::Init() called." << endl; | |
fa3d1cc7 | 217 | |
6d4d27f2 | 218 | FinishGeometry(); |
219 | ||
b9d0a01d | 220 | } |
221 | ||
bf3aa28e | 222 | void TFluka::FinishGeometry() { |
6d4d27f2 | 223 | // |
224 | // Build-up table with region to medium correspondance | |
225 | // | |
226 | char tmp[5]; | |
227 | ||
bf3aa28e | 228 | if (fVerbosityLevel >=3) |
229 | cout << "==> TFluka::FinishGeometry() called." << endl; | |
230 | ||
6d4d27f2 | 231 | // fGeometryManager->Ggclos(); |
bf3aa28e | 232 | |
6d4d27f2 | 233 | FGeometryInit* flugg = FGeometryInit::GetInstance(); |
234 | ||
235 | fMediaByRegion = new Int_t[fNVolumes+2]; | |
236 | for (Int_t i = 0; i < fNVolumes; i++) | |
237 | { | |
238 | FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]); | |
239 | TString volName = vol->GetName(); | |
240 | Int_t media = vol->GetMedium(); | |
241 | printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media); | |
242 | strcpy(tmp, volName.Data()); | |
243 | tmp[4] = '\0'; | |
244 | flugg->SetMediumFromName(tmp, media); | |
27b2f7fe | 245 | } |
6d4d27f2 | 246 | |
247 | flugg->BuildMediaMap(); | |
27b2f7fe | 248 | |
bf3aa28e | 249 | if (fVerbosityLevel >=3) |
250 | cout << "<== TFluka::FinishGeometry() called." << endl; | |
251 | } | |
252 | ||
253 | void TFluka::BuildPhysics() { | |
254 | if (fVerbosityLevel >=3) | |
255 | cout << "==> TFluka::BuildPhysics() called." << endl; | |
256 | ||
257 | ||
258 | if (fVerbosityLevel >=3) | |
259 | cout << "<== TFluka::BuildPhysics() called." << endl; | |
260 | } | |
261 | ||
b9d0a01d | 262 | void TFluka::ProcessEvent() { |
263 | if (fVerbosityLevel >=3) | |
264 | cout << "==> TFluka::ProcessEvent() called." << endl; | |
265 | ||
266 | if (fVerbosityLevel >=3) | |
267 | cout << "<== TFluka::ProcessEvent() called." << endl; | |
268 | } | |
269 | ||
bf3aa28e | 270 | |
b9d0a01d | 271 | void TFluka::ProcessRun(Int_t nevent) { |
272 | if (fVerbosityLevel >=3) | |
273 | cout << "==> TFluka::ProcessRun(" << nevent << ") called." | |
274 | << endl; | |
275 | ||
276 | if (fVerbosityLevel >=2) { | |
277 | cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl; | |
278 | cout << "\t* Calling flukam again..." << endl; | |
279 | } | |
280 | fApplication->GeneratePrimaries(); | |
281 | EPISOR.lsouit = true; | |
bf3aa28e | 282 | flukam(1); |
b9d0a01d | 283 | |
284 | if (fVerbosityLevel >=3) | |
285 | cout << "<== TFluka::ProcessRun(" << nevent << ") called." | |
286 | << endl; | |
287 | } | |
288 | ||
bf3aa28e | 289 | //_____________________________________________________________________________ |
290 | // methods for building/management of geometry | |
291 | //____________________________________________________________________________ | |
292 | // functions from GCONS | |
293 | void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z, | |
294 | Float_t &dens, Float_t &radl, Float_t &absl, | |
295 | Float_t* ubuf, Int_t& nbuf) { | |
296 | // | |
297 | fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf); | |
298 | } | |
299 | ||
300 | void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z, | |
301 | Double_t &dens, Double_t &radl, Double_t &absl, | |
302 | Double_t* ubuf, Int_t& nbuf) { | |
303 | // | |
304 | fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf); | |
305 | } | |
306 | ||
307 | // detector composition | |
308 | void TFluka::Material(Int_t& kmat, const char* name, Double_t a, | |
309 | Double_t z, Double_t dens, Double_t radl, Double_t absl, | |
310 | Float_t* buf, Int_t nwbuf) { | |
311 | // | |
312 | fGeometryManager | |
313 | ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); | |
314 | } | |
315 | void TFluka::Material(Int_t& kmat, const char* name, Double_t a, | |
316 | Double_t z, Double_t dens, Double_t radl, Double_t absl, | |
317 | Double_t* buf, Int_t nwbuf) { | |
318 | // | |
319 | fGeometryManager | |
320 | ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); | |
321 | } | |
322 | ||
323 | void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, | |
324 | Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) { | |
325 | // | |
326 | fGeometryManager | |
327 | ->Mixture(kmat, name, a, z, dens, nlmat, wmat); | |
328 | } | |
329 | void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, | |
330 | Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) { | |
331 | // | |
332 | fGeometryManager | |
333 | ->Mixture(kmat, name, a, z, dens, nlmat, wmat); | |
334 | } | |
335 | ||
336 | void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, | |
337 | Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, | |
338 | Double_t stemax, Double_t deemax, Double_t epsil, | |
339 | Double_t stmin, Float_t* ubuf, Int_t nbuf) { | |
340 | // | |
341 | fGeometryManager | |
342 | ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, | |
343 | epsil, stmin, ubuf, nbuf); | |
344 | } | |
345 | void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, | |
346 | Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, | |
347 | Double_t stemax, Double_t deemax, Double_t epsil, | |
348 | Double_t stmin, Double_t* ubuf, Int_t nbuf) { | |
349 | // | |
350 | fGeometryManager | |
351 | ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, | |
352 | epsil, stmin, ubuf, nbuf); | |
353 | } | |
354 | ||
355 | void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, | |
356 | Double_t thetaY, Double_t phiY, Double_t thetaZ, | |
357 | Double_t phiZ) { | |
358 | // | |
359 | fGeometryManager | |
360 | ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); | |
361 | } | |
362 | ||
363 | void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) { | |
364 | // | |
365 | fGeometryManager->Gstpar(itmed, param, parval); | |
366 | } | |
367 | ||
368 | // functions from GGEOM | |
369 | Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed, | |
370 | Float_t *upar, Int_t np) { | |
371 | // | |
6d4d27f2 | 372 | // fVolumeMediaMap[TString(name)] = nmed; |
373 | TClonesArray &lvols = *fVolumeMediaMap; | |
374 | new(lvols[fNVolumes++]) | |
375 | FlukaVolume(name, nmed); | |
376 | return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); | |
bf3aa28e | 377 | } |
378 | Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed, | |
379 | Double_t *upar, Int_t np) { | |
380 | // | |
6d4d27f2 | 381 | TClonesArray &lvols = *fVolumeMediaMap; |
382 | new(lvols[fNVolumes++]) | |
383 | FlukaVolume(name, nmed); | |
384 | ||
385 | return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); | |
bf3aa28e | 386 | } |
387 | ||
388 | void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, | |
389 | Int_t iaxis) { | |
390 | // | |
6d4d27f2 | 391 | fGeometryManager->Gsdvn(name, mother, ndiv, iaxis); |
bf3aa28e | 392 | } |
393 | ||
394 | void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, | |
395 | Int_t iaxis, Double_t c0i, Int_t numed) { | |
396 | // | |
6d4d27f2 | 397 | TClonesArray &lvols = *fVolumeMediaMap; |
398 | new(lvols[fNVolumes++]) | |
399 | FlukaVolume(name, numed); | |
400 | fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); | |
bf3aa28e | 401 | } |
402 | ||
403 | void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, | |
404 | Int_t iaxis, Int_t numed, Int_t ndvmx) { | |
6d4d27f2 | 405 | // |
406 | TClonesArray &lvols = *fVolumeMediaMap; | |
407 | new(lvols[fNVolumes++]) | |
408 | FlukaVolume(name, numed); | |
409 | fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx); | |
bf3aa28e | 410 | } |
411 | ||
412 | void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, | |
413 | Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { | |
414 | // | |
6d4d27f2 | 415 | TClonesArray &lvols = *fVolumeMediaMap; |
416 | new(lvols[fNVolumes++]) | |
417 | FlukaVolume(name, numed); | |
418 | fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); | |
bf3aa28e | 419 | } |
420 | ||
421 | void TFluka::Gsord(const char *name, Int_t iax) { | |
422 | // | |
423 | fGeometryManager->Gsord(name, iax); | |
424 | } | |
425 | ||
426 | void TFluka::Gspos(const char *name, Int_t nr, const char *mother, | |
427 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
428 | const char *konly) { | |
429 | // | |
430 | fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly); | |
431 | } | |
432 | ||
433 | void TFluka::Gsposp(const char *name, Int_t nr, const char *mother, | |
434 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
435 | const char *konly, Float_t *upar, Int_t np) { | |
436 | // | |
437 | fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); | |
438 | } | |
439 | void TFluka::Gsposp(const char *name, Int_t nr, const char *mother, | |
440 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
441 | const char *konly, Double_t *upar, Int_t np) { | |
442 | // | |
443 | fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); | |
444 | } | |
445 | ||
446 | void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) { | |
447 | // | |
448 | fGeometryManager->Gsbool(onlyVolName, manyVolName); | |
449 | } | |
450 | ||
451 | void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov, | |
452 | Float_t *absco, Float_t *effic, Float_t *rindex) { | |
453 | // | |
454 | fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex); | |
455 | } | |
456 | void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov, | |
457 | Double_t *absco, Double_t *effic, Double_t *rindex) { | |
458 | // | |
459 | fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex); | |
460 | } | |
461 | ||
462 | // Euclid | |
463 | void TFluka::WriteEuclid(const char* fileName, const char* topVol, | |
464 | Int_t number, Int_t nlevel) { | |
465 | // | |
466 | fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel); | |
467 | } | |
468 | ||
469 | ||
470 | ||
27b2f7fe | 471 | //_____________________________________________________________________________ |
472 | // methods needed by the stepping | |
473 | //____________________________________________________________________________ | |
6d4d27f2 | 474 | |
27b2f7fe | 475 | Int_t TFluka::GetMedium() const { |
6d4d27f2 | 476 | FGeometryInit* flugg = FGeometryInit::GetInstance(); |
477 | return flugg->GetMedium(fCurrentFlukaRegion); | |
27b2f7fe | 478 | } |
bf3aa28e | 479 | |
480 | ||
481 | ||
482 | //____________________________________________________________________________ | |
483 | // ID <--> PDG transformations | |
b9d0a01d | 484 | //_____________________________________________________________________________ |
485 | Int_t TFluka::IdFromPDG(Int_t pdg) const | |
486 | { | |
487 | // | |
f9cb2fec | 488 | // Return Fluka code from PDG and pseudo ENDF code |
489 | ||
490 | // MCIHAD() goes from pdg to fluka internal. | |
491 | Int_t intfluka = mcihad(pdg); | |
492 | // KPTOIP array goes from internal to official | |
493 | return GetFlukaKPTOIP(intfluka); | |
b9d0a01d | 494 | } |
495 | ||
b9d0a01d | 496 | Int_t TFluka::PDGFromId(Int_t id) const |
497 | { | |
498 | // | |
f9cb2fec | 499 | // Return PDG code and pseudo ENDF code from Fluka code |
b9d0a01d | 500 | |
f9cb2fec | 501 | //IPTOKP array goes from official to internal |
502 | Int_t intfluka = GetFlukaIPTOKP(id); | |
503 | //MPKDHA() goes from internal to PDG | |
504 | return mpdgha(intfluka); | |
6d4d27f2 | 505 | } |
506 | ||
fa3d1cc7 | 507 | //_____________________________________________________________________________ |
508 | // methods for step management | |
509 | //____________________________________________________________________________ | |
510 | // | |
511 | // dynamic properties | |
512 | // | |
513 | void TFluka::TrackPosition(TLorentzVector& position) const | |
514 | { | |
515 | // Return the current position in the master reference frame of the | |
516 | // track being transported | |
517 | // TRACKR.atrack = age of the particle | |
518 | // TRACKR.xtrack = x-position of the last point | |
519 | // TRACKR.ytrack = y-position of the last point | |
520 | // TRACKR.ztrack = z-position of the last point | |
521 | position.SetX(TRACKR.xtrack[TRACKR.ntrack]); | |
522 | position.SetY(TRACKR.ytrack[TRACKR.ntrack]); | |
523 | position.SetZ(TRACKR.ztrack[TRACKR.ntrack]); | |
524 | position.SetT(TRACKR.atrack); | |
525 | } | |
526 | ||
527 | void TFluka::TrackMomentum(TLorentzVector& momentum) const | |
528 | { | |
529 | // Return the direction and the momentum (GeV/c) of the track | |
530 | // currently being transported | |
531 | // TRACKR.ptrack = momentum of the particle (not always defined, if | |
532 | // < 0 must be obtained from etrack) | |
533 | // TRACKR.cx,y,ztrck = direction cosines of the current particle | |
534 | // TRACKR.etrack = total energy of the particle | |
535 | // TRACKR.jtrack = identity number of the particle | |
536 | // PAPROP.am[TRACKR.jtrack] = particle mass in gev | |
537 | if (TRACKR.ptrack >= 0) { | |
538 | momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck); | |
539 | momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck); | |
540 | momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck); | |
541 | momentum.SetE(TRACKR.etrack); | |
542 | return; | |
543 | } | |
544 | else { | |
545 | Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]); | |
546 | momentum.SetPx(p*TRACKR.cxtrck); | |
547 | momentum.SetPy(p*TRACKR.cytrck); | |
548 | momentum.SetPz(p*TRACKR.cztrck); | |
549 | momentum.SetE(TRACKR.etrack); | |
550 | return; | |
551 | } | |
552 | } | |
553 | ||
554 | Double_t TFluka::TrackStep() const | |
555 | { | |
556 | // Return the length in centimeters of the current step | |
557 | // TRACKR.ctrack = total curved path | |
558 | return TRACKR.ctrack; | |
559 | } | |
560 | ||
561 | Double_t TFluka::TrackLength() const | |
562 | { | |
563 | // It is wrong | |
564 | // should be the sum of all steps starting from the beginning of the track | |
565 | // for the time being returns only the length in centimeters of the current step | |
566 | return TRACKR.ctrack; | |
567 | } | |
568 | ||
569 | Double_t TFluka::TrackTime() const | |
570 | { | |
571 | // Return the current time of flight of the track being transported | |
572 | // TRACKR.atrack = age of the particle | |
573 | return TRACKR.atrack; | |
574 | } | |
575 | ||
576 | Double_t TFluka::Edep() const | |
577 | { | |
578 | // Energy deposition | |
579 | // if TRACKR.ntrack = 0, TRACKR.mtrack = 0: | |
580 | // -->local energy deposition (the value and the point are not recorded in TRACKR) | |
581 | // but in the variable "rull" of the procedure "endraw.cxx" | |
582 | // if TRACKR.ntrack > 0, TRACKR.mtrack = 0: | |
583 | // -->no energy loss along the track | |
584 | // if TRACKR.ntrack > 0, TRACKR.mtrack > 0: | |
585 | // -->energy loss distributed along the track | |
586 | // TRACKR.dtrack = energy deposition of the jth deposition even | |
587 | if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0) | |
588 | return fRull; | |
589 | else { | |
590 | Double_t sum = 0; | |
591 | for ( Int_t j=0;j<TRACKR.mtrack;j++) { | |
592 | sum +=TRACKR.dtrack[j]; | |
593 | } | |
594 | return sum; | |
595 | } | |
596 | } | |
597 | ||
598 | Int_t TFluka::TrackPid() const | |
599 | { | |
600 | // Return the id of the particle transported | |
601 | // TRACKR.jtrack = identity number of the particle | |
602 | return PDGFromId(TRACKR.jtrack); | |
603 | } | |
604 | ||
605 | Double_t TFluka::TrackCharge() const | |
606 | { | |
607 | // Return charge of the track currently transported | |
608 | // PAPROP.ichrge = electric charge of the particle | |
609 | return PAPROP.ichrge[TRACKR.jtrack+6]; | |
610 | } | |
611 | ||
612 | Double_t TFluka::TrackMass() const | |
613 | { | |
614 | // PAPROP.am = particle mass in GeV | |
615 | return PAPROP.am[TRACKR.jtrack+6]; | |
616 | } | |
617 | ||
618 | Double_t TFluka::Etot() const | |
619 | { | |
620 | // TRACKR.etrack = total energy of the particle | |
621 | return TRACKR.etrack; | |
622 | } | |
623 | ||
624 | // | |
625 | // track status | |
626 | // | |
627 | Bool_t TFluka::IsNewTrack() const | |
628 | { | |
629 | // ???????????????, | |
630 | // True if the track is not at the boundary of the current volume | |
631 | // Not true in some cases in bxdraw - to be solved | |
632 | return 1; | |
633 | } | |
634 | ||
635 | Bool_t TFluka::IsTrackInside() const | |
636 | { | |
637 | // True if the track is not at the boundary of the current volume | |
638 | // In Fluka a step is always inside one kind of material | |
639 | // If the step would go behind the region of one material, | |
640 | // it will be shortened to reach only the boundary. | |
641 | // Therefore IsTrackInside() is always true. | |
642 | // Not true in some cases in bxdraw - to be solved | |
643 | return 1; | |
644 | } | |
645 | ||
646 | Bool_t TFluka::IsTrackEntering() const | |
647 | { | |
648 | // True if this is the first step of the track in the current volume | |
649 | // Boundary- (X) crossing | |
650 | // Icode = 19: boundary crossing - call from Kaskad | |
651 | // Icode = 29: boundary crossing - call from Emfsco | |
652 | // Icode = 39: boundary crossing - call from Kasneu | |
653 | // Icode = 49: boundary crossing - call from Kashea | |
654 | // Icode = 59: boundary crossing - call from Kasoph | |
655 | if (fIcode == 19 || | |
656 | fIcode == 29 || | |
657 | fIcode == 39 || | |
658 | fIcode == 49 || | |
659 | fIcode == 59) return 1; | |
660 | else return 0; | |
661 | } | |
662 | ||
663 | Bool_t TFluka::IsTrackExiting() const | |
664 | { | |
665 | // True if this is the last step of the track in the current volume | |
666 | // Boundary- (X) crossing | |
667 | // Icode = 19: boundary crossing - call from Kaskad | |
668 | // Icode = 29: boundary crossing - call from Emfsco | |
669 | // Icode = 39: boundary crossing - call from Kasneu | |
670 | // Icode = 49: boundary crossing - call from Kashea | |
671 | // Icode = 59: boundary crossing - call from Kasoph | |
672 | if (fIcode == 19 || | |
673 | fIcode == 29 || | |
674 | fIcode == 39 || | |
675 | fIcode == 49 || | |
676 | fIcode == 59) return 1; | |
677 | else return 0; | |
678 | } | |
679 | ||
680 | Bool_t TFluka::IsTrackOut() const | |
681 | { | |
682 | // True if the track is out of the setup | |
683 | // means escape | |
684 | // Icode = 14: escape - call from Kaskad | |
685 | // Icode = 23: escape - call from Emfsco | |
686 | // Icode = 32: escape - call from Kasneu | |
687 | // Icode = 40: escape - call from Kashea | |
688 | // Icode = 51: escape - call from Kasoph | |
689 | if (fIcode == 14 || | |
690 | fIcode == 23 || | |
691 | fIcode == 32 || | |
692 | fIcode == 40 || | |
693 | fIcode == 51) return 1; | |
694 | else return 0; | |
695 | } | |
696 | ||
697 | Bool_t TFluka::IsTrackDisappeared() const | |
698 | { | |
699 | // means all inelastic interactions and decays | |
700 | // fIcode from usdraw | |
701 | if (fIcode == 101 || // inelastic interaction | |
702 | fIcode == 102 || // particle decay | |
703 | fIcode == 214 || // in-flight annihilation | |
704 | fIcode == 215 || // annihilation at rest | |
705 | fIcode == 217 || // pair production | |
706 | fIcode == 221) return 1; | |
707 | else return 0; | |
708 | } | |
709 | ||
710 | Bool_t TFluka::IsTrackStop() const | |
711 | { | |
712 | // True if the track energy has fallen below the threshold | |
713 | // means stopped by signal or below energy threshold | |
714 | // Icode = 12: stopping particle - call from Kaskad | |
715 | // Icode = 15: time kill - call from Kaskad | |
716 | // Icode = 21: below threshold, iarg=1 - call from Emfsco | |
717 | // Icode = 22: below threshold, iarg=2 - call from Emfsco | |
718 | // Icode = 24: time kill - call from Emfsco | |
719 | // Icode = 31: below threshold - call from Kasneu | |
720 | // Icode = 33: time kill - call from Kasneu | |
721 | // Icode = 41: time kill - call from Kashea | |
722 | // Icode = 52: time kill - call from Kasoph | |
723 | if (fIcode == 12 || | |
724 | fIcode == 15 || | |
725 | fIcode == 21 || | |
726 | fIcode == 22 || | |
727 | fIcode == 24 || | |
728 | fIcode == 31 || | |
729 | fIcode == 33 || | |
730 | fIcode == 41 || | |
731 | fIcode == 52) return 1; | |
732 | else return 0; | |
733 | } | |
734 | ||
735 | Bool_t TFluka::IsTrackAlive() const | |
736 | { | |
737 | // means not disappeared or not out | |
738 | if (IsTrackDisappeared() || IsTrackOut() ) return 0; | |
739 | else return 1; | |
740 | } | |
741 | ||
742 | // | |
743 | // secondaries | |
744 | // | |
745 | ||
746 | Int_t TFluka::NSecondaries() const | |
747 | // Number of secondary particles generated in the current step | |
748 | // fIcode = 100 = elastic interaction | |
749 | // fIcode = 101 = inelastic interaction | |
750 | // fIcode = 102 = particle decay | |
751 | // fIcode = 103 = delta ray | |
752 | // fIcode = 104 = pair production | |
753 | // fIcode = 105 = bremsstrahlung | |
754 | { | |
755 | if (fIcode >= 100 && fIcode <= 105) | |
756 | return FINUC.np + FHEAVY.npheav; | |
757 | else | |
758 | return -1; | |
759 | } | |
760 | ||
761 | void TFluka::GetSecondary(Int_t isec, Int_t& particleId, | |
762 | TLorentzVector& position, TLorentzVector& momentum) | |
763 | // | |
764 | // fIcode = 100 = elastic interaction | |
765 | // fIcode = 101 = inelastic interaction | |
766 | // fIcode = 102 = particle decay | |
767 | // fIcode = 103 = delta ray | |
768 | // fIcode = 104 = pair production | |
769 | // fIcode = 105 = bremsstrahlung | |
770 | { | |
771 | if (fIcode >= 100 && fIcode <= 105) { | |
772 | if (isec >= 0 && isec < FINUC.np) { | |
773 | particleId = PDGFromId(FINUC.kpart[isec]); | |
774 | position.SetX(fXsco); | |
775 | position.SetY(fYsco); | |
776 | position.SetZ(fZsco); | |
777 | position.SetT(FINUC.agesec[isec]); | |
778 | momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]); | |
779 | momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]); | |
780 | momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]); | |
781 | momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]); | |
782 | } | |
783 | if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) { | |
784 | Int_t jsec = isec - FINUC.np; | |
785 | particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!! | |
786 | position.SetX(fXsco); | |
787 | position.SetY(fYsco); | |
788 | position.SetZ(fZsco); | |
789 | position.SetT(FHEAVY.agheav[jsec]); | |
790 | momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]); | |
791 | momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]); | |
792 | momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]); | |
793 | if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) | |
794 | momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]); | |
795 | else if (FHEAVY.tkheav[jsec] > 6) | |
796 | momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!! | |
797 | } | |
798 | } | |
799 | } | |
800 | ||
801 | //TMCProcess ProdProcess(Int_t isec) const | |
802 | // Name of the process that has produced the secondary particles | |
803 | // in the current step | |
804 | //{ | |
805 | // will come from FINUC when called from USDRAW | |
806 | //} | |
807 | ||
808 | //Int_t StepProcesses(TArrayI &proc) const | |
809 | // Return processes active in the current step | |
810 | //{ | |
811 | //ck = total energy of the particl ???????????????? | |
812 | //} | |
813 | ||
814 | ||
815 | // =============================================================== | |
816 | void TFluka::FutoTest() | |
817 | { | |
818 | Int_t icode, mreg, newreg, particleId; | |
819 | // Int_t medium; | |
820 | Double_t rull, xsco, ysco, zsco; | |
821 | TLorentzVector position, momentum; | |
822 | icode = GetIcode(); | |
823 | if (icode == 0) { | |
824 | cout << " icode=" << icode << endl; | |
825 | /* | |
826 | cout << "TLorentzVector positionX=" << position.X() | |
827 | << "positionY=" << position.Y() | |
828 | << "positionZ=" << position.Z() | |
829 | << "timeT=" << position.T() << endl; | |
830 | cout << "TLorentzVector momentumX=" << momentum.X() | |
831 | << "momentumY=" << momentum.Y() | |
832 | << "momentumZ=" << momentum.Z() | |
833 | << "energyE=" << momentum.E() << endl; | |
834 | cout << "TrackPid=" << TrackPid() << endl; | |
835 | */ | |
836 | } | |
837 | ||
838 | else if (icode > 0 && icode <= 5) { | |
839 | // mgdraw | |
840 | mreg = GetMreg(); | |
841 | // medium = GetMedium(); | |
842 | cout << " icode=" << icode | |
843 | << " mreg=" << mreg | |
844 | // << " medium=" << medium | |
845 | << endl; | |
846 | TrackPosition(position); | |
847 | TrackMomentum(momentum); | |
848 | cout << "TLorentzVector positionX=" << position.X() | |
849 | << "positionY=" << position.Y() | |
850 | << "positionZ=" << position.Z() | |
851 | << "timeT=" << position.T() << endl; | |
852 | cout << "TLorentzVector momentumX=" << momentum.X() | |
853 | << "momentumY=" << momentum.Y() | |
854 | << "momentumZ=" << momentum.Z() | |
855 | << "energyE=" << momentum.E() << endl; | |
856 | cout << "TrackStep=" << TrackStep() << endl; | |
857 | cout << "TrackLength=" << TrackLength() << endl; | |
858 | cout << "TrackTime=" << TrackTime() << endl; | |
859 | cout << "Edep=" << Edep() << endl; | |
860 | cout << "TrackPid=" << TrackPid() << endl; | |
861 | cout << "TrackCharge=" << TrackCharge() << endl; | |
862 | cout << "TrackMass=" << TrackMass() << endl; | |
863 | cout << "Etot=" << Etot() << endl; | |
864 | cout << "IsNewTrack=" << IsNewTrack() << endl; | |
865 | cout << "IsTrackInside=" << IsTrackInside() << endl; | |
866 | cout << "IsTrackEntering=" << IsTrackEntering() << endl; | |
867 | cout << "IsTrackExiting=" << IsTrackExiting() << endl; | |
868 | cout << "IsTrackOut=" << IsTrackOut() << endl; | |
869 | cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl; | |
870 | cout << "IsTrackAlive=" << IsTrackAlive() << endl; | |
871 | } | |
872 | ||
873 | else if((icode >= 10 && icode <= 15) || | |
874 | (icode >= 20 && icode <= 24) || | |
875 | (icode >= 30 && icode <= 33) || | |
876 | (icode >= 40 && icode <= 41) || | |
877 | (icode >= 50 && icode <= 52)) { | |
878 | // endraw | |
879 | mreg = GetMreg(); | |
880 | // medium = GetMedium(); | |
881 | rull = GetRull(); | |
882 | xsco = GetXsco(); | |
883 | ysco = GetYsco(); | |
884 | zsco = GetZsco(); | |
885 | cout << " icode=" << icode | |
886 | << " mreg=" << mreg | |
887 | // << " medium=" << medium | |
888 | << " rull=" << rull | |
889 | << " xsco=" << xsco | |
890 | << " ysco=" << ysco | |
891 | << " zsco=" << zsco << endl; | |
892 | TrackPosition(position); | |
893 | TrackMomentum(momentum); | |
894 | cout << "Edep=" << Edep() << endl; | |
895 | cout << "Etot=" << Etot() << endl; | |
896 | cout << "TrackPid=" << TrackPid() << endl; | |
897 | cout << "TrackCharge=" << TrackCharge() << endl; | |
898 | cout << "TrackMass=" << TrackMass() << endl; | |
899 | cout << "IsTrackOut=" << IsTrackOut() << endl; | |
900 | cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl; | |
901 | cout << "IsTrackStop=" << IsTrackStop() << endl; | |
902 | cout << "IsTrackAlive=" << IsTrackAlive() << endl; | |
903 | } | |
904 | ||
905 | else if((icode >= 100 && icode <= 105) || | |
906 | (icode == 208) || | |
907 | (icode == 210) || | |
908 | (icode == 212) || | |
909 | (icode >= 214 && icode <= 215) || | |
910 | (icode == 217) || | |
911 | (icode == 219) || | |
912 | (icode == 221) || | |
913 | (icode == 225) || | |
914 | (icode == 300) || | |
915 | (icode == 400)) { | |
916 | // usdraw | |
917 | mreg = GetMreg(); | |
918 | // medium = GetMedium(); | |
919 | xsco = GetXsco(); | |
920 | ysco = GetYsco(); | |
921 | zsco = GetZsco(); | |
922 | cout << " icode=" << icode | |
923 | << " mreg=" << mreg | |
924 | // << " medium=" << medium | |
925 | << " xsco=" << xsco | |
926 | << " ysco=" << ysco | |
927 | << " zsco=" << zsco << endl; | |
928 | cout << "TrackPid=" << TrackPid() << endl; | |
929 | cout << "NSecondaries=" << NSecondaries() << endl; | |
930 | for (Int_t isec=0; isec< NSecondaries(); isec++) { | |
931 | //void TFluka::GetSecondary(Int_t isec, Int_t& particleId, | |
932 | // TLorentzVector& position, TLorentzVector& momentum) | |
933 | TFluka::GetSecondary(isec, particleId, position, momentum); | |
934 | cout << "TLorentzVector positionX=" << position.X() | |
935 | << "positionY=" << position.Y() | |
936 | << "positionZ=" << position.Z() | |
937 | << "timeT=" << position.T() << endl; | |
938 | cout << "TLorentzVector momentumX=" << momentum.X() | |
939 | << "momentumY=" << momentum.Y() | |
940 | << "momentumZ=" << momentum.Z() | |
941 | << "energyE=" << momentum.E() << endl; | |
942 | cout << "TrackPid=" << particleId << endl; | |
943 | ||
944 | } | |
945 | } | |
946 | ||
947 | else if((icode == 19) || | |
948 | (icode == 29) || | |
949 | (icode == 39) || | |
950 | (icode == 49) || | |
951 | (icode == 59)) { | |
952 | mreg = GetMreg(); | |
953 | // medium = GetMedium(); | |
954 | newreg = GetNewreg(); | |
955 | xsco = GetXsco(); | |
956 | ysco = GetYsco(); | |
957 | zsco = GetZsco(); | |
958 | cout << " icode=" << icode | |
959 | << " mreg=" << mreg | |
960 | // << " medium=" << medium | |
961 | << " newreg=" << newreg | |
962 | << " xsco=" << xsco | |
963 | << " ysco=" << ysco | |
964 | << " zsco=" << zsco << endl; | |
965 | } | |
966 | // | |
967 | // ==================================================================== | |
968 | // | |
969 | ||
970 | ||
971 | ||
972 | } // end of FutoTest |