Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/polarisation/src/G4PolarizedComptonModel.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/electromagnetic/polarisation/src/G4PolarizedComptonModel.cc (Version 11.3.0) and /processes/electromagnetic/polarisation/src/G4PolarizedComptonModel.cc (Version 10.1.p1)


  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 //
                                                   >>  27 // $Id: G4PolarizedComptonModel.cc 82755 2014-07-08 14:07:29Z gcosmo $
                                                   >>  28 //
 26 // -------------------------------------------     29 // -------------------------------------------------------------------
 27 //                                                 30 //
 28 // Geant4 Class file                           <<  31 // GEANT4 Class file
                                                   >>  32 //
 29 //                                                 33 //
 30 // File name:     G4PolarizedComptonModel          34 // File name:     G4PolarizedComptonModel
 31 //                                                 35 //
 32 // Author:        Andreas Schaelicke               36 // Author:        Andreas Schaelicke
                                                   >>  37 //
                                                   >>  38 // Creation date: 01.05.2005
                                                   >>  39 //
                                                   >>  40 // Modifications:
                                                   >>  41 // 18-07-06 use newly calculated cross sections (P. Starovoitov)
                                                   >>  42 // 21-08-05 update interface (A. Schaelicke)
                                                   >>  43 //
                                                   >>  44 // Class Description:
                                                   >>  45 //
                                                   >>  46 // -------------------------------------------------------------------
                                                   >>  47 //
                                                   >>  48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 33                                                    50 
 34 #include "G4PolarizedComptonModel.hh"              51 #include "G4PolarizedComptonModel.hh"
 35                                                << 
 36 #include "G4Exp.hh"                            << 
 37 #include "G4Log.hh"                            << 
 38 #include "G4ParticleChangeForGamma.hh"         << 
 39 #include "G4PhysicalConstants.hh"                  52 #include "G4PhysicalConstants.hh"
                                                   >>  53 #include "G4Electron.hh"
                                                   >>  54 #include "G4Gamma.hh"
                                                   >>  55 #include "Randomize.hh"
                                                   >>  56 #include "G4DataVector.hh"
                                                   >>  57 #include "G4ParticleChangeForGamma.hh"
                                                   >>  58 
                                                   >>  59 
                                                   >>  60 #include "G4StokesVector.hh"
 40 #include "G4PolarizationManager.hh"                61 #include "G4PolarizationManager.hh"
 41 #include "G4PolarizationHelper.hh"                 62 #include "G4PolarizationHelper.hh"
 42 #include "G4PolarizedComptonXS.hh"             <<  63 #include "G4PolarizedComptonCrossSection.hh"
 43 #include "G4StokesVector.hh"                   << 
 44 #include "G4SystemOfUnits.hh"                  << 
 45                                                    64 
 46 //....oooOO0OOooo........oooOO0OOooo........oo     65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  66 
 47 G4PolarizedComptonModel::G4PolarizedComptonMod     67 G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*,
 48                                                <<  68              const G4String& nam)
 49   : G4KleinNishinaCompton(nullptr, nam)        <<  69   : G4KleinNishinaCompton(0,nam),
 50   , fVerboseLevel(0)                           <<  70     verboseLevel(0)
 51 {                                                  71 {
 52   fCrossSectionCalculator = new G4PolarizedCom <<  72   crossSectionCalculator = new G4PolarizedComptonCrossSection();
 53   fBeamPolarization       = G4StokesVector::ZE << 
 54   fTargetPolarization     = G4StokesVector::ZE << 
 55 }                                                  73 }
 56                                                    74 
 57 //....oooOO0OOooo........oooOO0OOooo........oo     75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >>  76 
 58 G4PolarizedComptonModel::~G4PolarizedComptonMo     77 G4PolarizedComptonModel::~G4PolarizedComptonModel()
 59 {                                                  78 {
 60   delete fCrossSectionCalculator;              <<  79   delete crossSectionCalculator;
 61 }                                                  80 }
 62                                                    81 
 63 //....oooOO0OOooo........oooOO0OOooo........oo <<  82 G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom
 64 G4double G4PolarizedComptonModel::ComputeAsymm <<  83                        (G4double gammaEnergy, G4double /*Z*/)
 65                                                <<  84  
 66 {                                                  85 {
 67   G4double asymmetry = 0.0;                    <<  86   G4double asymmetry = 0.0 ;
 68                                                    87 
 69   G4double k0 = gammaEnergy / electron_mass_c2 <<  88   G4double k0 = gammaEnergy / electron_mass_c2 ;
 70   G4double k1 = 1. + 2. * k0;                  <<  89   G4double k1 = 1 + 2*k0 ;
 71                                                    90 
 72   asymmetry = -k0;                                 91   asymmetry = -k0;
 73   asymmetry *=                                 <<  92   asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
 74     (k0 + 1.) * sqr(k1) * G4Log(k1) - 2. * k0  <<  93   asymmetry /= ((k0 - 2.)*k0  -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);    
 75   asymmetry /= ((k0 - 2.) * k0 - 2.) * sqr(k1) <<  94 
 76                2. * k0 * (k0 * (k0 + 1.) * (k0 <<  95   // G4cout<<"energy = "<<GammaEnergy<<"  asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl;
 77                                                <<  96   if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl;
 78   if(asymmetry > 1.)                           << 
 79   {                                            << 
 80     G4ExceptionDescription ed;                 << 
 81     ed << "ERROR in G4PolarizedComptonModel::C << 
 82        << " asymmetry = " << asymmetry << "\n" << 
 83     G4Exception("G4PolarizedComptonModel::Comp << 
 84                 JustWarning, ed);              << 
 85   }                                            << 
 86                                                    97 
 87   return asymmetry;                                98   return asymmetry;
 88 }                                                  99 }
 89                                                   100 
 90 //....oooOO0OOooo........oooOO0OOooo........oo << 101 
 91 G4double G4PolarizedComptonModel::ComputeCross    102 G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom(
 92   const G4ParticleDefinition* pd, G4double kin << 103                                 const G4ParticleDefinition* pd,
 93   G4double cut, G4double emax)                 << 104                                       G4double kinEnergy, 
                                                   >> 105                                       G4double Z, 
                                                   >> 106                                       G4double A, 
                                                   >> 107                                       G4double cut,
                                                   >> 108                                       G4double emax)
 94 {                                                 109 {
 95   G4double xs = G4KleinNishinaCompton::Compute << 110   double xs = 
 96     pd, kinEnergy, Z, A, cut, emax);           << 111     G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy,
 97   G4double polzz = fBeamPolarization.p3() * fT << 112                   Z,A,cut,emax);
 98   if(polzz > 0.0)                              << 113   G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z();
 99   {                                            << 114   if (polzz > 0.0) {
100     G4double asym = ComputeAsymmetryPerAtom(ki << 115     G4double asym = ComputeAsymmetryPerAtom(kinEnergy, Z);  
101     xs *= (1. + polzz * asym);                 << 116     xs*=(1.+polzz*asym);
102   }                                               117   }
103   return xs;                                      118   return xs;
104 }                                                 119 }
105                                                   120 
106 //....oooOO0OOooo........oooOO0OOooo........oo    121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 122 
107 void G4PolarizedComptonModel::SampleSecondarie    123 void G4PolarizedComptonModel::SampleSecondaries(
108   std::vector<G4DynamicParticle*>* fvect, cons << 124                               std::vector<G4DynamicParticle*>* fvect,
109   const G4DynamicParticle* aDynamicGamma, G4do << 125                               const G4MaterialCutsCouple*,
                                                   >> 126             const G4DynamicParticle* aDynamicGamma,
                                                   >> 127             G4double, G4double)
110 {                                                 128 {
111   // do nothing below the threshold               129   // do nothing below the threshold
112   if(aDynamicGamma->GetKineticEnergy() <= LowE << 130   if(aDynamicGamma->GetKineticEnergy() <= LowEnergyLimit()) { return; }
113   {                                            << 
114     return;                                    << 
115   }                                            << 
116                                                   131 
117   const G4Track* aTrack       = fParticleChang << 132   const G4Track * aTrack = fParticleChange->GetCurrentTrack();
118   G4VPhysicalVolume* aPVolume = aTrack->GetVol << 133   G4VPhysicalVolume*  aPVolume  = aTrack->GetVolume();
119   G4LogicalVolume* aLVolume   = aPVolume->GetL << 134   G4LogicalVolume*    aLVolume  = aPVolume->GetLogicalVolume();
120                                                << 135 
121   if(fVerboseLevel >= 1)                       << 136   if (verboseLevel >= 1) {
122   {                                            << 137     G4cout<<"G4PolarizedComptonModel::SampleSecondaries in "
123     G4cout << "G4PolarizedComptonModel::Sample << 138           <<  aLVolume->GetName() <<G4endl;
124            << aLVolume->GetName() << G4endl;   << 
125   }                                               139   }
126   G4PolarizationManager* polarizationManager = << 140   G4PolarizationManager * polarizationManager = 
127     G4PolarizationManager::GetInstance();         141     G4PolarizationManager::GetInstance();
128                                                   142 
129   // obtain polarization of the beam              143   // obtain polarization of the beam
130   fBeamPolarization = G4StokesVector(aDynamicG << 144   theBeamPolarization =  aDynamicGamma->GetPolarization();
131   fBeamPolarization.SetPhoton();               << 145   theBeamPolarization.SetPhoton();
132                                                   146 
133   // obtain polarization of the media             147   // obtain polarization of the media
134   G4bool targetIsPolarized = polarizationManag    148   G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
135   fTargetPolarization = polarizationManager->G << 149   theTargetPolarization = 
                                                   >> 150     polarizationManager->GetVolumePolarization(aLVolume);
136                                                   151 
137   // if beam is linear polarized or target is  << 152   // if beam is linear polarized or target is transversely polarized 
138   // determine the angle to x-axis                153   // determine the angle to x-axis
139   // (assumes same PRF as in the polarization     154   // (assumes same PRF as in the polarization definition)
                                                   >> 155 
140   G4ThreeVector gamDirection0 = aDynamicGamma-    156   G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
141                                                   157 
142   // transfer fTargetPolarization              << 158   // transfere theTargetPolarization 
143   // into the gamma frame (problem electron is    159   // into the gamma frame (problem electron is at rest)
144   if(targetIsPolarized)                        << 160   if (targetIsPolarized) {
145   {                                            << 161     theTargetPolarization.rotateUz(gamDirection0);
146     fTargetPolarization.rotateUz(gamDirection0 << 
147   }                                               162   }
148   // The scattered gamma energy is sampled acc << 163   // The scattered gamma energy is sampled according to 
149   // Klein - Nishina formula.                     164   // Klein - Nishina formula.
150   // The random number techniques of Butcher & << 165   // The random number techniques of Butcher & Messel are used 
151   // (Nuc Phys 20(1960),15).                      166   // (Nuc Phys 20(1960),15).
152   // Note : Effects due to binding of atomic e << 167   // Note : Effects due to binding of atomic electrons are negliged.
153                                                << 168  
154   G4double gamEnergy0 = aDynamicGamma->GetKine    169   G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
155   G4double E0_m       = gamEnergy0 / electron_ << 170   G4double E0_m = gamEnergy0 / electron_mass_c2 ;
156                                                   171 
157   // sample the energy rate of the scattered g << 172   //
158   G4double epsilon, sint2;                     << 173   // sample the energy rate of the scattered gamma 
159   G4double onecost = 0.0;                      << 174   //
160   G4double Phi     = 0.0;                      << 175 
161   G4double greject = 1.0;                      << 176   G4double epsilon, epsilonsq, onecost, sint2, greject ;
162   G4double cosTeta = 1.0;                      << 177 
163   G4double sinTeta = 0.0;                      << 178   G4double eps0       = 1./(1. + 2.*E0_m);
164                                                << 179   G4double epsilon0sq = eps0*eps0;
165   G4double eps0       = 1. / (1. + 2. * E0_m); << 180   G4double alpha1     = - std::log(eps0);
166   G4double epsilon0sq = eps0 * eps0;           << 181   G4double alpha2     = 0.5*(1.- epsilon0sq);
167   G4double alpha1     = -G4Log(eps0);          << 182 
168   G4double alpha2     = alpha1 + 0.5 * (1. - e << 183   G4double polarization = 
169                                                << 184     theBeamPolarization.p3()*theTargetPolarization.p3();
170   G4double polarization = fBeamPolarization.p3 << 185 
171                                                << 186   G4int nloop = 0;
172   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 187   do {
173   G4int nloop                           = 0;   << 188     ++nloop;
174   G4bool end                            = fals << 189     // false interaction if too many iterations
175                                                << 190     if(nloop > 1000) { return; }
176   G4double rndm[3];                            << 191 
177                                                << 192     if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
178   do                                           << 193       epsilon   = std::exp(-alpha1*G4UniformRand());   // epsilon0**r
179   {                                            << 194       epsilonsq = epsilon*epsilon; 
180     do                                         << 195 
181     {                                          << 196     } else {
182       ++nloop;                                 << 197       epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
183       // false interaction if too many iterati << 198       epsilon   = std::sqrt(epsilonsq);
184       if(nloop > fLoopLim)                     << 199     }
185       {                                        << 200 
186         PrintWarning(aDynamicGamma, nloop, gre << 201     onecost = (1.- epsilon)/(epsilon*E0_m);
187                      "too many iterations");   << 202     sint2   = onecost*(2.-onecost);
188         return;                                << 203 
189       }                                        << 204     G4double gdiced = 2.*(1./epsilon+epsilon);
190                                                << 205     G4double gdist  = 1./epsilon + epsilon - sint2 
191       // 3 random numbers to sample scattering << 206       - polarization*(1./epsilon-epsilon)*(1.-onecost);
192       rndmEngineMod->flatArray(3, rndm);       << 207 
193                                                << 208     greject = gdist/gdiced;
194       if(alpha1 > alpha2 * rndm[0])            << 209 
195       {                                        << 210     if (greject>1) {
196         epsilon = G4Exp(-alpha1 * rndm[1]);    << 211       G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
197       }                                        << 212       <<" costh rejection does not work properly: "<<greject
198       else                                     << 213       <<G4endl;
199       {                                        << 214     }
200         epsilon = std::sqrt(epsilon0sq + (1. - << 215   } while (greject < G4UniformRand());
201       }                                        << 216  
202                                                << 217   //
203       onecost = (1. - epsilon) / (epsilon * E0 << 218   // scattered gamma angles. ( Z - axis along the parent gamma)
204       sint2   = onecost * (2. - onecost);      << 219   //
205                                                << 220 
206       G4double gdiced = 2. * (1. / epsilon + e << 221   G4double cosTeta = 1. - onecost; 
207       G4double gdist  = 1. / epsilon + epsilon << 222   G4double sinTeta = std::sqrt (sint2);
208                        polarization * (1. / ep << 223   G4double Phi;
209                                                << 224   do {
210       greject = gdist / gdiced;                << 225     ++nloop;
211                                                << 226     // false interaction if too many iterations
212       if(greject > 1.0)                        << 227     if(nloop > 1000) { return; }
213       {                                        << 228 
214         PrintWarning(aDynamicGamma, nloop, gre << 229     Phi = twopi * G4UniformRand();
215                      "theta majoranta wrong"); << 230     G4double gdiced = 1./epsilon + epsilon - sint2 
216       }                                        << 231        + std::abs(theBeamPolarization.p3())*
217       // Loop checking, 03-Aug-2015, Vladimir  << 232        ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3())
218     } while(greject < rndm[2]);                << 233   +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1()) 
219                                                << 234             + sqr(theTargetPolarization.p2()))))
220     // assuming phi loop successful            << 235        +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + 
221     end = true;                                << 236        sqr(theBeamPolarization.p2())));
222                                                << 237 
223     // scattered gamma angles. ( Z - axis alon << 238     G4double gdist = 1./epsilon + epsilon - sint2 
224     cosTeta = 1. - onecost;                    << 239        + theBeamPolarization.p3()*
225     sinTeta = std::sqrt(sint2);                << 240        ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()
226     do                                         << 241   +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+
227     {                                          << 242              std::sin(Phi)*theTargetPolarization.p2()))
228       ++nloop;                                 << 243        -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1()
229                                                << 244          +std::sin(2.*Phi)*theBeamPolarization.p2());
230       // 2 random numbers to sample scattering << 245     greject = gdist/gdiced;
231       rndmEngineMod->flatArray(2, rndm);       << 246 
232                                                << 247     if (greject>1.+1.e-10 || greject<0) {
233       // false interaction if too many iterati << 248       G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
234       Phi = twopi * rndm[0];                   << 249       <<" phi rejection does not work properly: "<<greject<<G4endl;
235       if(nloop > fLoopLim)                     << 250     }
236       {                                        << 251     if (greject<1.e-3) {
237         PrintWarning(aDynamicGamma, nloop, gre << 252       G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n"
238                      "too many iterations");   << 253       <<" phi rejection does not work properly: "<<greject<<"\n";
239         return;                                << 254       G4cout<<" greject="<<greject<<"  phi="<<Phi<<"   cost="<<cosTeta<<"\n";
240       }                                        << 255       G4cout<<" gdiced="<<gdiced<<"   gdist="<<gdist<<"\n";
241                                                << 256       G4cout<<" eps="<<epsilon<<"    1/eps="<<1./epsilon<<"\n";
242       G4double gdiced = 1. / epsilon + epsilon << 257     }
243                         std::abs(fBeamPolariza << 258      
244                           (std::abs((1. / epsi << 259   } while (greject < G4UniformRand());
245                                     fTargetPol << 260   G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), 
246                            (1. - epsilon) * si << 261     dirz = cosTeta;
247                              (std::sqrt(sqr(fT << 
248                                         sqr(fT << 
249                         sint2 * (std::sqrt(sqr << 
250                                            sqr << 
251                                                << 
252       G4double gdist =                         << 
253         1. / epsilon + epsilon - sint2 +       << 
254         fBeamPolarization.p3() *               << 
255           ((1. / epsilon - epsilon) * cosTeta  << 
256            (1. - epsilon) * sinTeta *          << 
257              (std::cos(Phi) * fTargetPolarizat << 
258               std::sin(Phi) * fTargetPolarizat << 
259         sint2 * (std::cos(2. * Phi) * fBeamPol << 
260                  std::sin(2. * Phi) * fBeamPol << 
261       greject = gdist / gdiced;                << 
262                                                << 
263       if(greject > 1.0)                        << 
264       {                                        << 
265         PrintWarning(aDynamicGamma, nloop, gre << 
266                      "phi majoranta wrong");   << 
267       }                                        << 
268                                                << 
269       if(greject < 1.e-3)                      << 
270       {                                        << 
271         PrintWarning(aDynamicGamma, nloop, gre << 
272                      "phi loop ineffective");  << 
273         // restart theta loop                  << 
274         end = false;                           << 
275         break;                                 << 
276       }                                        << 
277                                                << 
278       // Loop checking, 03-Aug-2015, Vladimir  << 
279     } while(greject < rndm[1]);                << 
280   } while(!end);                               << 
281   G4double dirx = sinTeta * std::cos(Phi);     << 
282   G4double diry = sinTeta * std::sin(Phi);     << 
283   G4double dirz = cosTeta;                     << 
284                                                   262 
                                                   >> 263   //
285   // update G4VParticleChange for the scattere    264   // update G4VParticleChange for the scattered gamma
286   G4ThreeVector gamDirection1(dirx, diry, dirz << 265   //
                                                   >> 266    
                                                   >> 267   G4ThreeVector gamDirection1 ( dirx,diry,dirz );
287   gamDirection1.rotateUz(gamDirection0);          268   gamDirection1.rotateUz(gamDirection0);
288   G4double gamEnergy1 = epsilon * gamEnergy0;  << 269   G4double gamEnergy1 = epsilon*gamEnergy0;
289                                                   270 
290   G4double edep = 0.0;                            271   G4double edep = 0.0;
291   if(gamEnergy1 > lowestSecondaryEnergy)       << 272   if(gamEnergy1 > lowestSecondaryEnergy) {
292   {                                            << 
293     fParticleChange->ProposeMomentumDirection(    273     fParticleChange->ProposeMomentumDirection(gamDirection1);
294     fParticleChange->SetProposedKineticEnergy(    274     fParticleChange->SetProposedKineticEnergy(gamEnergy1);
295   }                                            << 275   } else { 
296   else                                         << 
297   {                                            << 
298     fParticleChange->ProposeTrackStatus(fStopA    276     fParticleChange->ProposeTrackStatus(fStopAndKill);
299     fParticleChange->SetProposedKineticEnergy(    277     fParticleChange->SetProposedKineticEnergy(0.0);
300     edep = gamEnergy1;                            278     edep = gamEnergy1;
301   }                                               279   }
                                                   >> 280  
                                                   >> 281   // 
                                                   >> 282   // calculate Stokesvector of final state photon and electron
                                                   >> 283   //
                                                   >> 284   G4ThreeVector  nInteractionFrame = 
                                                   >> 285     G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0);
302                                                   286 
303   // calculate Stokes vector of final state ph << 287   // transfere theBeamPolarization and theTargetPolarization 
304   G4ThreeVector nInteractionFrame =            << 
305     G4PolarizationHelper::GetFrame(gamDirectio << 
306                                                << 
307   // transfer fBeamPolarization and fTargetPol << 
308   // into the interaction frame (note electron    288   // into the interaction frame (note electron is in gamma frame)
309   if(fVerboseLevel >= 1)                       << 289   if (verboseLevel>=1) {
310   {                                            << 290     G4cout << "========================================\n";
311     G4cout << "=============================== << 291     G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n";
312     G4cout << " nInteractionFrame = " << nInte << 292     G4cout << " GammaDirection0 = " <<gamDirection0<<"\n";
313     G4cout << " GammaDirection0 = " << gamDire << 293     G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
314     G4cout << " gammaPolarization = " << fBeam << 294     G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
315     G4cout << " electronPolarization = " << fT << 295   }
316   }                                            << 296 
317                                                << 297   theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
318   fBeamPolarization.InvRotateAz(nInteractionFr << 298   theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0);
319   fTargetPolarization.InvRotateAz(nInteraction << 299 
320                                                << 300   if (verboseLevel>=1) {
321   if(fVerboseLevel >= 1)                       << 301     G4cout << "----------------------------------------\n";
322   {                                            << 302     G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n";
323     G4cout << "------------------------------- << 303     G4cout << " electronPolarization = " <<theTargetPolarization<<"\n";
324     G4cout << " gammaPolarization = " << fBeam << 304     G4cout << "----------------------------------------\n";
325     G4cout << " electronPolarization = " << fT << 
326     G4cout << "------------------------------- << 
327   }                                               305   }
328                                                   306 
329   // initialize the polarization transfer matr    307   // initialize the polarization transfer matrix
330   fCrossSectionCalculator->Initialize(epsilon, << 308   crossSectionCalculator->Initialize(epsilon,E0_m,0.,
331                                       fTargetP << 309              theBeamPolarization,
332                                                << 310              theTargetPolarization,2);
333   if(gamEnergy1 > lowestSecondaryEnergy)       << 311   
334   {                                            << 312   if(gamEnergy1 > lowestSecondaryEnergy) {
                                                   >> 313  
335     // in interaction frame                       314     // in interaction frame
336     // calculate polarization transfer to the     315     // calculate polarization transfer to the photon (in interaction plane)
337     fFinalGammaPolarization = fCrossSectionCal << 316     finalGammaPolarization = crossSectionCalculator->GetPol2();
338     if(fVerboseLevel >= 1)                     << 317     if (verboseLevel>=1) {
339     {                                          << 318       G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
340       G4cout << " gammaPolarization1 = " << fF << 
341     }                                             319     }
342     fFinalGammaPolarization.SetPhoton();       << 320     finalGammaPolarization.SetPhoton();
343                                                   321 
344     // translate polarization into particle re    322     // translate polarization into particle reference frame
345     fFinalGammaPolarization.RotateAz(nInteract << 323     finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1);
346     if(fFinalGammaPolarization.mag() > 1. + 1. << 324     //store polarization vector
347     {                                          << 325     fParticleChange->ProposePolarization(finalGammaPolarization);
348       G4ExceptionDescription ed;               << 326     if (finalGammaPolarization.mag() > 1.+1.e-8){
349       ed << "ERROR in Polarizaed Compton Scatt << 327       G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
350       ed << "Polarization of final photon more << 328       G4cout<<"Polarization of final photon more than 100%"<<G4endl;
351       ed << fFinalGammaPolarization            << 329       G4cout<<finalGammaPolarization<<" mag = "
352          << " mag = " << fFinalGammaPolarizati << 330       <<finalGammaPolarization.mag()<<G4endl;
353       G4Exception("G4PolarizedComptonModel::Sa << 
354                   FatalException, ed);         << 
355     }                                             331     }
356     // store polarization vector               << 332     if (verboseLevel>=1) {
357     fParticleChange->ProposePolarization(fFina << 333       G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n";
358     if(fVerboseLevel >= 1)                     << 334       G4cout << " GammaDirection1 = " <<gamDirection1<<"\n";
359     {                                          << 
360       G4cout << " gammaPolarization1 = " << fF << 
361       G4cout << " GammaDirection1 = " << gamDi << 
362     }                                             335     }
363   }                                               336   }
364                                                   337 
                                                   >> 338   //
365   // kinematic of the scattered electron          339   // kinematic of the scattered electron
                                                   >> 340   //
366   G4double eKinEnergy = gamEnergy0 - gamEnergy    341   G4double eKinEnergy = gamEnergy0 - gamEnergy1;
367                                                   342 
368   if(eKinEnergy > lowestSecondaryEnergy)       << 343   if (eKinEnergy > lowestSecondaryEnergy) {
369   {                                            << 344   
370     G4ThreeVector eDirection =                 << 345     G4ThreeVector eDirection = 
371       gamEnergy0 * gamDirection0 - gamEnergy1  << 346       gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
372     eDirection = eDirection.unit();               347     eDirection = eDirection.unit();
373                                                   348 
374     finalElectronPolarization = fCrossSectionC << 349     finalElectronPolarization = crossSectionCalculator->GetPol3();
375     if(fVerboseLevel >= 1)                     << 350     if (verboseLevel>=1) {
376     {                                          << 351       G4cout << " electronPolarization1 = " 
377       G4cout << " electronPolarization1 = " << << 352        <<finalElectronPolarization<<"\n";
378              << G4endl;                        << 
379     }                                             353     }
380     // transfer into particle reference frame     354     // transfer into particle reference frame
381     finalElectronPolarization.RotateAz(nIntera << 355     finalElectronPolarization.RotateAz(nInteractionFrame,eDirection);
382     if(fVerboseLevel >= 1)                     << 356     if (verboseLevel>=1) {
383     {                                          << 357       G4cout << " electronPolarization1 = " 
384       G4cout << " electronPolarization1 = " << << 358        <<finalElectronPolarization<<"\n";
385              << G4endl << " ElecDirection = "  << 359       G4cout << " ElecDirection = " <<eDirection<<"\n";
386     }                                             360     }
387                                                   361 
388     // create G4DynamicParticle object for the    362     // create G4DynamicParticle object for the electron.
389     G4DynamicParticle* aElectron =             << 363     G4DynamicParticle* aElectron = 
390       new G4DynamicParticle(theElectron, eDire << 364       new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
391     // store polarization vector               << 365     //store polarization vector
392     if(finalElectronPolarization.mag() > 1. +  << 366     if (finalElectronPolarization.mag() > 1.+1.e-8){
393     {                                          << 367       G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl;
394       G4ExceptionDescription ed;               << 368       G4cout<<"Polarization of final electron more than 100%"<<G4endl;
395       ed << "ERROR in Polarized Compton Scatte << 369       G4cout<<finalElectronPolarization<<" mag = "
396       ed << "Polarization of final electron mo << 370       <<finalElectronPolarization.mag()<<G4endl;
397       ed << finalElectronPolarization          << 
398          << " mag = " << finalElectronPolariza << 
399       G4Exception("G4PolarizedComptonModel::Sa << 
400                   FatalException, ed);         << 
401     }                                             371     }
402     aElectron->SetPolarization(finalElectronPo    372     aElectron->SetPolarization(finalElectronPolarization.p1(),
403                                finalElectronPo << 373              finalElectronPolarization.p2(),
404                                finalElectronPo << 374              finalElectronPolarization.p3());
405     fvect->push_back(aElectron);                  375     fvect->push_back(aElectron);
406   }                                            << 376   } else {
407   else                                         << 377     edep += eKinEnergy;  
408   {                                            << 
409     edep += eKinEnergy;                        << 
410   }                                               378   }
411   // energy balance                               379   // energy balance
412   if(edep > 0.0)                               << 380   if(edep > 0.0) { 
413   {                                            << 
414     fParticleChange->ProposeLocalEnergyDeposit    381     fParticleChange->ProposeLocalEnergyDeposit(edep);
415   }                                               382   }
416 }                                                 383 }
417                                                   384 
418 //....oooOO0OOooo........oooOO0OOooo........oo    385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
419 void G4PolarizedComptonModel::PrintWarning(con << 386 
420                                            G4i << 387 
421                                            G4d << 
422                                            con << 
423 {                                              << 
424   G4ExceptionDescription ed;                   << 
425   ed << "Problem of scattering sampling: " <<  << 
426      << "Niter= " << nloop << " grej= " << gre << 
427      << " cos(theta)= " << 1.0 - onecos << " p << 
428      << "Gamma E(MeV)= " << dp->GetKineticEner << 
429      << " dir= " << dp->GetMomentumDirection() << 
430      << " pol= " << dp->GetPolarization();     << 
431   G4Exception("G4PolarizedComptonModel::Sample << 
432               JustWarning, ed, "");            << 
433 }                                              << 
434                                                   388