musrAtRestSpinRotation.cc

Go to the documentation of this file.

#include "musrAtRestSpinRotation.hh"

#include "G4ThreeVector.hh"
#include "G4MagneticField.hh"
#include "G4Track.hh"
#include "G4FieldManager.hh"
#include <iomanip>
#include "globals.hh"
#include "G4ParticleChange.hh"
#include "G4ios.hh"
#include "G4Transform3D.hh"
#include "G4UnitsTable.hh"
#include "G4TransportationManager.hh" 
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4VParticleChange.hh"
#include "G4ProcessVector.hh"
#include "G4ProcessManager.hh"

musrAtRestSpinRotation::musrAtRestSpinRotation(const G4String& processName)
  : G4VRestProcess (processName)
{
 G4cout << GetProcessName() << " is created "<< G4endl;
 pointer = this;
}


musrAtRestSpinRotation::~musrAtRestSpinRotation() 
{;} 


musrAtRestSpinRotation*  musrAtRestSpinRotation::pointer=0;
musrAtRestSpinRotation*  musrAtRestSpinRotation::GetInstance()
{
  return pointer;
}




G4VParticleChange* musrAtRestSpinRotation::AtRestDoIt(const G4Track& theTrack, const G4Step& aStep)
{
  
  theParticleChange.Initialize(theTrack);  
  theParticleChange.ProposeTrackStatus(fAlive);
  // G4cout<<"AT REST::: globaltime  "<< theTrack.GetGlobalTime()/ns <<G4endl;// getchar();
  theParticleChange.ProposePosition(theTrack.GetPosition());
  theParticleChange.ProposeMomentumDirection(theTrack.GetMomentumDirection());
  theParticleChange.ProposeEnergy(theTrack.GetKineticEnergy());
  
  // tao :: Get Time
  itime = theTrack.GetProperTime();
  G4double gtime = theTrack.GetGlobalTime();
  ftime = theTrack.GetDynamicParticle()->GetPreAssignedDecayProperTime(); 


  deltatime = ftime - itime;
  theParticleChange.ProposeGlobalTime(deltatime + itime -gtime);
  theParticleChange.ProposeProperTime(deltatime);
  
  polar = aStep.GetTrack()->GetPolarization();


  
  
  //  if(aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetFieldManager())
  if(G4TransportationManager::GetTransportationManager()->GetFieldManager())
    {
      //      G4FieldManager *fMgr = theTrack.GetVolume()->GetLogicalVolume()->GetFieldManager();
      G4FieldManager *fMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
      


      if(!fMgr->DoesFieldChangeEnergy())//then we have a magnetic field
        {
          

          //G4cout<<"AT REST::: MAGNETIC FIELD HERE" <<G4endl;// getchar();
          
          
          // tao :: Get Field
          point[0]=theTrack.GetPosition().x();
          point[1]=theTrack.GetPosition().y();
          point[2]=theTrack.GetPosition().z();
              
          mfield = fMgr->GetDetectorField();       

          mfield->GetFieldValue(point,B);
          

          //G4cout <<"AT REST::: MAGNETIC FIELD B="<< B[0] <<" " << B[1] <<" " << B[2] <<G4endl;
          

              
#ifdef G4SRVERBOSE
          G4cout <<"AT REST::: TIME= proper: "<< itime <<" ;  global: " << gtime  <<" decay: " << ftime <<G4endl;
#endif
          
          
          G4ThreeVector magField(B[0],B[1],B[2]);
          RotateSpin(aStep,magField,deltatime);
          
          
#ifdef G4SRVERBOSE
          G4cout<<"AT REST::: spin rotated";
#endif
          
        } 


          theParticleChange.ProposePolarization(polar);
      

      
#ifdef G4SRVERBOSE
      theParticleChange.DumpInfo();
#endif
      
      
      return &theParticleChange;
      
      // then the AtRestUpdateStep method is called cf LEMuSRParticleChangeForSpinRotation (SR)
      
    }
  
  else 
    {
      //  G4cout <<"No MAGNETIC FIELD do not call AtRestPrecession :: RotateSpin!!! "; 
      return &theParticleChange;



    }
  
  

  
}



void musrAtRestSpinRotation::RotateSpin( const G4Step& aStep, G4ThreeVector B, G4double deltatime )
{
  
  
  G4Transform3D Spin_rotation;
 
  G4double Bnorm = sqrt(sqr(B[0])  + sqr(B[1]) +sqr(B[2]) );


  //  G4cout<< "AT REST::: PARAMETERS\n" 
  //    << "Magnetic Field Norm  : " << G4BestUnit(Bnorm,"Magnetic flux density") <<"\n";
    
  G4double omega,q,a,fqz;
  
  q= aStep.GetTrack()->GetDefinition()->GetPDGCharge();
  a= 1.165922e-3;
  fqz = 8.5062e+7*rad/(s*kilogauss);
  
#ifdef G4SRVERBOSE
  G4cout<< "AT REST::: PARAMETERS\n" 
        << "Charge  : " << q <<"\n";
#endif
  
  omega= - (fqz)*(1.+a) * Bnorm;

#ifdef G4SRVERBOSE
  G4cout<< "AT REST::: PARAMETERS\n" 
        << "Frequency: " << G4BestUnit(omega/(2*M_PI*rad),"Frequency") <<"\n";
#endif
   
  
  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<< "AT REST::: PARAMETERS\n" 
  << "Initial spin  : " << spin <<"\n"
  << "Delta time    : " << deltatime <<"\n"
  << "Rotation angle: " << rotation_angle/rad <<"\n"
  << "New spin      : " << newspin <<"\n"
  << "Checked norms : " << x <<" " <<y <<" \n"
  << G4endl; 
  #endif
  
  polar = newspin;
  
}

---