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