Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/xrays/src/G4ForwardXrayTR.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/xrays/src/G4ForwardXrayTR.cc (Version 11.3.0) and /processes/electromagnetic/xrays/src/G4ForwardXrayTR.cc (Version 4.0)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 //
                                                   >>  24 // $Id: G4ForwardXrayTR.cc,v 1.7 2001/10/24 16:40:54 maire Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-04-00 $
                                                   >>  26 //
                                                   >>  27 // G4ForwardXrayTR class -- implementation file
                                                   >>  28 
                                                   >>  29 // GEANT 4 class implementation file --- Copyright CERN 1995
                                                   >>  30 // CERN Geneva Switzerland
                                                   >>  31 
                                                   >>  32 // For information related to this code, please, contact
                                                   >>  33 // CERN, CN Division, ASD Group
 26 // History:                                        34 // History:
 27 // 1st version 11.09.97 V. Grichine (Vladimir.     35 // 1st version 11.09.97 V. Grichine (Vladimir.Grichine@cern.ch )
 28 // 2nd version 17.12.97 V. Grichine                36 // 2nd version 17.12.97 V. Grichine
 29 // 17-09-01, migration of Materials to pure ST <<  37 // 17-09-01, migration of Materials to pure STL (mma) 
 30 // 10-03-03, migration to "cut per region" (V. << 
 31 // 03.06.03, V.Ivanchenko fix compilation warn << 
 32                                                    38 
 33 #include "G4ForwardXrayTR.hh"                  << 
 34                                                    39 
 35 #include "globals.hh"                          <<  40 #include <math.h>
 36 #include "G4Gamma.hh"                          <<  41 
 37 #include "G4GeometryTolerance.hh"              <<  42 // #include "G4ios.hh"
                                                   >>  43 // #include <fstream.h>
                                                   >>  44 // #include <stdlib.h>
                                                   >>  45 
                                                   >>  46 #include "G4ForwardXrayTR.hh"
 38 #include "G4Material.hh"                           47 #include "G4Material.hh"
 39 #include "G4PhysicalConstants.hh"              <<  48 #include "globals.hh"
 40 #include "G4PhysicsLinearVector.hh"            << 
 41 #include "G4PhysicsLogVector.hh"               << 
 42 #include "G4PhysicsTable.hh"                       49 #include "G4PhysicsTable.hh"
 43 #include "G4PhysicsVector.hh"                      50 #include "G4PhysicsVector.hh"
 44 #include "G4Poisson.hh"                        <<  51 #include "G4PhysicsLinearVector.hh"
 45 #include "G4ProductionCutsTable.hh"            <<  52 #include "G4PhysicsLogVector.hh"
 46 #include "G4SystemOfUnits.hh"                  <<  53 
 47 #include "G4PhysicsModelCatalog.hh"            <<  54 // Table initialization
                                                   >>  55 
                                                   >>  56 // G4PhysicsTable* G4ForwardXrayTR::fAngleDistrTable  = NULL ;
                                                   >>  57 // G4PhysicsTable* G4ForwardXrayTR::fEnergyDistrTable = NULL ;
                                                   >>  58 
                                                   >>  59 
                                                   >>  60 // Initialization of local constants
                                                   >>  61 
                                                   >>  62 G4int    G4ForwardXrayTR::fSympsonNumber =  100       ;
                                                   >>  63 
                                                   >>  64 G4double G4ForwardXrayTR::fTheMinEnergyTR   =    1.0*keV  ;
                                                   >>  65 G4double G4ForwardXrayTR::fTheMaxEnergyTR   =  100.0*keV  ;
                                                   >>  66 G4double G4ForwardXrayTR::fTheMaxAngle    =      1.0e-3   ;
                                                   >>  67 G4double G4ForwardXrayTR::fTheMinAngle    =      5.0e-6   ;
                                                   >>  68 G4int    G4ForwardXrayTR::fBinTR            =   50        ;
                                                   >>  69 
                                                   >>  70 G4double G4ForwardXrayTR::fMinProtonTkin = 100.0*GeV  ;
                                                   >>  71 G4double G4ForwardXrayTR::fMaxProtonTkin = 100.0*TeV  ;
                                                   >>  72 G4int    G4ForwardXrayTR::fTotBin        =  50        ;
                                                   >>  73 // Proton energy vector initialization
                                                   >>  74 
                                                   >>  75 G4PhysicsLogVector* G4ForwardXrayTR::
                                                   >>  76 fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin,
                                                   >>  77                                              fMaxProtonTkin,
                                                   >>  78                                                     fTotBin  ) ;
                                                   >>  79 
                                                   >>  80 G4double G4ForwardXrayTR::fPlasmaCof = 4.0*pi*fine_structure_const*
                                                   >>  81                                        hbarc*hbarc*hbarc/electron_mass_c2 ;
                                                   >>  82 
                                                   >>  83 G4double G4ForwardXrayTR::fCofTR     = fine_structure_const/pi ;
                                                   >>  84 
                                                   >>  85 
                                                   >>  86 /*   ************************************************************************
                                                   >>  87 
                                                   >>  88 
                                                   >>  89 ///////////////////////////////////////////////////////////////////////
                                                   >>  90 //
                                                   >>  91 // Constructor for preparation tables with angle and energy TR distributions
                                                   >>  92 // in all materials involved in test program. Lorentz factors correspond to
                                                   >>  93 // kinetic energies of protons between 100*GeV and 100*TeV, ~ 10^2-10^5
                                                   >>  94 //
                                                   >>  95 // Recommended only for use in applications with 
                                                   >>  96 // few light materials involved                     !!!!!!!!!!!!!!
                                                   >>  97 
                                                   >>  98 G4ForwardXrayTR::G4ForwardXrayTR()
                                                   >>  99   : G4TransitionRadiation("XrayTR")
                                                   >> 100 {
                                                   >> 101   G4int iMat, jMat, iTkin, iTR, iPlace ;
                                                   >> 102   static 
                                                   >> 103   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
                                                   >> 104   G4int numOfMat = G4Material::GetNumberOfMaterials();
                                                   >> 105   fGammaCutInKineticEnergy = new G4double[numOfMat] ;
                                                   >> 106   fGammaCutInKineticEnergy = fPtrGamma->GetEnergyCuts() ;
                                                   >> 107   fMatIndex1 = -1 ;
                                                   >> 108   fMatIndex2 = -1 ;
                                                   >> 109   fAngleDistrTable  = new G4PhysicsTable(numOfMat*(numOfMat - 1)*fTotBin) ;
                                                   >> 110   fEnergyDistrTable = new G4PhysicsTable(numOfMat*(numOfMat - 1)*fTotBin) ;
                                                   >> 111 
                                                   >> 112         G4PhysicsLogVector* aVector = new G4PhysicsLogVector(fMinProtonTkin,
                                                   >> 113                                                              fMaxProtonTkin,
                                                   >> 114                                                                     fTotBin  ) ;
                                                   >> 115 
                                                   >> 116   for(iMat=0;iMat<numOfMat;iMat++) // loop over pairs of different materials
                                                   >> 117   {
                                                   >> 118     for(jMat=0;jMat<numOfMat;jMat++)  // transition iMat -> jMat !!!
                                                   >> 119     {
                                                   >> 120       if(iMat == jMat)   continue ;      // no TR !!
                                                   >> 121       else
                                                   >> 122       {
                                                   >> 123         const G4Material* mat1 = (*theMaterialTable)[iMat] ;
                                                   >> 124         const G4Material* mat2 = (*theMaterialTable)[jMat] ;
                                                   >> 125 
                                                   >> 126         fSigma1 = fPlasmaCof*(mat1->GetElectronDensity()) ;
                                                   >> 127         fSigma2 = fPlasmaCof*(mat2->GetElectronDensity()) ;
                                                   >> 128 
                                                   >> 129 //        fGammaTkinCut = fGammaCutInKineticEnergy[jMat] ; // TR photon in jMat !
                                                   >> 130           fGammaTkinCut = 0.0 ;
                                                   >> 131 
                                                   >> 132         if(fGammaTkinCut > fTheMinEnergyTR) // setting of min/max TR energies 
                                                   >> 133   {
                                                   >> 134           fMinEnergyTR = fGammaTkinCut ;
                                                   >> 135   }
                                                   >> 136         else
                                                   >> 137   {
                                                   >> 138           fMinEnergyTR = fTheMinEnergyTR ;
                                                   >> 139   }
                                                   >> 140         if(fGammaTkinCut > fTheMaxEnergyTR)
                                                   >> 141   {
                                                   >> 142           fMaxEnergyTR = 2.0*fGammaTkinCut ;    // usually very low TR rate 
                                                   >> 143   }
                                                   >> 144         else
                                                   >> 145   {
                                                   >> 146           fMaxEnergyTR = fTheMaxEnergyTR ;
                                                   >> 147   }
                                                   >> 148         for(iTkin=0;iTkin<fTotBin;iTkin++)      // Lorentz factor loop
                                                   >> 149   {
                                                   >> 150           G4PhysicsLogVector* 
                                                   >> 151                     energyVector = new G4PhysicsLogVector(fMinEnergyTR,
                                                   >> 152                                                              fMaxEnergyTR,
                                                   >> 153                                                                    fBinTR  ) ;
                                                   >> 154           G4PhysicsLinearVector* 
                                                   >> 155                      angleVector = new G4PhysicsLinearVector(        0.0,
                                                   >> 156                                                              fMaxThetaTR,
                                                   >> 157                                                                   fBinTR  ) ;
                                                   >> 158           G4double energySum = 0.0 ;
                                                   >> 159           G4double angleSum  = 0.0 ;
                                                   >> 160           fGamma = 1.0 +   (aVector->GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
                                                   >> 161           fMaxThetaTR = 10000.0/(fGamma*fGamma) ;
                                                   >> 162           if(fMaxThetaTR > fTheMaxAngle)
                                                   >> 163           { 
                                                   >> 164             fMaxThetaTR = fTheMaxAngle ;
                                                   >> 165     }
                                                   >> 166           else
                                                   >> 167     {
                                                   >> 168             if(fMaxThetaTR < fTheMinAngle)
                                                   >> 169       {
                                                   >> 170               fMaxThetaTR = fTheMinAngle ;
                                                   >> 171       }
                                                   >> 172     }
                                                   >> 173           energyVector->PutValue(fBinTR-1,energySum) ;
                                                   >> 174           angleVector->PutValue(fBinTR-1,angleSum)   ;
                                                   >> 175 
                                                   >> 176           for(iTR=fBinTR-2;iTR>=0;iTR--)
                                                   >> 177     {
                                                   >> 178             energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR),
                                                   >> 179                                         energyVector->GetLowEdgeEnergy(iTR+1)) ; 
                                                   >> 180 
                                                   >> 181             angleSum  += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR),
                                                   >> 182                                          angleVector->GetLowEdgeEnergy(iTR+1)) ;
                                                   >> 183             energyVector->PutValue(iTR,energySum) ;
                                                   >> 184             angleVector->PutValue(iTR,angleSum)   ;
                                                   >> 185     }
                                                   >> 186           if(jMat < iMat)
                                                   >> 187     {
                                                   >> 188             iPlace = (iMat*(numOfMat-1)+jMat)*fTotBin+iTkin ;
                                                   >> 189     }
                                                   >> 190           else   // jMat > iMat right part of matrices (jMat-1) ! 
                                                   >> 191     {
                                                   >> 192             iPlace = (iMat*(numOfMat-1)+jMat-1)*fTotBin+iTkin ;
                                                   >> 193     } 
                                                   >> 194           fEnergyDistrTable->insertAt(iPlace,energyVector) ;
                                                   >> 195           fAngleDistrTable->insertAt(iPlace,angleVector) ;
                                                   >> 196   }    //                      iTkin    
                                                   >> 197       }      //         jMat != iMat
                                                   >> 198     }        //     jMat
                                                   >> 199   }          // iMat
                                                   >> 200 }
                                                   >> 201 
                                                   >> 202 
                                                   >> 203 ****************************************************************  */
                                                   >> 204 
 48                                                   205 
 49 //////////////////////////////////////////////    206 //////////////////////////////////////////////////////////////////////
 50 //                                                207 //
 51 // Constructor for creation of physics tables  << 208 // Constructor for creation of physics tables (angle and energy TR 
 52 // distributions) for a couple of selected mat    209 // distributions) for a couple of selected materials.
 53 //                                                210 //
 54 // Recommended for use in applications with ma    211 // Recommended for use in applications with many materials involved,
 55 // when only few (usually couple) materials ar    212 // when only few (usually couple) materials are interested for generation
 56 // of TR on the interface between them            213 // of TR on the interface between them
 57 G4ForwardXrayTR::G4ForwardXrayTR(const G4Strin << 214 
 58                                  const G4Strin << 215 
 59                                  const G4Strin << 216 G4ForwardXrayTR::
 60   : G4TransitionRadiation(processName)         << 217 G4ForwardXrayTR( const G4String& matName1,   //  G4Material* pMat1,
 61 {                                              << 218      const G4String& matName2,    //  G4Material* pMat2,
 62   secID = G4PhysicsModelCatalog::GetModelID("m << 219                  const G4String& processName                          )
 63   fPtrGamma                = nullptr;          << 220   :        G4TransitionRadiation(processName) 
 64   fGammaCutInKineticEnergy = nullptr;          << 221 {
 65   fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR  << 222   //  fMatIndex1 = pMat1->GetIndex() ;
 66   fGamma = fSigma1 = fSigma2 = 0.0;            << 223   //  fMatIndex2 = pMat2->GetIndex() ;
 67   fAngleDistrTable           = nullptr;        << 
 68   fEnergyDistrTable          = nullptr;        << 
 69   fMatIndex1 = fMatIndex2 = 0;                 << 
 70                                                << 
 71   // Proton energy vector initialization       << 
 72   fProtonEnergyVector =                        << 
 73     new G4PhysicsLogVector(fMinProtonTkin, fMa << 
 74   G4int iMat;                                     224   G4int iMat;
 75   const G4ProductionCutsTable* theCoupleTable  << 225   static 
 76     G4ProductionCutsTable::GetProductionCutsTa << 226   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
 77   G4int numOfCouples = (G4int)theCoupleTable-> << 
 78                                                   227 
 79   G4bool build = true;                         << 228   G4int numOfMat = G4Material::GetNumberOfMaterials();
 80                                                   229 
 81   for(iMat = 0; iMat < numOfCouples; ++iMat)   << 230   for(iMat=0;iMat<numOfMat;iMat++)    // check first material name
 82   {                                               231   {
 83     const G4MaterialCutsCouple* couple =       << 232     if( matName1 == (*theMaterialTable)[iMat]->GetName() )
 84       theCoupleTable->GetMaterialCutsCouple(iM << 
 85     if(matName1 == couple->GetMaterial()->GetN << 
 86     {                                             233     {
 87       fMatIndex1 = couple->GetIndex();         << 234       fMatIndex1 = (*theMaterialTable)[iMat]->GetIndex() ;
 88       break;                                   << 235       break ;
 89     }                                             236     }
 90   }                                               237   }
 91   if(iMat == numOfCouples)                     << 238   if(iMat == numOfMat)
 92   {                                               239   {
 93     G4Exception("G4ForwardXrayTR::G4ForwardXra << 240     G4Exception("Invalid first material name in G4ForwardXrayTR constructor") ;
 94                 JustWarning,                   << 
 95                 "Invalid first material name i << 
 96     build = false;                             << 
 97   }                                               241   }
 98                                                   242 
 99   if(build)                                    << 243   for(iMat=0;iMat<numOfMat;iMat++)    // check second material name
100   {                                               244   {
101     for(iMat = 0; iMat < numOfCouples; ++iMat) << 245     if( matName2 == (*theMaterialTable)[iMat]->GetName() )
102     {                                             246     {
103       const G4MaterialCutsCouple* couple =     << 247       fMatIndex2 = (*theMaterialTable)[iMat]->GetIndex() ;
104         theCoupleTable->GetMaterialCutsCouple( << 248       break ;
105       if(matName2 == couple->GetMaterial()->Ge << 
106       {                                        << 
107         fMatIndex2 = couple->GetIndex();       << 
108         break;                                 << 
109       }                                        << 
110     }                                          << 
111     if(iMat == numOfCouples)                   << 
112     {                                          << 
113       G4Exception(                             << 
114         "G4ForwardXrayTR::G4ForwardXrayTR", "F << 
115         "Invalid second material name in G4For << 
116       build = false;                           << 
117     }                                             249     }
118   }                                               250   }
119   if(build)                                    << 251   if(iMat == numOfMat)
120   {                                               252   {
121     BuildXrayTRtables();                       << 253     G4Exception("Invalid second material name in G4ForwardXrayTR constructor") ;
122   }                                               254   }
                                                   >> 255   //  G4cout<<"G4ForwardXray constructor is called"<<G4endl ;
                                                   >> 256   BuildXrayTRtables() ;
123 }                                                 257 }
124                                                   258 
125 //////////////////////////////////////////////    259 /////////////////////////////////////////////////////////////////////////
                                                   >> 260 //
126 // Constructor used by X-ray transition radiat    261 // Constructor used by X-ray transition radiation parametrisation models
127 G4ForwardXrayTR::G4ForwardXrayTR(const G4Strin << 262 
128   : G4TransitionRadiation(processName)         << 263 G4ForwardXrayTR::
                                                   >> 264 G4ForwardXrayTR( const G4String& processName  )
                                                   >> 265   :        G4TransitionRadiation(processName) 
129 {                                                 266 {
130   fPtrGamma                = nullptr;          << 267   ;
131   fGammaCutInKineticEnergy = nullptr;          << 
132   fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR  << 
133   fGamma = fSigma1 = fSigma2 = 0.0;            << 
134   fAngleDistrTable           = nullptr;        << 
135   fEnergyDistrTable          = nullptr;        << 
136   fMatIndex1 = fMatIndex2 = 0;                 << 
137                                                << 
138   // Proton energy vector initialization       << 
139   fProtonEnergyVector =                        << 
140     new G4PhysicsLogVector(fMinProtonTkin, fMa << 
141 }                                                 268 }
142                                                   269 
                                                   >> 270 
143 //////////////////////////////////////////////    271 //////////////////////////////////////////////////////////////////////
                                                   >> 272 //
144 // Destructor                                     273 // Destructor
145 G4ForwardXrayTR::~G4ForwardXrayTR()            << 274 //
146 {                                              << 
147   delete fAngleDistrTable;                     << 
148   delete fEnergyDistrTable;                    << 
149   delete fProtonEnergyVector;                  << 
150 }                                              << 
151                                                << 
152 void G4ForwardXrayTR::ProcessDescription(std:: << 
153 {                                              << 
154   out << "Simulation of forward X-ray transiti << 
155          "relativistic charged particles cross << 
156          "two materials.\n";                   << 
157 }                                              << 
158                                                   275 
159 G4double G4ForwardXrayTR::GetMeanFreePath(cons << 276 G4ForwardXrayTR::~G4ForwardXrayTR()
160                                           G4Fo << 
161 {                                                 277 {
162   *condition = Forced;                         << 278   ;
163   return DBL_MAX;  // so TR doesn't limit mean << 
164 }                                                 279 }
165                                                   280 
166 //////////////////////////////////////////////    281 //////////////////////////////////////////////////////////////////////////////
                                                   >> 282 //
167 // Build physics tables for energy and angular    283 // Build physics tables for energy and angular distributions of X-ray TR photon
                                                   >> 284 
168 void G4ForwardXrayTR::BuildXrayTRtables()         285 void G4ForwardXrayTR::BuildXrayTRtables()
169 {                                                 286 {
170   G4int iMat, jMat, iTkin, iTR, iPlace;        << 287   G4int iMat, jMat, iTkin, iTR, iPlace ;
171   const G4ProductionCutsTable* theCoupleTable  << 288   static 
172     G4ProductionCutsTable::GetProductionCutsTa << 289   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
173   G4int numOfCouples = (G4int)theCoupleTable-> << 290 
                                                   >> 291   G4int numOfMat = G4Material::GetNumberOfMaterials() ;
                                                   >> 292 
                                                   >> 293   fGammaCutInKineticEnergy = new G4double[numOfMat] ;
                                                   >> 294   fGammaCutInKineticEnergy = fPtrGamma->GetEnergyCuts() ;
174                                                   295 
175   fGammaCutInKineticEnergy = theCoupleTable->G << 296   fAngleDistrTable  = new G4PhysicsTable(2*fTotBin) ;
                                                   >> 297   fEnergyDistrTable = new G4PhysicsTable(2*fTotBin) ;
176                                                   298 
177   fAngleDistrTable  = new G4PhysicsTable(2 * f << 
178   fEnergyDistrTable = new G4PhysicsTable(2 * f << 
179                                                   299 
180   for(iMat = 0; iMat < numOfCouples;           << 300   for(iMat=0;iMat<numOfMat;iMat++)     // loop over pairs of different materials
181       ++iMat)  // loop over pairs of different << 
182   {                                               301   {
183     if(iMat != fMatIndex1 && iMat != fMatIndex << 302     if( iMat != fMatIndex1 && iMat != fMatIndex2 ) continue ;
184       continue;                                << 
185                                                   303 
186     for(jMat = 0; jMat < numOfCouples; ++jMat) << 304     for(jMat=0;jMat<numOfMat;jMat++)  // transition iMat -> jMat !!!
187     {                                             305     {
188       if(iMat == jMat || (jMat != fMatIndex1 & << 306       if( iMat == jMat || ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )   
189       {                                        << 307       { 
190         continue;                              << 308         continue ; 
191       }                                        << 309       } 
192       else                                        310       else
193       {                                           311       {
194         const G4MaterialCutsCouple* iCouple =  << 312         const G4Material* mat1 = (*theMaterialTable)[iMat] ;
195           theCoupleTable->GetMaterialCutsCoupl << 313         const G4Material* mat2 = (*theMaterialTable)[jMat] ;
196         const G4MaterialCutsCouple* jCouple =  << 
197           theCoupleTable->GetMaterialCutsCoupl << 
198         const G4Material* mat1 = iCouple->GetM << 
199         const G4Material* mat2 = jCouple->GetM << 
200                                                   314 
201         fSigma1 = fPlasmaCof * (mat1->GetElect << 315         fSigma1 = fPlasmaCof*(mat1->GetElectronDensity()) ;
202         fSigma2 = fPlasmaCof * (mat2->GetElect << 316         fSigma2 = fPlasmaCof*(mat2->GetElectronDensity()) ;
203                                                   317 
204         fGammaTkinCut = 0.0;                   << 318         // fGammaTkinCut = fGammaCutInKineticEnergy[jMat] ; // TR photon in jMat !
205                                                   319 
206         if(fGammaTkinCut > fTheMinEnergyTR)  / << 320         fGammaTkinCut = 0.0 ;
207         {                                      << 321 
208           fMinEnergyTR = fGammaTkinCut;        << 322         if(fGammaTkinCut > fTheMinEnergyTR)    // setting of min/max TR energies 
209         }                                      << 323   {
                                                   >> 324           fMinEnergyTR = fGammaTkinCut ;
                                                   >> 325   }
210         else                                      326         else
211         {                                      << 327   {
212           fMinEnergyTR = fTheMinEnergyTR;      << 328           fMinEnergyTR = fTheMinEnergyTR ;
213         }                                      << 329   }
214         if(fGammaTkinCut > fTheMaxEnergyTR)       330         if(fGammaTkinCut > fTheMaxEnergyTR)
215         {                                      << 331   {
216           fMaxEnergyTR = 2.0 * fGammaTkinCut;  << 332           fMaxEnergyTR = 2.0*fGammaTkinCut ;    // usually very low TR rate 
217         }                                      << 333   }
218         else                                      334         else
219         {                                      << 335   {
220           fMaxEnergyTR = fTheMaxEnergyTR;      << 336           fMaxEnergyTR = fTheMaxEnergyTR ;
221         }                                      << 337   }
222         for(iTkin = 0; iTkin < fTotBin; ++iTki << 338         for(iTkin=0;iTkin<fTotBin;iTkin++)      // Lorentz factor loop
223         {                                      << 339   {
224           auto energyVector =                  << 340           G4PhysicsLogVector* 
225             new G4PhysicsLogVector(fMinEnergyT << 341                     energyVector = new G4PhysicsLogVector( fMinEnergyTR,
226                                                << 342                                                            fMaxEnergyTR,
227           fGamma = 1.0 + (fProtonEnergyVector- << 343                                                            fBinTR         ) ;
228                           proton_mass_c2);     << 
229                                                   344 
230           fMaxThetaTR = 10000.0 / (fGamma * fG << 345           fGamma = 1.0 +   (fProtonEnergyVector->
                                                   >> 346                             GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
231                                                   347 
                                                   >> 348           fMaxThetaTR = 10000.0/(fGamma*fGamma) ;
                                                   >> 349  
232           if(fMaxThetaTR > fTheMaxAngle)          350           if(fMaxThetaTR > fTheMaxAngle)
233           {                                    << 351           { 
234             fMaxThetaTR = fTheMaxAngle;        << 352             fMaxThetaTR = fTheMaxAngle ;
235           }                                    << 353     }
236           else                                    354           else
237           {                                    << 355     {
238             if(fMaxThetaTR < fTheMinAngle)        356             if(fMaxThetaTR < fTheMinAngle)
239             {                                  << 357       {
240               fMaxThetaTR = fTheMinAngle;      << 358               fMaxThetaTR = fTheMinAngle ;
241             }                                  << 359       }
242           }                                    << 360     }
243           auto angleVector =                   << 361    // G4cout<<G4endl<<"fGamma = "<<fGamma<<"  fMaxThetaTR = "<<fMaxThetaTR<<G4endl ;
244             new G4PhysicsLinearVector(0.0, fMa << 362           G4PhysicsLinearVector* 
245           G4double energySum = 0.0;            << 363                      angleVector = new G4PhysicsLinearVector(        0.0,
246           G4double angleSum  = 0.0;            << 364                                                              fMaxThetaTR,
247                                                << 365                                                                   fBinTR  ) ;
248           energyVector->PutValue(fBinTR - 1, e << 366           G4double energySum = 0.0 ;
249           angleVector->PutValue(fBinTR - 1, an << 367           G4double angleSum  = 0.0 ;
250                                                << 368 
251           for(iTR = fBinTR - 2; iTR >= 0; --iT << 369           energyVector->PutValue(fBinTR-1,energySum) ;
252           {                                    << 370           angleVector->PutValue(fBinTR-1,angleSum)   ;
253             energySum +=                       << 371 
254               fCofTR * EnergySum(energyVector- << 372           for(iTR=fBinTR-2;iTR>=0;iTR--)
255                                  energyVector- << 373     {
256                                                << 374             energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR),
257             angleSum +=                        << 375                                         energyVector->GetLowEdgeEnergy(iTR+1)) ; 
258               fCofTR * AngleSum(angleVector->G << 376 
259                                 angleVector->G << 377             angleSum  += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR),
260                                                << 378                                          angleVector->GetLowEdgeEnergy(iTR+1)) ;
261             energyVector->PutValue(iTR, energy << 379 
262             angleVector->PutValue(iTR, angleSu << 380             energyVector->PutValue(iTR,energySum) ;
263           }                                    << 381             angleVector ->PutValue(iTR,angleSum)   ;
                                                   >> 382     }
                                                   >> 383     // G4cout<<"sumE = "<<energySum<<" ; sumA = "<<angleSum<<G4endl ;
264                                                   384 
265           if(jMat < iMat)                         385           if(jMat < iMat)
266           {                                    << 386     {
267             iPlace = fTotBin + iTkin;          << 387             iPlace = fTotBin+iTkin ;   // (iMat*(numOfMat-1)+jMat)*
268           }                                    << 388     }
269           else  // jMat > iMat right part of m << 389           else   // jMat > iMat right part of matrices (jMat-1) ! 
270           {                                    << 390     {
271             iPlace = iTkin;                    << 391             iPlace = iTkin ;  // (iMat*(numOfMat-1)+jMat-1)*fTotBin+
272           }                                    << 392     } 
273           fEnergyDistrTable->insertAt(iPlace,  << 393           fEnergyDistrTable->insertAt(iPlace,energyVector) ;
274           fAngleDistrTable->insertAt(iPlace, a << 394           fAngleDistrTable->insertAt(iPlace,angleVector) ;
275         }  //                      iTkin       << 395   }    //                      iTkin    
276       }    //         jMat != iMat             << 396       }      //         jMat != iMat
277     }      //     jMat                         << 397     }        //     jMat
278   }        // iMat                             << 398   }          // iMat
                                                   >> 399   //  G4cout<<"G4ForwardXrayTR::BuildXrayTRtables have been called"<<G4endl ;
279 }                                                 400 }
280                                                   401 
281 //////////////////////////////////////////////    402 ///////////////////////////////////////////////////////////////////////
282 //                                                403 //
283 // This function returns the spectral and angl    404 // This function returns the spectral and angle density of TR quanta
284 // in X-ray energy region generated forward wh    405 // in X-ray energy region generated forward when a relativistic
285 // charged particle crosses interface between     406 // charged particle crosses interface between two materials.
286 // The high energy small theta approximation i    407 // The high energy small theta approximation is applied.
287 // (matter1 -> matter2)                           408 // (matter1 -> matter2)
288 // varAngle =2* (1 - std::cos(Theta)) or appro << 409 // varAngle =2* (1 - cos(Theta)) or approximately = Theta*Theta
289 //                                                410 //
290 G4double G4ForwardXrayTR::SpectralAngleTRdensi << 411 
291                                                << 412 G4double
292 {                                              << 413 G4ForwardXrayTR::SpectralAngleTRdensity( G4double energy,
293   G4double formationLength1, formationLength2; << 414                                          G4double varAngle ) const
294   formationLength1 =                           << 415 {
295     1.0 / (1.0 / (fGamma * fGamma) + fSigma1 / << 416   G4double  formationLength1, formationLength2 ;
296   formationLength2 =                           << 417   formationLength1 = 1.0/
297     1.0 / (1.0 / (fGamma * fGamma) + fSigma2 / << 418   (1.0/(fGamma*fGamma)
298   return (varAngle / energy) * (formationLengt << 419   + fSigma1/(energy*energy)
299          (formationLength1 - formationLength2) << 420   + varAngle) ;
                                                   >> 421   formationLength2 = 1.0/
                                                   >> 422   (1.0/(fGamma*fGamma)
                                                   >> 423   + fSigma2/(energy*energy)
                                                   >> 424   + varAngle) ;
                                                   >> 425   return (varAngle/energy)*(formationLength1 - formationLength2)
                                                   >> 426               *(formationLength1 - formationLength2)  ;
                                                   >> 427 
300 }                                                 428 }
301                                                   429 
                                                   >> 430 
302 //////////////////////////////////////////////    431 //////////////////////////////////////////////////////////////////
                                                   >> 432 //
303 // Analytical formula for angular density of X    433 // Analytical formula for angular density of X-ray TR photons
304 G4double G4ForwardXrayTR::AngleDensity(G4doubl << 434 // 
                                                   >> 435 
                                                   >> 436 G4double G4ForwardXrayTR::AngleDensity( G4double energy,
                                                   >> 437                                         G4double varAngle ) const
305 {                                                 438 {
306   G4double x, x2, c, d, f, a2, b2, a4, b4;     << 439   G4double x, x2, a, b, c, d, f, a2, b2, a4, b4 ;
307   G4double cof1, cof2, cof3;                   << 440   G4double cof1, cof2, cof3 ;
308   x    = 1.0 / energy;                         << 441   x = 1.0/energy ;
309   x2   = x * x;                                << 442   x2 = x*x ;
310   c    = 1.0 / fSigma1;                        << 443   c = 1.0/fSigma1 ;
311   d    = 1.0 / fSigma2;                        << 444   d = 1.0/fSigma2 ;
312   f    = (varAngle + 1.0 / (fGamma * fGamma)); << 445   f = (varAngle + 1.0/(fGamma*fGamma)) ;
313   a2   = c * f;                                << 446   a2 = c*f ;
314   b2   = d * f;                                << 447   b2 = d*f ;
315   a4   = a2 * a2;                              << 448   a4 = a2*a2 ;
316   b4   = b2 * b2;                              << 449   b4 = b2*b2 ;
317   cof1 = c * c * (0.5 / (a2 * (x2 + a2)) + 0.5 << 450   a = sqrt(a2) ;
318   cof3 = d * d * (0.5 / (b2 * (x2 + b2)) + 0.5 << 451   b = sqrt(b2) ;
319   cof2 = -c * d *                              << 452   cof1 = c*c*(0.5/(a2*(x2 +a2)) +0.5*log(x2/(x2 +a2))/a4) ;
320          (std::log(x2 / (x2 + b2)) / b2 - std: << 453   cof3 = d*d*(0.5/(b2*(x2 +b2)) +0.5*log(x2/(x2 +b2))/b4) ;
321          (a2 - b2);                            << 454   cof2 = -c*d*(log(x2/(x2 +b2))/b2 - log(x2/(x2 +a2))/a2)/(a2 - b2)   ;
322   return -varAngle * (cof1 + cof2 + cof3);     << 455   return -varAngle*(cof1 + cof2 + cof3) ;
323 }                                                 456 }
324                                                   457 
325 //////////////////////////////////////////////    458 /////////////////////////////////////////////////////////////////////
                                                   >> 459 //
326 // Definite integral of X-ray TR spectral-angl    460 // Definite integral of X-ray TR spectral-angle density from energy1
327 // to energy2                                     461 // to energy2
328 G4double G4ForwardXrayTR::EnergyInterval(G4dou << 462 //
329                                          G4dou << 463 
                                                   >> 464 G4double G4ForwardXrayTR::EnergyInterval( G4double energy1,
                                                   >> 465                                           G4double energy2,
                                                   >> 466                                           G4double varAngle ) const
330 {                                                 467 {
331   return AngleDensity(energy2, varAngle) - Ang << 468   return     AngleDensity(energy2,varAngle)
                                                   >> 469            - AngleDensity(energy1,varAngle) ;
332 }                                                 470 }
333                                                   471 
334 //////////////////////////////////////////////    472 //////////////////////////////////////////////////////////////////////
                                                   >> 473 //
335 // Integral angle distribution of X-ray TR pho    474 // Integral angle distribution of X-ray TR photons based on analytical
336 // formula for angle density                      475 // formula for angle density
337 G4double G4ForwardXrayTR::AngleSum(G4double va << 476 //
                                                   >> 477 
                                                   >> 478 G4double G4ForwardXrayTR::AngleSum( G4double varAngle1,
                                                   >> 479                                     G4double varAngle2     )   const
338 {                                                 480 {
339   G4int i;                                     << 481   G4int i ;
340   G4double h, sumEven = 0.0, sumOdd = 0.0;     << 482   G4double h , sumEven = 0.0 , sumOdd = 0.0 ;
341   h = 0.5 * (varAngle2 - varAngle1) / fSympson << 483   h = 0.5*(varAngle2 - varAngle1)/fSympsonNumber ;
342   for(i = 1; i < fSympsonNumber; ++i)          << 484   for(i=1;i<fSympsonNumber;i++)
343   {                                            << 485   {
344     sumEven +=                                 << 486    sumEven += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1 + 2*i*h   ) ;
345       EnergyInterval(fMinEnergyTR, fMaxEnergyT << 487    sumOdd  += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,
346     sumOdd +=                                  << 488                                                    varAngle1 + (2*i - 1)*h ) ;
347       EnergyInterval(fMinEnergyTR, fMaxEnergyT << 489   }
348   }                                            << 490   sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,
349   sumOdd += EnergyInterval(fMinEnergyTR, fMaxE << 491                         varAngle1 + (2*fSympsonNumber - 1)*h    ) ;
350                            varAngle1 + (2 * fS << 492 
351                                                << 493   return h*(EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1)
352   return h *                                   << 494           + EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle2)
353          (EnergyInterval(fMinEnergyTR, fMaxEne << 495           + 4.0*sumOdd + 2.0*sumEven                          )/3.0 ;
354           EnergyInterval(fMinEnergyTR, fMaxEne << 
355           2.0 * sumEven) /                     << 
356          3.0;                                  << 
357 }                                                 496 }
358                                                   497 
359 //////////////////////////////////////////////    498 /////////////////////////////////////////////////////////////////////
                                                   >> 499 //
360 // Analytical Expression for   spectral densit    500 // Analytical Expression for   spectral density of Xray TR photons
361 // x = 2*(1 - std::cos(Theta)) ~ Theta^2       << 501 // x = 2*(1 - cos(Theta)) ~ Theta^2
362 G4double G4ForwardXrayTR::SpectralDensity(G4do << 502 //
                                                   >> 503 
                                                   >> 504 G4double G4ForwardXrayTR::SpectralDensity( G4double energy,
                                                   >> 505                                            G4double      x  ) const
363 {                                                 506 {
364   G4double a, b;                               << 507   G4double a, b ;
365   a = 1.0 / (fGamma * fGamma) + fSigma1 / (ene << 508   a =  1.0/(fGamma*fGamma)
366   b = 1.0 / (fGamma * fGamma) + fSigma2 / (ene << 509      + fSigma1/(energy*energy)  ;
367   return ((a + b) * std::log((x + b) / (x + a) << 510   b =  1.0/(fGamma*fGamma)
368           b / (x + b)) /                       << 511      + fSigma2/(energy*energy)  ;
369          energy;                               << 512   return ( (a + b)*log((x + b)/(x + a))/(a - b)
                                                   >> 513           + a/(x + a) + b/(x + b) )/energy ;
                                                   >> 514 
370 }                                                 515 }
371                                                   516 
372 //////////////////////////////////////////////    517 ////////////////////////////////////////////////////////////////////
                                                   >> 518 //
373 //  The spectral density in some angle interva    519 //  The spectral density in some angle interval from varAngle1 to
374 //  varAngle2                                     520 //  varAngle2
375 G4double G4ForwardXrayTR::AngleInterval(G4doub << 521 //
376                                         G4doub << 522 
                                                   >> 523 G4double G4ForwardXrayTR::AngleInterval( G4double    energy,
                                                   >> 524                                          G4double varAngle1,
                                                   >> 525                                          G4double varAngle2   ) const
377 {                                                 526 {
378   return SpectralDensity(energy, varAngle2) -  << 527   return     SpectralDensity(energy,varAngle2)
379          SpectralDensity(energy, varAngle1);   << 528            - SpectralDensity(energy,varAngle1) ;
380 }                                                 529 }
381                                                   530 
382 //////////////////////////////////////////////    531 ////////////////////////////////////////////////////////////////////
                                                   >> 532 //
383 // Integral spectral distribution of X-ray TR     533 // Integral spectral distribution of X-ray TR photons based on
384 // analytical formula for spectral density        534 // analytical formula for spectral density
385 G4double G4ForwardXrayTR::EnergySum(G4double e << 535 //
                                                   >> 536 
                                                   >> 537 G4double G4ForwardXrayTR::EnergySum( G4double energy1,
                                                   >> 538                                      G4double energy2     )   const
386 {                                                 539 {
387   G4int i;                                     << 540   G4int i ;
388   G4double h, sumEven = 0.0, sumOdd = 0.0;     << 541   G4double h , sumEven = 0.0 , sumOdd = 0.0 ;
389   h = 0.5 * (energy2 - energy1) / fSympsonNumb << 542   h = 0.5*(energy2 - energy1)/fSympsonNumber ;
390   for(i = 1; i < fSympsonNumber; ++i)          << 543   for(i=1;i<fSympsonNumber;i++)
391   {                                            << 544   {
392     sumEven += AngleInterval(energy1 + 2 * i * << 545    sumEven += AngleInterval(energy1 + 2*i*h,0.0,fMaxThetaTR);
393     sumOdd += AngleInterval(energy1 + (2 * i - << 546    sumOdd  += AngleInterval(energy1 + (2*i - 1)*h,0.0,fMaxThetaTR) ;
394   }                                            << 547   }
395   sumOdd +=                                    << 548   sumOdd += AngleInterval(energy1 + (2*fSympsonNumber - 1)*h,
396     AngleInterval(energy1 + (2 * fSympsonNumbe << 549                           0.0,fMaxThetaTR) ;
397                                                << 550 
398   return h *                                   << 551   return h*(  AngleInterval(energy1,0.0,fMaxThetaTR)
399          (AngleInterval(energy1, 0.0, fMaxThet << 552             + AngleInterval(energy2,0.0,fMaxThetaTR)
400           AngleInterval(energy2, 0.0, fMaxThet << 553             + 4.0*sumOdd + 2.0*sumEven                          )/3.0 ;
401           2.0 * sumEven) /                     << 
402          3.0;                                  << 
403 }                                                 554 }
404                                                   555 
405 //////////////////////////////////////////////    556 /////////////////////////////////////////////////////////////////////////
                                                   >> 557 //
406 // PostStepDoIt function for creation of forwa    558 // PostStepDoIt function for creation of forward X-ray photons in TR process
407 // on boundary between two materials with real << 559 // on boubndary between two materials with really different plasma energies
                                                   >> 560 //
                                                   >> 561 
408 G4VParticleChange* G4ForwardXrayTR::PostStepDo    562 G4VParticleChange* G4ForwardXrayTR::PostStepDoIt(const G4Track& aTrack,
409                                                << 563                         const G4Step&  aStep)
410 {                                                 564 {
411   aParticleChange.Initialize(aTrack);             565   aParticleChange.Initialize(aTrack);
412   G4int iMat, jMat, iTkin, iPlace, numOfTR, iT << 566   //  G4cout<<"call G4ForwardXrayTR::PostStepDoIt"<<G4endl ;
                                                   >> 567   G4int iMat, jMat, iTkin, iPlace, numOfMat, numOfTR, iTR, iTransfer ;
413                                                   568 
414   G4double energyPos, anglePos, energyTR, thet << 569   G4double energyPos, anglePos, energyTR, theta, phi, dirX, dirY, dirZ ;
415   G4double W, W1, W2, E1, E2;                  << 570   G4double W, W1, W2, E1, E2 ;
                                                   >> 571   static 
                                                   >> 572   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
                                                   >> 573   numOfMat = G4Material::GetNumberOfMaterials() ;
416                                                   574 
417   G4StepPoint* pPreStepPoint  = aStep.GetPreSt    575   G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
418   G4StepPoint* pPostStepPoint = aStep.GetPostS    576   G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
419   G4double tol =                               << 
420     0.5 * G4GeometryTolerance::GetInstance()-> << 
421                                                   577 
422   if(pPostStepPoint->GetStepStatus() != fGeomB << 578   if (pPostStepPoint->GetStepStatus() != fGeomBoundary)
423   {                                               579   {
424     return G4VDiscreteProcess::PostStepDoIt(aT    580     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
425   }                                               581   }
426   if(aTrack.GetStepLength() <= tol)            << 582   if (aTrack.GetStepLength() <= kCarTolerance*0.5)
                                                   >> 583 
427   {                                               584   {
428     return G4VDiscreteProcess::PostStepDoIt(aT    585     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
429   }                                               586   }
430   // Arrived at boundary, so begin to try TR   << 587   // Come on boundary, so begin to try TR
431                                                   588 
432   const G4MaterialCutsCouple* iCouple = pPreSt << 589   iMat = pPreStepPoint ->GetPhysicalVolume()->
433                                           ->Ge << 590        GetLogicalVolume()->GetMaterial()->GetIndex() ;
434                                           ->Ge << 591   jMat = pPostStepPoint->GetPhysicalVolume()->
435   const G4MaterialCutsCouple* jCouple = pPostS << 592        GetLogicalVolume()->GetMaterial()->GetIndex() ;
436                                           ->Ge << 
437                                           ->Ge << 
438   const G4Material* iMaterial = iCouple->GetMa << 
439   const G4Material* jMaterial = jCouple->GetMa << 
440   iMat                        = iCouple->GetIn << 
441   jMat                        = jCouple->GetIn << 
442                                                   593 
443   // The case of equal or approximate (in term    594   // The case of equal or approximate (in terms of plasma energy) materials
444   // No TR photons ?!                             595   // No TR photons ?!
445                                                   596 
446   if(iMat == jMat ||                           << 597   if (     iMat == jMat
447      ((fMatIndex1 >= 0 && fMatIndex2 >= 0) &&  << 598       || (    (fMatIndex1 >= 0 && fMatIndex1 >= 0) 
448       (iMat != fMatIndex1 && iMat != fMatIndex << 599            && ( iMat != fMatIndex1 && iMat != fMatIndex2 )
449       (jMat != fMatIndex1 && jMat != fMatIndex << 600            && ( jMat != fMatIndex1 && jMat != fMatIndex2 )  )
450                                                << 601 
451      || iMaterial->GetState() == jMaterial->Ge << 602       || (*theMaterialTable)[iMat]->GetState() ==
452                                                << 603          (*theMaterialTable)[jMat]->GetState()
453      || (iMaterial->GetState() == kStateSolid  << 604   
454          jMaterial->GetState() == kStateLiquid << 605       ||(    (*theMaterialTable)[iMat]->GetState() == kStateSolid
455                                                << 606           && (*theMaterialTable)[jMat]->GetState() == kStateLiquid )
456      || (iMaterial->GetState() == kStateLiquid << 607  
457          jMaterial->GetState() == kStateSolid) << 608       ||(    (*theMaterialTable)[iMat]->GetState() == kStateLiquid
                                                   >> 609           && (*theMaterialTable)[jMat]->GetState() == kStateSolid  )   )
458   {                                               610   {
459     return G4VDiscreteProcess::PostStepDoIt(aT << 611     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep) ;
460   }                                               612   }
461                                                   613 
462   const G4DynamicParticle* aParticle = aTrack.    614   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
463   G4double charge = aParticle->GetDefinition()    615   G4double charge = aParticle->GetDefinition()->GetPDGCharge();
464                                                << 616  
465   if(charge == 0.0)  // Uncharged particle doe << 617   if(charge == 0.0) // Uncharged particle doesn't Generate TR photons 
466   {                                               618   {
467     return G4VDiscreteProcess::PostStepDoIt(aT    619     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
468   }                                               620   }
469   // Now we are ready to Generate TR photons      621   // Now we are ready to Generate TR photons
470                                                   622 
471   G4double chargeSq  = charge * charge;        << 623   G4double chargeSq = charge*charge ;
472   G4double kinEnergy = aParticle->GetKineticEn << 624   G4double kinEnergy     = aParticle->GetKineticEnergy() ;
473   G4double massRatio =                         << 625   G4double massRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ;
474     proton_mass_c2 / aParticle->GetDefinition( << 626   G4double TkinScaled = kinEnergy*massRatio ;
475   G4double TkinScaled = kinEnergy * massRatio; << 627   for(iTkin=0;iTkin<fTotBin;iTkin++)
476   for(iTkin = 0; iTkin < fTotBin; ++iTkin)     << 
477   {                                               628   {
478     if(TkinScaled < fProtonEnergyVector->GetLo << 629     if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) // <= ?
479     {                                             630     {
480       break;                                   << 631       break ;
481     }                                          << 632     }     
482   }                                               633   }
483   if(jMat < iMat)                                 634   if(jMat < iMat)
484   {                                               635   {
485     iPlace = fTotBin + iTkin - 1;              << 636     iPlace = fTotBin + iTkin - 1 ; // (iMat*(numOfMat - 1) + jMat)*
486   }                                               637   }
487   else                                            638   else
488   {                                               639   {
489     iPlace = iTkin - 1;                        << 640     iPlace = iTkin - 1 ;  // (iMat*(numOfMat - 1) + jMat - 1)*fTotBin + 
490   }                                               641   }
                                                   >> 642   //  G4PhysicsVector*  energyVector1 = (*fEnergyDistrTable)(iPlace)     ;
                                                   >> 643   //  G4PhysicsVector*  energyVector2 = (*fEnergyDistrTable)(iPlace + 1) ;
491                                                   644 
492   G4ParticleMomentum particleDir = aParticle-> << 645   //  G4PhysicsVector*   angleVector1 = (*fAngleDistrTable)(iPlace)      ;
                                                   >> 646   //  G4PhysicsVector*   angleVector2 = (*fAngleDistrTable)(iPlace + 1)  ;
493                                                   647 
494   if(iTkin == fTotBin)  // TR plato, try from  << 648   G4ParticleMomentum particleDir = aParticle->GetMomentumDirection() ;
                                                   >> 649 
                                                   >> 650   if(iTkin == fTotBin)                 // TR plato, try from left
495   {                                               651   {
496     numOfTR = (G4int)G4Poisson(                << 652  // G4cout<<iTkin<<" mean TR number = "<<( (*(*fEnergyDistrTable)(iPlace))(0) +
497       ((*(*fEnergyDistrTable)(iPlace))(0) + (* << 653  //                                   (*(*fAngleDistrTable)(iPlace))(0) )
498       chargeSq * 0.5);                         << 654  //      *chargeSq*0.5<<G4endl ;
                                                   >> 655 
                                                   >> 656     numOfTR = RandPoisson::shoot( ( (*(*fEnergyDistrTable)(iPlace))(0) +
                                                   >> 657                                     (*(*fAngleDistrTable)(iPlace))(0) )
                                                   >> 658                                     *chargeSq*0.5 ) ;
499     if(numOfTR == 0)                              659     if(numOfTR == 0)
500     {                                             660     {
501       return G4VDiscreteProcess::PostStepDoIt(    661       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
502     }                                             662     }
503     else                                          663     else
504     {                                             664     {
                                                   >> 665       // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl ;
                                                   >> 666 
505       aParticleChange.SetNumberOfSecondaries(n    667       aParticleChange.SetNumberOfSecondaries(numOfTR);
506                                                   668 
507       for(iTR = 0; iTR < numOfTR; ++iTR)       << 669       for(iTR=0;iTR<numOfTR;iTR++)
508       {                                           670       {
509         energyPos = (*(*fEnergyDistrTable)(iPl << 671         energyPos = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand() ;
510         for(iTransfer = 0; iTransfer < fBinTR  << 672         for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
511         {                                      << 673   {
512           if(energyPos >= (*(*fEnergyDistrTabl << 674           if(energyPos >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break ;
513             break;                             << 675   }
514         }                                      << 676         energyTR = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
515         energyTR = (*fEnergyDistrTable)(iPlace << 677 
516                                                << 678   // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl ;
517         kinEnergy -= energyTR;                 << 679 
518         aParticleChange.ProposeEnergy(kinEnerg << 680         kinEnergy -= energyTR ; 
519                                                << 681         aParticleChange.SetEnergyChange(kinEnergy);
520         anglePos = (*(*fAngleDistrTable)(iPlac << 682 
521         for(iTransfer = 0; iTransfer < fBinTR  << 683         anglePos = (*(*fAngleDistrTable)(iPlace))(0)*G4UniformRand() ;
522         {                                      << 684         for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
523           if(anglePos > (*(*fAngleDistrTable)( << 685   {
524             break;                             << 686           if(anglePos > (*(*fAngleDistrTable)(iPlace))(iTransfer)) break ;
525         }                                      << 687   }
526         theta = std::sqrt(                     << 688         theta = sqrt((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1)) ;
527           (*fAngleDistrTable)(iPlace)->GetLowE << 689 
528                                                << 690   // G4cout<<iTransfer<<" :  theta = "<<theta<<G4endl ;
529         phi  = twopi * G4UniformRand();        << 691 
530         dirX = std::sin(theta) * std::cos(phi) << 692         phi = twopi*G4UniformRand() ;
531         dirY = std::sin(theta) * std::sin(phi) << 693         dirX = sin(theta)*cos(phi)  ;
532         dirZ = std::cos(theta);                << 694         dirY = sin(theta)*sin(phi)  ;
533         G4ThreeVector directionTR(dirX, dirY,  << 695         dirZ = cos(theta)           ;
534         directionTR.rotateUz(particleDir);     << 696         G4ThreeVector directionTR(dirX,dirY,dirZ) ;
535         auto aPhotonTR = new G4DynamicParticle << 697         directionTR.rotateUz(particleDir) ;
536                                                << 698         G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
537   // Create the G4Track                        << 699                                                              directionTR,
538   auto aSecondaryTrack = new G4Track(aPhotonTR << 700                                                              energyTR     ) ;
539   aSecondaryTrack->SetTouchableHandle(aStep.Ge << 701         aParticleChange.AddSecondary(aPhotonTR) ;
540   aSecondaryTrack->SetParentID(aTrack.GetTrack << 
541   aSecondaryTrack->SetCreatorModelID(secID);   << 
542   aParticleChange.AddSecondary(aSecondaryTrack << 
543       }                                           702       }
544     }                                             703     }
545   }                                               704   }
546   else                                            705   else
547   {                                               706   {
548     if(iTkin == 0)  // Tkin is too small, negl << 707     if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
549     {                                             708     {
550       return G4VDiscreteProcess::PostStepDoIt(    709       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
551     }                                          << 710     } 
552     else  // general case: Tkin between two ve << 711     else          // general case: Tkin between two vectors of the material
553     {                                             712     {
554       E1 = fProtonEnergyVector->GetLowEdgeEner << 713       E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 
555       E2 = fProtonEnergyVector->GetLowEdgeEner << 714       E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin)     ;
556       W  = 1.0 / (E2 - E1);                    << 715        W = 1.0/(E2 - E1) ;
557       W1 = (E2 - TkinScaled) * W;              << 716       W1 = (E2 - TkinScaled)*W ;
558       W2 = (TkinScaled - E1) * W;              << 717       W2 = (TkinScaled - E1)*W ;
559                                                << 718 
560       numOfTR = (G4int)G4Poisson((((*(*fEnergy << 719   // G4cout<<iTkin<<" mean TR number = "<<(((*(*fEnergyDistrTable)(iPlace))(0)+
561                                    (*(*fAngleD << 720   // (*(*fAngleDistrTable)(iPlace))(0))*W1 + 
562                                     W1 +       << 721   //                                ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
563                                   ((*(*fEnergy << 722   // (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
564                                    (*(*fAngleD << 723   //                                    *chargeSq*0.5<<G4endl ;
565                                     W2) *      << 724 
566                                  chargeSq * 0. << 725       numOfTR = RandPoisson::shoot((((*(*fEnergyDistrTable)(iPlace))(0)+
                                                   >> 726 (*(*fAngleDistrTable)(iPlace))(0))*W1 + 
                                                   >> 727                                     ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
                                                   >> 728 (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
                                                   >> 729                                     *chargeSq*0.5 ) ;
567       if(numOfTR == 0)                            730       if(numOfTR == 0)
568       {                                           731       {
569         return G4VDiscreteProcess::PostStepDoI    732         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
570       }                                           733       }
571       else                                        734       else
572       {                                           735       {
                                                   >> 736         // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl ;
                                                   >> 737 
573         aParticleChange.SetNumberOfSecondaries    738         aParticleChange.SetNumberOfSecondaries(numOfTR);
574         for(iTR = 0; iTR < numOfTR; ++iTR)     << 739         for(iTR=0;iTR<numOfTR;iTR++)
575         {                                         740         {
576           energyPos = ((*(*fEnergyDistrTable)( << 741           energyPos = ((*(*fEnergyDistrTable)(iPlace))(0)*W1+
577                        (*(*fEnergyDistrTable)( << 742                        (*(*fEnergyDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand() ;
578                       G4UniformRand();         << 743           for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
579           for(iTransfer = 0; iTransfer < fBinT << 744       {
580           {                                    << 745             if(energyPos >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1+
581             if(energyPos >=                    << 746                        (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer)*W2)) break ;
582                ((*(*fEnergyDistrTable)(iPlace) << 747       }
583                 (*(*fEnergyDistrTable)(iPlace  << 748           energyTR = ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer))*W1+
584               break;                           << 749                ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer))*W2 ;
585           }                                    << 750 
586           energyTR =                           << 751     // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl ;
587             ((*fEnergyDistrTable)(iPlace)->Get << 752 
588             ((*fEnergyDistrTable)(iPlace + 1)- << 753           kinEnergy -= energyTR ; 
589               W2;                              << 754           aParticleChange.SetEnergyChange(kinEnergy);
590                                                << 755 
591           kinEnergy -= energyTR;               << 756           anglePos = ((*(*fAngleDistrTable)(iPlace))(0)*W1+
592           aParticleChange.ProposeEnergy(kinEne << 757                        (*(*fAngleDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand() ;
593                                                << 758           for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
594           anglePos = ((*(*fAngleDistrTable)(iP << 759     {
595                       (*(*fAngleDistrTable)(iP << 760             if(anglePos > ((*(*fAngleDistrTable)(iPlace))(iTransfer)*W1+
596                      G4UniformRand();          << 761                       (*(*fAngleDistrTable)(iPlace + 1))(iTransfer)*W2)) break ;
597           for(iTransfer = 0; iTransfer < fBinT << 762     }
598           {                                    << 763           theta = sqrt(((*fAngleDistrTable)(iPlace)->
599             if(anglePos > ((*(*fAngleDistrTabl << 764                         GetLowEdgeEnergy(iTransfer-1))*W1+
600                            (*(*fAngleDistrTabl << 765                   ((*fAngleDistrTable)(iPlace + 1)->
601               break;                           << 766                         GetLowEdgeEnergy(iTransfer-1))*W2) ;
602           }                                    << 767 
603           theta = std::sqrt(                   << 768     // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl ;
604             ((*fAngleDistrTable)(iPlace)->GetL << 769 
605               W1 +                             << 770           phi = twopi*G4UniformRand() ;
606             ((*fAngleDistrTable)(iPlace + 1)-> << 771           dirX = sin(theta)*cos(phi)  ;
607               W2);                             << 772           dirY = sin(theta)*sin(phi)  ;
608                                                << 773           dirZ = cos(theta)           ;
609           phi  = twopi * G4UniformRand();      << 774           G4ThreeVector directionTR(dirX,dirY,dirZ) ;
610           dirX = std::sin(theta) * std::cos(ph << 775           directionTR.rotateUz(particleDir) ;
611           dirY = std::sin(theta) * std::sin(ph << 776           G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
612           dirZ = std::cos(theta);              << 777                                                                directionTR,
613           G4ThreeVector directionTR(dirX, dirY << 778                                                                energyTR     ) ;
614           directionTR.rotateUz(particleDir);   << 779           aParticleChange.AddSecondary(aPhotonTR) ;
615           auto aPhotonTR =                     << 
616             new G4DynamicParticle(G4Gamma::Gam << 
617                                                << 
618     // Create the G4Track                      << 
619     G4Track* aSecondaryTrack = new G4Track(aPh << 
620     aSecondaryTrack->SetTouchableHandle(aStep. << 
621     aSecondaryTrack->SetParentID(aTrack.GetTra << 
622     aSecondaryTrack->SetCreatorModelID(secID); << 
623     aParticleChange.AddSecondary(aSecondaryTra << 
624         }                                         780         }
625       }                                           781       }
626     }                                             782     }
627   }                                               783   }
628   return &aParticleChange;                     << 784   return &aParticleChange ;
629 }                                                 785 }
630                                                   786 
631 //////////////////////////////////////////////    787 ////////////////////////////////////////////////////////////////////////////
                                                   >> 788 //
632 // Test function for checking of PostStepDoIt     789 // Test function for checking of PostStepDoIt random preparation of TR photon
633 // energy                                         790 // energy
634 G4double G4ForwardXrayTR::GetEnergyTR(G4int iM << 791 //
                                                   >> 792 
                                                   >> 793 G4double 
                                                   >> 794 G4ForwardXrayTR::GetEnergyTR(G4int iMat, G4int jMat, G4int iTkin) const
635 {                                                 795 {
636   G4int iPlace, numOfTR, iTR, iTransfer;       << 796   G4int  iPlace, numOfMat, numOfTR, iTR, iTransfer ;
637   G4double energyTR = 0.0;  // return this val << 797   G4double energyTR = 0.0 ; // return this value for no TR photons
638   G4double energyPos;                          << 798   G4double energyPos  ;
639   G4double W1, W2;                             << 799   G4double  W1, W2;
640                                                << 800 
641   const G4ProductionCutsTable* theCoupleTable  << 801   static 
642     G4ProductionCutsTable::GetProductionCutsTa << 802   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
643   G4int numOfCouples = (G4int)theCoupleTable-> << 803   numOfMat = G4Material::GetNumberOfMaterials() ;
                                                   >> 804 
644                                                   805 
645   // The case of equal or approximate (in term    806   // The case of equal or approximate (in terms of plasma energy) materials
646   // No TR photons ?!                             807   // No TR photons ?!
647                                                   808 
648   const G4MaterialCutsCouple* iCouple =        << 
649     theCoupleTable->GetMaterialCutsCouple(iMat << 
650   const G4MaterialCutsCouple* jCouple =        << 
651     theCoupleTable->GetMaterialCutsCouple(jMat << 
652   const G4Material* iMaterial = iCouple->GetMa << 
653   const G4Material* jMaterial = jCouple->GetMa << 
654                                                   809 
655   if(iMat == jMat                              << 810   if (     iMat == jMat
656                                                << 
657      || iMaterial->GetState() == jMaterial->Ge << 
658                                                << 
659      || (iMaterial->GetState() == kStateSolid  << 
660          jMaterial->GetState() == kStateLiquid << 
661                                                << 
662      || (iMaterial->GetState() == kStateLiquid << 
663          jMaterial->GetState() == kStateSolid) << 
664                                                   811 
                                                   >> 812       || (*theMaterialTable)[iMat]->GetState() ==
                                                   >> 813          (*theMaterialTable)[jMat]->GetState()
                                                   >> 814   
                                                   >> 815       ||(    (*theMaterialTable)[iMat]->GetState() == kStateSolid
                                                   >> 816           && (*theMaterialTable)[jMat]->GetState() == kStateLiquid )
                                                   >> 817  
                                                   >> 818       ||(    (*theMaterialTable)[iMat]->GetState() == kStateLiquid
                                                   >> 819           && (*theMaterialTable)[jMat]->GetState() == kStateSolid  )   )
665   {                                               820   {
666     return energyTR;                           << 821     return energyTR ;
667   }                                               822   }
668                                                   823 
669   if(jMat < iMat)                                 824   if(jMat < iMat)
670   {                                               825   {
671     iPlace = (iMat * (numOfCouples - 1) + jMat << 826     iPlace = (iMat*(numOfMat - 1) + jMat)*fTotBin + iTkin - 1 ;
672   }                                               827   }
673   else                                            828   else
674   {                                               829   {
675     iPlace = (iMat * (numOfCouples - 1) + jMat << 830     iPlace = (iMat*(numOfMat - 1) + jMat - 1)*fTotBin + iTkin - 1 ;
676   }                                               831   }
677   G4PhysicsVector* energyVector1 = (*fEnergyDi << 832   G4PhysicsVector*  energyVector1 = (*fEnergyDistrTable)(iPlace)     ;
678   G4PhysicsVector* energyVector2 = (*fEnergyDi << 833   G4PhysicsVector*  energyVector2 = (*fEnergyDistrTable)(iPlace + 1) ;
679                                                   834 
680   if(iTkin == fTotBin)  // TR plato, try from  << 835   if(iTkin == fTotBin)                 // TR plato, try from left
681   {                                               836   {
682     numOfTR = (G4int)G4Poisson((*energyVector1 << 837     numOfTR = RandPoisson::shoot( (*energyVector1)(0)  ) ;
683     if(numOfTR == 0)                              838     if(numOfTR == 0)
684     {                                             839     {
685       return energyTR;                         << 840       return energyTR ;
686     }                                             841     }
687     else                                          842     else
688     {                                             843     {
689       for(iTR = 0; iTR < numOfTR; ++iTR)       << 844       for(iTR=0;iTR<numOfTR;iTR++)
690       {                                           845       {
691         energyPos = (*energyVector1)(0) * G4Un << 846         energyPos = (*energyVector1)(0)*G4UniformRand() ;
692         for(iTransfer = 0; iTransfer < fBinTR  << 847         for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
693         {                                      << 848   {
694           if(energyPos >= (*energyVector1)(iTr << 849           if(energyPos >= (*energyVector1)(iTransfer)) break ;
695             break;                             << 850   }
696         }                                      << 851         energyTR += energyVector1->GetLowEdgeEnergy(iTransfer) ;
697         energyTR += energyVector1->GetLowEdgeE << 
698       }                                           852       }
699     }                                             853     }
700   }                                               854   }
701   else                                            855   else
702   {                                               856   {
703     if(iTkin == 0)  // Tkin is too small, negl << 857     if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
704     {                                             858     {
705       return energyTR;                         << 859       return energyTR ;
706     }                                          << 860     } 
707     else  // general case: Tkin between two ve << 861     else          // general case: Tkin between two vectors of the material
708     {     // use trivial mean half/half        << 862     {             // use trivial mean half/half
709       W1      = 0.5;                           << 863       W1 = 0.5 ; 
710       W2      = 0.5;                           << 864       W2 = 0.5 ;
711       numOfTR = (G4int)G4Poisson((*energyVecto << 865      numOfTR = RandPoisson::shoot( (*energyVector1)(0)*W1 +
                                                   >> 866                                    (*energyVector2)(0)*W2  ) ;
712       if(numOfTR == 0)                            867       if(numOfTR == 0)
713       {                                           868       {
714         return energyTR;                       << 869         return energyTR ;
715       }                                           870       }
716       else                                        871       else
717       {                                           872       {
718         G4cout << "It is still OK in GetEnergy << 873   G4cout<<"It is still OK in GetEnergyTR(int,int,int)"<<G4endl;
719         for(iTR = 0; iTR < numOfTR; ++iTR)     << 874         for(iTR=0;iTR<numOfTR;iTR++)
720         {                                         875         {
721           energyPos = ((*energyVector1)(0) * W << 876           energyPos = ((*energyVector1)(0)*W1+
722                       G4UniformRand();         << 877                        (*energyVector2)(0)*W2)*G4UniformRand() ;
723           for(iTransfer = 0; iTransfer < fBinT << 878           for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
724           {                                    << 879       {
725             if(energyPos >= ((*energyVector1)( << 880             if(energyPos >= ((*energyVector1)(iTransfer)*W1+
726                              (*energyVector2)( << 881                              (*energyVector2)(iTransfer)*W2)) break ;
727               break;                           << 882       }
728           }                                    << 883           energyTR += (energyVector1->GetLowEdgeEnergy(iTransfer))*W1+
729           energyTR += (energyVector1->GetLowEd << 884                       (energyVector2->GetLowEdgeEnergy(iTransfer))*W2 ;
730                       (energyVector2->GetLowEd << 885 
731         }                                         886         }
732       }                                           887       }
733     }                                             888     }
734   }                                               889   }
735                                                   890 
736   return energyTR;                             << 891   return energyTR   ;
737 }                                                 892 }
738                                                   893 
739 //////////////////////////////////////////////    894 ////////////////////////////////////////////////////////////////////////////
                                                   >> 895 //
740 // Test function for checking of PostStepDoIt     896 // Test function for checking of PostStepDoIt random preparation of TR photon
741 // theta angle relative to particle direction     897 // theta angle relative to particle direction
742 G4double G4ForwardXrayTR::GetThetaTR(G4int, G4 << 898 //
743                                                << 899  
744 G4int G4ForwardXrayTR::GetSympsonNumber() { re << 900        
                                                   >> 901 G4double 
                                                   >> 902 G4ForwardXrayTR::GetThetaTR(G4int iMat, G4int jMat, G4int iTkin) const     
                                                   >> 903 {
                                                   >> 904   G4double theta = 0.0 ;
745                                                   905 
746 G4int G4ForwardXrayTR::GetBinTR() { return fBi << 906   return theta   ;
                                                   >> 907 }
747                                                   908 
748 G4double G4ForwardXrayTR::GetMinProtonTkin() { << 
749                                                   909 
750 G4double G4ForwardXrayTR::GetMaxProtonTkin() { << 
751                                                   910 
752 G4int G4ForwardXrayTR::GetTotBin() { return fT << 911 // end of G4ForwardXrayTR implementation file 
                                                   >> 912 //
                                                   >> 913 ///////////////////////////////////////////////////////////////////////////
753                                                   914