diff --git a/geant4/LEMuSR/src/LEMuSRMUONIUMScatt.cc b/geant4/LEMuSR/src/LEMuSRMUONIUMScatt.cc
index 06c5eaf..ef7bfa7 100644
--- a/geant4/LEMuSR/src/LEMuSRMUONIUMScatt.cc
+++ b/geant4/LEMuSR/src/LEMuSRMUONIUMScatt.cc
@@ -4,6 +4,10 @@
 #include "G4TransportationManager.hh"
 #include "Randomize.hh"
 #include "G4ProductionCutsTable.hh"
+#include "G4FieldManager.hh"
+#include "G4Field.hh"
+
+
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     
@@ -26,11 +30,11 @@ G4VParticleChange* LEMuSRMUONIUMScatt::PostStepDoIt(
   fParticleChange.Initialize(trackData);
   
   // tao :: Get Time
-  G4double itime = trackData.GetProperTime();
-  G4double gtime = trackData.GetGlobalTime();
-  G4double ftime = trackData.GetDynamicParticle()->GetPreAssignedDecayProperTime(); 
+  itime = trackData.GetProperTime();
+  gtime = trackData.GetGlobalTime();
+  ftime = trackData.GetDynamicParticle()->GetPreAssignedDecayProperTime(); 
   
-  G4double deltatime = ftime - itime;
+  deltatime = ftime - itime;
   fParticleChange.ProposeGlobalTime(deltatime + itime -gtime);
   
   // set position momentum energy
@@ -46,15 +50,18 @@ G4VParticleChange* LEMuSRMUONIUMScatt::PostStepDoIt(
       //      fParticleChange.ProposeEnergy(0);
    fParticleChange.ProposePosition(trackData.GetStep()->GetPreStepPoint()->GetPosition());
     
-      fParticleChange.ProposePolarization(trackData.GetPolarization());
+   //      fParticleChange.ProposePolarization(trackData.GetPolarization());
       //      G4cout<<" at rest "<< trackData.GetStep()->GetPreStepPoint()->GetPosition() <<G4endl;
   fParticleChange.ProposeTrackStatus(fStopButAlive) ;
       
-    }  
+    }   
   else
     {
-      fParticleChange.ProposePolarization(trackData.GetPolarization());
+      //     fParticleChange.ProposePolarization(trackData.GetPolarization());
     }
+
+  polar = RotateSpinIfMag(trackData, aStep);
+  fParticleChange.ProposePolarization(polar);
   //   fParticleChange.DumpInfo();
   return &fParticleChange;
 
@@ -78,7 +85,8 @@ G4bool LEMuSRMUONIUMScatt::CheckCondition( const G4Step& aStep)
   return condition;
  
 }  
-  
+
+ 
 
 G4double LEMuSRMUONIUMScatt:: GetMeanFreePath(const G4Track& ,
 			   G4double  ,
@@ -123,3 +131,134 @@ void LEMuSRMUONIUMScatt::PrepareSecondary(const G4Track& track)
 
 
 ;}
+
+
+G4ThreeVector LEMuSRMUONIUMScatt::RotateSpinIfMag( const G4Track& theTrack, const G4Step& aStep)
+{
+  G4ThreeVector theSpin = theTrack.GetPolarization() ;
+
+  if(theTrack.GetKineticEnergy()!=0&&aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetFieldManager())
+    {
+      
+      
+      G4FieldManager *fMgr = theTrack.GetVolume()->GetLogicalVolume()->GetFieldManager();
+      
+      
+      //  if(!fMgr->DoesFieldChangeEnergy())//then we have a magnetic field
+      if(fMgr->FieldHasMagComponent())//then we have a magnetic field
+	{
+	  
+#ifdef G4SRVERBOSE
+	  G4cout<<"IN FLIGHT::: MAGNETIC FIELD HERE" <<G4endl;// getchar();
+#endif
+	  
+	  // tao :: Get Field
+	  point[0]=theTrack.GetPosition().x();
+	  point[1]=theTrack.GetPosition().y();
+	  point[2]=theTrack.GetPosition().z();
+	 
+
+	  const G4Field* mfield;
+	  mfield = fMgr->GetDetectorField();       
+	  
+	  mfield->GetFieldValue(point,B);
+	  
+#ifdef G4SRVERBOSE
+	  G4cout <<"IN FLIGHT::: MAGNETIC FIELD B="<< B[0]/gauss <<" G, " << B[1]/gauss <<" G, " << B[2]/gauss<<" G" <<G4endl;
+	  
+#endif
+	      
+#ifdef G4SRVERBOSE
+	  G4cout <<"IN FLIGHT::: TIME= proper: "<< itime <<" ;  global: " << gtime  <<" decay: " << ftime <<G4endl;
+#endif
+	  
+	  
+	  G4ThreeVector magField(B[0],B[1],B[2]);
+	theSpin= RotateSpin(aStep,magField,deltatime);
+	  
+	  
+#ifdef G4SRVERBOSE
+	  G4cout<<"IN FLIGHT::: spin rotated";
+#endif
+
+	}
+    }
+
+	  return theSpin;
+}
+
+
+  
+
+
+
+G4ThreeVector LEMuSRMUONIUMScatt::RotateSpin( const G4Step& aStep, G4ThreeVector B, G4double deltatime )
+{
+  
+  
+  G4Transform3D Spin_rotation;
+ 
+  G4double Bnorm = sqrt(sqr(B[0])  + sqr(B[1]) +sqr(B[2]) );
+
+#ifdef G4SRVERBOSE
+  G4cout<< "IN FLIGHT::: PARAMETERS\n" 
+	<< "Magnetic Field Norm  : " << G4BestUnit(Bnorm,"Magnetic flux density") <<"\n";
+#endif
+    
+  G4double omega,q,a,fqz;
+  G4double gamma;
+  q= aStep.GetTrack()->GetDefinition()->GetPDGCharge();
+  a= 1.165922e-3;
+  fqz = 8.5062e+7*rad/(s*kilogauss);
+  
+  gamma = aStep.GetTrack()->GetDefinition()->GetGyromagneticRatio()*rad;
+
+  //  G4cout<<  fqz*(s*tesla)<<G4endl;
+  //  G4cout<<  gamma*(s*tesla)<<G4endl;
+ 
+
+
+#ifdef G4SRVERBOSE
+  G4cout<< "IN FLIGHT::: PARAMETERS\n" 
+	<< "Charge  : " << q <<"\n";
+#endif
+  
+  //  omega= - (fqz)*(1.+a) * Bnorm;
+  omega= - (gamma) * Bnorm;
+
+#ifdef G4SRVERBOSE
+  G4cout<< "IN FLIGHT::: PARAMETERS\n" 
+	<< "Frequency: " << G4BestUnit(fqz*gauss/(2*M_PI*rad),"Frequency") <<"\n";
+  G4cout<< "IN FLIGHT::: PARAMETERS\n" 
+	<< "FrequencyG: " << G4BestUnit(gamma*gauss/(2*M_PI*rad),"Frequency") <<"\n";  
+#endif
+   
+  
+  G4double rotation_angle = deltatime*omega; 
+  
+
+  Spin_rotation= G4Rotate3D(rotation_angle,B/Bnorm);
+  
+  HepVector3D spin = aStep.GetTrack()->GetPolarization();
+  HepVector3D newspin;
+  newspin = Spin_rotation*spin;
+  
+  G4double x,y,z,alpha;
+  x = sqrt(spin*spin);
+  y = sqrt(newspin*newspin);
+  z = spin*newspin/x/y;
+  alpha = acos(z);
+  
+#ifdef G4SRVERBOSE
+  G4cout<< "IN FLIGHT::: PARAMETERS\n" 
+	<< "Initial spin  : " << spin <<"\n"
+	<< "Delta time    : " << deltatime <<"\n"
+	<< "Rotation angle: " << rotation_angle/(M_PI*rad) <<"\n"
+	<< "New spin      : " << newspin <<"\n"
+	<< "Checked norms : " << x <<" " <<y <<" \n"
+	<< G4endl; 
+#endif
+  
+  return newspin; 
+}
+