PDGFromId: protection against intfluka = 0.
[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   // IPTOKP array goes from official to internal
504
505     Int_t intfluka = GetFlukaIPTOKP(id);
506
507     if (! intfluka) {
508         printf("\n Warning PDGFromId: internal id of %d is zero \n", id);
509         return -1;
510     }
511     
512     //MPKDHA() goes from internal to PDG
513     return mpdgha(intfluka);
514 }
515
516 //_____________________________________________________________________________
517 // methods for step management
518 //____________________________________________________________________________ 
519 //
520 // dynamic properties
521 //
522 void TFluka::TrackPosition(TLorentzVector& position) const
523 {
524 // Return the current position in the master reference frame of the
525 // track being transported
526 // TRACKR.atrack = age of the particle
527 // TRACKR.xtrack = x-position of the last point
528 // TRACKR.ytrack = y-position of the last point
529 // TRACKR.ztrack = z-position of the last point
530   position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
531   position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
532   position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
533   position.SetT(TRACKR.atrack);
534 }
535
536 void TFluka::TrackMomentum(TLorentzVector& momentum) const
537 {
538 // Return the direction and the momentum (GeV/c) of the track
539 // currently being transported
540 // TRACKR.ptrack = momentum of the particle (not always defined, if
541 //               < 0 must be obtained from etrack) 
542 // TRACKR.cx,y,ztrck = direction cosines of the current particle
543 // TRACKR.etrack = total energy of the particle
544 // TRACKR.jtrack = identity number of the particle
545 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
546   if (TRACKR.ptrack >= 0) {
547     momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
548     momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
549     momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
550     momentum.SetE(TRACKR.etrack);
551     return;
552   }
553   else {
554     Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
555     momentum.SetPx(p*TRACKR.cxtrck);
556     momentum.SetPy(p*TRACKR.cytrck);
557     momentum.SetPz(p*TRACKR.cztrck);
558     momentum.SetE(TRACKR.etrack);
559     return;
560   }
561 }
562
563 Double_t TFluka::TrackStep() const
564 {
565 // Return the length in centimeters of the current step
566 // TRACKR.ctrack = total curved path
567     return TRACKR.ctrack;
568 }
569
570 Double_t TFluka::TrackLength() const
571 {
572 // It is wrong
573 // should be the sum of all steps starting from the beginning of the track
574 // for the time being returns only the length in centimeters of the current step
575     return TRACKR.ctrack;
576 }
577
578 Double_t TFluka::TrackTime() const
579 {
580 // Return the current time of flight of the track being transported
581 // TRACKR.atrack = age of the particle
582   return TRACKR.atrack;
583 }
584
585 Double_t TFluka::Edep() const
586 {
587 // Energy deposition
588 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
589 // -->local energy deposition (the value and the point are not recorded in TRACKR)
590 //    but in the variable "rull" of the procedure "endraw.cxx"
591 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
592 // -->no energy loss along the track
593 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
594 // -->energy loss distributed along the track
595 // TRACKR.dtrack = energy deposition of the jth deposition even
596   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
597     return fRull;
598   else {
599     Double_t sum = 0;
600     for ( Int_t j=0;j<TRACKR.mtrack;j++) {
601       sum +=TRACKR.dtrack[j];  
602     }
603     return sum;
604   }
605 }
606
607 Int_t TFluka::TrackPid() const
608 {
609 // Return the id of the particle transported
610 // TRACKR.jtrack = identity number of the particle
611   return PDGFromId(TRACKR.jtrack);
612 }
613
614 Double_t TFluka::TrackCharge() const
615 {
616 // Return charge of the track currently transported
617 // PAPROP.ichrge = electric charge of the particle
618   return PAPROP.ichrge[TRACKR.jtrack+6];
619 }
620
621 Double_t TFluka::TrackMass() const
622 {
623 // PAPROP.am = particle mass in GeV
624   return PAPROP.am[TRACKR.jtrack+6];
625 }
626
627 Double_t TFluka::Etot() const
628 {
629 // TRACKR.etrack = total energy of the particle
630   return TRACKR.etrack;
631 }
632
633 //
634 // track status
635 //
636 Bool_t   TFluka::IsNewTrack() const
637 {
638 // ???????????????,
639 // True if the track is not at the boundary of the current volume
640 // Not true in some cases in bxdraw - to be solved
641   return 1;
642 }
643
644 Bool_t   TFluka::IsTrackInside() const
645 {
646 // True if the track is not at the boundary of the current volume
647 // In Fluka a step is always inside one kind of material
648 // If the step would go behind the region of one material,
649 // it will be shortened to reach only the boundary.
650 // Therefore IsTrackInside() is always true.
651 // Not true in some cases in bxdraw - to be solved
652   return 1;
653 }
654
655 Bool_t   TFluka::IsTrackEntering() const
656 {
657 // True if this is the first step of the track in the current volume
658 // Boundary- (X) crossing
659 // Icode = 19: boundary crossing - call from Kaskad
660 // Icode = 29: boundary crossing - call from Emfsco
661 // Icode = 39: boundary crossing - call from Kasneu
662 // Icode = 49: boundary crossing - call from Kashea
663 // Icode = 59: boundary crossing - call from Kasoph
664   if (fIcode == 19 ||
665       fIcode == 29 ||
666       fIcode == 39 ||
667       fIcode == 49 ||
668       fIcode == 59) return 1;
669   else return 0;
670 }
671
672 Bool_t   TFluka::IsTrackExiting() const
673 {
674 // True if this is the last step of the track in the current volume
675 // Boundary- (X) crossing
676 // Icode = 19: boundary crossing - call from Kaskad
677 // Icode = 29: boundary crossing - call from Emfsco
678 // Icode = 39: boundary crossing - call from Kasneu
679 // Icode = 49: boundary crossing - call from Kashea
680 // Icode = 59: boundary crossing - call from Kasoph
681   if (fIcode == 19 ||
682       fIcode == 29 ||
683       fIcode == 39 ||
684       fIcode == 49 ||
685       fIcode == 59) return 1;
686   else return 0;
687 }
688
689 Bool_t   TFluka::IsTrackOut() const
690 {
691 // True if the track is out of the setup
692 // means escape
693 // Icode = 14: escape - call from Kaskad
694 // Icode = 23: escape - call from Emfsco
695 // Icode = 32: escape - call from Kasneu
696 // Icode = 40: escape - call from Kashea
697 // Icode = 51: escape - call from Kasoph
698   if (fIcode == 14 ||
699       fIcode == 23 ||
700       fIcode == 32 ||
701       fIcode == 40 ||
702       fIcode == 51) return 1;
703   else return 0;
704 }
705
706 Bool_t   TFluka::IsTrackDisappeared() const
707 {
708 // means all inelastic interactions and decays
709 // fIcode from usdraw
710   if (fIcode == 101 || // inelastic interaction
711       fIcode == 102 || // particle decay
712       fIcode == 214 || // in-flight annihilation
713       fIcode == 215 || // annihilation at rest
714       fIcode == 217 || // pair production
715       fIcode == 221) return 1;
716   else return 0;
717 }
718
719 Bool_t   TFluka::IsTrackStop() const
720 {
721 // True if the track energy has fallen below the threshold
722 // means stopped by signal or below energy threshold
723 // Icode = 12: stopping particle       - call from Kaskad
724 // Icode = 15: time kill               - call from Kaskad
725 // Icode = 21: below threshold, iarg=1 - call from Emfsco
726 // Icode = 22: below threshold, iarg=2 - call from Emfsco
727 // Icode = 24: time kill               - call from Emfsco
728 // Icode = 31: below threshold         - call from Kasneu
729 // Icode = 33: time kill               - call from Kasneu
730 // Icode = 41: time kill               - call from Kashea
731 // Icode = 52: time kill               - call from Kasoph
732   if (fIcode == 12 ||
733       fIcode == 15 ||
734       fIcode == 21 ||
735       fIcode == 22 ||
736       fIcode == 24 ||
737       fIcode == 31 ||
738       fIcode == 33 ||
739       fIcode == 41 ||
740       fIcode == 52) return 1;
741   else return 0;
742 }
743
744 Bool_t   TFluka::IsTrackAlive() const
745 {
746 // means not disappeared or not out
747   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
748   else return 1;
749 }
750
751 //
752 // secondaries
753 //
754
755 Int_t TFluka::NSecondaries() const
756 // Number of secondary particles generated in the current step
757 // fIcode = 100 = elastic interaction
758 // fIcode = 101 = inelastic interaction
759 // fIcode = 102 = particle decay
760 // fIcode = 103 = delta ray
761 // fIcode = 104 = pair production
762 // fIcode = 105 = bremsstrahlung
763 {
764   if (fIcode >= 100 && fIcode <= 105)
765     return FINUC.np + FHEAVY.npheav;
766   else 
767     return -1;
768 }
769
770 void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
771                 TLorentzVector& position, TLorentzVector& momentum)
772 // 
773 // fIcode = 100 = elastic interaction
774 // fIcode = 101 = inelastic interaction
775 // fIcode = 102 = particle decay
776 // fIcode = 103 = delta ray
777 // fIcode = 104 = pair production
778 // fIcode = 105 = bremsstrahlung
779 {
780   if (fIcode >= 100 && fIcode <= 105) {
781     if (isec >= 0 && isec < FINUC.np) {
782       particleId = PDGFromId(FINUC.kpart[isec]);
783       position.SetX(fXsco);
784       position.SetY(fYsco);
785       position.SetZ(fZsco);
786       position.SetT(FINUC.agesec[isec]);
787       momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
788       momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
789       momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
790       momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
791     }
792     if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
793       Int_t jsec = isec - FINUC.np;
794       particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
795       position.SetX(fXsco);
796       position.SetY(fYsco);
797       position.SetZ(fZsco);
798       position.SetT(FHEAVY.agheav[jsec]);
799       momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
800       momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
801       momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
802       if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
803         momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
804       else if (FHEAVY.tkheav[jsec] > 6)
805         momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
806     }
807   }
808 }
809
810 //TMCProcess ProdProcess(Int_t isec) const
811 // Name of the process that has produced the secondary particles
812 // in the current step
813 //{
814 // will come from FINUC when called from USDRAW
815 //}
816
817 //Int_t StepProcesses(TArrayI &proc) const
818 // Return processes active in the current step
819 //{
820 //ck = total energy of the particl ???????????????? 
821 //}
822
823
824 // ===============================================================
825 void TFluka::FutoTest() 
826 {
827   Int_t icode, mreg, newreg, particleId;
828 //  Int_t medium;
829   Double_t rull, xsco, ysco, zsco;
830   TLorentzVector position, momentum;
831   icode = GetIcode();
832   if (icode == 0) {
833     cout << " icode=" << icode << endl;
834     /*
835     cout << "TLorentzVector positionX=" << position.X()
836        << "positionY=" << position.Y()
837        << "positionZ=" << position.Z()
838        << "timeT=" << position.T() << endl;
839     cout << "TLorentzVector momentumX=" << momentum.X()
840        << "momentumY=" << momentum.Y()
841        << "momentumZ=" << momentum.Z()
842        << "energyE=" << momentum.E() << endl;
843     cout << "TrackPid=" << TrackPid() << endl;
844     */
845   }
846
847   else if (icode > 0 && icode <= 5) {
848 // mgdraw
849     mreg = GetMreg();
850 //    medium = GetMedium();
851     cout << " icode=" << icode
852          << " mreg=" << mreg
853 //       << " medium=" << medium
854          << endl;
855   TrackPosition(position);
856   TrackMomentum(momentum);
857   cout << "TLorentzVector positionX=" << position.X()
858        << "positionY=" << position.Y()
859        << "positionZ=" << position.Z()
860        << "timeT=" << position.T() << endl;
861   cout << "TLorentzVector momentumX=" << momentum.X()
862        << "momentumY=" << momentum.Y()
863        << "momentumZ=" << momentum.Z()
864        << "energyE=" << momentum.E() << endl;
865   cout << "TrackStep=" << TrackStep() << endl;
866   cout << "TrackLength=" << TrackLength() << endl;
867   cout << "TrackTime=" << TrackTime() << endl;
868   cout << "Edep=" << Edep() << endl;
869   cout << "TrackPid=" << TrackPid() << endl;
870   cout << "TrackCharge=" << TrackCharge() << endl;
871   cout << "TrackMass=" << TrackMass() << endl;
872   cout << "Etot=" << Etot() << endl;
873   cout << "IsNewTrack=" << IsNewTrack() << endl;
874   cout << "IsTrackInside=" << IsTrackInside() << endl;
875   cout << "IsTrackEntering=" << IsTrackEntering() << endl;
876   cout << "IsTrackExiting=" << IsTrackExiting() << endl;
877   cout << "IsTrackOut=" << IsTrackOut() << endl;
878   cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
879   cout << "IsTrackAlive=" << IsTrackAlive() << endl;
880   }
881
882   else if((icode >= 10 && icode <= 15) ||
883           (icode >= 20 && icode <= 24) ||
884           (icode >= 30 && icode <= 33) ||
885           (icode >= 40 && icode <= 41) ||
886           (icode >= 50 && icode <= 52)) {
887 // endraw
888     mreg = GetMreg();
889 //    medium = GetMedium();
890     rull = GetRull();
891     xsco = GetXsco();
892     ysco = GetYsco();
893     zsco = GetZsco();
894     cout << " icode=" << icode
895          << " mreg=" << mreg
896 //       << " medium=" << medium
897          << " rull=" << rull
898          << " xsco=" << xsco
899          << " ysco=" << ysco
900          << " zsco=" << zsco << endl;
901   TrackPosition(position);
902   TrackMomentum(momentum);
903   cout << "Edep=" << Edep() << endl;
904   cout << "Etot=" << Etot() << endl;
905   cout << "TrackPid=" << TrackPid() << endl;
906   cout << "TrackCharge=" << TrackCharge() << endl;
907   cout << "TrackMass=" << TrackMass() << endl;
908   cout << "IsTrackOut=" << IsTrackOut() << endl;
909   cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
910   cout << "IsTrackStop=" << IsTrackStop() << endl;
911   cout << "IsTrackAlive=" << IsTrackAlive() << endl;
912   }
913
914   else if((icode >= 100 && icode <= 105) ||
915            (icode == 208) ||
916            (icode == 210) ||
917            (icode == 212) ||
918            (icode >= 214 && icode <= 215) ||
919            (icode == 217) ||
920            (icode == 219) ||
921            (icode == 221) ||
922            (icode == 225) ||
923            (icode == 300) ||
924            (icode == 400)) {
925 // usdraw
926     mreg = GetMreg();
927 //    medium = GetMedium();
928     xsco = GetXsco();
929     ysco = GetYsco();
930     zsco = GetZsco();
931     cout << " icode=" << icode
932          << " mreg=" << mreg
933 //       << " medium=" << medium
934          << " xsco=" << xsco
935          << " ysco=" << ysco
936          << " zsco=" << zsco << endl;
937     cout << "TrackPid=" << TrackPid() << endl;
938     cout << "NSecondaries=" << NSecondaries() << endl;
939     for (Int_t isec=0; isec< NSecondaries(); isec++) {
940 //void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
941 //                 TLorentzVector& position, TLorentzVector& momentum)
942       TFluka::GetSecondary(isec, particleId, position, momentum);
943       cout << "TLorentzVector positionX=" << position.X()
944            << "positionY=" << position.Y()
945            << "positionZ=" << position.Z()
946            << "timeT=" << position.T() << endl;
947       cout << "TLorentzVector momentumX=" << momentum.X()
948            << "momentumY=" << momentum.Y()
949            << "momentumZ=" << momentum.Z()
950            << "energyE=" << momentum.E() << endl;
951       cout << "TrackPid=" << particleId << endl;
952
953     }
954   }
955
956   else if((icode == 19) ||
957           (icode == 29) ||
958           (icode == 39) ||
959           (icode == 49) ||
960           (icode == 59)) {
961     mreg = GetMreg();
962 //    medium = GetMedium();
963     newreg = GetNewreg();
964     xsco = GetXsco();
965     ysco = GetYsco();
966     zsco = GetZsco();
967     cout << " icode=" << icode
968          << " mreg=" << mreg
969 //       << " medium=" << medium
970          << " newreg=" << newreg
971          << " xsco=" << xsco
972          << " ysco=" << ysco
973          << " zsco=" << zsco << endl;
974   }
975 //
976 // ====================================================================
977 //
978
979   
980
981 } // end of FutoTest