diff --git a/src/program/example-config.json b/src/program/example-config.json
index 391791b398cfcbbc39e0316fa89a6fe24e660b7b..b0c815f38eab0599726fedfbf8edc350382178f1 100644
--- a/src/program/example-config.json
+++ b/src/program/example-config.json
@@ -1,33 +1,36 @@
 {
   "mc" : {
-    "type": "pythia6m",
-    "write_lund": true,
+    "type": "jlab",
+    "write_lund": "false",
     "generator": {
       "type": "pythia6",
       "target": {
+        "mode": "fixt",
         "A": 1,
         "Z": 1, 
         "nucleon": "p"
       },
       "beam": {
-        "particle": "ELEC",
-        "energy": 11
+        "particle": "e+",
+        "energy": 27.6
       },
-      "Q2": {
-        "min": 0.9,
-        "max": 60
-      },
-      "x": {
-        "min": 0.002,
-        "max": 0.99
-      },
-      "y": {
-        "min": 0.05,
-        "max": 0.95 
-      },
-      "E": {
-        "min": 8.2,
-        "max": 11
+      "limits": {
+        "Q2": {
+          "min": 0.9,
+          "max": 60
+        },
+        "x": {
+          "min": 0.002,
+          "max": 0.99
+        },
+        "y": {
+          "min": 0.05,
+          "max": 0.95 
+        },
+        "E": {
+          "min": 8.2,
+          "max": 11
+        }
       }
     }
   }
diff --git a/src/pythia6m/fmotion/fmotion.F b/src/pythia6m/fmotion/fmotion.F
index e25e424510c28354ac8c76d9d8cee1c4e5529e28..ed850bb2fba174f0c3db4c0c05a29ca39be420dd 100644
--- a/src/pythia6m/fmotion/fmotion.F
+++ b/src/pythia6m/fmotion/fmotion.F
@@ -1,15 +1,19 @@
-      subroutine fmotion (x,y,z)
+C p is the initial 4-momentum (5-component) of the beam without fermi
+C momentum. The routine will pick a proper fermi-momentum in the
+C Rest frame and then boost back to follow the beam.
+      subroutine fmotion (p)
 
       implicit none
 
-      real FM_table(1000),step_size_FM,FMlimit
+      double precision FM_table(1000),step_size_FM,FMlimit
       integer FM_nb
       common/FMvar/FM_table,step_size_FM,FMlimit,FM_nb
-
+C Actual fermi momentum variables
       integer         i
-      real            ThFM,PhiFM,Rd,Kf
-      real            rndm
-      DOUBLE PRECISION x,y,z
+      double precision            ThFM,PhiFM,Rd,Kf
+      double precision            rndm
+      double precision p(5)
+      double precision k(5)
 
       double precision pi
       parameter(pi=3.1415926535898d0)
@@ -28,8 +32,49 @@ ccc Pick a random position in the FM distribution
 ccc Convert i into momentum and smooth the distribution
       Kf =  (i - rndm(0)) * step_size_FM
 
-      x = cos(PhiFM)*sin(ThFM)*Kf
-      y = sin(PhiFM)*sin(ThFM)*Kf
-      z = cos(ThFM)*Kf
+      k(1) = cos(PhiFM)*sin(ThFM)*Kf
+      k(2) = sin(PhiFM)*sin(ThFM)*Kf
+      k(3) = cos(ThFM)*Kf
+      k(4) = sqrt(k(1)**2 + k(2)**2 + k(3)**2 + p(5)**2)
+      k(5) = p(5)
+
+ccc Boost to lab frame
+      call boost_to_lab(k, p)
+      p(1) = k(1)
+      p(2) = k(2)
+      p(3) = k(3)
+      p(4) = k(4)
+      p(5) = k(5)
+
+      end
+
+
+C     Variables we want to boost: p[5]
+C     Lab 4-momentum: p_lab[5]
+C     Sets p to the boosted variables
+      subroutine boost_to_lab (p, p_lab)
+
+      implicit none
+
+C     (px, py, pz, E, m) --> same as in Pythia
+      double precision p(5), p_lab(5)
+      double precision beta(3), btot, ptot
+      double precision gamma, beta_p, kappa
+
+      beta(1) = p_lab(1) / p_lab(4)
+      beta(2) = p_lab(2) / p_lab(4)
+      beta(3) = p_lab(3) / p_lab(4)
+      btot = sqrt(beta(1)**2 + beta(2)**2 + beta(3)**2)
+
+C Only proceed if our boost is meaningfully large
+      if (btot.gt.1d-10) then
+        gamma = 1/sqrt(1D0-btot**2)
+        beta_p = beta(1) * p(1) + beta(2) * p(2) + beta(3) * p(3)
+        kappa = gamma * (gamma * beta_p / (1D0 + gamma) + p(4))
+        p(1) = p(1) + kappa * beta(1)
+        p(2) = p(2) + kappa * beta(2)
+        p(3) = p(3) + kappa * beta(3)
+        p(4) = gamma * (p(4) + beta_p)
+      endif
 
       end
diff --git a/src/pythia6m/fmotion/fmotion.h b/src/pythia6m/fmotion/fmotion.h
index 8f6b765d2b6bec1ff5d7842cc2ffc2c620c5a0a5..0d53bdb8b00811ae9b01af420db5b612f7db4a90 100644
--- a/src/pythia6m/fmotion/fmotion.h
+++ b/src/pythia6m/fmotion/fmotion.h
@@ -5,7 +5,7 @@
 extern "C" {
 #endif
 
-void initfmtable_(int*, int*, float*, int*);
+void initfmtable_(int*, int*, double*, int*);
 #define INIT_FM(z, a, fmco, fmd) initfmtable_(&z, &a, &fmco, &fmd)
 
 #ifdef __cplusplus
diff --git a/src/pythia6m/fmotion/initfmotion.F b/src/pythia6m/fmotion/initfmotion.F
index 914a9c16a3809ff9dd7747f911041e8e1f9e220e..d39e85371f401467de3c5fdf69de9fcc40762c87 100644
--- a/src/pythia6m/fmotion/initfmotion.F
+++ b/src/pythia6m/fmotion/initfmotion.F
@@ -2,7 +2,7 @@
       implicit none
 
 ccccc Table for FM distribution
-      real FM_table(1000),step_size_FM,FMlimit
+      double precision FM_table(1000),step_size_FM,FMlimit
       integer FM_nb
       common/FMvar/FM_table,step_size_FM,FMlimit,FM_nb
 
@@ -12,7 +12,7 @@ ccccc Table for FM distribution
 
       integer i,Z,A
       integer iZ,iA,FMBS,FMD
-      real FMCO
+      double precision FMCO
       double precision rhoCS,mom,momfm,proba,ptot,step
 
 c     TODO force good iZ iA combination
diff --git a/src/pythia6m/gmc/gen_set.h b/src/pythia6m/gmc/gen_set.h
index 0c9e894470bb16795f3687f4bbcaa1ab8f476abd..78cc4baacb2a5321bbfca0cb03459c9be30a8204 100644
--- a/src/pythia6m/gmc/gen_set.h
+++ b/src/pythia6m/gmc/gen_set.h
@@ -9,15 +9,16 @@ struct GeneratorSettings {
   int PhotoType;    // user control switch (MSEL, 9.2)
   int DecayPiK;     // Have LUND decay pi+-, K+- and K0_L
   int DecayLamKs;   // Have LUND decay Lambda, Sigma, Xi, Omega, K0_S
+  int dummy;
   int radgenGenLUT; // generate LUT for pol. rad. corr.
   int radgenUseLUT; // use LUT for pol. rad. corr.
   int enableFM;     // activate or not the FM
   int FMNBins;      // Number of bin for FM discretization
-  float FMCutOff;   // Max FM (in GeV)
-  char PyModel[4];  // l/N (DIS) or gamma/l pN (GVMD,REAL,DIS,RAD)
-  char GenRad[4];   // radiative effects (NO,POL)
-  char FStruct[4];  // structure function parameterization
-  char R[4];        // param. of R=sL/sT (1990, 1998; or 0 for R=0)
+  double FMCutOff;   // Max FM (in GeV)
+  char PyModel[8];  // l/N (DIS) or gamma/l pN (GVMD,REAL,DIS,RAD)
+  char GenRad[8];   // radiative effects (NO,POL)
+  char FStruct[8];  // structure function parameterization
+  char R[8];        // param. of R=sL/sT (1990, 1998; or 0 for R=0)
 };
 extern struct GeneratorSettings common_gen_set_;
 #define glb_gen common_gen_set_
diff --git a/src/pythia6m/gmc/gen_set.inc b/src/pythia6m/gmc/gen_set.inc
index a6d9bc3de59d39c56f4cd70dfe6f9cd69cf3ad3b..f939b6d0f18e8532d9251981022ab4b4f000f04f 100644
--- a/src/pythia6m/gmc/gen_set.inc
+++ b/src/pythia6m/gmc/gen_set.inc
@@ -2,6 +2,7 @@
      +  genSet_PhotoType,
      +  genSet_DecayPiK,
      +  genSet_DecayLamKs,
+     +  genSet_dummy,
      +  genSet_GenLUT,
      +  genSet_UseLUT,
      +  genSet_enableFM,
@@ -15,11 +16,11 @@
 
       integer
      +  genSet_PhotoType,
-     +  genSet_fastMC,
+     +  genSet_dummy,
      +  genSet_enableFM,
      +  genSet_FMNBins
 
-      real*4
+      double precision
      +  genSet_FMCutOff
 
       logical
@@ -28,7 +29,7 @@
      +  genSet_GenLUT,
      +  genSet_UseLUT
 
-      character*4
+      character*8
      +  genSet_PyModel,
      +  genSet_GenRad,
      +  genSet_FStruct,
diff --git a/src/pythia6m/gmc/gmc_event.F b/src/pythia6m/gmc/gmc_event.F
index 4c8949f960cde77423197c0b227ed98d557771d9..dc7f813d65d9654d6d3c16762cfedbcb04a45b3a 100644
--- a/src/pythia6m/gmc/gmc_event.F
+++ b/src/pythia6m/gmc/gmc_event.F
@@ -21,7 +21,8 @@
 
         integer nevgen, i, iscat, nlines
         logical lcut, ok
-        real rlu, pbsgen
+        double precision rlu, pbsgen
+        double precision ptarg(5)
 
 ! Initialize event generation counter
 
@@ -33,8 +34,18 @@
 10      nevgen = nevgen+1
 
 ! recalculate nucleon momentum if required 
+
 100     if (genSet_enableFM.ne.0) then
-          call fmotion(P(2,1), P(2,2), P(2,3))
+          ptarg(1) = 0
+          ptarg(2) = 0
+          ptarg(3) = mcSet_TargetMomentum
+          ptarg(4) = sqrt(mcSet_TargetMomentum**2 + mcSet_TargetMass**2)
+          ptarg(5) = mcSet_TargetMass
+          call fmotion(ptarg)
+          P(2,1) = 0
+          P(2,2) = 0
+          P(2,3) = ptarg(3)
+
         endif
 ! generate a new beam photon if needed
         if (mcSet_BeamParType(1:1).eq.'G') then
diff --git a/src/pythia6m/gmc/gmc_init.F b/src/pythia6m/gmc/gmc_init.F
index 568b31b28a150b32728a4382b7a4c571ae541cac..782684fa17974301222a2a8c656d9b72e0b1ff41 100644
--- a/src/pythia6m/gmc/gmc_init.F
+++ b/src/pythia6m/gmc/gmc_init.F
@@ -17,11 +17,12 @@
 #include "py6radg.inc"
 
         double precision pbeam
+        double precision targmass
         integer j
         character*80    cline
         character*8     cBeam, cTarget, cBeamG
         logical         ok
-        real            rndm
+        double precision rndm
         character*8     cFrame
 
 ! ... force block data modules to be read
@@ -70,15 +71,15 @@ C        mstp(1) = genSet_PyMaxFl
 
         pbeam = mcSet_EneBeam
 
-        if (genSet_enableFM.ne.0) then
+        if (genSet_enableFM.ne.0.or.mcSet_Collider.ne.0) then
           p(1,1) = 0
           p(1,2) = 0
           p(1,3) = pbeam
 
-! use nucleon at rest for initialization
+! use nucleon at central beam momentum for initialization (even with FM)
           p(2,1) = 0
           p(2,2) = 0
-          p(2,3) = 0
+          p(2,3) = mcSet_TargetMomentum
 
           cFrame = '3mom'
 
@@ -119,10 +120,10 @@ C        mstp(1) = genSet_PyMaxFl
 ! ... for realistic bremsstrahlung photons, run with MSTP(142) = 2 to
 ! ... correctly account for the BS cross section
 
-        if ((mcSet_BeamParType.eq.'ELEC').and.
+        if ((mcSet_BeamParType(1:3).eq.'ELE').and.
      &      (genSet_PyModel(1:3).eq.'DIS')) then
           cBeam = 'e-'
-        elseif ((mcSet_BeamParType.eq.'POSI').and. 
+        elseif ((mcSet_BeamParType(1:3).eq.'POS').and. 
      &          (genSet_PyModel(1:3).eq.'DIS')) then
           cBeam = 'e+'
         elseif ((mcSet_BeamParType(1:3).eq.'GAM').and. 
@@ -131,15 +132,15 @@ C        mstp(1) = genSet_PyMaxFl
           mstp(142) = 2
           parp(2) = 1.0
         elseif (genSet_PyModel(1:4).eq.'GVMD') then
-          if (mcSet_BeamParType.eq.'ELEC') then
+          if (mcSet_BeamParType(1:3).eq.'ELE') then
             cBeam = 'gamma/e-'
-          elseif (mcSet_BeamParType.eq.'POSI') then
+          elseif (mcSet_BeamParType(1:3).eq.'POS') then
             cBeam = 'gamma/e+'
           endif
         elseif (genSet_PyModel(1:3).eq.'RAD') then
-          if (mcSet_BeamParType.eq.'ELEC') then
+          if (mcSet_BeamParType(1:3).eq.'ELE') then
             cBeam = 'gamma/e-'
-          elseif (mcSet_BeamParType.eq.'POSI') then
+          elseif (mcSet_BeamParType(1:3).eq.'POS') then
             cBeam = 'gamma/e+'
           endif
 C switch to turn on the rad.corrections in pythia
diff --git a/src/pythia6m/gmc/mcRadCor.h b/src/pythia6m/gmc/mcRadCor.h
index ec99a0b48dfa98f9c30c9a428687717fbd3273db..bf3f1a7dfd606257792a34acca71f91149bc4913 100644
--- a/src/pythia6m/gmc/mcRadCor.h
+++ b/src/pythia6m/gmc/mcRadCor.h
@@ -9,28 +9,29 @@ extern "C" {
         int mcRadCor;
         int ID;
         char cType[4];
-        float xTrue;
-        float yTrue;
-        float nuTrue;
-        float Q2True;
-        float W2True;
-        float thetaBrems;
-        float phiBrems;
-        float sigRad;
-        float sigCor;
-        float sigCorErr;
-        float tailIne;
-        float tailEla;
-        float tailCoh;
-        float vacuum;
-        float vertex;
-        float small;
-        float redFac;
-        float EBrems;
-        float pxTrue;           // true lepton variables
-        float pyTrue;           //
-        float pzTrue;           //
-        float ETrue;            //
+        int dummy;
+        double xTrue;
+        double yTrue;
+        double nuTrue;
+        double Q2True;
+        double W2True;
+        double thetaBrems;
+        double phiBrems;
+        double sigRad;
+        double sigCor;
+        double sigCorErr;
+        double tailIne;
+        double tailEla;
+        double tailCoh;
+        double vacuum;
+        double vertex;
+        double small;
+        double redFac;
+        double EBrems;
+        double pxTrue;           // true lepton variables
+        double pyTrue;           //
+        double pzTrue;           //
+        double ETrue;            //
         int mcRadCor_9999;
     };
     extern RadCorData mcRadCor_;
diff --git a/src/pythia6m/gmc/mcRadCor.inc b/src/pythia6m/gmc/mcRadCor.inc
index 0f6f847a535b227aa08dda5b9aa1ce90cef893d6..885ee9074d0924fd60925afdf416e31d65c5add0 100644
--- a/src/pythia6m/gmc/mcRadCor.inc
+++ b/src/pythia6m/gmc/mcRadCor.inc
@@ -1,8 +1,9 @@
       INTEGER                   mcRadCor,
-     +                          mcRadCor_9999
+     +                          mcRadCor_9999,
+     +                          mcRadCor_dummy
       INTEGER                   mcRadCor_ID
       CHARACTER*4               mcRadCor_cType
-      REAL                      mcRadCor_XTrue,
+      DOUBLE PRECISION          mcRadCor_XTrue,
      +                          mcRadCor_YTrue,
      +                          mcRadCor_NuTrue,
      +                          mcRadCor_Q2True,
@@ -29,6 +30,7 @@
       COMMON /mcRadCor/         mcRadCor,
      +                          mcRadCor_ID,
      +                          mcRadCor_cType,
+     +                          mcRadCor_dummy,
      +                          mcRadCor_XTrue,
      +                          mcRadCor_YTrue,
      +                          mcRadCor_NuTrue,
diff --git a/src/pythia6m/gmc/mc_set.h b/src/pythia6m/gmc/mc_set.h
index 01923255dbe0fb99ff1c7340443dad4907f7a6c0..26de27116a084fa2665eeaf66bf083556ebd336d 100644
--- a/src/pythia6m/gmc/mc_set.h
+++ b/src/pythia6m/gmc/mc_set.h
@@ -6,22 +6,26 @@ extern "C" {
 #endif
 
 struct MCSettings {
-  int RunNo;             // run number
-  int NEvents;           // number of events to generate
-  int TarA;              // target atomic number
-  int TarZ;              // target charge
-  float EneBeam;         // beam energy [GeV]
-  float PBValue;         // Beam Polarization [-1,1]
-  float PTValue;         // Target Polarization [-1,1]
-  float Q2Min;           // min Q2
-  float Q2Max;           // max Q2
-  float xMin;            // min x
-  float xMax;            // max x
-  float yMin;            // min y
-  float yMax;            // max y
-  float EMin;            // min E for a BS beam
-  float EMax;            // max E for a BS beam
-  float RL;              // radiation length (for a photon beam)
+  int RunNo;    // run number
+  int NEvents;  // number of events to generate
+  int TarA;     // target atomic number
+  int TarZ;     // target charge
+  int Collider; // 0: fixed target, 1: collider
+  int dummy;
+  double EneBeam;        // beam energy [GeV]
+  double TargetMomentum; // target momentum [GeV]
+  double TargetMass;     // Target nucleon mass [GeV]
+  double PBValue;        // Beam Polarization [-1,1]
+  double PTValue;        // Target Polarization [-1,1]
+  double Q2Min;          // min Q2
+  double Q2Max;          // max Q2
+  double xMin;           // min x
+  double xMax;           // max x
+  double yMin;           // min y
+  double yMax;           // max y
+  double EMin;           // min E for a BS beam
+  double EMax;           // max E for a BS beam
+  double RL;             // radiation length (for a photon beam)
   char BeamParType[4];   // beam particle type
   char PBeam[4];         // Beam polarization type (L)
   char PTarget[4];       // target polarization type (T,L,DT)
diff --git a/src/pythia6m/gmc/mc_set.inc b/src/pythia6m/gmc/mc_set.inc
index 3a875bf991ce54785461af88ebeb7a5fbae5be6f..3de3e7e244e19df5cd13455e1eea8252e51f5a1d 100644
--- a/src/pythia6m/gmc/mc_set.inc
+++ b/src/pythia6m/gmc/mc_set.inc
@@ -3,7 +3,11 @@
      +  mcSet_NEvents,
      +  mcSet_TarA,
      +  mcSet_TarZ,
+     +  mcSet_Collider,
+     +  mcSet_dummy,
      +  mcSet_EneBeam,
+     +  mcSet_TargetMomentum,
+     +  mcSet_TargetMass,
      +  mcSet_PBValue,
      +  mcSet_PTValue,
      +  mcSet_Q2Min,
@@ -25,10 +29,14 @@
      +  mcSet_RunNo,
      +  mcSet_NEvents,
      +  mcSet_TarA,
-     +  mcSet_TarZ
+     +  mcSet_TarZ,
+     +  mcSet_Collider,
+     +  mcSet_dummy
 
-      real
+      double precision
      +  mcSet_EneBeam,
+     +  mcSet_TargetMomentum,
+     +  mcSet_TargetMass,
      +  mcSet_PBValue,
      +  mcSet_PTValue,
      +  mcSet_Q2Min,
diff --git a/src/pythia6m/gmc/mcevnt.h b/src/pythia6m/gmc/mcevnt.h
index dc6e367c060336da441d447c2214623e4422cb4e..e1c0bbd436dd8e59637cfd10eadf3bcf619347df 100644
--- a/src/pythia6m/gmc/mcevnt.h
+++ b/src/pythia6m/gmc/mcevnt.h
@@ -16,27 +16,28 @@ extern "C" {
         // event generation
         int process;    // PYTHIA process number (MSTI(1))
         int scatIdx;    // PYTHIA index of the scattered lepton
-        float weight;   // event weight
-        float Q2Gen;    // vertex kinematics
-        float nuGen;    //
-        float xGen;     //
-        float yGen;     //
-        float W2Gen;    //
-        float thetaGen; // scattered lepton
-        float phiGen;   //
-        float EGen;     //
-        float pGen;     //
-        float pxGen;    // and its 3-vector
-        float pyGen;    //
-        float pzGen;    //
-        float vxGen;    // vertex position
-        float vyGen;    //
-        float vzGen;    //
+        int dummy;
+        double weight;   // event weight
+        double Q2Gen;    // vertex kinematics
+        double nuGen;    //
+        double xGen;     //
+        double yGen;     //
+        double W2Gen;    //
+        double thetaGen; // scattered lepton
+        double phiGen;   //
+        double EGen;     //
+        double pGen;     //
+        double pxGen;    // and its 3-vector
+        double pyGen;    //
+        double pzGen;    //
+        double vxGen;    // vertex position
+        double vyGen;    //
+        double vzGen;    //
         // timing
-        float timePrev; // time before current event
-        float timeNew;  // time after current event
-        float timeMax;  // maximum event generation time
-        float timeTot;  // total time spent on event generation
+        double timePrev; // time before current event
+        double timeNew;  // time after current event
+        double timeMax;  // maximum event generation time
+        double timeTot;  // total time spent on event generation
     };
     extern struct EventData mcevnt_;
 #define glb_event mcevnt_
diff --git a/src/pythia6m/gmc/mcevnt.inc b/src/pythia6m/gmc/mcevnt.inc
index b480ea87438d276b3d341eb20e05f837f78c742a..bc4e1fb9a7be2f6e99fc7995a7051030428565ae 100644
--- a/src/pythia6m/gmc/mcevnt.inc
+++ b/src/pythia6m/gmc/mcevnt.inc
@@ -15,11 +15,12 @@
      +		mcoldievgen,		! cumulative ievgen BEFORE current ev
      +		mcnerror,		! Monte Carlo errors
      +      mcprocess,      ! PYTHIA process number
-     +      mciscat         ! PYTHIA index of the scattered lepton
+     +      mciscat,        ! PYTHIA index of the scattered lepton
+     +      dummy
 
 ! Event generation 
 
-	real	weight,					! event weight
+	double precision	weight,					! event weight
      +		genq2, gennu, genx, geny, genw2,	! vertex kinematics
      +		genthe, genphi, geneprim, genpprim,	! scattered lepton
      +		genpx, genpy, genpz,			! scat lepton 3-vector
@@ -27,7 +28,7 @@
 
 ! Timing
 
-	real	timold,		! time before current event
+	double precision	timold,		! time before current event
      +		timnew,		! time after current event
      +		timmax,		! maximum time per event we've seen so far
      +		timtot		! total integrated time spent on events
@@ -39,7 +40,7 @@
      +      beamiltype,
      +		mcievent, mcievgen, mcoldievgen, 
      +		mcnerror,
-     +      mcprocess, mciscat,
+     +      mcprocess, mciscat, dummy,
      +		weight,
      +		genq2, gennu, genx, geny, genw2, 
      +		genthe, genphi, geneprim, genpprim,
diff --git a/src/pythia6m/gmc/random.F b/src/pythia6m/gmc/random.F
index fb082975ccb825f5958a9096a3d19b41f3908119..5061c4377a4c82055d723823a92faaac37d71578 100644
--- a/src/pythia6m/gmc/random.F
+++ b/src/pythia6m/gmc/random.F
@@ -35,7 +35,7 @@
 
         integer           nseed1, nseed2, nseq
         integer           i, iseed1, iseed2, iseq, ilux
-        real              r
+        double precision              r
         character*(*)     chopt
         character*1       c1opt
 
@@ -73,11 +73,11 @@
 !-----------------------------------------------------------------
 ! Replace the obsolete CERNLIB RNDM functions
 
-      real function rndm (dummy)
+      double precision function rndm (dummy)
 
         implicit none
 
-        real  dummy, r
+        double precision  dummy, r
 
         call ranlux(r,1)
 
@@ -90,7 +90,7 @@
 
         implicit none
 
-        real dummy, r
+        double precision dummy, r
         integer i
 
         equivalence (r,i)
@@ -109,7 +109,7 @@
         implicit none
 
         integer n
-        real  r(n)
+        double precision  r(n)
 
         call ranlux(r,n)
 
@@ -124,7 +124,7 @@
 
         implicit none
 
-        real  a, b, r(2)
+        double precision  a, b, r(2)
         external nran
 
         call rnormx(r,2,nran)
@@ -144,7 +144,7 @@
         implicit none
 
         integer n
-        real r(n)
+        double precision r(n)
         external nran
 
         call rnormx(r, n, nran)
@@ -155,11 +155,11 @@
 !-----------------------------------------------------------------
 ! Replace the F77 RANF
 
-      real function ranf (dummy)
+      double precision function ranf (dummy)
 
         implicit none
 
-        real   dummy, r
+        double precision   dummy, r
 
         call ranlux(r,1)    
 
@@ -171,12 +171,12 @@
 !-----------------------------------------------------------------
 ! Replace the JETSET random number generator
 
-      real function rlu(idummy)
+      double precision function rlu(idummy)
 
         implicit none
 
         integer idummy
-        real r
+        double precision r
 
         call ranlux(r,1)
 
@@ -193,7 +193,7 @@
         implicit none
 
         integer idummy
-        real r
+        double precision r
 
         call ranlux(r,1)
 
@@ -210,7 +210,7 @@
         implicit none
 
         integer n
-        real r(n)
+        double precision r(n)
 
         call ranlux(r,n)
 
diff --git a/src/pythia6m/gmc/random.h b/src/pythia6m/gmc/random.h
index 03575ab02ed35884d8dd4973f1dcf1afabeeaf84..a32640e306acf8cc2911152920a886d27bd7f4e6 100644
--- a/src/pythia6m/gmc/random.h
+++ b/src/pythia6m/gmc/random.h
@@ -42,35 +42,35 @@ extern "C" {
     rndmq_(&nseed1, &nseed2, &nseq, chopt, strlen(chopt))
 
     // Replace the obsolete CERNLIB RNDM functions
-    float rndm_(float* dummy);
-    int irndm_(float* dummy);
+    double rndm_(double* dummy);
+    int irndm_(double* dummy);
 #define RNDM(dummy) rndm_(&dummy)
 #define IRNDM(dummy) irndm_(&dummy)
 
     // Replace the obsolete CERNLIB NRAN subroutine
-    void nran_(float* array, int* length);
+    void nran_(double* array, int* length);
 #define NRAN(array, length) nran_(array, &length)
 
     // Replace the obsolete CERNLIB RANNOR subroutine
-    void rannor_(float* a, float* b);
+    void rannor_(double* a, double* b);
 #define RANNOR(a, b) rannor_(&a, &b)
 
     // Replace the obsolete CERNLIB RNORML subroutine
-    void rnorml_(float* array, int* length);
+    void rnorml_(double* array, int* length);
 #define RNORML(array, length) rnorml_(array, &length)
 
     // Replace the F77 RANF
-    float ranf_(float* dummy);
+    double ranf_(double* dummy);
 #define RANF(dummy) ranf_(&dummy)
 
     // Replace the JETSET random number generator
-    float rlu_(int* dummy);
+    double rlu_(int* dummy);
     double drlu_(int* dummy);
 #define RLU(dummy) rlu_(&dummy);
 #define DRLU(dummy) drlu_(&dummy);
 
     // Replace the DIVONNE random number generator
-    void ranums_(float* array, int* length);
+    void ranums_(double* array, int* length);
 #define RANUMS(array, length) ranums_(array, &length)
 
 #ifdef __cplusplus
diff --git a/src/pythia6m/interface/pythia6m_generator.cc b/src/pythia6m/interface/pythia6m_generator.cc
index 93cef7c76be5d1a54d0c9637c8ee1876333ece03..6ae1a78be0c4c41417c20966aa221b777fb8e3b7 100644
--- a/src/pythia6m/interface/pythia6m_generator.cc
+++ b/src/pythia6m/interface/pythia6m_generator.cc
@@ -1,6 +1,5 @@
 #include "pythia6m_generator.hh"
 
-
 #include <cstring>
 #include <pythia6m/core/exception.hh>
 #include <pythia6m/core/logger.hh>
@@ -19,7 +18,8 @@ MCSettings glb_mc;
 EventData glb_event;
 
 namespace pythia6m {
-pythia6m_generator::pythia6m_generator(const configuration& conf, const string_path& path)
+pythia6m_generator::pythia6m_generator(const configuration& conf,
+                                       const string_path& path)
     : conf_(conf, path) {
   LOG_INFO("PYTHIA6M", "Initializing the generator");
   // configuration: first GMC, then the generator (depends on e.g. beam settings
@@ -38,7 +38,7 @@ pythia6m_generator::pythia6m_generator(const configuration& conf, const string_p
   init_gmc();
 }
 
-pythia6m_generator::~pythia6m_generator() { 
+pythia6m_generator::~pythia6m_generator() {
   GMC_RUNSTAT();
   LOG_INFO("PYTHIA6M",
            "TOTAL Generated events: " + std::to_string(glb_event.nEvGen));
@@ -87,7 +87,7 @@ void pythia6m_generator::conf_generator() const {
     sprintf(glb_gen.PyModel, "REAL");
   } else {
     sprintf(glb_gen.PyModel, "RAD");
-    //sprintf(glb_gen.PyModel, "GVMD");
+    // sprintf(glb_gen.PyModel, "GVMD");
   }
 }
 // =============================================================================
@@ -100,15 +100,34 @@ void pythia6m_generator::conf_gmc() const {
   glb_mc.TarA = conf_.get<int>("target/Z");
   glb_mc.TargetNucleon[0] = conf_.get<char>("target/nucleon");
   glb_mc.EneBeam = conf_.get<double>("beam/energy");
-
+  if (glb_mc.TargetNucleon[0] == 'p') {
+    glb_mc.TargetMass = .93827208816;
+  } else if (glb_mc.TargetNucleon[0] = 'n') {
+    glb_mc.TargetMass = .93956542052;
+  } else {
+    LOG_ERROR("conf_gmc", "Valid target nucleons are 'p' or 'n'");
+    throw conf_.value_error("target/nucleon");
+  }
+  if (conf_.get<std::string>("target/mode") == "fixt") {
+    glb_mc.Collider = false;
+    glb_mc.TargetMomentum = 0.;
+  } else if (conf_.get<std::string>("target/mode") == "collider") {
+    glb_mc.Collider = false;
+    // energy per nucleon
+    const double E = conf_.get<double>("target/energy") / glb_mc.TarA;
+    glb_mc.TargetMomentum = sqrt(E * E - glb_mc.TargetMass * glb_mc.TargetMass);
+  } else {
+    LOG_ERROR("conf_gmc", "Valid target modes are 'fixt' or 'collider'");
+    throw conf_.value_error("target/nucleon");
+  }
   // beam particle
   std::string beam = conf_.get<std::string>("beam/particle");
-  if (beam != "ELEC" && beam != "POSI" && beam != "GAM") {
-    LOG_ERROR("conf_gmc", "Valid beam particles are 'ELEC', 'POSI' or 'GAM'");
+  if (beam != "e-" && beam != "e+" && beam != "gamma") {
+    LOG_ERROR("conf_gmc", "Valid beam particles are 'e-', 'e+' or 'gamma'");
     throw conf_.value_error("beam/particle");
   }
-  sprintf(glb_mc.BeamParType, "%s", beam.c_str());
-  if (beam == "GAM") {
+  if (beam == "gamma") {
+    sprintf(glb_mc.BeamParType, "GAM");
     // real photon setup
     LOG_INFO("pythia6", "Setting up MC for photoproduction");
     glb_mc.Q2Min = 0;
@@ -117,18 +136,23 @@ void pythia6m_generator::conf_gmc() const {
     glb_mc.xMax = 1.;
     glb_mc.yMin = 0.;
     glb_mc.yMax = 1.;
-    glb_mc.EMin = conf_.get<double>("E/min");
-    glb_mc.EMax = fmin(conf_.get<double>("E/max"), glb_mc.EneBeam);
-    glb_mc.RL = conf_.get<double>("beam/RL");
+    glb_mc.EMin = conf_.get<double>("limits/E/min");
+    glb_mc.EMax = fmin(conf_.get<double>("limits/E/max"), glb_mc.EneBeam);
+    glb_mc.RL = conf_.get<double>("limits/beam/RL");
   } else {
     // electron or positron
+    if (beam == "e-") {
+      sprintf(glb_mc.BeamParType, "ELE");
+    } else {
+      sprintf(glb_mc.BeamParType, "POS");
+    }
     LOG_INFO("pythia6", "Setting up MC for electroproduction");
-    glb_mc.Q2Min = conf_.get<double>("Q2/min");
-    glb_mc.Q2Max = conf_.get<double>("Q2/max");
-    glb_mc.xMin = conf_.get<double>("x/min");
-    glb_mc.xMax = conf_.get<double>("x/max");
-    glb_mc.yMin = conf_.get<double>("y/min");
-    glb_mc.yMax = conf_.get<double>("y/max");
+    glb_mc.Q2Min = conf_.get<double>("limits/Q2/min");
+    glb_mc.Q2Max = conf_.get<double>("limits/Q2/max");
+    glb_mc.xMin = conf_.get<double>("limits/x/min");
+    glb_mc.xMax = conf_.get<double>("limits/x/max");
+    glb_mc.yMin = conf_.get<double>("limits/y/min");
+    glb_mc.yMax = conf_.get<double>("limits/y/max");
     glb_mc.EMin = 0.;
     glb_mc.EMax = 10000000.;
     glb_mc.RL = 0;
@@ -148,9 +172,9 @@ void pythia6m_generator::conf_gmc() const {
 void pythia6m_generator::init_event() const {
   LOG_INFO("init", "Initializing the event common blocks");
   // beam particle
-  if (!strncmp(glb_mc.BeamParType, "POSI", 4)) {
+  if (!strncmp(glb_mc.BeamParType, "POS", 3)) {
     glb_event.beamLType = -11;
-  } else if (!strncmp(glb_mc.BeamParType, "ELEC", 4)) {
+  } else if (!strncmp(glb_mc.BeamParType, "ELE", 3)) {
     glb_event.beamLType = 11;
   } else if (!strncmp(glb_mc.BeamParType, "GAM", 3)) {
     glb_event.beamLType = 22;
@@ -209,4 +233,4 @@ void pythia6m_generator::init_random() const {
   RNDMQ(dummy1, dummy2, glb_mc.RunNo, command);
 }
 
-} // ns pythia6m
+} // namespace pythia6m
diff --git a/src/pythia6m/pythia6/pyr.F b/src/pythia6m/pythia6/pyr.F
index 4d5d9423e7c8ed417f2bec37850bd16e01da8535..b7f2a4432e26b013d425d3ac49573c6b6ce5bf07 100644
--- a/src/pythia6m/pythia6/pyr.F
+++ b/src/pythia6m/pythia6/pyr.F
@@ -21,7 +21,7 @@
 	implicit none
 
 	integer	idummy
-	real	r
+	double precision	r
 
 	call ranlux(r,1)
 	pyr = r
diff --git a/src/pythia6m/radgen/F b/src/pythia6m/radgen/F
new file mode 100644
index 0000000000000000000000000000000000000000..0cfb4f4b1a84f546959278846be745a7ee406ade
--- /dev/null
+++ b/src/pythia6m/radgen/F
@@ -0,0 +1 @@
+[?1049h[?1h=[?25l"radgen" [New File]                                                                                           1 ¬                                                                                                       ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                           ~                                                                                                            NORMAL   develop  radgen                                               [unix]  100% ☰    0/1  :  1                    [?25h[?25l:q
[?1l>[?25h[?1049l
\ No newline at end of file
diff --git a/src/pythia6m/radgen/density.inc b/src/pythia6m/radgen/density.inc
index 91194ee3eaa9a32e17fb58247fcab4b52096f22e..fdca302f621231fe10693e834501f2cbe79176a9 100644
--- a/src/pythia6m/radgen/density.inc
+++ b/src/pythia6m/radgen/density.inc
@@ -5,9 +5,9 @@
       parameter(ntpho=44)
       parameter(nt=ntpho)
 
-      real*4 denstk,width
-      real*4 densdis,widdis
-      real*4 denspho,widpho
+      real*8 denstk,width
+      real*8 densdis,widdis
+      real*8 denspho,widpho
 * make these arrays large enough for all cases
       common/density/ntx,nty,denstk(nt,nt,245,3),width(nt,nt,7,3)
      &              ,densdis(ntdis,ntdis,245,3),widdis(ntdis,ntdis,7,3)
diff --git a/src/pythia6m/radgen/radgen.F b/src/pythia6m/radgen/radgen.F
index ae8e38c0c7b9c4969b41ca42b92a045e6441a4c2..6ea32f96d2b647267ff5e3995e0ee55b57283221 100644
--- a/src/pythia6m/radgen/radgen.F
+++ b/src/pythia6m/radgen/radgen.F
@@ -45,7 +45,7 @@ C **********************************************************************
 *        q2,nu : uncorrected kinematics of scattered lepton
 *        phil  : azimuthal angle of scattered lepton
 *  OUTPUT :
-*	 PhRAD(4) : real rad. photon 4 vector
+*	 PhRAD(4) : double precision rad. photon 4 vector
 *	 q2tr	  : true Q2
 *	 Utr	  : true U
 ***********************************************************
@@ -57,11 +57,11 @@ C **********************************************************************
 #include "radgenkeys.inc"
 
 C eca 041106 to allow that g* theta is really small
-C      real PhRAD(4)
+C      double precision PhRAD(4)
 
       double precision phrad(4)
-      REAL q2tr,anutr
-      REAL e1,q2l,nul,phil,xs,ys,weight
+      double precision q2tr,anutr
+      double precision e1,q2l,nul,phil,xs,ys,weight
       INTEGER itagen
 *
 
@@ -87,12 +87,12 @@ C      real PhRAD(4)
 
       if(itagen.ne.0.and.abs(ixytest).ne.2)then
 	iphi=0
-	call mpolrad(e1,ys,xs,1.,plrun,pnrun,itagen)
+	call mpolrad(e1,ys,xs,1.D0,plrun,pnrun,itagen)
 	if(ixytest.eq.1)ixytest=-2
       endif
       call radgam_pol(e1,ys,xs,phil,ixytest,itagen,q2tr,anutr)
       
-* set kinematics of possibly radiated real gamma
+* set kinematics of possibly radiated double precision gamma
       phrad(1)=dplabg(1)
       phrad(2)=dplabg(2)
       phrad(3)=dplabg(3)
@@ -165,6 +165,7 @@ C **********************************************************************
 
       subroutine radgam_pol(erad,yrad,xrad,genphi,
      +     ixytest,itagen,q2tr,anutr)
+      implicit real*8(a-c,e-h,o-z)
 
 #include "double.inc"
 #include "phiout.inc"
@@ -176,7 +177,7 @@ C **********************************************************************
 #include "density.inc"
 
       real*8 podinl
-      real temp                               ! Temp variable for a denom
+      double precision temp                               ! Temp variable for a denom
 
 *
 *     selection of mass dmj and mass index iadmj
@@ -326,7 +327,7 @@ C	 pause
         deg=dom
         dthg=dtk
         dphig=dphk
-        sigma_total=sngl(tbor+tine+tpro+tnuc)
+        sigma_total=(tbor+tine+tpro+tnuc)
         
         q2tr=y+rrout*taout
         anutr=sx/ap-dom
@@ -451,7 +452,7 @@ C **********************************************************************
 #include "density.inc"
 #include "py6radg.inc"
 
-      real e1curr,yscurr,xscurr,uncurr,plcurr,pncurr
+      double precision e1curr,yscurr,xscurr,uncurr,plcurr,pncurr
 
       dimension tai(5),si(2,3),si0(2,3),tls(2,3,4)
 
@@ -473,7 +474,7 @@ c        negative polarization is twice as large as positive one
       call conk2(e1,xs,ys,1)
 c
 c
-c delta is factorizing part of virtual and real leptonic bremsstrahlung
+c delta is factorizing part of virtual and double precision leptonic bremsstrahlung
 c
       call deltas(delta,delinf,tr)
 
@@ -826,7 +827,7 @@ C **********************************************************************
 
       subroutine deltas(delta,delinf,tr)
 c
-c delta is factorizing part of virtual and real leptonic bremsstrahlung
+c delta is factorizing part of virtual and double precision leptonic bremsstrahlung
 c
       implicit real*8(a-h,o-z)
 
@@ -1384,7 +1385,7 @@ c
 #include "radgenkeys.inc"
 #include "py6radg.inc"
 
-      real forgetp,forintp
+      double precision forgetp,forintp
       external forgetp,forintp
 
       dimension sfm(8)
@@ -1557,7 +1558,7 @@ c
           b4=0d0
         elseif((ire.eq.14).or.(ire.eq.20).or.(ire.eq.84).or.
      +          (ire.eq.131))then
-          ff=forgetp(sngl(t))
+          ff=forgetp(t)
 c         if (t.lt.0.002) then 
 c           ff=forgetp(sngl(t))
 c           print *,t,forgetp(sngl(t))
@@ -1803,13 +1804,13 @@ C DGAM2=(RP/SQRT(3/2))**2 = 0.66666666            for He4
       end
 
 
-      REAL FUNCTION FORGETp(Q2)
+      double precision FUNCTION FORGETp(Q2)
 C     GET FORM FACTOR BY INTERPOLATION OUT OF TABLE FFNUC
 
       implicit none
 #include "tailcom.inc"
       integer iq
-      real q2
+      double precision q2
 
       FORGETp=0.
 C     LOWER BIN OF Q2
@@ -1827,9 +1828,9 @@ C     of the charge distribution
 
       implicit none
 #include "tailcom.inc"
-      real q2max,rmin,rmax,eps,hbc
-      real q2
-      real forintp,gauss
+      double precision q2max,rmin,rmax,eps,hbc
+      double precision q2
+      double precision forintp,gauss
       external forintp
       integer nqbin1,iq
 
@@ -1858,7 +1859,7 @@ C
       END
 
 
-        REAL FUNCTION FORINTp(R)
+      double precision FUNCTION FORINTp(R)
 C       INTEGRAND OF FOURIER TRANSFORMATION INTEGRAL
 C       THREE-DIMENSIONAL INTEGRATION IS REDUCED TO
 C       AN INTEGRATION OVER THE RADIUS BY MULTIPLYING
@@ -1866,7 +1867,7 @@ C       WITH A BESSEL FKT.
 
       implicit none
 #include "tailcom.inc"
-      real r,qr,forrhop
+      double precision r,qr,forrhop
 
         IF (QFOR.GT.1.E-5) THEN
            QR=QFOR*R
@@ -1879,7 +1880,7 @@ c        print *,'forintp=',forintp
         END
 
 
-        REAL FUNCTION FORRHOp(R)
+        double precision FUNCTION FORRHOp(R)
 C       charge distribution as fkt. of radius
 C     generalized 2-parameter fermi for not too light nuclei
 C     Ref: T.de Forest and J.D.Walecka, Adv.in Phys.15(1966)1  -- page 21
@@ -1889,7 +1890,7 @@ C         z = 2.4 / (4*ln(3))
       implicit none
 #include "tailcom.inc"
 #include "cmpcom.inc"
-       real r,zpara,cpara,wpara
+       double precision r,zpara,cpara,wpara
 
        if (ire.eq.14) then
 C - 3 parameter Fermi charge distr of 14N
@@ -2193,14 +2194,17 @@ C **********************************************************************
 
       subroutine itafun(ys,xs,pl,pn,ixytest,itagen)
 
+      implicit real*8 (a-h,o-z)
 #include "cmpcom.inc"
 #include "radgen.inc"
 #include "density.inc"
 #include "xytabcom.inc"
 
+      integer ixytest,itagen
+
       dimension xym(2),na(2),a(2*nt),wrk(nt,nt)
       dimension xx(nt),work(nt),w(nt),d(nt),xwork1(nt),xwork2(nt)
-	real prevy
+	double precision prevy
 
       do nny=nty,1,-1
         if(ys.gt.y(nny))goto 10
@@ -2361,6 +2365,7 @@ C
 C **********************************************************************
 
       subroutine xytabl(tname,e1,plrun,pnrun,ixytest,ire)
+      implicit real*8(a-h,o-z)
 
 #include "cmpcom.inc"
 #include "radgen.inc"
@@ -2410,7 +2415,7 @@ C **********************************************************************
             ylim=(bmc2-bmp**2)/(2.*bmp*e1*(1d0-x(ix)))
             ys=y(iy)
             if(ys.lt.ylim ) ys=ylim
-            call mpolrad(e1,ys,x(ix),1.,plrun,pnrun,-1)
+            call mpolrad(e1,ys,x(ix),1.D0,plrun,pnrun,-1)
             tbor_u(ix,iy)=tbor
             sig1g_u(ix,iy)=sig1g
             sigrad_u(ix,iy)=sigrad
diff --git a/src/pythia6m/radgen/radgen.inc b/src/pythia6m/radgen/radgen.inc
index 31feab1099e212c9dafb2e3634678db799e765da..996908939507e88300e7e876c1ab54d667ec4932 100644
--- a/src/pythia6m/radgen/radgen.inc
+++ b/src/pythia6m/radgen/radgen.inc
@@ -1,4 +1,4 @@
-      real radgen_xmin, radgen_xmax, radgen_ymin, radgen_ymax
+      double precision radgen_xmin, radgen_xmax, radgen_ymin, radgen_ymax
 
       parameter (radgen_xmin=1.0e-09)
       parameter (radgen_xmax=0.99)
@@ -10,9 +10,9 @@
      +     ,dsts,dcts
      +     ,taout,rrout,dsitkm,dcvtkm,ddetkm,dsigmr,drcurr,ddeler
 
-      integer ntk,nrr,itkcur,iphi,ndxtkm
+      integer ntk,nrr,itkcur,iphi,ndxtkm,rcdummy
 
-      real sigradu, sigradp, sig1gu, sig1gp
+      double precision sigradu, sigradp, sig1gu, sig1gp
 
       common /rgencom/sigrad,tine,tnuc,tpro,tbor,demin
      +     ,sig1g,sigcor,vac,vertex,small,redfac
@@ -20,7 +20,7 @@
      +     ,phipoi,taout,rrout
      +     ,dsitkm(400,3),dcvtkm(400,3),ddetkm(400,3)
      +     ,dsigmr(200,400),drcurr(200,400),ddeler(200,400)
-     +     ,ntk,nrr ,itkcur,iphi,ndxtkm(3)
+     +     ,ntk,nrr ,itkcur,iphi,ndxtkm(3),rcdummy
      +     ,sigradu, sigradp, sig1gu, sig1gp
 
       
diff --git a/src/pythia6m/radgen/radgen_event.F b/src/pythia6m/radgen/radgen_event.F
index 2244b9d17b12864c2a4b7aca68c5ab31f85e728e..98e1146a5ebcbf27cf9cb8c26d666282f713b54d 100644
--- a/src/pythia6m/radgen/radgen_event.F
+++ b/src/pythia6m/radgen/radgen_event.F
@@ -16,10 +16,10 @@ C#include "partap.inc"
 
       integer iphot
 C eca 041106 to allow that g* theta is really small
-C      real PhRAD(4)
+C      double precision PhRAD(4)
       double precision phrad(4)
-      real q2true,nutrue,radweight
-      real nu,q2,phi
+      double precision q2true,nutrue,radweight
+      double precision nu,q2,phi
 
 ! calculate radiative corrections
       nu=gennu
@@ -44,23 +44,23 @@ C      real PhRAD(4)
       mcRadCor_W2True=amp2 - q2true + 2.*amp*nutrue
 *     print *,amp,amp2,mcRadCor_XTrue,nutrue,q2true,mcRadCor_W2True
 
-! ... kinematics of real photon
+! ... kinematics of double precision photon
 
-      mcRadCor_EBrems=sngl(phrad(4))
+      mcRadCor_EBrems=phrad(4)
       mcRadCor_ThetaBrems=0.
       mcRadCor_PhiBrems=0.
 
       if(phrad(4).gt.0.) then
-       mcRadCor_ThetaBrems = sngl(acos(phrad(3)/phrad(4)))
+       mcRadCor_ThetaBrems = acos(phrad(3)/phrad(4))
       endif
 
       if (.not.(phrad(1).eq.0..and.phrad(2).eq.0.)) then
-        mcRadCor_PhiBrems = sngl(atan2(phrad(2),phrad(1)))
+        mcRadCor_PhiBrems = atan2(phrad(2),phrad(1))
         if (mcRadCor_PhiBrems.lt.0.)
      +       mcRadCor_PhiBrems = mcRadCor_PhiBrems + twopi
       endif
 
-c...if we would like to have the TSAI system angles for the real gamma than
+c...if we would like to have the TSAI system angles for the double precision gamma than
 c     mcRadCor_ThetaBrems = dthg
 c     mcRadCor_PhiBrems = dphig
 
diff --git a/src/pythia6m/radgen/radgen_init.F b/src/pythia6m/radgen/radgen_init.F
index d9337258a2a3ff8a5e070dd919a02f5cbbd87b5d..bc41947795c7225ad221e770252b6dbff85d07bb 100644
--- a/src/pythia6m/radgen/radgen_init.F
+++ b/src/pythia6m/radgen/radgen_init.F
@@ -10,7 +10,7 @@
 #include "radgen.inc"
 #include "radgenkeys.inc"
 
-      real ebeam,polbeam,poltarget
+      double precision ebeam,polbeam,poltarget
       logical bUseLUT,bGenLUT
 C      INTEGER hermesfile
       CHARACTER*256 tname,fname
diff --git a/src/pythia6m/radgen/radgenkeys.inc b/src/pythia6m/radgen/radgenkeys.inc
index 9fd8f5a0954f9a72b825107164396ba9e384c926..ccbc17f448ab2fe3b5731e7ed6874c8dc8adcaa3 100644
--- a/src/pythia6m/radgen/radgenkeys.inc
+++ b/src/pythia6m/radgen/radgenkeys.inc
@@ -1,4 +1,4 @@
 
       integer ixytest, kill_elas_res
-      real plrun,pnrun
+      double precision plrun,pnrun
       common/radgenkeys/plrun,pnrun,ixytest,kill_elas_res
diff --git a/src/pythia6m/radgen/tailcom.inc b/src/pythia6m/radgen/tailcom.inc
index a2855a304353e70d9d90d78a67be8d8c43c20b7e..142861bf10687de209e00ae02d4f28b484fa44c6 100644
--- a/src/pythia6m/radgen/tailcom.inc
+++ b/src/pythia6m/radgen/tailcom.inc
@@ -1,4 +1,4 @@
-      real    pl,pn
+      double precision    pl,pn
       integer ita,isf1,isf2,isf3,ire
       double precision qfor,q2bin,ffnuc,un,qn
       integer          nqbin
diff --git a/src/pythia6m/radgen/xytabcom.inc b/src/pythia6m/radgen/xytabcom.inc
index a182ae5575952a19925ab0757341f5d41b8fffaa..7ff13f88821e4a653e06ba1729ea9c9a3d762675 100644
--- a/src/pythia6m/radgen/xytabcom.inc
+++ b/src/pythia6m/radgen/xytabcom.inc
@@ -1,4 +1,4 @@
-      real x,y
+      double precision x,y
      +    ,sig1g_u,sigrad_u
      +    ,tbor_u,tine_u,tnuc_u,tpro_u
      +    ,sig1g_p,sigrad_p
diff --git a/src/pythia6m/radiator/pbsgen.F b/src/pythia6m/radiator/pbsgen.F
index 9a3f51b11b0602d273b533fc357a7d6e8439e7d1..8bbbda78d06cfc0389f3d2f61965f45c4eb4ef1e 100644
--- a/src/pythia6m/radiator/pbsgen.F
+++ b/src/pythia6m/radiator/pbsgen.F
@@ -7,11 +7,11 @@ C...  Does not actually evaluate the bremsstrahlung spectrum, this has
 C...  to be done through a PBSINT call. This secondary call is performed
 C...  during the PYEVWT step (when MSTP(142) = 2)
  
-      REAL FUNCTION PBSGEN(PBSEMIN, PBSEMAX)
+      DOUBLE PRECISION FUNCTION PBSGEN(PBSEMIN, PBSEMAX)
  
 C...Double precision and integer declarations.
       IMPLICIT NONE
-      REAL PBSEMIN, PBSEMAX
+      DOUBLE PRECISION PBSEMIN, PBSEMAX
 
       CALL RANLUX(PBSGEN,1)
       PBSGEN = (PBSEMIN + PBSGEN * (PBSEMAX - PBSEMIN))