SetTrackIsExiting, SetTrackIsEntering, SetTrackIsInside added.
[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$
24969d13 18Revision 1.15 2003/02/18 16:12:17 morsch
19Protect mpdgha against negative argument.
20
6015a930 21Revision 1.14 2003/02/18 12:47:59 morsch
22Gmtod and Gdtom added.
23
dc37cac6 24Revision 1.13 2003/01/31 14:01:51 morsch
25Major update on
26- Getters related to geometry.
27- Communication with run manager (event steering)
28
b0d8df96 29Revision 1.12 2003/01/31 12:18:53 morsch
30Corrected indices. (E. Futo)
31
c230803a 32Revision 1.9 2002/12/06 12:41:29 morsch
33Mess from last merge cleaned up.
34
a7e55c80 35Revision 1.8 2002/12/06 12:28:44 morsch
36Region to media mapping corrected and improved.
37
6d4d27f2 38Revision 1.7 2002/12/06 12:21:32 morsch
39User stepping methods added (E. Futo)
40
fa3d1cc7 41Revision 1.6 2002/11/21 18:40:06 iglez2
42Media handling added
43
27b2f7fe 44Revision 1.5 2002/11/07 17:59:10 iglez2
45Included the geometry through geant4_vmc/FLUGG
46
bf3aa28e 47Revision 1.4 2002/11/04 16:00:46 iglez2
48The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
49
f9cb2fec 50Revision 1.3 2002/10/22 15:12:14 alibrary
51Introducing Riostream.h
52
eae0fe66 53Revision 1.2 2002/10/14 14:57:40 hristov
54Merging the VirtualMC branch to the main development branch (HEAD)
55
b9d0a01d 56Revision 1.1.2.8 2002/10/08 16:33:17 iglez2
57LSOUIT is set to true before the second call to flukam.
58
59Revision 1.1.2.7 2002/10/08 09:30:37 iglez2
60Solved stupid missing ;
61
62Revision 1.1.2.6 2002/10/07 13:40:22 iglez2
63First implementations of the PDG <--> Fluka Id conversion routines
64
65Revision 1.1.2.5 2002/09/26 16:26:03 iglez2
66Added verbosity
67Call to gAlice->Generator()->Generate()
68
69Revision 1.1.2.4 2002/09/26 13:22:23 iglez2
70Naive implementation of ProcessRun and ProcessEvent
71Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
72
73Revision 1.1.2.3 2002/09/20 15:35:51 iglez2
74Modification of LFDRTR. Value is passed to FLUKA !!!
75
76Revision 1.1.2.2 2002/09/18 14:34:44 iglez2
77Revised version with all pure virtual methods implemented
78
79Revision 1.1.2.1 2002/07/24 08:49:41 alibrary
80Adding TFluka to VirtualMC
81
82Revision 1.1 2002/07/05 13:10:07 morsch
83First commit of Fluka interface.
84
03ca248b 85*/
86
eae0fe66 87#include <Riostream.h>
b9d0a01d 88
6d4d27f2 89#include "TClonesArray.h"
03ca248b 90#include "TFluka.h"
b9d0a01d 91#include "TCallf77.h" //For the fortran calls
92#include "Fdblprc.h" //(DBLPRC) fluka common
b9d0a01d 93#include "Fepisor.h" //(EPISOR) fluka common
fa3d1cc7 94#include "Ffinuc.h" //(FINUC) fluka common
95#include "Fiounit.h" //(IOUNIT) fluka common
96#include "Fpaprop.h" //(PAPROP) fluka common
f9cb2fec 97#include "Fpart.h" //(PART) fluka common
fa3d1cc7 98#include "Ftrackr.h" //(TRACKR) fluka common
6d4d27f2 99#include "Fpaprop.h" //(PAPROP) fluka common
fa3d1cc7 100#include "Ffheavy.h" //(FHEAVY) fluka common
b9d0a01d 101
fa3d1cc7 102#include "TVirtualMC.h"
bf3aa28e 103#include "TG4GeometryManager.h" //For the geometry management
104#include "TG4DetConstruction.h" //For the detector construction
105
106#include "FGeometryInit.hh"
fa3d1cc7 107#include "TLorentzVector.h"
6d4d27f2 108#include "FlukaVolume.h"
bf3aa28e 109
b9d0a01d 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_
f9cb2fec 115# define mcihad mcihad_
116# define mpdgha mpdgha_
b9d0a01d 117#else
118# define flukam FLUKAM
119# define fluka_openinp FLUKA_OPENINP
120# define fluka_closeinp FLUKA_CLOSEINP
f9cb2fec 121# define mcihad MCIHAD
122# define mpdgha MPDGHA
b9d0a01d 123#endif
124
125extern "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&);
f9cb2fec 133 int type_of_call mcihad(const int&);
134 int type_of_call mpdgha(const int&);
b9d0a01d 135}
136
137//
138// Class implementation for ROOT
139//
03ca248b 140ClassImp(TFluka)
b9d0a01d 141
142//
bf3aa28e 143//----------------------------------------------------------------------------
144// TFluka constructors and destructors.
b9d0a01d 145//____________________________________________________________________________
146TFluka::TFluka()
147 :TVirtualMC(),
148 fVerbosityLevel(0),
bf3aa28e 149 fInputFileName(""),
27b2f7fe 150 fDetector(0),
151 fCurrentFlukaRegion(-1)
b9d0a01d 152{
153 //
154 // Default constructor
155 //
156}
157
b9d0a01d 158TFluka::TFluka(const char *title, Int_t verbosity)
159 :TVirtualMC("TFluka",title),
160 fVerbosityLevel(verbosity),
bf3aa28e 161 fInputFileName(""),
24969d13 162 fTrackIsEntering(0),
163 fTrackIsExiting(0),
27b2f7fe 164 fDetector(0),
165 fCurrentFlukaRegion(-1)
b9d0a01d 166{
167 if (fVerbosityLevel >=3)
168 cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
169
bf3aa28e 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
b9d0a01d 187 if (fVerbosityLevel >=3)
188 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
6d4d27f2 189
190 fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
191 fNVolumes = 0;
192 fMediaByRegion = 0;
b9d0a01d 193}
194
bf3aa28e 195TFluka::~TFluka() {
196 if (fVerbosityLevel >=3)
197 cout << "==> TFluka::~TFluka() destructor called." << endl;
198
199 delete fGeometryManager;
6d4d27f2 200 fVolumeMediaMap->Delete();
201 delete fVolumeMediaMap;
202
bf3aa28e 203
204 if (fVerbosityLevel >=3)
205 cout << "<== TFluka::~TFluka() destructor called." << endl;
206}
207
208//
209//_____________________________________________________________________________
210// TFluka control methods
b9d0a01d 211//____________________________________________________________________________
212void 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;
bf3aa28e 228 flukam(1);
b9d0a01d 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;
fa3d1cc7 236
6d4d27f2 237 FinishGeometry();
238
b9d0a01d 239}
240
bf3aa28e 241void TFluka::FinishGeometry() {
6d4d27f2 242//
243// Build-up table with region to medium correspondance
244//
245 char tmp[5];
246
bf3aa28e 247 if (fVerbosityLevel >=3)
248 cout << "==> TFluka::FinishGeometry() called." << endl;
249
6d4d27f2 250// fGeometryManager->Ggclos();
bf3aa28e 251
6d4d27f2 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';
b0d8df96 263 flugg->SetMediumFromName(tmp, media, i+1);
264 fMediaByRegion[i] = media;
27b2f7fe 265 }
6d4d27f2 266
267 flugg->BuildMediaMap();
27b2f7fe 268
bf3aa28e 269 if (fVerbosityLevel >=3)
270 cout << "<== TFluka::FinishGeometry() called." << endl;
271}
272
273void 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
b9d0a01d 282void TFluka::ProcessEvent() {
283 if (fVerbosityLevel >=3)
284 cout << "==> TFluka::ProcessEvent() called." << endl;
b0d8df96 285 fApplication->GeneratePrimaries();
286 EPISOR.lsouit = true;
287 flukam(1);
b9d0a01d 288 if (fVerbosityLevel >=3)
289 cout << "<== TFluka::ProcessEvent() called." << endl;
290}
291
bf3aa28e 292
b9d0a01d 293void 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 }
b0d8df96 302 fApplication->InitGeometry();
303 fApplication->BeginEvent();
304 ProcessEvent();
305 fApplication->FinishEvent();
b9d0a01d 306 if (fVerbosityLevel >=3)
307 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
308 << endl;
b0d8df96 309
b9d0a01d 310}
311
bf3aa28e 312//_____________________________________________________________________________
313// methods for building/management of geometry
314//____________________________________________________________________________
315// functions from GCONS
316void 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
323void 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
331void 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}
338void 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
346void 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}
352void 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
359void 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}
368void 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
378void 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
386void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
387//
388 fGeometryManager->Gstpar(itmed, param, parval);
389}
390
391// functions from GGEOM
392Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
393 Float_t *upar, Int_t np) {
394//
6d4d27f2 395// fVolumeMediaMap[TString(name)] = nmed;
b0d8df96 396 printf("TFluka::Gsvolu() name = %s, nmed = %d\n", name, nmed);
397
6d4d27f2 398 TClonesArray &lvols = *fVolumeMediaMap;
399 new(lvols[fNVolumes++])
400 FlukaVolume(name, nmed);
401 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
bf3aa28e 402}
403Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
404 Double_t *upar, Int_t np) {
405//
6d4d27f2 406 TClonesArray &lvols = *fVolumeMediaMap;
407 new(lvols[fNVolumes++])
408 FlukaVolume(name, nmed);
409
410 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
bf3aa28e 411}
412
413void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
414 Int_t iaxis) {
415//
b0d8df96 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);
6d4d27f2 422 fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
bf3aa28e 423}
424
425void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
426 Int_t iaxis, Double_t c0i, Int_t numed) {
427//
6d4d27f2 428 TClonesArray &lvols = *fVolumeMediaMap;
429 new(lvols[fNVolumes++])
430 FlukaVolume(name, numed);
431 fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
bf3aa28e 432}
433
434void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
435 Int_t iaxis, Int_t numed, Int_t ndvmx) {
6d4d27f2 436//
437 TClonesArray &lvols = *fVolumeMediaMap;
438 new(lvols[fNVolumes++])
439 FlukaVolume(name, numed);
440 fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
bf3aa28e 441}
442
443void 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//
6d4d27f2 446 TClonesArray &lvols = *fVolumeMediaMap;
447 new(lvols[fNVolumes++])
448 FlukaVolume(name, numed);
449 fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
bf3aa28e 450}
451
452void TFluka::Gsord(const char *name, Int_t iax) {
453//
454 fGeometryManager->Gsord(name, iax);
455}
456
457void 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
464void 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}
470void 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
477void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
478//
479 fGeometryManager->Gsbool(onlyVolName, manyVolName);
480}
481
482void 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}
487void 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
494void 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
27b2f7fe 502//_____________________________________________________________________________
503// methods needed by the stepping
504//____________________________________________________________________________
6d4d27f2 505
27b2f7fe 506Int_t TFluka::GetMedium() const {
b0d8df96 507//
508// Get the medium number for the current fluka region
509//
6d4d27f2 510 FGeometryInit* flugg = FGeometryInit::GetInstance();
511 return flugg->GetMedium(fCurrentFlukaRegion);
27b2f7fe 512}
bf3aa28e 513
514
515
516//____________________________________________________________________________
517// ID <--> PDG transformations
b9d0a01d 518//_____________________________________________________________________________
519Int_t TFluka::IdFromPDG(Int_t pdg) const
520{
521 //
f9cb2fec 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);
b9d0a01d 528}
529
b9d0a01d 530Int_t TFluka::PDGFromId(Int_t id) const
531{
532 //
f9cb2fec 533 // Return PDG code and pseudo ENDF code from Fluka code
c230803a 534
bc021b12 535 //IPTOKP array goes from official to internal
b0d8df96 536 if (id == 0) {
537 printf("PDGFromId: Error id = 0");
538 return -1;
539 }
540
bc021b12 541 Int_t intfluka = GetFlukaIPTOKP(id);
b0d8df96 542 if (intfluka == 0) {
543 printf("PDGFromId: Error intfluka = 0");
544 return -1;
6015a930 545 } else if (intfluka < 0) {
546 printf("PDGFromId: Error intfluka < 0");
547 return -1;
b0d8df96 548 }
6015a930 549 printf("mpdgha called with %d %d \n", id, intfluka);
bc021b12 550 return mpdgha(intfluka);
6d4d27f2 551}
552
fa3d1cc7 553//_____________________________________________________________________________
554// methods for step management
555//____________________________________________________________________________
bc021b12 556//
557// set methods
558//
559void TFluka::SetMaxStep(Double_t)
560{
561// SetMaxStep is dummy procedure in TFluka !
562 cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
563}
564
565void TFluka::SetMaxNStep(Int_t)
566{
567// SetMaxNStep is dummy procedure in TFluka !
568 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
569}
570
571void TFluka::SetUserDecay(Int_t)
572{
573// SetUserDecay is dummy procedure in TFluka !
574 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
575}
576
fa3d1cc7 577//
578// dynamic properties
579//
580void 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);
24969d13 592//
593//
594//
595
fa3d1cc7 596}
597
598void 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 {
b8b430a9 616 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
fa3d1cc7 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
625Double_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
632Double_t TFluka::TrackLength() const
633{
b8b430a9 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
fa3d1cc7 639// for the time being returns only the length in centimeters of the current step
b8b430a9 640 Double_t sum = 0;
641 for ( Int_t j=0;j<TRACKR.ntrack;j++) {
642 sum +=TRACKR.ttrack[j];
643 }
644 return sum;
fa3d1cc7 645}
646
647Double_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
654Double_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
676Int_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
683Double_t TFluka::TrackCharge() const
684{
685// Return charge of the track currently transported
686// PAPROP.ichrge = electric charge of the particle
bc021b12 687// TRACKR.jtrack = identity number of the particle
fa3d1cc7 688 return PAPROP.ichrge[TRACKR.jtrack+6];
689}
690
691Double_t TFluka::TrackMass() const
692{
693// PAPROP.am = particle mass in GeV
bc021b12 694// TRACKR.jtrack = identity number of the particle
fa3d1cc7 695 return PAPROP.am[TRACKR.jtrack+6];
696}
697
698Double_t TFluka::Etot() const
699{
700// TRACKR.etrack = total energy of the particle
701 return TRACKR.etrack;
702}
703
704//
705// track status
706//
707Bool_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
715Bool_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
726Bool_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 ||
24969d13 739 fIcode == 59 ||
740 fTrackIsEntering) return 1;
fa3d1cc7 741 else return 0;
742}
743
744Bool_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 ||
24969d13 757 fIcode == 59 ||
758 fTrackIsExiting) return 1;
fa3d1cc7 759 else return 0;
760}
761
762Bool_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
779Bool_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
792Bool_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
817Bool_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
828Int_t TFluka::NSecondaries() const
829// Number of secondary particles generated in the current step
bc021b12 830// FINUC.np = number of secondaries except light and heavy ions
b8b430a9 831// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
fa3d1cc7 832{
bc021b12 833 return FINUC.np + FHEAVY.npheav;
fa3d1cc7 834}
835
836void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
837 TLorentzVector& position, TLorentzVector& momentum)
fa3d1cc7 838{
bc021b12 839 if (isec >= 0 && isec < FINUC.np) {
b8b430a9 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
bc021b12 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]);
fa3d1cc7 868 }
bc021b12 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 }
fa3d1cc7 885}
886
bc021b12 887TMCProcess TFluka::ProdProcess(Int_t isec) const
fa3d1cc7 888// Name of the process that has produced the secondary particles
889// in the current step
bc021b12 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;
b0d8df96 911// const TMCProcess kIpPCerenkov = kPCerenkov;
912// const TMCProcess kIpPSynchrotron = kPSynchrotron;
bc021b12 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}
fa3d1cc7 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
b0d8df96 949Int_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
958const 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
969Int_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
992Int_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
1004Int_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
1023const 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
1035const 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
1047Int_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
dc37cac6 1060void 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
1086void 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
1111void 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 }
1132void 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 }
b0d8df96 1153
fa3d1cc7 1154// ===============================================================
1155void 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;
24969d13 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]);
fa3d1cc7 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++) {
fa3d1cc7 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
b0d8df96 1322