]>
Commit | Line | Data |
---|---|---|
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$ | |
18 | Revision 1.15 2003/02/18 16:12:17 morsch | |
19 | Protect mpdgha against negative argument. | |
20 | ||
21 | Revision 1.14 2003/02/18 12:47:59 morsch | |
22 | Gmtod and Gdtom added. | |
23 | ||
24 | Revision 1.13 2003/01/31 14:01:51 morsch | |
25 | Major update on | |
26 | - Getters related to geometry. | |
27 | - Communication with run manager (event steering) | |
28 | ||
29 | Revision 1.12 2003/01/31 12:18:53 morsch | |
30 | Corrected indices. (E. Futo) | |
31 | ||
32 | Revision 1.9 2002/12/06 12:41:29 morsch | |
33 | Mess from last merge cleaned up. | |
34 | ||
35 | Revision 1.8 2002/12/06 12:28:44 morsch | |
36 | Region to media mapping corrected and improved. | |
37 | ||
38 | Revision 1.7 2002/12/06 12:21:32 morsch | |
39 | User stepping methods added (E. Futo) | |
40 | ||
41 | Revision 1.6 2002/11/21 18:40:06 iglez2 | |
42 | Media handling added | |
43 | ||
44 | Revision 1.5 2002/11/07 17:59:10 iglez2 | |
45 | Included the geometry through geant4_vmc/FLUGG | |
46 | ||
47 | Revision 1.4 2002/11/04 16:00:46 iglez2 | |
48 | The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent. | |
49 | ||
50 | Revision 1.3 2002/10/22 15:12:14 alibrary | |
51 | Introducing Riostream.h | |
52 | ||
53 | Revision 1.2 2002/10/14 14:57:40 hristov | |
54 | Merging the VirtualMC branch to the main development branch (HEAD) | |
55 | ||
56 | Revision 1.1.2.8 2002/10/08 16:33:17 iglez2 | |
57 | LSOUIT is set to true before the second call to flukam. | |
58 | ||
59 | Revision 1.1.2.7 2002/10/08 09:30:37 iglez2 | |
60 | Solved stupid missing ; | |
61 | ||
62 | Revision 1.1.2.6 2002/10/07 13:40:22 iglez2 | |
63 | First implementations of the PDG <--> Fluka Id conversion routines | |
64 | ||
65 | Revision 1.1.2.5 2002/09/26 16:26:03 iglez2 | |
66 | Added verbosity | |
67 | Call to gAlice->Generator()->Generate() | |
68 | ||
69 | Revision 1.1.2.4 2002/09/26 13:22:23 iglez2 | |
70 | Naive implementation of ProcessRun and ProcessEvent | |
71 | Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init() | |
72 | ||
73 | Revision 1.1.2.3 2002/09/20 15:35:51 iglez2 | |
74 | Modification of LFDRTR. Value is passed to FLUKA !!! | |
75 | ||
76 | Revision 1.1.2.2 2002/09/18 14:34:44 iglez2 | |
77 | Revised version with all pure virtual methods implemented | |
78 | ||
79 | Revision 1.1.2.1 2002/07/24 08:49:41 alibrary | |
80 | Adding TFluka to VirtualMC | |
81 | ||
82 | Revision 1.1 2002/07/05 13:10:07 morsch | |
83 | First commit of Fluka interface. | |
84 | ||
85 | */ | |
86 | ||
87 | #include <Riostream.h> | |
88 | ||
89 | #include "TClonesArray.h" | |
90 | #include "TFluka.h" | |
91 | #include "TCallf77.h" //For the fortran calls | |
92 | #include "Fdblprc.h" //(DBLPRC) fluka common | |
93 | #include "Fepisor.h" //(EPISOR) fluka common | |
94 | #include "Ffinuc.h" //(FINUC) fluka common | |
95 | #include "Fiounit.h" //(IOUNIT) fluka common | |
96 | #include "Fpaprop.h" //(PAPROP) fluka common | |
97 | #include "Fpart.h" //(PART) fluka common | |
98 | #include "Ftrackr.h" //(TRACKR) fluka common | |
99 | #include "Fpaprop.h" //(PAPROP) fluka common | |
100 | #include "Ffheavy.h" //(FHEAVY) fluka common | |
101 | ||
102 | #include "TVirtualMC.h" | |
103 | #include "TG4GeometryManager.h" //For the geometry management | |
104 | #include "TG4DetConstruction.h" //For the detector construction | |
105 | ||
106 | #include "FGeometryInit.hh" | |
107 | #include "TLorentzVector.h" | |
108 | #include "FlukaVolume.h" | |
109 | ||
110 | // Fluka methods that may be needed. | |
111 | #ifndef WIN32 | |
112 | # define flukam flukam_ | |
113 | # define fluka_openinp fluka_openinp_ | |
114 | # define fluka_closeinp fluka_closeinp_ | |
115 | # define mcihad mcihad_ | |
116 | # define mpdgha mpdgha_ | |
117 | #else | |
118 | # define flukam FLUKAM | |
119 | # define fluka_openinp FLUKA_OPENINP | |
120 | # define fluka_closeinp FLUKA_CLOSEINP | |
121 | # define mcihad MCIHAD | |
122 | # define mpdgha MPDGHA | |
123 | #endif | |
124 | ||
125 | extern "C" | |
126 | { | |
127 | // | |
128 | // Prototypes for FLUKA functions | |
129 | // | |
130 | void type_of_call flukam(const int&); | |
131 | void type_of_call fluka_openinp(const int&, DEFCHARA); | |
132 | void type_of_call fluka_closeinp(const int&); | |
133 | int type_of_call mcihad(const int&); | |
134 | int type_of_call mpdgha(const int&); | |
135 | } | |
136 | ||
137 | // | |
138 | // Class implementation for ROOT | |
139 | // | |
140 | ClassImp(TFluka) | |
141 | ||
142 | // | |
143 | //---------------------------------------------------------------------------- | |
144 | // TFluka constructors and destructors. | |
145 | //____________________________________________________________________________ | |
146 | TFluka::TFluka() | |
147 | :TVirtualMC(), | |
148 | fVerbosityLevel(0), | |
149 | fInputFileName(""), | |
150 | fDetector(0), | |
151 | fCurrentFlukaRegion(-1) | |
152 | { | |
153 | // | |
154 | // Default constructor | |
155 | // | |
156 | } | |
157 | ||
158 | TFluka::TFluka(const char *title, Int_t verbosity) | |
159 | :TVirtualMC("TFluka",title), | |
160 | fVerbosityLevel(verbosity), | |
161 | fInputFileName(""), | |
162 | fTrackIsEntering(0), | |
163 | fTrackIsExiting(0), | |
164 | fDetector(0), | |
165 | fCurrentFlukaRegion(-1) | |
166 | { | |
167 | if (fVerbosityLevel >=3) | |
168 | cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl; | |
169 | ||
170 | ||
171 | // create geometry manager | |
172 | if (fVerbosityLevel >=2) | |
173 | cout << "\t* Creating G4 Geometry manager..." << endl; | |
174 | fGeometryManager = new TG4GeometryManager(); | |
175 | if (fVerbosityLevel >=2) | |
176 | cout << "\t* Creating G4 Detector..." << endl; | |
177 | fDetector = new TG4DetConstruction(); | |
178 | FGeometryInit* geominit = FGeometryInit::GetInstance(); | |
179 | if (geominit) | |
180 | geominit->setDetConstruction(fDetector); | |
181 | else { | |
182 | cerr << "ERROR: Could not create FGeometryInit!" << endl; | |
183 | cerr << " Exiting!!!" << endl; | |
184 | abort(); | |
185 | } | |
186 | ||
187 | if (fVerbosityLevel >=3) | |
188 | cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl; | |
189 | ||
190 | fVolumeMediaMap = new TClonesArray("FlukaVolume",1000); | |
191 | fNVolumes = 0; | |
192 | fMediaByRegion = 0; | |
193 | } | |
194 | ||
195 | TFluka::~TFluka() { | |
196 | if (fVerbosityLevel >=3) | |
197 | cout << "==> TFluka::~TFluka() destructor called." << endl; | |
198 | ||
199 | delete fGeometryManager; | |
200 | fVolumeMediaMap->Delete(); | |
201 | delete fVolumeMediaMap; | |
202 | ||
203 | ||
204 | if (fVerbosityLevel >=3) | |
205 | cout << "<== TFluka::~TFluka() destructor called." << endl; | |
206 | } | |
207 | ||
208 | // | |
209 | //_____________________________________________________________________________ | |
210 | // TFluka control methods | |
211 | //____________________________________________________________________________ | |
212 | void TFluka::Init() { | |
213 | if (fVerbosityLevel >=3) | |
214 | cout << "==> TFluka::Init() called." << endl; | |
215 | ||
216 | if (fVerbosityLevel >=2) | |
217 | cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F') | |
218 | << ") in fluka..." << endl; | |
219 | GLOBAL.lfdrtr = true; | |
220 | ||
221 | if (fVerbosityLevel >=2) | |
222 | cout << "\t* Opening file " << fInputFileName << endl; | |
223 | const char* fname = fInputFileName; | |
224 | fluka_openinp(lunin, PASSCHARA(fname)); | |
225 | ||
226 | if (fVerbosityLevel >=2) | |
227 | cout << "\t* Calling flukam..." << endl; | |
228 | flukam(1); | |
229 | ||
230 | if (fVerbosityLevel >=2) | |
231 | cout << "\t* Closing file " << fInputFileName << endl; | |
232 | fluka_closeinp(lunin); | |
233 | ||
234 | if (fVerbosityLevel >=3) | |
235 | cout << "<== TFluka::Init() called." << endl; | |
236 | ||
237 | FinishGeometry(); | |
238 | ||
239 | } | |
240 | ||
241 | void TFluka::FinishGeometry() { | |
242 | // | |
243 | // Build-up table with region to medium correspondance | |
244 | // | |
245 | char tmp[5]; | |
246 | ||
247 | if (fVerbosityLevel >=3) | |
248 | cout << "==> TFluka::FinishGeometry() called." << endl; | |
249 | ||
250 | // fGeometryManager->Ggclos(); | |
251 | ||
252 | FGeometryInit* flugg = FGeometryInit::GetInstance(); | |
253 | ||
254 | fMediaByRegion = new Int_t[fNVolumes+2]; | |
255 | for (Int_t i = 0; i < fNVolumes; i++) | |
256 | { | |
257 | FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]); | |
258 | TString volName = vol->GetName(); | |
259 | Int_t media = vol->GetMedium(); | |
260 | printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media); | |
261 | strcpy(tmp, volName.Data()); | |
262 | tmp[4] = '\0'; | |
263 | flugg->SetMediumFromName(tmp, media, i+1); | |
264 | fMediaByRegion[i] = media; | |
265 | } | |
266 | ||
267 | flugg->BuildMediaMap(); | |
268 | ||
269 | if (fVerbosityLevel >=3) | |
270 | cout << "<== TFluka::FinishGeometry() called." << endl; | |
271 | } | |
272 | ||
273 | void TFluka::BuildPhysics() { | |
274 | if (fVerbosityLevel >=3) | |
275 | cout << "==> TFluka::BuildPhysics() called." << endl; | |
276 | ||
277 | ||
278 | if (fVerbosityLevel >=3) | |
279 | cout << "<== TFluka::BuildPhysics() called." << endl; | |
280 | } | |
281 | ||
282 | void TFluka::ProcessEvent() { | |
283 | if (fVerbosityLevel >=3) | |
284 | cout << "==> TFluka::ProcessEvent() called." << endl; | |
285 | fApplication->GeneratePrimaries(); | |
286 | EPISOR.lsouit = true; | |
287 | flukam(1); | |
288 | if (fVerbosityLevel >=3) | |
289 | cout << "<== TFluka::ProcessEvent() called." << endl; | |
290 | } | |
291 | ||
292 | ||
293 | void TFluka::ProcessRun(Int_t nevent) { | |
294 | if (fVerbosityLevel >=3) | |
295 | cout << "==> TFluka::ProcessRun(" << nevent << ") called." | |
296 | << endl; | |
297 | ||
298 | if (fVerbosityLevel >=2) { | |
299 | cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl; | |
300 | cout << "\t* Calling flukam again..." << endl; | |
301 | } | |
302 | fApplication->InitGeometry(); | |
303 | fApplication->BeginEvent(); | |
304 | ProcessEvent(); | |
305 | fApplication->FinishEvent(); | |
306 | if (fVerbosityLevel >=3) | |
307 | cout << "<== TFluka::ProcessRun(" << nevent << ") called." | |
308 | << endl; | |
309 | ||
310 | } | |
311 | ||
312 | //_____________________________________________________________________________ | |
313 | // methods for building/management of geometry | |
314 | //____________________________________________________________________________ | |
315 | // functions from GCONS | |
316 | void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z, | |
317 | Float_t &dens, Float_t &radl, Float_t &absl, | |
318 | Float_t* ubuf, Int_t& nbuf) { | |
319 | // | |
320 | fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf); | |
321 | } | |
322 | ||
323 | void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z, | |
324 | Double_t &dens, Double_t &radl, Double_t &absl, | |
325 | Double_t* ubuf, Int_t& nbuf) { | |
326 | // | |
327 | fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf); | |
328 | } | |
329 | ||
330 | // detector composition | |
331 | void TFluka::Material(Int_t& kmat, const char* name, Double_t a, | |
332 | Double_t z, Double_t dens, Double_t radl, Double_t absl, | |
333 | Float_t* buf, Int_t nwbuf) { | |
334 | // | |
335 | fGeometryManager | |
336 | ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); | |
337 | } | |
338 | void TFluka::Material(Int_t& kmat, const char* name, Double_t a, | |
339 | Double_t z, Double_t dens, Double_t radl, Double_t absl, | |
340 | Double_t* buf, Int_t nwbuf) { | |
341 | // | |
342 | fGeometryManager | |
343 | ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); | |
344 | } | |
345 | ||
346 | void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, | |
347 | Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) { | |
348 | // | |
349 | fGeometryManager | |
350 | ->Mixture(kmat, name, a, z, dens, nlmat, wmat); | |
351 | } | |
352 | void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, | |
353 | Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) { | |
354 | // | |
355 | fGeometryManager | |
356 | ->Mixture(kmat, name, a, z, dens, nlmat, wmat); | |
357 | } | |
358 | ||
359 | void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, | |
360 | Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, | |
361 | Double_t stemax, Double_t deemax, Double_t epsil, | |
362 | Double_t stmin, Float_t* ubuf, Int_t nbuf) { | |
363 | // | |
364 | fGeometryManager | |
365 | ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, | |
366 | epsil, stmin, ubuf, nbuf); | |
367 | } | |
368 | void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, | |
369 | Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, | |
370 | Double_t stemax, Double_t deemax, Double_t epsil, | |
371 | Double_t stmin, Double_t* ubuf, Int_t nbuf) { | |
372 | // | |
373 | fGeometryManager | |
374 | ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, | |
375 | epsil, stmin, ubuf, nbuf); | |
376 | } | |
377 | ||
378 | void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, | |
379 | Double_t thetaY, Double_t phiY, Double_t thetaZ, | |
380 | Double_t phiZ) { | |
381 | // | |
382 | fGeometryManager | |
383 | ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); | |
384 | } | |
385 | ||
386 | void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) { | |
387 | // | |
388 | fGeometryManager->Gstpar(itmed, param, parval); | |
389 | } | |
390 | ||
391 | // functions from GGEOM | |
392 | Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed, | |
393 | Float_t *upar, Int_t np) { | |
394 | // | |
395 | // fVolumeMediaMap[TString(name)] = nmed; | |
396 | printf("TFluka::Gsvolu() name = %s, nmed = %d\n", name, nmed); | |
397 | ||
398 | TClonesArray &lvols = *fVolumeMediaMap; | |
399 | new(lvols[fNVolumes++]) | |
400 | FlukaVolume(name, nmed); | |
401 | return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); | |
402 | } | |
403 | Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed, | |
404 | Double_t *upar, Int_t np) { | |
405 | // | |
406 | TClonesArray &lvols = *fVolumeMediaMap; | |
407 | new(lvols[fNVolumes++]) | |
408 | FlukaVolume(name, nmed); | |
409 | ||
410 | return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); | |
411 | } | |
412 | ||
413 | void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, | |
414 | Int_t iaxis) { | |
415 | // | |
416 | // The medium of the daughter is the one of the mother | |
417 | Int_t volid = TFluka::VolId(mother); | |
418 | Int_t med = TFluka::VolId2Mate(volid); | |
419 | TClonesArray &lvols = *fVolumeMediaMap; | |
420 | new(lvols[fNVolumes++]) | |
421 | FlukaVolume(name, med); | |
422 | fGeometryManager->Gsdvn(name, mother, ndiv, iaxis); | |
423 | } | |
424 | ||
425 | void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, | |
426 | Int_t iaxis, Double_t c0i, Int_t numed) { | |
427 | // | |
428 | TClonesArray &lvols = *fVolumeMediaMap; | |
429 | new(lvols[fNVolumes++]) | |
430 | FlukaVolume(name, numed); | |
431 | fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); | |
432 | } | |
433 | ||
434 | void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, | |
435 | Int_t iaxis, Int_t numed, Int_t ndvmx) { | |
436 | // | |
437 | TClonesArray &lvols = *fVolumeMediaMap; | |
438 | new(lvols[fNVolumes++]) | |
439 | FlukaVolume(name, numed); | |
440 | fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx); | |
441 | } | |
442 | ||
443 | void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, | |
444 | Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { | |
445 | // | |
446 | TClonesArray &lvols = *fVolumeMediaMap; | |
447 | new(lvols[fNVolumes++]) | |
448 | FlukaVolume(name, numed); | |
449 | fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); | |
450 | } | |
451 | ||
452 | void TFluka::Gsord(const char *name, Int_t iax) { | |
453 | // | |
454 | fGeometryManager->Gsord(name, iax); | |
455 | } | |
456 | ||
457 | void TFluka::Gspos(const char *name, Int_t nr, const char *mother, | |
458 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
459 | const char *konly) { | |
460 | // | |
461 | fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly); | |
462 | } | |
463 | ||
464 | void TFluka::Gsposp(const char *name, Int_t nr, const char *mother, | |
465 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
466 | const char *konly, Float_t *upar, Int_t np) { | |
467 | // | |
468 | fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); | |
469 | } | |
470 | void TFluka::Gsposp(const char *name, Int_t nr, const char *mother, | |
471 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
472 | const char *konly, Double_t *upar, Int_t np) { | |
473 | // | |
474 | fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); | |
475 | } | |
476 | ||
477 | void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) { | |
478 | // | |
479 | fGeometryManager->Gsbool(onlyVolName, manyVolName); | |
480 | } | |
481 | ||
482 | void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov, | |
483 | Float_t *absco, Float_t *effic, Float_t *rindex) { | |
484 | // | |
485 | fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex); | |
486 | } | |
487 | void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov, | |
488 | Double_t *absco, Double_t *effic, Double_t *rindex) { | |
489 | // | |
490 | fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex); | |
491 | } | |
492 | ||
493 | // Euclid | |
494 | void TFluka::WriteEuclid(const char* fileName, const char* topVol, | |
495 | Int_t number, Int_t nlevel) { | |
496 | // | |
497 | fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel); | |
498 | } | |
499 | ||
500 | ||
501 | ||
502 | //_____________________________________________________________________________ | |
503 | // methods needed by the stepping | |
504 | //____________________________________________________________________________ | |
505 | ||
506 | Int_t TFluka::GetMedium() const { | |
507 | // | |
508 | // Get the medium number for the current fluka region | |
509 | // | |
510 | FGeometryInit* flugg = FGeometryInit::GetInstance(); | |
511 | return flugg->GetMedium(fCurrentFlukaRegion); | |
512 | } | |
513 | ||
514 | ||
515 | ||
516 | //____________________________________________________________________________ | |
517 | // ID <--> PDG transformations | |
518 | //_____________________________________________________________________________ | |
519 | Int_t TFluka::IdFromPDG(Int_t pdg) const | |
520 | { | |
521 | // | |
522 | // Return Fluka code from PDG and pseudo ENDF code | |
523 | ||
524 | // MCIHAD() goes from pdg to fluka internal. | |
525 | Int_t intfluka = mcihad(pdg); | |
526 | // KPTOIP array goes from internal to official | |
527 | return GetFlukaKPTOIP(intfluka); | |
528 | } | |
529 | ||
530 | Int_t TFluka::PDGFromId(Int_t id) const | |
531 | { | |
532 | // | |
533 | // Return PDG code and pseudo ENDF code from Fluka code | |
534 | ||
535 | //IPTOKP array goes from official to internal | |
536 | if (id == 0) { | |
537 | printf("PDGFromId: Error id = 0"); | |
538 | return -1; | |
539 | } | |
540 | ||
541 | Int_t intfluka = GetFlukaIPTOKP(id); | |
542 | if (intfluka == 0) { | |
543 | printf("PDGFromId: Error intfluka = 0"); | |
544 | return -1; | |
545 | } else if (intfluka < 0) { | |
546 | printf("PDGFromId: Error intfluka < 0"); | |
547 | return -1; | |
548 | } | |
549 | printf("mpdgha called with %d %d \n", id, intfluka); | |
550 | return mpdgha(intfluka); | |
551 | } | |
552 | ||
553 | //_____________________________________________________________________________ | |
554 | // methods for step management | |
555 | //____________________________________________________________________________ | |
556 | // | |
557 | // set methods | |
558 | // | |
559 | void TFluka::SetMaxStep(Double_t) | |
560 | { | |
561 | // SetMaxStep is dummy procedure in TFluka ! | |
562 | cout << "SetMaxStep is dummy procedure in TFluka !" << endl; | |
563 | } | |
564 | ||
565 | void TFluka::SetMaxNStep(Int_t) | |
566 | { | |
567 | // SetMaxNStep is dummy procedure in TFluka ! | |
568 | cout << "SetMaxNStep is dummy procedure in TFluka !" << endl; | |
569 | } | |
570 | ||
571 | void TFluka::SetUserDecay(Int_t) | |
572 | { | |
573 | // SetUserDecay is dummy procedure in TFluka ! | |
574 | cout << "SetUserDecay is dummy procedure in TFluka !" << endl; | |
575 | } | |
576 | ||
577 | // | |
578 | // dynamic properties | |
579 | // | |
580 | void TFluka::TrackPosition(TLorentzVector& position) const | |
581 | { | |
582 | // Return the current position in the master reference frame of the | |
583 | // track being transported | |
584 | // TRACKR.atrack = age of the particle | |
585 | // TRACKR.xtrack = x-position of the last point | |
586 | // TRACKR.ytrack = y-position of the last point | |
587 | // TRACKR.ztrack = z-position of the last point | |
588 | position.SetX(TRACKR.xtrack[TRACKR.ntrack]); | |
589 | position.SetY(TRACKR.ytrack[TRACKR.ntrack]); | |
590 | position.SetZ(TRACKR.ztrack[TRACKR.ntrack]); | |
591 | position.SetT(TRACKR.atrack); | |
592 | // | |
593 | // | |
594 | // | |
595 | ||
596 | } | |
597 | ||
598 | void TFluka::TrackMomentum(TLorentzVector& momentum) const | |
599 | { | |
600 | // Return the direction and the momentum (GeV/c) of the track | |
601 | // currently being transported | |
602 | // TRACKR.ptrack = momentum of the particle (not always defined, if | |
603 | // < 0 must be obtained from etrack) | |
604 | // TRACKR.cx,y,ztrck = direction cosines of the current particle | |
605 | // TRACKR.etrack = total energy of the particle | |
606 | // TRACKR.jtrack = identity number of the particle | |
607 | // PAPROP.am[TRACKR.jtrack] = particle mass in gev | |
608 | if (TRACKR.ptrack >= 0) { | |
609 | momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck); | |
610 | momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck); | |
611 | momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck); | |
612 | momentum.SetE(TRACKR.etrack); | |
613 | return; | |
614 | } | |
615 | else { | |
616 | Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]); | |
617 | momentum.SetPx(p*TRACKR.cxtrck); | |
618 | momentum.SetPy(p*TRACKR.cytrck); | |
619 | momentum.SetPz(p*TRACKR.cztrck); | |
620 | momentum.SetE(TRACKR.etrack); | |
621 | return; | |
622 | } | |
623 | } | |
624 | ||
625 | Double_t TFluka::TrackStep() const | |
626 | { | |
627 | // Return the length in centimeters of the current step | |
628 | // TRACKR.ctrack = total curved path | |
629 | return TRACKR.ctrack; | |
630 | } | |
631 | ||
632 | Double_t TFluka::TrackLength() const | |
633 | { | |
634 | // Still wrong !!! | |
635 | // This is the sum of substeps !!! | |
636 | // TRACKR.ctrack = total curved path of the current step | |
637 | // Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field | |
638 | // The sum of all step length starting from the beginning of the track | |
639 | // for the time being returns only the length in centimeters of the current step | |
640 | Double_t sum = 0; | |
641 | for ( Int_t j=0;j<TRACKR.ntrack;j++) { | |
642 | sum +=TRACKR.ttrack[j]; | |
643 | } | |
644 | return sum; | |
645 | } | |
646 | ||
647 | Double_t TFluka::TrackTime() const | |
648 | { | |
649 | // Return the current time of flight of the track being transported | |
650 | // TRACKR.atrack = age of the particle | |
651 | return TRACKR.atrack; | |
652 | } | |
653 | ||
654 | Double_t TFluka::Edep() const | |
655 | { | |
656 | // Energy deposition | |
657 | // if TRACKR.ntrack = 0, TRACKR.mtrack = 0: | |
658 | // -->local energy deposition (the value and the point are not recorded in TRACKR) | |
659 | // but in the variable "rull" of the procedure "endraw.cxx" | |
660 | // if TRACKR.ntrack > 0, TRACKR.mtrack = 0: | |
661 | // -->no energy loss along the track | |
662 | // if TRACKR.ntrack > 0, TRACKR.mtrack > 0: | |
663 | // -->energy loss distributed along the track | |
664 | // TRACKR.dtrack = energy deposition of the jth deposition even | |
665 | if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0) | |
666 | return fRull; | |
667 | else { | |
668 | Double_t sum = 0; | |
669 | for ( Int_t j=0;j<TRACKR.mtrack;j++) { | |
670 | sum +=TRACKR.dtrack[j]; | |
671 | } | |
672 | return sum; | |
673 | } | |
674 | } | |
675 | ||
676 | Int_t TFluka::TrackPid() const | |
677 | { | |
678 | // Return the id of the particle transported | |
679 | // TRACKR.jtrack = identity number of the particle | |
680 | return PDGFromId(TRACKR.jtrack); | |
681 | } | |
682 | ||
683 | Double_t TFluka::TrackCharge() const | |
684 | { | |
685 | // Return charge of the track currently transported | |
686 | // PAPROP.ichrge = electric charge of the particle | |
687 | // TRACKR.jtrack = identity number of the particle | |
688 | return PAPROP.ichrge[TRACKR.jtrack+6]; | |
689 | } | |
690 | ||
691 | Double_t TFluka::TrackMass() const | |
692 | { | |
693 | // PAPROP.am = particle mass in GeV | |
694 | // TRACKR.jtrack = identity number of the particle | |
695 | return PAPROP.am[TRACKR.jtrack+6]; | |
696 | } | |
697 | ||
698 | Double_t TFluka::Etot() const | |
699 | { | |
700 | // TRACKR.etrack = total energy of the particle | |
701 | return TRACKR.etrack; | |
702 | } | |
703 | ||
704 | // | |
705 | // track status | |
706 | // | |
707 | Bool_t TFluka::IsNewTrack() const | |
708 | { | |
709 | // ???????????????, | |
710 | // True if the track is not at the boundary of the current volume | |
711 | // Not true in some cases in bxdraw - to be solved | |
712 | return 1; | |
713 | } | |
714 | ||
715 | Bool_t TFluka::IsTrackInside() const | |
716 | { | |
717 | // True if the track is not at the boundary of the current volume | |
718 | // In Fluka a step is always inside one kind of material | |
719 | // If the step would go behind the region of one material, | |
720 | // it will be shortened to reach only the boundary. | |
721 | // Therefore IsTrackInside() is always true. | |
722 | // Not true in some cases in bxdraw - to be solved | |
723 | return 1; | |
724 | } | |
725 | ||
726 | Bool_t TFluka::IsTrackEntering() const | |
727 | { | |
728 | // True if this is the first step of the track in the current volume | |
729 | // Boundary- (X) crossing | |
730 | // Icode = 19: boundary crossing - call from Kaskad | |
731 | // Icode = 29: boundary crossing - call from Emfsco | |
732 | // Icode = 39: boundary crossing - call from Kasneu | |
733 | // Icode = 49: boundary crossing - call from Kashea | |
734 | // Icode = 59: boundary crossing - call from Kasoph | |
735 | if (fIcode == 19 || | |
736 | fIcode == 29 || | |
737 | fIcode == 39 || | |
738 | fIcode == 49 || | |
739 | fIcode == 59 || | |
740 | fTrackIsEntering) return 1; | |
741 | else return 0; | |
742 | } | |
743 | ||
744 | Bool_t TFluka::IsTrackExiting() const | |
745 | { | |
746 | // True if this is the last step of the track in the current volume | |
747 | // Boundary- (X) crossing | |
748 | // Icode = 19: boundary crossing - call from Kaskad | |
749 | // Icode = 29: boundary crossing - call from Emfsco | |
750 | // Icode = 39: boundary crossing - call from Kasneu | |
751 | // Icode = 49: boundary crossing - call from Kashea | |
752 | // Icode = 59: boundary crossing - call from Kasoph | |
753 | if (fIcode == 19 || | |
754 | fIcode == 29 || | |
755 | fIcode == 39 || | |
756 | fIcode == 49 || | |
757 | fIcode == 59 || | |
758 | fTrackIsExiting) return 1; | |
759 | else return 0; | |
760 | } | |
761 | ||
762 | Bool_t TFluka::IsTrackOut() const | |
763 | { | |
764 | // True if the track is out of the setup | |
765 | // means escape | |
766 | // Icode = 14: escape - call from Kaskad | |
767 | // Icode = 23: escape - call from Emfsco | |
768 | // Icode = 32: escape - call from Kasneu | |
769 | // Icode = 40: escape - call from Kashea | |
770 | // Icode = 51: escape - call from Kasoph | |
771 | if (fIcode == 14 || | |
772 | fIcode == 23 || | |
773 | fIcode == 32 || | |
774 | fIcode == 40 || | |
775 | fIcode == 51) return 1; | |
776 | else return 0; | |
777 | } | |
778 | ||
779 | Bool_t TFluka::IsTrackDisappeared() const | |
780 | { | |
781 | // means all inelastic interactions and decays | |
782 | // fIcode from usdraw | |
783 | if (fIcode == 101 || // inelastic interaction | |
784 | fIcode == 102 || // particle decay | |
785 | fIcode == 214 || // in-flight annihilation | |
786 | fIcode == 215 || // annihilation at rest | |
787 | fIcode == 217 || // pair production | |
788 | fIcode == 221) return 1; | |
789 | else return 0; | |
790 | } | |
791 | ||
792 | Bool_t TFluka::IsTrackStop() const | |
793 | { | |
794 | // True if the track energy has fallen below the threshold | |
795 | // means stopped by signal or below energy threshold | |
796 | // Icode = 12: stopping particle - call from Kaskad | |
797 | // Icode = 15: time kill - call from Kaskad | |
798 | // Icode = 21: below threshold, iarg=1 - call from Emfsco | |
799 | // Icode = 22: below threshold, iarg=2 - call from Emfsco | |
800 | // Icode = 24: time kill - call from Emfsco | |
801 | // Icode = 31: below threshold - call from Kasneu | |
802 | // Icode = 33: time kill - call from Kasneu | |
803 | // Icode = 41: time kill - call from Kashea | |
804 | // Icode = 52: time kill - call from Kasoph | |
805 | if (fIcode == 12 || | |
806 | fIcode == 15 || | |
807 | fIcode == 21 || | |
808 | fIcode == 22 || | |
809 | fIcode == 24 || | |
810 | fIcode == 31 || | |
811 | fIcode == 33 || | |
812 | fIcode == 41 || | |
813 | fIcode == 52) return 1; | |
814 | else return 0; | |
815 | } | |
816 | ||
817 | Bool_t TFluka::IsTrackAlive() const | |
818 | { | |
819 | // means not disappeared or not out | |
820 | if (IsTrackDisappeared() || IsTrackOut() ) return 0; | |
821 | else return 1; | |
822 | } | |
823 | ||
824 | // | |
825 | // secondaries | |
826 | // | |
827 | ||
828 | Int_t TFluka::NSecondaries() const | |
829 | // Number of secondary particles generated in the current step | |
830 | // FINUC.np = number of secondaries except light and heavy ions | |
831 | // FHEAVY.npheav = number of secondaries for light and heavy secondary ions | |
832 | { | |
833 | return FINUC.np + FHEAVY.npheav; | |
834 | } | |
835 | ||
836 | void TFluka::GetSecondary(Int_t isec, Int_t& particleId, | |
837 | TLorentzVector& position, TLorentzVector& momentum) | |
838 | { | |
839 | if (isec >= 0 && isec < FINUC.np) { | |
840 | // more fine condition depending on icode | |
841 | // icode = 100 ? | |
842 | // icode = 101 OK | |
843 | // icode = 102 OK | |
844 | // icode = 103 ? | |
845 | // icode = 104 ? | |
846 | // icode = 105 ? | |
847 | // icode = 208 ? | |
848 | // icode = 210 ? | |
849 | // icode = 212 ? | |
850 | // icode = 214 OK | |
851 | // icode = 215 OK | |
852 | // icode = 219 ? | |
853 | // icode = 221 OK | |
854 | // icode = 225 ? | |
855 | // icode = 300 ? | |
856 | // icode = 400 ? | |
857 | ||
858 | particleId = PDGFromId(FINUC.kpart[isec]); | |
859 | position.SetX(fXsco); | |
860 | position.SetY(fYsco); | |
861 | position.SetZ(fZsco); | |
862 | position.SetT(TRACKR.atrack); | |
863 | // position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem. | |
864 | momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]); | |
865 | momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]); | |
866 | momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]); | |
867 | momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]); | |
868 | } | |
869 | if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) { | |
870 | Int_t jsec = isec - FINUC.np; | |
871 | particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!! | |
872 | position.SetX(fXsco); | |
873 | position.SetY(fYsco); | |
874 | position.SetZ(fZsco); | |
875 | position.SetT(TRACKR.atrack); | |
876 | // position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem. | |
877 | momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]); | |
878 | momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]); | |
879 | momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]); | |
880 | if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) | |
881 | momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]); | |
882 | else if (FHEAVY.tkheav[jsec] > 6) | |
883 | momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!! | |
884 | } | |
885 | } | |
886 | ||
887 | TMCProcess TFluka::ProdProcess(Int_t isec) const | |
888 | // Name of the process that has produced the secondary particles | |
889 | // in the current step | |
890 | { | |
891 | const TMCProcess kIpNoProc = kPNoProcess; | |
892 | const TMCProcess kIpPDecay = kPDecay; | |
893 | const TMCProcess kIpPPair = kPPair; | |
894 | //const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton; | |
895 | //const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton; | |
896 | const TMCProcess kIpPCompton = kPCompton; | |
897 | const TMCProcess kIpPPhotoelectric = kPPhotoelectric; | |
898 | const TMCProcess kIpPBrem = kPBrem; | |
899 | //const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy; | |
900 | //const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron; | |
901 | const TMCProcess kIpPDeltaRay = kPDeltaRay; | |
902 | //const TMCProcess kIpPMoller = kPMoller; | |
903 | //const TMCProcess kIpPBhabha = kPBhabha; | |
904 | const TMCProcess kIpPAnnihilation = kPAnnihilation; | |
905 | //const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight; | |
906 | //const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest; | |
907 | const TMCProcess kIpPHadronic = kPHadronic; | |
908 | const TMCProcess kIpPMuonNuclear = kPMuonNuclear; | |
909 | const TMCProcess kIpPPhotoFission = kPPhotoFission; | |
910 | const TMCProcess kIpPRayleigh = kPRayleigh; | |
911 | // const TMCProcess kIpPCerenkov = kPCerenkov; | |
912 | // const TMCProcess kIpPSynchrotron = kPSynchrotron; | |
913 | ||
914 | Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11; | |
915 | if (fIcode == 102) return kIpPDecay; | |
916 | else if (fIcode == 104 || fIcode == 217) return kIpPPair; | |
917 | //else if (fIcode == 104) return kIpPairFromPhoton; | |
918 | //else if (fIcode == 217) return kIpPPairFromVirtualPhoton; | |
919 | else if (fIcode == 219) return kIpPCompton; | |
920 | else if (fIcode == 221) return kIpPPhotoelectric; | |
921 | else if (fIcode == 105 || fIcode == 208) return kIpPBrem; | |
922 | //else if (fIcode == 105) return kIpPBremFromHeavy; | |
923 | //else if (fIcode == 208) return kPBremFromElectronOrPositron; | |
924 | else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay; | |
925 | else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay; | |
926 | //else if (fIcode == 210) return kIpPMoller; | |
927 | //else if (fIcode == 212) return kIpPBhabha; | |
928 | else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation; | |
929 | //else if (fIcode == 214) return kIpPAnnihilInFlight; | |
930 | //else if (fIcode == 215) return kIpPAnnihilAtRest; | |
931 | else if (fIcode == 101) return kIpPHadronic; | |
932 | else if (fIcode == 101) { | |
933 | if (!mugamma) return kIpPHadronic; | |
934 | else if (TRACKR.jtrack == 7) return kIpPPhotoFission; | |
935 | else return kIpPMuonNuclear; | |
936 | } | |
937 | else if (fIcode == 225) return kIpPRayleigh; | |
938 | // Fluka codes 100, 300 and 400 still to be investigasted | |
939 | else return kIpNoProc; | |
940 | } | |
941 | ||
942 | //Int_t StepProcesses(TArrayI &proc) const | |
943 | // Return processes active in the current step | |
944 | //{ | |
945 | //ck = total energy of the particl ???????????????? | |
946 | //} | |
947 | ||
948 | ||
949 | Int_t TFluka::VolId2Mate(Int_t id) const | |
950 | { | |
951 | // | |
952 | // Returns the material number for a given volume ID | |
953 | // | |
954 | printf("VolId2Mate %d %d\n", id, fMediaByRegion[id]); | |
955 | return fMediaByRegion[id-1]; | |
956 | } | |
957 | ||
958 | const char* TFluka::VolName(Int_t id) const | |
959 | { | |
960 | // | |
961 | // Returns the volume name for a given volume ID | |
962 | // | |
963 | FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[id-1]); | |
964 | const char* name = vol->GetName(); | |
965 | printf("VolName %d %s \n", id, name); | |
966 | return name; | |
967 | } | |
968 | ||
969 | Int_t TFluka::VolId(const Text_t* volName) const | |
970 | { | |
971 | // | |
972 | // Converts from volume name to volume ID. | |
973 | // Time consuming. (Only used during set-up) | |
974 | // Could be replaced by hash-table | |
975 | // | |
976 | char tmp[5]; | |
977 | Int_t i =0; | |
978 | for (i = 0; i < fNVolumes; i++) | |
979 | { | |
980 | FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]); | |
981 | TString name = vol->GetName(); | |
982 | strcpy(tmp, name.Data()); | |
983 | tmp[4] = '\0'; | |
984 | if (!strcmp(tmp, volName)) break; | |
985 | } | |
986 | i++; | |
987 | ||
988 | return i; | |
989 | } | |
990 | ||
991 | ||
992 | Int_t TFluka::CurrentVolID(Int_t& copyNo) const | |
993 | { | |
994 | // | |
995 | // Return the logical id and copy number corresponding to the current fluka region | |
996 | // | |
997 | int ir = fCurrentFlukaRegion; | |
998 | int id = (FGeometryInit::GetInstance())->CurrentVolID(ir, copyNo); | |
999 | printf("CurrentVolID: %d %d %d \n", ir, id, copyNo); | |
1000 | return id; | |
1001 | ||
1002 | } | |
1003 | ||
1004 | Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const | |
1005 | { | |
1006 | // | |
1007 | // Return the logical id and copy number of off'th mother | |
1008 | // corresponding to the current fluka region | |
1009 | // | |
1010 | if (off == 0) | |
1011 | return CurrentVolID(copyNo); | |
1012 | ||
1013 | int ir = fCurrentFlukaRegion; | |
1014 | int id = (FGeometryInit::GetInstance())->CurrentVolOffID(ir, off, copyNo); | |
1015 | ||
1016 | printf("CurrentVolOffID: %d %d %d \n", ir, id, copyNo); | |
1017 | if (id == -1) | |
1018 | printf("CurrentVolOffID: Warning Mother not found !!!\n"); | |
1019 | return id; | |
1020 | } | |
1021 | ||
1022 | ||
1023 | const char* TFluka::CurrentVolName() const | |
1024 | { | |
1025 | // | |
1026 | // Return the current volume name | |
1027 | // | |
1028 | Int_t copy; | |
1029 | Int_t id = TFluka::CurrentVolID(copy); | |
1030 | const char* name = TFluka::VolName(id); | |
1031 | printf("CurrentVolumeName: %d %s \n", fCurrentFlukaRegion, name); | |
1032 | return name; | |
1033 | } | |
1034 | ||
1035 | const char* TFluka::CurrentVolOffName(Int_t off) const | |
1036 | { | |
1037 | // | |
1038 | // Return the volume name of the off'th mother of the current volume | |
1039 | // | |
1040 | Int_t copy; | |
1041 | Int_t id = TFluka::CurrentVolOffID(off, copy); | |
1042 | const char* name = TFluka::VolName(id); | |
1043 | printf("CurrentVolumeOffName: %d %s \n", fCurrentFlukaRegion, name); | |
1044 | return name; | |
1045 | } | |
1046 | ||
1047 | Int_t TFluka::CurrentMaterial(Float_t &a, Float_t &z, | |
1048 | Float_t &dens, Float_t &radl, Float_t &absl) const | |
1049 | { | |
1050 | // | |
1051 | // Return the current medium number | |
1052 | // | |
1053 | Int_t copy; | |
1054 | Int_t id = TFluka::CurrentVolID(copy); | |
1055 | Int_t med = TFluka::VolId2Mate(id); | |
1056 | printf("CurrentMaterial: %d %d \n", fCurrentFlukaRegion, med); | |
1057 | return med; | |
1058 | } | |
1059 | ||
1060 | void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag) | |
1061 | { | |
1062 | // Transforms a position from the world reference frame | |
1063 | // to the current volume reference frame. | |
1064 | // | |
1065 | // Geant3 desription: | |
1066 | // ================== | |
1067 | // Computes coordinates XD (in DRS) | |
1068 | // from known coordinates XM in MRS | |
1069 | // The local reference system can be initialized by | |
1070 | // - the tracking routines and GMTOD used in GUSTEP | |
1071 | // - a call to GMEDIA(XM,NUMED) | |
1072 | // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) | |
1073 | // (inverse routine is GDTOM) | |
1074 | // | |
1075 | // If IFLAG=1 convert coordinates | |
1076 | // IFLAG=2 convert direction cosinus | |
1077 | // | |
1078 | // --- | |
1079 | Double_t xmD[3], xdD[3]; | |
1080 | xmD[0] = xm[0]; xmD[1] = xm[1]; xmD[2] = xm[2]; | |
1081 | (FGeometryInit::GetInstance())->Gmtod(xmD, xdD, iflag); | |
1082 | xd[0] = xdD[0]; xd[1] = xdD[1]; xd[2] = xdD[2]; | |
1083 | } | |
1084 | ||
1085 | ||
1086 | void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag) | |
1087 | { | |
1088 | // Transforms a position from the world reference frame | |
1089 | // to the current volume reference frame. | |
1090 | // | |
1091 | // Geant3 desription: | |
1092 | // ================== | |
1093 | // Computes coordinates XD (in DRS) | |
1094 | // from known coordinates XM in MRS | |
1095 | // The local reference system can be initialized by | |
1096 | // - the tracking routines and GMTOD used in GUSTEP | |
1097 | // - a call to GMEDIA(XM,NUMED) | |
1098 | // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) | |
1099 | // (inverse routine is GDTOM) | |
1100 | // | |
1101 | // If IFLAG=1 convert coordinates | |
1102 | // IFLAG=2 convert direction cosinus | |
1103 | // | |
1104 | // --- | |
1105 | Double_t xmD[3], xdD[3]; | |
1106 | xdD[0] = xd[0]; xdD[1] = xd[1]; xdD[2] = xd[2]; | |
1107 | (FGeometryInit::GetInstance())->Gdtom(xmD, xdD, iflag); | |
1108 | xm[0] = xmD[0]; xm[1] = xmD[1]; xm[2] = xmD[2]; | |
1109 | } | |
1110 | ||
1111 | void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag) | |
1112 | { | |
1113 | // Transforms a position from the current volume reference frame | |
1114 | // to the world reference frame. | |
1115 | // | |
1116 | // Geant3 desription: | |
1117 | // ================== | |
1118 | // Computes coordinates XM (Master Reference System | |
1119 | // knowing the coordinates XD (Detector Ref System) | |
1120 | // The local reference system can be initialized by | |
1121 | // - the tracking routines and GDTOM used in GUSTEP | |
1122 | // - a call to GSCMED(NLEVEL,NAMES,NUMBER) | |
1123 | // (inverse routine is GMTOD) | |
1124 | // | |
1125 | // If IFLAG=1 convert coordinates | |
1126 | // IFLAG=2 convert direction cosinus | |
1127 | // | |
1128 | // --- | |
1129 | ||
1130 | ||
1131 | } | |
1132 | void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag) | |
1133 | { | |
1134 | // Transforms a position from the current volume reference frame | |
1135 | // to the world reference frame. | |
1136 | // | |
1137 | // Geant3 desription: | |
1138 | // ================== | |
1139 | // Computes coordinates XM (Master Reference System | |
1140 | // knowing the coordinates XD (Detector Ref System) | |
1141 | // The local reference system can be initialized by | |
1142 | // - the tracking routines and GDTOM used in GUSTEP | |
1143 | // - a call to GSCMED(NLEVEL,NAMES,NUMBER) | |
1144 | // (inverse routine is GMTOD) | |
1145 | // | |
1146 | // If IFLAG=1 convert coordinates | |
1147 | // IFLAG=2 convert direction cosinus | |
1148 | // | |
1149 | // --- | |
1150 | ||
1151 | (FGeometryInit::GetInstance())->Gdtom(xm, xd, iflag); | |
1152 | } | |
1153 | ||
1154 | // =============================================================== | |
1155 | void TFluka::FutoTest() | |
1156 | { | |
1157 | Int_t icode, mreg, newreg, particleId; | |
1158 | // Int_t medium; | |
1159 | Double_t rull, xsco, ysco, zsco; | |
1160 | TLorentzVector position, momentum; | |
1161 | icode = GetIcode(); | |
1162 | if (icode == 0) { | |
1163 | cout << " icode=" << icode << endl; | |
1164 | /* | |
1165 | cout << "TLorentzVector positionX=" << position.X() | |
1166 | << "positionY=" << position.Y() | |
1167 | << "positionZ=" << position.Z() | |
1168 | << "timeT=" << position.T() << endl; | |
1169 | cout << "TLorentzVector momentumX=" << momentum.X() | |
1170 | << "momentumY=" << momentum.Y() | |
1171 | << "momentumZ=" << momentum.Z() | |
1172 | << "energyE=" << momentum.E() << endl; | |
1173 | cout << "TrackPid=" << TrackPid() << endl; | |
1174 | */ | |
1175 | } | |
1176 | ||
1177 | else if (icode > 0 && icode <= 5) { | |
1178 | // mgdraw | |
1179 | mreg = GetMreg(); | |
1180 | // medium = GetMedium(); | |
1181 | cout << " icode=" << icode | |
1182 | << " mreg=" << mreg | |
1183 | // << " medium=" << medium | |
1184 | << endl; | |
1185 | TrackPosition(position); | |
1186 | TrackMomentum(momentum); | |
1187 | cout << "TLorentzVector positionX=" << position.X() | |
1188 | << "positionY=" << position.Y() | |
1189 | << "positionZ=" << position.Z() | |
1190 | << "timeT=" << position.T() << endl; | |
1191 | cout << "TLorentzVector momentumX=" << momentum.X() | |
1192 | << "momentumY=" << momentum.Y() | |
1193 | << "momentumZ=" << momentum.Z() | |
1194 | << "energyE=" << momentum.E() << endl; | |
1195 | cout << "TrackStep=" << TrackStep() << endl; | |
1196 | cout << "TrackLength=" << TrackLength() << endl; | |
1197 | cout << "TrackTime=" << TrackTime() << endl; | |
1198 | cout << "Edep=" << Edep() << endl; | |
1199 | cout << "TrackPid=" << TrackPid() << endl; | |
1200 | cout << "TrackCharge=" << TrackCharge() << endl; | |
1201 | cout << "TrackMass=" << TrackMass() << endl; | |
1202 | cout << "Etot=" << Etot() << endl; | |
1203 | cout << "IsNewTrack=" << IsNewTrack() << endl; | |
1204 | cout << "IsTrackInside=" << IsTrackInside() << endl; | |
1205 | cout << "IsTrackEntering=" << IsTrackEntering() << endl; | |
1206 | cout << "IsTrackExiting=" << IsTrackExiting() << endl; | |
1207 | cout << "IsTrackOut=" << IsTrackOut() << endl; | |
1208 | cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl; | |
1209 | cout << "IsTrackAlive=" << IsTrackAlive() << endl; | |
1210 | ||
1211 | Float_t x = position.X(); | |
1212 | Float_t y = position.Y(); | |
1213 | Float_t z = position.Z(); | |
1214 | Float_t xm[3]; | |
1215 | Float_t xd[3]; | |
1216 | xm[0] = x; xm[1] = y; xm[2] = z; | |
1217 | printf("Global trackPosition: %f %f %f \n", x, y, z); | |
1218 | Gmtod(xm, xd, 1); | |
1219 | printf("Local trackPosition: %f %f %f \n", xd[0], xd[1], xd[2]); | |
1220 | Gdtom(xd, xm, 1); | |
1221 | printf("New trackPosition: %f %f %f \n", xm[0], xm[1], xm[2]); | |
1222 | } | |
1223 | ||
1224 | else if((icode >= 10 && icode <= 15) || | |
1225 | (icode >= 20 && icode <= 24) || | |
1226 | (icode >= 30 && icode <= 33) || | |
1227 | (icode >= 40 && icode <= 41) || | |
1228 | (icode >= 50 && icode <= 52)) { | |
1229 | // endraw | |
1230 | mreg = GetMreg(); | |
1231 | // medium = GetMedium(); | |
1232 | rull = GetRull(); | |
1233 | xsco = GetXsco(); | |
1234 | ysco = GetYsco(); | |
1235 | zsco = GetZsco(); | |
1236 | cout << " icode=" << icode | |
1237 | << " mreg=" << mreg | |
1238 | // << " medium=" << medium | |
1239 | << " rull=" << rull | |
1240 | << " xsco=" << xsco | |
1241 | << " ysco=" << ysco | |
1242 | << " zsco=" << zsco << endl; | |
1243 | TrackPosition(position); | |
1244 | TrackMomentum(momentum); | |
1245 | cout << "Edep=" << Edep() << endl; | |
1246 | cout << "Etot=" << Etot() << endl; | |
1247 | cout << "TrackPid=" << TrackPid() << endl; | |
1248 | cout << "TrackCharge=" << TrackCharge() << endl; | |
1249 | cout << "TrackMass=" << TrackMass() << endl; | |
1250 | cout << "IsTrackOut=" << IsTrackOut() << endl; | |
1251 | cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl; | |
1252 | cout << "IsTrackStop=" << IsTrackStop() << endl; | |
1253 | cout << "IsTrackAlive=" << IsTrackAlive() << endl; | |
1254 | } | |
1255 | ||
1256 | else if((icode >= 100 && icode <= 105) || | |
1257 | (icode == 208) || | |
1258 | (icode == 210) || | |
1259 | (icode == 212) || | |
1260 | (icode >= 214 && icode <= 215) || | |
1261 | (icode == 217) || | |
1262 | (icode == 219) || | |
1263 | (icode == 221) || | |
1264 | (icode == 225) || | |
1265 | (icode == 300) || | |
1266 | (icode == 400)) { | |
1267 | // usdraw | |
1268 | mreg = GetMreg(); | |
1269 | // medium = GetMedium(); | |
1270 | xsco = GetXsco(); | |
1271 | ysco = GetYsco(); | |
1272 | zsco = GetZsco(); | |
1273 | cout << " icode=" << icode | |
1274 | << " mreg=" << mreg | |
1275 | // << " medium=" << medium | |
1276 | << " xsco=" << xsco | |
1277 | << " ysco=" << ysco | |
1278 | << " zsco=" << zsco << endl; | |
1279 | cout << "TrackPid=" << TrackPid() << endl; | |
1280 | cout << "NSecondaries=" << NSecondaries() << endl; | |
1281 | for (Int_t isec=0; isec< NSecondaries(); isec++) { | |
1282 | TFluka::GetSecondary(isec, particleId, position, momentum); | |
1283 | cout << "TLorentzVector positionX=" << position.X() | |
1284 | << "positionY=" << position.Y() | |
1285 | << "positionZ=" << position.Z() | |
1286 | << "timeT=" << position.T() << endl; | |
1287 | cout << "TLorentzVector momentumX=" << momentum.X() | |
1288 | << "momentumY=" << momentum.Y() | |
1289 | << "momentumZ=" << momentum.Z() | |
1290 | << "energyE=" << momentum.E() << endl; | |
1291 | cout << "TrackPid=" << particleId << endl; | |
1292 | ||
1293 | } | |
1294 | } | |
1295 | ||
1296 | else if((icode == 19) || | |
1297 | (icode == 29) || | |
1298 | (icode == 39) || | |
1299 | (icode == 49) || | |
1300 | (icode == 59)) { | |
1301 | mreg = GetMreg(); | |
1302 | // medium = GetMedium(); | |
1303 | newreg = GetNewreg(); | |
1304 | xsco = GetXsco(); | |
1305 | ysco = GetYsco(); | |
1306 | zsco = GetZsco(); | |
1307 | cout << " icode=" << icode | |
1308 | << " mreg=" << mreg | |
1309 | // << " medium=" << medium | |
1310 | << " newreg=" << newreg | |
1311 | << " xsco=" << xsco | |
1312 | << " ysco=" << ysco | |
1313 | << " zsco=" << zsco << endl; | |
1314 | } | |
1315 | // | |
1316 | // ==================================================================== | |
1317 | // | |
1318 | ||
1319 | ||
1320 | ||
1321 | } // end of FutoTest | |
1322 |