Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointhIonisationModel.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/G4AdjointhIonisationModel.cc (Version 11.3.0) and /processes/electromagnetic/adjoint/src/G4AdjointhIonisationModel.cc (Version 9.6.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                                                <<  26 // $Id$
                                                   >>  27 //
 27 #include "G4AdjointhIonisationModel.hh"            28 #include "G4AdjointhIonisationModel.hh"
 28                                                    29 
                                                   >>  30 #include "G4PhysicalConstants.hh"
                                                   >>  31 #include "G4SystemOfUnits.hh"
 29 #include "G4AdjointCSManager.hh"                   32 #include "G4AdjointCSManager.hh"
                                                   >>  33 #include "G4Integrator.hh"
                                                   >>  34 #include "G4TrackStatus.hh"
                                                   >>  35 #include "G4ParticleChange.hh"
 30 #include "G4AdjointElectron.hh"                    36 #include "G4AdjointElectron.hh"
 31 #include "G4AdjointProton.hh"                      37 #include "G4AdjointProton.hh"
                                                   >>  38 #include "G4AdjointInterpolator.hh"
 32 #include "G4BetheBlochModel.hh"                    39 #include "G4BetheBlochModel.hh"
 33 #include "G4BraggModel.hh"                         40 #include "G4BraggModel.hh"
 34 #include "G4NistManager.hh"                    << 
 35 #include "G4ParticleChange.hh"                 << 
 36 #include "G4PhysicalConstants.hh"              << 
 37 #include "G4Proton.hh"                             41 #include "G4Proton.hh"
 38 #include "G4SystemOfUnits.hh"                  <<  42 #include "G4NistManager.hh"
 39 #include "G4TrackStatus.hh"                    << 
 40                                                    43 
 41 //////////////////////////////////////////////     44 ////////////////////////////////////////////////////////////////////////////////
 42 G4AdjointhIonisationModel::G4AdjointhIonisatio <<  45 //
 43   : G4VEmAdjointModel("Adjoint_hIonisation")   <<  46 G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* projectileDefinition):
 44 {                                              <<  47   G4VEmAdjointModel("Adjoint_hIonisation")
 45   fUseMatrix               = true;             <<  48 { 
 46   fUseMatrixPerElement     = true;             <<  49 
 47   fApplyCutInRange         = true;             <<  50 
 48   fOneMatrixForAllElements = true;             << 
 49   fSecondPartSameType      = false;            << 
 50                                                << 
 51   // The direct EM Model is taken as BetheBloc << 
 52   // computation of the differential cross sec << 
 53   // The Bragg model could be used as an alter << 
 54   // differential cross section                << 
 55                                                << 
 56   fDirectModel              = new G4BetheBloch << 
 57   fBraggDirectEMModel       = new G4BraggModel << 
 58   fAdjEquivDirectSecondPart = G4AdjointElectro << 
 59   fDirectPrimaryPart        = pDef;            << 
 60                                                << 
 61   if(pDef == G4Proton::Proton())               << 
 62   {                                            << 
 63     fAdjEquivDirectPrimPart = G4AdjointProton: << 
 64   }                                            << 
 65                                                    51 
                                                   >>  52   UseMatrix =true;
                                                   >>  53   UseMatrixPerElement = true;
                                                   >>  54   ApplyCutInRange = true;
                                                   >>  55   UseOnlyOneMatrixForAllElements = true; 
                                                   >>  56   CS_biasing_factor =1.;
                                                   >>  57   second_part_of_same_type =false;
                                                   >>  58   
                                                   >>  59   //The direct EM Modfel is taken has BetheBloch it is only used for the computation 
                                                   >>  60   // of the differential cross section.
                                                   >>  61   //The Bragg model could be used  as an alternative as  it offers the same differential cross section
                                                   >>  62   
                                                   >>  63   theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
                                                   >>  64   theBraggDirectEMModel = new G4BraggModel(projectileDefinition); 
                                                   >>  65   theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
                                                   >>  66   
                                                   >>  67   theDirectPrimaryPartDef = projectileDefinition;
                                                   >>  68   theAdjEquivOfDirectPrimPartDef = 0;
                                                   >>  69   if (projectileDefinition == G4Proton::Proton()) {
                                                   >>  70     theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
                                                   >>  71   
                                                   >>  72   }
                                                   >>  73   
 66   DefineProjectileProperty();                      74   DefineProjectileProperty();
 67 }                                                  75 }
 68                                                << 
 69 //////////////////////////////////////////////     76 ////////////////////////////////////////////////////////////////////////////////
 70 G4AdjointhIonisationModel::~G4AdjointhIonisati <<  77 //
                                                   >>  78 G4AdjointhIonisationModel::~G4AdjointhIonisationModel()
                                                   >>  79 {;}
                                                   >>  80 
 71                                                    81 
 72 //////////////////////////////////////////////     82 ////////////////////////////////////////////////////////////////////////////////
 73 void G4AdjointhIonisationModel::SampleSecondar <<  83 //
 74   const G4Track& aTrack, G4bool isScatProjToPr <<  84 void G4AdjointhIonisationModel::SampleSecondaries(const G4Track& aTrack,
 75   G4ParticleChange* fParticleChange)           <<  85                                 G4bool IsScatProjToProjCase,
 76 {                                              <<  86         G4ParticleChange* fParticleChange)
 77   if(!fUseMatrix)                              <<  87 { 
 78     return RapidSampleSecondaries(aTrack, isSc <<  88  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
 79                                                <<  89  
 80   const G4DynamicParticle* theAdjointPrimary = <<  90  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
 81                                                <<  91   
 82   // Elastic inverse scattering                <<  92  //Elastic inverse scattering 
 83   G4double adjointPrimKinEnergy = theAdjointPr <<  93  //---------------------------------------------------------
 84   G4double adjointPrimP         = theAdjointPr <<  94  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
 85                                                <<  95  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
 86   if(adjointPrimKinEnergy > GetHighEnergyLimit <<  96 
 87   {                                            <<  97  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
 88     return;                                    <<  98   return;
 89   }                                            <<  99  }
                                                   >> 100  
                                                   >> 101  //Sample secondary energy
                                                   >> 102  //-----------------------
                                                   >> 103  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
                                                   >> 104  CorrectPostStepWeight(fParticleChange, 
                                                   >> 105             aTrack.GetWeight(), 
                                                   >> 106             adjointPrimKinEnergy,
                                                   >> 107             projectileKinEnergy,
                                                   >> 108             IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
                                                   >> 109 
                                                   >> 110      
                                                   >> 111  //Kinematic: 
                                                   >> 112  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
                                                   >> 113  // him part of its  energy
                                                   >> 114  //----------------------------------------------------------------------------------------
                                                   >> 115  
                                                   >> 116  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
                                                   >> 117  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
                                                   >> 118  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 
                                                   >> 119  
                                                   >> 120  
                                                   >> 121  
                                                   >> 122  //Companion
                                                   >> 123  //-----------
                                                   >> 124  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
                                                   >> 125  if (IsScatProjToProjCase) {
                                                   >> 126     companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
                                                   >> 127  } 
                                                   >> 128  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
                                                   >> 129  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;  
                                                   >> 130  
                                                   >> 131  
                                                   >> 132  //Projectile momentum
                                                   >> 133  //--------------------
                                                   >> 134  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
                                                   >> 135  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
                                                   >> 136  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
                                                   >> 137  G4double phi =G4UniformRand()*2.*3.1415926;
                                                   >> 138  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
                                                   >> 139  projectileMomentum.rotateUz(dir_parallel);
                                                   >> 140  
                                                   >> 141  
                                                   >> 142  
                                                   >> 143  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
                                                   >> 144   fParticleChange->ProposeTrackStatus(fStopAndKill);
                                                   >> 145   fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
                                                   >> 146   //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
                                                   >> 147  }
                                                   >> 148  else {
                                                   >> 149   fParticleChange->ProposeEnergy(projectileKinEnergy);
                                                   >> 150   fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
                                                   >> 151  }  
                                                   >> 152       
                                                   >> 153 
                                                   >> 154   
 90                                                   155 
 91   // Sample secondary energy                   << 
 92   G4double projectileKinEnergy =               << 
 93     SampleAdjSecEnergyFromCSMatrix(adjointPrim << 
 94   CorrectPostStepWeight(fParticleChange, aTrac << 
 95                         adjointPrimKinEnergy,  << 
 96                         isScatProjToProj);     << 
 97   // Caution!!! this weight correction should  << 
 98                                                << 
 99   // Kinematic:                                << 
100   // we consider a two body elastic scattering << 
101   // the projectile knock on an e- at rest and << 
102   G4double projectileM0          = fAdjEquivDi << 
103   G4double projectileTotalEnergy = projectileM << 
104   G4double projectileP2 =                      << 
105     projectileTotalEnergy * projectileTotalEne << 
106                                                << 
107   // Companion                                 << 
108   G4double companionM0 = fAdjEquivDirectPrimPa << 
109   if(isScatProjToProj)                         << 
110   {                                            << 
111     companionM0 = fAdjEquivDirectSecondPart->G << 
112   }                                            << 
113   G4double companionTotalEnergy =              << 
114     companionM0 + projectileKinEnergy - adjoin << 
115   G4double companionP2 =                       << 
116     companionTotalEnergy * companionTotalEnerg << 
117                                                << 
118   // Projectile momentum                       << 
119   G4double P_parallel =                        << 
120     (adjointPrimP * adjointPrimP + projectileP << 
121     (2. * adjointPrimP);                       << 
122   G4double P_perp = std::sqrt(projectileP2 - P << 
123   G4ThreeVector dir_parallel = theAdjointPrima << 
124   G4double phi               = G4UniformRand() << 
125   G4ThreeVector projectileMomentum =           << 
126     G4ThreeVector(P_perp * std::cos(phi), P_pe << 
127   projectileMomentum.rotateUz(dir_parallel);   << 
128                                                << 
129   if(!isScatProjToProj)                        << 
130   {  // kill the primary and add a secondary   << 
131     fParticleChange->ProposeTrackStatus(fStopA << 
132     fParticleChange->AddSecondary(             << 
133       new G4DynamicParticle(fAdjEquivDirectPri << 
134   }                                            << 
135   else                                         << 
136   {                                            << 
137     fParticleChange->ProposeEnergy(projectileK << 
138     fParticleChange->ProposeMomentumDirection( << 
139   }                                            << 
140 }                                                 156 }
141                                                   157 
142 //////////////////////////////////////////////    158 ////////////////////////////////////////////////////////////////////////////////
143 void G4AdjointhIonisationModel::RapidSampleSec << 159 //
144   const G4Track& aTrack, G4bool isScatProjToPr << 160 void G4AdjointhIonisationModel::RapidSampleSecondaries(const G4Track& aTrack,
145   G4ParticleChange* fParticleChange)           << 161                        G4bool IsScatProjToProjCase,
146 {                                              << 162                  G4ParticleChange* fParticleChange)
147   const G4DynamicParticle* theAdjointPrimary = << 163 { 
148   DefineCurrentMaterial(aTrack.GetMaterialCuts << 164 
149                                                << 165  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
150   G4double adjointPrimKinEnergy = theAdjointPr << 166  DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
151   G4double adjointPrimP         = theAdjointPr << 167  
152                                                << 168  
153   if(adjointPrimKinEnergy > GetHighEnergyLimit << 169  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
154   {                                            << 170  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
155     return;                                    << 171  
156   }                                            << 172  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
                                                   >> 173   return;
                                                   >> 174  }
                                                   >> 175   
                                                   >> 176  G4double projectileKinEnergy =0.;
                                                   >> 177  G4double eEnergy=0.;
                                                   >> 178  G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
                                                   >> 179  if (!IsScatProjToProjCase){//1/E^2 distribution
                                                   >> 180   
                                                   >> 181   eEnergy=adjointPrimKinEnergy;
                                                   >> 182   G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
                                                   >> 183         G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
                                                   >> 184   if (Emin>=Emax) return;
                                                   >> 185   G4double a=1./Emax;
                                                   >> 186   G4double b=1./Emin;
                                                   >> 187   newCS=newCS*(b-a)/eEnergy;
                                                   >> 188   
                                                   >> 189   projectileKinEnergy =1./(b- (b-a)*G4UniformRand()); 
                                                   >> 190   
                                                   >> 191   
                                                   >> 192  }
                                                   >> 193  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
                                                   >> 194   G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
                                                   >> 195   if (Emin>=Emax) return;
                                                   >> 196   G4double diff1=Emin-adjointPrimKinEnergy;
                                                   >> 197   G4double diff2=Emax-adjointPrimKinEnergy;
                                                   >> 198   
                                                   >> 199   G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
                                                   >> 200   G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
                                                   >> 201   /*G4double f31=diff1/Emin;
                                                   >> 202   G4double f32=diff2/Emax/f31;*/
                                                   >> 203   G4double t3=2.*std::log(Emax/Emin);
                                                   >> 204   G4double sum_t=t1+t2+t3;
                                                   >> 205   newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
                                                   >> 206   G4double t=G4UniformRand()*sum_t;
                                                   >> 207   if (t <=t1 ){
                                                   >> 208     G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
                                                   >> 209     projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
                                                   >> 210     
                                                   >> 211   }
                                                   >> 212   else if (t <=t2 )  {
                                                   >> 213     G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
                                                   >> 214     projectileKinEnergy =1./(1./Emin-q);
                                                   >> 215   }
                                                   >> 216   else {  
                                                   >> 217     projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
                                                   >> 218     
                                                   >> 219   }
                                                   >> 220   eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
                                                   >> 221   
                                                   >> 222   
                                                   >> 223  }
                                                   >> 224 
                                                   >> 225  
                                                   >> 226 
                                                   >> 227  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy; 
                                                   >> 228   
                                                   >> 229   
                                                   >> 230                 
                                                   >> 231  //Weight correction
                                                   >> 232  //-----------------------
                                                   >> 233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
                                                   >> 234  G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
                                                   >> 235  
                                                   >> 236  //G4cout<<w_corr<<G4endl;
                                                   >> 237  w_corr*=newCS/lastCS;
                                                   >> 238  //G4cout<<w_corr<<G4endl;
                                                   >> 239  //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
                                                   >> 240  //Here we consider the true  diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
                                                   >> 241  
                                                   >> 242  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
                                                   >> 243  w_corr*=diffCS/diffCS_perAtom_Used;
                                                   >> 244  //G4cout<<w_corr<<G4endl;
                                                   >> 245      
                                                   >> 246  G4double new_weight = aTrack.GetWeight()*w_corr;
                                                   >> 247  fParticleChange->SetParentWeightByProcess(false);
                                                   >> 248  fParticleChange->SetSecondaryWeightByProcess(false);
                                                   >> 249  fParticleChange->ProposeParentWeight(new_weight);
                                                   >> 250  
                                                   >> 251  
                                                   >> 252  
                                                   >> 253  
                                                   >> 254  //Kinematic: 
                                                   >> 255  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest  and gives
                                                   >> 256  // him part of its  energy
                                                   >> 257  //----------------------------------------------------------------------------------------
                                                   >> 258  
                                                   >> 259  G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
                                                   >> 260  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
                                                   >> 261  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 
                                                   >> 262  
                                                   >> 263  
                                                   >> 264  
                                                   >> 265  //Companion
                                                   >> 266  //-----------
                                                   >> 267  G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
                                                   >> 268  if (IsScatProjToProjCase) {
                                                   >> 269     companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
                                                   >> 270  } 
                                                   >> 271  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
                                                   >> 272  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;  
                                                   >> 273  
                                                   >> 274  
                                                   >> 275  //Projectile momentum
                                                   >> 276  //--------------------
                                                   >> 277  G4double  P_parallel = (adjointPrimP*adjointPrimP +  projectileP2 - companionP2)/(2.*adjointPrimP);
                                                   >> 278  G4double  P_perp = std::sqrt( projectileP2 -  P_parallel*P_parallel);
                                                   >> 279  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
                                                   >> 280  G4double phi =G4UniformRand()*2.*3.1415926;
                                                   >> 281  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
                                                   >> 282  projectileMomentum.rotateUz(dir_parallel);
                                                   >> 283  
                                                   >> 284  
                                                   >> 285  
                                                   >> 286  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
                                                   >> 287   fParticleChange->ProposeTrackStatus(fStopAndKill);
                                                   >> 288   fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
                                                   >> 289   //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
                                                   >> 290  }
                                                   >> 291  else {
                                                   >> 292   fParticleChange->ProposeEnergy(projectileKinEnergy);
                                                   >> 293   fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
                                                   >> 294  }  
                                                   >> 295       
157                                                   296 
158   G4double projectileKinEnergy = 0.;           << 297  
159   G4double eEnergy             = 0.;           << 298  
160   G4double newCS =                             << 
161     fCurrentMaterial->GetElectronDensity() * t << 
162   if(!isScatProjToProj)                        << 
163   {  // 1/E^2 distribution                     << 
164                                                << 
165     eEnergy       = adjointPrimKinEnergy;      << 
166     G4double Emax = GetSecondAdjEnergyMaxForPr << 
167     G4double Emin = GetSecondAdjEnergyMinForPr << 
168     if(Emin >= Emax)                           << 
169       return;                                  << 
170     G4double a = 1. / Emax;                    << 
171     G4double b = 1. / Emin;                    << 
172     newCS      = newCS * (b - a) / eEnergy;    << 
173                                                   299 
174     projectileKinEnergy = 1. / (b - (b - a) *  << 
175   }                                            << 
176   else                                         << 
177   {                                            << 
178     G4double Emax =                            << 
179       GetSecondAdjEnergyMaxForScatProjToProj(a << 
180     G4double Emin =                            << 
181       GetSecondAdjEnergyMinForScatProjToProj(a << 
182     if(Emin >= Emax)                           << 
183       return;                                  << 
184     G4double diff1 = Emin - adjointPrimKinEner << 
185     G4double diff2 = Emax - adjointPrimKinEner << 
186                                                << 
187     G4double t1    = adjointPrimKinEnergy * (1 << 
188     G4double t2    = adjointPrimKinEnergy * (1 << 
189     G4double t3    = 2. * std::log(Emax / Emin << 
190     G4double sum_t = t1 + t2 + t3;             << 
191     newCS      = newCS * sum_t / adjointPrimKi << 
192     G4double t = G4UniformRand() * sum_t;      << 
193     if(t <= t1)                                << 
194     {                                          << 
195       G4double q          = G4UniformRand() *  << 
196       projectileKinEnergy = adjointPrimKinEner << 
197     }                                          << 
198     else if(t <= t2)                           << 
199     {                                          << 
200       G4double q          = G4UniformRand() *  << 
201       projectileKinEnergy = 1. / (1. / Emin -  << 
202     }                                          << 
203     else                                       << 
204     {                                          << 
205       projectileKinEnergy = Emin * std::pow(Em << 
206     }                                          << 
207     eEnergy = projectileKinEnergy - adjointPri << 
208   }                                            << 
209                                                   300 
210   G4double diffCS_perAtom_Used = twopi_mc2_rcl << 
211                                  projectileKin << 
212                                  eEnergy / eEn << 
213                                                << 
214   // Weight correction                         << 
215   // First w_corr is set to the ratio between  << 
216   G4double w_corr =                            << 
217     G4AdjointCSManager::GetAdjointCSManager()- << 
218                                                << 
219   w_corr *= newCS / fLastCS;                   << 
220   // Then another correction is needed due to  << 
221   // differential CS has been used rather than << 
222   // direct model. Here we consider the true d << 
223   // numerical differentiation over Tcut of th << 
224                                                << 
225   G4double diffCS =                            << 
226     DiffCrossSectionPerAtomPrimToSecond(projec << 
227   w_corr *= diffCS / diffCS_perAtom_Used;      << 
228                                                << 
229   if (isScatProjToProj && fTcutSecond>0.005) w << 
230                                                << 
231   G4double new_weight = aTrack.GetWeight() * w << 
232   fParticleChange->SetParentWeightByProcess(fa << 
233   fParticleChange->SetSecondaryWeightByProcess << 
234   fParticleChange->ProposeParentWeight(new_wei << 
235                                                << 
236   // Kinematic:                                << 
237   // we consider a two body elastic scattering << 
238   // the projectile knocks on an e- at rest an << 
239   G4double projectileM0          = fAdjEquivDi << 
240   G4double projectileTotalEnergy = projectileM << 
241   G4double projectileP2 =                      << 
242     projectileTotalEnergy * projectileTotalEne << 
243                                                << 
244   // Companion                                 << 
245   G4double companionM0 = fAdjEquivDirectPrimPa << 
246   if(isScatProjToProj)                         << 
247   {                                            << 
248     companionM0 = fAdjEquivDirectSecondPart->G << 
249   }                                            << 
250   G4double companionTotalEnergy =              << 
251     companionM0 + projectileKinEnergy - adjoin << 
252   G4double companionP2 =                       << 
253     companionTotalEnergy * companionTotalEnerg << 
254                                                << 
255   // Projectile momentum                       << 
256   G4double P_parallel =                        << 
257     (adjointPrimP * adjointPrimP + projectileP << 
258     (2. * adjointPrimP);                       << 
259   G4double P_perp = std::sqrt(projectileP2 - P << 
260   G4ThreeVector dir_parallel = theAdjointPrima << 
261   G4double phi               = G4UniformRand() << 
262   G4ThreeVector projectileMomentum =           << 
263     G4ThreeVector(P_perp * std::cos(phi), P_pe << 
264   projectileMomentum.rotateUz(dir_parallel);   << 
265                                                << 
266   if(!isScatProjToProj)                        << 
267   {  // kill the primary and add a secondary   << 
268     fParticleChange->ProposeTrackStatus(fStopA << 
269     fParticleChange->AddSecondary(             << 
270       new G4DynamicParticle(fAdjEquivDirectPri << 
271   }                                            << 
272   else                                         << 
273   {                                            << 
274     fParticleChange->ProposeEnergy(projectileK << 
275     fParticleChange->ProposeMomentumDirection( << 
276   }                                            << 
277 }                                              << 
278                                                   301 
                                                   >> 302 } 
                                                   >> 303     
279 //////////////////////////////////////////////    304 ////////////////////////////////////////////////////////////////////////////////
                                                   >> 305 //
280 G4double G4AdjointhIonisationModel::DiffCrossS    306 G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
281   G4double kinEnergyProj, G4double kinEnergyPr << 307                                       G4double kinEnergyProj, 
282 {  // Probably here the Bragg Model should be  << 308                                       G4double kinEnergyProd,
283    // kinEnergyProj/nuc < 2 MeV                << 309               G4double Z, 
284                                                << 310                                       G4double A)
285   G4double dSigmadEprod = 0.;                  << 311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
286   G4double Emax_proj    = GetSecondAdjEnergyMa << 312   
287   G4double Emin_proj    = GetSecondAdjEnergyMi << 313 
288                                                << 314  
289   // the produced particle should have a kinet << 315  G4double dSigmadEprod=0;
290   // projectile                                << 316  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
291   if(kinEnergyProj > Emin_proj && kinEnergyPro << 317  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
292   {                                            << 318  
293     G4double Tmax = kinEnergyProj;             << 319  
294     G4double E1   = kinEnergyProd;             << 320  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 
295     //1.0006 factor seems to give the best dif << 321   G4double Tmax=kinEnergyProj;
296     G4double E2   = kinEnergyProd *1.0006;     << 322   
297     G4double sigma1, sigma2;                   << 323         G4double E1=kinEnergyProd;
298     if(kinEnergyProj > 2. * MeV)               << 324   G4double E2=kinEnergyProd*1.000001;
299     {                                          << 325   G4double dE=(E2-E1);
300       sigma1 = fDirectModel->ComputeCrossSecti << 326   G4double sigma1,sigma2;
301         fDirectPrimaryPart, kinEnergyProj, Z,  << 327   if (kinEnergyProj >2.*MeV){
302       sigma2 = fDirectModel->ComputeCrossSecti << 328       sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
303         fDirectPrimaryPart, kinEnergyProj, Z,  << 329       sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
304     }                                          << 330         }
305     else                                       << 331   else {
306     {                                          << 332       sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
307       sigma1 = fBraggDirectEMModel->ComputeCro << 333       sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
308         fDirectPrimaryPart, kinEnergyProj, Z,  << 334   }
309       sigma2 = fBraggDirectEMModel->ComputeCro << 335   
310         fDirectPrimaryPart, kinEnergyProj, Z,  << 336   
311     }                                          << 337   dSigmadEprod=(sigma1-sigma2)/dE;
312                                                << 338   if (dSigmadEprod>1.) {
313     dSigmadEprod = (sigma1 - sigma2) / (E2 - E << 339     G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
314     if(dSigmadEprod > 1.)                      << 340     G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
315     {                                          << 341     G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
316       G4cout << "sigma1 " << kinEnergyProj / M << 342     
317              << '\t' << sigma1 << G4endl;      << 343   }
318       G4cout << "sigma2 " << kinEnergyProj / M << 344   
319              << '\t' << sigma2 << G4endl;      << 345    
320       G4cout << "dsigma " << kinEnergyProj / M << 346   
321              << '\t' << dSigmadEprod << G4endl << 347    //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
322     }                                          << 348    //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
323                                                << 349    //to test the rejection of a secondary
324     // correction of differential cross sectio << 350    //-------------------------
325     // the suppression of particle at secondar << 351    
326     // Bloch Model. This correction consists o << 352    //Source code taken from    G4BetheBlochModel::SampleSecondaries
327     // function used to test the rejection of  << 353    
328     // from G4BetheBlochModel::SampleSecondari << 354    G4double deltaKinEnergy = kinEnergyProd;
329     G4double deltaKinEnergy = kinEnergyProd;   << 355    
330                                                << 356    //Part of the taken code
331     // projectile formfactor - suppression of  << 357    //----------------------
332     // delta-electron production at high energ << 358    
333     G4double x = fFormFact * deltaKinEnergy;   << 359    
334     if(x > 1.e-6)                              << 360    
335     {                                          << 361    // projectile formfactor - suppresion of high energy
336       G4double totEnergy = kinEnergyProj + fMa << 362          // delta-electron production at high energy
337       G4double etot2     = totEnergy * totEner << 363    G4double x = formfact*deltaKinEnergy;
338       G4double beta2 = kinEnergyProj * (kinEne << 364           if(x > 1.e-6) {
339       G4double f     = 1.0 - beta2 * deltaKinE << 365     
340       G4double f1    = 0.0;                    << 366     
341       if(0.5 == fSpin)                         << 367     G4double totEnergy     = kinEnergyProj + mass;
342       {                                        << 368           G4double etot2         = totEnergy*totEnergy;
343         f1 = 0.5 * deltaKinEnergy * deltaKinEn << 369     G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
344         f += f1;                               << 370     G4double f;
345       }                                        << 371     G4double f1 = 0.0;
346       G4double x1 = 1.0 + x;                   << 372     f = 1.0 - beta2*deltaKinEnergy/Tmax;
347       G4double gg = 1.0 / (x1 * x1);           << 373         if( 0.5 == spin ) {
348       if(0.5 == fSpin)                         << 374             f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
349       {                                        << 375             f += f1;
350         G4double x2 = 0.5 * electron_mass_c2 * << 376         }
351         gg *= (1.0 + fMagMoment2 * (x2 - f1 /  << 377     G4double x1 = 1.0 + x;
352       }                                        << 378         G4double gg  = 1.0/(x1*x1);
353       if(gg > 1.0)                             << 379         if( 0.5 == spin ) {
354       {                                        << 380             G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
355         G4cout << "### G4BetheBlochModel in Ad << 381             gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
356                << G4endl;                      << 382         }
357         gg = 1.;                               << 383         if(gg > 1.0) {
358       }                                        << 384             G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
359       dSigmadEprod *= gg;                      << 385             << G4endl;
360     }                                          << 386       gg=1.;
                                                   >> 387     }
                                                   >> 388     //G4cout<<"gg"<<gg<<G4endl;
                                                   >> 389     dSigmadEprod*=gg;   
                                                   >> 390        }
                                                   >> 391    
361   }                                               392   }
362                                                   393 
363   return dSigmadEprod;                         << 394  return dSigmadEprod; 
364 }                                                 395 }
365                                                   396 
366 ////////////////////////////////////////////// << 397 
                                                   >> 398 
                                                   >> 399 //////////////////////////////////////////////////////////////////////////////////////////////
                                                   >> 400 //
367 void G4AdjointhIonisationModel::DefineProjecti    401 void G4AdjointhIonisationModel::DefineProjectileProperty()
368 {                                              << 402 {   
369   // Slightly modified code taken from G4Bethe << 403     //Slightly modified code taken from G4BetheBlochModel::SetParticle
370   G4String pname = fDirectPrimaryPart->GetPart << 404     //------------------------------------------------
371                                                << 405     G4String pname = theDirectPrimaryPartDef->GetParticleName();
372   fMass           = fDirectPrimaryPart->GetPDG << 406     if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
373   fSpin           = fDirectPrimaryPart->GetPDG << 407   pname != "deuteron" && pname != "triton") {
374   fMassRatio      = electron_mass_c2 / fMass;  << 408       isIon = true;
375   fOnePlusRatio2  = (1. + fMassRatio) * (1. +  << 
376   fOneMinusRatio2 = (1. - fMassRatio) * (1. -  << 
377   G4double magmom = fDirectPrimaryPart->GetPDG << 
378                     (0.5 * eplus * hbar_Planck << 
379   fMagMoment2 = magmom * magmom - 1.0;         << 
380   fFormFact   = 0.0;                           << 
381   if(fDirectPrimaryPart->GetLeptonNumber() ==  << 
382   {                                            << 
383     G4double x = 0.8426 * GeV;                 << 
384     if(fSpin == 0.0 && fMass < GeV)            << 
385     {                                          << 
386       x = 0.736 * GeV;                         << 
387     }                                          << 
388     else if(fMass > GeV)                       << 
389     {                                          << 
390       x /= G4NistManager::Instance()->GetZ13(f << 
391     }                                             409     }
392     fFormFact = 2.0 * electron_mass_c2 / (x *  << 410     
393   }                                            << 411     mass = theDirectPrimaryPartDef->GetPDGMass();
                                                   >> 412     mass_ratio_projectile = proton_mass_c2/theDirectPrimaryPartDef->GetPDGMass();;
                                                   >> 413     spin = theDirectPrimaryPartDef->GetPDGSpin();
                                                   >> 414     G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
                                                   >> 415     chargeSquare = q*q;
                                                   >> 416     ratio = electron_mass_c2/mass;
                                                   >> 417     ratio2 = ratio*ratio;
                                                   >> 418     one_plus_ratio_2=(1+ratio)*(1+ratio);
                                                   >> 419     one_minus_ratio_2=(1-ratio)*(1-ratio);
                                                   >> 420     G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
                                                   >> 421       *mass/(0.5*eplus*hbar_Planck*c_squared);
                                                   >> 422     magMoment2 = magmom*magmom - 1.0;
                                                   >> 423     formfact = 0.0;
                                                   >> 424     if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
                                                   >> 425       G4double x = 0.8426*GeV;
                                                   >> 426       if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
                                                   >> 427       else if(mass > GeV) {
                                                   >> 428         x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
                                                   >> 429   //  tlimit = 51.2*GeV*A13[iz]*A13[iz];
                                                   >> 430       }
                                                   >> 431       formfact = 2.0*electron_mass_c2/(x*x);
                                                   >> 432       tlimit   = 2.0/formfact;
                                                   >> 433    }
394 }                                                 434 }
395                                                   435 
396 //////////////////////////////////////////////    436 ////////////////////////////////////////////////////////////////////////////////
397 G4double G4AdjointhIonisationModel::AdjointCro << 437 //
398   const G4MaterialCutsCouple* aCouple, G4doubl << 438 G4double G4AdjointhIonisationModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
399   G4bool isScatProjToProj)                     << 439                      G4double primEnergy,
400 {                                              << 440                      G4bool IsScatProjToProjCase)
401   if(fUseMatrix)                               << 441 { 
402     return G4VEmAdjointModel::AdjointCrossSect << 442   if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
403                                                << 
404   DefineCurrentMaterial(aCouple);                 443   DefineCurrentMaterial(aCouple);
                                                   >> 444     
                                                   >> 445   
                                                   >> 446   G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
                                                   >> 447   
                                                   >> 448   if (!IsScatProjToProjCase ){
                                                   >> 449     G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
                                                   >> 450     G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
                                                   >> 451   if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
                                                   >> 452     Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
                                                   >> 453   } 
                                                   >> 454   else Cross=0.;
                                                   >> 455     
                                                   >> 456   
                                                   >> 457 
                                                   >> 458   
                                                   >> 459   
                                                   >> 460   
                                                   >> 461   }
                                                   >> 462   else {
                                                   >> 463     G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
                                                   >> 464   G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
                                                   >> 465   G4double diff1=Emin_proj-primEnergy;
                                                   >> 466   G4double diff2=Emax_proj-primEnergy;
                                                   >> 467   G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
                                                   >> 468   //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
                                                   >> 469   G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
                                                   >> 470   Cross*=(t1+t2);
405                                                   471 
406   G4double Cross =                             << 472     
407     fCurrentMaterial->GetElectronDensity() * t << 
408                                                << 
409   if(!isScatProjToProj)                        << 
410   {                                            << 
411     G4double Emax_proj = GetSecondAdjEnergyMax << 
412     G4double Emin_proj = GetSecondAdjEnergyMin << 
413     if(Emax_proj > Emin_proj && primEnergy > f << 
414     {                                          << 
415       Cross *= (1. / Emin_proj - 1. / Emax_pro << 
416     }                                          << 
417     else                                       << 
418       Cross = 0.;                              << 
419   }                                            << 
420   else                                         << 
421   {                                            << 
422     G4double Emax_proj = GetSecondAdjEnergyMax << 
423     G4double Emin_proj =                       << 
424       GetSecondAdjEnergyMinForScatProjToProj(p << 
425     G4double diff1 = Emin_proj - primEnergy;   << 
426     G4double diff2 = Emax_proj - primEnergy;   << 
427     G4double t1 =                              << 
428       (1. / diff1 + 1. / Emin_proj - 1. / diff << 
429     G4double t2 =                              << 
430       2. * std::log(Emax_proj / Emin_proj) / p << 
431     Cross *= (t1 + t2);                        << 
432   }                                               473   }
433   fLastCS = Cross;                             << 474   lastCS =Cross;
434   return Cross;                                << 475   return Cross; 
435 }                                                 476 }
436                                                << 
437 //////////////////////////////////////////////    477 //////////////////////////////////////////////////////////////////////////////
438 G4double G4AdjointhIonisationModel::GetSecondA << 478 //
439   G4double primAdjEnergy)                      << 479 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
440 {                                              << 480 { 
441   G4double Tmax = primAdjEnergy * fOnePlusRati << 481   G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
442                   (fOneMinusRatio2 - 2. * fMas << 
443   return Tmax;                                    482   return Tmax;
444 }                                                 483 }
445                                                << 
446 //////////////////////////////////////////////    484 //////////////////////////////////////////////////////////////////////////////
447 G4double G4AdjointhIonisationModel::GetSecondA << 485 //
448   G4double primAdjEnergy, G4double tcut)       << 486 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
449 {                                              << 487 { return PrimAdjEnergy+Tcut;
450   return primAdjEnergy + tcut;                 << 
451 }                                                 488 }
452                                                << 
453 //////////////////////////////////////////////    489 //////////////////////////////////////////////////////////////////////////////
454 G4double G4AdjointhIonisationModel::GetSecondA << 490 //        
455 {                                              << 491 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
456   return GetHighEnergyLimit();                 << 492 { return HighEnergyLimit;
457 }                                                 493 }
458                                                << 
459 //////////////////////////////////////////////    494 //////////////////////////////////////////////////////////////////////////////
460 G4double G4AdjointhIonisationModel::GetSecondA << 495 //
461   G4double primAdjEnergy)                      << 496 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
462 {                                              << 497 {  G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
463   G4double Tmin =                              << 498    return Tmin;
464     (2. * primAdjEnergy - 4. * fMass +         << 
465      std::sqrt(4. * primAdjEnergy * primAdjEne << 
466                8. * primAdjEnergy * fMass * (1 << 
467     4.;                                        << 
468   return Tmin;                                 << 
469 }                                                 499 }
470                                                   500