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