]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TGeant3/TGeant3.cxx
No instance of AliDecayer created to initialize fDecayer.
[u/mrichter/AliRoot.git] / TGeant3 / TGeant3.cxx
index ba8f55f0e6b3046240320c18003914ffcfc1bdfc..498ee8c44b197892c627049c6a9a4cde41ddcf22 100644 (file)
@@ -1,3 +1,91 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+Revision 1.34  2000/09/07 12:12:01  morsch
+Comment inside comment removed.
+
+Revision 1.33  2000/09/06 16:03:42  morsch
+Set ExternalDecayer, Decayer  and SetForceDecay methods added.
+Gspart calls for charmed and bottom hadrons added.
+Decay mode definitions for charmed and beauty hadrons have been taken out.
+They will be  handled by an external decayer.
+
+Revision 1.32  2000/08/24 16:28:53  hristov
+TGeant3::IsNewTrack corrected by F.Carminati
+
+Revision 1.31  2000/07/13 16:19:10  fca
+Mainly coding conventions + some small bug fixes
+
+Revision 1.30  2000/07/12 08:56:30  fca
+Coding convention correction and warning removal
+
+Revision 1.29  2000/07/11 18:24:59  fca
+Coding convention corrections + few minor bug fixes
+
+Revision 1.28  2000/06/29 10:51:55  morsch
+Add some charmed and bottom baryons to the particle list (TDatabasePDG). This
+is needed by Hijing. Should be part of a future review of TDatabasePDG.
+
+Revision 1.27  2000/06/21 17:40:15  fca
+Adding possibility to set ISTRA, PAI model
+
+Revision 1.26  2000/05/16 13:10:41  fca
+New method IsNewTrack and fix for a problem in Father-Daughter relations
+
+Revision 1.25  2000/04/07 11:12:35  fca
+G4 compatibility changes
+
+Revision 1.24  2000/02/28 21:03:57  fca
+Some additions to improve the compatibility with G4
+
+Revision 1.23  2000/02/23 16:25:25  fca
+AliVMC and AliGeant3 classes introduced
+ReadEuclid moved from AliRun to AliModule
+
+Revision 1.22  2000/01/18 15:40:13  morsch
+Interface to GEANT3 routines GFTMAT, GBRELM and GPRELM added
+Define geant particle type 51: Feedback Photon with Cherenkov photon properties.
+
+Revision 1.21  2000/01/17 19:41:17  fca
+Add SetERAN function
+
+Revision 1.20  2000/01/12 11:29:27  fca
+Close material file
+
+Revision 1.19  1999/12/17 09:03:12  fca
+Introduce a names array
+
+Revision 1.18  1999/11/26 16:55:39  fca
+Reimplement CurrentVolName() to avoid memory leaks
+
+Revision 1.17  1999/11/03 16:58:28  fca
+Correct source of address violation in creating character string
+
+Revision 1.16  1999/11/03 13:17:08  fca
+Have ProdProcess return const char*
+
+Revision 1.15  1999/10/26 06:04:50  fca
+Introduce TLorentzVector in AliMC::GetSecondary. Thanks to I.Hrivnacova
+
+Revision 1.14  1999/09/29 09:24:30  fca
+Introduction of the Copyright and cvs Log
+
+*/
+
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
 //  Interface Class to the Geant3.21 MonteCarlo                              //
 #include "ctype.h" 
 #include <TDatabasePDG.h>
 #include "AliCallf77.h" 
+#include "AliDecayer.h" 
  
 #ifndef WIN32 
 # define gzebra  gzebra_ 
 # define gfmate  gfmate_ 
 # define gfpart  gfpart_ 
 # define gftmed  gftmed_ 
+# define gftmat  gftmat_ 
 # define gmate   gmate_ 
 # define gpart   gpart_ 
 # define gsdk    gsdk_ 
 # define gsxyz   gsxyz_ 
 # define gtrack  gtrack_ 
 # define gtreve  gtreve_ 
-# define gtreve_root  gtreve_root_ 
+# define gtreveroot  gtreveroot_ 
 # define grndm   grndm_ 
 # define grndmq  grndmq_ 
 # define gdtom   gdtom_ 
 # define gckmat  gckmat_
 # define geditv  geditv_
 # define mzdrop  mzdrop_
+
+# define ertrak  ertrak_
+# define ertrgo  ertrgo_
  
 # define setbomb setbomb_
 # define setclip setclip_
 # define gcomad gcomad_
 
+# define gbrelm gbrelm_
+# define gprelm gprelm_
 #else 
 # define gzebra  GZEBRA 
 # define grfile  GRFILE 
 # define gfmate  GFMATE 
 # define gfpart  GFPART 
 # define gftmed  GFTMED 
+# define gftmat  GFTMAT
 # define gmate   GMATE 
 # define gpart   GPART 
 # define gsdk    GSDK 
 # define gsxyz   GSXYZ 
 # define gtrack  GTRACK 
 # define gtreve  GTREVE 
-# define gtreve_root  GTREVE_ROOT
+# define gtreveroot  GTREVEROOT
 # define grndm   GRNDM
 # define grndmq  GRNDMQ
 # define gdtom   GDTOM 
 # define geditv  GEDITV
 # define mzdrop  MZDROP 
 
+# define ertrak  ERTRAK
+# define ertrgo  ERTRGO
 # define setbomb SETBOMB
 # define setclip SETCLIP
 # define gcomad  GCOMAD
  
+# define gbrelm GBRELM
+# define gprelm GPRELM
+
 #endif 
 
 //____________________________________________________________________________ 
@@ -271,7 +373,7 @@ extern "C"
 
   void type_of_call gtreve(); 
 
-  void type_of_call gtreve_root(); 
+  void type_of_call gtreveroot(); 
 
   void type_of_call grndm(Float_t *, const Int_t &); 
 
@@ -307,6 +409,10 @@ extern "C"
                           Float_t &, Float_t &, Float_t &, Float_t &,
                           Float_t &, Float_t &, Float_t *, Int_t * DEFCHARL); 
 
+  void type_of_call gftmat(const Int_t&, const Int_t&, DEFCHARD, const Int_t&,
+                          Float_t*, Float_t*
+                          ,Float_t *, Int_t & DEFCHARL);
+
   void type_of_call gsmate(const Int_t&, DEFCHARD, Float_t &, Float_t &,
                           Float_t &, Float_t &, Float_t &, Float_t *,
                           Int_t & DEFCHARL); 
@@ -409,12 +515,21 @@ extern "C"
   void type_of_call setclip(DEFCHARD, Float_t &,Float_t &,Float_t &,Float_t &,
                            Float_t &, Float_t & DEFCHARL); 
   void type_of_call gcomad(DEFCHARD, Int_t*& DEFCHARL); 
+
+  void type_of_call ertrak(const Float_t *const x1, const Float_t *const p1,
+                          const Float_t *x2, const Float_t *p2,
+                          const Int_t &ipa, DEFCHARD DEFCHARL);
+
+  void type_of_call ertrgo();
+  
+    float type_of_call gbrelm(const Float_t &z, const Float_t& t, const Float_t& cut);
+    float type_of_call gprelm(const Float_t &z, const Float_t& t, const Float_t& cut);
 }
 
 //
 // Geant3 global pointer
 //
-static Int_t defSize = 600;
+static const Int_t kDefSize = 600;
 
 ClassImp(TGeant3) 
  
@@ -447,6 +562,7 @@ TGeant3::TGeant3(const char *title, Int_t nwgeant)
   //
   // Zero number of particles
   fNPDGCodes=0;
+  fDecayer=0;
 } 
 
 //____________________________________________________________________________ 
@@ -470,11 +586,11 @@ void TGeant3::DefaultRange()
   //
   // Set range of current drawing pad to 20x20 cm
   //
-  if (!higz) {
-    new THIGZ(defSize); 
+  if (!gHigz) {
+    new THIGZ(kDefSize); 
     gdinit();
   }
-  higz->Range(0,0,20,20);
+  gHigz->Range(0,0,20,20);
 }
 
 //____________________________________________________________________________ 
@@ -483,8 +599,8 @@ void TGeant3::InitHIGZ()
   //
   // Initialise HIGZ
   //
-  if (!higz) {
-    new THIGZ(defSize); 
+  if (!gHigz) {
+    new THIGZ(kDefSize); 
     gdinit();
   }
 }
@@ -501,6 +617,7 @@ void TGeant3::LoadAddress()
   gcomad(PASSCHARD("GCBANK"),(int*&) fGcbank  PASSCHARL("GCBANK"));
   gcomad(PASSCHARD("GCLINK"),(int*&) fGclink  PASSCHARL("GCLINK"));
   gcomad(PASSCHARD("GCCUTS"),(int*&) fGccuts  PASSCHARL("GCCUTS"));
+  gcomad(PASSCHARD("GCMULO"),(int*&) fGcmulo  PASSCHARL("GCMULO"));
   gcomad(PASSCHARD("GCFLAG"),(int*&) fGcflag  PASSCHARL("GCFLAG"));
   gcomad(PASSCHARD("GCKINE"),(int*&) fGckine  PASSCHARL("GCKINE"));
   gcomad(PASSCHARD("GCKING"),(int*&) fGcking  PASSCHARL("GCKING"));
@@ -514,10 +631,18 @@ void TGeant3::LoadAddress()
   gcomad(PASSCHARD("GCNUM"), (int*&) fGcnum   PASSCHARL("GCNUM"));
   gcomad(PASSCHARD("GCSETS"),(int*&) fGcsets  PASSCHARL("GCSETS"));
   gcomad(PASSCHARD("GCPHYS"),(int*&) fGcphys  PASSCHARL("GCPHYS"));
+  gcomad(PASSCHARD("GCPHLT"),(int*&) fGcphlt  PASSCHARL("GCPHLT"));
   gcomad(PASSCHARD("GCOPTI"),(int*&) fGcopti  PASSCHARL("GCOPTI"));
   gcomad(PASSCHARD("GCTLIT"),(int*&) fGctlit  PASSCHARL("GCTLIT"));
   gcomad(PASSCHARD("GCVDMA"),(int*&) fGcvdma  PASSCHARL("GCVDMA"));
 
+  // Commons for GEANE
+  gcomad(PASSCHARD("ERTRIO"),(int*&) fErtrio  PASSCHARL("ERTRIO"));
+  gcomad(PASSCHARD("EROPTS"),(int*&) fEropts  PASSCHARL("EROPTS"));
+  gcomad(PASSCHARD("EROPTC"),(int*&) fEroptc  PASSCHARL("EROPTC"));
+  gcomad(PASSCHARD("ERWORK"),(int*&) fErwork  PASSCHARL("ERWORK"));
+
+  // Variables for ZEBRA store
   gcomad(PASSCHARD("IQ"), addr  PASSCHARL("IQ"));
   fZiq = addr;
   gcomad(PASSCHARD("LQ"), addr  PASSCHARL("LQ"));
@@ -535,6 +660,13 @@ void TGeant3::GeomIter()
   fNextVol=fGcvolu->nlevel;
 }
 
+//____________________________________________________________________________ 
+void TGeant3::FinishGeometry()
+{
+  //Close the geometry structure
+  Ggclos();
+}
+  
 //____________________________________________________________________________ 
 Int_t TGeant3::NextVolUp(Text_t *name, Int_t &copy)
 {
@@ -546,16 +678,21 @@ Int_t TGeant3::NextVolUp(Text_t *name, Int_t &copy)
   fNextVol--;
   if(fNextVol>=0) {
     gname=fGcvolu->names[fNextVol];
-    strncpy(name,(char *) &gname, 4);
-    name[4]='\0';
     copy=fGcvolu->number[fNextVol];
     i=fGcvolu->lvolum[fNextVol];
+    name = fVolNames[i-1];
     if(gname == fZiq[fGclink->jvolum+i]) return i;
     else printf("GeomTree: Volume %s not found in bank\n",name);
   }
   return 0;
 }
 
+//_____________________________________________________________________________
+void TGeant3::BuildPhysics()
+{
+  Gphysi();
+}
+
 //_____________________________________________________________________________
 Int_t TGeant3::CurrentVolID(Int_t &copy) const
 {
@@ -603,17 +740,13 @@ const char* TGeant3::CurrentVolName() const
   // Returns the current volume name
   //
   Int_t i, gname;
-  char *name;
   if( (i=fGcvolu->nlevel-1) < 0 ) {
     Warning("CurrentVolName","Stack depth %d\n",fGcvolu->nlevel);
   } else {
     gname=fGcvolu->names[i];
-    name = new char[5];
-    strncpy(name,(char *) &gname, 4);
-    name[4]='\0';
     i=fGcvolu->lvolum[i];   
-    if(gname == fZiq[fGclink->jvolum+i]) return name;
-    else Warning("CurrentVolName","Volume %4s not found\n",name);
+    if(gname == fZiq[fGclink->jvolum+i]) return fVolNames[i-1];
+    else Warning("CurrentVolName","Volume %4s not found\n",(char*) &gname);
   }
   return 0;
 }
@@ -627,18 +760,14 @@ const char* TGeant3::CurrentVolOffName(Int_t off) const
   // if name=0 no name is returned
   //
   Int_t i, gname;
-  char *name;
   if( (i=fGcvolu->nlevel-off-1) < 0 ) {
     Warning("CurrentVolOffName",
            "Offset requested %d but stack depth %d\n",off,fGcvolu->nlevel);
   } else {
     gname=fGcvolu->names[i];
-    name = new char[5];
-    strncpy(name,(char *) &gname, 4);
-    name[4]='\0';
     i=fGcvolu->lvolum[i];    
-    if(gname == fZiq[fGclink->jvolum+i]) return name;
-    else Warning("CurrentVolOffName","Volume %4s not found\n",name);
+    if(gname == fZiq[fGclink->jvolum+i]) return fVolNames[i-1];
+    else Warning("CurrentVolOffName","Volume %4s not found\n",(char*)&gname);
   }
   return 0;
 }
@@ -648,7 +777,7 @@ Int_t TGeant3::IdFromPDG(Int_t pdg) const
 {
   //
   // Return Geant3 code from PDG and pseudo ENDF code
-
+  //
   for(Int_t i=0;i<fNPDGCodes;++i)
     if(pdg==fPDGCode[i]) return i;
   return -1;
@@ -657,6 +786,9 @@ Int_t TGeant3::IdFromPDG(Int_t pdg) const
 //_____________________________________________________________________________
 Int_t TGeant3::PDGFromId(Int_t id) const 
 {
+  //
+  // Return PDG code and pseudo ENDF code from Geant3 code
+  //
   if(id>0 && id<fNPDGCodes) return fPDGCode[id];
   else return -1;
 }
@@ -737,10 +869,10 @@ void TGeant3::DefineParticles()
   fPDGCode[fNPDGCodes++]=213;   // 42 = RHO+
 
   Gspart(43, "RHO -", 4, 0.768, -1., 4.353e-24);
-  fPDGCode[fNPDGCodes++]=-213;   // 40 = RHO-
+  fPDGCode[fNPDGCodes++]=-213;   // 43 = RHO-
 
   Gspart(44, "RHO 0", 3, 0.768, 0., 4.353e-24);
-  fPDGCode[fNPDGCodes++]=113;   // 37 = D0
+  fPDGCode[fNPDGCodes++]=113;   //  44 = RHO0
 
   //
   // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
@@ -754,28 +886,93 @@ void TGeant3::DefineParticles()
 
   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
 
-  const Double_t autogev=0.9314943228;
-  const Double_t hslash = 1.0545726663e-27;
-  const Double_t erggev = 1/1.6021773349e-3;
-  const Double_t hshgev = hslash*erggev;
-  const Double_t yearstosec = 3600*24*365.25;
+  const Double_t kAu2Gev=0.9314943228;
+  const Double_t khSlash = 1.0545726663e-27;
+  const Double_t kErg2Gev = 1/1.6021773349e-3;
+  const Double_t khShGev = khSlash*kErg2Gev;
+  const Double_t kYear2Sec = 3600*24*365.25;
+//
+// Bottom mesons
+// mass and life-time from PDG
+  pdgDB->AddParticle("B(s)*0","B(s)*0",
+                    5.4163, kTRUE, 0.047, +0.,"Meson",  533);
+
+  pdgDB->AddParticle("B(s)*0 bar","B(s)*0 bar",
+                    5.4163, kTRUE, 0.047, -0.,"Meson", -533);
+
+// Charmed baryons
+// 
+// value for mass used by Hijing
+  pdgDB->AddParticle("Sigma(c)*+","Sigma(c)*+",
+                    2.4536, kTRUE, -1., +1.,"Baryon",  4214);
+
+  pdgDB->AddParticle("Sigma(c)*-","Sigma(c)*-",
+                    2.4536, kTRUE, -1., -1.,"Baryon", -4214);
+// equivalent to 4312 ? Hijing uses m=2.55
+  pdgDB->AddParticle("Xsi(c)0","Xsi(c)0",
+                    2.4703, kTRUE, -1., +0.,"Baryon",  4132);
+
+  pdgDB->AddParticle("Xsi(c)0 bar","Xsi(c)0 bar",
+                    2.4703, kTRUE, -1., -0.,"Baryon", -4132);
+// equivalent to 4322 ? Hijing uses m=2.55
+  pdgDB->AddParticle("Xi(c)+","Xi(c)+",
+                    2.4656, kFALSE, -1., +1.,"Baryon",  4232);
+  
+  pdgDB->AddParticle("Xi(c)-","Xi(c)-",
+                    2.4656, kFALSE, -1., -1.,"Baryon", -4232);
+// mass values from Hijing
+
+  pdgDB->AddParticle("Xsi(c)*0","Xsi(c)*0",
+                    2.63, kTRUE, -1., +0.,"Baryon",  4314);
 
+  pdgDB->AddParticle("Xsi(c)*0 bar","Xsi(c)*0 bar",
+                    2.63, kTRUE, -1., -0.,"Baryon", -4314);
 
-  pdgDB->AddParticle("Deuteron","Deuteron",2*autogev+8.071e-3,kTRUE,
+  pdgDB->AddParticle("Xsi(c)*+","Xsi(c)*+",
+                    2.63, kTRUE, -1., +1.,"Baryon",  4324);
+
+  pdgDB->AddParticle("Xsi(c)*-","Xsi(c)*-",
+                    2.63, kTRUE, -1., -1.,"Baryon", -4324);
+
+// pdg mass value, Hijing uses m=2.73.
+  pdgDB->AddParticle("Omega(c)0","Omega(c)0",
+                    2.7040, kFALSE, khShGev/0.064e-12, +0.,"Baryon",  4332);
+  
+  pdgDB->AddParticle("Omega(c)0 bar","Omega(c)0 bar",
+                    2.7040, kFALSE, khShGev/0.064e-12, -0.,"Baryon", -4332);
+// mass value from Hijing
+  pdgDB->AddParticle("Omega(c)*0","Omega(c)*0",
+                    2.8000, kFALSE, -1., +0.,"Baryon",  4334);
+  
+  pdgDB->AddParticle("Omega(c)*0 bar","Omega(c)*0",
+                    2.8000, kFALSE, -1., -0.,"Baryon", -4334);
+
+// Bottom baryons
+//
+// mass value from Hijing
+  pdgDB->AddParticle("Sigma(b)*+","Sigma(b)*+",
+                    5.8100, kFALSE, -1., +1.,"Baryon", 5224);
+
+  pdgDB->AddParticle("Sigma(b)*-","Sigma(b)*-",
+                    5.8100, kFALSE, -1., -1.,"Baryon", -5224);
+//
+//
+  pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
                     0,1,"Ion",kion+10020);
   fPDGCode[fNPDGCodes++]=kion+10020;   // 45 = Deuteron
 
-  pdgDB->AddParticle("Triton","Triton",3*autogev+14.931e-3,kFALSE,
-                    hshgev/(12.33*yearstosec),1,"Ion",kion+10030);
+  pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
+                    khShGev/(12.33*kYear2Sec),1,"Ion",kion+10030);
   fPDGCode[fNPDGCodes++]=kion+10030;   // 46 = Triton
 
-  pdgDB->AddParticle("Alpha","Alpha",4*autogev+2.424e-3,kTRUE,
-                    hshgev/(12.33*yearstosec),2,"Ion",kion+20040);
+  pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
+                    khShGev/(12.33*kYear2Sec),2,"Ion",kion+20040);
   fPDGCode[fNPDGCodes++]=kion+20040;   // 47 = Alpha
 
   fPDGCode[fNPDGCodes++]=0;   // 48 = geantino mapped to rootino
 
-  pdgDB->AddParticle("HE3","HE3",3*autogev+14.931e-3,kFALSE,
+  pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
                     0,2,"Ion",kion+20030);
   fPDGCode[fNPDGCodes++]=kion+20030;   // 49 = HE3
 
@@ -783,6 +980,67 @@ void TGeant3::DefineParticles()
                     0,0,"Special",kspe+50);
   fPDGCode[fNPDGCodes++]=kspe+50;   // 50 = Cherenkov
 
+  Gspart(51, "FeedbackPhoton", 7, 0., 0.,1.e20 );
+  pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
+                    0,0,"Special",kspe+51);
+  fPDGCode[fNPDGCodes++]=kspe+51;   // 51 = FeedbackPhoton
+  Gspart(52, "Lambda_c+", 4, 2.2849, +1., 2.06e-13);
+  fPDGCode[fNPDGCodes++]=4122;   //52 = Lambda_c+
+
+  Gspart(53, "Lambda_c-", 4, 2.2849, -1., 2.06e-13);
+  fPDGCode[fNPDGCodes++]=-4122; //53 = Lambda_c-  
+
+  Gspart(54, "D_s+", 4, 1.9685, +1., 4.67e-13);
+  fPDGCode[fNPDGCodes++]=431;   //54 = D_s+
+
+  Gspart(55, "D_s-", 4, 1.9685, -1., 4.67e-13);
+  fPDGCode[fNPDGCodes++]=-431; //55 = D_s-
+
+  Gspart(56, "Tau+", 5, 1.77705, +1., 2.9e-13);
+  fPDGCode[fNPDGCodes++]=15;   //56 = Tau+
+
+  Gspart(57, "Tau-", 5, 1.77705, -1., 2.9e-13);
+  fPDGCode[fNPDGCodes++]=-15;  //57 = Tau-  
+
+  Gspart(58, "B0",     3, 5.2792, +0., 1.56e-12);
+  fPDGCode[fNPDGCodes++]=511;   //58 = B0
+
+  Gspart(59, "B0 bar", 3, 5.2792, -0., 1.56e-12);
+  fPDGCode[fNPDGCodes++]=-511;  //58 = B0bar
+
+  Gspart(60, "B+",     4, 5.2789, +1., 1.65e-12);
+  fPDGCode[fNPDGCodes++]=521;   //60 = B+
+
+  Gspart(61, "B-",     4, 5.2789, -1., 1.65e-12);
+  fPDGCode[fNPDGCodes++]=-521;  //61 = B-
+
+  Gspart(62, "Bs",     3, 5.3693, +0., 1.54e-12);
+  fPDGCode[fNPDGCodes++]=521;   //62 = B_s
+
+  Gspart(63, "Bs bar", 3, 5.3693, -0., 1.54e-12);
+  fPDGCode[fNPDGCodes++]=-521;  //63 = B_s bar
+
+  Gspart(64, "Lambda_b",     3, 5.624, +0., 1.24e-12);
+  fPDGCode[fNPDGCodes++]=5122;   //64 = Lambda_b
+
+  Gspart(65, "Lambda_b bar", 3, 5.624, -0., 1.24e-12);
+  fPDGCode[fNPDGCodes++]=-5122;  //65 = Lambda_b bar
+
+  Gspart(66, "J/Psi", 3.09688, 3, 0., 0.);
+  fPDGCode[fNPDGCodes++]=443;       // 66 = J/Psi
+
+  Gspart(67, "Psi Prime", 3, 3.686, 0., 0.);
+  fPDGCode[fNPDGCodes++]=20443;    // 67 = Psi prime
+
+  Gspart(68, "Upsilon(1S)", 9.46037, 3, 0., 0.);
+  fPDGCode[fNPDGCodes++]=553;       // 68 = Upsilon(1S)
+
+  Gspart(69, "Upsilon(2S)", 10.0233, 3, 0., 0.);
+  fPDGCode[fNPDGCodes++]=20553;       // 69 = Upsilon(2S)
+
+  Gspart(70, "Upsilon(3S)", 10.3553, 3, 0., 0.);
+  fPDGCode[fNPDGCodes++]=30553;       // 70 = Upsilon(3S)
+
 /* --- Define additional decay modes --- */
 /* --- omega(783) --- */
     for (kz = 0; kz < 6; ++kz) {
@@ -815,6 +1073,7 @@ void TGeant3::DefineParticles()
     mode[4] = 1701;
     Gsdk(ipa, bratio, mode);
 /* --- D+ --- */
+    /*
     for (kz = 0; kz < 6; ++kz) {
        bratio[kz] = 0.;
        mode[kz] = 0;
@@ -829,7 +1088,9 @@ void TGeant3::DefineParticles()
     mode[2] = 111208;
     mode[3] = 110809;
     Gsdk(ipa, bratio, mode);
+    */
 /* --- D- --- */
+    /*
     for (kz = 0; kz < 6; ++kz) {
        bratio[kz] = 0.;
        mode[kz] = 0;
@@ -844,7 +1105,9 @@ void TGeant3::DefineParticles()
     mode[2] = 121109;
     mode[3] = 120908;
     Gsdk(ipa, bratio, mode);
+    */
 /* --- D0 --- */
+    /*
     for (kz = 0; kz < 6; ++kz) {
        bratio[kz] = 0.;
        mode[kz] = 0;
@@ -857,7 +1120,9 @@ void TGeant3::DefineParticles()
     mode[1] = 1208;
     mode[2] = 1112;
     Gsdk(ipa, bratio, mode);
+    */
 /* --- Anti D0 --- */
+    /*
     for (kz = 0; kz < 6; ++kz) {
        bratio[kz] = 0.;
        mode[kz] = 0;
@@ -870,6 +1135,7 @@ void TGeant3::DefineParticles()
     mode[1] = 1109;
     mode[2] = 1112;
     Gsdk(ipa, bratio, mode);
+    */
 /* --- rho+ --- */
     for (kz = 0; kz < 6; ++kz) {
        bratio[kz] = 0.;
@@ -920,7 +1186,7 @@ void TGeant3::DefineParticles()
 }
 
 //_____________________________________________________________________________
-Int_t TGeant3::VolId(Text_t *name) const
+Int_t TGeant3::VolId(const Text_t *name) const
 {
   //
   // Return the unique numeric identifier for volume name
@@ -948,13 +1214,96 @@ const char* TGeant3::VolName(Int_t id) const
   //
   // Return the volume name given the volume identifier
   //
-  static char name[5];
   if(id<1 || id > fGcnum->nvolum || fGclink->jvolum<=0) 
-    strcpy(name,"NULL");
+    return fVolNames[fGcnum->nvolum];
   else
-    strncpy(name,(char *)&fZiq[fGclink->jvolum+id],4);
-  name[4]='\0';
-  return name;
+    return fVolNames[id-1];
+}
+
+//_____________________________________________________________________________
+void    TGeant3::SetCut(const char* cutName, Float_t cutValue)
+{
+  //
+  // Set transport cuts for particles
+  //
+  if(!strcmp(cutName,"CUTGAM")) 
+    fGccuts->cutgam=cutValue; 
+  else if(!strcmp(cutName,"CUTGAM")) 
+    fGccuts->cutele=cutValue; 
+  else if(!strcmp(cutName,"CUTELE")) 
+    fGccuts->cutneu=cutValue; 
+  else if(!strcmp(cutName,"CUTHAD")) 
+    fGccuts->cuthad=cutValue; 
+  else if(!strcmp(cutName,"CUTMUO")) 
+    fGccuts->cutmuo=cutValue; 
+  else if(!strcmp(cutName,"BCUTE")) 
+    fGccuts->bcute=cutValue; 
+  else if(!strcmp(cutName,"BCUTM")) 
+    fGccuts->bcutm=cutValue; 
+  else if(!strcmp(cutName,"DCUTE")) 
+    fGccuts->dcute=cutValue; 
+  else if(!strcmp(cutName,"DCUTM")) 
+    fGccuts->dcutm=cutValue; 
+  else if(!strcmp(cutName,"PPCUTM")) 
+    fGccuts->ppcutm=cutValue; 
+  else if(!strcmp(cutName,"TOFMAX")) 
+    fGccuts->tofmax=cutValue; 
+  else Warning("SetCut","Cut %s not implemented\n",cutName);
+}
+
+//_____________________________________________________________________________
+void    TGeant3::SetProcess(const char* flagName, Int_t flagValue)
+{
+  //
+  // Set thresholds for different processes
+  //
+  if(!strcmp(flagName,"PAIR")) 
+    fGcphys->ipair=flagValue;
+  else if(!strcmp(flagName,"COMP")) 
+    fGcphys->icomp=flagValue;
+  else if(!strcmp(flagName,"PHOT")) 
+    fGcphys->iphot=flagValue;
+  else if(!strcmp(flagName,"PFIS")) 
+    fGcphys->ipfis=flagValue;
+  else if(!strcmp(flagName,"DRAY")) 
+    fGcphys->idray=flagValue;
+  else if(!strcmp(flagName,"ANNI")) 
+    fGcphys->ianni=flagValue;
+  else if(!strcmp(flagName,"BREM")) 
+    fGcphys->ibrem=flagValue;
+  else if(!strcmp(flagName,"HADR")) 
+    fGcphys->ihadr=flagValue;
+  else if(!strcmp(flagName,"MUNU")) 
+    fGcphys->imunu=flagValue;
+  else if(!strcmp(flagName,"DCAY")) 
+    fGcphys->idcay=flagValue;
+  else if(!strcmp(flagName,"LOSS")) 
+    fGcphys->iloss=flagValue;
+  else if(!strcmp(flagName,"MULS")) 
+    fGcphys->imuls=flagValue;
+  else if(!strcmp(flagName,"RAYL")) 
+    fGcphys->irayl=flagValue;
+  else if(!strcmp(flagName,"STRA")) 
+    fGcphlt->istra=flagValue;
+  else if(!strcmp(flagName,"SYNC")) 
+    fGcphlt->isync=flagValue;
+  else  Warning("SetFlag","Flag %s not implemented\n",flagName);
+}
+
+//_____________________________________________________________________________
+Float_t TGeant3::Xsec(char* reac, Float_t /* energy */, 
+                     Int_t part, Int_t /* mate */)
+{
+  //
+  // Calculate X-sections -- dummy for the moment
+  //
+  if(!strcmp(reac,"PHOT"))
+  {
+    if(part!=22) {
+      Error("Xsec","Can calculate photoelectric only for photons\n");
+    }
+  }
+  return 0;
 }
 
 //_____________________________________________________________________________
@@ -1038,6 +1387,15 @@ Float_t TGeant3::TrackLength() const
   return fGctrak->sleng;
 }
 
+//_____________________________________________________________________________
+Bool_t TGeant3::IsNewTrack() const
+{
+  //
+  // True if the track is not at the boundary of the current volume
+  //
+  return (fGctrak->sleng==0);
+}
+
 //_____________________________________________________________________________
 Bool_t TGeant3::IsTrackInside() const
 {
@@ -1102,19 +1460,20 @@ Int_t   TGeant3::CurrentEvent() const
 }
 
 //_____________________________________________________________________________
-void    TGeant3::ProdProcess(char* proc) const
+const char* TGeant3::ProdProcess() const
 {
   //
   // Name of the process that has produced the secondary particles
   // in the current step
   //
-  const Int_t ipmec[13] = { 5,6,7,8,9,10,11,12,21,23,25,105,108 };
+  static char proc[5];
+  const Int_t kIpMec[13] = { 5,6,7,8,9,10,11,12,21,23,25,105,108 };
   Int_t mec, km, im;
   //
   if(fGcking->ngkine>0) {
     for (km = 0; km < fGctrak->nmec; ++km) {
       for (im = 0; im < 13; ++im) {
-       if (fGctrak->lmec[km] == ipmec[im]) {
+       if (fGctrak->lmec[km] == kIpMec[im]) {
          mec = fGctrak->lmec[km];
          if (0 < mec && mec < 31) {
            strncpy(proc,(char *)&fGctrak->namec[mec - 1],4);
@@ -1122,16 +1481,18 @@ void    TGeant3::ProdProcess(char* proc) const
            strncpy(proc,(char *)&fGctpol->namec1[mec - 101],4);
          }
          proc[4]='\0';
-         return;
+         return proc;
        }
       }
     }
     strcpy(proc,"UNKN");
   } else strcpy(proc,"NONE");
+  return proc;
 }
 
 //_____________________________________________________________________________
-void    TGeant3::GetSecondary(Int_t isec, Int_t& ipart, Float_t* x, Float_t* p)
+void    TGeant3::GetSecondary(Int_t isec, Int_t& ipart, 
+                             TLorentzVector &x, TLorentzVector &p)
 {
   //
   // Get the parameters of the secondary track number isec produced
@@ -1156,6 +1517,9 @@ void    TGeant3::GetSecondary(Int_t isec, Int_t& ipart, Float_t* x, Float_t* p)
 //_____________________________________________________________________________
 void TGeant3::InitLego()
 {
+  //
+  // Set switches for lego transport
+  //
   SetSWIT(4,0);
   SetDEBU(0,0,0);  //do not print a message 
 }
@@ -1208,40 +1572,6 @@ Float_t TGeant3::MaxStep() const
   return fGctmed->stemax;
 }
 
-//_____________________________________________________________________________
-void TGeant3::SetColors()
-{
-  //
-  // Set the colors for all the volumes
-  // this is done sequentially for all volumes
-  // based on the number of their medium
-  //
-  Int_t kv, icol;
-  Int_t jvolum=fGclink->jvolum;
-  //Int_t jtmed=fGclink->jtmed;
-  //Int_t jmate=fGclink->jmate;
-  Int_t nvolum=fGcnum->nvolum;
-  char name[5];
-  //
-  //    Now for all the volumes
-  for(kv=1;kv<=nvolum;kv++) {
-    //     Get the tracking medium
-    Int_t itm=Int_t (fZq[fZlq[jvolum-kv]+4]);
-    //     Get the material
-    //Int_t ima=Int_t (fZq[fZlq[jtmed-itm]+6]);
-    //     Get z
-    //Float_t z=fZq[fZlq[jmate-ima]+7];
-    //     Find color number
-    //icol = Int_t(z)%6+2;
-    //icol = 17+Int_t(z*150./92.);
-    //icol = kv%6+2;
-    icol = itm%6+2;
-    strncpy(name,(char*)&fZiq[jvolum+kv],4);
-    name[4]='\0';
-    Gsatt(name,"COLO",icol);
-  }
-}
-
 //_____________________________________________________________________________
 void TGeant3::SetMaxStep(Float_t maxstep)
 {
@@ -1411,31 +1741,6 @@ void TGeant3::Matrix(Int_t& krot, Float_t thex, Float_t phix, Float_t they,
   gsrotm(krot, thex, phix, they, phiy, thez, phiz);
 }
 
-//_____________________________________________________________________________
-void TGeant3::GetParticle(const Int_t pdg, char *name, Float_t &mass) const
-{
-  Int_t ipart = IdFromPDG(pdg);
-  if(ipart<0) {
-    printf("Particle %d not in geant\n",pdg);
-    name=new char[7];
-    strcpy(name,"Unknown");
-    mass=-1;
-    return;
-  }
-  //
-  // Return name and mass of particle code ipart
-  // Geant321 conventions
-  //
-  Int_t hname[6];
-  Int_t jpart=fGclink->jpart;
-  Int_t jpa=fZlq[jpart-ipart];
-  //
-  for(Int_t i=1; i<6; i++) hname[i-1]=fZiq[jpa+i];
-  hname[5]=0;
-  strncpy(name,(char *)hname, 21);
-  mass=fZq[jpa+7];
-}
-
 //_____________________________________________________________________________
 Int_t TGeant3::GetMedium() const
 {
@@ -1515,9 +1820,9 @@ void  TGeant3::Gpcxyz()
   //
   //    Print track and volume parameters at current point
   //
-  gpcxyz(); 
+    
+    gpcxyz(); 
 } 
 //_____________________________________________________________________________
 void  TGeant3::Ggclos() 
 { 
@@ -1541,6 +1846,14 @@ void  TGeant3::Ggclos()
   //   through the routine GHCLOS. 
   //
   ggclos(); 
+  // Create internal list of volumes
+  fVolNames = new char[fGcnum->nvolum+1][5];
+  Int_t i;
+  for(i=0; i<fGcnum->nvolum; ++i) {
+    strncpy(fVolNames[i], (char *) &fZiq[fGclink->jvolum+i+1], 4);
+    fVolNames[i][4]='\0';
+  }
+  strcpy(fVolNames[fGcnum->nvolum],"NULL");
 } 
  
 //_____________________________________________________________________________
@@ -1578,7 +1891,7 @@ void  TGeant3::Gtrig()
 { 
   //
   // Steering function to process one event
-  //
+  //  
   gtrig(); 
 } 
  
@@ -1645,7 +1958,8 @@ void  TGeant3::Gfpart(Int_t ipart, char *name, Int_t &itrtyp,
   //
   Float_t *ubuf=0; 
   Int_t   nbuf; 
-  gfpart(ipart, PASSCHARD(name), itrtyp, amass, charge, tlife, ubuf, nbuf
+  Int_t igpart = IdFromPDG(ipart);
+  gfpart(igpart, PASSCHARD(name), itrtyp, amass, charge, tlife, ubuf, nbuf
         PASSCHARL(name)); 
 } 
  
@@ -1660,7 +1974,38 @@ void  TGeant3::Gftmed(Int_t numed, char *name, Int_t &nmat, Int_t &isvol,
   //
   gftmed(numed, PASSCHARD(name), nmat, isvol, ifield, fieldm, tmaxfd, stemax,  
          deemax, epsil, stmin, ubuf, nbuf PASSCHARL(name)); 
+}
+
+ void  TGeant3::Gftmat(Int_t imate, Int_t ipart, char *chmeca, Int_t kdim,
+                     Float_t* tkin, Float_t* value, Float_t* pcut,
+                     Int_t &ixst)
+{ 
+  //
+  // Return parameters for tracking medium NUMED
+  //
+  gftmat(imate, ipart, PASSCHARD(chmeca), kdim,
+        tkin, value, pcut, ixst PASSCHARL(chmeca));
+
 } 
+
+//_____________________________________________________________________________
+Float_t TGeant3::Gbrelm(Float_t z, Float_t t, Float_t bcut)
+{
+  //
+  // To calculate energy loss due to soft muon BREMSSTRAHLUNG
+  //
+  return gbrelm(z,t,bcut);
+}
+
+//_____________________________________________________________________________
+Float_t TGeant3::Gprelm(Float_t z, Float_t t, Float_t bcut)
+{
+  //
+  // To calculate DE/DX in GeV*barn/atom for direct pair production by muons
+  //
+  return gprelm(z,t,bcut);
+}
  
 //_____________________________________________________________________________
 void  TGeant3::Gmate() 
@@ -2037,12 +2382,12 @@ void  TGeant3::Gtreve()
 } 
 
 //_____________________________________________________________________________
-void  TGeant3::Gtreve_root() 
+void  TGeant3::GtreveRoot() 
 { 
   //
   //   Controls tracking of all particles belonging to the current event
   //
-  gtreve_root(); 
+  gtreveroot(); 
 } 
 
 //_____________________________________________________________________________
@@ -2595,7 +2940,7 @@ void TGeant3::Gdopen(Int_t iview)
   //  with solid colours can now be stored in a view bank or in 'PICTURE FILES'
   //
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   gdopen(iview);
 }
  
@@ -2700,7 +3045,7 @@ void TGeant3::Gdraw(const char *name,Float_t theta, Float_t phi, Float_t psi,
   //  string for the NAME of the volume can be found using the command DTREE).
   //
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   char vname[5];
   Vname(name,vname);
   if (fGcvdma->raytra != 1) {
@@ -2731,7 +3076,7 @@ void TGeant3::Gdrawc(const char *name,Int_t axis, Float_t cut,Float_t u0,
   //  the CVOL/BOX function.
   //  
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   char vname[5];
   Vname(name,vname);
   gdrawc(PASSCHARD(vname), axis,cut,u0,v0,ul,vl PASSCHARL(vname)); 
@@ -2760,7 +3105,7 @@ void TGeant3::Gdrawx(const char *name,Float_t cutthe, Float_t cutphi,
   //  The resulting picture is seen from the viewing angles theta,phi.
   //
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   char vname[5];
   Vname(name,vname);
   gdrawx(PASSCHARD(vname), cutthe,cutphi,cutval,theta,phi,u0,v0,ul,vl
@@ -2831,7 +3176,7 @@ void TGeant3::Gdspec(const char *name)
   //  volume.
   //  
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   char vname[5];
   Vname(name,vname);
   gdspec(PASSCHARD(vname) PASSCHARL(vname)); 
@@ -2844,12 +3189,12 @@ void TGeant3::DrawOneSpec(const char *name)
   //  Function called when one double-clicks on a volume name
   //  in a TPavelabel drawn by Gdtree.
   //
-  THIGZ *higzSave = higz;
+  THIGZ *higzSave = gHigz;
   higzSave->SetName("higzSave");
   THIGZ *higzSpec = (THIGZ*)gROOT->FindObject("higzSpec");
-  //printf("DrawOneSpec, higz=%x, higzSpec=%x\n",higz,higzSpec);
-  if (higzSpec) higz     = higzSpec;
-  else          higzSpec = new THIGZ(defSize);
+  //printf("DrawOneSpec, gHigz=%x, higzSpec=%x\n",gHigz,higzSpec);
+  if (higzSpec) gHigz     = higzSpec;
+  else          higzSpec = new THIGZ(kDefSize);
   higzSpec->SetName("higzSpec");
   higzSpec->cd();
   higzSpec->Clear();
@@ -2859,7 +3204,7 @@ void TGeant3::DrawOneSpec(const char *name)
   higzSpec->Update();
   higzSave->cd();
   higzSave->SetName("higz");
-  higz = higzSave;
+  gHigz = higzSave;
 } 
 
 //_____________________________________________________________________________
@@ -2879,11 +3224,11 @@ void TGeant3::Gdtree(const char *name,Int_t levmax, Int_t isel)
   //    - drawing tree of parent
   //  
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   char vname[5];
   Vname(name,vname);
   gdtree(PASSCHARD(vname), levmax, isel PASSCHARL(vname)); 
-  higz->fPname = "";
+  gHigz->SetPname("");
 } 
 
 //_____________________________________________________________________________
@@ -2897,7 +3242,7 @@ void TGeant3::GdtreeParent(const char *name,Int_t levmax, Int_t isel)
   //  This function draws the logical tree of the parent of name.
   //  
   InitHIGZ();
-  higz->Clear();
+  gHigz->Clear();
   // Scan list of volumes in JVOLUM
   char vname[5];
   Int_t gname, i, jvo, in, nin, jin, num;
@@ -2913,7 +3258,7 @@ void TGeant3::GdtreeParent(const char *name,Int_t levmax, Int_t isel)
        strncpy(vname,(char*)&fZiq[fGclink->jvolum+i],4);
        vname[4] = 0;           
        gdtree(PASSCHARD(vname), levmax, isel PASSCHARL(vname)); 
-       higz->fPname = "";
+       gHigz->SetPname("");
        return;
       }
     }
@@ -3135,6 +3480,20 @@ void TGeant3::SetDRAY(Int_t par)
   fGcphys->idray = par;
 }
  
+//_____________________________________________________________________________
+void TGeant3::SetERAN(Float_t ekmin, Float_t ekmax, Int_t nekbin)
+{
+  //
+  //  To control cross section tabulations
+  //   ekmin = minimum kinetic energy in GeV
+  //   ekmax = maximum kinetic energy in GeV
+  //   nekbin = number of logatithmic bins (<200)
+  //  
+  fGcmulo->ekmin = ekmin;
+  fGcmulo->ekmax = ekmax;
+  fGcmulo->nekbin = nekbin;
+}
 //_____________________________________________________________________________
 void TGeant3::SetHADR(Int_t par)
 {
@@ -3274,6 +3633,18 @@ void TGeant3::SetRAYL(Int_t par)
   fGcphys->irayl = par;
 }
  
+//_____________________________________________________________________________
+void TGeant3::SetSTRA(Int_t par)
+{
+  //
+  //  To control energy loss fluctuations
+  //  with the PhotoAbsorption Ionisation model.
+  //   par =0 no Straggling.
+  //       =1 Straggling yes => no Delta rays.
+  //  
+  fGcphlt->istra = par;
+}
 //_____________________________________________________________________________
 void TGeant3::SetSWIT(Int_t sw, Int_t val)
 {
@@ -3337,6 +3708,56 @@ void TGeant3::Vname(const char *name, char *vname)
   vname[4] = 0;      
 }
  
+//______________________________________________________________________________
+void TGeant3::Ertrgo()
+{
+  //
+  // Perform the tracking of the track Track parameters are in VECT
+  //
+  ertrgo();
+}
+
+//______________________________________________________________________________
+void TGeant3::Ertrak(const Float_t *const x1, const Float_t *const p1, 
+                       const Float_t *x2, const Float_t *p2,
+                       Int_t ipa,  Option_t *chopt)
+{
+  //************************************************************************
+  //*                                                                      *
+  //*          Perform the tracking of the track from point X1 to          *
+  //*                    point X2                                          *
+  //*          (Before calling this routine the user should also provide   *
+  //*                    the input informations in /EROPTS/ and /ERTRIO/   *
+  //*                    using subroutine EUFIL(L/P/V)                     *
+  //*                 X1       - Starting coordinates (Cartesian)          *
+  //*                 P1       - Starting 3-momentum  (Cartesian)          *
+  //*                 X2       - Final coordinates    (Cartesian)          *
+  //*                 P2       - Final 3-momentum     (Cartesian)          *
+  //*                 IPA      - Particle code (a la GEANT) of the track   *
+  //*                                                                      *
+  //*                 CHOPT                                                *
+  //*                     'B'   'Backward tracking' - i.e. energy loss     *
+  //*                                        added to the current energy   *
+  //*                     'E'   'Exact' calculation of errors assuming     *
+  //*                                        helix (i.e. pathlength not    *
+  //*                                        assumed as infinitesimal)     *
+  //*                     'L'   Tracking upto prescribed Lengths reached   *
+  //*                     'M'   'Mixed' prediction (not yet coded)         *
+  //*                     'O'   Tracking 'Only' without calculating errors *
+  //*                     'P'   Tracking upto prescribed Planes reached    *
+  //*                     'V'   Tracking upto prescribed Volumes reached   *
+  //*                     'X'   Tracking upto prescribed Point approached  *
+  //*                                                                      *
+  //*                Interface with GEANT :                                *
+  //*             Track parameters are in /CGKINE/ and /GCTRAK/            *
+  //*                                                                      *
+  //*          ==>Called by : USER                                         *
+  //*             Authors   M.Maire, E.Nagy  ********//*                     *
+  //*                                                                      *
+  //************************************************************************
+  ertrak(x1,p1,x2,p2,ipa,PASSCHARD(chopt) PASSCHARL(chopt));
+}
+        
 //_____________________________________________________________________________
 void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
                          Int_t number, Int_t nlevel)
@@ -3366,7 +3787,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
   //     this procedure fails, ordering starts from 1.
   //     Arrays IOTMED and IOMATE are used for this procedure
   //
-  const char shape[][5]={"BOX ","TRD1","TRD2","TRAP","TUBE","TUBS","CONE",
+  const char kShape[][5]={"BOX ","TRD1","TRD2","TRAP","TUBE","TUBS","CONE",
                         "CONS","SPHE","PARA","PGON","PCON","ELTU","HYPE",
                         "GTRA","CTUB"};
   Int_t i, end, itm, irm, jrm, k, nmed;
@@ -3395,15 +3816,15 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
   Int_t itmedc, nmatc, isvolc, ifieldc, nwbufc, isvol, nmat, ifield, nwbuf;
   Float_t fieldmc, tmaxfdc, stemaxc, deemaxc, epsilc, stminc, fieldm;
   Float_t tmaxf, stemax, deemax, epsil, stmin;
-  const char *f10000="!\n%s\n!\n";
+  const char *k10000="!\n%s\n!\n";
   //Open the input file
   end=strlen(filnam);
   for(i=0;i<end;i++) if(filnam[i]=='.') {
     end=i;
     break;
   }
-  filext=new char[end+4];
-  filetme=new char[end+4];
+  filext=new char[end+5];
+  filetme=new char[end+5];
   strncpy(filext,filnam,end);
   strncpy(filetme,filnam,end);
   //
@@ -3596,7 +4017,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
   // *** write down the tracking medium definition
   //
   strcpy(card,"!       Tracking medium");
-  fprintf(lun,f10000,card);
+  fprintf(lun,k10000,card);
   //
   for(itm=1;itm<=fGcnum->ntmed;itm++) {
     if (iws[iadtmd+itm]>0) {
@@ -3626,7 +4047,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
       //* *** write down the rotation matrix
   //*
   strcpy(card,"!       Reperes");
-  fprintf(lun,f10000,card);
+  fprintf(lun,k10000,card);
   //
   for(irm=1;irm<=fGcnum->nrotm;irm++) {
     if (iws[iadrot+irm]>0) {
@@ -3640,7 +4061,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
   //* *** write down the volume definition
   //*
   strcpy(card,"!       Volumes");
-  fprintf(lun,f10000,card);
+  fprintf(lun,k10000,card);
   //*
   for(ivstak=1;ivstak<=nvstak;ivstak++) {
     ivo = iws[ivstak];
@@ -3654,7 +4075,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
       if (npar>0) {
        if (ivstak>1) for(i=0;i<npar;i++) par[i]=fZq[jvo+7+i];
        Gckpar (ish,npar,par);
-       fprintf(lun,"VOLU '%4s' '%4s' %3d %3d\n",name,shape[ish-1],iotmed[nmed],npar);
+       fprintf(lun,"VOLU '%4s' '%4s' %3d %3d\n",name,kShape[ish-1],iotmed[nmed],npar);
        for(i=0;i<(npar-1)/6+1;i++) {
          fprintf(lun,"     ");
          left=npar-i*6;
@@ -3662,14 +4083,14 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
          fprintf(lun,"\n");
        }
       } else {
-       fprintf(lun,"VOLU '%4s' '%4s' %3d %3d\n",name,shape[ish-1],iotmed[nmed],npar);
+       fprintf(lun,"VOLU '%4s' '%4s' %3d %3d\n",name,kShape[ish-1],iotmed[nmed],npar);
       }
     }
   }
   //*
   //* *** write down the division of volumes
   //*
-  fprintf(lun,f10000,"!       Divisions");
+  fprintf(lun,k10000,"!       Divisions");
   for(ivstak=1;ivstak<=nvstak;ivstak++) {
     ivo = TMath::Abs(iws[ivstak]);
     jvo  = fZlq[fGclink->jvolum-ivo];
@@ -3707,7 +4128,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
   //*
   //* *** write down the the positionnement of volumes
   //*
-  fprintf(lun,f10000,"!       Positionnements\n");
+  fprintf(lun,k10000,"!       Positionnements\n");
   //
   for(ivstak = 1;ivstak<=nvstak;ivstak++) {
     ivo = TMath::Abs(iws[ivstak]);
@@ -3788,6 +4209,7 @@ void TGeant3::WriteEuclid(const char* filnam, const char* topvol,
     }
   }
   fprintf(lun,"END\n");
+  fclose(lun);
   printf(" *** GWEUCL *** file: %s is now written out\n",filext);
   printf(" *** GWEUCL *** file: %s is now written out\n",filetme);
   // Clean up