G4DecayWithSpin.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // ------------------------------------------------------------
00027 //      GEANT 4 class header file
00028 //
00029 //      History:
00030 //      17 August 2004  P. Gumplinger, T. MacPhail
00031 //      11 April  2008  Kamil Sedlak (PSI), Toni Shiroka (PSI)
00032 // ------------------------------------------------------------
00033 //
00034 #include "G4DecayWithSpin.hh"
00035 
00036 #include "G4Step.hh"
00037 #include "G4Track.hh"
00038 #include "G4DecayTable.hh"
00039 #include "G4MuonDecayChannelWithSpin.hh"
00040 
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4Vector3D.hh"
00044 
00045 #include "G4TransportationManager.hh"
00046 #include "G4PropagatorInField.hh"
00047 #include "G4FieldManager.hh"
00048 #include "G4Field.hh"
00049 
00050 #include "G4Transform3D.hh"
00051 
00052 G4DecayWithSpin::G4DecayWithSpin(const G4String& processName):G4Decay(processName)
00053 {
00054   // set Process Sub Type   
00055   SetProcessSubType(static_cast<int>(DECAY_WithSpin));
00056 
00057 }
00058 
00059 G4DecayWithSpin::~G4DecayWithSpin(){}
00060 
00061 G4VParticleChange* G4DecayWithSpin::DecayIt(const G4Track& aTrack, const G4Step& aStep)
00062 {
00063 
00064 // get particle
00065   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00066   const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00067 
00068 // get parent_polarization
00069   G4ThreeVector parent_polarization = aParticle->GetPolarization();
00070 
00071   if(parent_polarization == G4ThreeVector(0,0,0))
00072   {
00073     // Generate random polarization direction
00074 
00075     G4double cost = 1. - 2.*G4UniformRand();
00076     G4double sint = std::sqrt((1.-cost)*(1.+cost));
00077 
00078     G4double phi = twopi*G4UniformRand();
00079     G4double sinp = std::sin(phi);
00080     G4double cosp = std::cos(phi);
00081 
00082     G4double px = sint*cosp;
00083     G4double py = sint*sinp;
00084     G4double pz = cost;
00085 
00086     parent_polarization.setX(px);
00087     parent_polarization.setY(py);
00088     parent_polarization.setZ(pz);
00089 
00090   }else{
00091 
00092     G4FieldManager* fieldMgr = aStep.GetTrack()->GetVolume()->
00093                                      GetLogicalVolume()->GetFieldManager();
00094 
00095     if (!fieldMgr) {
00096        G4TransportationManager *transportMgr =
00097                          G4TransportationManager::GetTransportationManager();
00098        G4PropagatorInField* fFieldPropagator = 
00099                                         transportMgr->GetPropagatorInField();
00100        if (fFieldPropagator) fieldMgr = 
00101                                   fFieldPropagator->GetCurrentFieldManager();
00102     }
00103 
00104     const G4Field* field = NULL;
00105     if(fieldMgr)field = fieldMgr->GetDetectorField();
00106 
00107     if (field) {
00108 
00109        G4double point[4];
00110        point[0] = (aStep.GetPostStepPoint()->GetPosition())[0];
00111        point[1] = (aStep.GetPostStepPoint()->GetPosition())[1];
00112        point[2] = (aStep.GetPostStepPoint()->GetPosition())[2];
00113        point[3] = aTrack.GetGlobalTime();
00114 
00115        G4double fieldValue[6];
00116        field -> GetFieldValue(point,fieldValue);
00117 
00118        G4ThreeVector B(fieldValue[0],fieldValue[1],fieldValue[2]);
00119 
00120        // Call the spin precession only for non-zero mag. field
00121        if (B.mag2() > 0.) parent_polarization = 
00122                                  Spin_Precession(aStep,B,fRemainderLifeTime);
00123 
00124     }
00125   }
00126 
00127 // decay table
00128   G4DecayTable *decaytable = aParticleDef->GetDecayTable();
00129 
00130   if (decaytable) {
00131      G4MuonDecayChannelWithSpin *decaychannel;
00132      decaychannel = (G4MuonDecayChannelWithSpin*)decaytable->SelectADecayChannel();
00133      if (decaychannel) decaychannel->SetPolarization(parent_polarization);
00134   }
00135 
00136   G4ParticleChangeForDecay* pParticleChangeForDecay;
00137 
00138   pParticleChangeForDecay = (G4ParticleChangeForDecay*)G4Decay::DecayIt(aTrack,aStep);
00139 
00140   pParticleChangeForDecay->ProposePolarization(parent_polarization);
00141 
00142   return pParticleChangeForDecay;
00143 
00144 }
00145 
00146 G4ThreeVector G4DecayWithSpin::Spin_Precession( const G4Step& aStep,
00147                                         G4ThreeVector B, G4double deltatime )
00148 {
00149   G4double Bnorm = std::sqrt(sqr(B[0])  + sqr(B[1]) +sqr(B[2]) );
00150 
00151   G4double q = aStep.GetTrack()->GetDefinition()->GetPDGCharge();
00152   G4double a = 1.165922e-3;
00153   G4double s_omega = 8.5062e+7*rad/(s*kilogauss);
00154 
00155   G4double omega = -(q*s_omega)*(1.+a) * Bnorm;
00156 
00157   G4double rotationangle = deltatime * omega;
00158 
00159   G4Transform3D SpinRotation = G4Rotate3D(rotationangle,B.unit());
00160 
00161   G4Vector3D Spin = aStep.GetTrack() -> GetPolarization();
00162 
00163   G4Vector3D newSpin = SpinRotation * Spin;
00164 
00165 #ifdef G4VERBOSE
00166   if (GetVerboseLevel()>2) {
00167     G4double normspin = std::sqrt(Spin*Spin);
00168     G4double normnewspin = std::sqrt(newSpin*newSpin);
00169     //G4double cosalpha = Spin*newSpin/normspin/normnewspin;
00170     //G4double alpha = std::acos(cosalpha);
00171 
00172     G4cout << "AT REST::: PARAMETERS " << G4endl;
00173     G4cout << "Initial spin  : " << Spin  << G4endl;
00174     G4cout << "Delta time    : " << deltatime  << G4endl;
00175     G4cout << "Rotation angle: " << rotationangle/rad  << G4endl;
00176     G4cout << "New spin      : " << newSpin  << G4endl;
00177     G4cout << "Checked norms : " << normspin <<" " << normnewspin << G4endl;
00178   }
00179 #endif
00180 
00181   return newSpin;
00182 
00183 }

Generated on Mon May 27 17:48:00 2013 for Geant4 by  doxygen 1.4.7