From 040457f56e65286a970ecdfd8691c29e6aa7f342 Mon Sep 17 00:00:00 2001
From: Thomas <thomas@pegasus.(none)>
Date: Mon, 11 Aug 2014 22:42:46 +0200
Subject: [PATCH] Begin to adapt for Geant4.10

---
 include/F04GlobalField.hh             |   2 +-
 include/musrPrimaryGeneratorAction.hh |   2 +-
 include/musrRootOutput.hh             |  32 +--
 include/musrRunAction.hh              |   2 +-
 src/musrDetectorConstruction.cc       | 350 +++++++++++++-------------
 src/musrUniformField.cc               |   4 +-
 6 files changed, 196 insertions(+), 196 deletions(-)

diff --git a/include/F04GlobalField.hh b/include/F04GlobalField.hh
index c0d4704..a153ef6 100644
--- a/include/F04GlobalField.hh
+++ b/include/F04GlobalField.hh
@@ -118,7 +118,7 @@ public:
   void SetStepper();
 
   /// Set the minimum step length
-  void SetMinStep(G4double s) { minStep = s;  G4cout<<"F04GlobalField::SetMinStep: minStep set to "<<minStep/mm<<" mm."<<G4endl; }
+  void SetMinStep(G4double s) { minStep = s;  G4cout<<"F04GlobalField::SetMinStep: minStep set to "<<minStep/CLHEP::mm<<" mm."<<G4endl; }
 
   /// Set the delta chord length
   void SetDeltaChord(G4double s) { deltaChord = s; }
diff --git a/include/musrPrimaryGeneratorAction.hh b/include/musrPrimaryGeneratorAction.hh
index 6bcf77a..62c87f0 100644
--- a/include/musrPrimaryGeneratorAction.hh
+++ b/include/musrPrimaryGeneratorAction.hh
@@ -82,7 +82,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
   //    void SetOrReadTheRandomNumberSeeds(G4int eventID);
     void SetOrReadTheRandomNumberSeeds(G4Event* anEvent);
     void SetTurtleMomentumBite (G4ThreeVector smearingParam)
-                               {turtleMomentumBite=true; turtleMomentumP0=smearingParam[0]*MeV; turtleSmearingFactor=smearingParam[1]*0.01;}
+                               {turtleMomentumBite=true; turtleMomentumP0=smearingParam[0]*CLHEP::MeV; turtleSmearingFactor=smearingParam[1]*0.01;}
     void SetTurtleMomentumScalingFactor(G4double momentumScaling) {turtleMomentumScalingFactor=momentumScaling;}
     void SetPrimaryParticule(G4String particleName);
     static G4String GetPrimaryName();                
diff --git a/include/musrRootOutput.hh b/include/musrRootOutput.hh
index aea30bb..e2dbcf9 100644
--- a/include/musrRootOutput.hh
+++ b/include/musrRootOutput.hh
@@ -54,20 +54,20 @@ class musrRootOutput  {
 
   // Getting variables (just for debugging)
     G4double GetDecayPositionZ() {return muDecayPosZ;};
-    G4double GetDecayTime()      {return muDecayTime*microsecond;};
-    G4double GetTimeInTarget()   {return muTargetTime*microsecond;};
+    G4double GetDecayTime()      {return muDecayTime*CLHEP::microsecond;};
+    G4double GetTimeInTarget()   {return muTargetTime*CLHEP::microsecond;};
 
   // Setting variables common to the whole event:
     void SetRunID          (G4int id) {runID = id;};
     void SetEventID        (G4int id) {eventID = id;};
-    void SetTimeToNextEvent(G4double deltaT) {timeToNextEvent = deltaT/microsecond;}
+    void SetTimeToNextEvent(G4double deltaT) {timeToNextEvent = deltaT/CLHEP::microsecond;}
     void SetDecayDetectorID (std::string detectorName) {muDecayDetID = SensDetectorMapping[detectorName];};
     Int_t GetDecayDetectorID() {return  muDecayDetID;};
-    void SetBField       (G4double F[6]) {B_t[0]=F[0]/tesla; B_t[1]=F[1]/tesla; B_t[2]=F[2]/tesla; 
-                                          B_t[3]=F[3]/tesla; B_t[4]=F[4]/tesla; B_t[5]=F[5]/tesla;};
+    void SetBField       (G4double F[6]) {B_t[0]=F[0]/CLHEP::tesla; B_t[1]=F[1]/CLHEP::tesla; B_t[2]=F[2]/CLHEP::tesla; 
+                                          B_t[3]=F[3]/CLHEP::tesla; B_t[4]=F[4]/CLHEP::tesla; B_t[5]=F[5]/CLHEP::tesla;};
     void SetDecayPolarisation (G4ThreeVector pol) {muDecayPolX=pol.x(); muDecayPolY=pol.y(); muDecayPolZ=pol.z();};
-    void SetDecayPosition (G4ThreeVector pos) {muDecayPosX=pos.x()/mm; muDecayPosY=pos.y()/mm; 
-                                               muDecayPosZ=pos.z()/mm;};
+    void SetDecayPosition (G4ThreeVector pos) {muDecayPosX=pos.x()/CLHEP::mm; muDecayPosY=pos.y()/CLHEP::mm; 
+                                               muDecayPosZ=pos.z()/CLHEP::mm;};
     void SetEventWeight  (G4double w) {weight *= w;}
     void SetDetectorInfo (G4int nDetectors, G4int ID, G4int particleID, G4double edep, 
                           G4double edep_el, G4double edep_pos, 
@@ -93,7 +93,7 @@ class musrRootOutput  {
 
     void SetInitialMuonParameters(G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz, 
 				  G4double xpolaris, G4double ypolaris, G4double zpolaris, G4double particleTime) {
-      muIniTime=particleTime/microsecond;
+      muIniTime=particleTime/CLHEP::microsecond;
       muIniPosX=x;  muIniPosY=y;  muIniPosZ=z;
       muIniMomX=px; muIniMomY=py; muIniMomZ=pz;
       muIniPolX=xpolaris; muIniPolY=ypolaris; muIniPolZ=zpolaris; 
@@ -107,24 +107,24 @@ class musrRootOutput  {
     }
 
     void SetPolInTarget(G4ThreeVector pol) {muTargetPolX=pol.x(); muTargetPolY=pol.y(); muTargetPolZ=pol.z();}
-    void SetTimeInTarget(G4double time) {muTargetTime = time/microsecond;}
-    void SetMomentumInTarget(G4ThreeVector mom) {muTargetMomX=(mom.x())/MeV; muTargetMomY=(mom.y())/MeV; muTargetMomZ=(mom.z())/MeV;}
+    void SetTimeInTarget(G4double time) {muTargetTime = time/CLHEP::microsecond;}
+    void SetMomentumInTarget(G4ThreeVector mom) {muTargetMomX=(mom.x())/CLHEP::MeV; muTargetMomY=(mom.y())/CLHEP::MeV; muTargetMomZ=(mom.z())/CLHEP::MeV;}
     void SetPolInM0(G4ThreeVector pol) {muM0PolX=pol.x(); muM0PolY=pol.y(); muM0PolZ=pol.z();}
-    void SetTimeInM0(G4double time) {muM0Time = time/microsecond;}
+    void SetTimeInM0(G4double time) {muM0Time = time/CLHEP::microsecond;}
     void SetPolInM1(G4ThreeVector pol) {muM1PolX=pol.x(); muM1PolY=pol.y(); muM1PolZ=pol.z();}
-    void SetTimeInM1(G4double time) {muM1Time = time/microsecond;}
+    void SetTimeInM1(G4double time) {muM1Time = time/CLHEP::microsecond;}
     void SetPolInM2(G4ThreeVector pol) {muM2PolX=pol.x(); muM2PolY=pol.y(); muM2PolZ=pol.z();}
-    void SetTimeInM2(G4double time) {muM2Time = time/microsecond;}
+    void SetTimeInM2(G4double time) {muM2Time = time/CLHEP::microsecond;}
     void SetInitialPositronMomentum(G4ThreeVector mom) {posIniMomx=mom.x();  posIniMomy=mom.y(); posIniMomz=mom.z();}
     void SetNOptPhot(G4int value) {nOptPhot=value;}
     void SetPhotDetTime(G4double time);
-    void SetDecayTime(G4double time) {muDecayTime=time/microsecond;}
+    void SetDecayTime(G4double time) {muDecayTime=time/CLHEP::microsecond;}
     void SetNrFieldNomVal(G4int n) {nFieldNomVal = n;}
     void SetFieldNomVal(G4int i, G4double value);
     G4int GetNrOfVolumes() {return det_nMax;}
     void SetBFieldIntegral(G4double BxInt,G4double ByInt,G4double BzInt,G4double BzInt1,G4double BzInt2,G4double BzInt3) {
-      BxIntegral=BxInt/m/tesla; ByIntegral=ByInt/m/tesla; BzIntegral=BzInt/m/tesla;
-      BzIntegral1=BzInt1/m/tesla;BzIntegral2=BzInt2/mm;BzIntegral3=BzInt3/mm;
+      BxIntegral=BxInt/CLHEP::m/CLHEP::tesla; ByIntegral=ByInt/CLHEP::m/CLHEP::tesla; BzIntegral=BzInt/CLHEP::m/CLHEP::tesla;
+      BzIntegral1=BzInt1/CLHEP::m/CLHEP::tesla;BzIntegral2=BzInt2/CLHEP::mm;BzIntegral3=BzInt3/CLHEP::mm;
     }
     void StoreGeantParameter(Int_t i, Double_t value) {
       if (i<maxNGeantParameters) { GeantParametersD[i]=value; }
diff --git a/include/musrRunAction.hh b/include/musrRunAction.hh
index e8bdf91..3a90b67 100644
--- a/include/musrRunAction.hh
+++ b/include/musrRunAction.hh
@@ -26,7 +26,7 @@
 
 #include "G4UserRunAction.hh"
 #include "globals.hh"
-#include "musrRootOutput.hh"
+//#include "musrRootOutput.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc
index 4c67088..828719b 100644
--- a/src/musrDetectorConstruction.cc
+++ b/src/musrDetectorConstruction.cc
@@ -158,11 +158,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	double pp1, pp2, pp3, pp4=0;
 	sscanf(&line[0],"%*s %*s %s %lf %lf %lf %lf",matrixName,&pp1,&pp2,&pp3,&pp4);
 	if (pp4==0) {
-	  G4RotationMatrix* pRotMatrix = new G4RotationMatrix(pp1*deg,pp2*deg,pp3*deg);  // pp1=phi, pp2=theta, pp3=psi
+	  G4RotationMatrix* pRotMatrix = new G4RotationMatrix(pp1*CLHEP::deg,pp2*CLHEP::deg,pp3*CLHEP::deg);  // pp1=phi, pp2=theta, pp3=psi
 	  pointerToRotationMatrix[matrixName]=pRotMatrix;
 	}
 	else {
-	  G4RotationMatrix* pRotMatrix = new G4RotationMatrix(G4ThreeVector(pp1,pp2,pp3),pp4*deg);
+	  G4RotationMatrix* pRotMatrix = new G4RotationMatrix(G4ThreeVector(pp1,pp2,pp3),pp4*CLHEP::deg);
 	  pointerToRotationMatrix[matrixName]=pRotMatrix;
 	}
       }
@@ -223,14 +223,14 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  solid = new G4Tubs(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  solid = new G4Tubs(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	}
 	else if (strcmp(tmpString2,"cons")==0){
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  solid = new G4Cons(solidName,x1*mm,x2*mm,x3*mm,x4*mm,x5*mm,x6*deg,x7*deg);
+	  solid = new G4Cons(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::mm,x5*CLHEP::mm,x6*CLHEP::deg,x7*CLHEP::deg);
 	}
 	else if (strcmp(tmpString2,"box")==0){
 	  //	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %s %lf %lf %lf %s %s",
@@ -239,21 +239,21 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
 	  //	  G4cout<<"DEBUG KAMIL: xxx1="<<xxx1<<",  xxx2="<<xxx2<<",  xxx3="<<xxx3<<G4endl;
-	  solid = new G4Box(solidName,x1*mm,x2*mm,x3*mm);
+	  solid = new G4Box(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
 	}
         else if ((strcmp(tmpString2,"trd")==0)||(strcmp(tmpString2,"trd90y")==0)) {
           sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
                  name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
           sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
           solidName+=name;
-          solid = new G4Trd(solidName,x1*mm,x2*mm,x3*mm,x4*mm,x5*mm);
+          solid = new G4Trd(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::mm,x5*CLHEP::mm);
         }
 	else if (strcmp(tmpString2,"sphere")==0){
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  solid = new G4Sphere(solidName,x1*mm,x2*mm,x3*deg,x4*deg,x5*deg,x6*deg);
+	  solid = new G4Sphere(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::deg,x4*CLHEP::deg,x5*CLHEP::deg,x6*CLHEP::deg);
 	}
 	else if ((strcmp(tmpString2,"polyconeA")==0)||(strcmp(tmpString2,"polyconeB")==0)){
 	  G4bool polyconeA=true; if (strcmp(tmpString2,"polyconeB")==0) polyconeA=false;
@@ -273,7 +273,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  }
 	  solidName+=name;
 	  G4cout<<name<<G4endl;
-	  //	  solid = new G4Sphere(solidName,x1*mm,x2*mm,x3*deg,x4*deg,x5*deg,x6*deg);
+	  //	  solid = new G4Sphere(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::deg,x4*CLHEP::deg,x5*CLHEP::deg,x6*CLHEP::deg);
 	  //	  double* zPLANE = pointerToArray[zPlane];
 	  iterArray = pointerToArray.find(zPlane);
 	  if (iterArray==pointerToArray.end()) {  // array does not exist
@@ -317,10 +317,10 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  }
 	  //
 	  if (polyconeA) {
-	    solid = new G4Polycone(solidName,x1*deg,x2*deg,numZPlanes,zPLANE,rINNER,rOUTER);
+	    solid = new G4Polycone(solidName,x1*CLHEP::deg,x2*CLHEP::deg,numZPlanes,zPLANE,rINNER,rOUTER);
 	  }
 	  else {  //  polyconeB
-	   solid = new G4Polycone(solidName,x1*deg,x2*deg,numZPlanes,zPLANE,rINNER);
+	   solid = new G4Polycone(solidName,x1*CLHEP::deg,x2*CLHEP::deg,numZPlanes,zPLANE,rINNER);
 	  } 
 	}
 	else if (strcmp(tmpString2,"para")==0){  // NOT YET TESTED
@@ -328,23 +328,23 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  solid = new G4Para(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg,x6*deg);
+	  solid = new G4Para(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg,x6*CLHEP::deg);
 	}
 	else if (strcmp(tmpString2,"torus")==0){  // added JSL. Torus piece, Rmin (bore), Rmax (outer), Rtorus (swept), pSPhi (start angle), pDPhi (swept angle)
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  solid = new G4Torus(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  solid = new G4Torus(solidName,x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	}
         else if (strcmp(tmpString2,"cylpart")==0){  // Volume introduced by Pavel Bakule on 12 May 2009
           sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %s %lf %lf %lf %s %s",
                  name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix);
           sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
           solidName+=name;
-          G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-          G4Box* solidSubtractedBox = new G4Box("solidSubtractedBox",x2*mm,(x2-x4)*mm,x3*mm+roundingErr);
-          G4Tubs* solidCylinder = new G4Tubs("solidCylinder",0.*mm,x2*mm,x3*mm,0.*deg,180.*deg);
+          G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+          G4Box* solidSubtractedBox = new G4Box("solidSubtractedBox",x2*CLHEP::mm,(x2-x4)*CLHEP::mm,x3*CLHEP::mm+roundingErr);
+          G4Tubs* solidCylinder = new G4Tubs("solidCylinder",0.*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,0.*CLHEP::deg,180.*CLHEP::deg);
           solid = new G4SubtractionSolid(solidName, solidCylinder, solidSubtractedBox);
         }
 	else if (strcmp(tmpString2,"uprofile")==0){
@@ -356,11 +356,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
 	  
-	  G4double roundingErr=0.01*mm;
-	  G4Box* box1 = new G4Box("Box1",x1*mm,x2*mm,x3*mm);
-	  G4Box* box2 = new G4Box("Box2",(x1-x4)*mm,(x2-x4/2.+roundingErr)*mm,x3*mm);
+	  G4double roundingErr=0.01*CLHEP::mm;
+	  G4Box* box1 = new G4Box("Box1",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4Box* box2 = new G4Box("Box2",(x1-x4)*CLHEP::mm,(x2-x4/2.+roundingErr)*CLHEP::mm,x3*CLHEP::mm);
 	  G4RotationMatrix rot(0,0,0);
-	  G4ThreeVector zTrans(0,x4/2.*mm,0);
+	  G4ThreeVector zTrans(0,x4/2.*CLHEP::mm,0);
 	  G4Transform3D transform(rot,zTrans);
 	  solid =  new G4SubtractionSolid(solidName, box1, box2, transform);
 	}
@@ -374,13 +374,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
 
-	  G4Box* box1 = new G4Box("Box1",x1*mm,x2*mm,x3*mm);
-	  G4Box* box2 = new G4Box("Box2",(x1+1.0)*mm,x4*sqrt(2.)*mm,x4*sqrt(2.)*mm);
-	  G4RotationMatrix rot(0,45*deg,0);
-	  G4ThreeVector zTransA(0,x2*mm,x3*mm);
-	  G4ThreeVector zTransB(0,-x2*mm,x3*mm);
-	  G4ThreeVector zTransC(0,x2*mm,-x3*mm);
-	  G4ThreeVector zTransD(0,-x2*mm,-x3*mm);
+	  G4Box* box1 = new G4Box("Box1",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4Box* box2 = new G4Box("Box2",(x1+1.0)*CLHEP::mm,x4*sqrt(2.)*CLHEP::mm,x4*sqrt(2.)*CLHEP::mm);
+	  G4RotationMatrix rot(0,45*CLHEP::deg,0);
+	  G4ThreeVector zTransA(0,x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4ThreeVector zTransB(0,-x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4ThreeVector zTransC(0,x2*CLHEP::mm,-x3*CLHEP::mm);
+	  G4ThreeVector zTransD(0,-x2*CLHEP::mm,-x3*CLHEP::mm);
 	  G4Transform3D transformA(rot,zTransA);
 	  G4Transform3D transformB(rot,zTransB);
 	  G4Transform3D transformC(rot,zTransC);
@@ -398,15 +398,15 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
 	  
-	  G4Box* box1 = new G4Box("Box1",x1*mm,4.7*mm,130*mm);
-	  G4Box* box2 = new G4Box("Box2",(x1-0.4)*mm,3*mm,130*mm);
+	  G4Box* box1 = new G4Box("Box1",x1*CLHEP::mm,4.7*CLHEP::mm,130*CLHEP::mm);
+	  G4Box* box2 = new G4Box("Box2",(x1-0.4)*CLHEP::mm,3*CLHEP::mm,130*CLHEP::mm);
 	  G4RotationMatrix rot(0,0,0);
-	  G4ThreeVector zTrans(0,1.3*mm,0);
+	  G4ThreeVector zTrans(0,1.3*CLHEP::mm,0);
 	  G4Transform3D transform(rot,zTrans);
 	  G4SubtractionSolid* solid_1 =  new G4SubtractionSolid("Drzak", box1, box2, transform);
 	  if (x2!=0) {
-	    G4Box* box3 = new G4Box("Box3",12.5*mm,4.5*mm,4*mm);
-	    G4ThreeVector zTrans2(0,(-4.7-x2)*mm,0);
+	    G4Box* box3 = new G4Box("Box3",12.5*CLHEP::mm,4.5*CLHEP::mm,4*CLHEP::mm);
+	    G4ThreeVector zTrans2(0,(-4.7-x2)*CLHEP::mm,0);
 	    G4Transform3D transform2(rot,zTrans2);
 	    solid = new G4UnionSolid("solidName", solid_1, box3, transform2);
 	  }
@@ -434,11 +434,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",0,x1*mm,x6/2*mm+roundingErr,x4*deg,x5*deg);
-	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",0,x2*mm,x3*mm,x4*deg,x5*deg);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",0,x1*CLHEP::mm,x6/2*CLHEP::mm+roundingErr,x4*CLHEP::deg,x5*CLHEP::deg);
+	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",0,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	  G4RotationMatrix rot(0,0,0);
-	  G4ThreeVector zTrans(0,0,(x6/2.-x3)*mm);
+	  G4ThreeVector zTrans(0,0,(x6/2.-x3)*CLHEP::mm);
 	  G4Transform3D transform(rot,zTrans);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetTube, transform);
 	}
@@ -461,11 +461,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",0,x6*mm,x7/2*mm+roundingErr,x4*deg,x5*deg);
-	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",0,x6*CLHEP::mm,x7/2*CLHEP::mm+roundingErr,x4*CLHEP::deg,x5*CLHEP::deg);
+	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	  G4RotationMatrix rot(0,0,0);
-	  G4ThreeVector zTrans(0,0,(x7/2.-x3)*mm);
+	  G4ThreeVector zTrans(0,0,(x7/2.-x3)*CLHEP::mm);
 	  G4Transform3D transform(rot,zTrans);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetTube, transform);
 	}
@@ -492,11 +492,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
 		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",x6*mm-roundingErr,x7*mm+roundingErr,x8/2*mm+roundingErr,x4*deg,x5*deg);
-	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Tubs* solidInnerDetTube = new G4Tubs("SolidInnerDetTube",x6*CLHEP::mm-roundingErr,x7*CLHEP::mm+roundingErr,x8/2*CLHEP::mm+roundingErr,x4*CLHEP::deg,x5*CLHEP::deg);
+	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	  G4RotationMatrix rot(0,0,0);
-	  G4ThreeVector zTrans(0,0,(x8/2.-x3)*mm);
+	  G4ThreeVector zTrans(0,0,(x8/2.-x3)*CLHEP::mm);
 	  G4Transform3D transform(rot,zTrans);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetTube, transform);
 	}
@@ -506,9 +506,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Box* solidInnerDetBox = new G4Box("SolidInnerDetBox",x1*mm,x1*mm,x3*mm+roundingErr);
-	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Box* solidInnerDetBox = new G4Box("SolidInnerDetBox",x1*CLHEP::mm,x1*CLHEP::mm,x3*CLHEP::mm+roundingErr);
+	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",0.*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetBox);
 	}
 	else if (strcmp(tmpString2,"tubsbox2")==0){
@@ -519,9 +519,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
 		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Box* solidInnerDetBox2 = new G4Box("SolidInnerDetBox2",x1*mm,x2*mm,x3*mm+roundingErr);
-	  G4Tubs* solidOuterDetTube2 = new G4Tubs("SolidOuterDetTube2",x4*mm,x5*mm,x6*mm,x7*deg,x8*deg);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Box* solidInnerDetBox2 = new G4Box("SolidInnerDetBox2",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm+roundingErr);
+	  G4Tubs* solidOuterDetTube2 = new G4Tubs("SolidOuterDetTube2",x4*CLHEP::mm,x5*CLHEP::mm,x6*CLHEP::mm,x7*CLHEP::deg,x8*CLHEP::deg);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube2, solidInnerDetBox2);
 	}
 	else if (strcmp(tmpString2,"boxbox")==0){
@@ -532,9 +532,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
 		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Box* solidInnerDetBox = new G4Box("SolidInnerDetBox",x1*mm+roundingErr,x2*mm+roundingErr,x3*mm+roundingErr);
-	  G4Box* solidOuterDetBox = new G4Box("SolidOuterDetBox",x4*mm,x5*mm,x6*mm);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Box* solidInnerDetBox = new G4Box("SolidInnerDetBox",x1*CLHEP::mm+roundingErr,x2*CLHEP::mm+roundingErr,x3*CLHEP::mm+roundingErr);
+	  G4Box* solidOuterDetBox = new G4Box("SolidOuterDetBox",x4*CLHEP::mm,x5*CLHEP::mm,x6*CLHEP::mm);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetBox, solidInnerDetBox);
 	}
 	else if (strcmp(tmpString2,"GPSforward")==0){
@@ -545,9 +545,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Box* GPS_SolidDetBox = new G4Box("GPS_SolidDetBox",x1*mm,x2*mm,x3*mm);
-	  G4Cons* GPS_SolidDetCons = new G4Cons("GPS_SolidDetCons",0.,x4*mm,0.,x5*mm,x3*mm+roundingErr,0.,360.);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Box* GPS_SolidDetBox = new G4Box("GPS_SolidDetBox",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4Cons* GPS_SolidDetCons = new G4Cons("GPS_SolidDetCons",0.,x4*CLHEP::mm,0.,x5*CLHEP::mm,x3*CLHEP::mm+roundingErr,0.,360.);
 	  solid =  new G4SubtractionSolid(solidName, GPS_SolidDetBox, GPS_SolidDetCons);
 	}
 	else if (strcmp(tmpString2,"GPSbackward")==0){
@@ -558,9 +558,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Box* GPS_SolidDetBox2 = new G4Box("GPS_SolidDetBox2",x1*mm,x2*mm,x3*mm);
-	  G4Trd* GPS_SolidDetTrd2 = new G4Trd("GPS_SolidDetTrd2",x4*mm,x5*mm,x4*mm,x5*mm,x3*mm+roundingErr);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Box* GPS_SolidDetBox2 = new G4Box("GPS_SolidDetBox2",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4Trd* GPS_SolidDetTrd2 = new G4Trd("GPS_SolidDetTrd2",x4*CLHEP::mm,x5*CLHEP::mm,x4*CLHEP::mm,x5*CLHEP::mm,x3*CLHEP::mm+roundingErr);
 	  solid =  new G4SubtractionSolid(solidName, GPS_SolidDetBox2, GPS_SolidDetTrd2);
 	}
 	else if (strcmp(tmpString2,"GPSbackwardVeto")==0){
@@ -572,9 +572,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
 		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Trd* GPS_OuterTrd = new G4Trd("GPS_OuterTrd",x1*mm,x2*mm,x3*mm,x4*mm,x5*mm);
-	  G4Trd* GPS_InnerTrd = new G4Trd("GPS_InnerTrd",x6*mm,x7*mm,x8*mm,x9*mm,x5*mm+roundingErr);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Trd* GPS_OuterTrd = new G4Trd("GPS_OuterTrd",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::mm,x5*CLHEP::mm);
+	  G4Trd* GPS_InnerTrd = new G4Trd("GPS_InnerTrd",x6*CLHEP::mm,x7*CLHEP::mm,x8*CLHEP::mm,x9*CLHEP::mm,x5*CLHEP::mm+roundingErr);
 	  solid =  new G4SubtractionSolid(solidName,GPS_OuterTrd,GPS_InnerTrd);
 	}
 	else if (strcmp(tmpString2,"tubeWithWindows")==0){
@@ -591,13 +591,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
 		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Tubs* solidShield = new G4Tubs("SolidShield",x1*mm,x2*mm,x3*mm,0.,360.);
-	  G4Tubs* solidWindF = new G4Tubs("SolidWindF",0.,x4*mm,x2*mm,0.,360.);
-	  G4Tubs* solidWindB = new G4Tubs("SolidWindB",0.,x6*mm,x2*mm,0.,360.);
-	  G4RotationMatrix* rot = new G4RotationMatrix(0,90*deg,0);
-	  G4ThreeVector zTransF(0.,x2*mm+roundingErr,(x5-x3)*mm);
-	  G4ThreeVector zTransB(0.,-x2*mm+roundingErr,(x7-x3)*mm);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Tubs* solidShield = new G4Tubs("SolidShield",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,0.,360.);
+	  G4Tubs* solidWindF = new G4Tubs("SolidWindF",0.,x4*CLHEP::mm,x2*CLHEP::mm,0.,360.);
+	  G4Tubs* solidWindB = new G4Tubs("SolidWindB",0.,x6*CLHEP::mm,x2*CLHEP::mm,0.,360.);
+	  G4RotationMatrix* rot = new G4RotationMatrix(0,90*CLHEP::deg,0);
+	  G4ThreeVector zTransF(0.,x2*CLHEP::mm+roundingErr,(x5-x3)*CLHEP::mm);
+	  G4ThreeVector zTransB(0.,-x2*CLHEP::mm+roundingErr,(x7-x3)*CLHEP::mm);
 	  G4SubtractionSolid* solidF1 =  new G4SubtractionSolid("solidF1", solidShield, solidWindF, rot, zTransF);
 	  solid =  new G4SubtractionSolid(solidName, solidF1, solidWindB, rot, zTransB);
 	}
@@ -610,11 +610,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",
 		 &posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4Box* solidDetBox = new G4Box("SolidDetBox",x1*mm,x2*mm,x3*mm);
-	  G4Box* solidHole   = new G4Box("SolidDetBox",x11*mm,x2*mm+0.1,x12*mm);
-	  G4Tubs* solidDetTube = new G4Tubs("SolidDetTube",0.,x4*mm,x5*mm,x6*deg,x7*deg);
+	  G4Box* solidDetBox = new G4Box("SolidDetBox",x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
+	  G4Box* solidHole   = new G4Box("SolidDetBox",x11*CLHEP::mm,x2*CLHEP::mm+0.1,x12*CLHEP::mm);
+	  G4Tubs* solidDetTube = new G4Tubs("SolidDetTube",0.,x4*CLHEP::mm,x5*CLHEP::mm,x6*CLHEP::deg,x7*CLHEP::deg);
 	  G4RotationMatrix* yRot = new G4RotationMatrix();
-	  G4ThreeVector zTrans(-x8*mm,-x9*mm,-x10*mm);
+	  G4ThreeVector zTrans(-x8*CLHEP::mm,-x9*CLHEP::mm,-x10*CLHEP::mm);
 	  G4SubtractionSolid* solidPart1 =  new G4SubtractionSolid("solidPart1", solidDetBox, solidDetTube, yRot, zTrans);
 	  solid =  new G4SubtractionSolid(solidName,solidPart1,solidHole);
 	}
@@ -625,15 +625,15 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,orientation,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Box* solidDetBox = new G4Box("SolidDetBox",x2*mm,x2*mm,x3*mm+roundingErr);
-	  G4Tubs* solidDetTube = new G4Tubs("SolidDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Box* solidDetBox = new G4Box("SolidDetBox",x2*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm+roundingErr);
+	  G4Tubs* solidDetTube = new G4Tubs("SolidDetTube",0.*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm,x4*CLHEP::deg,x5*CLHEP::deg);
 	  G4RotationMatrix* yRot = new G4RotationMatrix;
 	  G4ThreeVector zTrans;
-	  if (strcmp(orientation,"D")==0) zTrans=G4ThreeVector((x1-x2)*mm,(-x1-x2)*mm,0);
-	  else if (strcmp(orientation,"U")==0) zTrans=G4ThreeVector((x2-x1)*mm,(x1+x2)*mm,0);
-	  else if (strcmp(orientation,"R")==0) zTrans=G4ThreeVector((-x1-x2)*mm,(x2-x1)*mm,0);
-	  else if (strcmp(orientation,"L")==0) zTrans=G4ThreeVector((x1+x2)*mm,(x1-x2)*mm,0);
+	  if (strcmp(orientation,"D")==0) zTrans=G4ThreeVector((x1-x2)*CLHEP::mm,(-x1-x2)*CLHEP::mm,0);
+	  else if (strcmp(orientation,"U")==0) zTrans=G4ThreeVector((x2-x1)*CLHEP::mm,(x1+x2)*CLHEP::mm,0);
+	  else if (strcmp(orientation,"R")==0) zTrans=G4ThreeVector((-x1-x2)*CLHEP::mm,(x2-x1)*CLHEP::mm,0);
+	  else if (strcmp(orientation,"L")==0) zTrans=G4ThreeVector((x1+x2)*CLHEP::mm,(x1-x2)*CLHEP::mm,0);
 	  else {
 	    G4cout<<"Unknown orientation of the tubsboxsegm volume!!"
 		  <<"  orientation="<<orientation<<G4endl;
@@ -643,27 +643,27 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	}
 	else if (strcmp(tmpString2,"GPDsampleHolderA")==0){
 	  // First part of the GPD sample holder, where posx, posy, posz = centre of the whole (long) tube 
-	  //                                                      (=111.25mm below the centre of the holes)
+	  //                                                      (=111.25CLHEP::mm below the centre of the holes)
 	  sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID);
 	  solidName+=name;
-	  G4double upperPartHalfHeight = 65*mm;
-	  G4double lowerPartHalfHeight = (287.5-65)/2*mm;
+	  G4double upperPartHalfHeight = 65*CLHEP::mm;
+	  G4double lowerPartHalfHeight = (287.5-65)/2*CLHEP::mm;
 	  G4double halfHeight = upperPartHalfHeight+lowerPartHalfHeight;
-	  G4double innerR = 37.*mm; G4double outerR = 37.5*mm;
+	  G4double innerR = 37.*CLHEP::mm; G4double outerR = 37.5*CLHEP::mm;
 	  
-	  //	  G4Box* solidDetBox = new G4Box("SolidDetBox",x2*mm,x2*mm,x3*mm+roundingErr);
+	  //	  G4Box* solidDetBox = new G4Box("SolidDetBox",x2*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm+roundingErr);
 	  G4Tubs* solidGPDTubeA1 = new G4Tubs("SolidGPDTubeA1",innerR,outerR,halfHeight,0,360);
-	  G4Box*  solidGPDBoxA2  = new G4Box ("SolidGPDBoxA2",10*mm,38*mm, 70*mm);
-	  G4Box*  solidGPDBoxA5  = new G4Box ("solidGPDBoxA5",5*mm,40*mm,30*mm); 
+	  G4Box*  solidGPDBoxA2  = new G4Box ("SolidGPDBoxA2",10*CLHEP::mm,38*CLHEP::mm, 70*CLHEP::mm);
+	  G4Box*  solidGPDBoxA5  = new G4Box ("solidGPDBoxA5",5*CLHEP::mm,40*CLHEP::mm,30*CLHEP::mm); 
 
 	  G4RotationMatrix* yRotA12 = new G4RotationMatrix();
-	  G4ThreeVector zTransA12( 30*mm,0, (111.25+20)*mm);
-	  G4ThreeVector zTransA13(-30*mm,0, (111.25+20)*mm);
+	  G4ThreeVector zTransA12( 30*CLHEP::mm,0, (111.25+20)*CLHEP::mm);
+	  G4ThreeVector zTransA13(-30*CLHEP::mm,0, (111.25+20)*CLHEP::mm);
 	  G4SubtractionSolid* solidA12 =  new G4SubtractionSolid("solidA12", solidGPDTubeA1, solidGPDBoxA2, yRotA12, zTransA12);
 	  G4SubtractionSolid* solidA123 =  new G4SubtractionSolid("solidA123", solidA12 , solidGPDBoxA2, yRotA12, zTransA13);
-	  G4ThreeVector zTransA5(0,0,111.25*mm);
+	  G4ThreeVector zTransA5(0,0,111.25*CLHEP::mm);
 	  solid =  new G4SubtractionSolid(solidName, solidA123, solidGPDBoxA5, yRotA12, zTransA5);
 	}
 	else if (strcmp(tmpString2,"GPDmHolder")==0){
@@ -672,12 +672,12 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID);
 	  solidName+=name;
-	  G4Box* solidGPDcutOut = new G4Box ("solidGPDcutOut",x1*mm,x1*mm,(x3+0.01)*mm);
-	  G4Box* solidGPDmHolder = new G4Box ("solidGPDmHolder",sqrt((double)2)*x1*mm,x2*mm,x3*mm);
-	  //	  G4RotationMatrix* rot = new G4RotationMatrix(45*deg,0,0);
-	  G4RotationMatrix* rot = new G4RotationMatrix(G4ThreeVector(0,0,1),45*deg);
+	  G4Box* solidGPDcutOut = new G4Box ("solidGPDcutOut",x1*CLHEP::mm,x1*CLHEP::mm,(x3+0.01)*CLHEP::mm);
+	  G4Box* solidGPDmHolder = new G4Box ("solidGPDmHolder",sqrt((double)2)*x1*CLHEP::mm,x2*CLHEP::mm,x3*CLHEP::mm);
+	  //	  G4RotationMatrix* rot = new G4RotationMatrix(45*CLHEP::deg,0,0);
+	  G4RotationMatrix* rot = new G4RotationMatrix(G4ThreeVector(0,0,1),45*CLHEP::deg);
 	  //	  G4RotationMatrix* rot = new G4RotationMatrix();
-	  G4ThreeVector  trans(0,x2*mm,0);
+	  G4ThreeVector  trans(0,x2*CLHEP::mm,0);
 	  solid =  new G4SubtractionSolid(solidName,solidGPDmHolder,solidGPDcutOut,rot,trans);
 	  //	  solid =  new G4SubtractionSolid(solidName,solidGPDcutOut,solidGPDmHolder,rot,trans);
 	}
@@ -689,23 +689,23 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
 	  solidName+=name;
-	  //G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  G4Tubs* solidMainBore       = new G4Tubs("solidMainBore"      ,0.*mm,x1*mm,x4*mm,0*deg,360*deg);
-	  G4Tubs* solidHorizCrossBore = new G4Tubs("solidHorizCrossBore",0.*mm,x2*mm,x5*mm,0*deg,360*deg);
-	  G4Tubs* solidBottomBore     = new G4Tubs("solidBottomBore"    ,0.*mm,x2*mm,x5/2*mm,0*deg,360*deg);
-	  G4Tubs* solidTopBore        = new G4Tubs("solidTopBore"       ,0.*mm,x3*mm,x5/2*mm,0*deg,360*deg);
-	  G4RotationMatrix* rotx = new G4RotationMatrix(G4ThreeVector(0.,1.,0.),90.*deg);
-	  G4RotationMatrix* roty = new G4RotationMatrix(G4ThreeVector(1.,0.,0.),90.*deg);
+	  //G4double roundingErr=0.01*CLHEP::mm;  // to avoid some displaying problems of the subtracted volumes
+	  G4Tubs* solidMainBore       = new G4Tubs("solidMainBore"      ,0.*CLHEP::mm,x1*CLHEP::mm,x4*CLHEP::mm,0*CLHEP::deg,360*CLHEP::deg);
+	  G4Tubs* solidHorizCrossBore = new G4Tubs("solidHorizCrossBore",0.*CLHEP::mm,x2*CLHEP::mm,x5*CLHEP::mm,0*CLHEP::deg,360*CLHEP::deg);
+	  G4Tubs* solidBottomBore     = new G4Tubs("solidBottomBore"    ,0.*CLHEP::mm,x2*CLHEP::mm,x5/2*CLHEP::mm,0*CLHEP::deg,360*CLHEP::deg);
+	  G4Tubs* solidTopBore        = new G4Tubs("solidTopBore"       ,0.*CLHEP::mm,x3*CLHEP::mm,x5/2*CLHEP::mm,0*CLHEP::deg,360*CLHEP::deg);
+	  G4RotationMatrix* rotx = new G4RotationMatrix(G4ThreeVector(0.,1.,0.),90.*CLHEP::deg);
+	  G4RotationMatrix* roty = new G4RotationMatrix(G4ThreeVector(1.,0.,0.),90.*CLHEP::deg);
 	  G4ThreeVector nowhere( 0,0,0);
-	  G4ThreeVector UpOffset( 0., x5/2.*mm, 0.);
-	  G4ThreeVector DownOffset( 0., -x5/2.*mm, 0.);
+	  G4ThreeVector UpOffset( 0., x5/2.*CLHEP::mm, 0.);
+	  G4ThreeVector DownOffset( 0., -x5/2.*CLHEP::mm, 0.);
 	  G4UnionSolid* partXZonly = new G4UnionSolid("partXZonly", solidMainBore, solidHorizCrossBore, rotx, nowhere);
 	  G4UnionSolid* partXZBonly = new G4UnionSolid("partXZBonly", partXZonly, solidBottomBore, roty, DownOffset);
 	  solid =  new G4UnionSolid(solidName, partXZBonly, solidTopBore, roty, UpOffset);
 	}
 	else ReportGeometryProblem(line);
 
-	G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm);
+	G4ThreeVector position = G4ThreeVector (posx*CLHEP::mm,posy*CLHEP::mm,posz*CLHEP::mm);
 	//	G4cout << "New volume: "<<tmpString2<<" "<<name<<", solid geometry parameters="<<x1<<" "<<x2<<" "<<x3<<" "<<x4<<" "<<x5
 	//	       << ", position="<<position<<", mother="<<mothersName<<G4endl;
 
@@ -737,10 +737,10 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	if (strcmp(sensitiveDet,"dead")) {G4cout<<"     sensitive:   ";}
 	else {G4cout<<" non-sensitive:   ";}
 		 
-        // For trapezoid "trd90y" rotate the trapezoid by 90 degrees in y direction
+        // For trapezoid "trd90y" rotate the trapezoid by 90 CLHEP::degrees in y direction
 	// (in addition to the requested rotation)
 	if (strcmp(tmpString2,"trd90y")==0) {
-	  pRot->rotateY(90.0*deg);
+	  pRot->rotateY(90.0*CLHEP::deg);
         }
 
 	G4FieldManager* pFieldMan = pointerToField[actualFieldName];
@@ -990,7 +990,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  }
 	  else if (strcmp(varName,"signalSeparationTime")==0) {
 	    sscanf(&line[0],"%*s %*s %*s %lf",&fVarValue);
-	    myMusrScintSD -> Set_OPSA_SignalSeparationTime(fVarValue*nanosecond);
+	    myMusrScintSD -> Set_OPSA_SignalSeparationTime(fVarValue*CLHEP::nanosecond);
 	  }
 	  else if (strcmp(varName,"photonFractions")==0) {
 	    double a, b, c, d;
@@ -1006,7 +1006,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	    int nBins;
 	    double min, max;
 	    sscanf(&line[0],"%*s %*s %*s %d %lf %lf",&nBins,&min,&max);
-	    myMusrScintSD -> SetOPSAhistoBinning(nBins,min*nanosecond,max*nanosecond);
+	    myMusrScintSD -> SetOPSAhistoBinning(nBins,min*CLHEP::nanosecond,max*CLHEP::nanosecond);
 	  }
 	  else if (strcmp(varName,"pulseShapeArray")==0) {
 	    char fileName[500];
@@ -1027,7 +1027,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  else if (strcmp(varName,"SetAPDcellsTimeVariationSigma")==0) {
 	    double sigma;
 	    sscanf(&line[0],"%*s %*s %*s %lf",&sigma);
-	    myMusrScintSD ->SetAPDcellsTimeVariationSigma(sigma*nanosecond);
+	    myMusrScintSD ->SetAPDcellsTimeVariationSigma(sigma*CLHEP::nanosecond);
 	  }
 	  else if (strcmp(varName,"SetAPDcrossTalk")==0) {
 	    double crossTalkProb;
@@ -1074,11 +1074,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  	  
 	  // Construct the field
 	  musrTabulatedElementField*  myElementTableField = 
-	    new musrTabulatedElementField(fieldInputFileName, fieldTableType, fieldValue*tesla, logVol, position);
+	    new musrTabulatedElementField(fieldInputFileName, fieldTableType, fieldValue*CLHEP::tesla, logVol, position);
 	  myElementTableField->SetElementFieldName(tmpString2);
 	  if (fieldNrOfSteps>0) {
 	    //cks The following line might require some correction for the electric field
-	    myElementTableField->SetEventNrDependentField(fieldValue*tesla,fieldValueFinal*tesla,fieldNrOfSteps);
+	    myElementTableField->SetEventNrDependentField(fieldValue*CLHEP::tesla,fieldValueFinal*CLHEP::tesla,fieldNrOfSteps);
 	  }
 	  //	    FieldList* fields = F04GlobalField::getObject()->getFields();
 	  //	    if (fields) {
@@ -1104,8 +1104,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	    musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
 	  }
 	  G4double fieldValue_tmp[6] = {
-	  fieldValue[0]*tesla,        fieldValue[1]*tesla,        fieldValue[2]*tesla,
-	  fieldValue[3]*(kilovolt/mm),fieldValue[4]*(kilovolt/mm),fieldValue[5]*(kilovolt/mm)};
+	  fieldValue[0]*CLHEP::tesla,        fieldValue[1]*CLHEP::tesla,        fieldValue[2]*CLHEP::tesla,
+	  fieldValue[3]*(CLHEP::kilovolt/CLHEP::mm),fieldValue[4]*(CLHEP::kilovolt/CLHEP::mm),fieldValue[5]*(CLHEP::kilovolt/CLHEP::mm)};
 	  
 	  musrUniformField* myElementUniformField = new musrUniformField(fieldValue_tmp, half_x, half_y, half_z, logVol, position);
 	  myElementUniformField->SetElementFieldName(tmpString2);
@@ -1123,10 +1123,10 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		    logicalVolumeName);
 	    musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
 	  }
-	  musrQuadrupole* myMusrQuadrupole = new musrQuadrupole(halfLength*mm,fieldRadius*mm,gradientValue*(tesla/m),fringeFactor,logVol,position);
+	  musrQuadrupole* myMusrQuadrupole = new musrQuadrupole(halfLength*CLHEP::mm,fieldRadius*CLHEP::mm,gradientValue*(CLHEP::tesla/CLHEP::m),fringeFactor,logVol,position);
 	  myMusrQuadrupole->SetElementFieldName(tmpString2);
 	  if (gradientNrOfSteps>0) {
-	    myMusrQuadrupole->SetEventNrDependentField(gradientValue*(tesla/m),gradientValueFinal*(tesla/m),gradientNrOfSteps);
+	    myMusrQuadrupole->SetEventNrDependentField(gradientValue*(CLHEP::tesla/CLHEP::m),gradientValueFinal*(CLHEP::tesla/CLHEP::m),gradientNrOfSteps);
 	  }
 	}
 
@@ -1148,11 +1148,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	    char parameterName[100];
 	    double parameterValue;
 	    sscanf(&line[0],"%*s %*s %*s %s %lf",parameterName,&parameterValue);
-	    if (strcmp(parameterName,"SetDeltaIntersection")==0){ fieldMgr->SetDeltaIntersection(parameterValue*mm); }
-	    else if (strcmp(parameterName,"SetDeltaOneStep")==0){ fieldMgr->SetDeltaOneStep(parameterValue*mm); }
+	    if (strcmp(parameterName,"SetDeltaIntersection")==0){ fieldMgr->SetDeltaIntersection(parameterValue*CLHEP::mm); }
+	    else if (strcmp(parameterName,"SetDeltaOneStep")==0){ fieldMgr->SetDeltaOneStep(parameterValue*CLHEP::mm); }
 	    else if (strcmp(parameterName,"SetMinimumEpsilonStep")==0){ fieldMgr->SetMinimumEpsilonStep(parameterValue); }
 	    else if (strcmp(parameterName,"SetMaximumEpsilonStep")==0){ fieldMgr->SetMaximumEpsilonStep(parameterValue); }
-	    else if (strcmp(parameterName,"SetLargestAcceptableStep")==0) { propagMgr->SetLargestAcceptableStep(parameterValue*mm); }
+	    else if (strcmp(parameterName,"SetLargestAcceptableStep")==0) { propagMgr->SetLargestAcceptableStep(parameterValue*CLHEP::mm); }
 	    else if (strcmp(parameterName,"SetMaxLoopCount")==0) {propagMgr->SetMaxLoopCount(int(parameterValue)); }
 	    else {
 	      G4cout<<"musrDetectorConstruction.cc: ERROR: Unknown parameterName \""
@@ -1170,8 +1170,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	    G4cout<<"                          Can not print the accuracy parameters of the magnetic field."<<G4endl;
 	  }
 	  else { 
-	    G4cout<<"GetDeltaIntersection()    = "<<fieldMgr->GetDeltaIntersection()/mm<<" mm"<<G4endl;
-	    G4cout<<"GetDeltaOneStep()         = "<<fieldMgr->GetDeltaOneStep()/mm<<" mm"<<G4endl;
+	    G4cout<<"GetDeltaIntersection()    = "<<fieldMgr->GetDeltaIntersection()/CLHEP::mm<<" CLHEP::mm"<<G4endl;
+	    G4cout<<"GetDeltaOneStep()         = "<<fieldMgr->GetDeltaOneStep()/CLHEP::mm<<" CLHEP::mm"<<G4endl;
 	    G4cout<<"GetMinimumEpsilonStep()   = "<<fieldMgr->GetMinimumEpsilonStep()<<G4endl;
 	    G4cout<<"GetMaximumEpsilonStep()   = "<<fieldMgr->GetMaximumEpsilonStep()<<G4endl;
 	  }
@@ -1181,7 +1181,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	    G4cout<<"                          Can not print the accuracy parameters of the magnetic field."<<G4endl;
 	  }
 	  else {
-	    G4cout<<"GetLargestAcceptableStep()= "<<propagMgr->GetLargestAcceptableStep()/mm<<" mm"<<G4endl;
+	    G4cout<<"GetLargestAcceptableStep()= "<<propagMgr->GetLargestAcceptableStep()/CLHEP::mm<<" CLHEP::mm"<<G4endl;
 	    G4cout<<"GetMaxLoopCount()         = "<<propagMgr->GetMaxLoopCount()<<G4endl;
 	  }	    
 	}
@@ -1225,11 +1225,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %lf %lf", &ustepMax, &utrakMax, &utimeMax, &uekinMin, &urangMin);
 	G4UserLimits* myUserLimits = new G4UserLimits();
 	G4cout<<"musrDetectorConstruction.cc:  G4UserLimits in "<<tmpString2<<": ";
-	if (ustepMax>0) {myUserLimits->SetMaxAllowedStep(ustepMax*mm);    G4cout<<"ustepMax = "<<ustepMax<<" mm, ";}
-	if (utrakMax>0) {myUserLimits->SetUserMaxTrackLength(utrakMax*mm);G4cout<<"utrakMax = "<<utrakMax<<" mm, ";}
- 	if (utimeMax>0) {myUserLimits->SetUserMaxTime(utimeMax*ns);       G4cout<<"utimeMax = "<<utimeMax<<" ns, ";}
- 	if (uekinMin>0) {myUserLimits->SetUserMinEkine(uekinMin*MeV);     G4cout<<"uekinMin = "<<uekinMin<<" MeV, ";}
- 	if (urangMin>0) {myUserLimits->SetUserMinRange(urangMin*mm);      G4cout<<"urangMin = "<<urangMin<<" mm, ";}
+	if (ustepMax>0) {myUserLimits->SetMaxAllowedStep(ustepMax*CLHEP::mm);    G4cout<<"ustepMax = "<<ustepMax<<" CLHEP::mm, ";}
+	if (utrakMax>0) {myUserLimits->SetUserMaxTrackLength(utrakMax*CLHEP::mm);G4cout<<"utrakMax = "<<utrakMax<<" CLHEP::mm, ";}
+ 	if (utimeMax>0) {myUserLimits->SetUserMaxTime(utimeMax*CLHEP::ns);       G4cout<<"utimeMax = "<<utimeMax<<" ns, ";}
+ 	if (uekinMin>0) {myUserLimits->SetUserMinEkine(uekinMin*CLHEP::MeV);     G4cout<<"uekinMin = "<<uekinMin<<" MeV, ";}
+ 	if (urangMin>0) {myUserLimits->SetUserMinRange(urangMin*CLHEP::mm);      G4cout<<"urangMin = "<<urangMin<<" CLHEP::mm, ";}
 	G4cout<<G4endl;
 	pLogVol->SetUserLimits(myUserLimits);
       }
@@ -1321,14 +1321,14 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		 << tmpString3 <<"\" not found!"<<G4endl<<"S T O P    F O R C E D"<<G4endl;
 	  ReportGeometryProblem(line);
 	}
-	G4cout<<"Mass of the detector "<<tmpString2<<" is "<<massVol->GetMass()/kg<<" kg."<<G4endl;
+	G4cout<<"Mass of the detector "<<tmpString2<<" is "<<massVol->GetMass()/CLHEP::kg<<" kg."<<G4endl;
       }
       
 
       else if (strcmp(tmpString1,"signalSeparationTime")==0){
 	double timeSeparation;
 	sscanf(&line[0],"%*s %*s %lf",&timeSeparation);
-        musrParameters::signalSeparationTime = timeSeparation*nanosecond;
+        musrParameters::signalSeparationTime = timeSeparation*CLHEP::nanosecond;
       }
 
       else if (strcmp(tmpString1,"maximumRunTimeAllowed")==0){   // in seconds
@@ -1590,50 +1590,50 @@ void musrDetectorConstruction::DefineMaterials()
   G4Element* Mg = man->FindOrBuildElement("Mg");
 
   // compounds required for MCP Macor
-  G4Material* MgO = new G4Material("MgO", 3.60*g/cm3, ncomponents=2);
+  G4Material* MgO = new G4Material("MgO", 3.60*CLHEP::g/CLHEP::cm3, ncomponents=2);
   MgO->AddElement(Mg, natoms=1);
   MgO->AddElement(O,  natoms=1);
 
-  G4Material* SiO2 = new G4Material("SiO2", 2.533*g/cm3, ncomponents=2); // quartz
+  G4Material* SiO2 = new G4Material("SiO2", 2.533*CLHEP::g/CLHEP::cm3, ncomponents=2); // quartz
   SiO2->AddElement(O,  natoms=2);
   SiO2->AddElement(Si, natoms=1);
 
-  G4Material* Al2O3 = new G4Material("Al2O3", 3.985*g/cm3, ncomponents=2); // saphire
+  G4Material* Al2O3 = new G4Material("Al2O3", 3.985*CLHEP::g/CLHEP::cm3, ncomponents=2); // saphire
   Al2O3->AddElement (Al, natoms=2);
   Al2O3->AddElement (O,  natoms=3);
  
-  G4Material* K2O = new G4Material("K2O", 2.350*g/cm3, ncomponents=2);
+  G4Material* K2O = new G4Material("K2O", 2.350*CLHEP::g/CLHEP::cm3, ncomponents=2);
   K2O->AddElement(O, natoms=1);
   K2O->AddElement(K, natoms=2);
   
-  G4Material* B2O3 = new G4Material("B2O3", 2.550*g/cm3, ncomponents=2);
+  G4Material* B2O3 = new G4Material("B2O3", 2.550*CLHEP::g/CLHEP::cm3, ncomponents=2);
   B2O3->AddElement (B, natoms=2);
   B2O3->AddElement (O, natoms=3);
 
   G4Material* Sci = 
-    new G4Material("Scintillator", density= 1.032*g/cm3, ncomponents=2);
+    new G4Material("Scintillator", density= 1.032*CLHEP::g/CLHEP::cm3, ncomponents=2);
   Sci->AddElement(C, natoms=9);
   Sci->AddElement(H, natoms=10);
 
   G4Material* Myl = 
-    new G4Material("Mylar", density= 1.397*g/cm3, ncomponents=3);
+    new G4Material("Mylar", density= 1.397*CLHEP::g/CLHEP::cm3, ncomponents=3);
   Myl->AddElement(C, natoms=10);
   Myl->AddElement(H, natoms= 8);
   Myl->AddElement(O, natoms= 4);
 
   // Brass
-  G4Material* brass = new G4Material("Brass", density= 8.40*g/cm3, ncomponents=2);
-  brass -> AddElement(Zn, fractionmass = 30*perCent);
-  brass -> AddElement(Cu, fractionmass = 70*perCent);
+  G4Material* brass = new G4Material("Brass", density= 8.40*CLHEP::g/CLHEP::cm3, ncomponents=2);
+  brass -> AddElement(Zn, fractionmass = 30*CLHEP::perCent);
+  brass -> AddElement(Cu, fractionmass = 70*CLHEP::perCent);
 
   // Stainless steel
-  G4Material* steel = new G4Material("Steel", density= 7.93*g/cm3, ncomponents=3);
+  G4Material* steel = new G4Material("Steel", density= 7.93*CLHEP::g/CLHEP::cm3, ncomponents=3);
   steel->AddElement(Ni, fractionmass=0.11);
   steel->AddElement(Cr, fractionmass=0.18);
   steel->AddElement(Fe, fractionmass=0.71);
 
   G4Material* macor= // Macor (used in the MCP detector)
-    new G4Material("Macor", density=2.52*g/cm3, ncomponents=5);
+    new G4Material("Macor", density=2.52*CLHEP::g/CLHEP::cm3, ncomponents=5);
   macor->AddMaterial(SiO2, fractionmass=0.470); // quartz
   macor->AddMaterial(MgO,  fractionmass=0.180);
   macor->AddMaterial(Al2O3,fractionmass=0.170); // saphire
@@ -1641,7 +1641,7 @@ void musrDetectorConstruction::DefineMaterials()
   macor->AddMaterial(B2O3, fractionmass=0.075);
     
   G4Material* mcpglass = // Glass of the Multi Channel Plate
-    new G4Material("MCPglass", density=2.0*g/cm3, ncomponents=9);
+    new G4Material("MCPglass", density=2.0*CLHEP::g/CLHEP::cm3, ncomponents=9);
   mcpglass->AddElement(Pb, fractionmass= 0.480);
   mcpglass->AddElement(O,  fractionmass= 0.258);
   mcpglass->AddElement(Si, fractionmass= 0.182);
@@ -1657,42 +1657,42 @@ void musrDetectorConstruction::DefineMaterials()
 //
 
   G4Material* Air = 
-    new G4Material("Air"  , density= 1.290*mg/cm3, ncomponents=2);
+    new G4Material("Air"  , density= 1.290*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   Air->AddElement(N, fractionmass=0.7);
   Air->AddElement(O, fractionmass=0.3);
 
   G4Material* Air1mbar = 
-    new G4Material("Air1mbar"  , density= 1.290e-3*mg/cm3, ncomponents=2);
+    new G4Material("Air1mbar"  , density= 1.290e-3*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   Air1mbar->AddElement(N, fractionmass=0.7);
   Air1mbar->AddElement(O, fractionmass=0.3);
 
   G4Material* AirE1mbar = 
-    new G4Material("AirE1mbar"  , density= 1.290e-4*mg/cm3, ncomponents=2);
+    new G4Material("AirE1mbar"  , density= 1.290e-4*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   AirE1mbar->AddElement(N, fractionmass=0.7);
   AirE1mbar->AddElement(O, fractionmass=0.3);
 
   G4Material* AirE2mbar = 
-    new G4Material("AirE2mbar"  , density= 1.290e-5*mg/cm3, ncomponents=2);
+    new G4Material("AirE2mbar"  , density= 1.290e-5*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   AirE2mbar->AddElement(N, fractionmass=0.7);
   AirE2mbar->AddElement(O, fractionmass=0.3);
 
   G4Material* AirE3mbar = 
-    new G4Material("AirE3mbar"  , density= 1.290e-6*mg/cm3, ncomponents=2);
+    new G4Material("AirE3mbar"  , density= 1.290e-6*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   AirE3mbar->AddElement(N, fractionmass=0.7);
   AirE3mbar->AddElement(O, fractionmass=0.3);
 
   G4Material* AirE4mbar = 
-    new G4Material("AirE4mbar"  , density= 1.290e-7*mg/cm3, ncomponents=2);
+    new G4Material("AirE4mbar"  , density= 1.290e-7*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   AirE4mbar->AddElement(N, fractionmass=0.7);
   AirE4mbar->AddElement(O, fractionmass=0.3);
 
   G4Material* AirE5mbar = 
-    new G4Material("AirE5mbar"  , density= 1.290e-8*mg/cm3, ncomponents=2);
+    new G4Material("AirE5mbar"  , density= 1.290e-8*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   AirE5mbar->AddElement(N, fractionmass=0.7);
   AirE5mbar->AddElement(O, fractionmass=0.3);
 
   G4Material* AirE6mbar = 
-    new G4Material("AirE6mbar"  , density= 1.290e-9*mg/cm3, ncomponents=2);
+    new G4Material("AirE6mbar"  , density= 1.290e-9*CLHEP::mg/CLHEP::cm3, ncomponents=2);
   AirE6mbar->AddElement(N, fractionmass=0.7);
   AirE6mbar->AddElement(O, fractionmass=0.3);
 
@@ -1701,34 +1701,34 @@ void musrDetectorConstruction::DefineMaterials()
 //
 
 //  G4Material* Vacuum =
-  new G4Material("Vacuum", z=1., a=1.01*g/mole,density= universe_mean_density,
-		   kStateGas, 2.73*kelvin, 3.e-18*pascal);
+  new G4Material("Vacuum", z=1., a=1.01*CLHEP::g/CLHEP::mole,density= CLHEP::universe_mean_density,
+		   kStateGas, 2.73*CLHEP::kelvin, 3.e-18*CLHEP::pascal);
   
-  new G4Material("ArgonGas", z= 18., a= 39.95*g/mole, density= 0.00000000001*mg/cm3);
+  new G4Material("ArgonGas", z= 18., a= 39.95*CLHEP::g/CLHEP::mole, density= 0.00000000001*CLHEP::mg/CLHEP::cm3);
 
-  new G4Material("HeliumGas5mbar", z=2., a=4.002602*g/mole, density= 0.00000088132*g/cm3);
-  new G4Material("HeliumGas6mbar", z=2., a=4.002602*g/mole, density= 0.000001057584*g/cm3);
-  new G4Material("HeliumGas7mbar", z=2., a=4.002602*g/mole, density= 0.000001233848*g/cm3);
-  new G4Material("HeliumGas8mbar", z=2., a=4.002602*g/mole, density= 0.000001410112*g/cm3);
-  new G4Material("HeliumGas9mbar", z=2., a=4.002602*g/mole, density= 0.000001586376*g/cm3);
-  new G4Material("HeliumGas10mbar",z=2., a=4.002602*g/mole, density= 0.00000176264*g/cm3);
-  new G4Material("HeliumGas11mbar",z=2., a=4.002602*g/mole, density= 0.000001938904*g/cm3);
-  new G4Material("HeliumGas12mbar",z=2., a=4.002602*g/mole, density= 0.000002115168*g/cm3);
-  new G4Material("HeliumGas13mbar",z=2., a=4.002602*g/mole, density= 0.000002291432*g/cm3);
-  new G4Material("HeliumGas14mbar",z=2., a=4.002602*g/mole, density= 0.000002467696*g/cm3);
-  new G4Material("HeliumGas15mbar",z=2., a=4.002602*g/mole, density= 0.00000264396*g/cm3);
-  new G4Material("HeliumGas20mbar",z=2., a=4.002602*g/mole, density= 0.00000352528*g/cm3);
-  new G4Material("HeliumGas30mbar",z=2., a=4.002602*g/mole, density= 0.00000528792*g/cm3);
-  new G4Material("HeliumGas40mbar",z=2., a=4.002602*g/mole, density= 0.00000705056*g/cm3);
-  new G4Material("HeliumGas50mbar",z=2., a=4.002602*g/mole, density= 0.00000881320*g/cm3);
-  new G4Material("HeliumGas60mbar",z=2., a=4.002602*g/mole, density= 0.00001057584*g/cm3);
-  new G4Material("HeliumGas70mbar",z=2., a=4.002602*g/mole, density= 0.00001233848*g/cm3);
-  new G4Material("HeliumGas80mbar",z=2., a=4.002602*g/mole, density= 0.00001410112*g/cm3);
-  new G4Material("HeliumGas90mbar",z=2., a=4.002602*g/mole, density= 0.00001586376*g/cm3);
-  new G4Material("HeliumGas100mbar",z=2.,a=4.002602*g/mole, density= 0.00001762640*g/cm3);
-  new G4Material("HeliumGasSat4K",z=2.,  a=4.002602*g/mole, density= 0.016891*g/cm3); // saturated vapour above liquid at 4.2K (JSL)
-  new G4Material("HeliumGas5mbar4K",z=2.,a=4.002602*g/mole, density= 0.016891*5.0/1013.0*g/cm3); // typical cold exchange gas, 4.2K and 5 mbar
-  new G4Material("HeliumGas2mbar4K",z=2.,a=4.002602*g/mole, density= 0.016891*2.0/1013.0*g/cm3); // typical cold exchange gas, 4.2K and 5 mbar
+  new G4Material("HeliumGas5mbar", z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000088132*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas6mbar", z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000001057584*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas7mbar", z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000001233848*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas8mbar", z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000001410112*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas9mbar", z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000001586376*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas10mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000176264*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas11mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000001938904*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas12mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000002115168*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas13mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000002291432*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas14mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.000002467696*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas15mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000264396*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas20mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000352528*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas30mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000528792*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas40mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000705056*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas50mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00000881320*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas60mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00001057584*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas70mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00001233848*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas80mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00001410112*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas90mbar",z=2., a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00001586376*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGas100mbar",z=2.,a=4.002602*CLHEP::g/CLHEP::mole, density= 0.00001762640*CLHEP::g/CLHEP::cm3);
+  new G4Material("HeliumGasSat4K",z=2.,  a=4.002602*CLHEP::g/CLHEP::mole, density= 0.016891*CLHEP::g/CLHEP::cm3); // saturated vapour above liquid at 4.2K (JSL)
+  new G4Material("HeliumGas5mbar4K",z=2.,a=4.002602*CLHEP::g/CLHEP::mole, density= 0.016891*5.0/1013.0*CLHEP::g/CLHEP::cm3); // typical cold exchange gas, 4.2K and 5 mbar
+  new G4Material("HeliumGas2mbar4K",z=2.,a=4.002602*CLHEP::g/CLHEP::mole, density= 0.016891*2.0/1013.0*CLHEP::g/CLHEP::cm3); // typical cold exchange gas, 4.2K and 5 mbar
 
   if (musrParameters::boolG4OpticalPhotons) {
 
diff --git a/src/musrUniformField.cc b/src/musrUniformField.cc
index 28b799b..f1f9aee 100644
--- a/src/musrUniformField.cc
+++ b/src/musrUniformField.cc
@@ -45,8 +45,8 @@ musrUniformField::musrUniformField(G4double EMF[6], G4double half_X, G4double ha
   G4String volName = lv->GetName().substr(4);
   G4cout << "\n ---> EM field in volume " << volName << " set to:" << G4endl;
   printf   ("      B = (%0.3g, %0.3g, %0.3g) T,  E = (%0.3g, %0.3g, %0.3g) kV/mm\n",
-  EMF[0]/tesla,         EMF[1]/tesla,         EMF[2]/tesla, 
-  EMF[3]/(kilovolt/mm), EMF[4]/(kilovolt/mm), EMF[5]/(kilovolt/mm));
+  EMF[0]/CLHEP::tesla,         EMF[1]/CLHEP::tesla,         EMF[2]/CLHEP::tesla, 
+  EMF[3]/(CLHEP::kilovolt/CLHEP::mm), EMF[4]/(CLHEP::kilovolt/CLHEP::mm), EMF[5]/(CLHEP::kilovolt/CLHEP::mm));
   G4cout << "-----------------------------------------------------------" << G4endl;
 }