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