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