3c0d40b99e74bce3fade723ba0183714df509082
[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.17  2003/06/05 10:22:57  morsch
19 All printout under verbosity level control.
20
21 Revision 1.16  2003/03/26 13:30:35  morsch
22 SetTrackIsExiting, SetTrackIsEntering, SetTrackIsInside added.
23
24 Revision 1.15  2003/02/18 16:12:17  morsch
25 Protect  mpdgha against negative argument.
26
27 Revision 1.14  2003/02/18 12:47:59  morsch
28 Gmtod and Gdtom added.
29
30 Revision 1.13  2003/01/31 14:01:51  morsch
31 Major update on
32 - Getters related to geometry.
33 - Communication with run manager (event steering)
34
35 Revision 1.12  2003/01/31 12:18:53  morsch
36 Corrected indices. (E. Futo)
37
38 Revision 1.9  2002/12/06 12:41:29  morsch
39 Mess from last merge cleaned up.
40
41 Revision 1.8  2002/12/06 12:28:44  morsch
42 Region to media mapping corrected and improved.
43
44 Revision 1.7  2002/12/06 12:21:32  morsch
45 User stepping methods added (E. Futo)
46
47 Revision 1.6  2002/11/21 18:40:06  iglez2
48 Media handling added
49
50 Revision 1.5  2002/11/07 17:59:10  iglez2
51 Included the geometry through geant4_vmc/FLUGG
52
53 Revision 1.4  2002/11/04 16:00:46  iglez2
54 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
55
56 Revision 1.3  2002/10/22 15:12:14  alibrary
57 Introducing Riostream.h
58
59 Revision 1.2  2002/10/14 14:57:40  hristov
60 Merging the VirtualMC branch to the main development branch (HEAD)
61
62 Revision 1.1.2.8  2002/10/08 16:33:17  iglez2
63 LSOUIT is set to true before the second call to flukam.
64
65 Revision 1.1.2.7  2002/10/08 09:30:37  iglez2
66 Solved stupid missing ;
67
68 Revision 1.1.2.6  2002/10/07 13:40:22  iglez2
69 First implementations of the PDG <--> Fluka Id conversion routines
70
71 Revision 1.1.2.5  2002/09/26 16:26:03  iglez2
72 Added verbosity
73 Call to gAlice->Generator()->Generate()
74
75 Revision 1.1.2.4  2002/09/26 13:22:23  iglez2
76 Naive implementation of ProcessRun and ProcessEvent
77 Opening/Closing of input file (sInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
78
79 Revision 1.1.2.3  2002/09/20 15:35:51  iglez2
80 Modification of LFDRTR. Value is passed to FLUKA !!!
81
82 Revision 1.1.2.2  2002/09/18 14:34:44  iglez2
83 Revised version with all pure virtual methods implemented
84
85 Revision 1.1.2.1  2002/07/24 08:49:41  alibrary
86 Adding TFluka to VirtualMC
87
88 Revision 1.1  2002/07/05 13:10:07  morsch
89 First commit of Fluka interface.
90
91 */
92
93 #include <Riostream.h>
94
95 #include "TClonesArray.h"
96 #include "TFluka.h"
97 #include "TCallf77.h"      //For the fortran calls
98 #include "Fdblprc.h"       //(DBLPRC) fluka common
99 #include "Fepisor.h"       //(EPISOR) fluka common
100 #include "Ffinuc.h"        //(FINUC) fluka common
101 #include "Fiounit.h"       //(IOUNIT) fluka common
102 #include "Fpaprop.h"       //(PAPROP) fluka common
103 #include "Fpart.h"         //(PART)   fluka common
104 #include "Ftrackr.h"       //(TRACKR) fluka common
105 #include "Fpaprop.h"       //(PAPROP) fluka common
106 #include "Ffheavy.h"       //(FHEAVY) fluka common
107
108 #include "TVirtualMC.h"
109 #include "TG4GeometryManager.h" //For the geometry management
110 #include "TG4DetConstruction.h" //For the detector construction
111
112 #include "FGeometryInit.hh"
113 #include "TLorentzVector.h"
114 #include "FlukaVolume.h"
115
116 // Fluka methods that may be needed.
117 #ifndef WIN32 
118 # define flukam  flukam_
119 # define fluka_openinp fluka_openinp_
120 # define fluka_closeinp fluka_closeinp_
121 # define mcihad mcihad_
122 # define mpdgha mpdgha_
123 #else 
124 # define flukam  FLUKAM
125 # define fluka_openinp FLUKA_OPENINP
126 # define fluka_closeinp FLUKA_CLOSEINP
127 # define mcihad MCIHAD
128 # define mpdgha MPDGHA
129 #endif
130
131 extern "C" 
132 {
133   //
134   // Prototypes for FLUKA functions
135   //
136   void type_of_call flukam(const int&);
137   void type_of_call fluka_openinp(const int&, DEFCHARA);
138   void type_of_call fluka_closeinp(const int&);
139   int  type_of_call mcihad(const int&);
140   int  type_of_call mpdgha(const int&);
141 }
142
143 //
144 // Class implementation for ROOT
145 //
146 ClassImp(TFluka)
147
148 //
149 //----------------------------------------------------------------------------
150 // TFluka constructors and destructors.
151 //____________________________________________________________________________ 
152 TFluka::TFluka()
153   :TVirtualMC(),
154    fVerbosityLevel(0),
155    sInputFileName(""),
156    fDetector(0),
157    fCurrentFlukaRegion(-1)
158
159   //
160   // Default constructor
161   //
162
163  
164 TFluka::TFluka(const char *title, Int_t verbosity)
165   :TVirtualMC("TFluka",title),
166    fVerbosityLevel(verbosity),
167    sInputFileName(""),
168    fTrackIsEntering(0),
169    fTrackIsExiting(0),
170    fDetector(0),
171    fCurrentFlukaRegion(-1)
172 {
173   if (fVerbosityLevel >=3)
174     cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
175
176   
177   // create geometry manager
178   if (fVerbosityLevel >=2)
179     cout << "\t* Creating G4 Geometry manager..." << endl;
180   fGeometryManager = new TG4GeometryManager();
181   if (fVerbosityLevel >=2)
182     cout << "\t* Creating G4 Detector..." << endl;
183   fDetector = new TG4DetConstruction();
184   FGeometryInit* geominit = FGeometryInit::GetInstance();
185   if (geominit)
186     geominit->setDetConstruction(fDetector);
187   else {
188     cerr << "ERROR: Could not create FGeometryInit!" << endl;
189     cerr << "       Exiting!!!" << endl;
190     abort();
191   }
192
193   if (fVerbosityLevel >=3)
194     cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
195
196   fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
197   fNVolumes      = 0;
198   fMediaByRegion = 0;
199 }
200
201 TFluka::~TFluka() {
202   if (fVerbosityLevel >=3)
203     cout << "==> TFluka::~TFluka() destructor called." << endl;
204
205   delete fGeometryManager;
206   fVolumeMediaMap->Delete();
207   delete  fVolumeMediaMap;
208   
209
210   if (fVerbosityLevel >=3)
211     cout << "<== TFluka::~TFluka() destructor called." << endl;
212 }
213
214 //
215 //_____________________________________________________________________________
216 // TFluka control methods
217 //____________________________________________________________________________ 
218 void TFluka::Init() {
219
220   if (fVerbosityLevel >=3)
221     cout << "==> TFluka::Init() called." << endl;
222
223     cout << "\t* InitPhysics() - Prepare input file to be called" << endl; 
224     InitPhysics(); // prepare input file
225     cout << "\t* InitPhysics() - Prepare input file called" << endl; 
226
227   if (fVerbosityLevel >=2)
228     cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
229          << ") in fluka..." << endl;
230   GLOBAL.lfdrtr = true;
231
232   if (fVerbosityLevel >=2)
233     cout << "\t* Opening file " << sInputFileName << endl;
234   const char* fname = sInputFileName;
235   fluka_openinp(lunin, PASSCHARA(fname));
236
237   if (fVerbosityLevel >=2)
238     cout << "\t* Calling flukam..." << endl;
239   flukam(1);
240
241   if (fVerbosityLevel >=2)
242     cout << "\t* Closing file " << sInputFileName << endl;
243   fluka_closeinp(lunin);
244
245   FinishGeometry();
246
247   if (fVerbosityLevel >=3)
248     cout << "<== TFluka::Init() called." << endl;
249
250 }
251
252 void TFluka::FinishGeometry() {
253 //
254 // Build-up table with region to medium correspondance
255 //
256     char tmp[5];
257     
258   if (fVerbosityLevel >=3)
259     cout << "==> TFluka::FinishGeometry() called." << endl;
260
261 //  fGeometryManager->Ggclos();
262
263   FGeometryInit* flugg = FGeometryInit::GetInstance();  
264   
265   fMediaByRegion = new Int_t[fNVolumes+2];
266   for (Int_t i = 0; i < fNVolumes; i++)
267   {
268       FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
269       TString volName = vol->GetName();
270       Int_t   media   = vol->GetMedium();
271       if (fVerbosityLevel >= 3)
272       printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
273       strcpy(tmp, volName.Data());
274       tmp[4] = '\0';
275       flugg->SetMediumFromName(tmp, media, i+1);
276       fMediaByRegion[i] = media;
277   }
278
279   flugg->BuildMediaMap();
280   
281   if (fVerbosityLevel >=3)
282     cout << "<== TFluka::FinishGeometry() called." << endl;
283
284
285 void TFluka::BuildPhysics() {
286   if (fVerbosityLevel >=3)
287     cout << "==> TFluka::BuildPhysics() called." << endl;
288
289
290   if (fVerbosityLevel >=3)
291     cout << "<== TFluka::BuildPhysics() called." << endl;
292 }  
293
294 void TFluka::ProcessEvent() {
295   if (fVerbosityLevel >=3)
296     cout << "==> TFluka::ProcessEvent() called." << endl;
297   fApplication->GeneratePrimaries();
298   EPISOR.lsouit = true;
299   flukam(1);
300   if (fVerbosityLevel >=3)
301     cout << "<== TFluka::ProcessEvent() called." << endl;
302 }
303
304
305 void TFluka::ProcessRun(Int_t nevent) {
306   if (fVerbosityLevel >=3)
307     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
308          << endl;
309
310   if (fVerbosityLevel >=2) {
311     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
312     cout << "\t* Calling flukam again..." << endl;
313   }
314   fApplication->InitGeometry();
315   fApplication->BeginEvent();
316   ProcessEvent();
317   fApplication->FinishEvent();
318   if (fVerbosityLevel >=3)
319     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
320          << endl;
321
322 }
323
324 //_____________________________________________________________________________
325 // methods for building/management of geometry
326 //____________________________________________________________________________ 
327 // functions from GCONS 
328 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,  
329                     Float_t &dens, Float_t &radl, Float_t &absl,
330                     Float_t* ubuf, Int_t& nbuf) {
331 //
332   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
333
334
335 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,  
336                     Double_t &dens, Double_t &radl, Double_t &absl,
337                     Double_t* ubuf, Int_t& nbuf) {
338 //
339   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
340
341
342 // detector composition
343 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
344                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
345                       Float_t* buf, Int_t nwbuf) {
346 //
347   fGeometryManager
348     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); 
349
350 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
351                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
352                       Double_t* buf, Int_t nwbuf) {
353 //
354   fGeometryManager
355     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); 
356
357
358 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
359                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
360 //
361    fGeometryManager
362      ->Mixture(kmat, name, a, z, dens, nlmat, wmat); 
363
364 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
365                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
366 //
367    fGeometryManager
368      ->Mixture(kmat, name, a, z, dens, nlmat, wmat); 
369
370
371 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
372                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
373                     Double_t stemax, Double_t deemax, Double_t epsil, 
374                     Double_t stmin, Float_t* ubuf, Int_t nbuf) { 
375   //
376   fGeometryManager
377     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
378              epsil, stmin, ubuf, nbuf);
379
380 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
381                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
382                     Double_t stemax, Double_t deemax, Double_t epsil, 
383                     Double_t stmin, Double_t* ubuf, Int_t nbuf) { 
384   //
385   fGeometryManager
386     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
387              epsil, stmin, ubuf, nbuf);
388
389
390 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
391                     Double_t thetaY, Double_t phiY, Double_t thetaZ, 
392                     Double_t phiZ) {
393 //                   
394   fGeometryManager
395     ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
396
397
398 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
399 //
400   fGeometryManager->Gstpar(itmed, param, parval); 
401 }    
402
403 // functions from GGEOM 
404 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
405                      Float_t *upar, Int_t np)  {
406 //
407 //  fVolumeMediaMap[TString(name)] = nmed;
408     if (fVerbosityLevel >= 3)
409     printf("TFluka::Gsvolu() name = %s, nmed = %d\n", name, nmed);
410     
411     TClonesArray &lvols = *fVolumeMediaMap;
412     new(lvols[fNVolumes++]) 
413         FlukaVolume(name, nmed);
414     return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); 
415 }
416 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
417                      Double_t *upar, Int_t np)  {
418 //
419     TClonesArray &lvols = *fVolumeMediaMap;
420     new(lvols[fNVolumes++]) 
421         FlukaVolume(name, nmed);
422
423     return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); 
424 }
425  
426 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
427                    Int_t iaxis) {
428 //
429 //  The medium of the daughter is the one of the mother
430     Int_t volid = TFluka::VolId(mother);
431     Int_t med   = TFluka::VolId2Mate(volid);
432     TClonesArray &lvols = *fVolumeMediaMap;
433     new(lvols[fNVolumes++]) 
434         FlukaVolume(name, med);
435     fGeometryManager->Gsdvn(name, mother, ndiv, iaxis); 
436
437
438 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
439                     Int_t iaxis, Double_t c0i, Int_t numed) {
440 //
441     TClonesArray &lvols = *fVolumeMediaMap;
442     new(lvols[fNVolumes++]) 
443         FlukaVolume(name, numed);
444     fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
445
446
447 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
448                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
449 //      
450     TClonesArray &lvols = *fVolumeMediaMap;
451     new(lvols[fNVolumes++]) 
452         FlukaVolume(name, numed);               
453     fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
454
455
456 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
457                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { 
458 //
459     TClonesArray &lvols = *fVolumeMediaMap;
460     new(lvols[fNVolumes++]) 
461         FlukaVolume(name, numed);
462     fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
463
464
465 void TFluka::Gsord(const char *name, Int_t iax) {
466 //
467   fGeometryManager->Gsord(name, iax); 
468
469
470 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
471                    Double_t x, Double_t y, Double_t z, Int_t irot, 
472                    const char *konly) {
473 //
474   fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly); 
475
476
477 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
478                     Double_t x, Double_t y, Double_t z, Int_t irot,
479                     const char *konly, Float_t *upar, Int_t np)  {
480   //
481   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
482
483 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
484                     Double_t x, Double_t y, Double_t z, Int_t irot,
485                     const char *konly, Double_t *upar, Int_t np)  {
486   //
487   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
488
489
490 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
491 //
492   fGeometryManager->Gsbool(onlyVolName, manyVolName);
493 }
494
495 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
496                          Float_t *absco, Float_t *effic, Float_t *rindex) {
497 //
498   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
499 }  
500 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
501                          Double_t *absco, Double_t *effic, Double_t *rindex) {
502 //
503   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
504 }  
505     
506 // Euclid
507 void TFluka::WriteEuclid(const char* fileName, const char* topVol, 
508                           Int_t number, Int_t nlevel) {
509 //
510   fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel); 
511
512
513
514
515 //_____________________________________________________________________________
516 // methods needed by the stepping
517 //____________________________________________________________________________ 
518
519 Int_t TFluka::GetMedium() const {
520 //
521 //  Get the medium number for the current fluka region
522 //
523     FGeometryInit* flugg = FGeometryInit::GetInstance();  
524     return flugg->GetMedium(fCurrentFlukaRegion);
525 }
526
527
528
529 //____________________________________________________________________________ 
530 // particle table usage
531 // ID <--> PDG transformations
532 //_____________________________________________________________________________
533 Int_t TFluka::IdFromPDG(Int_t pdg) const 
534 {
535   //
536   // Return Fluka code from PDG and pseudo ENDF code
537
538   // MCIHAD() goes from pdg to fluka internal.
539   Int_t intfluka = mcihad(pdg);
540   // KPTOIP array goes from internal to official
541   return GetFlukaKPTOIP(intfluka);
542 }
543
544 Int_t TFluka::PDGFromId(Int_t id) const 
545 {
546   //
547   // Return PDG code and pseudo ENDF code from Fluka code
548
549   //IPTOKP array goes from official to internal
550     if (id == 0) {
551       if (fVerbosityLevel >= 1)
552         printf("PDGFromId: Error id = 0");
553         return -1;
554     }
555     
556   Int_t intfluka = GetFlukaIPTOKP(id);
557     if (intfluka == 0) {
558       if (fVerbosityLevel >= 1)
559         printf("PDGFromId: Error intfluka = 0");
560         return -1;
561     } else if (intfluka < 0) {
562       if (fVerbosityLevel >= 1)
563         printf("PDGFromId: Error intfluka < 0");
564         return -1;
565     }
566     if (fVerbosityLevel >= 3)
567     printf("mpdgha called with %d %d \n", id, intfluka);
568   return mpdgha(intfluka);
569 }
570
571 //_____________________________________________________________________________
572 // methods for physics management
573 //____________________________________________________________________________ 
574 //
575 // set methods
576 //
577
578 void TFluka::SetProcess(const char* flagName, Int_t flagValue)
579 {
580   Int_t i;
581   if (iNbOfProc < 100) {
582     for (i=0; i<iNbOfProc; i++) {
583       if (strcmp(&sProcessFlag[i][0],flagName) == 0) {
584         iProcessValue[iNbOfProc] = flagValue;
585         goto fin;
586       }
587     }
588     strcpy(&sProcessFlag[iNbOfProc][0],flagName);
589     iProcessValue[iNbOfProc++] = flagValue;
590   }
591   else
592     cout << "Nb of SetProcess calls exceeds 100 - ignored" << endl;
593 fin:
594   iNbOfProc = iNbOfProc;
595 }
596
597 void TFluka::SetCut(const char* cutName, Double_t cutValue)
598 {
599   Int_t i;
600   if (iNbOfCut < 100) {
601     for (i=0; i<iNbOfCut; i++) {
602       if (strcmp(&sCutFlag[i][0],cutName) == 0) {
603         fCutValue[iNbOfCut] = cutValue;
604         goto fin;
605       }
606     }
607     strcpy(&sCutFlag[iNbOfCut][0],cutName);
608     fCutValue[iNbOfCut++] = cutValue;
609   }
610   else
611     cout << "Nb of SetCut calls exceeds 100 - ignored" << endl;
612 fin:
613   iNbOfCut = iNbOfCut;
614 }
615
616 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
617 {
618   printf("WARNING: Xsec not yet implemented !\n"); return -1.;
619 }
620
621
622 void TFluka::InitPhysics()
623 {
624 // Last material number taken from the "corealice.inp" file, presently 31
625 // !!! it should be available from Flugg !!!
626   Float_t fLastMaterial = 31.0;
627   Float_t fLastRegion   = 692.;
628   
629 // construct file names
630   TString sAliceInp = getenv("ALICE_ROOT");
631   sAliceInp +="/TFluka/input/";
632   TString sAliceCoreInp = sAliceInp;
633   sAliceInp += GetInputFileName();
634   sAliceCoreInp += GetCoreInputFileName();
635   ifstream AliceCoreInp(sAliceCoreInp.Data());
636   ofstream AliceInp(sAliceInp.Data());
637
638 // copy core input file until (not included) START card
639   Char_t sLine[255];
640   Float_t fEventsPerRun;
641   while (AliceCoreInp.getline(sLine,255)) {
642     if (strncmp(sLine,"START",5) != 0)
643       AliceInp << sLine << endl;
644     else {
645       sscanf(sLine+10,"%10f",&fEventsPerRun);
646       goto fin;
647     }
648   } //end of while
649
650 fin:
651 // in G3 the process control values meaning can be different for
652 // different processes, but for most of them is:
653 //   0  process is not activated
654 //   1  process is activated WITH generation of secondaries
655 //   2  process is activated WITHOUT generation of secondaries
656 // if process does not generate secondaries => 1 same as 2
657 //
658 // Exceptions:
659 //   MULS:  also 3
660 //   LOSS:  also 3, 4
661 //   RAYL:  only 0,1
662 //   HADR:  may be > 2
663 //
664  
665 // Loop over number of SetProcess calls  
666   AliceInp << "*----------------------------------------------------------------------------- "; 
667   AliceInp << endl;
668   AliceInp << "*----- The following data are generated from SetProcess and SetCut calls ----- "; 
669   AliceInp << endl;
670   AliceInp << "*----------------------------------------------------------------------------- "; 
671     AliceInp << endl;
672   for (Int_t i=0; i<iNbOfProc; i++) {
673
674     // annihilation
675     // G3 default value: 1
676     // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
677     // Particles: e+
678     // Physics:   EM
679     // gMC ->SetProcess("ANNI",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. ANNH-THR
680     if ((strncmp(&sProcessFlag[i][0],"ANNI",4) == 0) && iProcessValue[i] == 1) {
681       AliceInp << "*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0."; 
682       AliceInp << endl;
683       AliceInp << "*Generated from call: SetProcess('ANNI',1);"; 
684       AliceInp << endl;
685       AliceInp << setw(10) << "EMFCUT    "; 
686       AliceInp << setiosflags(ios::scientific) << setprecision(5);
687       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
688       AliceInp << setw(10) << -1.0; // kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
689       AliceInp << setw(10) << 0.0;  // not used
690       AliceInp << setw(10) << 0.0;  // not used
691       AliceInp << setw(10) << 3.0;  // lower bound of the material indices in which the respective thresholds apply
692       AliceInp << setw(10) << setprecision(2);
693       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
694       AliceInp << setprecision(1);
695       AliceInp << setw(10) << 1.0;  // step length in assigning indices
696       AliceInp << setw(8)  << "ANNH-THR"; 
697       AliceInp << endl;
698     }
699     
700     // bremsstrahlung
701     // G3 default value: 1
702     // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
703     //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
704     //               G4LowEnergyBremstrahlung
705     // Particles: e-/e+; mu+/mu-
706     // Physics:   EM
707     // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
708                                  // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
709     else if ((strncmp(&sProcessFlag[i][0],"BREM",4) == 0) && iProcessValue[i] == 1) {
710       AliceInp << "*Bremsstrahlung by muons and charged hadrons is activated"; 
711       AliceInp << endl;
712       AliceInp << "*Generated from call: SetProcess('BREM',1);"; 
713       AliceInp << endl;
714       AliceInp << setw(10) << "PAIRBREM  "; 
715       AliceInp << setiosflags(ios::scientific) << setprecision(5);
716       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
717       AliceInp << setw(10) << 2.0; // bremsstrahlung by muons and charged hadrons is activated
718       AliceInp << setw(10) << 0.0; // e+, e- kinetic energy threshold (in GeV) for explicit pair production. A value of 0.0 is meaningful.
719       AliceInp << setw(10) << 0.0; // no explicit bremsstrahlung production is simulated
720       AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
721       AliceInp << setw(10) << setprecision(2);
722       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
723       AliceInp << endl;
724       // for e+ and e-
725       AliceInp << "*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0."; 
726       AliceInp << endl;
727       AliceInp << "*Generated from call: SetProcess('BREM',1);"; 
728       AliceInp << endl;
729       AliceInp << setw(10) << "EMFCUT    "; 
730       AliceInp << setiosflags(ios::scientific) << setprecision(5);
731       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
732       AliceInp << setw(10) << -1.0; // kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
733       AliceInp << setw(10) << 0.0;  // not used
734       AliceInp << setw(10) << 0.0;  // not used
735       AliceInp << setw(10) << 3.0;  // lower bound of the material indices in which the respective thresholds apply
736       AliceInp << setw(10) << setprecision(2);
737       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
738       AliceInp << setprecision(1);
739       AliceInp << setw(10) << 1.0; // step length in assigning indices
740       AliceInp << setw(8)  << "ELPO-THR"; 
741       AliceInp << endl;
742     }
743     
744     // Compton scattering
745     // G3 default value: 1
746     // G4 processes: G4ComptonScattering,
747     //               G4LowEnergyCompton,
748     //               G4PolarizedComptonScattering
749     // Particles: gamma
750     //                                                                               // Physics:   EM
751     // gMC ->SetProcess("COMP",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. PHOT-THR
752     else if ((strncmp(&sProcessFlag[i][0],"COMP",4) == 0) && iProcessValue[i] == 1) {
753       AliceInp << "*Energy threshold (GeV) for Compton scattering - resets to default=0."; 
754       AliceInp << endl;
755       AliceInp << "*Generated from call: SetProcess('COMP',1);"; 
756       AliceInp << endl;
757       AliceInp << setw(10) << "EMFCUT    "; 
758       AliceInp << setiosflags(ios::scientific) << setprecision(5);
759       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
760       AliceInp << setw(10) << -1.0; // energy threshold (GeV) for Compton scattering - resets to default=0.
761       AliceInp << setw(10) << 0.0;  // not used
762       AliceInp << setw(10) << 0.0;  // not used
763       AliceInp << setw(10) << 3.0;  // lower bound of the material indices in which the respective thresholds apply
764       AliceInp << setprecision(2);
765       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
766       AliceInp << setprecision(1);
767       AliceInp << setw(10) << 1.0; // step length in assigning indices
768       AliceInp << setw(8) << "PHOT-THR"; 
769       AliceInp << endl;
770     }
771
772     // decay
773     // G3 default value: 1
774     // G4 process: G4Decay
775     // 
776     // Particles: all which decay is applicable for
777     // Physics:   General
778     //gMC ->SetProcess("DCAY",1); // not available
779     else if ((strncmp(&sProcessFlag[i][0],"DCAY",4) == 0) && iProcessValue[i] == 1) 
780       cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not avaliable!" << endl;
781       
782     // delta-ray
783     // G3 default value: 2
784     // !! G4 treats delta rays in different way
785     // G4 processes: G4eIonisation/G4IeIonization,
786     //               G4MuIonisation/G4IMuIonization,
787     //               G4hIonisation/G4IhIonisation
788     //                                                                               // Particles: charged
789     // Physics:   EM
790     // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0.  0. 3. lastmat 0.
791     else if ((strncmp(&sProcessFlag[i][0],"DRAY",4) == 0) && iProcessValue[i] == 0) {
792       AliceInp << "*Kinetic energy threshold (GeV) for delta ray production"; 
793       AliceInp << endl;
794       AliceInp << "*Generated from call: SetProcess('DRAY',1);"; 
795       AliceInp << endl;
796       AliceInp << setw(10) << "DELTARAY  "; 
797       AliceInp << setiosflags(ios::scientific) << setprecision(5);
798       AliceInp << setw(10) << 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
799       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
800       AliceInp << setw(10) << 0.0; // ignored
801       AliceInp << setw(10) << 0.0; // ignored
802       AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
803       AliceInp << setw(10) << setprecision(2);
804       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
805       AliceInp << setprecision(1);
806       AliceInp << setw(10) << 1.0; // step length in assigning indices
807       AliceInp << endl;
808     }
809     
810     // muon nuclear interaction
811     // G3 default value: 0
812     // G4 processes: G4MuNuclearInteraction,
813     // G4MuonMinusCaptureAtRest
814     // 
815     // Particles: mu
816     // Physics:   Not set
817     // gMC ->SetProcess("MUNU",1); // MUPHOTON  1.   0.  0. 3. lastmat
818     else if ((strncmp(&sProcessFlag[i][0],"MUNU",4) == 0) && iProcessValue[i] == 1) {
819       AliceInp << "*Muon nuclear interactions with production of secondary hadrons"; 
820       AliceInp << endl;
821       AliceInp << "*Generated from call: SetProcess('MUNU',1);"; 
822       AliceInp << endl;
823       AliceInp << setw(10) << "MUPHOTON  "; 
824       AliceInp << setiosflags(ios::scientific) << setprecision(5);
825       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
826       AliceInp << setw(10) << 1.0;  // full simulation of muon nuclear interactions and production of secondary hadrons
827       AliceInp << setw(10) << 0.0; // ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
828       AliceInp << setw(10) << 0.0; // fraction of rho-like interactions ( must be < 1) - Default = 0.75.
829       AliceInp << setprecision(1);
830       AliceInp << setw(10) << 3.0;  // lower bound of the material indices in which the respective thresholds apply
831       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
832       AliceInp << setprecision(1);
833       AliceInp << setw(10) << 1.0;  // step length in assigning indices
834       AliceInp << endl;
835     }
836
837     // muon nuclear interaction
838     // G3 default value: 0
839     // G4 processes: G4MuNuclearInteraction,
840     // G4MuonMinusCaptureAtRest
841     // 
842     // Particles: mu
843     // Physics:   Not set
844     // gMC ->SetProcess("MUNU",1); // MUPHOTON  1.   0.  0. 3. lastmat
845     else if ((strncmp(&sProcessFlag[i][0],"MUNU",4) == 0) && iProcessValue[i] == 2) {
846       AliceInp << "*Muon nuclear interactions without production of secondary hadrons"; 
847       AliceInp << endl;
848       AliceInp << "*Generated from call: SetProcess('MUNU',2);"; 
849       AliceInp << endl;
850       AliceInp << setw(10) << "MUPHOTON  "; 
851       AliceInp << setiosflags(ios::scientific) << setprecision(5);
852       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
853       AliceInp << setw(10) << 2.0; // full simulation of muon nuclear interactions and production of secondary hadrons
854       AliceInp << setw(10) << 0.0; // ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
855       AliceInp << setw(10) << 0.0; // fraction of rho-like interactions ( must be < 1) - Default = 0.75.
856       AliceInp << setprecision(1);
857       AliceInp << setw(10) << 3.0;  // lower bound of the material indices in which the respective thresholds apply
858       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
859       AliceInp << setprecision(1);
860       AliceInp << setw(10) << 1.0;  // step length in assigning indices
861       AliceInp << endl;
862     }
863
864   // pair production
865     // G3 default value: 1
866     // G4 processes: G4GammaConversion,
867     //               G4MuPairProduction/G4IMuPairProduction
868     //               G4LowEnergyGammaConversion
869     // Particles: gamma, mu
870     // Physics:   EM
871   // gMC ->SetProcess("PAIR",1); // PAIRBREM  1.   0.  0. 3. lastmat
872                                  // EMFCUT    0.   0. -1. 3. lastmat 0. PHOT-THR
873     else if ((strncmp(&sProcessFlag[i][0],"PAIR",4) == 0) && iProcessValue[i] == 1) {
874       AliceInp << "*Pair production by muons and charged hadrons is activated"; 
875       AliceInp << endl;
876       AliceInp << "*Generated from call: SetProcess('PAIR',1);"; 
877       AliceInp << endl;
878       AliceInp << setw(10) << "PAIRBREM  "; 
879       AliceInp << setiosflags(ios::scientific) << setprecision(5);
880       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
881       AliceInp << setw(10) << 1.0; // pair production by muons and charged hadrons is activated
882       AliceInp << setw(10) << 0.0; // e+, e- kinetic energy threshold (in GeV) for explicit pair production.
883       AliceInp << setw(10) << 0.0; // no explicit bremsstrahlung production is simulated
884       AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
885       AliceInp << setprecision(2);
886       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
887       AliceInp << endl;
888       // for e+ and e-
889       AliceInp << "*Pair production by electrons is activated"; 
890       AliceInp << endl;
891       AliceInp << "*Generated from call: SetProcess('PAIR',1);"; 
892       AliceInp << endl;
893       AliceInp << setw(10) << "EMFCUT    "; 
894       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
895       AliceInp << setw(10) << 0.0;  // ignored
896       AliceInp << setw(10) << 0.0;  // ignored
897       AliceInp << setw(10) << -1.0; // resets to default=0.
898       AliceInp << setw(10) << 3.0;  // lower bound of the material indices in which the respective thresholds apply
899       AliceInp << setprecision(2);
900       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
901       AliceInp << setprecision(1);
902       AliceInp << setw(10) << 1.0;  // step length in assigning indices
903       AliceInp << setw(8) << "PHOT-THR"; 
904       AliceInp << endl;
905     }
906
907     // photofission
908     // G3 default value: 0
909     // G4 process: ??
910     //
911     // Particles: gamma
912     // Physics:   ??
913     // gMC ->SetProcess("PFIS",0); // PHOTONUC -1.   0.  0. 3. lastmat 0.
914     else if ((strncmp(&sProcessFlag[i][0],"PFIS",4) == 0) && iProcessValue[i] == 0) {
915       AliceInp << "*No photonuclear interactions";
916       AliceInp << endl;
917       AliceInp << "*Generated from call: SetProcess('PFIS',0);"; 
918       AliceInp << endl;
919       AliceInp << setw(10) << "PHOTONUC  "; 
920       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
921       AliceInp << setw(10) << -1.0; // no photonuclear interactions
922       AliceInp << setw(10) << 0.0;  // not used
923       AliceInp << setw(10) << 0.0;  // not used
924       AliceInp << setw(10) << 3.0;  // upper bound of the material indices in which the respective thresholds apply
925       AliceInp << setprecision(2); 
926       AliceInp << setw(10) << fLastMaterial;
927       AliceInp << setprecision(1);  // upper bound of the material indices in which the respective thresholds apply
928       AliceInp << setprecision(1);
929       AliceInp << setw(10) << 1.0;  // step length in assigning indices
930       AliceInp << endl;
931     }
932
933     else if ((strncmp(&sProcessFlag[i][0],"PFIS",4) == 0) && iProcessValue[i] == 1) {
934       AliceInp << "*Photon nuclear interactions are activated at all energies";
935       AliceInp << endl;
936       AliceInp << "*Generated from call: SetProcess('PFIS',1);"; 
937       AliceInp << endl;
938       AliceInp << setw(10) << "PHOTONUC  "; 
939       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
940       AliceInp << setw(10) << 1.0; // photonuclear interactions are activated at all energies
941       AliceInp << setw(10) << 0.0; // not used
942       AliceInp << setw(10) << 0.0; // not used
943       AliceInp << setprecision(2); 
944       AliceInp << setw(10) << 3.0; // upper bound of the material indices in which the respective thresholds apply
945       AliceInp << setw(10) << fLastMaterial;
946       AliceInp << setprecision(1); // upper bound of the material indices in which the respective thresholds apply
947       AliceInp << setprecision(1);
948       AliceInp << setw(10) << 1.0; // step length in assigning indices
949       AliceInp << endl;
950     }
951
952     // photo electric effect
953     // G3 default value: 1
954     // G4 processes: G4PhotoElectricEffect
955     //               G4LowEnergyPhotoElectric
956     // Particles: gamma
957     // Physics:   EM
958     // gMC ->SetProcess("PHOT",1); // EMFCUT    0.  -1.  0. 3. lastmat 0. PHOT-THR
959     else if ((strncmp(&sProcessFlag[i][0],"PHOT",4) == 0) && iProcessValue[i] == 1) {
960       AliceInp << "*Photo electric effect is activated"; 
961       AliceInp << endl;
962       AliceInp << "*Generated from call: SetProcess('PHOT',1);"; 
963       AliceInp << endl;
964       AliceInp << setw(10) << "EMFCUT    "; 
965       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
966       AliceInp << setw(10) << 0.0;  // ignored
967       AliceInp << setw(10) << -1.0; // resets to default=0.
968       AliceInp << setw(10) << 0.0;  // ignored
969       AliceInp << setw(10) << 3.0;  // upper bound of the material indices in which the respective thresholds apply
970       AliceInp << setprecision(2);
971       AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
972       AliceInp << setprecision(1);
973       AliceInp << setw(10) << 1.0;  // step length in assigning indices
974       AliceInp << setw(8) << "PHOT-THR"; 
975       AliceInp << endl;
976     }
977
978     else { // processes not yet treated
979     //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
980             
981     // Cerenkov photon generation
982     // G3 default value: 0
983     // G4 process: G4Cerenkov
984     // 
985     // Particles: charged
986     // Physics:   Optical
987     //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
988             
989     //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
990
991     // hadronic process
992     // G3 default value: 1
993     // G4 processes: all defined by TG4PhysicsConstructorHadron
994     //  
995     // Particles: hadrons
996     // Physics:   Hadron
997     // gMC ->SetProcess("HADR",1); // ??? hadronic process
998
999     // light photon absorption
1000     // it is turned on when Cerenkov process is turned on
1001     // G3 default value: 0
1002     // G4 process: G4OpAbsorption, G4OpBoundaryProcess
1003     // 
1004     // Particles: optical photon
1005     // Physics:   Optical
1006     // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
1007
1008     // energy loss
1009     // G3 default value: 2
1010     // G4 processes: G4eIonisation/G4IeIonization,
1011     //               G4MuIonisation/G4IMuIonization,
1012     //               G4hIonisation/G4IhIonisation
1013     // 
1014     // Particles: charged
1015     // Physics:   EM
1016     // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
1017         
1018     // multiple scattering
1019     // G3 default value: 1
1020     // G4 process: G4MultipleScattering/G4IMultipleScattering
1021     // 
1022     // Particles: charged
1023     // Physics:   EM
1024     // gMC ->SetProcess("MULS",1); // ??? MULSOPT  ? multiple scattering
1025
1026     // Rayleigh scattering
1027     // G3 default value: 0
1028     // G4 process: G4OpRayleigh
1029     // 
1030     // Particles: optical photon
1031     // Physics:   Optical
1032     //xx gMC ->SetProcess("RAYL",1);
1033         
1034     //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
1035
1036     // synchrotron radiation in magnetic field
1037     // G3 default value: 0
1038     // G4 process: G4SynchrotronRadiation
1039     // 
1040     // Particles: ??
1041     // Physics:   Not set
1042     //xx gMC ->SetProcess("SYNC",1); // ??? synchrotron radiation generation
1043
1044       cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not yet implemented!" << endl;
1045     }
1046   } //end of loop number of SetProcess calls
1047
1048  
1049 // Loop over number of SetCut calls  
1050   for (Int_t i=0; i<iNbOfCut; i++) {
1051
1052     // gammas
1053     // G4 particles: "gamma"
1054     // G3 default value: 0.001 GeV
1055     //gMC ->SetCut("CUTGAM",cut); // cut for gammas
1056     if (strncmp(&sCutFlag[i][0],"CUTGAM",6) == 0) {
1057       AliceInp << "*Cut for gamma"; 
1058       AliceInp << endl;
1059       AliceInp << "*Generated from call: SetCut('CUTGAM',cut);"; 
1060       AliceInp << endl;
1061       AliceInp << setw(10) << "PART-THR  "; 
1062       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1063       AliceInp << setw(10) << -fCutValue[i];
1064       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1065       AliceInp << setw(10) << 7.0;
1066       AliceInp << endl;
1067     }
1068
1069     // electrons
1070     // G4 particles: "e-"
1071     // ?? positrons
1072     // G3 default value: 0.001 GeV
1073     //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
1074     else if (strncmp(&sCutFlag[i][0],"CUTELE",6) == 0) {
1075       AliceInp << "*Cut for electrons"; 
1076       AliceInp << endl;
1077       AliceInp << "*Generated from call: SetCut('CUTELE',cut);"; 
1078       AliceInp << endl;
1079       AliceInp << setw(10) << "PART-THR  "; 
1080       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1081       AliceInp << setw(10) << -fCutValue[i];
1082       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1083       AliceInp << setw(10) << 3.0;
1084       AliceInp << setw(10) << 4.0;
1085       AliceInp << setw(10) << 1.0;
1086       AliceInp << endl;
1087     }
1088
1089     // neutral hadrons
1090     // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
1091     // G3 default value: 0.01 GeV
1092     //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
1093     else if (strncmp(&sCutFlag[i][0],"CUTNEU",6) == 0) {
1094       AliceInp << "*Cut for neutral hadrons"; 
1095       AliceInp << endl;
1096       AliceInp << "*Generated from call: SetCut('CUTNEU',cut);"; 
1097       AliceInp << endl;
1098       AliceInp << setw(10) << "PART-THR  "; 
1099       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1100       AliceInp << setw(10) << -fCutValue[i];
1101       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1102       AliceInp << setw(10) << 8.0; // Neutron
1103       AliceInp << setw(10) << 9.0; // Antineutron
1104       AliceInp << endl;
1105       AliceInp << setw(10) << "PART-THR  "; 
1106       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1107       AliceInp << setw(10) << -fCutValue[i];
1108       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1109       AliceInp << setw(10) << 12.0; // Kaon zero long
1110       AliceInp << setw(10) << 12.0; // Kaon zero long
1111       AliceInp << endl;
1112       AliceInp << setw(10) << "PART-THR  "; 
1113       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1114       AliceInp << setw(10) << -fCutValue[i];
1115       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1116       AliceInp << setw(10) << 17.0; // Lambda, 18=Antilambda
1117       AliceInp << setw(10) << 19.0; // Kaon zero short
1118       AliceInp << endl;
1119       AliceInp << setw(10) << "PART-THR  "; 
1120       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1121       AliceInp << setw(10) << -fCutValue[i];
1122       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1123       AliceInp << setw(10) << 22.0; // Sigma zero, Pion zero, Kaon zero
1124       AliceInp << setw(10) << 25.0; // Antikaon zero
1125       AliceInp << endl;
1126       AliceInp << setw(10) << "PART-THR  "; 
1127       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1128       AliceInp << setw(10) << -fCutValue[i];
1129       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1130       AliceInp << setw(10) << 32.0; // Antisigma zero
1131       AliceInp << setw(10) << 32.0; // Antisigma zero
1132       AliceInp << endl;
1133       AliceInp << setw(10) << "PART-THR  "; 
1134       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1135       AliceInp << setw(10) << -fCutValue[i];
1136       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1137       AliceInp << setw(10) << 34.0; // Xi zero
1138       AliceInp << setw(10) << 35.0; // AntiXi zero
1139       AliceInp << endl;
1140       AliceInp << setw(10) << "PART-THR  "; 
1141       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1142       AliceInp << setw(10) << -fCutValue[i];
1143       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1144       AliceInp << setw(10) << 47.0; // D zero
1145       AliceInp << setw(10) << 48.0; // AntiD zero
1146       AliceInp << endl;
1147       AliceInp << setw(10) << "PART-THR  "; 
1148       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1149       AliceInp << setw(10) << -fCutValue[i];
1150       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1151       AliceInp << setw(10) << 53.0; // Xi_c zero
1152       AliceInp << setw(10) << 53.0; // Xi_c zero
1153       AliceInp << endl;
1154       AliceInp << setw(10) << "PART-THR  "; 
1155       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1156       AliceInp << setw(10) << -fCutValue[i];
1157       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1158       AliceInp << setw(10) << 55.0; // Xi'_c zero
1159       AliceInp << setw(10) << 56.0; // Omega_c zero
1160       AliceInp << endl;
1161       AliceInp << setw(10) << "PART-THR  "; 
1162       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1163       AliceInp << setw(10) << -fCutValue[i];
1164       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1165       AliceInp << setw(10) << 59.0; // AntiXi_c zero
1166       AliceInp << setw(10) << 59.0; // AntiXi_c zero
1167       AliceInp << endl;
1168       AliceInp << setw(10) << "PART-THR  "; 
1169       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1170       AliceInp << setw(10) << -fCutValue[i];
1171       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1172       AliceInp << setw(10) << 61.0; // AntiXi'_c zero
1173       AliceInp << setw(10) << 62.0; // AntiOmega_c zero
1174       AliceInp << endl;
1175     }
1176
1177     // charged hadrons
1178     // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
1179     // G3 default value: 0.01 GeV
1180     //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
1181     else if (strncmp(&sCutFlag[i][0],"CUTHAD",6) == 0) {
1182       AliceInp << "*Cut for charged hadrons"; 
1183       AliceInp << endl;
1184       AliceInp << "*Generated from call: SetCut('CUTHAD',cut);"; 
1185       AliceInp << endl;
1186       AliceInp << setw(10) << "PART-THR  "; 
1187       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1188       AliceInp << setw(10) << -fCutValue[i];
1189       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1190       AliceInp << setw(10) << 1.0; // Proton
1191       AliceInp << setw(10) << 2.0; // Antiproton
1192       AliceInp << endl;
1193       AliceInp << setw(10) << "PART-THR  "; 
1194       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1195       AliceInp << setw(10) << -fCutValue[i];
1196       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1197       AliceInp << setw(10) << 13.0; // Positive Pion, Negative Pion, Positive Kaon
1198       AliceInp << setw(10) << 16.0; // Negative Kaon
1199       AliceInp << endl;
1200       AliceInp << setw(10) << "PART-THR  "; 
1201       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1202       AliceInp << setw(10) << -fCutValue[i];
1203       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1204       AliceInp << setw(10) << 20.0; // Negative Sigma
1205       AliceInp << setw(10) << 16.0; // Positive Sigma
1206       AliceInp << endl;
1207       AliceInp << setw(10) << "PART-THR  "; 
1208       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1209       AliceInp << setw(10) << -fCutValue[i];
1210       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1211       AliceInp << setw(10) << 31.0; // Antisigma minus
1212       AliceInp << setw(10) << 33.0; // Antisigma plus
1213       AliceInp << setprecision(1);
1214       AliceInp << setw(10) << 2.0;  // step length
1215       AliceInp << endl;
1216       AliceInp << setw(10) << "PART-THR  "; 
1217       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1218       AliceInp << setw(10) << -fCutValue[i];
1219       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1220       AliceInp << setw(10) << 36.0; // Negative Xi, Positive Xi, Omega minus
1221       AliceInp << setw(10) << 39.0; // Antiomega
1222       AliceInp << endl;
1223       AliceInp << setw(10) << "PART-THR  "; 
1224       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1225       AliceInp << setw(10) << -fCutValue[i];
1226       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1227       AliceInp << setw(10) << 45.0; // D plus
1228       AliceInp << setw(10) << 46.0; // D minus
1229       AliceInp << endl;
1230       AliceInp << setw(10) << "PART-THR  "; 
1231       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1232       AliceInp << setw(10) << -fCutValue[i];
1233       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1234       AliceInp << setw(10) << 49.0; // D_s plus, D_s minus, Lambda_c plus
1235       AliceInp << setw(10) << 52.0; // Xi_c plus
1236       AliceInp << endl;
1237       AliceInp << setw(10) << "PART-THR  "; 
1238       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1239       AliceInp << setw(10) << -fCutValue[i];
1240       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1241       AliceInp << setw(10) << 54.0; // Xi'_c plus
1242       AliceInp << setw(10) << 60.0; // AntiXi'_c minus
1243       AliceInp << setprecision(1);
1244       AliceInp << setw(10) << 6.0;  // step length
1245       AliceInp << endl;
1246       AliceInp << setw(10) << "PART-THR  "; 
1247       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1248       AliceInp << setw(10) << -fCutValue[i];
1249       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1250       AliceInp << setw(10) << 57.0; // Antilambda_c minus
1251       AliceInp << setw(10) << 58.0; // AntiXi_c minus
1252       AliceInp << endl;
1253     }
1254
1255     // muons
1256     // G4 particles: "mu+", "mu-"
1257     // G3 default value: 0.01 GeV
1258     //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
1259     else if (strncmp(&sCutFlag[i][0],"CUTMUO",6) == 0) {
1260       AliceInp << "*Cut for muons"; 
1261       AliceInp << endl;
1262       AliceInp << "*Generated from call: SetCut('CUTMUO',cut);"; 
1263       AliceInp << endl;
1264       AliceInp << setw(10) << "PART-THR  "; 
1265       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1266       AliceInp << setw(10) << -fCutValue[i];
1267       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1268       AliceInp << setprecision(2);
1269       AliceInp << setw(10) << 10.0;
1270       AliceInp << setw(10) << 11.0;
1271       AliceInp << endl;
1272     }
1273
1274     // electron bremsstrahlung
1275     // G4 particles: "gamma"
1276     // G3 default value: CUTGAM=0.001 GeV
1277     //gMC ->SetCut("BCUTE",cut);  // cut for electron bremsstrahlung
1278     else if (strncmp(&sCutFlag[i][0],"BCUTE",5) == 0) {
1279       AliceInp << "*Cut for electron bremsstrahlung"; 
1280       AliceInp << endl;
1281       AliceInp << "*Generated from call: SetCut('BCUTE',cut);"; 
1282       AliceInp << endl;
1283       AliceInp << setw(10) << "EMFCUT    "; 
1284       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1285       AliceInp << setw(10) << -fCutValue[i];
1286       AliceInp << setw(10) << setiosflags(ios::fixed);
1287       AliceInp << setw(10) << setprecision(1);
1288       AliceInp << setw(10) << 0.0; // photon cut-off is unchanged
1289       AliceInp << setw(10) << 0.0; // ignored
1290       AliceInp << setw(10) << 2.0;
1291       AliceInp << setprecision(4);
1292       AliceInp << setw(10) << fLastRegion; // upper bound of the material indices in which the respective thresholds apply
1293       AliceInp << setprecision(1);
1294       AliceInp << setw(10) << 1.0; // step length in assigning indices
1295       AliceInp << endl;
1296     }
1297
1298     // muon and hadron bremsstrahlung
1299     // G4 particles: "gamma"
1300     // G3 default value: CUTGAM=0.001 GeV
1301     //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung ????????????
1302     else if (strncmp(&sCutFlag[i][0],"BCUTM",5) == 0) {
1303       AliceInp << "*Cut for muon and hadron bremsstrahlung ????????????"; 
1304       AliceInp << endl;
1305       AliceInp << "*Generated from call: SetCut('BCUTM',cut);"; 
1306       AliceInp << endl;
1307       AliceInp << setw(10) << "PAIRBREM  "; 
1308       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1309       AliceInp << setw(10) << -fCutValue[i];
1310       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1311       AliceInp << setw(10) << 0.0;
1312       AliceInp << setw(10) << 2.0;
1313       AliceInp << setw(10) << 2.0;
1314       AliceInp << setw(10) << 2.0;
1315       AliceInp << setw(10) << 1.0;
1316       AliceInp << endl;
1317     }
1318
1319     // delta-rays by electrons
1320     // G4 particles: "e-"
1321     // G3 default value: 10**4 GeV
1322     //gMC ->SetCut("DCUTE",cut);  // cut for deltarays by electrons ???????????????
1323     else if (strncmp(&sCutFlag[i][0],"DCUTE",5) == 0) {
1324       AliceInp << "*Cut for deltarays by electrons ????????????"; 
1325       AliceInp << endl;
1326       AliceInp << "*Generated from call: SetCut('DCUTE',cut);"; 
1327       AliceInp << endl;
1328       AliceInp << setw(10) << "EMFCUT    "; 
1329       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1330       AliceInp << setw(10) << -fCutValue[i];
1331       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1332       AliceInp << setw(10) << 0.0;
1333       AliceInp << setw(10) << 0.0;
1334       AliceInp << setw(10) << 2.0;
1335       AliceInp << setprecision(4);
1336       AliceInp << setw(10) << fLastRegion;
1337       AliceInp << setprecision(1);
1338       AliceInp << setw(10) << 1.0;
1339       AliceInp << endl;
1340     }
1341     
1342     // delta-rays by muons
1343     // G4 particles: "e-"
1344     // G3 default value: 10**4 GeV
1345     //gMC ->SetCut("DCUTM",cut);  // cut for deltarays by muons
1346     else if (strncmp(&sCutFlag[i][0],"DCUTM",5) == 0) {
1347       AliceInp << "*Cut for deltarays by muons ????????????"; 
1348       AliceInp << endl;
1349       AliceInp << "*Generated from call: SetCut('DCUTM',cut);"; 
1350       AliceInp << endl;
1351       AliceInp << setw(10) << "DELTARAY  "; 
1352       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1353       AliceInp << setw(10) << fCutValue[i];
1354       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1355       AliceInp << setw(10) << 0.0;
1356       AliceInp << setw(10) << 0.0;
1357       AliceInp << setw(10) << 3.0;
1358       AliceInp << setprecision(2);
1359       AliceInp << setw(10) << fLastMaterial;
1360       AliceInp << setprecision(1);
1361       AliceInp << setw(10) << 1.0;
1362       AliceInp << endl;
1363     }
1364     
1365     // direct pair production by muons
1366     // G4 particles: "e-", "e+"
1367     // G3 default value: 0.01 GeV
1368     //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons ?????????????????????????
1369     else if (strncmp(&sCutFlag[i][0],"PPCUTM",6) == 0) {
1370       AliceInp << "*Total energy cut for direct pair prod. by muons ????????????"; 
1371       AliceInp << endl;
1372       AliceInp << "*Generated from call: SetCut('PPCUTM',cut);"; 
1373       AliceInp << endl;
1374       AliceInp << setw(10) << "PAIRBREM  "; 
1375       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1376       AliceInp << setw(10) << -fCutValue[i];
1377       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1378       AliceInp << setw(10) << 0.0;
1379       AliceInp << setw(10) << 2.0;
1380       AliceInp << setw(10) << 2.0;
1381       AliceInp << setw(10) << 2.0;
1382       AliceInp << setw(10) << 1.0;
1383       AliceInp << endl;
1384     }
1385
1386     // time of flight cut in seconds
1387     // G4 particles: all
1388     // G3 default value: 0.01 GeV
1389     //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
1390     else if (strncmp(&sCutFlag[i][0],"TOFMAX",6) == 0) {
1391       AliceInp << "*Time of flight cuts in seconds"; 
1392       AliceInp << endl;
1393       AliceInp << "*Generated from call: SetCut('TOFMAX',tofmax);"; 
1394       AliceInp << endl;
1395       AliceInp << setw(10) << "TIME-CUT  "; 
1396       AliceInp << setiosflags(ios::scientific) << setprecision(5);
1397       AliceInp << setw(10) << fCutValue[i]*1.e9;
1398       AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1399       AliceInp << setw(10) << 0.0;
1400       AliceInp << setw(10) << 0.0;
1401       AliceInp << setw(10) << -6.0; // lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1402       AliceInp << setprecision(2);
1403       AliceInp << setw(10) << 64.0; // upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1404       AliceInp << setprecision(1);
1405       AliceInp << setw(10) << 1.0; // step length in assigning numbers
1406       AliceInp << endl;
1407     }
1408
1409     else {
1410       cout << "SetCut for flag=" << &sCutFlag[i][0] << " value=" << fCutValue[i] << " not yet implemented!" << endl;
1411     }
1412   } //end of loop over SeCut calls
1413     
1414 // Add START and STOP card
1415    AliceInp << setw(10) << "START     "; 
1416    AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint);
1417    AliceInp << setw(10) << fEventsPerRun;
1418    AliceInp << endl;
1419    AliceInp << setw(10) << "STOP      "; 
1420    AliceInp << endl;
1421
1422 }
1423
1424 //_____________________________________________________________________________
1425 // methods for step management
1426 //____________________________________________________________________________ 
1427 //
1428 // set methods
1429 //
1430 void TFluka::SetMaxStep(Double_t)
1431 {
1432 // SetMaxStep is dummy procedure in TFluka !
1433   if (fVerbosityLevel >=3)
1434   cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
1435 }
1436
1437 void TFluka::SetMaxNStep(Int_t)
1438 {
1439 // SetMaxNStep is dummy procedure in TFluka !
1440   if (fVerbosityLevel >=3)
1441   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1442 }
1443
1444 void TFluka::SetUserDecay(Int_t)
1445 {
1446 // SetUserDecay is dummy procedure in TFluka !
1447   if (fVerbosityLevel >=3)
1448   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1449 }
1450
1451 //
1452 // dynamic properties
1453 //
1454 void TFluka::TrackPosition(TLorentzVector& position) const
1455 {
1456 // Return the current position in the master reference frame of the
1457 // track being transported
1458 // TRACKR.atrack = age of the particle
1459 // TRACKR.xtrack = x-position of the last point
1460 // TRACKR.ytrack = y-position of the last point
1461 // TRACKR.ztrack = z-position of the last point
1462   Int_t caller = GetCaller();
1463   if (caller == 1 || caller == 3 || caller == 6) { //bxdraw,endraw,usdraw
1464     position.SetX(GetXsco());
1465     position.SetY(GetYsco());
1466     position.SetZ(GetZsco());
1467     position.SetT(TRACKR.atrack);
1468   }
1469   else if (caller == 4) { // mgdraw
1470     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1471     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1472     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1473     position.SetT(TRACKR.atrack);
1474   }
1475   else if (caller == 5) { // sodraw
1476     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1477     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1478     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1479     position.SetT(0);
1480   }
1481   else
1482     Warning("TrackPosition","position not available");
1483 }
1484
1485 //
1486 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1487 {
1488 // Return the current position in the master reference frame of the
1489 // track being transported
1490 // TRACKR.atrack = age of the particle
1491 // TRACKR.xtrack = x-position of the last point
1492 // TRACKR.ytrack = y-position of the last point
1493 // TRACKR.ztrack = z-position of the last point
1494   Int_t caller = GetCaller();
1495   if (caller == 1 || caller == 3 || caller == 6) { //bxdraw,endraw,usdraw
1496     x = GetXsco();
1497     y = GetYsco();
1498     z = GetZsco();
1499   }
1500   else if (caller == 4) { // mgdraw
1501     x = TRACKR.xtrack[TRACKR.ntrack];
1502     y = TRACKR.ytrack[TRACKR.ntrack];
1503     z = TRACKR.ztrack[TRACKR.ntrack];
1504   }
1505   else if (caller == 5) { // sodraw
1506     x = TRACKR.xtrack[TRACKR.ntrack];
1507     y = TRACKR.ytrack[TRACKR.ntrack];
1508     z = TRACKR.ztrack[TRACKR.ntrack];
1509   }
1510   else
1511     Warning("TrackPosition","position not available");
1512 }
1513
1514 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1515 {
1516 // Return the direction and the momentum (GeV/c) of the track
1517 // currently being transported
1518 // TRACKR.ptrack = momentum of the particle (not always defined, if
1519 //               < 0 must be obtained from etrack) 
1520 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1521 // TRACKR.etrack = total energy of the particle
1522 // TRACKR.jtrack = identity number of the particle
1523 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1524   Int_t caller = GetCaller();
1525   if (caller != 2) { // not eedraw 
1526     if (TRACKR.ptrack >= 0) {
1527       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1528       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1529       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1530       momentum.SetE(TRACKR.etrack);
1531       return;
1532     }
1533     else {
1534       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1535       momentum.SetPx(p*TRACKR.cxtrck);
1536       momentum.SetPy(p*TRACKR.cytrck);
1537       momentum.SetPz(p*TRACKR.cztrck);
1538       momentum.SetE(TRACKR.etrack);
1539       return;
1540     }
1541   }
1542   else
1543     Warning("TrackMomentum","momentum not available");
1544 }
1545
1546 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1547 {
1548 // Return the direction and the momentum (GeV/c) of the track
1549 // currently being transported
1550 // TRACKR.ptrack = momentum of the particle (not always defined, if
1551 //               < 0 must be obtained from etrack) 
1552 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1553 // TRACKR.etrack = total energy of the particle
1554 // TRACKR.jtrack = identity number of the particle
1555 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1556   Int_t caller = GetCaller();
1557   if (caller != 2) { // not eedraw 
1558     if (TRACKR.ptrack >= 0) {
1559       px = TRACKR.ptrack*TRACKR.cxtrck;
1560       py = TRACKR.ptrack*TRACKR.cytrck;
1561       pz = TRACKR.ptrack*TRACKR.cztrck;
1562       e = TRACKR.etrack;
1563       return;
1564     }
1565     else {
1566       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1567       px = p*TRACKR.cxtrck;
1568       py = p*TRACKR.cytrck;
1569       pz = p*TRACKR.cztrck;
1570       e = TRACKR.etrack;
1571       return;
1572     }
1573   }
1574   else
1575     Warning("TrackMomentum","momentum not available");
1576 }
1577
1578 Double_t TFluka::TrackStep() const
1579 {
1580 // Return the length in centimeters of the current step
1581 // TRACKR.ctrack = total curved path
1582   Int_t caller = GetCaller();
1583   if (caller == 1 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
1584     return 0.0;
1585   else if (caller == 4) //mgdraw
1586     return TRACKR.ctrack;
1587   else
1588     return -1.0;
1589 }
1590
1591 Double_t TFluka::TrackLength() const
1592 {
1593 // Still wrong !!!
1594 // This is the sum of substeps !!!
1595 // TRACKR.ctrack = total curved path of the current step
1596 // Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
1597 // The sum of all step length starting from the beginning of the track
1598 // for the time being returns only the length in centimeters of the current step
1599   Double_t sum = 0;
1600   Int_t caller = GetCaller();
1601   if (caller == 1 || caller == 3 || caller == 4 || caller == 6) { //bxdraw,endraw,mgdraw,usdraw
1602     for ( Int_t j=0;j<TRACKR.ntrack;j++) {
1603       sum +=TRACKR.ttrack[j];
1604     }
1605     return sum;
1606   }
1607   else 
1608     return -1.0;
1609 }
1610
1611 Double_t TFluka::TrackTime() const
1612 {
1613 // Return the current time of flight of the track being transported
1614 // TRACKR.atrack = age of the particle
1615   Int_t caller = GetCaller();
1616   if (caller == 1 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1617     return TRACKR.atrack;
1618   else 
1619     return -1;
1620 }
1621
1622 Double_t TFluka::Edep() const
1623 {
1624 // Energy deposition
1625 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1626 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1627 //    but in the variable "rull" of the procedure "endraw.cxx"
1628 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1629 // -->no energy loss along the track
1630 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1631 // -->energy loss distributed along the track
1632 // TRACKR.dtrack = energy deposition of the jth deposition even
1633   Double_t sum = 0;
1634   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1635     sum +=TRACKR.dtrack[j];  
1636   }
1637   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1638     return fRull + sum;
1639   else {
1640     return sum;
1641   }
1642 }
1643
1644 Int_t TFluka::TrackPid() const
1645 {
1646 // Return the id of the particle transported
1647 // TRACKR.jtrack = identity number of the particle
1648   Int_t caller = GetCaller();
1649   if (caller != 2)  // not eedraw 
1650     return PDGFromId(TRACKR.jtrack);
1651   else
1652     return -1000;
1653 }
1654
1655 Double_t TFluka::TrackCharge() const
1656 {
1657 // Return charge of the track currently transported
1658 // PAPROP.ichrge = electric charge of the particle
1659 // TRACKR.jtrack = identity number of the particle
1660   Int_t caller = GetCaller();
1661   if (caller != 2)  // not eedraw 
1662     return PAPROP.ichrge[TRACKR.jtrack+6];
1663   else
1664     return -1000.0;
1665 }
1666
1667 Double_t TFluka::TrackMass() const
1668 {
1669 // PAPROP.am = particle mass in GeV
1670 // TRACKR.jtrack = identity number of the particle
1671   Int_t caller = GetCaller();
1672   if (caller != 2)  // not eedraw 
1673     return PAPROP.am[TRACKR.jtrack+6];
1674   else
1675     return -1000.0;
1676 }
1677
1678 Double_t TFluka::Etot() const
1679 {
1680 // TRACKR.etrack = total energy of the particle
1681   Int_t caller = GetCaller();
1682   if (caller != 2)  // not eedraw
1683     return TRACKR.etrack;
1684   else
1685     return -1000.0;
1686 }
1687
1688 //
1689 // track status
1690 //
1691 Bool_t   TFluka::IsNewTrack() const
1692 {
1693 // ???????????????,
1694 // True if the track is not at the boundary of the current volume
1695 // Not true in some cases in bxdraw - to be solved
1696   Int_t caller = GetCaller();
1697   if (caller == 1)
1698     return 1; // how to handle double step ?????????????
1699   else
1700     return 0; // ??????????????
1701 }
1702
1703 Bool_t   TFluka::IsTrackInside() const
1704 {
1705 // True if the track is not at the boundary of the current volume
1706 // In Fluka a step is always inside one kind of material
1707 // If the step would go behind the region of one material,
1708 // it will be shortened to reach only the boundary.
1709 // Therefore IsTrackInside() is always true.
1710   Int_t caller = GetCaller();
1711   if (caller == 1)  // bxdraw
1712     return 0;
1713   else
1714     return 1;
1715 }
1716
1717 Bool_t   TFluka::IsTrackEntering() const
1718 {
1719 // True if this is the first step of the track in the current volume
1720   Int_t caller = GetCaller();
1721   if (caller == 11 || caller == 4)  // bxdraw entering
1722     return 1;
1723   else return 0;
1724 }
1725
1726 Bool_t   TFluka::IsTrackExiting() const
1727 {
1728   Int_t caller = GetCaller();
1729   if (caller == 12)  // bxdraw exiting
1730     return 1;
1731   else return 0;
1732 }
1733
1734 Bool_t   TFluka::IsTrackOut() const
1735 {
1736 // True if the track is out of the setup
1737 // means escape
1738 // Icode = 14: escape - call from Kaskad
1739 // Icode = 23: escape - call from Emfsco
1740 // Icode = 32: escape - call from Kasneu
1741 // Icode = 40: escape - call from Kashea
1742 // Icode = 51: escape - call from Kasoph
1743   if (iIcode == 14 ||
1744       iIcode == 23 ||
1745       iIcode == 32 ||
1746       iIcode == 40 ||
1747       iIcode == 51) return 1;
1748   else return 0;
1749 }
1750
1751 Bool_t   TFluka::IsTrackDisappeared() const
1752 {
1753 // means all inelastic interactions and decays
1754 // iIcode from usdraw
1755   if (iIcode == 101 || // inelastic interaction
1756       iIcode == 102 || // particle decay
1757       iIcode == 214 || // in-flight annihilation
1758       iIcode == 215 || // annihilation at rest
1759       iIcode == 217 || // pair production
1760       iIcode == 221) return 1;
1761   else return 0;
1762 }
1763
1764 Bool_t   TFluka::IsTrackStop() const
1765 {
1766 // True if the track energy has fallen below the threshold
1767 // means stopped by signal or below energy threshold
1768 // Icode = 12: stopping particle       - call from Kaskad
1769 // Icode = 15: time kill               - call from Kaskad
1770 // Icode = 21: below threshold, iarg=1 - call from Emfsco
1771 // Icode = 22: below threshold, iarg=2 - call from Emfsco
1772 // Icode = 24: time kill               - call from Emfsco
1773 // Icode = 31: below threshold         - call from Kasneu
1774 // Icode = 33: time kill               - call from Kasneu
1775 // Icode = 41: time kill               - call from Kashea
1776 // Icode = 52: time kill               - call from Kasoph
1777   if (iIcode == 12 ||
1778       iIcode == 15 ||
1779       iIcode == 21 ||
1780       iIcode == 22 ||
1781       iIcode == 24 ||
1782       iIcode == 31 ||
1783       iIcode == 33 ||
1784       iIcode == 41 ||
1785       iIcode == 52) return 1;
1786   else return 0;
1787 }
1788
1789 Bool_t   TFluka::IsTrackAlive() const
1790 {
1791 // means not disappeared or not out
1792   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1793   else return 1;
1794 }
1795
1796 //
1797 // secondaries
1798 //
1799
1800 Int_t TFluka::NSecondaries() const
1801 // Number of secondary particles generated in the current step
1802 // FINUC.np = number of secondaries except light and heavy ions
1803 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
1804 {
1805   Int_t caller = GetCaller();
1806   if (caller == 6)  // valid only after usdraw
1807     return FINUC.np + FHEAVY.npheav;
1808   else
1809     return 0;
1810 } // end of NSecondaries
1811
1812 void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1813                 TLorentzVector& position, TLorentzVector& momentum)
1814 {
1815   Int_t caller = GetCaller();
1816   if (caller == 6) {  // valid only after usdraw
1817     if (isec >= 0 && isec < FINUC.np) {
1818       particleId = PDGFromId(FINUC.kpart[isec]);
1819       position.SetX(fXsco);
1820       position.SetY(fYsco);
1821       position.SetZ(fZsco);
1822       position.SetT(TRACKR.atrack);
1823 //    position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem.
1824       momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
1825       momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
1826       momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
1827       momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
1828     }
1829     else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
1830       Int_t jsec = isec - FINUC.np;
1831       particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1832       position.SetX(fXsco);
1833       position.SetY(fYsco);
1834       position.SetZ(fZsco);
1835       position.SetT(TRACKR.atrack);
1836 //    position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem.
1837       momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1838       momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1839       momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1840       if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
1841         momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1842       else if (FHEAVY.tkheav[jsec] > 6)
1843         momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1844     }
1845     else
1846       Warning("GetSecondary","isec out of range");
1847   }
1848   else
1849     Warning("GetSecondary","no secondaries available");
1850 } // end of GetSecondary
1851
1852 TMCProcess TFluka::ProdProcess(Int_t isec) const
1853 // Name of the process that has produced the secondary particles
1854 // in the current step
1855 {
1856     const TMCProcess kIpNoProc = kPNoProcess;
1857     const TMCProcess kIpPDecay = kPDecay;
1858     const TMCProcess kIpPPair = kPPair;
1859 //  const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
1860 //  const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
1861     const TMCProcess kIpPCompton = kPCompton;
1862     const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
1863     const TMCProcess kIpPBrem = kPBrem;
1864 //  const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
1865 //  const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
1866     const TMCProcess kIpPDeltaRay = kPDeltaRay;
1867 //  const TMCProcess kIpPMoller = kPMoller;
1868 //  const TMCProcess kIpPBhabha = kPBhabha;
1869     const TMCProcess kIpPAnnihilation = kPAnnihilation;
1870 //  const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
1871 //  const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
1872     const TMCProcess kIpPHadronic = kPHadronic;
1873     const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
1874     const TMCProcess kIpPPhotoFission = kPPhotoFission;
1875     const TMCProcess kIpPRayleigh = kPRayleigh;
1876 //  const TMCProcess kIpPCerenkov = kPCerenkov;
1877 //  const TMCProcess kIpPSynchrotron = kPSynchrotron;
1878
1879     Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
1880     if (iIcode == 102) return kIpPDecay;
1881     else if (iIcode == 104 || iIcode == 217) return kIpPPair;
1882 //  else if (iIcode == 104) return kIpPairFromPhoton;
1883 //  else if (iIcode == 217) return kIpPPairFromVirtualPhoton;
1884     else if (iIcode == 219) return kIpPCompton;
1885     else if (iIcode == 221) return kIpPPhotoelectric;
1886     else if (iIcode == 105 || iIcode == 208) return kIpPBrem;
1887 //  else if (iIcode == 105) return kIpPBremFromHeavy;
1888 //  else if (iIcode == 208) return kPBremFromElectronOrPositron;
1889     else if (iIcode == 103 || iIcode == 400) return kIpPDeltaRay;
1890     else if (iIcode == 210 || iIcode == 212) return kIpPDeltaRay;
1891 //  else if (iIcode == 210) return kIpPMoller;
1892 //  else if (iIcode == 212) return kIpPBhabha;
1893     else if (iIcode == 214 || iIcode == 215) return kIpPAnnihilation;
1894 //  else if (iIcode == 214) return kIpPAnnihilInFlight;
1895 //  else if (iIcode == 215) return kIpPAnnihilAtRest;
1896     else if (iIcode == 101) return kIpPHadronic;
1897     else if (iIcode == 101) {
1898       if (!mugamma) return kIpPHadronic;
1899       else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
1900       else return kIpPMuonNuclear;
1901     }
1902     else if (iIcode == 225) return kIpPRayleigh;
1903 // Fluka codes 100, 300 and 400 still to be investigasted
1904     else return kIpNoProc;
1905 }
1906
1907 //Int_t StepProcesses(TArrayI &proc) const
1908 // Return processes active in the current step
1909 //{
1910 //ck = total energy of the particl ???????????????? 
1911 //}
1912
1913
1914 Int_t TFluka::VolId2Mate(Int_t id) const
1915 {
1916 //
1917 // Returns the material number for a given volume ID
1918 //
1919     if (fVerbosityLevel >= 3)
1920     printf("VolId2Mate %d %d\n", id, fMediaByRegion[id]); 
1921     return fMediaByRegion[id-1];
1922 }
1923
1924 const char* TFluka::VolName(Int_t id) const
1925 {
1926 //
1927 // Returns the volume name for a given volume ID
1928 //
1929     FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[id-1]);
1930     const char* name = vol->GetName();
1931     if (fVerbosityLevel >= 3)
1932     printf("VolName %d %s \n", id, name);
1933     return name;
1934 }
1935
1936 Int_t TFluka::VolId(const Text_t* volName) const
1937 {
1938 //
1939 // Converts from volume name to volume ID.
1940 // Time consuming. (Only used during set-up)
1941 // Could be replaced by hash-table
1942 //
1943     char tmp[5];
1944     Int_t i =0;
1945     for (i = 0; i < fNVolumes; i++)
1946   {
1947       FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
1948       TString name = vol->GetName();
1949       strcpy(tmp, name.Data());
1950       tmp[4] = '\0';
1951       if (!strcmp(tmp, volName)) break;
1952   }
1953     i++;
1954
1955     return i;
1956 }
1957
1958
1959 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
1960 {
1961 //
1962 // Return the logical id and copy number corresponding to the current fluka region
1963 //
1964     int ir = fCurrentFlukaRegion;
1965     int id = (FGeometryInit::GetInstance())->CurrentVolID(ir, copyNo);
1966     if (fVerbosityLevel >= 3)
1967     printf("CurrentVolID: %d %d %d \n", ir, id, copyNo); 
1968     return id;
1969
1970
1971
1972 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
1973 {
1974 //
1975 // Return the logical id and copy number of off'th mother 
1976 // corresponding to the current fluka region
1977 //
1978     if (off == 0) 
1979         return CurrentVolID(copyNo);
1980     
1981     int ir = fCurrentFlukaRegion;
1982     int id = (FGeometryInit::GetInstance())->CurrentVolOffID(ir, off, copyNo);
1983     if (fVerbosityLevel >= 3)
1984     printf("CurrentVolOffID: %d %d %d \n", ir, id, copyNo); 
1985     if (id == -1) 
1986         if (fVerbosityLevel >= 0)
1987         printf("CurrentVolOffID: Warning Mother not found !!!\n"); 
1988     return id;
1989 }
1990
1991
1992 const char* TFluka::CurrentVolName() const
1993 {
1994 //
1995 // Return the current volume name
1996 //
1997     Int_t copy;
1998     Int_t id = TFluka::CurrentVolID(copy);
1999     const char* name = TFluka::VolName(id);
2000     if (fVerbosityLevel >= 3)
2001     printf("CurrentVolumeName: %d %s \n", fCurrentFlukaRegion,  name); 
2002     return name;
2003 }
2004
2005 const char* TFluka::CurrentVolOffName(Int_t off) const
2006 {
2007 //
2008 // Return the volume name of the off'th mother of the current volume
2009 //
2010     Int_t copy;
2011     Int_t id = TFluka::CurrentVolOffID(off, copy);
2012     const char* name = TFluka::VolName(id);
2013     if (fVerbosityLevel >= 3)
2014     printf("CurrentVolumeOffName: %d %s \n", fCurrentFlukaRegion,  name); 
2015     return name;
2016 }
2017
2018 Int_t TFluka::CurrentMaterial(Float_t &a, Float_t &z, 
2019                       Float_t &dens, Float_t &radl, Float_t &absl) const
2020 {
2021 //
2022 //  Return the current medium number
2023 //
2024     Int_t copy;
2025     Int_t id  =  TFluka::CurrentVolID(copy);
2026     Int_t med =  TFluka::VolId2Mate(id);
2027     if (fVerbosityLevel >= 3)
2028     printf("CurrentMaterial: %d %d \n", fCurrentFlukaRegion,  med); 
2029     return med;
2030 }
2031
2032 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2033     {
2034 // Transforms a position from the world reference frame
2035 // to the current volume reference frame.
2036 //
2037 //  Geant3 desription:
2038 //  ==================
2039 //       Computes coordinates XD (in DRS) 
2040 //       from known coordinates XM in MRS 
2041 //       The local reference system can be initialized by
2042 //         - the tracking routines and GMTOD used in GUSTEP
2043 //         - a call to GMEDIA(XM,NUMED)
2044 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) 
2045 //             (inverse routine is GDTOM) 
2046 //
2047 //        If IFLAG=1  convert coordinates 
2048 //           IFLAG=2  convert direction cosinus
2049 //
2050 // ---
2051         Double_t xmD[3], xdD[3];        
2052         xmD[0] = xm[0]; xmD[1] = xm[1]; xmD[2] = xm[2]; 
2053         (FGeometryInit::GetInstance())->Gmtod(xmD, xdD, iflag);
2054         xd[0] = xdD[0]; xd[1] = xdD[1]; xd[2] = xdD[2]; 
2055     }
2056
2057   
2058 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2059     {
2060 // Transforms a position from the world reference frame
2061 // to the current volume reference frame.
2062 //
2063 //  Geant3 desription:
2064 //  ==================
2065 //       Computes coordinates XD (in DRS) 
2066 //       from known coordinates XM in MRS 
2067 //       The local reference system can be initialized by
2068 //         - the tracking routines and GMTOD used in GUSTEP
2069 //         - a call to GMEDIA(XM,NUMED)
2070 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) 
2071 //             (inverse routine is GDTOM) 
2072 //
2073 //        If IFLAG=1  convert coordinates 
2074 //           IFLAG=2  convert direction cosinus
2075 //
2076 // ---
2077         Double_t xmD[3], xdD[3];        
2078         xdD[0] = xd[0]; xdD[1] = xd[1]; xdD[2] = xd[2]; 
2079         (FGeometryInit::GetInstance())->Gdtom(xmD, xdD, iflag);
2080         xm[0] = xmD[0]; xm[1] = xmD[1]; xm[2] = xmD[2]; 
2081     }
2082
2083 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2084     {
2085 // Transforms a position from the current volume reference frame
2086 // to the world reference frame.
2087 //
2088 //  Geant3 desription:
2089 //  ==================
2090 //  Computes coordinates XM (Master Reference System
2091 //  knowing the coordinates XD (Detector Ref System)
2092 //  The local reference system can be initialized by
2093 //    - the tracking routines and GDTOM used in GUSTEP
2094 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2095 //        (inverse routine is GMTOD)
2096 // 
2097 //   If IFLAG=1  convert coordinates
2098 //      IFLAG=2  convert direction cosinus
2099 //
2100 // ---
2101
2102
2103     }
2104 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2105     {
2106 // Transforms a position from the current volume reference frame
2107 // to the world reference frame.
2108 //
2109 //  Geant3 desription:
2110 //  ==================
2111 //  Computes coordinates XM (Master Reference System
2112 //  knowing the coordinates XD (Detector Ref System)
2113 //  The local reference system can be initialized by
2114 //    - the tracking routines and GDTOM used in GUSTEP
2115 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2116 //        (inverse routine is GMTOD)
2117 // 
2118 //   If IFLAG=1  convert coordinates
2119 //      IFLAG=2  convert direction cosinus
2120 //
2121 // ---
2122
2123         (FGeometryInit::GetInstance())->Gdtom(xm, xd, iflag);
2124     }
2125
2126 // ===============================================================
2127 void TFluka::FutoTest() 
2128 {
2129     Int_t icode, mreg, newreg, particleId;
2130     Double_t rull, xsco, ysco, zsco;
2131     TLorentzVector position, momentum;
2132     icode = GetIcode();
2133     if (icode == 0) {
2134         if (fVerbosityLevel >=3)
2135             cout << " icode=" << icode << endl;
2136     } else if (icode > 0 && icode <= 5) {
2137 // mgdraw
2138         mreg = GetMreg();
2139         if (fVerbosityLevel >=3)
2140             cout << " icode=" << icode
2141                  << " mreg=" << mreg
2142                  << endl;
2143         TrackPosition(position);
2144         TrackMomentum(momentum);
2145         if (fVerbosityLevel >=3) {
2146             cout << "TLorentzVector positionX=" << position.X()
2147                  << "positionY=" << position.Y()
2148                  << "positionZ=" << position.Z()
2149                  << "timeT=" << position.T() << endl;
2150             cout << "TLorentzVector momentumX=" << momentum.X()
2151                  << "momentumY=" << momentum.Y()
2152                  << "momentumZ=" << momentum.Z()
2153                  << "energyE=" << momentum.E() << endl;
2154             cout << "TrackStep=" << TrackStep() << endl;
2155             cout << "TrackLength=" << TrackLength() << endl;
2156             cout << "TrackTime=" << TrackTime() << endl;
2157             cout << "Edep=" << Edep() << endl;
2158             cout << "TrackPid=" << TrackPid() << endl;
2159             cout << "TrackCharge=" << TrackCharge() << endl;
2160             cout << "TrackMass=" << TrackMass() << endl;
2161             cout << "Etot=" << Etot() << endl;
2162             cout << "IsNewTrack=" << IsNewTrack() << endl;
2163             cout << "IsTrackInside=" << IsTrackInside() << endl;
2164             cout << "IsTrackEntering=" << IsTrackEntering() << endl;
2165             cout << "IsTrackExiting=" << IsTrackExiting() << endl;
2166             cout << "IsTrackOut=" << IsTrackOut() << endl;
2167             cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
2168             cout << "IsTrackAlive=" << IsTrackAlive() << endl;
2169         }
2170         
2171         Float_t x = position.X();
2172         Float_t y = position.Y();
2173         Float_t z = position.Z();
2174         Float_t xm[3];
2175         Float_t xd[3];
2176         xm[0] = x; xm[1] = y; xm[2] = z;
2177         if (fVerbosityLevel >= 3)
2178             printf("Global trackPosition: %f %f %f \n", x, y, z);
2179         Gmtod(xm, xd, 1);
2180         if (fVerbosityLevel >= 3)
2181             printf("Local trackPosition: %f %f %f \n", xd[0], xd[1], xd[2]);
2182         Gdtom(xd, xm, 1);
2183         if (fVerbosityLevel >= 3)
2184             printf("New trackPosition: %f %f %f \n", xm[0], xm[1], xm[2]);
2185     } else if((icode >= 10 && icode <= 15) ||
2186               (icode >= 20 && icode <= 24) ||
2187               (icode >= 30 && icode <= 33) ||
2188               (icode >= 40 && icode <= 41) ||
2189               (icode >= 50 && icode <= 52)) {
2190 // endraw
2191         mreg = GetMreg();
2192         rull = GetRull();
2193         xsco = GetXsco();
2194         ysco = GetYsco();
2195         zsco = GetZsco();
2196         
2197         if (fVerbosityLevel >=3) {     
2198             cout << " icode=" << icode
2199                  << " mreg=" << mreg
2200                  << " rull=" << rull
2201                  << " xsco=" << xsco
2202                  << " ysco=" << ysco
2203                  << " zsco=" << zsco << endl;
2204         }
2205         TrackPosition(position);
2206         TrackMomentum(momentum);
2207         if (fVerbosityLevel >=3) {
2208             cout << "Edep=" << Edep() << endl;
2209             cout << "Etot=" << Etot() << endl;
2210             cout << "TrackPid=" << TrackPid() << endl;
2211             cout << "TrackCharge=" << TrackCharge() << endl;
2212             cout << "TrackMass=" << TrackMass() << endl;
2213             cout << "IsTrackOut=" << IsTrackOut() << endl;
2214             cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
2215             cout << "IsTrackStop=" << IsTrackStop() << endl;
2216             cout << "IsTrackAlive=" << IsTrackAlive() << endl;
2217         }
2218     } else if((icode >= 100 && icode <= 105) ||
2219               (icode == 208) ||
2220               (icode == 210) ||
2221               (icode == 212) ||
2222               (icode >= 214 && icode <= 215) ||
2223               (icode == 217) ||
2224               (icode == 219) ||
2225               (icode == 221) ||
2226               (icode == 225) ||
2227               (icode == 300) ||
2228               (icode == 400)) {
2229 // usdraw
2230           mreg = GetMreg();
2231           xsco = GetXsco();
2232           ysco = GetYsco();
2233           zsco = GetZsco();
2234           
2235           if (fVerbosityLevel >=3) {
2236               cout << " icode=" << icode
2237                    << " mreg=" << mreg
2238                    << " xsco=" << xsco
2239                    << " ysco=" << ysco
2240                    << " zsco=" << zsco << endl;
2241               cout << "TrackPid=" << TrackPid() << endl;
2242               cout << "NSecondaries=" << NSecondaries() << endl;
2243           }
2244           
2245           for (Int_t isec=0; isec< NSecondaries(); isec++) {
2246               TFluka::GetSecondary(isec, particleId, position, momentum);
2247               if (fVerbosityLevel >=3) {
2248                   cout << "TLorentzVector positionX=" << position.X()
2249                        << "positionY=" << position.Y()
2250                        << "positionZ=" << position.Z()
2251                        << "timeT=" << position.T() << endl;
2252                   cout << "TLorentzVector momentumX=" << momentum.X()
2253                        << "momentumY=" << momentum.Y()
2254                        << "momentumZ=" << momentum.Z()
2255                        << "energyE=" << momentum.E() << endl;
2256                   cout << "TrackPid=" << particleId << endl;
2257               }
2258           }
2259     } else if((icode == 19) ||
2260               (icode == 29) ||
2261               (icode == 39) ||
2262               (icode == 49) ||
2263               (icode == 59)) {
2264         mreg = GetMreg();
2265         newreg = GetNewreg();
2266         xsco = GetXsco();
2267         ysco = GetYsco();
2268         zsco = GetZsco();
2269         if (fVerbosityLevel >=3) {
2270             cout << " icode=" << icode
2271                  << " mreg=" << mreg
2272                  << " newreg=" << newreg
2273                  << " xsco=" << xsco
2274                  << " ysco=" << ysco
2275                  << " zsco=" << zsco << endl;
2276         }
2277     }
2278 } // end of FutoTest
2279