Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.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/adjoint/src/G4AdjointBremsstrahlungModel.cc (Version 11.3.0) and /processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc (Version 9.5.p2)


  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                                                <<  26 // $Id: G4AdjointBremsstrahlungModel.cc,v 1.6 2010-11-11 11:51:56 ldesorgh Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
                                                   >>  28 //
 27 #include "G4AdjointBremsstrahlungModel.hh"         29 #include "G4AdjointBremsstrahlungModel.hh"
 28                                                << 
 29 #include "G4AdjointCSManager.hh"                   30 #include "G4AdjointCSManager.hh"
                                                   >>  31 #include "G4Integrator.hh"
                                                   >>  32 #include "G4TrackStatus.hh"
                                                   >>  33 #include "G4ParticleChange.hh"
 30 #include "G4AdjointElectron.hh"                    34 #include "G4AdjointElectron.hh"
 31 #include "G4AdjointGamma.hh"                       35 #include "G4AdjointGamma.hh"
 32 #include "G4Electron.hh"                           36 #include "G4Electron.hh"
 33 #include "G4EmModelManager.hh"                 << 
 34 #include "G4Gamma.hh"                          << 
 35 #include "G4ParticleChange.hh"                 << 
 36 #include "G4PhysicalConstants.hh"              << 
 37 #include "G4SeltzerBergerModel.hh"             << 
 38 #include "G4SystemOfUnits.hh"                  << 
 39 #include "G4TrackStatus.hh"                    << 
 40                                                    37 
 41 ////////////////////////////////////////////// <<  38 #include "G4Timer.hh"
 42 G4AdjointBremsstrahlungModel::G4AdjointBremsst <<  39 //#include "G4PenelopeBremsstrahlungModel.hh"
 43   : G4VEmAdjointModel("AdjointeBremModel")     << 
 44 {                                              << 
 45   fDirectModel = aModel;                       << 
 46   Initialize();                                << 
 47 }                                              << 
 48                                                    40 
 49 ////////////////////////////////////////////// << 
 50 G4AdjointBremsstrahlungModel::G4AdjointBremsst << 
 51   : G4VEmAdjointModel("AdjointeBremModel")     << 
 52 {                                              << 
 53   fDirectModel = new G4SeltzerBergerModel();   << 
 54   Initialize();                                << 
 55 }                                              << 
 56                                                    41 
 57 //////////////////////////////////////////////     42 ////////////////////////////////////////////////////////////////////////////////
 58 void G4AdjointBremsstrahlungModel::Initialize( <<  43 //
 59 {                                              <<  44 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
                                                   >>  45  G4VEmAdjointModel("AdjointeBremModel"),
                                                   >>  46   MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi)
                                                   >>  47 { 
 60   SetUseMatrix(false);                             48   SetUseMatrix(false);
 61   SetUseMatrixPerElement(false);                   49   SetUseMatrixPerElement(false);
 62                                                <<  50   
 63   fEmModelManagerForFwdModels = new G4EmModelM <<  51   theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel");
 64   fEmModelManagerForFwdModels->AddEmModel(1, f <<  52   theDirectEMModel=theDirectStdBremModel;
                                                   >>  53   theEmModelManagerForFwdModels = new G4EmModelManager();
                                                   >>  54   isDirectModelInitialised = false;
                                                   >>  55   G4VEmFluctuationModel* f=0;
                                                   >>  56   G4Region* r=0;
                                                   >>  57   theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r);
                                                   >>  58  // theDirectPenelopeBremModel =0;
                                                   >>  59   
 65   SetApplyCutInRange(true);                        60   SetApplyCutInRange(true);
                                                   >>  61   highKinEnergy= 100.*TeV;
                                                   >>  62   lowKinEnergy = 1.0*keV;
 66                                                    63 
 67   fElectron = G4Electron::Electron();          <<  64   lastCZ =0.;
 68   fGamma    = G4Gamma::Gamma();                << 
 69                                                    65 
 70   fAdjEquivDirectPrimPart   = G4AdjointElectro <<  66   
 71   fAdjEquivDirectSecondPart = G4AdjointGamma:: <<  67   theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
 72   fDirectPrimaryPart        = fElectron;       <<  68   theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
 73   fSecondPartSameType       = false;           <<  69   theDirectPrimaryPartDef=G4Electron::Electron();
                                                   >>  70   second_part_of_same_type=false;
                                                   >>  71   
                                                   >>  72   /*UsePenelopeModel=false;
                                                   >>  73   if (UsePenelopeModel) {
                                                   >>  74     G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem");
                                                   >>  75   theEmModelManagerForFwdModels = new G4EmModelManager();
                                                   >>  76     isPenelopeModelInitialised = false;
                                                   >>  77   G4VEmFluctuationModel* f=0;
                                                   >>  78   G4Region* r=0;
                                                   >>  79   theDirectEMModel=thePenelopeModel;
                                                   >>  80   theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r);
                                                   >>  81   }
                                                   >>  82   */  
                                                   >>  83   
 74                                                    84 
 75   fCSManager = G4AdjointCSManager::GetAdjointC <<  85   
 76 }                                                  86 }
 77                                                << 
 78 //////////////////////////////////////////////     87 ////////////////////////////////////////////////////////////////////////////////
                                                   >>  88 //
 79 G4AdjointBremsstrahlungModel::~G4AdjointBremss     89 G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
 80 {                                              <<  90 {if (theDirectStdBremModel) delete theDirectStdBremModel;
 81   if(fEmModelManagerForFwdModels)              <<  91  if (theEmModelManagerForFwdModels) delete theEmModelManagerForFwdModels;
 82     delete fEmModelManagerForFwdModels;        << 
 83 }                                                  92 }
 84                                                    93 
 85 //////////////////////////////////////////////     94 ////////////////////////////////////////////////////////////////////////////////
 86 void G4AdjointBremsstrahlungModel::SampleSecon <<  95 //
 87   const G4Track& aTrack, G4bool isScatProjToPr <<  96 void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
 88   G4ParticleChange* fParticleChange)           <<  97                        G4bool IsScatProjToProjCase,
                                                   >>  98                  G4ParticleChange* fParticleChange)
 89 {                                                  99 {
 90   if(!fUseMatrix)                              << 100  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 
 91     return RapidSampleSecondaries(aTrack, isSc << 
 92                                                << 
 93   const G4DynamicParticle* theAdjointPrimary = << 
 94   DefineCurrentMaterial(aTrack.GetMaterialCuts << 
 95                                                << 
 96   G4double adjointPrimKinEnergy   = theAdjoint << 
 97   G4double adjointPrimTotalEnergy = theAdjoint << 
 98                                                   101 
 99   if(adjointPrimKinEnergy > GetHighEnergyLimit << 102  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
100   {                                            << 103  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
101     return;                                    << 104  
102   }                                            << 105  
                                                   >> 106  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
                                                   >> 107  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
                                                   >> 108  
                                                   >> 109  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
                                                   >> 110   return;
                                                   >> 111  }
                                                   >> 112   
                                                   >> 113   G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
                                                   >> 114                 IsScatProjToProjCase);
                                                   >> 115  //Weight correction
                                                   >> 116  //-----------------------             
                                                   >> 117  CorrectPostStepWeight(fParticleChange, 
                                                   >> 118            aTrack.GetWeight(), 
                                                   >> 119            adjointPrimKinEnergy,
                                                   >> 120            projectileKinEnergy,
                                                   >> 121            IsScatProjToProjCase); 
                                                   >> 122  
                                                   >> 123  
                                                   >> 124  //Kinematic
                                                   >> 125  //---------
                                                   >> 126  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
                                                   >> 127  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
                                                   >> 128  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 
                                                   >> 129  G4double projectileP = std::sqrt(projectileP2);
                                                   >> 130  
                                                   >> 131  
                                                   >> 132  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
                                                   >> 133  //------------------------------------------------
                                                   >> 134   G4double u;
                                                   >> 135   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
103                                                   136 
104   G4double projectileKinEnergy =               << 137   if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
105     SampleAdjSecEnergyFromCSMatrix(adjointPrim << 138      else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
106                                                   139 
107   // Weight correction                         << 140   G4double theta = u*electron_mass_c2/projectileTotalEnergy;
108   CorrectPostStepWeight(fParticleChange, aTrac << 
109                         adjointPrimKinEnergy,  << 
110                         isScatProjToProj);     << 
111                                                << 
112   // Kinematic                                 << 
113   G4double projectileM0          = fAdjEquivDi << 
114   G4double projectileTotalEnergy = projectileM << 
115   G4double projectileP2 =                      << 
116     projectileTotalEnergy * projectileTotalEne << 
117   G4double projectileP = std::sqrt(projectileP << 
118                                                   141 
119   // Angle of the gamma direction with the pro << 142   G4double sint = std::sin(theta);
120   // G4eBremsstrahlungModel                    << 143   G4double cost = std::cos(theta);
121   G4double u;                                  << 144 
122   if(0.25 > G4UniformRand())                   << 145   G4double phi = twopi * G4UniformRand() ;
123     u = -std::log(G4UniformRand() * G4UniformR << 146   
124   else                                         << 147   G4ThreeVector projectileMomentum;
125     u = -std::log(G4UniformRand() * G4UniformR << 148   projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
126                                                << 149   if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
127   G4double theta = u * electron_mass_c2 / proj << 150     G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
128   G4double sint  = std::sin(theta);            << 151   G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
129   G4double cost  = std::cos(theta);            << 152   G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
130                                                << 153   G4double sint1 =  std::sqrt(1.-cost1*cost1);
131   G4double phi = twopi * G4UniformRand();      << 154   projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
132                                                << 155   
133   G4ThreeVector projectileMomentum =           << 
134     G4ThreeVector(std::cos(phi) * sint, std::s << 
135     projectileP;  // gamma frame               << 
136   if(isScatProjToProj)                         << 
137   {  // the adjoint primary is the scattered e << 
138     G4ThreeVector gammaMomentum =              << 
139       (projectileTotalEnergy - adjointPrimTota << 
140       G4ThreeVector(0., 0., 1.);               << 
141     G4ThreeVector dirProd = projectileMomentum << 
142     G4double cost1        = std::cos(dirProd.a << 
143     G4double sint1        = std::sqrt(1. - cos << 
144     projectileMomentum =                       << 
145       G4ThreeVector(std::cos(phi) * sint1, std << 
146       projectileP;                             << 
147   }                                               156   }
148                                                << 157   
149   projectileMomentum.rotateUz(theAdjointPrimar    158   projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
150                                                << 159  
151   if(!isScatProjToProj)                        << 160  
152   {  // kill the primary and add a secondary   << 161  
153     fParticleChange->ProposeTrackStatus(fStopA << 162   if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
154     fParticleChange->AddSecondary(             << 163   fParticleChange->ProposeTrackStatus(fStopAndKill);
155       new G4DynamicParticle(fAdjEquivDirectPri << 164   fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
156   }                                            << 165   }
157   else                                         << 166   else {
158   {                                            << 167   fParticleChange->ProposeEnergy(projectileKinEnergy);
159     fParticleChange->ProposeEnergy(projectileK << 168   fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
160     fParticleChange->ProposeMomentumDirection( << 169   
161   }                                            << 170   } 
162 }                                              << 171 } 
163                                                << 
164 //////////////////////////////////////////////    172 ////////////////////////////////////////////////////////////////////////////////
165 void G4AdjointBremsstrahlungModel::RapidSample << 173 //
166   const G4Track& aTrack, G4bool isScatProjToPr << 174 void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack,
167   G4ParticleChange* fParticleChange)           << 175                        G4bool IsScatProjToProjCase,
168 {                                              << 176                  G4ParticleChange* fParticleChange)
169   const G4DynamicParticle* theAdjointPrimary = << 177 { 
170   DefineCurrentMaterial(aTrack.GetMaterialCuts << 178 
171                                                << 179  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
172   G4double adjointPrimKinEnergy   = theAdjoint << 180  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
173   G4double adjointPrimTotalEnergy = theAdjoint << 181  
174                                                << 182  
175   if(adjointPrimKinEnergy > GetHighEnergyLimit << 183  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
176   {                                            << 184  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
177     return;                                    << 185  
178   }                                            << 186  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
                                                   >> 187   return;
                                                   >> 188  }
                                                   >> 189   
                                                   >> 190  G4double projectileKinEnergy =0.;
                                                   >> 191  G4double gammaEnergy=0.;
                                                   >> 192  G4double diffCSUsed=0.; 
                                                   >> 193  if (!IsScatProjToProjCase){
                                                   >> 194   gammaEnergy=adjointPrimKinEnergy;
                                                   >> 195   G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
                                                   >> 196         G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
                                                   >> 197   if (Emin>=Emax) return;
                                                   >> 198   projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
                                                   >> 199   diffCSUsed=lastCZ/projectileKinEnergy;
                                                   >> 200   
                                                   >> 201  }
                                                   >> 202  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
                                                   >> 203   G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
                                                   >> 204   if (Emin>=Emax) return;
                                                   >> 205   G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
                                                   >> 206   G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
                                                   >> 207   //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl;
                                                   >> 208   projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
                                                   >> 209   gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
                                                   >> 210   diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
                                                   >> 211   
                                                   >> 212  }
                                                   >> 213   
                                                   >> 214   
                                                   >> 215   
                                                   >> 216                 
                                                   >> 217  //Weight correction
                                                   >> 218  //-----------------------
                                                   >> 219  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
                                                   >> 220  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
                                                   >> 221 
                                                   >> 222  //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model
                                                   >> 223  //Here we consider the true  diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term.
                                                   >> 224  //Basically any other differential CS   diffCS could be used here (example Penelope). 
                                                   >> 225  
                                                   >> 226  G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
                                                   >> 227  w_corr*=diffCS/diffCSUsed;
                                                   >> 228      
                                                   >> 229  G4double new_weight = aTrack.GetWeight()*w_corr;
                                                   >> 230  fParticleChange->SetParentWeightByProcess(false);
                                                   >> 231  fParticleChange->SetSecondaryWeightByProcess(false);
                                                   >> 232  fParticleChange->ProposeParentWeight(new_weight);
                                                   >> 233  
                                                   >> 234  //Kinematic
                                                   >> 235  //---------
                                                   >> 236  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
                                                   >> 237  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
                                                   >> 238  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 
                                                   >> 239  G4double projectileP = std::sqrt(projectileP2);
                                                   >> 240  
                                                   >> 241  
                                                   >> 242  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
                                                   >> 243  //------------------------------------------------
                                                   >> 244   G4double u;
                                                   >> 245   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
179                                                   246 
180   G4double projectileKinEnergy = 0.;           << 247   if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
181   G4double gammaEnergy         = 0.;           << 248      else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
182   G4double diffCSUsed          = 0.;           << 
183   if(!isScatProjToProj)                        << 
184   {                                            << 
185     gammaEnergy   = adjointPrimKinEnergy;      << 
186     G4double Emax = GetSecondAdjEnergyMaxForPr << 
187     G4double Emin = GetSecondAdjEnergyMinForPr << 
188     if(Emin >= Emax)                           << 
189       return;                                  << 
190     projectileKinEnergy = Emin * std::pow(Emax << 
191     diffCSUsed          = fCsBiasingFactor * f << 
192   }                                            << 
193   else                                         << 
194   {                                            << 
195     G4double Emax =                            << 
196       GetSecondAdjEnergyMaxForScatProjToProj(a << 
197     G4double Emin =                            << 
198       GetSecondAdjEnergyMinForScatProjToProj(a << 
199     if(Emin >= Emax)                           << 
200       return;                                  << 
201     G4double f1 = (Emin - adjointPrimKinEnergy << 
202     G4double f2 = (Emax - adjointPrimKinEnergy << 
203     projectileKinEnergy =                      << 
204       adjointPrimKinEnergy / (1. - f1 * std::p << 
205     gammaEnergy = projectileKinEnergy - adjoin << 
206     diffCSUsed =                               << 
207       fLastCZ * adjointPrimKinEnergy / project << 
208   }                                            << 
209                                                   249 
210   // Weight correction:                        << 250   G4double theta = u*electron_mass_c2/projectileTotalEnergy;
211   // First w_corr is set to the ratio between  << 
212   // if this has to be done in the model.      << 
213   // For the case of forced interaction this w << 
214   // the forced interaction.  It is important  << 
215   // creation of the secondary                 << 
216   G4double w_corr = fOutsideWeightFactor;      << 
217   if(fInModelWeightCorr)                       << 
218   {                                            << 
219     w_corr = fCSManager->GetPostStepWeightCorr << 
220   }                                            << 
221                                                   251 
222   // Then another correction is needed due to  << 252   G4double sint = std::sin(theta);
223   // differential CS has been used rather than << 253   G4double cost = std::cos(theta);
224   // direct model Here we consider the true di << 254 
225   // numerical differentiation over Tcut of th << 255   G4double phi = twopi * G4UniformRand() ;
226   // Migdal term. Basically any other differen << 256   
227   // (example Penelope).                       << 257   G4ThreeVector projectileMomentum;
228   G4double diffCS = DiffCrossSectionPerVolumeP << 258   projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
229     fCurrentMaterial, projectileKinEnergy, gam << 259   if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
230   w_corr *= diffCS / diffCSUsed;               << 260     G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
231                                                << 261   G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
232   G4double new_weight = aTrack.GetWeight() * w << 262   G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
233   fParticleChange->SetParentWeightByProcess(fa << 263   G4double sint1 =  std::sqrt(1.-cost1*cost1);
234   fParticleChange->SetSecondaryWeightByProcess << 264   projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
235   fParticleChange->ProposeParentWeight(new_wei << 265   
236                                                << 
237   // Kinematic                                 << 
238   G4double projectileM0          = fAdjEquivDi << 
239   G4double projectileTotalEnergy = projectileM << 
240   G4double projectileP2 =                      << 
241     projectileTotalEnergy * projectileTotalEne << 
242   G4double projectileP = std::sqrt(projectileP << 
243                                                << 
244   // Use the angular model of the forward mode << 
245   // Dummy dynamic particle to use the model   << 
246   G4DynamicParticle* aDynPart =                << 
247     new G4DynamicParticle(fElectron, G4ThreeVe << 
248                                                << 
249   // Get the element from the direct model     << 
250   const G4Element* elm = fDirectModel->SelectR << 
251     fCurrentCouple, fElectron, projectileKinEn << 
252   G4int Z         = elm->GetZasInt();          << 
253   G4double energy = aDynPart->GetTotalEnergy() << 
254   G4ThreeVector projectileMomentum =           << 
255     fDirectModel->GetAngularDistribution()->Sa << 
256                                             fC << 
257   G4double phi = projectileMomentum.getPhi();  << 
258                                                << 
259   if(isScatProjToProj)                         << 
260   {  // the adjoint primary is the scattered e << 
261     G4ThreeVector gammaMomentum =              << 
262       (projectileTotalEnergy - adjointPrimTota << 
263       G4ThreeVector(0., 0., 1.);               << 
264     G4ThreeVector dirProd = projectileMomentum << 
265     G4double cost1        = std::cos(dirProd.a << 
266     G4double sint1        = std::sqrt(1. - cos << 
267     projectileMomentum =                       << 
268       G4ThreeVector(std::cos(phi) * sint1, std << 
269       projectileP;                             << 
270   }                                               266   }
271                                                << 267   
272   projectileMomentum.rotateUz(theAdjointPrimar    268   projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
273                                                << 269  
274   if(!isScatProjToProj)                        << 270  
275   {  // kill the primary and add a secondary   << 271  
276     fParticleChange->ProposeTrackStatus(fStopA << 272   if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
277     fParticleChange->AddSecondary(             << 273   fParticleChange->ProposeTrackStatus(fStopAndKill);
278       new G4DynamicParticle(fAdjEquivDirectPri << 274   fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
279   }                                            << 275   }
280   else                                         << 276   else {
281   {                                            << 277   fParticleChange->ProposeEnergy(projectileKinEnergy);
282     fParticleChange->ProposeEnergy(projectileK << 278   fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
283     fParticleChange->ProposeMomentumDirection( << 279   
284   }                                            << 280   } 
285 }                                              << 281 } 
                                                   >> 282 ////////////////////////////////////////////////////////////////////////////////
                                                   >> 283 //
                                                   >> 284 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
                                                   >> 285                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
                                                   >> 286                                       G4double kinEnergyProd // kinetic energy of the secondary particle 
                                                   >> 287               )
                                                   >> 288 {if (!isDirectModelInitialised) {
                                                   >> 289   theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
                                                   >> 290   isDirectModelInitialised =true;
                                                   >> 291  }
                                                   >> 292 
                                                   >> 293  return  DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
                                                   >> 294                                                      kinEnergyProj, 
                                                   >> 295                                                      kinEnergyProd);
                                                   >> 296  /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
                                                   >> 297                                                      kinEnergyProj, 
                                                   >> 298                                                      kinEnergyProd);*/                                 
                                                   >> 299 }             
286                                                   300 
287 //////////////////////////////////////////////    301 ////////////////////////////////////////////////////////////////////////////////
288 G4double G4AdjointBremsstrahlungModel::DiffCro << 302 //
289   const G4Material* aMaterial,                 << 303 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1(
290   G4double kinEnergyProj,  // kin energy of pr << 304                 const G4Material* aMaterial,
291   G4double kinEnergyProd   // kinetic energy o << 305                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
292 )                                              << 306                                       G4double kinEnergyProd // kinetic energy of the secondary particle 
                                                   >> 307               )
293 {                                                 308 {
294   if(!fIsDirectModelInitialised)               << 309  G4double dCrossEprod=0.;
295   {                                            << 310  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
296     fEmModelManagerForFwdModels->Initialise(fE << 311  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
297     fIsDirectModelInitialised = true;          << 312  
298   }                                            << 313  
299   return G4VEmAdjointModel::DiffCrossSectionPe << 314  //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
300     aMaterial, kinEnergyProj, kinEnergyProd);  << 315  //This is what is applied in the discrete standard model before the  rejection test  that make a correction
                                                   >> 316  //The application of the same rejection function is not possible here.
                                                   >> 317  //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the 
                                                   >> 318  // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and 
                                                   >> 319  // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. 
                                                   >> 320  // In the future we plan to use the brem secondary spectra from the G4Penelope implementation 
                                                   >> 321  
                                                   >> 322  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
                                                   >> 323   G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
                                                   >> 324   dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
                                                   >> 325  }
                                                   >> 326  return dCrossEprod;
                                                   >> 327   
301 }                                                 328 }
302                                                   329 
303 //////////////////////////////////////////////    330 ////////////////////////////////////////////////////////////////////////////////
304 G4double G4AdjointBremsstrahlungModel::Adjoint << 331 //
305   const G4MaterialCutsCouple* aCouple, G4doubl << 332 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2(
306   G4bool isScatProjToProj)                     << 333                 const G4Material* material,
                                                   >> 334                                       G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction 
                                                   >> 335                                       G4double kinEnergyProd // kinetic energy of the secondary particle 
                                                   >> 336               )
307 {                                                 337 {
308   static constexpr G4double maxEnergy = 100. * << 338  //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor 
309   // 2.78.. == std::exp(1.)                    << 339   //used in the direct model
310   if(!fIsDirectModelInitialised)               << 340  
311   {                                            << 341  G4double dCrossEprod=0.;
312     fEmModelManagerForFwdModels->Initialise(fE << 342  
313     fIsDirectModelInitialised = true;          << 343  const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 344  const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
                                                   >> 345  G4double dum=0.;
                                                   >> 346  G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
                                                   >> 347  G4double dE=E2-E1;
                                                   >> 348  for (size_t i=0; i<material->GetNumberOfElements(); i++) { 
                                                   >> 349   G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
                                                   >> 350   G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
                                                   >> 351   dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
                                                   >> 352    
                                                   >> 353  }
                                                   >> 354  
                                                   >> 355  //Now the Migdal correction
                                                   >> 356  
                                                   >> 357  G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
                                                   >> 358  G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
                                                   >> 359                                              *(material->GetElectronDensity());
                                                   >> 360                  
                                                   >> 361  
                                                   >> 362  G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct
                                                   >> 363                     //model is different than the one used in the secondary sampling by a
                                                   >> 364                     //factor (1.+kp2) To be checked!
                                                   >> 365  
                                                   >> 366  dCrossEprod*=MigdalFactor;
                                                   >> 367  return dCrossEprod;
                                                   >> 368   
                                                   >> 369 }
                                                   >> 370 ////////////////////////////////////////////////////////////////////////////////
                                                   >> 371 //
                                                   >> 372 G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
                                                   >> 373                      G4double primEnergy,
                                                   >> 374                      G4bool IsScatProjToProjCase)
                                                   >> 375 { if (!isDirectModelInitialised) {
                                                   >> 376     theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
                                                   >> 377   isDirectModelInitialised =true;
314   }                                               378   }
315   if(fUseMatrix)                               << 379   if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
316     return G4VEmAdjointModel::AdjointCrossSect << 
317                                                << 
318   DefineCurrentMaterial(aCouple);                 380   DefineCurrentMaterial(aCouple);
319   G4double Cross = 0.;                         << 381   G4double Cross=0.;
320   // this gives the constant above             << 382   lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
321   fLastCZ = fDirectModel->CrossSectionPerVolum << 383   
322     aCouple->GetMaterial(), fDirectPrimaryPart << 384   if (!IsScatProjToProjCase ){
323                                                << 385     G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
324   if(!isScatProjToProj)                        << 386     G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
325   {                                            << 387   if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj);
326     G4double Emax_proj = GetSecondAdjEnergyMax << 388   }
327     G4double Emin_proj = GetSecondAdjEnergyMin << 389   else {
328     if(Emax_proj > Emin_proj && primEnergy > f << 390     G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
329       Cross = fCsBiasingFactor * fLastCZ * std << 391   G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
330   }                                            << 392   if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
331   else                                         << 393     
332   {                                            << 394   }
333     G4double Emax_proj = GetSecondAdjEnergyMax << 395   return Cross; 
334     G4double Emin_proj =                       << 396 }              
335       GetSecondAdjEnergyMinForScatProjToProj(p << 397 
336     if(Emax_proj > Emin_proj)                  << 398 G4double G4AdjointBremsstrahlungModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
337       Cross = fLastCZ * std::log((Emax_proj -  << 399                      G4double primEnergy,
338                                  Emax_proj / ( << 400                      G4bool IsScatProjToProjCase)
339   }                                            << 401 { 
340   return Cross;                                << 402   return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
                                                   >> 403   lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
                                                   >> 404   return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
                                                   >> 405     
341 }                                                 406 }
                                                   >> 407 
                                                   >> 408 
                                                   >> 409 
342                                                   410