Geant4  10.01.p02
G4PolarizationHelper.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4PolarizationHelper.cc 68046 2013-03-13 14:31:38Z gcosmo $
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4PolarizationHelper
32 //
33 // Author: Andreas Schaelicke
34 //
35 // Creation date: 12.08.2006
36 //
37 // Modifications:
38 //
39 // Class Description:
40 //
41 // Provides some basic polarization transformation routines.
42 //
43 #include "G4PolarizationHelper.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4StokesVector.hh"
46 
47 
49 {
50  G4ThreeVector normal = (mom1.cross(mom2)).unit();
51  return normal;
52  // return 1./normal.mag()*normal;
53 }
54 
56 {
57  // compare also G4ThreeVector::rotateUz()
58 
59  if (uZ.x()==0. && uZ.y()==0.) {
60  return G4ThreeVector(0.,1.,0.);
61  }
62 
63  G4double invPerp = 1./std::sqrt(sqr(uZ.x())+sqr(uZ.y()));
64  return G4ThreeVector(-uZ.y()*invPerp,uZ.x()*invPerp,0);
65 }
66 
68 {
69  // compare also G4ThreeVector::rotateUz()
70 
71  if (uZ.x()==0. && uZ.y()==0.) {
72  if (uZ.z()>=0.) return G4ThreeVector(1.,0.,0.);
73  return G4ThreeVector(-1.,0.,0.);
74  }
75 
76  G4double perp = std::sqrt(sqr(uZ.x())+sqr(uZ.y()));
77  G4double invPerp = uZ.z()/perp;
78  return G4ThreeVector(uZ.x()*invPerp,uZ.y()*invPerp,-perp);
79 }
80 
82 {
83  G4double phi =2.*pi*G4UniformRand();
84  G4ThreeVector normal = std::cos(phi)*GetParticleFrameX(mom1)
85  + std::sin(phi)*G4PolarizationHelper::GetParticleFrameY(mom1);
86  return normal;
87 }
88 
89 
91 {
92  // compare also G4ThreeVector::rotateUz()
93 
94  if (uZ.x()==0. && uZ.y()==0.) {
95  if (uZ.z()>=0.) return spin;
96  return G4ThreeVector(-spin.x(),spin.y(),-spin.z());
97  }
98 
99  G4double perp = std::sqrt(sqr(uZ.x())+sqr(uZ.y()));
100  G4double invPerp = 1./perp;
101 
102  G4ThreeVector uX(uZ.x()*uZ.z()*invPerp,uZ.y()*uZ.z()*invPerp,-perp);
103  G4ThreeVector uY(-uZ.y()*invPerp,uZ.x()*invPerp,0);
104 
105  return G4ThreeVector(spin*uX,spin*uY,spin*uZ);
106 }
107 
109 {
110  G4double theta=0.;
111  G4cout<<"========================================\n\n";
112  for (G4int i=0; i<=10; ++i) {
113  theta=pi*i/10.;
114  G4ThreeVector zAxis = G4ThreeVector(std::sin(theta),0.,std::cos(theta));
115  if (i==5) zAxis = G4ThreeVector(1.,0.,0.);
116  if (i==10) zAxis = G4ThreeVector(0.,0.,-1.);
117  G4ThreeVector yAxis = GetParticleFrameY(zAxis);
118 
119  G4cout<<zAxis<<" "<<zAxis.mag()<<"\n";
120  G4cout<<yAxis<<" "<<yAxis.mag()<<"\n";
121  G4ThreeVector xAxis = yAxis.cross(zAxis);
122  G4cout<<xAxis<<" "<<xAxis.mag()<<"\n\n";
123  }
124 
125  G4cout<<"========================================\n\n";
126 
127  for (G4int i=0; i<=10; ++i) {
128  theta=pi*i/10.;
129  G4ThreeVector zAxis = G4ThreeVector(0.,std::sin(theta),std::cos(theta));
130  if (i==5) zAxis = G4ThreeVector(0.,1.,0.);
131  if (i==10) zAxis = G4ThreeVector(0.,0.,-1.);
132  G4ThreeVector yAxis = GetParticleFrameY(zAxis);
133 
134  G4cout<<zAxis<<" "<<zAxis.mag()<<"\n";
135  G4cout<<yAxis<<" "<<yAxis.mag()<<"\n";
136  G4ThreeVector xAxis = yAxis.cross(zAxis);
137  G4cout<<xAxis<<" "<<xAxis.mag()<<"\n\n";
138 
139  G4cout<<"spat : "<<xAxis*yAxis.cross(zAxis)<<"\n\n";
140  }
141  G4cout<<"========================================\n\n";
142 }
143 
145 {
146  // check transformation procedure for polarisation transfer
147  // calculation in scattering processes
148  // a) transfer target polarisation in beam particle reference frame (PRF)
149  // b) calc correct asymmetry w.r.t. scattering plane
150  // c) determine incomming polarisation in interaction frame (IF)
151  // d) transfer outgoing polarisation from IF to PRF
152  G4cout<<"========================================\n\n";
153 
154  G4double theta=0.;
155 
156  G4ThreeVector dir0=G4ThreeVector(0.,0.,1.);
157  G4ThreeVector dir2=G4ThreeVector(std::sin(theta),0.,std::cos(theta));
158 
159  G4StokesVector pol0=G4ThreeVector(0.,0.,1.);
160  G4StokesVector pol1=G4ThreeVector(0.,0.,1.);
161 
162  pol1.rotateUz(dir0);
163 
164  G4cout<<"========================================\n\n";
165 
166 
167 }
static void TestPolarizationTransformations()
static G4ThreeVector GetRandomFrame(const G4ThreeVector &)
CLHEP::Hep3Vector G4ThreeVector
const G4double pi
int G4int
Definition: G4Types.hh:78
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
static void TestInteractionFrame()
static G4ThreeVector GetSpinInPRF(const G4ThreeVector &uZ, const G4ThreeVector &spin)
static G4ThreeVector GetParticleFrameY(const G4ThreeVector &)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
static G4ThreeVector GetParticleFrameX(const G4ThreeVector &)