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