diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
index f098c28047c0eb8385690b733db3184c05a11196..527d89ef1d0c1f8aaf5bbc40ee893882ebd2dcc1 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
@@ -62,7 +62,7 @@ G4int nEvts = 100;
 
 G4double thickness = 0.3; //in mm
 
-G4String physList = "FTFP_BERT";
+G4String dimu_xsec_fac = "1000";
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -91,6 +91,12 @@ int main(int argc,char** argv) {
           thickness = G4double(atof(argv[i]));
           ++i;
        }
+       if(G4String(argv[i]) == "--scale")
+       {
+          ++i;
+          dimu_xsec_fac = G4String(argv[i]);
+          ++i;
+       }
 
   }
  
@@ -130,7 +136,7 @@ int main(int argc,char** argv) {
 
   // Increase Gamma2MuMu XSec to 1000x
   
-  UI->ApplyCommand("/testem/phys/SetGammaToMuPairFac  1000");  
+  UI->ApplyCommand("/testem/phys/SetGammaToMuPairFac "+dimu_xsec_fac);  
 
   //Start Simulation
   UI->ApplyCommand("/run/beamOn");
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh
index 45b3517a0b36112036fa32b232cc9d6dc19db575..4fa954952553b7d69242590e46bec0729dbb81eb 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh
@@ -47,18 +47,16 @@ class PrimaryGeneratorMessenger;
 class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
 {
   public:
-    PrimaryGeneratorAction(DetectorConstruction*, G4double, G4int);
+    PrimaryGeneratorAction(DetectorConstruction*, G4double, G4int);    
    ~PrimaryGeneratorAction();
 
   public:
-    void SetRndmBeam(G4double val)  {fRndmBeam = val;}
-    void GeneratePrimaries(G4Event*) override;
+    virtual void GeneratePrimaries(G4Event*);
 
   private:
     G4ParticleGun*             fParticleGun;
     DetectorConstruction*      fDetector;
-    G4double                   fRndmBeam;
-    PrimaryGeneratorMessenger* fGunMessenger;
+    PrimaryGeneratorMessenger* fGunMessenger;     
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh
index 0e08a95f28cd60f38bc9d2d1c00b5cb80c298f6c..0743dd38b3a7ce4bf8cfaa83177d4dca0fbc90a3 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh
@@ -48,13 +48,10 @@ class PrimaryGeneratorMessenger: public G4UImessenger
   public:
     PrimaryGeneratorMessenger(PrimaryGeneratorAction*);
    ~PrimaryGeneratorMessenger();
-
-    void SetNewValue(G4UIcommand*, G4String) override;
-
+    
   private:
     PrimaryGeneratorAction*    fAction;
-    G4UIdirectory*             fGunDir;
-    G4UIcmdWithADouble*        fRndmCmd;
+    G4UIdirectory*             fGunDir;     
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh
index bc835afb8cfb0839d84cb29627f3f2d20481d16c..96913fd0d0103ce58777ae963890838ee375b1c7 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh
@@ -36,6 +36,8 @@
 
 #include "Geant4/G4UserSteppingAction.hh"
 #include "Geant4/globals.hh"
+#include "Geant4/G4ParticleDefinition.hh"
+#include "Geant4/G4ParticleTypes.hh"
 
 class RunAction;
 
@@ -47,12 +49,17 @@ public:
   SteppingAction(RunAction*);
  ~SteppingAction();
 
-  void UserSteppingAction(const G4Step*) override;
-
+  virtual void UserSteppingAction(const G4Step*);
+    
 private:
   RunAction* fRunAction;
   G4double   fMuonMass;
   G4double   fEMass;
+
+  G4MuonPlus*  muplus_def = G4MuonPlus::MuonPlusDefinition();
+  G4Positron* eplus_def  = G4Positron::PositronDefinition();
+
+  int tuple_index    = -1;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
index 617f8faed7d85ac949f1f609d845e10565013be9..9f30f0c97ba2e9c66fe5e8f4854d0e479c012597 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
@@ -14,6 +14,7 @@ parser = argparse.ArgumentParser("DILEPTONPROD")
 parser.add_argument("--nevts", default=10000, help="Number of Events")
 parser.add_argument("--energies", default="[1000]", help="Energies List")
 parser.add_argument("--thickness", default="[0.3]", help="List of thicknesses")
+parser.add_argument("--scale", default="100000", help="Scale DiMuon Conversion Cross Section by factor")
 parser.add_argument("--debug", default=False, action="store_true", help="Debug Mode")
 
 args = parser.parse_args()
@@ -41,10 +42,14 @@ except TypeError:
 
 for energy in energies:
     for t in thickness:
-        _logger.info("%s\n\nRunning Gamma->l+ l- with options:\n\n\tNEvents   :\t%s\n\tEnergy    :\t%s GeV\n\tThickness :\t%s mm\n", 
-         '='*40, args.nevts, energy, t)
+        _logger.info("%s\n\nRunning Gamma->l+ l- Conversion in Aluminium with options:\n"+
+                     "\n\tNumber of Events         :\t%s "+
+                     "\n\tBeam Energy              :\t%s GeV"+
+                     "\n\tTarget Thickness         :\t%s mm"+
+                     "\n\tDiMuon XSec Scale Factor :\t%s\n", 
+         '='*40, args.nevts, energy, t, args.scale)
         _logger.info('='*40)
-        subprocess.check_call("G4GammaToDiLeptonConversionTest.exe --nevts {} --energy {} --thickness {}".format(args.nevts, energy, t), shell=True)
+        subprocess.check_call("G4GammaToDiLeptonConversionTest.exe --nevts {} --energy {} --thickness {} --scale {}".format(args.nevts, energy, t, args.scale), shell=True)
 
 
 _logger.info("Creating Single Combined Plots ROOT File...")
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc
index 47af5bebc905a4012d758aee36f5011423607aa9..9141a2f4366cd26aa2270fc9b31c05a62a961d89 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc
@@ -41,7 +41,6 @@
 #include "Geant4/G4ParticleTable.hh"
 #include "Geant4/G4ParticleDefinition.hh"
 #include "Geant4/G4SystemOfUnits.hh"
-#include "Geant4/Randomize.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -50,7 +49,6 @@ PrimaryGeneratorAction::PrimaryGeneratorAction(
 :G4VUserPrimaryGeneratorAction(),                                              
  fParticleGun(0),
  fDetector(DC),
- fRndmBeam(0),       
  fGunMessenger(0)
 {
   fParticleGun  = new G4ParticleGun(nevts);
@@ -78,18 +76,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 {
   //this function is called at the begining of event
   //
-  G4double x0 = -0.5*(fDetector->GetBoxDepth());
-  G4double y0 = 0.*cm, z0 = 0.*cm;
-    
-  //randomize the beam, if requested.
-  //
-  if (fRndmBeam > 0.) 
-    {
-      G4double rbeam = 0.5*(fDetector->GetBoxDepth())*fRndmBeam;
-      y0 = (2*G4UniformRand()-1.)*rbeam;
-      z0 = (2*G4UniformRand()-1.)*rbeam;
-    }
-  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));  
+  fParticleGun->SetParticlePosition(G4ThreeVector(-0.5*(fDetector->GetBoxDepth()),0.*cm,0.*cm));  
   fParticleGun->GeneratePrimaryVertex(anEvent); 
 }
 
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc
index 1ef2718eb908437bd3aa4df6d0c6e8dd08d10205..e778f31a20b21254b911393ad7e1befd8ebe8267 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc
@@ -42,36 +42,18 @@
 PrimaryGeneratorMessenger::PrimaryGeneratorMessenger(
                                              PrimaryGeneratorAction* Gun)
 :G4UImessenger(),fAction(Gun),
- fGunDir(0),     
- fRndmCmd(0)
+ fGunDir(0)
 { 
   fGunDir = new G4UIdirectory("/testem/gun/");
   fGunDir->SetGuidance("gun control");
 
-  fRndmCmd = new G4UIcmdWithADouble("/testem/gun/rndm",this);
-  fRndmCmd->SetGuidance("random lateral extension on the beam");
-  fRndmCmd->SetGuidance("in fraction of 0.5*sizeYZ");
-  fRndmCmd->SetParameterName("rBeam",false);
-  fRndmCmd->SetRange("rBeam>=0.&&rBeam<=1.");
-  fRndmCmd->AvailableForStates(G4State_Idle);  
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 PrimaryGeneratorMessenger::~PrimaryGeneratorMessenger()
 {
-  delete fRndmCmd;
   delete fGunDir;  
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void PrimaryGeneratorMessenger::SetNewValue(G4UIcommand* command,
-                                               G4String newValue)
-{ 
-  if (command == fRndmCmd)
-   {fAction->SetRndmBeam(fRndmCmd->GetNewDoubleValue(newValue));}   
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
index 532af81c4fef02d5654d58f15e1757f28dca6439..a76c886147c84e74113921dd41a1f5a09b848d7e 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
@@ -47,8 +47,6 @@
 #include "Geant4/G4RunManager.hh"
 #include "Geant4/G4MuonMinus.hh"
 
-#include "Geant4/Randomize.hh"
-
 #include <sstream>
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -107,6 +105,12 @@ void RunAction::BeginOfRunAction(const G4Run* aRun)
   fAnalysis->CreateNtupleDColumn("DiMuon_ThetaPlus");
   fAnalysis->CreateNtupleDColumn("DiMuon_ThetaMinus");
   fAnalysis->CreateNtupleDColumn("DiMuon_InvMass");
+  fAnalysis->CreateNtupleDColumn("DiMuon_PX_Plus");
+  fAnalysis->CreateNtupleDColumn("DiMuon_PY_Plus");
+  fAnalysis->CreateNtupleDColumn("DiMuon_PZ_Plus");
+  fAnalysis->CreateNtupleDColumn("DiMuon_PX_Minus");
+  fAnalysis->CreateNtupleDColumn("DiMuon_PY_Minus");
+  fAnalysis->CreateNtupleDColumn("DiMuon_PZ_Minus");
   fAnalysis->FinishNtuple();
  
   fAnalysis->CreateNtuple("DiElectron", "DiElectron Production");
@@ -115,6 +119,12 @@ void RunAction::BeginOfRunAction(const G4Run* aRun)
   fAnalysis->CreateNtupleDColumn("DiElectron_ThetaPlus");
   fAnalysis->CreateNtupleDColumn("DiElectron_ThetaMinus");
   fAnalysis->CreateNtupleDColumn("DiElectron_InvMass");
+  fAnalysis->CreateNtupleDColumn("DiElectron_PX_Plus");
+  fAnalysis->CreateNtupleDColumn("DiElectron_PY_Plus");
+  fAnalysis->CreateNtupleDColumn("DiElectron_PZ_Plus");
+  fAnalysis->CreateNtupleDColumn("DiElectron_PX_Minus");
+  fAnalysis->CreateNtupleDColumn("DiElectron_PY_Minus");
+  fAnalysis->CreateNtupleDColumn("DiElectron_PZ_Minus");
   fAnalysis->FinishNtuple();
 
   G4cout << "\n----> NTuple file is opened in " << fileName << G4endl;   
@@ -138,12 +148,6 @@ void RunAction::CountProcesses(G4String procName)
 
 void RunAction::EndOfRunAction(const G4Run*)
 {
-  // show Rndm status : not needed
-  //
-  //CLHEP::HepRandom::showEngineStatus();
-
-  //total number of process calls
-  //
   G4double countTot = 0.;
   G4cout << "\n Number of process calls --->";
   for (size_t i=0; i< fProcCounter->size();i++) {
diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
index df2417a07a9a557a44b85dea7c919b7d9ec60481..3ad44d9545de01f1830a41cbadf5340a163c019f 100644
--- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
+++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
@@ -35,7 +35,6 @@
 #include "RunAction.hh"
 #include "Geant4/G4SteppingManager.hh"
 #include "Geant4/G4VProcess.hh"
-#include "Geant4/G4ParticleTypes.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -67,117 +66,64 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep)
  G4ThreeVector PGamma  = PrePoint->GetMomentum();
     
  G4double      Eplus(0), Eminus(0), ETot(0);
- G4ThreeVector Pplus   , Pminus,    PTot;
+ G4ThreeVector Pplus   , Pminus,    PTot; 
  const G4TrackVector* secondary = fpSteppingManager->GetSecondary();
- if(processName == "GammaToMuPair")
- {
-     for (size_t lp=0; lp<(*secondary).size(); lp++) {
-       
-       if ((*secondary)[lp]->GetDefinition()==G4MuonPlus::MuonPlusDefinition()) {
-         Eplus  = (*secondary)[lp]->GetTotalEnergy();
-         Pplus  = (*secondary)[lp]->GetMomentum();
-       } else {
-         Eminus = (*secondary)[lp]->GetTotalEnergy();
-         Pminus = (*secondary)[lp]->GetMomentum();                 
-       }
-
-       ETot += (*secondary)[lp]->GetTotalEnergy();
-       PTot += (*secondary)[lp]->GetMomentum();
-       
-     }
-                   
-     G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma;
-     G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus);
-     G4double GammaPlus = Eplus/fMuonMass;
-     G4double GammaMinus= Eminus/fMuonMass;
-     
-     G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
-
-     if(0.0 == thetaPlus || 0.0 == thetaMinus) {
-       G4cout << "SteppingAction for GammaToMuPair: "
-              << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus
-              << " gamPlus= " << GammaPlus << " gamMinus= " <<  GammaMinus
-              << "  " << thetaPlus *GammaPlus - thetaMinus*GammaMinus << G4endl;
-       return;
-     }
-     analysisManager->FillNtupleDColumn(0, 0, xPlus);
-     analysisManager->FillNtupleDColumn(0, 1, xMinus);
-     analysisManager->FillNtupleDColumn(0, 2, thetaPlus);
-     analysisManager->FillNtupleDColumn(0, 3, thetaMinus);
-     analysisManager->FillNtupleDColumn(0, 4, sqrt(pow(ETot,2)-pow(PTot.getR(),2)));
-     analysisManager->AddNtupleRow(0);
-
-     /*analysisManager->FillH1(1,1./(1.+std::pow(thetaPlus*GammaPlus,2)));
-     analysisManager->FillH1(2,std::log10(thetaPlus*GammaPlus));
-
-     analysisManager->FillH1(3,std::log10(thetaMinus*GammaMinus));
-     analysisManager->FillH1(4,std::log10(std::fabs(thetaPlus *GammaPlus
-                                                  -thetaMinus*GammaMinus)));
-     
-     analysisManager->FillH1(5,xPlus);
-     analysisManager->FillH1(6,xMinus);
-
-     analysisManager->FillH1(18, sqrt(pow(ETot,2)-pow(PTot.getR(),2)));
-     analysisManager->FillH1(19, thetaPlus);
-     analysisManager->FillH1(20, thetaMinus);*/
- }
- else if(processName == "conv")
- {
-     for (size_t lp=0; lp<(*secondary).size(); lp++) {
-       
-       if ((*secondary)[lp]->GetDefinition()==G4Positron::PositronDefinition()) {
-         Eplus  = (*secondary)[lp]->GetTotalEnergy();
-         Pplus  = (*secondary)[lp]->GetMomentum();
-       } else {
-         Eminus = (*secondary)[lp]->GetTotalEnergy();
-         Pminus = (*secondary)[lp]->GetMomentum();                 
-       }
-
-       ETot += (*secondary)[lp]->GetTotalEnergy();
-       PTot += (*secondary)[lp]->GetMomentum();
-       
-     }
-                   
-     G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma;
-     G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus);
-     G4double GammaPlus = Eplus/fEMass;
-     G4double GammaMinus= Eminus/fEMass;
-     
-     G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
-
-     if(0.0 == thetaPlus || 0.0 == thetaMinus) {
-       G4cout << "SteppingAction for GammaToEPair: "
-              << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus
-              << " gamPlus= " << GammaPlus << " gamMinus= " <<  GammaMinus
-              << "  " << thetaPlus *GammaPlus - thetaMinus*GammaMinus << G4endl;
-       return;
-     }
-     analysisManager->FillNtupleDColumn(1, 0, xPlus);
-     analysisManager->FillNtupleDColumn(1, 1, xMinus);
-     analysisManager->FillNtupleDColumn(1, 2, thetaPlus);
-     analysisManager->FillNtupleDColumn(1, 3, thetaMinus);
-     analysisManager->FillNtupleDColumn(1, 4, sqrt(pow(ETot,2)-pow(PTot.getR(),2)));
-     analysisManager->AddNtupleRow(1);
-
-     /*analysisManager->FillH1(21,1./(1.+std::pow(thetaPlus*GammaPlus,2)));
-     analysisManager->FillH1(22,std::log10(thetaPlus*GammaPlus));
-
-     analysisManager->FillH1(23,std::log10(thetaMinus*GammaMinus));
-     analysisManager->FillH1(24,std::log10(std::fabs(thetaPlus *GammaPlus
-                                                  -thetaMinus*GammaMinus)));
-     
-     analysisManager->FillH1(25,xPlus);
-     analysisManager->FillH1(26,xMinus);
+ 
+ G4double particle_mass = -1;
+ tuple_index = -1;
 
-     analysisManager->FillH1(27, sqrt(pow(ETot,2)-pow(PTot.getR(),2)));
-     analysisManager->FillH1(28, thetaPlus);
-     analysisManager->FillH1(29, thetaMinus);*/
+ if(processName == "conv")
+ {  
+     particle_mass = fEMass; 
+     tuple_index = 1;
  }
  else
- {
-      G4cout << "Unrecognised Process!" << G4endl; //Should not get to here!
+ {  
+     particle_mass = fMuonMass; 
+     tuple_index = 0;
+ }
 
+ for (size_t lp=0; lp<(*secondary).size(); lp++) {
+   if ((*secondary)[lp]->GetDefinition()== muplus_def || (*secondary)[lp]->GetDefinition()== eplus_def) {
+     Eplus  = (*secondary)[lp]->GetTotalEnergy();
+     Pplus  = (*secondary)[lp]->GetMomentum();
+   } else {
+     Eminus = (*secondary)[lp]->GetTotalEnergy();
+     Pminus = (*secondary)[lp]->GetMomentum();
+   }
+ 
+   ETot += (*secondary)[lp]->GetTotalEnergy();
+   PTot += (*secondary)[lp]->GetMomentum();
+   
  }
+               
+ G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma;
+ G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus);
+ G4double GammaPlus = Eplus/particle_mass;
+ G4double GammaMinus= Eminus/particle_mass;
+ 
+ G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
+ 
+ if(0.0 == thetaPlus || 0.0 == thetaMinus) {
+   G4cout << "SteppingAction for GammaToMuPair: "
+          << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus
+          << " gamPlus= " << GammaPlus << " gamMinus= " <<  GammaMinus
+          << "  " << thetaPlus *GammaPlus - thetaMinus*GammaMinus << G4endl;
+   return;
+ }
+ analysisManager->FillNtupleDColumn(tuple_index, 0, xPlus);
+ analysisManager->FillNtupleDColumn(tuple_index, 1, xMinus);
+ analysisManager->FillNtupleDColumn(tuple_index, 2, thetaPlus);
+ analysisManager->FillNtupleDColumn(tuple_index, 3, thetaMinus);
+ analysisManager->FillNtupleDColumn(tuple_index, 4, sqrt(pow(ETot,2)-pow(PTot.getR(),2)));
+ analysisManager->FillNtupleDColumn(tuple_index, 5, Pplus.x());
+ analysisManager->FillNtupleDColumn(tuple_index, 6, Pplus.y());
+ analysisManager->FillNtupleDColumn(tuple_index, 7, Pplus.z());
+ analysisManager->FillNtupleDColumn(tuple_index, 8, Pminus.x());
+ analysisManager->FillNtupleDColumn(tuple_index, 9, Pminus.y());
+ analysisManager->FillNtupleDColumn(tuple_index, 10, Pminus.z());
+ analysisManager->AddNtupleRow(tuple_index);
+ 
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......