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