ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StokesVector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4StokesVector.cc
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 //
27 // GEANT4 Class file
28 //
29 //
30 // File name: G4StokesVector
31 //
32 // Author: Andreas Schaelicke
33 //
34 // Creation date: 01.05.2005
35 //
36 // Modifications:
37 //
38 // Class Description:
39 //
40 // Provides Stokesvector representation employed in polarized
41 // processes.
42 //
43 #include "G4StokesVector.hh"
44 #include "G4PolarizationHelper.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "Randomize.hh"
47 
55 
58 {
59 }
60 
63 {
64 }
65 
67 {
68  return *this==ZERO;
69 }
70 
71 void G4StokesVector::RotateAz(G4ThreeVector nInteractionFrame,
72  G4ThreeVector particleDirection)
73 {
74  G4ThreeVector yParticleFrame =
75  G4PolarizationHelper::GetParticleFrameY(particleDirection);
76 
77 
78  G4double cosphi=yParticleFrame*nInteractionFrame;
79  if (cosphi>(1.+1.e-8) || cosphi<(-1.-1.e-8)) {
80  G4cout<<" warning G4StokesVector::RotateAz cosphi>1 or cosphi<-1\n"
81  <<" cosphi="<<cosphi<<"\n"
82  <<" zAxis="<<particleDirection<<" ("<<particleDirection.mag()<<")\n"
83  <<" yAxis="<<yParticleFrame<<" ("<<yParticleFrame.mag()<<")\n"
84  <<" nAxis="<<nInteractionFrame<<" ("
85  <<nInteractionFrame.mag()<<")"<<G4endl;
86  }
87  if (cosphi>1.) cosphi=1.;
88  else if (cosphi<-1.) cosphi=-1.;
89 
90 // G4cout<<" cosphi="<<cosphi<<"\n"
91 // <<" zAxis="<<particleDirection<<" ("<<particleDirection.mag()<<")\n"
92 // <<" yAxis="<<yParticleFrame<<" ("<<yParticleFrame.mag()<<","<<(yParticleFrame*particleDirection)<<")\n"
93 // <<" nAxis="<<nInteractionFrame<<" ("
94 // <<nInteractionFrame.mag()<<")"<<G4endl;
95 
96  // G4double hel=sgn(cross(yParticleFrame*nInteractionFrame)*zInteractionFrame);
97  // Why not particleDirection instead of zInteractionFrame ???!!!
98  // -> is the same, since SYSIN is called with p1, and p2 as first parameter!
99  G4double hel=(yParticleFrame.cross(nInteractionFrame)*particleDirection)>0?1.:-1.;
100 
101  G4double sinphi=hel*std::sqrt(1.-cosphi*cosphi);
102  // G4cout<<" sin2 + cos2 -1 = "<<(sinphi*sinphi+cosphi*cosphi-1)<<"\n";
103 
104  RotateAz(cosphi,sinphi);
105 }
106 
107 
109  G4ThreeVector particleDirection)
110 {
111  // note if incomming particle is on z-axis,
112  // we might encounter some nummerical problems, since
113  // nInteratonFrame and yParticleFrame are actually (almost) the same momentum
114  // and the normalization is only good to 10^-12 !
115 
116  G4ThreeVector yParticleFrame =
117  G4PolarizationHelper::GetParticleFrameY(particleDirection);
118  G4double cosphi=yParticleFrame*nInteractionFrame;
119 
120  if (cosphi>1.+1.e-8 || cosphi<-1.-1.e-8) {
121  G4cout<<" warning G4StokesVector::RotateAz cosphi>1 or cosphi<-1\n";
122  }
123  if (cosphi>1.) cosphi=1.;
124  else if (cosphi<-1.)cosphi=-1.;
125 
126  // check sign once more!
127  G4double hel=(yParticleFrame.cross(nInteractionFrame)*particleDirection)>0?1.:-1.;
128  G4double sinphi=hel*std::sqrt(std::fabs(1.-cosphi*cosphi));
129  RotateAz(cosphi,-sinphi);
130 }
131 
133 {
134  if (!isPhoton) {
135  G4double xsi1= cosphi*p1() + sinphi*p2();
136  G4double xsi2= -sinphi*p1() + cosphi*p2();
137  setX(xsi1);
138  setY(xsi2);
139  return;
140  }
141 
142  G4double sin2phi=2.*cosphi*sinphi;
143  G4double cos2phi=cosphi*cosphi-sinphi*sinphi;
144 
145  G4double xsi1= cos2phi*p1() + sin2phi*p2();
146  G4double xsi2= -sin2phi*p1() + cos2phi*p2();
147  setX(xsi1);
148  setY(xsi2);
149 }
150 
152 {
153  G4double bet=getPhi();
154  if (isPhoton) { bet *= 0.5; }
155  return bet;
156 }
157 
159 {
160  G4double costheta=2.*G4UniformRand()-1.;
161  G4double sintheta=std::sqrt(1.-costheta*costheta);
162  G4double aphi =2.*pi*G4UniformRand();
163  setX(std::sin(aphi)*sintheta);
164  setY(std::cos(aphi)*sintheta);
165  setZ(costheta);
166 }
167 
169 {
170  if (G4UniformRand()>0.5) setX(1.);
171  else setX(-1.);
172  setY(0.);
173  setZ(0.);
174 }
175 
177 {
178  setX(0.);
179  if (G4UniformRand()>0.5) setY(1.);
180  else setY(-1.);
181  setZ(0.);
182 }
183 
185 {
186  setX(0.);
187  setY(0.);
188  if (G4UniformRand()>0.5) setZ(1.);
189  else setZ(-1.);
190 }
191 
193 {
194  setZ(-z());
195 }
196 
198 {
199  // delta x = sqrt[ ( <x^2> - <x>^2 )/(n-1) ]
200  G4StokesVector mean=(1./n)*(*this);
201  return G4StokesVector((1./(n-1.)*((1./n)*sum2 - mean.PolSqr()))).PolSqrt();
202 }
203 
205  {return G4ThreeVector(b.x()!=0. ? x()/b.x() : 11111.,
206  b.y()!=0. ? y()/b.y() : 11111.,
207  b.z()!=0. ? z()/b.z() : 11111.);}