Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/decay/src/G4DecayWithSpin.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/decay/src/G4DecayWithSpin.cc (Version 11.3.0) and /processes/decay/src/G4DecayWithSpin.cc (Version 9.6.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // -------------------------------------------     26 // ------------------------------------------------------------
 27 //      GEANT 4 class header file                  27 //      GEANT 4 class header file
 28 //                                                 28 //
 29 //      History:                                   29 //      History:
 30 //      17 August 2004  P. Gumplinger, T. MacP     30 //      17 August 2004  P. Gumplinger, T. MacPhail
 31 //      11 April  2008  Kamil Sedlak (PSI), To     31 //      11 April  2008  Kamil Sedlak (PSI), Toni Shiroka (PSI)
 32 // -------------------------------------------     32 // ------------------------------------------------------------
 33 //                                                 33 //
 34 #include "G4DecayWithSpin.hh"                      34 #include "G4DecayWithSpin.hh"
 35                                                    35 
 36 #include "G4Step.hh"                               36 #include "G4Step.hh"
 37 #include "G4Track.hh"                              37 #include "G4Track.hh"
 38 #include "G4DecayTable.hh"                         38 #include "G4DecayTable.hh"
 39 #include "G4MuonDecayChannelWithSpin.hh"           39 #include "G4MuonDecayChannelWithSpin.hh"
 40                                                    40 
 41 #include "G4PhysicalConstants.hh"                  41 #include "G4PhysicalConstants.hh"
 42 #include "G4SystemOfUnits.hh"                      42 #include "G4SystemOfUnits.hh"
 43 #include "G4Vector3D.hh"                           43 #include "G4Vector3D.hh"
 44                                                    44 
 45 #include "G4TransportationManager.hh"              45 #include "G4TransportationManager.hh"
 46 #include "G4PropagatorInField.hh"                  46 #include "G4PropagatorInField.hh"
 47 #include "G4FieldManager.hh"                       47 #include "G4FieldManager.hh"
 48 #include "G4Field.hh"                              48 #include "G4Field.hh"
 49                                                    49 
 50 #include "G4Transform3D.hh"                        50 #include "G4Transform3D.hh"
 51                                                    51 
 52 G4DecayWithSpin::G4DecayWithSpin(const G4Strin     52 G4DecayWithSpin::G4DecayWithSpin(const G4String& processName):G4Decay(processName)
 53 {                                                  53 {
 54   // set Process Sub Type                          54   // set Process Sub Type   
 55   SetProcessSubType(static_cast<int>(DECAY_Wit     55   SetProcessSubType(static_cast<int>(DECAY_WithSpin));
 56                                                    56 
 57 }                                                  57 }
 58                                                    58 
 59 G4DecayWithSpin::~G4DecayWithSpin(){}              59 G4DecayWithSpin::~G4DecayWithSpin(){}
 60                                                    60 
 61 G4VParticleChange* G4DecayWithSpin::PostStepDo <<  61 G4VParticleChange* G4DecayWithSpin::DecayIt(const G4Track& aTrack, const G4Step& aStep)
 62 {                                                  62 {
 63   if ( (aTrack.GetTrackStatus() == fStopButAli << 
 64        (aTrack.GetTrackStatus() == fStopAndKil << 
 65     fParticleChangeForDecay.Initialize(aTrack) << 
 66     return &fParticleChangeForDecay;           << 
 67   }                                            << 
 68                                                    63 
 69 // get particle                                    64 // get particle
 70   const G4DynamicParticle* aParticle = aTrack.     65   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
 71   const G4ParticleDefinition* aParticleDef = a     66   const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
 72                                                    67 
 73 // get parent_polarization                         68 // get parent_polarization
 74   G4ThreeVector parent_polarization = aParticl     69   G4ThreeVector parent_polarization = aParticle->GetPolarization();
 75                                                    70 
 76   if(parent_polarization == G4ThreeVector(0,0, <<  71   if(parent_polarization == G4ThreeVector(0,0,0))
                                                   >>  72   {
 77     // Generate random polarization direction      73     // Generate random polarization direction
 78     G4double cost = 1. - 2.*G4UniformRand();   << 
 79     G4double sint = std::sqrt((1.-cost)*(1.+co << 
 80                                                << 
 81     G4double phi = twopi*G4UniformRand();      << 
 82     G4double sinp = std::sin(phi);             << 
 83     G4double cosp = std::cos(phi);             << 
 84                                                << 
 85     G4double px = sint*cosp;                   << 
 86     G4double py = sint*sinp;                   << 
 87     G4double pz = cost;                        << 
 88                                                << 
 89     parent_polarization.setX(px);              << 
 90     parent_polarization.setY(py);              << 
 91     parent_polarization.setZ(pz);              << 
 92   }                                            << 
 93                                                << 
 94 // decay table                                 << 
 95   G4DecayTable *decaytable = aParticleDef->Get << 
 96   if (decaytable != nullptr) {                 << 
 97     for (G4int ip=0; ip<decaytable->entries(); << 
 98       decaytable->GetDecayChannel(ip)->SetPola << 
 99     }                                          << 
100   }                                            << 
101                                                << 
102   G4ParticleChangeForDecay* pParticleChangeFor << 
103   pParticleChangeForDecay = (G4ParticleChangeF << 
104   pParticleChangeForDecay->ProposePolarization << 
105                                                << 
106   return pParticleChangeForDecay;              << 
107 }                                              << 
108                                                << 
109 G4VParticleChange* G4DecayWithSpin::AtRestDoIt << 
110 {                                              << 
111                                                << 
112 // get particle                                << 
113   const G4DynamicParticle* aParticle = aTrack. << 
114   const G4ParticleDefinition* aParticleDef = a << 
115                                                << 
116 // get parent_polarization                     << 
117   G4ThreeVector parent_polarization = aParticl << 
118                                                    74 
119   if(parent_polarization == G4ThreeVector(0,0, << 
120     // Generate random polarization direction  << 
121     G4double cost = 1. - 2.*G4UniformRand();       75     G4double cost = 1. - 2.*G4UniformRand();
122     G4double sint = std::sqrt((1.-cost)*(1.+co     76     G4double sint = std::sqrt((1.-cost)*(1.+cost));
123                                                    77 
124     G4double phi = twopi*G4UniformRand();          78     G4double phi = twopi*G4UniformRand();
125     G4double sinp = std::sin(phi);                 79     G4double sinp = std::sin(phi);
126     G4double cosp = std::cos(phi);                 80     G4double cosp = std::cos(phi);
127                                                    81 
128     G4double px = sint*cosp;                       82     G4double px = sint*cosp;
129     G4double py = sint*sinp;                       83     G4double py = sint*sinp;
130     G4double pz = cost;                            84     G4double pz = cost;
131                                                    85 
132     parent_polarization.setX(px);                  86     parent_polarization.setX(px);
133     parent_polarization.setY(py);                  87     parent_polarization.setY(py);
134     parent_polarization.setZ(pz);                  88     parent_polarization.setZ(pz);
135                                                    89 
136   }else{                                           90   }else{
137                                                    91 
138     G4FieldManager* fieldMgr = aStep.GetTrack(     92     G4FieldManager* fieldMgr = aStep.GetTrack()->GetVolume()->
139                                      GetLogica     93                                      GetLogicalVolume()->GetFieldManager();
140     if (fieldMgr == nullptr) {                 <<  94 
                                                   >>  95     if (!fieldMgr) {
141        G4TransportationManager *transportMgr =     96        G4TransportationManager *transportMgr =
142                          G4TransportationManag     97                          G4TransportationManager::GetTransportationManager();
143        G4PropagatorInField* fFieldPropagator =     98        G4PropagatorInField* fFieldPropagator = 
144                                         transp     99                                         transportMgr->GetPropagatorInField();
145        if (fFieldPropagator) fieldMgr =           100        if (fFieldPropagator) fieldMgr = 
146                                   fFieldPropag    101                                   fFieldPropagator->GetCurrentFieldManager();
147     }                                             102     }
148                                                   103 
149     const G4Field* field = nullptr;            << 104     const G4Field* field = NULL;
150     if (fieldMgr != nullptr) field = fieldMgr- << 105     if(fieldMgr)field = fieldMgr->GetDetectorField();
                                                   >> 106 
                                                   >> 107     if (field) {
151                                                   108 
152     if ( field != nullptr ) {                  << 
153        G4double point[4];                         109        G4double point[4];
154        point[0] = (aStep.GetPostStepPoint()->G    110        point[0] = (aStep.GetPostStepPoint()->GetPosition())[0];
155        point[1] = (aStep.GetPostStepPoint()->G    111        point[1] = (aStep.GetPostStepPoint()->GetPosition())[1];
156        point[2] = (aStep.GetPostStepPoint()->G    112        point[2] = (aStep.GetPostStepPoint()->GetPosition())[2];
157        point[3] = aTrack.GetGlobalTime();         113        point[3] = aTrack.GetGlobalTime();
158                                                   114 
159        G4double fieldValue[6] ={ 0., 0., 0., 0 << 115        G4double fieldValue[6];
160        field -> GetFieldValue(point,fieldValue    116        field -> GetFieldValue(point,fieldValue);
                                                   >> 117 
161        G4ThreeVector B(fieldValue[0],fieldValu    118        G4ThreeVector B(fieldValue[0],fieldValue[1],fieldValue[2]);
162                                                   119 
163        // Call the spin precession only for no    120        // Call the spin precession only for non-zero mag. field
164        if (B.mag2() > 0.) parent_polarization     121        if (B.mag2() > 0.) parent_polarization = 
165                                  Spin_Precessi    122                                  Spin_Precession(aStep,B,fRemainderLifeTime);
166                                                   123 
167     }                                             124     }
168   }                                               125   }
169                                                   126 
170 // decay table                                    127 // decay table
171   G4DecayTable *decaytable = aParticleDef->Get    128   G4DecayTable *decaytable = aParticleDef->GetDecayTable();
172   if ( decaytable != nullptr) {                << 129 
173     for (G4int ip=0; ip<decaytable->entries(); << 130   if (decaytable) {
174       decaytable->GetDecayChannel(ip)->SetPola << 131      G4MuonDecayChannelWithSpin *decaychannel;
175     }                                          << 132      decaychannel = (G4MuonDecayChannelWithSpin*)decaytable->SelectADecayChannel();
176   }                                            << 133      if (decaychannel) decaychannel->SetPolarization(parent_polarization);
                                                   >> 134   }
177                                                   135 
178   G4ParticleChangeForDecay* pParticleChangeFor    136   G4ParticleChangeForDecay* pParticleChangeForDecay;
                                                   >> 137 
179   pParticleChangeForDecay = (G4ParticleChangeF    138   pParticleChangeForDecay = (G4ParticleChangeForDecay*)G4Decay::DecayIt(aTrack,aStep);
                                                   >> 139 
180   pParticleChangeForDecay->ProposePolarization    140   pParticleChangeForDecay->ProposePolarization(parent_polarization);
181                                                   141 
182   return pParticleChangeForDecay;                 142   return pParticleChangeForDecay;
183                                                   143 
184 }                                                 144 }
185                                                   145 
186 G4ThreeVector G4DecayWithSpin::Spin_Precession    146 G4ThreeVector G4DecayWithSpin::Spin_Precession( const G4Step& aStep,
187                                         G4Thre    147                                         G4ThreeVector B, G4double deltatime )
188 {                                                 148 {
189   G4double Bnorm = std::sqrt(sqr(B[0])  + sqr(    149   G4double Bnorm = std::sqrt(sqr(B[0])  + sqr(B[1]) +sqr(B[2]) );
190                                                   150 
191   G4double q = aStep.GetTrack()->GetDefinition    151   G4double q = aStep.GetTrack()->GetDefinition()->GetPDGCharge();
192   G4double a = 1.165922e-3;                       152   G4double a = 1.165922e-3;
193   G4double s_omega = 8.5062e+7*rad/(s*kilogaus    153   G4double s_omega = 8.5062e+7*rad/(s*kilogauss);
194                                                   154 
195   G4double omega = -(q*s_omega)*(1.+a) * Bnorm    155   G4double omega = -(q*s_omega)*(1.+a) * Bnorm;
196                                                   156 
197   G4double rotationangle = deltatime * omega;     157   G4double rotationangle = deltatime * omega;
198                                                   158 
199   G4Transform3D SpinRotation = G4Rotate3D(rota    159   G4Transform3D SpinRotation = G4Rotate3D(rotationangle,B.unit());
200                                                   160 
201   G4Vector3D Spin = aStep.GetTrack() -> GetPol    161   G4Vector3D Spin = aStep.GetTrack() -> GetPolarization();
202                                                   162 
203   G4Vector3D newSpin = SpinRotation * Spin;       163   G4Vector3D newSpin = SpinRotation * Spin;
204                                                   164 
205 #ifdef G4VERBOSE                                  165 #ifdef G4VERBOSE
206   if (GetVerboseLevel()>2) {                      166   if (GetVerboseLevel()>2) {
207     G4double normspin = std::sqrt(Spin*Spin);     167     G4double normspin = std::sqrt(Spin*Spin);
208     G4double normnewspin = std::sqrt(newSpin*n    168     G4double normnewspin = std::sqrt(newSpin*newSpin);
209     //G4double cosalpha = Spin*newSpin/normspi    169     //G4double cosalpha = Spin*newSpin/normspin/normnewspin;
210     //G4double alpha = std::acos(cosalpha);       170     //G4double alpha = std::acos(cosalpha);
211                                                   171 
212     G4cout << "AT REST::: PARAMETERS " << G4en    172     G4cout << "AT REST::: PARAMETERS " << G4endl;
213     G4cout << "Initial spin  : " << Spin  << G    173     G4cout << "Initial spin  : " << Spin  << G4endl;
214     G4cout << "Delta time    : " << deltatime     174     G4cout << "Delta time    : " << deltatime  << G4endl;
215     G4cout << "Rotation angle: " << rotationan    175     G4cout << "Rotation angle: " << rotationangle/rad  << G4endl;
216     G4cout << "New spin      : " << newSpin  <    176     G4cout << "New spin      : " << newSpin  << G4endl;
217     G4cout << "Checked norms : " << normspin <    177     G4cout << "Checked norms : " << normspin <<" " << normnewspin << G4endl;
218   }                                               178   }
219 #endif                                            179 #endif
220                                                   180 
221   return newSpin;                                 181   return newSpin;
222                                                   182 
223 }                                              << 
224                                                << 
225 void G4DecayWithSpin::ProcessDescription(std:: << 
226 {                                              << 
227   outFile << GetProcessName()                  << 
228     << ": Decay of particles considering paren << 
229     << "kinematics of daughters are dertermine << 
230 }                                                 183 }
231                                                   184