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