TVirtualMC.h

Go to the documentation of this file.

// @(#)root/vmc:$Name:  $:$Id$
// Authors: Ivana Hrivnacova, Rene Brun, Federico Carminati 13/04/2002

/*************************************************************************
 * Copyright (C) 2006, Rene Brun and Fons Rademakers.                    *
 * Copyright (C) 2002, ALICE Experiment at CERN.                         *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef ROOT_TVirtualMC
#define ROOT_TVirtualMC

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//                                                                           //
//   Abstract Monte Carlo interface                                          //
//                                                                           //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include "TMCProcess.h"
#include "TMCParticleType.h"
#include "TMCOptical.h"
#include "TMCtls.h"
#include "TVirtualMCApplication.h"
#include "TVirtualMCStack.h"
#include "TVirtualMCDecayer.h"
#include "TVirtualMagField.h"
#include "TRandom.h"
#include "TString.h"
#include "TError.h"

class TLorentzVector;
class TGeoHMatrix;
class TArrayI;
class TArrayD;

class TVirtualMC : public TNamed {

public:
   /// Standard constructor
   ///
   /// isRootGeometrySupported = True if implementation of TVirtualMC
   ///        supports geometry defined with TGeo
   TVirtualMC(const char *name, const char *title,
              Bool_t isRootGeometrySupported = kFALSE);

   /// Default constructor
   TVirtualMC();

   /// Destructor
   virtual ~TVirtualMC();

   /// Static access method
   static TVirtualMC* GetMC();

   //
   // ------------------------------------------------
   // methods for building/management of geometry
   // ------------------------------------------------
   //

   /// Info about supporting geometry defined via Root
   virtual Bool_t IsRootGeometrySupported() const = 0;

   //
   // functions from GCONS
   // ------------------------------------------------
   //

   /// Define a material
   /// - kmat   number assigned to the material
   /// - name   material name
   /// - a      atomic mass in au
   /// - z      atomic number
   /// - dens   density in g/cm3
   /// - absl   absorption length in cm;
   ///               if >=0 it is ignored and the program
   ///               calculates it, if <0. -absl is taken
   /// - radl   radiation length in cm
   ///               if >=0 it is ignored and the program
   ///               calculates it, if <0. -radl is taken
   /// - buf    pointer to an array of user words
   /// - nwbuf  number of user words
   virtual void  Material(Int_t& kmat, const char* name, Double_t a,
                    Double_t z, Double_t dens, Double_t radl, Double_t absl,
                    Float_t* buf, Int_t nwbuf) = 0;

   /// The same as previous but in double precision
   virtual void  Material(Int_t& kmat, const char* name, Double_t a,
                     Double_t z, Double_t dens, Double_t radl, Double_t absl,
                     Double_t* buf, Int_t nwbuf) = 0;

   /// Define a mixture or a compound
   /// with a number kmat composed by the basic nlmat materials defined
   /// by arrays a, z and wmat
   ///
   /// If nlmat > 0 then wmat contains the proportion by
   /// weights of each basic material in the mixture.
   ///
   /// If nlmat < 0 then wmat contains the number of atoms
   /// of a given kind into the molecule of the compound.
   /// In this case, wmat in output is changed to relative
   /// weights.
   virtual void  Mixture(Int_t& kmat, const char *name, Float_t *a,
                     Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) = 0;

   /// The same as previous but in double precision
   virtual void  Mixture(Int_t& kmat, const char *name, Double_t *a,
                     Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) = 0;

   /// Define a medium.
   /// - kmed      tracking medium number assigned
   /// - name      tracking medium name
   /// - nmat      material number
   /// - isvol     sensitive volume flag
   /// - ifield    magnetic field:
   ///                  - ifield = 0 if no magnetic field;
   ///                  - ifield = -1 if user decision in guswim;
   ///                  - ifield = 1 if tracking performed with g3rkuta;
   ///                  - ifield = 2 if tracking performed with g3helix;
   ///                  - ifield = 3 if tracking performed with g3helx3.
   /// - fieldm    max. field value (kilogauss)
   /// - tmaxfd    max. angle due to field (deg/step)
   /// - stemax    max. step allowed
   /// - deemax    max. fraction of energy lost in a step
   /// - epsil     tracking precision (cm)
   /// - stmin     min. step due to continuous processes (cm)
   /// - ubuf      pointer to an array of user words
   /// - nbuf      number of user words
   virtual void  Medium(Int_t& kmed, const char *name, Int_t nmat,
                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
                     Double_t stemax, Double_t deemax, Double_t epsil,
                     Double_t stmin, Float_t* ubuf, Int_t nbuf) = 0;

   /// The same as previous but in double precision
   virtual void  Medium(Int_t& kmed, const char *name, Int_t nmat,
                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
                     Double_t stemax, Double_t deemax, Double_t epsil,
                     Double_t stmin, Double_t* ubuf, Int_t nbuf) = 0;

   /// Define a rotation matrix
   /// - krot     rotation matrix number assigned
   /// - thetaX   polar angle for axis X
   /// - phiX     azimuthal angle for axis X
   /// - thetaY   polar angle for axis Y
   /// - phiY     azimuthal angle for axis Y
   /// - thetaZ   polar angle for axis Z
   /// - phiZ     azimuthal angle for axis Z
   virtual void  Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
                     Double_t phiZ) = 0;

   /// Change the value of cut or mechanism param
   /// to a new value parval for tracking medium itmed.
   /// In Geant3, the  data  structure JTMED contains the standard tracking
   /// parameters (CUTS and flags to control the physics processes)  which
   /// are used  by default for all tracking media.
   /// It is possible to redefine individually with this function any of these
   /// parameters for a given tracking medium.
   /// - itmed   tracking medium number
   /// - param   is a character string (variable name)
   /// - parval  must be given as a floating point.
   virtual void  Gstpar(Int_t itmed, const char *param, Double_t parval) = 0;

   //
   // functions from GGEOM
   // ------------------------------------------------
   //

   /// Create a new volume
   /// - name   Volume name
   /// - shape  Volume type
   /// - nmed   Tracking medium number
   /// - np     Number of shape parameters
   /// - upar   Vector containing shape parameters
   virtual Int_t  Gsvolu(const char *name, const char *shape, Int_t nmed,
                          Float_t *upar, Int_t np) = 0;

   /// The same as previous but in double precision
   virtual Int_t  Gsvolu(const char *name, const char *shape, Int_t nmed,
                          Double_t *upar, Int_t np) = 0;

   /// Create a new volume by dividing an existing one.
   /// It divides a previously defined volume
   /// - name   Volume name
   /// - mother Mother volume name
   /// - ndiv   Number of divisions
   /// - iaxis  Axis value:
   ///               X,Y,Z of CAXIS will be translated to 1,2,3 for IAXIS.
   virtual void  Gsdvn(const char *name, const char *mother, Int_t ndiv,
                         Int_t iaxis) = 0;

   /// Create a new volume by dividing an existing one.
   /// Divide mother into ndiv divisions called name
   /// along axis iaxis starting at coordinate value c0i.
   /// The new volume created will be medium number numed.
   virtual void  Gsdvn2(const char *name, const char *mother, Int_t ndiv,
                         Int_t iaxis, Double_t c0i, Int_t numed) = 0;

   /// Create a new volume by dividing an existing one
   /// Divide mother into divisions called name along
   /// axis iaxis in steps of step. If not exactly divisible
   /// will make as many as possible and will center them
   /// with respect to the mother. Divisions will have medium
   /// number numed. If numed is 0, numed of mother is taken.
   /// ndvmx is the expected maximum number of divisions
   /// (If 0, no protection tests are performed in Geant3)
   virtual void  Gsdvt(const char *name, const char *mother, Double_t step,
                         Int_t iaxis, Int_t numed, Int_t ndvmx) = 0;

   /// Create a new volume by dividing an existing one
   /// Divides mother into divisions called name along
   /// axis iaxis starting at coordinate value c0 with step
   /// size step.
   /// The new volume created will have medium number numed.
   /// If numed is 0, numed of mother is taken.
   /// ndvmx is the expected maximum number of divisions
   /// (If 0, no protection tests are performed in Geant3)
   virtual void  Gsdvt2(const char *name, const char *mother, Double_t step,
                         Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) = 0;

   /// Flag volume name whose contents will have to be ordered
   /// along axis iax, by setting the search flag to -iax
   /// (Geant3 only)
   virtual void  Gsord(const char *name, Int_t iax) = 0;

   /// Position a volume into an existing one.
   /// It positions a previously defined volume in the mother.
   /// - name   Volume name
   /// - nr     Copy number of the volume
   /// - mother Mother volume name
   /// - x      X coord. of the volume in mother ref. sys.
   /// - y      Y coord. of the volume in mother ref. sys.
   /// - z      Z coord. of the volume in mother ref. sys.
   /// - irot   Rotation matrix number w.r.t. mother ref. sys.
   /// - konly  ONLY/MANY flag
   virtual void  Gspos(const char *name, Int_t nr, const char *mother,
                         Double_t x, Double_t y, Double_t z, Int_t irot,
                         const char *konly="ONLY") = 0;

   /// Place a copy of generic volume name with user number
   ///  nr inside mother, with its parameters upar(1..np)
   virtual void  Gsposp(const char *name, Int_t nr, const char *mother,
                         Double_t x, Double_t y, Double_t z, Int_t irot,
                         const char *konly, Float_t *upar, Int_t np) = 0;

   /// The same as previous but in double precision
   virtual void  Gsposp(const char *name, Int_t nr, const char *mother,
                         Double_t x, Double_t y, Double_t z, Int_t irot,
                         const char *konly, Double_t *upar, Int_t np) = 0;

   /// Helper function for resolving MANY.
   /// Specify the ONLY volume that overlaps with the
   /// specified MANY and has to be substracted.
   /// (Geant4 only)
   virtual void  Gsbool(const char* onlyVolName, const char* manyVolName) = 0;

   /// Define the tables for UV photon tracking in medium itmed.
   /// Please note that it is the user's responsibility to
   /// provide all the coefficients:
   /// - itmed       Tracking medium number
   /// - npckov      Number of bins of each table
   /// - ppckov      Value of photon momentum (in GeV)
   /// - absco       Absorption coefficients
   ///                     - dielectric: absorption length in cm
   ///                     - metals    : absorption fraction (0<=x<=1)
   /// - effic       Detection efficiency for UV photons
   /// - rindex      Refraction index (if=0 metal)
   virtual void  SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
                               Float_t *absco, Float_t *effic, Float_t *rindex) = 0;

   /// The same as previous but in double precision
   virtual void  SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
                               Double_t *absco, Double_t *effic, Double_t *rindex) = 0;

   //
   // functions for definition of surfaces
   // and material properties for optical physics
   // ------------------------------------------------
   //

   /// Define the optical surface
   /// - name           surface name
   /// - model          selection of model (see #EMCOpSurfaceModel values)
   /// - surfaceType    surface type (see #EMCOpSurfaceType values)
   /// - surfaceFinish  surface quality (see #EMCOpSurfaceType values)
   /// - sigmaAlpha     an unified model surface parameter
   /// (Geant4 only)
   virtual void  DefineOpSurface(const char* name,
                         EMCOpSurfaceModel model,
                         EMCOpSurfaceType surfaceType,
                         EMCOpSurfaceFinish surfaceFinish,
                         Double_t sigmaAlpha) = 0;

   /// Define the optical surface border
   /// - name        border surface name
   /// - vol1Name    first volume name
   /// - vol1CopyNo  first volume copy number
   /// - vol2Name    second volume name
   /// - vol2CopyNo  second volume copy number
   /// - opSurfaceName  name of optical surface which this border belongs to
   /// (Geant4 only)
   virtual void  SetBorderSurface(const char* name,
                         const char* vol1Name, int vol1CopyNo,
                         const char* vol2Name, int vol2CopyNo,
                         const char* opSurfaceName) = 0;

   /// Define the optical skin surface
   /// - name        skin surface name
   /// - volName     volume name
   /// - opSurfaceName  name of optical surface which this border belongs to
   /// (Geant4 only)
   virtual void  SetSkinSurface(const char* name,
                         const char* volName,
                         const char* opSurfaceName) = 0;

   /// Define material property via a table of values
   /// - itmed         tracking medium id
   /// - propertyName  property name
   /// - np            number of bins of the table
   /// - pp            value of photon momentum (in GeV)
   /// - values        property values
   /// (Geant4 only)
   virtual void  SetMaterialProperty(
                         Int_t itmed, const char* propertyName,
                         Int_t np, Double_t* pp, Double_t* values) = 0;

   /// Define material property via a value
   /// - itmed         tracking medium id
   /// - propertyName  property name
   /// - value         property value
   /// (Geant4 only)
   virtual void  SetMaterialProperty(
                         Int_t itmed, const char* propertyName,
                         Double_t value) = 0;

   /// Define optical surface property via a table of values
   /// - surfaceName   optical surface name
   /// - propertyName  property name
   /// - np            number of bins of the table
   /// - pp            value of photon momentum (in GeV)
   /// - values        property values
   /// (Geant4 only)
   virtual void  SetMaterialProperty(
                         const char* surfaceName, const char* propertyName,
                         Int_t np, Double_t* pp, Double_t* values) = 0;

   //
   // functions for access to geometry
   // ------------------------------------------------
   //

   /// Return the transformation matrix between the volume specified by
   /// the path volumePath and the top or master volume.
   virtual Bool_t GetTransformation(const TString& volumePath,
                         TGeoHMatrix& matrix) = 0;

   /// Return the name of the shape (shapeType)  and its parameters par
   /// for the volume specified by the path volumePath .
   virtual Bool_t GetShape(const TString& volumePath,
                         TString& shapeType, TArrayD& par) = 0;

   /// Return the material parameters for the material specified by
   /// the material Id
   virtual Bool_t GetMaterial(Int_t imat, TString& name,
                               Double_t& a, Double_t& z, Double_t& density,
                               Double_t& radl, Double_t& inter, TArrayD& par) = 0;

   /// Return the material parameters for the volume specified by
   /// the volumeName.
   virtual Bool_t GetMaterial(const TString& volumeName,
                               TString& name, Int_t& imat,
                               Double_t& a, Double_t& z, Double_t& density,
                               Double_t& radl, Double_t& inter, TArrayD& par) = 0;

   /// Return the medium parameters for the volume specified by the
   /// volumeName.
   virtual Bool_t GetMedium(const TString& volumeName,
                             TString& name, Int_t& imed,
                             Int_t& nmat, Int_t& isvol, Int_t& ifield,
                             Double_t& fieldm, Double_t& tmaxfd, Double_t& stemax,
                             Double_t& deemax, Double_t& epsil, Double_t& stmin,
                             TArrayD& par) = 0;

   /// Write out the geometry of the detector in EUCLID file format
   /// - filnam  file name - will be with the extension .euc                 *
   /// - topvol  volume name of the starting node
   /// - number  copy number of topvol (relevant for gsposp)
   /// - nlevel  number of  levels in the tree structure
   ///                to be written out, starting from topvol
   /// (Geant3 only)
   /// Deprecated
   virtual void  WriteEuclid(const char* filnam, const char* topvol,
                             Int_t number, Int_t nlevel) = 0;

   /// Set geometry from Root (built via TGeo)
   virtual void  SetRootGeometry() = 0;

   /// Activate the parameters defined in tracking media
   /// (DEEMAX, STMIN, STEMAX), which are, be default, ignored.
   /// In Geant4 case, only STEMAX is taken into account.
   /// In FLUKA, all tracking media parameters are ignored.
   virtual void SetUserParameters(Bool_t isUserParameters) = 0;

   //
   // get methods
   // ------------------------------------------------
   //

   /// Return the unique numeric identifier for volume name volName
   virtual Int_t VolId(const char* volName) const = 0;

   /// Return the volume name for a given volume identifier id
   virtual const char* VolName(Int_t id) const = 0;

   /// Return the unique numeric identifier for medium name mediumName
   virtual Int_t MediumId(const char* mediumName) const = 0;

   /// Return total number of volumes in the geometry
   virtual Int_t NofVolumes() const = 0;

   /// Return material number for a given volume id
   virtual Int_t VolId2Mate(Int_t id) const = 0;

   /// Return number of daughters of the volume specified by volName
   virtual Int_t NofVolDaughters(const char* volName) const = 0;

   /// Return the name of i-th daughter of the volume specified by volName
   virtual const char*  VolDaughterName(const char* volName, Int_t i) const = 0;

   /// Return the copyNo of i-th daughter of the volume specified by volName
   virtual Int_t        VolDaughterCopyNo(const char* volName, Int_t i) const = 0;

   //
   // ------------------------------------------------
   // methods for physics management
   // ------------------------------------------------
   //

   //
   // set methods
   // ------------------------------------------------
   //

   /// Set transport cuts for particles
   virtual Bool_t   SetCut(const char* cutName, Double_t cutValue) = 0;

   /// Set process control
   virtual Bool_t   SetProcess(const char* flagName, Int_t flagValue) = 0;

   /// Set a user defined particle
   /// Function is ignored if particle with specified pdg
   /// already exists and error report is printed.
   /// - pdg           PDG encoding
   /// - name          particle name
   /// - mcType        VMC Particle type
   /// - mass          mass [GeV]
   /// - charge        charge [eplus]
   /// - lifetime      time of life [s]
   /// - pType         particle type as in Geant4
   /// - width         width [GeV]
   /// - iSpin         spin
   /// - iParity       parity
   /// - iConjugation  conjugation
   /// - iIsospin      isospin
   /// - iIsospinZ     isospin - #rd component
   /// - gParity       gParity
   /// - lepton        lepton number
   /// - baryon        baryon number
   /// - stable        stability
   /// - shortlived    is shorlived?
   /// - subType       particle subType as in Geant4
   /// - antiEncoding  anti encoding
   /// - magMoment     magnetic moment
   /// - excitation    excitation energy [GeV]
   virtual Bool_t   DefineParticle(Int_t pdg, const char* name,
                        TMCParticleType mcType,
                        Double_t mass, Double_t charge, Double_t lifetime) = 0;

   /// Set a user defined particle
   /// Function is ignored if particle with specified pdg
   /// already exists and error report is printed.
   /// - pdg           PDG encoding
   /// - name          particle name
   /// - mcType        VMC Particle type
   /// - mass          mass [GeV]
   /// - charge        charge [eplus]
   /// - lifetime      time of life [s]
   /// - pType         particle type as in Geant4
   /// - width         width [GeV]
   /// - iSpin         spin
   /// - iParity       parity
   /// - iConjugation  conjugation
   /// - iIsospin      isospin
   /// - iIsospinZ     isospin - #rd component
   /// - gParity       gParity
   /// - lepton        lepton number
   /// - baryon        baryon number
   /// - stable        stability
   /// - shortlived    is shorlived?
   /// - subType       particle subType as in Geant4
   /// - antiEncoding  anti encoding
   /// - magMoment     magnetic moment
   /// - excitation    excitation energy [GeV]
   virtual Bool_t   DefineParticle(Int_t pdg, const char* name,
                        TMCParticleType mcType,
                        Double_t mass, Double_t charge, Double_t lifetime,
                        const TString& pType, Double_t width,
                        Int_t iSpin, Int_t iParity, Int_t iConjugation,
                        Int_t iIsospin, Int_t iIsospinZ, Int_t gParity,
                        Int_t lepton, Int_t baryon,
                        Bool_t stable, Bool_t shortlived = kFALSE,
                        const TString& subType = "",
                        Int_t antiEncoding = 0, Double_t magMoment = 0.0,
                        Double_t excitation = 0.0) = 0;

   /// Set a user defined ion.
   /// - name          ion name
   /// - Z             atomic number
   /// - A             atomic mass
   /// - Q             charge [eplus}
   /// - excitation    excitation energy [GeV]
   /// - mass          mass  [GeV] (if not specified by user, approximative
   ///                 mass is calculated)
   virtual Bool_t   DefineIon(const char* name, Int_t Z, Int_t A,
                        Int_t Q, Double_t excEnergy, Double_t mass = 0.) = 0;

   /// Set a user phase space decay for a particle
   /// -  pdg           particle PDG encoding
   /// -  bratios       the array with branching ratios (in %)
   /// -  mode[6][3]    the array with daughters particles PDG codes  for each
   ///                 decay channel
   virtual Bool_t   SetDecayMode(Int_t pdg, Float_t bratio[6], Int_t mode[6][3]) = 0;

   /// Calculate X-sections
   /// (Geant3 only)
   /// Deprecated
   virtual Double_t Xsec(char*, Double_t, Int_t, Int_t) = 0;

   //
   // particle table usage
   // ------------------------------------------------
   //

   /// Return MC specific code from a PDG and pseudo ENDF code (pdg)
   virtual Int_t   IdFromPDG(Int_t pdg) const =0;

   /// Return PDG code and pseudo ENDF code from MC specific code (id)
   virtual Int_t   PDGFromId(Int_t id) const =0;

   //
   // get methods
   // ------------------------------------------------
   //

   /// Return name of the particle specified by pdg.
   virtual TString   ParticleName(Int_t pdg) const = 0;

   /// Return mass of the particle specified by pdg.
   virtual Double_t  ParticleMass(Int_t pdg) const = 0;

   /// Return charge (in e units) of the particle specified by pdg.
   virtual Double_t  ParticleCharge(Int_t pdg) const = 0;

   /// Return life time of the particle specified by pdg.
   virtual Double_t  ParticleLifeTime(Int_t pdg) const = 0;

   /// Return VMC type of the particle specified by pdg.
   virtual TMCParticleType ParticleMCType(Int_t pdg) const = 0;
   //
   // ------------------------------------------------
   // methods for step management
   // ------------------------------------------------
   //

   //
   // action methods
   // ------------------------------------------------
   //

   /// Stop the transport of the current particle and skip to the next
   virtual void StopTrack() = 0;

   /// Stop simulation of the current event and skip to the next
   virtual void StopEvent() = 0;

   /// Stop simulation of the current event and set the abort run flag to true
   virtual void StopRun() = 0;

   //
   // set methods
   // ------------------------------------------------
   //

   /// Set the maximum step allowed till the particle is in the current medium
   virtual void SetMaxStep(Double_t) = 0;

   /// Set the maximum number of steps till the particle is in the current medium
   virtual void SetMaxNStep(Int_t) = 0;

   /// Force the decays of particles to be done with Pythia
   /// and not with the Geant routines.
   virtual void SetUserDecay(Int_t pdg) = 0;

   /// Force the decay time of the current particle
   virtual void ForceDecayTime(Float_t) = 0;

   //
   // tracking volume(s)
   // ------------------------------------------------
   //

   /// Return the current volume ID and copy number
   virtual Int_t    CurrentVolID(Int_t& copyNo) const =0;

   /// Return the current volume off upward in the geometrical tree
   /// ID and copy number
   virtual Int_t    CurrentVolOffID(Int_t off, Int_t& copyNo) const =0;

   /// Return the current volume name
   virtual const char* CurrentVolName() const =0;

   /// Return the current volume off upward in the geometrical tree
   /// name and copy number'
   /// if name=0 no name is returned
   virtual const char* CurrentVolOffName(Int_t off) const =0;

   /// Return the path in geometry tree for the current volume
   virtual const char* CurrentVolPath() = 0;

   /// If track is on a geometry boundary, fill the normal vector of the crossing
   /// volume surface and return true, return false otherwise
   virtual Bool_t   CurrentBoundaryNormal(
                       Double_t &x, Double_t &y, Double_t &z) const = 0;

   /// Return the parameters of the current material during transport
   virtual Int_t    CurrentMaterial(Float_t &a, Float_t &z,
                       Float_t &dens, Float_t &radl, Float_t &absl) const =0;

   //// Return the number of the current medium
   virtual Int_t    CurrentMedium() const = 0;
                         // new function (to replace GetMedium() const)

   /// Return the number of the current event
   virtual Int_t    CurrentEvent() const =0;

   /// Computes coordinates xd in daughter reference system
   /// from known coordinates xm in mother reference system.
   /// - xm    coordinates in mother reference system (input)
   /// - xd    coordinates in daughter reference system (output)
   /// - iflag
   ///   - IFLAG = 1  convert coordinates
   ///   - IFLAG = 2  convert direction cosines
   virtual void     Gmtod(Float_t* xm, Float_t* xd, Int_t iflag) = 0;

   /// The same as previous but in double precision
   virtual void     Gmtod(Double_t* xm, Double_t* xd, Int_t iflag) = 0;

   /// Computes coordinates xm in mother reference system
   /// from known coordinates xd in daughter reference system.
   /// - xd    coordinates in daughter reference system (input)
   /// - xm    coordinates in mother reference system (output)
   /// - iflag
   ///   - IFLAG = 1  convert coordinates
   ///   - IFLAG = 2  convert direction cosines
   virtual void     Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)= 0 ;

   /// The same as previous but in double precision
   virtual void     Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)= 0 ;

   /// Return the maximum step length in the current medium
   virtual Double_t MaxStep() const =0;

   /// Return the maximum number of steps allowed in the current medium
   virtual Int_t    GetMaxNStep() const = 0;

   //
   // get methods
   // tracking particle
   // dynamic properties
   // ------------------------------------------------
   //

   /// Return the current position in the master reference frame of the
   /// track being transported
   virtual void     TrackPosition(TLorentzVector& position) const =0;

   /// Return the current position in the master reference frame of the
   /// track being transported
   virtual void     TrackPosition(Double_t &x, Double_t &y, Double_t &z) const =0;

   /// Return the direction and the momentum (GeV/c) of the track
   /// currently being transported
   virtual void     TrackMomentum(TLorentzVector& momentum) const =0;

   /// Return the direction and the momentum (GeV/c) of the track
   /// currently being transported
   virtual void     TrackMomentum(Double_t &px, Double_t &py, Double_t &pz, Double_t &etot) const =0;

   /// Return the length in centimeters of the current step (in cm)
   virtual Double_t TrackStep() const =0;

   /// Return the length of the current track from its origin (in cm)
   virtual Double_t TrackLength() const =0;

   /// Return the current time of flight of the track being transported
   virtual Double_t TrackTime() const =0;

   /// Return the energy lost in the current step
   virtual Double_t Edep() const =0;

   //
   // get methods
   // tracking particle
   // static properties
   // ------------------------------------------------
   //

   /// Return the PDG of the particle transported
   virtual Int_t    TrackPid() const =0;

   /// Return the charge of the track currently transported
   virtual Double_t TrackCharge() const =0;

   /// Return the mass of the track currently transported
   virtual Double_t TrackMass() const =0;

   /// Return the total energy of the current track
   virtual Double_t Etot() const =0;

   //
   // get methods - track status
   // ------------------------------------------------
   //

   /// Return true when the track performs the first step
   virtual Bool_t   IsNewTrack() const =0;

   /// Return true if the track is not at the boundary of the current volume
   virtual Bool_t   IsTrackInside() const =0;

   /// Return true if this is the first step of the track in the current volume
   virtual Bool_t   IsTrackEntering() const =0;

   /// Return true if this is the last step of the track in the current volume
   virtual Bool_t   IsTrackExiting() const =0;

   /// Return true if the track is out of the setup
   virtual Bool_t   IsTrackOut() const =0;

   /// Return true if the current particle has disappeared
   /// either because it decayed or because it underwent
   /// an inelastic collision
   virtual Bool_t   IsTrackDisappeared() const =0;

   /// Return true if the track energy has fallen below the threshold
   virtual Bool_t   IsTrackStop() const =0;

   /// Return true if the current particle is alive and will continue to be
   /// transported
   virtual Bool_t   IsTrackAlive() const=0;

   //
   // get methods - secondaries
   // ------------------------------------------------
   //

   /// Return the number of secondary particles generated in the current step
   virtual Int_t    NSecondaries() const=0;

   /// Return the parameters of the secondary track number isec produced
   /// in the current step
   virtual void     GetSecondary(Int_t isec, Int_t& particleId,
                                 TLorentzVector& position, TLorentzVector& momentum) =0;

   /// Return the VMC code of the process that has produced the secondary
   /// particles in the current step
   virtual TMCProcess ProdProcess(Int_t isec) const =0;

   /// Return the array of the VMC code of the processes active in the current
   /// step
   virtual Int_t    StepProcesses(TArrayI &proc) const = 0;

   /// Return the information about the transport order needed by the stack
   virtual Bool_t   SecondariesAreOrdered() const = 0;


   //
   // ------------------------------------------------
   // Control methods
   // ------------------------------------------------
   //

   /// Initialize MC
   virtual void Init() = 0;

   /// Initialize MC physics
   virtual void BuildPhysics() = 0;

   /// Process one event
   /// Deprecated
   virtual void ProcessEvent() = 0;

   /// Process one  run and return true if run has finished successfully,
   /// return false in other cases (run aborted by user)
   virtual Bool_t ProcessRun(Int_t nevent) = 0;

   /// Set switches for lego transport
   virtual void InitLego() = 0;

   /// (In)Activate collecting TGeo tracks
   virtual void SetCollectTracks(Bool_t collectTracks) = 0;

   /// Return the info if collecting tracks is activated
   virtual Bool_t IsCollectTracks() const = 0;

   /// Return the info if multi-threading is supported/activated
   virtual Bool_t IsMT() const { return kFALSE; }

   //
   // ------------------------------------------------
   // Set methods
   // ------------------------------------------------
   //

   /// Set the particle stack
   virtual void SetStack(TVirtualMCStack* stack);

   /// Set the external decayer
   virtual void SetExternalDecayer(TVirtualMCDecayer* decayer);

   /// Set the random number generator
   virtual void SetRandom(TRandom* random);

   /// Set the magnetic field
   virtual void SetMagField(TVirtualMagField* field);

    //
    // ------------------------------------------------
    // Get methods
    // ------------------------------------------------
    //

    /// Return the particle stack
    TVirtualMCStack*   GetStack() const   { return fStack; }

    /// Return the external decayer
    TVirtualMCDecayer* GetDecayer() const { return fDecayer; }

    /// Return the random number generator
    TRandom*           GetRandom() const  { return fRandom; }

    /// Return the magnetic field
    TVirtualMagField*  GetMagField() const  { return fMagField; }

protected:
   TVirtualMCApplication* fApplication; //!< User MC application

private:
   TVirtualMC(const TVirtualMC &mc);
   TVirtualMC & operator=(const TVirtualMC &);

#if !defined(__CINT__)
   static TMCThreadLocal TVirtualMC*  fgMC; ///< Monte Carlo singleton instance
#else
   static                TVirtualMC*  fgMC; ///< Monte Carlo singleton instance
#endif

   TVirtualMCStack*    fStack;   //!< Particles stack
   TVirtualMCDecayer*  fDecayer; //!< External decayer
   TRandom*            fRandom;  //!< Random number generator
   TVirtualMagField*   fMagField;//!< Magnetic field

   ClassDef(TVirtualMC,1)  //Interface to Monte Carlo
};

#define gMC (TVirtualMC::GetMC())

#endif //ROOT_TVirtualMC