Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointCSManager.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/G4AdjointCSManager.cc (Version 11.3.0) and /processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc (Version 9.6.p3)


  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 // $Id: G4AdjointCSManager.cc 69844 2013-05-16 09:19:33Z gcosmo $
 26 //                                                 27 //
 27                                                    28 
                                                   >>  29 #include <fstream>
                                                   >>  30 #include <iomanip>
                                                   >>  31 
 28 #include "G4AdjointCSManager.hh"                   32 #include "G4AdjointCSManager.hh"
 29                                                    33 
                                                   >>  34 #include "G4PhysicalConstants.hh"
                                                   >>  35 #include "G4SystemOfUnits.hh"
 30 #include "G4AdjointCSMatrix.hh"                    36 #include "G4AdjointCSMatrix.hh"
 31 #include "G4AdjointElectron.hh"                << 
 32 #include "G4AdjointGamma.hh"                   << 
 33 #include "G4AdjointInterpolator.hh"                37 #include "G4AdjointInterpolator.hh"
 34 #include "G4AdjointProton.hh"                  <<  38 #include "G4AdjointCSMatrix.hh"
 35 #include "G4Electron.hh"                       <<  39 #include "G4VEmAdjointModel.hh"
 36 #include "G4Element.hh"                        << 
 37 #include "G4ElementTable.hh"                       40 #include "G4ElementTable.hh"
 38 #include "G4Gamma.hh"                          <<  41 #include "G4Element.hh"
 39 #include "G4ParticleDefinition.hh"                 42 #include "G4ParticleDefinition.hh"
 40 #include "G4PhysicalConstants.hh"              <<  43 #include "G4Element.hh"
 41 #include "G4PhysicsLogVector.hh"               << 
 42 #include "G4PhysicsTable.hh"                   << 
 43 #include "G4ProductionCutsTable.hh"            << 
 44 #include "G4Proton.hh"                         << 
 45 #include "G4SystemOfUnits.hh"                  << 
 46 #include "G4VEmAdjointModel.hh"                << 
 47 #include "G4VEmProcess.hh"                         44 #include "G4VEmProcess.hh"
 48 #include "G4VEnergyLossProcess.hh"                 45 #include "G4VEnergyLossProcess.hh"
                                                   >>  46 #include "G4PhysicsTable.hh" 
                                                   >>  47 #include "G4PhysicsLogVector.hh"
                                                   >>  48 #include "G4PhysicsTableHelper.hh"
                                                   >>  49 #include "G4Electron.hh"
                                                   >>  50 #include "G4Gamma.hh"
                                                   >>  51 #include "G4Proton.hh"
                                                   >>  52 #include "G4AdjointElectron.hh"
                                                   >>  53 #include "G4AdjointGamma.hh"
                                                   >>  54 #include "G4AdjointProton.hh"
                                                   >>  55 #include "G4ProductionCutsTable.hh"
                                                   >>  56 #include "G4ProductionCutsTable.hh"
 49                                                    57 
 50 G4ThreadLocal G4AdjointCSManager* G4AdjointCSM <<  58 G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
 51                                                << 
 52 constexpr G4double G4AdjointCSManager::fTmin;  << 
 53 constexpr G4double G4AdjointCSManager::fTmax;  << 
 54 constexpr G4int G4AdjointCSManager::fNbins;    << 
 55                                                << 
 56 //////////////////////////////////////////////     59 ///////////////////////////////////////////////////////
                                                   >>  60 //
 57 G4AdjointCSManager* G4AdjointCSManager::GetAdj     61 G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
 58 {                                              <<  62 { if(theInstance == 0) {
 59   if(fInstance == nullptr)                     <<  63     static G4AdjointCSManager ins;
 60   {                                            <<  64      theInstance = &ins;
 61     static G4ThreadLocalSingleton<G4AdjointCSM << 
 62     fInstance = inst.Instance();               << 
 63   }                                                65   }
 64   return fInstance;                            <<  66  return theInstance; 
 65 }                                                  67 }
 66                                                    68 
 67 //////////////////////////////////////////////     69 ///////////////////////////////////////////////////////
                                                   >>  70 //
 68 G4AdjointCSManager::G4AdjointCSManager()           71 G4AdjointCSManager::G4AdjointCSManager()
 69 {                                              <<  72 { CrossSectionMatrixesAreBuilt=false;
                                                   >>  73   TotalSigmaTableAreBuilt=false;
                                                   >>  74   theTotalForwardSigmaTableVector.clear();
                                                   >>  75   theTotalAdjointSigmaTableVector.clear();
                                                   >>  76   listOfForwardEmProcess.clear();
                                                   >>  77   listOfForwardEnergyLossProcess.clear();
                                                   >>  78   theListOfAdjointParticlesInAction.clear(); 
                                                   >>  79   EminForFwdSigmaTables.clear();
                                                   >>  80   EminForAdjSigmaTables.clear();
                                                   >>  81   EkinofFwdSigmaMax.clear();
                                                   >>  82   EkinofAdjSigmaMax.clear();
                                                   >>  83   listSigmaTableForAdjointModelScatProjToProj.clear();
                                                   >>  84   listSigmaTableForAdjointModelProdToProj.clear();
                                                   >>  85   Tmin=0.1*keV;
                                                   >>  86   Tmax=100.*TeV;
                                                   >>  87   nbins=320; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
                                                   >>  88   
 70   RegisterAdjointParticle(G4AdjointElectron::A     89   RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
 71   RegisterAdjointParticle(G4AdjointGamma::Adjo     90   RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
 72   RegisterAdjointParticle(G4AdjointProton::Adj     91   RegisterAdjointParticle(G4AdjointProton::AdjointProton());
 73 }                                              <<  92   
                                                   >>  93   verbose  = 1;
                                                   >>  94  
                                                   >>  95   lastPartDefForCS =0;
                                                   >>  96   LastEkinForCS =0;
                                                   >>  97   LastCSCorrectionFactor =1.;
                                                   >>  98   
                                                   >>  99   forward_CS_mode = true;
                                                   >> 100   
                                                   >> 101   currentParticleDef = 0;
                                                   >> 102   currentCouple =0;
                                                   >> 103   currentMaterial=0;
                                                   >> 104   lastMaterial=0;
                                                   >> 105 
                                                   >> 106   
                                                   >> 107   theAdjIon = 0;
                                                   >> 108   theFwdIon = 0;  
                                                   >> 109 
 74                                                   110 
                                                   >> 111  
                                                   >> 112  
                                                   >> 113 }
 75 //////////////////////////////////////////////    114 ///////////////////////////////////////////////////////
                                                   >> 115 //
 76 G4AdjointCSManager::~G4AdjointCSManager()         116 G4AdjointCSManager::~G4AdjointCSManager()
 77 {                                              << 117 {;
 78   for (auto& p : fAdjointCSMatricesForProdToPr << 
 79     for (auto p1 : p) {                        << 
 80       if (p1) {                                << 
 81         delete p1;                             << 
 82         p1 = nullptr;                          << 
 83       }                                        << 
 84     }                                          << 
 85     p.clear();                                 << 
 86   }                                            << 
 87   fAdjointCSMatricesForProdToProj.clear();     << 
 88                                                << 
 89   for (auto& p : fAdjointCSMatricesForScatProj << 
 90     for (auto p1 : p) {                        << 
 91       if (p1) {                                << 
 92         delete p1;                             << 
 93         p1 = nullptr;                          << 
 94       }                                        << 
 95     }                                          << 
 96     p.clear();                                 << 
 97   }                                            << 
 98   fAdjointCSMatricesForScatProjToProj.clear(); << 
 99                                                << 
100   for (auto p : fAdjointModels) {              << 
101     if (p) {                                   << 
102       delete p;                                << 
103       p = nullptr;                             << 
104     }                                          << 
105   }                                            << 
106   fAdjointModels.clear();                      << 
107                                                << 
108   for (auto p : fTotalAdjSigmaTable) {         << 
109     p->clearAndDestroy();                      << 
110     delete p;                                  << 
111     p = nullptr;                               << 
112   }                                            << 
113   fTotalAdjSigmaTable.clear();                 << 
114                                                << 
115   for (auto p : fSigmaTableForAdjointModelScat << 
116     p->clearAndDestroy();                      << 
117     delete p;                                  << 
118     p = nullptr;                               << 
119   }                                            << 
120   fSigmaTableForAdjointModelScatProjToProj.cle << 
121                                                << 
122   for (auto p : fSigmaTableForAdjointModelProd << 
123     p->clearAndDestroy();                      << 
124     delete p;                                  << 
125     p = nullptr;                               << 
126   }                                            << 
127   fSigmaTableForAdjointModelProdToProj.clear() << 
128                                                << 
129   for (auto p : fTotalFwdSigmaTable) {         << 
130     p->clearAndDestroy();                      << 
131     delete p;                                  << 
132     p = nullptr;                               << 
133   }                                            << 
134   fTotalFwdSigmaTable.clear();                 << 
135                                                << 
136   for (auto p : fForwardProcesses) {           << 
137     delete p;                                  << 
138     p = nullptr;                               << 
139   }                                            << 
140   fForwardProcesses.clear();                   << 
141                                                << 
142   for (auto p : fForwardLossProcesses) {       << 
143     delete p;                                  << 
144     p = nullptr;                               << 
145   }                                            << 
146   fForwardLossProcesses.clear();               << 
147 }                                                 118 }
148                                                << 
149 //////////////////////////////////////////////    119 ///////////////////////////////////////////////////////
150 std::size_t G4AdjointCSManager::RegisterEmAdjo << 120 //
151 {                                              << 121 size_t G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
152   fAdjointModels.push_back(aModel);            << 122 {listOfAdjointEMModel.push_back(aModel);
153   fSigmaTableForAdjointModelScatProjToProj.pus << 123  listSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable);
154   fSigmaTableForAdjointModelProdToProj.push_ba << 124  listSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable);
155   return fAdjointModels.size() - 1;            << 125  return listOfAdjointEMModel.size() -1;
                                                   >> 126  
156 }                                                 127 }
157                                                << 
158 //////////////////////////////////////////////    128 ///////////////////////////////////////////////////////
159 void G4AdjointCSManager::RegisterEmProcess(G4V << 129 //
160                                            G4P << 130 void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
161 {                                              << 131 { 
162   G4ParticleDefinition* anAdjPartDef =         << 132   G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
163     GetAdjointParticleEquivalent(aFwdPartDef); << 133   if (anAdjPartDef && aProcess){
164   if(anAdjPartDef && aProcess)                 << 134     RegisterAdjointParticle(anAdjPartDef);
165   {                                            << 135   G4int index=-1;
166     RegisterAdjointParticle(anAdjPartDef);     << 136   
167                                                << 137     for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
168     for(std::size_t i = 0; i < fAdjointParticl << 138       if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
169     {                                          << 139     }
170       if(anAdjPartDef->GetParticleName() ==    << 140     listOfForwardEmProcess[index]->push_back(aProcess);
171          fAdjointParticlesInAction[i]->GetPart << 
172         fForwardProcesses[i]->push_back(aProce << 
173     }                                          << 
174   }                                               141   }
175 }                                                 142 }
176                                                << 
177 //////////////////////////////////////////////    143 ///////////////////////////////////////////////////////
178 void G4AdjointCSManager::RegisterEnergyLossPro << 144 //
179   G4VEnergyLossProcess* aProcess, G4ParticleDe << 145 void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
180 {                                                 146 {
181   G4ParticleDefinition* anAdjPartDef =         << 147   G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
182     GetAdjointParticleEquivalent(aFwdPartDef); << 148   if (anAdjPartDef && aProcess){
183   if(anAdjPartDef && aProcess)                 << 149         RegisterAdjointParticle(anAdjPartDef);
184   {                                            << 150     G4int index=-1;
185     RegisterAdjointParticle(anAdjPartDef);     << 151     for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
186     for(std::size_t i = 0; i < fAdjointParticl << 152       if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
187     {                                          << 153     }
188       if(anAdjPartDef->GetParticleName() ==    << 154     listOfForwardEnergyLossProcess[index]->push_back(aProcess);
189          fAdjointParticlesInAction[i]->GetPart << 155    }
190                               fForwardLossProc << 
191                                                << 
192     }                                          << 
193   }                                            << 
194 }                                                 156 }
195                                                << 
196 //////////////////////////////////////////////    157 ///////////////////////////////////////////////////////
                                                   >> 158 //
197 void G4AdjointCSManager::RegisterAdjointPartic    159 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
198 {                                              << 160 {  G4int index=-1;
199   G4bool found = false;                        << 161    for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
200   for(auto p : fAdjointParticlesInAction)      << 162     if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
201   {                                            << 163    }
202     if(p->GetParticleName() == aPartDef->GetPa << 164   
203     {                                          << 165    if (index ==-1){
204       found = true;                            << 166     listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>());
205     }                                          << 167   theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable);
206   }                                            << 168   theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
207   if(!found)                                   << 169   listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
208   {                                            << 170   theListOfAdjointParticlesInAction.push_back(aPartDef);
209     fForwardLossProcesses.push_back(new std::v << 171   EminForFwdSigmaTables.push_back(std::vector<G4double> ());
210     fTotalFwdSigmaTable.push_back(new G4Physic << 172   EminForAdjSigmaTables.push_back(std::vector<G4double> ());
211     fTotalAdjSigmaTable.push_back(new G4Physic << 173   EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
212     fForwardProcesses.push_back(new std::vecto << 174   EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
213     fAdjointParticlesInAction.push_back(aPartD << 175   
214     fEminForFwdSigmaTables.push_back(std::vect << 176    }
215     fEminForAdjSigmaTables.push_back(std::vect << 
216     fEkinofFwdSigmaMax.push_back(std::vector<G << 
217     fEkinofAdjSigmaMax.push_back(std::vector<G << 
218   }                                            << 
219 }                                                 177 }
220                                                << 
221 //////////////////////////////////////////////    178 ///////////////////////////////////////////////////////
                                                   >> 179 //
222 void G4AdjointCSManager::BuildCrossSectionMatr    180 void G4AdjointCSManager::BuildCrossSectionMatrices()
223 {                                              << 181 { 
224   if(fCSMatricesBuilt)                         << 182   if (CrossSectionMatrixesAreBuilt) return;
225     return;                                    << 183     //Tcut, Tmax 
226   // The Tcut and Tmax matrices will be comput << 184       //The matrices will be computed probably just once
227   // When Tcut changes, some PhysicsTable will << 185         //When Tcut will change some PhysicsTable will be recomputed  
228   // for each MaterialCutCouple but not all th << 186         // for each MaterialCutCouple but not all the matrices  
229   // The Tcut defines a lower limit in the ene << 187         //The Tcut defines a lower limit in the energy of the Projectile before the scattering
230   // scattering. In the Projectile to Scattere << 188         //In the Projectile to Scattered Projectile case we have
231   //      E_ScatProj<E_Proj-Tcut               << 189         //      E_ScatProj<E_Proj-Tcut
232   // Therefore in the adjoint case we have     << 190         //Therefore in the adjoint case we have
233   //      Eproj> E_ScatProj+Tcut               << 191         //      Eproj> E_ScatProj+Tcut
234   // This implies that when computing the adjo << 192         //This implies that when computing the adjoint CS we should integrate over Epro
235   // Epro from E_ScatProj+Tcut to Emax         << 193         // from E_ScatProj+Tcut to Emax
236   // In the Projectile to Secondary case Tcut  << 194         //In the Projectile to Secondary case Tcut plays a role only in the fact that  
237   // Esecond should be greater than Tcut to ha << 195         // Esecond should be greater than Tcut to have the possibility to have any adjoint
238   // adjoint process.                          << 196         //process      
239   // To avoid recomputing matrices for all cha << 197         //To avoid to recompute the matrices for all changes of MaterialCutCouple
240   // we propose to compute the matrices only o << 198         //We propose to compute the matrices only once for the minimum possible Tcut and then
241   // and then to interpolate the probability f << 199         //to interpolate the probability for a new Tcut (implemented in G4VAdjointEmModel)
242   // G4VAdjointEmModel)                        << 200   
243                                                << 201   
244   fAdjointCSMatricesForScatProjToProj.clear(); << 202   theAdjointCSMatricesForScatProjToProj.clear();
245   fAdjointCSMatricesForProdToProj.clear();     << 203   theAdjointCSMatricesForProdToProj.clear();
246   const G4ElementTable* theElementTable   = G4 << 204   const G4ElementTable* theElementTable = G4Element::GetElementTable();
247   const G4MaterialTable* theMaterialTable = G4 << 205   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
248                                                << 206   
249   G4cout << "========== Computation of cross s << 207   G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
250             "models =========="                << 208   for (size_t i=0; i<listOfAdjointEMModel.size();i++){
251          << G4endl;                            << 209     G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
252   for(const auto& aModel : fAdjointModels)     << 210     G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
253   {                                            << 211     if (aModel->GetUseMatrix()){
254     G4cout << "Build adjoint cross section mat << 212       std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
255            << G4endl;                          << 213       std::vector<G4AdjointCSMatrix*>* aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
256     if(aModel->GetUseMatrix())                 << 214       aListOfMat1->clear();
257     {                                          << 215       aListOfMat2->clear();
258       std::vector<G4AdjointCSMatrix*>* aListOf << 216       if (aModel->GetUseMatrixPerElement()){
259         new std::vector<G4AdjointCSMatrix*>(); << 217         if (aModel->GetUseOnlyOneMatrixForAllElements()){
260       std::vector<G4AdjointCSMatrix*>* aListOf << 218             std::vector<G4AdjointCSMatrix*>
261         new std::vector<G4AdjointCSMatrix*>(); << 219             two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
262       if(aModel->GetUseMatrixPerElement())     << 220             aListOfMat1->push_back(two_matrices[0]);
263       {                                        << 221             aListOfMat2->push_back(two_matrices[1]);
264         if(aModel->GetUseOnlyOneMatrixForAllEl << 222         }
265         {                                      << 223         else {    
266           std::vector<G4AdjointCSMatrix*> two_ << 224           for (size_t j=0; j<theElementTable->size();j++){
267             BuildCrossSectionsModelAndElement( << 225             G4Element* anElement=(*theElementTable)[j];
268           aListOfMat1->push_back(two_matrices[ << 226             G4int Z = int(anElement->GetZ());
269           aListOfMat2->push_back(two_matrices[ << 227             G4int A = int(anElement->GetA());
270         }                                      << 228             std::vector<G4AdjointCSMatrix*>
271         else                                   << 229               two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
272         {                                      << 230             aListOfMat1->push_back(two_matrices[0]);
273           for(const auto& anElement : *theElem << 231             aListOfMat2->push_back(two_matrices[1]);
274           {                                    << 232           }
275             G4int Z = G4lrint(anElement->GetZ( << 233         } 
276             G4int A = G4lrint(anElement->GetN( << 234       }
277             std::vector<G4AdjointCSMatrix*> tw << 235       else { //Per material case
278               BuildCrossSectionsModelAndElemen << 236         for (size_t j=0; j<theMaterialTable->size();j++){
279             aListOfMat1->push_back(two_matrice << 237           G4Material* aMaterial=(*theMaterialTable)[j];
280             aListOfMat2->push_back(two_matrice << 238           std::vector<G4AdjointCSMatrix*>
281           }                                    << 239             two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
282         }                                      << 240           aListOfMat1->push_back(two_matrices[0]);
283       }                                        << 241           aListOfMat2->push_back(two_matrices[1]);
284       else                                     << 242         }
285       {  // Per material case                  << 243       
286         for(const auto& aMaterial : *theMateri << 244       }
287         {                                      << 245       theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
288           std::vector<G4AdjointCSMatrix*> two_ << 246       theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);  
289             BuildCrossSectionsModelAndMaterial << 247       aModel->SetCSMatrices(aListOfMat1, aListOfMat2);  
290           aListOfMat1->push_back(two_matrices[ << 248     }
291           aListOfMat2->push_back(two_matrices[ << 249     else {  G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
292         }                                      << 250       std::vector<G4AdjointCSMatrix*> two_empty_matrices;
293       }                                        << 251       theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
294       fAdjointCSMatricesForProdToProj.push_bac << 252       theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
295       fAdjointCSMatricesForScatProjToProj.push << 253       
296       aModel->SetCSMatrices(aListOfMat1, aList << 254     }   
297     }                                          << 255   }
298     else                                       << 256   G4cout<<"              All adjoint cross section matrices are computed!"<<G4endl;
299     {                                          << 257   G4cout<<"======================================================================"<<G4endl;
300       G4cout << "The model " << aModel->GetNam << 258   
301              << " does not use cross section m << 259   CrossSectionMatrixesAreBuilt = true;
302       std::vector<G4AdjointCSMatrix*> two_empt << 260 
303       fAdjointCSMatricesForProdToProj.push_bac << 
304       fAdjointCSMatricesForScatProjToProj.push << 
305     }                                          << 
306   }                                            << 
307   G4cout << "              All adjoint cross s << 
308          << G4endl;                            << 
309   G4cout                                       << 
310     << "====================================== << 
311     << G4endl;                                 << 
312                                                   261 
313   fCSMatricesBuilt = true;                     << 
314 }                                                 262 }
315                                                   263 
                                                   >> 264 
316 //////////////////////////////////////////////    265 ///////////////////////////////////////////////////////
                                                   >> 266 //
317 void G4AdjointCSManager::BuildTotalSigmaTables    267 void G4AdjointCSManager::BuildTotalSigmaTables()
318 {                                              << 268 { if (TotalSigmaTableAreBuilt) return;
319   if(fSigmaTableBuilt)                         << 
320     return;                                    << 
321                                                << 
322   const G4ProductionCutsTable* theCoupleTable  << 
323     G4ProductionCutsTable::GetProductionCutsTa << 
324                                                << 
325   // Prepare the Sigma table for all AdjointEM << 
326   for(std::size_t i = 0; i < fAdjointModels.si << 
327   {                                            << 
328     fSigmaTableForAdjointModelScatProjToProj[i << 
329     fSigmaTableForAdjointModelProdToProj[i]->c << 
330     for(std::size_t j = 0; j < theCoupleTable- << 
331     {                                          << 
332       fSigmaTableForAdjointModelScatProjToProj << 
333         new G4PhysicsLogVector(fTmin, fTmax, f << 
334       fSigmaTableForAdjointModelProdToProj[i]- << 
335         new G4PhysicsLogVector(fTmin, fTmax, f << 
336     }                                          << 
337   }                                            << 
338                                                << 
339   for(std::size_t i = 0; i < fAdjointParticles << 
340   {                                            << 
341     G4ParticleDefinition* thePartDef = fAdjoin << 
342     DefineCurrentParticle(thePartDef);         << 
343     fTotalFwdSigmaTable[i]->clearAndDestroy(); << 
344     fTotalAdjSigmaTable[i]->clearAndDestroy(); << 
345     fEminForFwdSigmaTables[i].clear();         << 
346     fEminForAdjSigmaTables[i].clear();         << 
347     fEkinofFwdSigmaMax[i].clear();             << 
348     fEkinofAdjSigmaMax[i].clear();             << 
349                                                << 
350     for(std::size_t j = 0; j < theCoupleTable- << 
351     {                                          << 
352       const G4MaterialCutsCouple* couple =     << 
353         theCoupleTable->GetMaterialCutsCouple( << 
354                                                << 
355       // make first the total fwd CS table for << 
356       G4PhysicsVector* aVector = new G4Physics << 
357       G4bool Emin_found        = false;        << 
358       G4double sigma_max       = 0.;           << 
359       G4double e_sigma_max     = 0.;           << 
360       for(std::size_t l = 0; l < fNbins; ++l)  << 
361       {                                        << 
362         G4double totCS = 0.;                   << 
363         G4double e     = aVector->Energy(l);   << 
364         for(std::size_t k = 0; k < fForwardPro << 
365         {                                      << 
366           totCS += (*fForwardProcesses[i])[k]- << 
367         }                                      << 
368         for(std::size_t k = 0; k < fForwardLos << 
369         {                                      << 
370           if(thePartDef == fAdjIon)            << 
371           {  // e is considered already as the << 
372             std::size_t mat_index = couple->Ge << 
373             G4VEmModel* currentModel =         << 
374               (*fForwardLossProcesses[i])[k]-> << 
375                                                << 
376             G4double chargeSqRatio = currentMo << 
377               fFwdIon, couple->GetMaterial(),  << 
378             (*fForwardLossProcesses[i])[k]->Se << 
379                                                << 
380           }                                    << 
381           G4double e1 = e / fMassRatio;        << 
382           totCS += (*fForwardLossProcesses[i]) << 
383         }                                      << 
384         aVector->PutValue(l, totCS);           << 
385         if(totCS > sigma_max)                  << 
386         {                                      << 
387           sigma_max   = totCS;                 << 
388           e_sigma_max = e;                     << 
389         }                                      << 
390         if(totCS > 0 && !Emin_found)           << 
391         {                                      << 
392           fEminForFwdSigmaTables[i].push_back( << 
393           Emin_found = true;                   << 
394         }                                      << 
395       }                                        << 
396                                                << 
397       fEkinofFwdSigmaMax[i].push_back(e_sigma_ << 
398                                                << 
399       if(!Emin_found)                          << 
400         fEminForFwdSigmaTables[i].push_back(fT << 
401                                                << 
402       fTotalFwdSigmaTable[i]->push_back(aVecto << 
403                                                << 
404       Emin_found                = false;       << 
405       sigma_max                 = 0;           << 
406       e_sigma_max               = 0.;          << 
407       G4PhysicsVector* aVector1 = new G4Physic << 
408       for(std::size_t eindex = 0; eindex < fNb << 
409       {                                        << 
410         G4double e     = aVector1->Energy(eind << 
411         G4double totCS = ComputeTotalAdjointCS << 
412           couple, thePartDef,                  << 
413           e * 0.9999999 / fMassRatio);  // fMa << 
414         aVector1->PutValue(eindex, totCS);     << 
415         if(totCS > sigma_max)                  << 
416         {                                      << 
417           sigma_max   = totCS;                 << 
418           e_sigma_max = e;                     << 
419         }                                      << 
420         if(totCS > 0 && !Emin_found)           << 
421         {                                      << 
422           fEminForAdjSigmaTables[i].push_back( << 
423           Emin_found = true;                   << 
424         }                                      << 
425       }                                        << 
426       fEkinofAdjSigmaMax[i].push_back(e_sigma_ << 
427       if(!Emin_found)                          << 
428         fEminForAdjSigmaTables[i].push_back(fT << 
429                                                   269 
430       fTotalAdjSigmaTable[i]->push_back(aVecto << 270   
431     }                                          << 271   const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 272  
                                                   >> 273  
                                                   >> 274  //Prepare the Sigma table for all AdjointEMModel, will be filled later on 
                                                   >> 275   for (size_t i=0; i<listOfAdjointEMModel.size();i++){
                                                   >> 276     listSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy();
                                                   >> 277   listSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy();
                                                   >> 278   for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
                                                   >> 279     listSigmaTableForAdjointModelScatProjToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
                                                   >> 280     listSigmaTableForAdjointModelProdToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
                                                   >> 281   }
                                                   >> 282   }
                                                   >> 283   
                                                   >> 284 
                                                   >> 285 
                                                   >> 286   for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
                                                   >> 287     G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
                                                   >> 288   DefineCurrentParticle(thePartDef);
                                                   >> 289     theTotalForwardSigmaTableVector[i]->clearAndDestroy();
                                                   >> 290   theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
                                                   >> 291   EminForFwdSigmaTables[i].clear();
                                                   >> 292   EminForAdjSigmaTables[i].clear();
                                                   >> 293   EkinofFwdSigmaMax[i].clear();
                                                   >> 294   EkinofAdjSigmaMax[i].clear();
                                                   >> 295   //G4cout<<thePartDef->GetParticleName();
                                                   >> 296   
                                                   >> 297     for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
                                                   >> 298     const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
                                                   >> 299     
                                                   >> 300     /*
                                                   >> 301     G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
                                                   >> 302     G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
                                                   >> 303     
                                                   >> 304     std::fstream FileOutputAdjCS(file_name1, std::ios::out);
                                                   >> 305     std::fstream FileOutputFwdCS(file_name2, std::ios::out);
                                                   >> 306     
                                                   >> 307     
                                                   >> 308     
                                                   >> 309     FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
                                                   >> 310     FileOutputAdjCS<<std::setprecision(6);
                                                   >> 311     FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
                                                   >> 312     FileOutputFwdCS<<std::setprecision(6);  
                                                   >> 313     */
                                                   >> 314         
                                                   >> 315 
                                                   >> 316     //make first the total fwd CS table for FwdProcess
                                                   >> 317     G4PhysicsVector* aVector =  new G4PhysicsLogVector(Tmin, Tmax, nbins);
                                                   >> 318     G4bool Emin_found=false;
                                                   >> 319     G4double sigma_max =0.;
                                                   >> 320     G4double e_sigma_max =0.;
                                                   >> 321     for(size_t l=0; l<aVector->GetVectorLength(); l++) { 
                                                   >> 322       G4double totCS=0.;
                                                   >> 323       G4double e=aVector->GetLowEdgeEnergy(l);
                                                   >> 324       for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
                                                   >> 325         totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
                                                   >> 326       }
                                                   >> 327       for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
                                                   >> 328         if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
                                                   >> 329                 size_t mat_index = couple->GetIndex();
                                                   >> 330           G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
                                                   >> 331             G4double chargeSqRatio =  currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
                                                   >> 332           (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
                                                   >> 333         }
                                                   >> 334         G4double e1=e/massRatio;
                                                   >> 335         totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
                                                   >> 336       }
                                                   >> 337       aVector->PutValue(l,totCS);
                                                   >> 338       if (totCS>sigma_max){
                                                   >> 339         sigma_max=totCS;
                                                   >> 340         e_sigma_max = e;
                                                   >> 341         
                                                   >> 342       }
                                                   >> 343       //FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
                                                   >> 344       
                                                   >> 345       if (totCS>0 && !Emin_found) {
                                                   >> 346         EminForFwdSigmaTables[i].push_back(e);
                                                   >> 347         Emin_found=true;
                                                   >> 348       }
                                                   >> 349       
                                                   >> 350     
                                                   >> 351     }
                                                   >> 352     //FileOutputFwdCS.close();
                                                   >> 353     
                                                   >> 354     EkinofFwdSigmaMax[i].push_back(e_sigma_max);
                                                   >> 355     
                                                   >> 356     
                                                   >> 357     if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
                                                   >> 358     
                                                   >> 359     theTotalForwardSigmaTableVector[i]->push_back(aVector);
                                                   >> 360     
                                                   >> 361     
                                                   >> 362     Emin_found=false;
                                                   >> 363     sigma_max=0;
                                                   >> 364     e_sigma_max =0.;
                                                   >> 365     G4PhysicsVector* aVector1 =  new G4PhysicsLogVector(Tmin, Tmax, nbins);
                                                   >> 366     for(eindex=0; eindex<aVector->GetVectorLength(); eindex++) { 
                                                   >> 367       G4double e=aVector->GetLowEdgeEnergy(eindex);
                                                   >> 368       G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
                                                   >> 369       aVector1->PutValue(eindex,totCS);
                                                   >> 370       if (totCS>sigma_max){
                                                   >> 371         sigma_max=totCS;
                                                   >> 372         e_sigma_max = e;
                                                   >> 373         
                                                   >> 374       }
                                                   >> 375       //FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
                                                   >> 376       if (totCS>0 && !Emin_found) {
                                                   >> 377         EminForAdjSigmaTables[i].push_back(e);
                                                   >> 378         Emin_found=true;
                                                   >> 379       }
                                                   >> 380     
                                                   >> 381     }
                                                   >> 382     //FileOutputAdjCS.close();
                                                   >> 383     EkinofAdjSigmaMax[i].push_back(e_sigma_max);
                                                   >> 384     if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
                                                   >> 385       
                                                   >> 386     theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
                                                   >> 387     
                                                   >> 388   }
432   }                                               389   }
433   fSigmaTableBuilt = true;                     << 390   TotalSigmaTableAreBuilt =true;
434 }                                              << 391    
435                                                << 
436 ////////////////////////////////////////////// << 
437 G4double G4AdjointCSManager::GetTotalAdjointCS << 
438   G4ParticleDefinition* aPartDef, G4double Eki << 
439   const G4MaterialCutsCouple* aCouple)         << 
440 {                                              << 
441   DefineCurrentMaterial(aCouple);              << 
442   DefineCurrentParticle(aPartDef);             << 
443   return (((*fTotalAdjSigmaTable[fCurrentParti << 
444             ->Value(Ekin * fMassRatio));       << 
445 }                                                 392 }
446                                                << 
447 //////////////////////////////////////////////    393 ///////////////////////////////////////////////////////
448 G4double G4AdjointCSManager::GetTotalForwardCS << 394 //
449   G4ParticleDefinition* aPartDef, G4double Eki << 395 G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
450   const G4MaterialCutsCouple* aCouple)         << 396                const G4MaterialCutsCouple* aCouple)
451 {                                              << 397 { DefineCurrentMaterial(aCouple);
452   DefineCurrentMaterial(aCouple);              << 398   DefineCurrentParticle(aPartDef);  
                                                   >> 399   G4bool b;
                                                   >> 400   return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
                                                   >> 401   
                                                   >> 402   
                                                   >> 403   
                                                   >> 404 }            
                                                   >> 405 ///////////////////////////////////////////////////////
                                                   >> 406 //             
                                                   >> 407 G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
                                                   >> 408                const G4MaterialCutsCouple* aCouple)
                                                   >> 409 { DefineCurrentMaterial(aCouple);
453   DefineCurrentParticle(aPartDef);                410   DefineCurrentParticle(aPartDef);
454   return (((*fTotalFwdSigmaTable[fCurrentParti << 411   G4bool b;
455             ->Value(Ekin * fMassRatio));       << 412   return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
                                                   >> 413   
                                                   >> 414   
456 }                                                 415 }
457                                                << 
458 //////////////////////////////////////////////    416 ///////////////////////////////////////////////////////
459 G4double G4AdjointCSManager::GetAdjointSigma(  << 417 //
460   G4double Ekin_nuc, std::size_t index_model,  << 418 G4double G4AdjointCSManager::GetAdjointSigma(G4double Ekin_nuc, size_t index_model,G4bool is_scat_proj_to_proj,
461   const G4MaterialCutsCouple* aCouple)         << 419                const G4MaterialCutsCouple* aCouple)
462 {                                              << 420 { DefineCurrentMaterial(aCouple);
463   DefineCurrentMaterial(aCouple);              << 421   G4bool b;
464   if(is_scat_proj_to_proj)                     << 422   if (is_scat_proj_to_proj) return (((*listSigmaTableForAdjointModelScatProjToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
465     return (((*fSigmaTableForAdjointModelScatP << 423   else return (((*listSigmaTableForAdjointModelProdToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
466                [fCurrentMatIndex])->Value(Ekin << 424 }                        
467   else                                         << 
468     return (                                   << 
469       ((*fSigmaTableForAdjointModelProdToProj[ << 
470         ->Value(Ekin_nuc));                    << 
471 }                                              << 
472                                                << 
473 //////////////////////////////////////////////    425 ///////////////////////////////////////////////////////
                                                   >> 426 //             
474 void G4AdjointCSManager::GetEminForTotalCS(G4P    427 void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
475                                            con << 428                const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
476                                            G4d << 429 { DefineCurrentMaterial(aCouple);
477                                            G4d << 
478 {                                              << 
479   DefineCurrentMaterial(aCouple);              << 
480   DefineCurrentParticle(aPartDef);                430   DefineCurrentParticle(aPartDef);
481   emin_adj = fEminForAdjSigmaTables[fCurrentPa << 431   emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
482              fMassRatio;                       << 432   emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
483   emin_fwd = fEminForFwdSigmaTables[fCurrentPa << 433   
484              fMassRatio;                       << 434   
                                                   >> 435   
485 }                                                 436 }
486                                                << 
487 //////////////////////////////////////////////    437 ///////////////////////////////////////////////////////
                                                   >> 438 //             
488 void G4AdjointCSManager::GetMaxFwdTotalCS(G4Pa    439 void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
489                                           cons << 440                const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
490                                           G4do << 441 { DefineCurrentMaterial(aCouple);
491                                           G4do << 
492 {                                              << 
493   DefineCurrentMaterial(aCouple);              << 
494   DefineCurrentParticle(aPartDef);                442   DefineCurrentParticle(aPartDef);
495   e_sigma_max = fEkinofFwdSigmaMax[fCurrentPar << 443   e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
496   sigma_max = ((*fTotalFwdSigmaTable[fCurrentP << 444   G4bool b;
497                 ->Value(e_sigma_max);          << 445   sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
498   e_sigma_max /= fMassRatio;                   << 446   e_sigma_max/=massRatio;
                                                   >> 447   
                                                   >> 448   
499 }                                                 449 }
500                                                << 
501 //////////////////////////////////////////////    450 ///////////////////////////////////////////////////////
                                                   >> 451 //             
502 void G4AdjointCSManager::GetMaxAdjTotalCS(G4Pa    452 void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
503                                           cons << 453                const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
504                                           G4do << 454 { DefineCurrentMaterial(aCouple);
505                                           G4do << 
506 {                                              << 
507   DefineCurrentMaterial(aCouple);              << 
508   DefineCurrentParticle(aPartDef);                455   DefineCurrentParticle(aPartDef);
509   e_sigma_max = fEkinofAdjSigmaMax[fCurrentPar << 456   e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
510   sigma_max = ((*fTotalAdjSigmaTable[fCurrentP << 457   G4bool b;
511                 ->Value(e_sigma_max);          << 458   sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
512   e_sigma_max /= fMassRatio;                   << 459   e_sigma_max/=massRatio;
513 }                                              << 460   
514                                                << 461   
515 ////////////////////////////////////////////// << 462 }            
516 G4double G4AdjointCSManager::GetCrossSectionCo << 463 ///////////////////////////////////////////////////////
517   G4ParticleDefinition* aPartDef, G4double Pre << 464 //             
518   const G4MaterialCutsCouple* aCouple, G4bool& << 465 G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
519 {                                              << 466                     G4double& fwd_TotCS)
520   static G4double lastEkin = 0.;               << 467 { G4double corr_fac = 1.;
521   static G4ParticleDefinition* lastPartDef;    << 468   if (forward_CS_mode) {
522                                                << 469     fwd_TotCS=PrefwdCS;
523   G4double corr_fac = 1.;                      << 470     if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
524   if(fForwardCSMode && aPartDef)               << 471     DefineCurrentMaterial(aCouple);
525   {                                            << 472     PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
526     if(lastEkin != PreStepEkin || aPartDef !=  << 473     PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
527        aCouple != fCurrentCouple)              << 474     LastEkinForCS = PreStepEkin;
528     {                                          << 475     lastPartDefForCS = aPartDef;
529       DefineCurrentMaterial(aCouple);          << 476     if (PrefwdCS >0. &&  PreadjCS >0.) {
530       G4double preadjCS = GetTotalAdjointCS(aP << 477       forward_CS_is_used = true;
531       G4double prefwdCS = GetTotalForwardCS(aP << 478       LastCSCorrectionFactor = PrefwdCS/PreadjCS;
532       lastEkin          = PreStepEkin;         << 479     }
533       lastPartDef       = aPartDef;            << 480     else {
534       if(prefwdCS > 0. && preadjCS > 0.)       << 481       forward_CS_is_used = false;
535       {                                        << 482       LastCSCorrectionFactor = 1.;
536         fForwardCSUsed          = true;        << 483       
537         fLastCSCorrectionFactor = prefwdCS / p << 484     }
538       }                                        << 485     
539       else                                     << 486   }
540       {                                        << 487   corr_fac =LastCSCorrectionFactor;
541         fForwardCSUsed          = false;       << 488   
542         fLastCSCorrectionFactor = 1.;          << 489   
543       }                                        << 490   
544     }                                          << 491   }  
545     corr_fac = fLastCSCorrectionFactor;        << 492   else {
546   }                                            << 493     forward_CS_is_used = false;
547   else                                         << 494     LastCSCorrectionFactor = 1.;
548   {                                            << 495   }
549     fForwardCSUsed          = false;           << 496   fwd_TotCS=PrefwdCS; 
550     fLastCSCorrectionFactor = 1.;              << 497   fwd_is_used = forward_CS_is_used;
551   }                                            << 498   return  corr_fac;
552   fwd_is_used = fForwardCSUsed;                << 499 }            
553   return corr_fac;                             << 500 ///////////////////////////////////////////////////////
554 }                                              << 501 //                       
                                                   >> 502 G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
                                                   >> 503                const G4MaterialCutsCouple* aCouple, G4double step_length)
                                                   >> 504 {  G4double corr_fac = 1.;
                                                   >> 505   //return corr_fac;
                                                   >> 506   //G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
                                                   >> 507   G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
                                                   >> 508   G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
                                                   >> 509   if (!forward_CS_is_used || pre_adjCS ==0. ||  after_fwdCS==0.) {
                                                   >> 510   forward_CS_is_used=false;
                                                   >> 511   G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
                                                   >> 512   corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
                                                   >> 513   LastCSCorrectionFactor = 1.;
                                                   >> 514   }
                                                   >> 515   else {
                                                   >> 516   LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
                                                   >> 517   } 
                                                   >> 518   
555                                                   519 
                                                   >> 520  
                                                   >> 521   return corr_fac; 
                                                   >> 522 }            
556 //////////////////////////////////////////////    523 ///////////////////////////////////////////////////////
557 G4double G4AdjointCSManager::GetContinuousWeig << 524 //                       
558   G4ParticleDefinition* aPartDef, G4double Pre << 525 G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
559   const G4MaterialCutsCouple* aCouple, G4doubl << 526 {//return 1.; 
560 {                                              << 527  return  1./LastCSCorrectionFactor;
561   G4double corr_fac    = 1.;                   << 528   
562   G4double after_fwdCS = GetTotalForwardCS(aPa << 529 }               
563   G4double pre_adjCS   = GetTotalAdjointCS(aPa << 
564   if(!fForwardCSUsed || pre_adjCS == 0. || aft << 
565   {                                            << 
566     G4double pre_fwdCS = GetTotalForwardCS(aPa << 
567     corr_fac *= std::exp((pre_adjCS - pre_fwdC << 
568     fLastCSCorrectionFactor = 1.;              << 
569   }                                            << 
570   else                                         << 
571   {                                            << 
572     fLastCSCorrectionFactor = after_fwdCS / pr << 
573   }                                            << 
574   return corr_fac;                             << 
575 }                                              << 
576                                                << 
577 //////////////////////////////////////////////    530 ///////////////////////////////////////////////////////
578 G4double G4AdjointCSManager::GetPostStepWeight << 531 //
579 {                                              << 532 G4double  G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
580   return 1. / fLastCSCorrectionFactor;         << 533                    G4VEmAdjointModel* aModel, 
581 }                                              << 534                    G4double PrimEnergy,
582                                                << 535                    G4double Tcut,
583 ////////////////////////////////////////////// << 536                    G4bool IsScatProjToProjCase,
584 G4double G4AdjointCSManager::ComputeAdjointCS( << 537                std::vector<G4double>& CS_Vs_Element)
585   G4Material* aMaterial, G4VEmAdjointModel* aM << 538 { 
586   G4double Tcut, G4bool isScatProjToProj, std: << 539   
587 {                                              << 540   G4double EminSec=0;
588   G4double EminSec = 0.;                       << 541   G4double EmaxSec=0;
589   G4double EmaxSec = 0.;                       << 542   
590                                                << 543   if (IsScatProjToProjCase){
591   static G4double lastPrimaryEnergy = 0.;      << 544   EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
592   static G4double lastTcut          = 0.;      << 545   EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
593   static G4Material* lastMaterial   = nullptr; << 546   }
594                                                << 547   else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
595   if(isScatProjToProj)                         << 548   EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
596   {                                            << 549   EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
597     EminSec = aModel->GetSecondAdjEnergyMinFor << 550   }
598     EmaxSec = aModel->GetSecondAdjEnergyMaxFor << 551   if (EminSec >= EmaxSec) return 0.;
599   }                                            << 552 
600   else if(PrimEnergy > Tcut || !aModel->GetApp << 553 
601   {                                            << 554   G4bool need_to_compute=false;
602     EminSec = aModel->GetSecondAdjEnergyMinFor << 555   if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
603     EmaxSec = aModel->GetSecondAdjEnergyMaxFor << 556     lastMaterial =aMaterial;
604   }                                            << 557     lastPrimaryEnergy = PrimEnergy;
605   if(EminSec >= EmaxSec)                       << 558     lastTcut=Tcut;
606     return 0.;                                 << 559   listOfIndexOfAdjointEMModelInAction.clear();
607                                                << 560   listOfIsScatProjToProjCase.clear();
608   G4bool need_to_compute = false;              << 561   lastAdjointCSVsModelsAndElements.clear();
609   if(aMaterial != lastMaterial || PrimEnergy ! << 562   need_to_compute=true;
610      Tcut != lastTcut)                         << 563   
611   {                                            << 564   }
612     lastMaterial      = aMaterial;             << 565   size_t ind=0;
613     lastPrimaryEnergy = PrimEnergy;            << 566   if (!need_to_compute){
614     lastTcut          = Tcut;                  << 567     need_to_compute=true;
615     fIndexOfAdjointEMModelInAction.clear();    << 568     for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){
616     fIsScatProjToProj.clear();                 << 569     size_t ind1=listOfIndexOfAdjointEMModelInAction[i];
617     fLastAdjointCSVsModelsAndElements.clear(); << 570     if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){
618     need_to_compute = true;                    << 571       need_to_compute=false;
619   }                                            << 572       CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind];
620                                                << 573     }
621   std::size_t ind = 0;                         << 574     ind++;
622   if(!need_to_compute)                         << 575   }
623   {                                            << 576   }
624     need_to_compute = true;                    << 577   
625     for(std::size_t i = 0; i < fIndexOfAdjoint << 578   if (need_to_compute){
626     {                                          << 579     size_t ind_model=0;
627       std::size_t ind1 = fIndexOfAdjointEMMode << 580   for (size_t i=0;i<listOfAdjointEMModel.size();i++){
628       if(aModel == fAdjointModels[ind1] &&     << 581     if (aModel == listOfAdjointEMModel[i]){
629          isScatProjToProj == fIsScatProjToProj << 582       ind_model=i;
630       {                                        << 583       break;
631         need_to_compute = false;               << 584     }
632         CS_Vs_Element   = fLastAdjointCSVsMode << 585   }
633       }                                        << 586   G4double Tlow=Tcut;
634       ++ind;                                   << 587   if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit();
635     }                                          << 588   listOfIndexOfAdjointEMModelInAction.push_back(ind_model); 
636   }                                            << 589   listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
637                                                << 590   CS_Vs_Element.clear();
638   if(need_to_compute)                          << 591   if (!aModel->GetUseMatrix()){
639   {                                            << 592     CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
640     std::size_t ind_model = 0;                 << 593                  
641     for(std::size_t i = 0; i < fAdjointModels. << 594   
642     {                                          << 595   }
643       if(aModel == fAdjointModels[i])          << 596   else if (aModel->GetUseMatrixPerElement()){
644       {                                        << 597       size_t n_el = aMaterial->GetNumberOfElements();
645         ind_model = i;                         << 598     if (aModel->GetUseOnlyOneMatrixForAllElements()){
646         break;                                 << 599       G4AdjointCSMatrix* theCSMatrix;
647       }                                        << 600       if (IsScatProjToProjCase){
648     }                                          << 601         theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0];
649     G4double Tlow = Tcut;                      << 602       }
650     if(!fAdjointModels[ind_model]->GetApplyCut << 603       else    theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0];
651       Tlow = fAdjointModels[ind_model]->GetLow << 604       G4double CS =0.;
652     fIndexOfAdjointEMModelInAction.push_back(i << 605       if (PrimEnergy > Tlow)
653     fIsScatProjToProj.push_back(isScatProjToPr << 606           CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
654     CS_Vs_Element.clear();                     << 607       G4double factor=0.;
655     if(!aModel->GetUseMatrix())                << 608       for (size_t i=0;i<n_el;i++){ //this could be computed only once
656     {                                          << 609         //size_t ind_el = aMaterial->GetElement(i)->GetIndex();
657       CS_Vs_Element.push_back(aModel->AdjointC << 610         factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
658         fCurrentCouple, PrimEnergy, isScatProj << 611       }
659     }                                          << 612       CS *=factor;
660     else if(aModel->GetUseMatrixPerElement())  << 613       CS_Vs_Element.push_back(CS);
661     {                                          << 614                 
662       std::size_t n_el = aMaterial->GetNumberO << 615     }
663       if(aModel->GetUseOnlyOneMatrixForAllElem << 616     else {
664       {                                        << 617       for (size_t i=0;i<n_el;i++){
665         G4AdjointCSMatrix* theCSMatrix;        << 618         size_t ind_el = aMaterial->GetElement(i)->GetIndex();
666         if(isScatProjToProj)                   << 619         //G4cout<<aMaterial->GetName()<<G4endl;
667         {                                      << 620         G4AdjointCSMatrix* theCSMatrix;
668           theCSMatrix = fAdjointCSMatricesForS << 621         if (IsScatProjToProjCase){
669         }                                      << 622           theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
670         else                                   << 623         }
671           theCSMatrix = fAdjointCSMatricesForP << 624         else    theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
672         G4double CS = 0.;                      << 625         G4double CS =0.;
673         if(PrimEnergy > Tlow)                  << 626         if (PrimEnergy > Tlow)
674           CS = ComputeAdjointCS(PrimEnergy, th << 627           CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
675         G4double factor = 0.;                  << 628         //G4cout<<CS<<G4endl;     
676         for(G4int i = 0; i < (G4int)n_el; ++i) << 629         CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i])); 
677         {  // this could be computed only once << 630       }
678           factor += aMaterial->GetElement(i)-> << 631     }   
679                     aMaterial->GetVecNbOfAtoms << 632     
680         }                                      << 633   }
681         CS *= factor;                          << 634   else {
682         CS_Vs_Element.push_back(CS);           << 635     size_t ind_mat = aMaterial->GetIndex();
683       }                                        << 636     G4AdjointCSMatrix* theCSMatrix;
684       else                                     << 637     if (IsScatProjToProjCase){
685       {                                        << 638       theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
686         for(G4int i = 0; i < (G4int)n_el; ++i) << 639     }
687         {                                      << 640     else    theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat];
688           std::size_t ind_el = aMaterial->GetE << 641     G4double CS =0.;
689           G4AdjointCSMatrix* theCSMatrix;      << 642     if (PrimEnergy > Tlow)
690           if(isScatProjToProj)                 << 643       CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
691           {                                    << 644     CS_Vs_Element.push_back(CS);              
692             theCSMatrix =                      << 645       
693               fAdjointCSMatricesForScatProjToP << 646     
694           }                                    << 647   }
695           else                                 << 648   lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
696             theCSMatrix = fAdjointCSMatricesFo << 649   
697           G4double CS = 0.;                    << 650   }
698           if(PrimEnergy > Tlow)                << 651   
699             CS = ComputeAdjointCS(PrimEnergy,  << 652   
700           CS_Vs_Element.push_back(CS *         << 653   G4double CS=0;
701                                   (aMaterial-> << 654   for (size_t i=0;i<CS_Vs_Element.size();i++){
702         }                                      << 655     CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
703       }                                        << 656     
704     }                                          << 
705     else                                       << 
706     {                                          << 
707       std::size_t ind_mat = aMaterial->GetInde << 
708       G4AdjointCSMatrix* theCSMatrix;          << 
709       if(isScatProjToProj)                     << 
710       {                                        << 
711         theCSMatrix = fAdjointCSMatricesForSca << 
712       }                                        << 
713       else                                     << 
714         theCSMatrix = fAdjointCSMatricesForPro << 
715       G4double CS = 0.;                        << 
716       if(PrimEnergy > Tlow)                    << 
717         CS = ComputeAdjointCS(PrimEnergy, theC << 
718       CS_Vs_Element.push_back(CS);             << 
719     }                                          << 
720     fLastAdjointCSVsModelsAndElements.push_bac << 
721   }                                            << 
722                                                << 
723   G4double CS = 0.;                            << 
724   for(const auto& cs_vs_el : CS_Vs_Element)    << 
725   {                                            << 
726     // We could put the progressive sum of the << 
727     // element itself                          << 
728     CS += cs_vs_el;                            << 
729   }                                               657   }
730   return CS;                                      658   return CS;
731 }                                              << 659 }             
732                                                << 
733 //////////////////////////////////////////////    660 ///////////////////////////////////////////////////////
734 G4Element* G4AdjointCSManager::SampleElementFr << 661 //  
735   G4Material* aMaterial, G4VEmAdjointModel* aM << 662 G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial,
736   G4double Tcut, G4bool isScatProjToProj)      << 663                          G4VEmAdjointModel* aModel,
737 {                                              << 664                          G4double PrimEnergy,
738   std::vector<G4double> CS_Vs_Element;         << 665                          G4double Tcut,
739   G4double CS    = ComputeAdjointCS(aMaterial, << 666                          G4bool IsScatProjToProjCase)
740                                  isScatProjToP << 667 { std::vector<G4double> CS_Vs_Element;
741   G4double SumCS = 0.;                         << 668   G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
742   std::size_t ind     = 0;                     << 669   G4double rand_var= G4UniformRand();
743   for(std::size_t i = 0; i < CS_Vs_Element.siz << 670   G4double SumCS=0.;
744   {                                            << 671   size_t ind=0;
745     SumCS += CS_Vs_Element[i];                 << 672   for (size_t i=0;i<CS_Vs_Element.size();i++){
746     if(G4UniformRand() <= SumCS / CS)          << 673   SumCS+=CS_Vs_Element[i];
747     {                                          << 674   if (rand_var<=SumCS/CS){
748       ind = i;                                 << 675     ind=i;
749       break;                                   << 676     break;
750     }                                          << 677   }
751   }                                               678   }
752                                                << 679  
753   return const_cast<G4Element*>(aMaterial->Get << 680   return const_cast<G4Element*>(aMaterial->GetElement(ind));
754 }                                              << 681  
755                                                << 682  
                                                   >> 683               
                                                   >> 684 }               
756 //////////////////////////////////////////////    685 ///////////////////////////////////////////////////////
757 G4double G4AdjointCSManager::ComputeTotalAdjoi << 686 //
758   const G4MaterialCutsCouple* aCouple, G4Parti << 687 G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple,
759   G4double Ekin)                               << 688                    G4ParticleDefinition* aPartDef,
760 {                                              << 689                    G4double Ekin)
761   G4double TotalCS = 0.;                       << 690 {
762                                                << 691  G4double TotalCS=0.;
763   DefineCurrentMaterial(aCouple);              << 692 
764                                                << 693  DefineCurrentMaterial(aCouple);
765   std::vector<G4double> CS_Vs_Element;         << 694 
766   G4double CS;                                 << 695   
767   G4VEmAdjointModel* adjModel = nullptr;       << 696  std::vector<G4double> CS_Vs_Element;
768   for(std::size_t i = 0; i < fAdjointModels.si << 697  G4double CS;
769   {                                            << 698  for (size_t i=0; i<listOfAdjointEMModel.size();i++){
770     G4double Tlow = 0.;                        << 699   
771     adjModel      = fAdjointModels[i];         << 700   G4double Tlow=0;
772     if(!adjModel->GetApplyCutInRange())        << 701   if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
773       Tlow = adjModel->GetLowEnergyLimit();    << 702   else {
774     else                                       << 703     G4ParticleDefinition* theDirSecondPartDef = 
775     {                                          << 704       GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
776       G4ParticleDefinition* theDirSecondPartDe << 705     size_t idx=56;
777         adjModel->GetAdjointEquivalentOfDirect << 706     if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
778       std::size_t idx = 56;                    << 707     else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
779       if(theDirSecondPartDef->GetParticleName( << 708     else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
780         idx = 0;                               << 709     if (idx <56) {
781       else if(theDirSecondPartDef->GetParticle << 710       const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
782         idx = 1;                               << 711       Tlow =(*aVec)[aCouple->GetIndex()];
783       else if(theDirSecondPartDef->GetParticle << 712     } 
784         idx = 2;                               << 713     
785       if(idx < 56)                             << 714   
786       {                                        << 715   }
787         const std::vector<G4double>* aVec =    << 716   if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
788           G4ProductionCutsTable::GetProduction << 717     if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
789             idx);                              << 718       CS=ComputeAdjointCS(currentMaterial,
790         Tlow = (*aVec)[aCouple->GetIndex()];   << 719                        listOfAdjointEMModel[i], 
791       }                                        << 720                        Ekin, Tlow,true,CS_Vs_Element);
792     }                                          << 721       TotalCS += CS;
793     if(Ekin <= adjModel->GetHighEnergyLimit()  << 722       (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,CS);             
794        Ekin >= adjModel->GetLowEnergyLimit())  << 723     }
795     {                                          << 724     if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
796       if(aPartDef ==                           << 725       CS = ComputeAdjointCS(currentMaterial,
797          adjModel->GetAdjointEquivalentOfDirec << 726                        listOfAdjointEMModel[i], 
798       {                                        << 727                        Ekin, Tlow,false, CS_Vs_Element);
799         CS = ComputeAdjointCS(fCurrentMaterial << 728       TotalCS += CS;
800                               CS_Vs_Element);  << 729       (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,CS);
801         TotalCS += CS;                         << 730     }
802         (*fSigmaTableForAdjointModelScatProjTo << 731     
803           ->PutValue(fNbins, CS);              << 732   }
804       }                                        << 733   else {
805       if(aPartDef ==                           << 734     (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,0.);
806          adjModel->GetAdjointEquivalentOfDirec << 735     (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,0.);
807       {                                        << 736     
808         CS = ComputeAdjointCS(fCurrentMaterial << 737   }
809                               CS_Vs_Element);  << 738  }
810         TotalCS += CS;                         << 739  return TotalCS;
811         (*fSigmaTableForAdjointModelProdToProj << 740     
812           fNbins, CS);                         << 741  
813       }                                        << 742 } 
814     }                                          << 
815     else                                       << 
816     {                                          << 
817       (*fSigmaTableForAdjointModelScatProjToPr << 
818         ->PutValue(fNbins, 0.);                << 
819       (*fSigmaTableForAdjointModelProdToProj[i << 
820         fNbins, 0.);                           << 
821     }                                          << 
822   }                                            << 
823   return TotalCS;                              << 
824 }                                              << 
825                                                << 
826 //////////////////////////////////////////////    743 ///////////////////////////////////////////////////////
                                                   >> 744 //
827 std::vector<G4AdjointCSMatrix*>                   745 std::vector<G4AdjointCSMatrix*>
828 G4AdjointCSManager::BuildCrossSectionsModelAnd << 746 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
829                                                << 747                         G4int nbin_pro_decade)
830                                                << 748 { 
831 {                                              << 749   G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
832   G4AdjointCSMatrix* theCSMatForProdToProjBack << 750   G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
833     new G4AdjointCSMatrix(false);              << 751  
834   G4AdjointCSMatrix* theCSMatForScatProjToProj << 752  
835     new G4AdjointCSMatrix(true);               << 753   //make the vector of primary energy of the adjoint particle, could try to make this just once ? 
836                                                << 754   
837   // make the vector of primary energy of the  << 755    G4double EkinMin =aModel->GetLowEnergyLimit();
838   G4double EkinMin        = aModel->GetLowEner << 756    G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
839   G4double EkinMaxForScat = aModel->GetHighEne << 757    G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
840   G4double EkinMaxForProd = aModel->GetHighEne << 758    if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
841   if(aModel->GetSecondPartOfSameType())        << 759    
842     EkinMaxForProd = EkinMaxForProd / 2.;      << 760     
843                                                << 761    //Product to projectile backward scattering
844   // Product to projectile backward scattering << 762    //-----------------------------------------
845   G4double dE = std::pow(10., 1. / nbin_pro_de << 763    G4double fE=std::pow(10.,1./nbin_pro_decade);
846   G4double E2 =                                << 764    G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
847     std::pow(10., double(int(std::log10(EkinMi << 765    G4double E1=EkinMin;
848                     nbin_pro_decade) / dE;     << 766    while (E1 <EkinMaxForProd){
849   G4double E1 = EkinMin;                       << 767     E1=std::max(EkinMin,E2);
850   // Loop checking, 07-Aug-2015, Vladimir Ivan << 768   E1=std::min(EkinMaxForProd,E1);
851   while(E1 < EkinMaxForProd)                   << 769   std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
852   {                                            << 770   if (aMat.size()>=2) {
853     E1 = std::max(EkinMin, E2);                << 771     std::vector< double>* log_ESecVec=aMat[0];
854     E1 = std::min(EkinMaxForProd, E1);         << 772     std::vector< double>* log_CSVec=aMat[1];
855     std::vector<std::vector<double>*> aMat =   << 773     G4double log_adjointCS=log_CSVec->back();
856       aModel->ComputeAdjointCrossSectionVector << 774     //normalise CSVec such that it becomes a probability vector
857                                                << 775     for (size_t j=0;j<log_CSVec->size();j++) {
858     if(aMat.size() >= 2)                       << 776             if (j==0) (*log_CSVec)[j] = 0.; 
859     {                                          << 777       else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
860       std::vector<double>* log_ESecVec = aMat[ << 778     } 
861       std::vector<double>* log_CSVec   = aMat[ << 779     (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
862       G4double log_adjointCS           = log_C << 780     theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
863       // normalise CSVec such that it becomes  << 781   } 
864       for(std::size_t j = 0; j < log_CSVec->si << 782     E1=E2;
865       {                                        << 783     E2*=fE;
866         if(j == 0)                             << 784    }
867           (*log_CSVec)[j] = 0.;                << 785    
868         else                                   << 786    //Scattered projectile to projectile backward scattering
869           (*log_CSVec)[j] =                    << 787    //-----------------------------------------
870             std::log(1. - std::exp((*log_CSVec << 788    
871       }                                        << 789    E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
872       (*log_CSVec)[log_CSVec->size() - 1] =    << 790    E1=EkinMin;
873         (*log_CSVec)[log_CSVec->size() - 2] -  << 791    while (E1 <EkinMaxForScat){
874       theCSMatForProdToProjBackwardScattering- << 792     E1=std::max(EkinMin,E2);
875         std::log(E1), log_adjointCS, log_ESecV << 793   E1=std::min(EkinMaxForScat,E1);
876     }                                          << 794   std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
877     E1 = E2;                                   << 795   if (aMat.size()>=2) {
878     E2 *= dE;                                  << 796     std::vector< double>* log_ESecVec=aMat[0];
879   }                                            << 797     std::vector< double>* log_CSVec=aMat[1];
880                                                << 798     G4double log_adjointCS=log_CSVec->back();
881   // Scattered projectile to projectile backwa << 799     //normalise CSVec such that it becomes a probability vector
882   E2 = std::pow(10., G4double(int(std::log10(E << 800     for (size_t j=0;j<log_CSVec->size();j++) {
883                        nbin_pro_decade) / dE;  << 801             if (j==0) (*log_CSVec)[j] = 0.; 
884   E1 = EkinMin;                                << 802       else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
885   // Loop checking, 07-Aug-2015, Vladimir Ivan << 803     } 
886   while(E1 < EkinMaxForScat)                   << 804     (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
887   {                                            << 805     theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
888     E1 = std::max(EkinMin, E2);                << 806     }
889     E1 = std::min(EkinMaxForScat, E1);         << 807   E1=E2;
890     std::vector<std::vector<G4double>*> aMat = << 808     E2*=fE;
891       aModel->ComputeAdjointCrossSectionVector << 809    }
892         E1, Z, A, nbin_pro_decade);            << 810    
893     if(aMat.size() >= 2)                       << 811    
894     {                                          << 
895       std::vector<G4double>* log_ESecVec = aMa << 
896       std::vector<G4double>* log_CSVec   = aMa << 
897       G4double log_adjointCS             = log << 
898       // normalise CSVec such that it becomes  << 
899       for(std::size_t j = 0; j < log_CSVec->si << 
900       {                                        << 
901         if(j == 0)                             << 
902           (*log_CSVec)[j] = 0.;                << 
903         else                                   << 
904           (*log_CSVec)[j] =                    << 
905             std::log(1. - std::exp((*log_CSVec << 
906       }                                        << 
907       (*log_CSVec)[log_CSVec->size() - 1] =    << 
908         (*log_CSVec)[log_CSVec->size() - 2] -  << 
909       theCSMatForScatProjToProjBackwardScatter << 
910         std::log(E1), log_adjointCS, log_ESecV << 
911     }                                          << 
912     E1 = E2;                                   << 
913     E2 *= dE;                                  << 
914   }                                            << 
915                                                << 
916   std::vector<G4AdjointCSMatrix*> res;            812   std::vector<G4AdjointCSMatrix*> res;
                                                   >> 813   res.clear();
917   res.push_back(theCSMatForProdToProjBackwardS    814   res.push_back(theCSMatForProdToProjBackwardScattering);
918   res.push_back(theCSMatForScatProjToProjBackw    815   res.push_back(theCSMatForScatProjToProjBackwardScattering);
                                                   >> 816   
                                                   >> 817 
                                                   >> 818 /*
                                                   >> 819   G4String file_name;
                                                   >> 820   std::stringstream astream;
                                                   >> 821   G4String str_Z;
                                                   >> 822   astream<<Z;
                                                   >> 823   astream>>str_Z;  
                                                   >> 824   theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt"); 
                                                   >> 825   theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
                                                   >> 826  
                                                   >> 827 */
919                                                   828 
                                                   >> 829   
920   return res;                                     830   return res;
                                                   >> 831   
                                                   >> 832   
921 }                                                 833 }
922                                                << 
923 //////////////////////////////////////////////    834 ///////////////////////////////////////////////////////
                                                   >> 835 //
924 std::vector<G4AdjointCSMatrix*>                   836 std::vector<G4AdjointCSMatrix*>
925 G4AdjointCSManager::BuildCrossSectionsModelAnd << 837 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel,
926   G4VEmAdjointModel* aModel, G4Material* aMate << 838                         G4Material* aMaterial,
927 {                                              << 839                         G4int nbin_pro_decade)
928   G4AdjointCSMatrix* theCSMatForProdToProjBack << 840 { 
929     new G4AdjointCSMatrix(false);              << 841   G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
930   G4AdjointCSMatrix* theCSMatForScatProjToProj << 842   G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
931     new G4AdjointCSMatrix(true);               << 843  
932                                                << 844  
933   G4double EkinMin        = aModel->GetLowEner << 845   //make the vector of primary energy of the adjoint particle, could try to make this just once ? 
934   G4double EkinMaxForScat = aModel->GetHighEne << 846   
935   G4double EkinMaxForProd = aModel->GetHighEne << 847    G4double EkinMin =aModel->GetLowEnergyLimit();
936   if(aModel->GetSecondPartOfSameType())        << 848    G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
937     EkinMaxForProd /= 2.;                      << 849    G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
938                                                << 850    if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
939   // Product to projectile backward scattering << 851    
940   G4double dE = std::pow(10., 1. / nbin_pro_de << 852     
941   G4double E2 =                                << 853    
942     std::pow(10., double(int(std::log10(EkinMi << 854    
943                     nbin_pro_decade) / dE;     << 855    
944   G4double E1 = EkinMin;                       << 856    
945   // Loop checking, 07-Aug-2015, Vladimir Ivan << 857    
946   while(E1 < EkinMaxForProd)                   << 858    //Product to projectile backward scattering
947   {                                            << 859    //-----------------------------------------
948     E1 = std::max(EkinMin, E2);                << 860    G4double fE=std::pow(10.,1./nbin_pro_decade);
949     E1 = std::min(EkinMaxForProd, E1);         << 861    G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
950     std::vector<std::vector<G4double>*> aMat = << 862    G4double E1=EkinMin;
951       aModel->ComputeAdjointCrossSectionVector << 863    while (E1 <EkinMaxForProd){
952         aMaterial, E1, nbin_pro_decade);       << 864     E1=std::max(EkinMin,E2);
953     if(aMat.size() >= 2)                       << 865   E1=std::min(EkinMaxForProd,E1);
954     {                                          << 866   std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
955       std::vector<G4double>* log_ESecVec = aMa << 867   if (aMat.size()>=2) {
956       std::vector<G4double>* log_CSVec   = aMa << 868     std::vector< double>* log_ESecVec=aMat[0];
957       G4double log_adjointCS             = log << 869     std::vector< double>* log_CSVec=aMat[1];
958                                                << 870     G4double log_adjointCS=log_CSVec->back();
959       // normalise CSVec such that it becomes  << 871   
960       for(std::size_t j = 0; j < log_CSVec->si << 872     //normalise CSVec such that it becomes a probability vector
961       {                                        << 873     for (size_t j=0;j<log_CSVec->size();j++) {
962         if(j == 0)                             << 874             //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
963           (*log_CSVec)[j] = 0.;                << 875       if (j==0) (*log_CSVec)[j] = 0.; 
964         else                                   << 876       else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
965           (*log_CSVec)[j] =                    << 877       //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
966             std::log(1. - std::exp((*log_CSVec << 878     } 
967       }                                        << 879     (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
968       (*log_CSVec)[log_CSVec->size() - 1] =    << 880     theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
969         (*log_CSVec)[log_CSVec->size() - 2] -  << 881   } 
970       theCSMatForProdToProjBackwardScattering- << 882   
971         std::log(E1), log_adjointCS, log_ESecV << 883   
972     }                                          << 884   
973                                                << 885     E1=E2;
974     E1 = E2;                                   << 886     E2*=fE;
975     E2 *= dE;                                  << 887    }
976   }                                            << 888    
977                                                << 889    //Scattered projectile to projectile backward scattering
978   // Scattered projectile to projectile backwa << 890    //-----------------------------------------
979   E2 =                                         << 891    
980     std::pow(10., G4double(G4int(std::log10(Ek << 892    E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
981                     nbin_pro_decade) /         << 893    E1=EkinMin;
982     dE;                                        << 894    while (E1 <EkinMaxForScat){
983   E1 = EkinMin;                                << 895     E1=std::max(EkinMin,E2);
984   while(E1 < EkinMaxForScat)                   << 896   E1=std::min(EkinMaxForScat,E1);
985   {                                            << 897   std::vector< std::vector< double>* >  aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
986     E1 = std::max(EkinMin, E2);                << 898   if (aMat.size()>=2) {
987     E1 = std::min(EkinMaxForScat, E1);         << 899     std::vector< double>* log_ESecVec=aMat[0];
988     std::vector<std::vector<G4double>*> aMat = << 900     std::vector< double>* log_CSVec=aMat[1];
989       aModel->ComputeAdjointCrossSectionVector << 901     G4double log_adjointCS=log_CSVec->back();
990         aMaterial, E1, nbin_pro_decade);       << 902   
991     if(aMat.size() >= 2)                       << 903     for (size_t j=0;j<log_CSVec->size();j++) {
992     {                                          << 904             //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
993       std::vector<G4double>* log_ESecVec = aMa << 905       if (j==0) (*log_CSVec)[j] = 0.; 
994       std::vector<G4double>* log_CSVec   = aMa << 906       else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
995       G4double log_adjointCS             = log << 907     //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
996                                                << 908  
997       for(std::size_t j = 0; j < log_CSVec->si << 909     } 
998       {                                        << 910     (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
999         if(j == 0)                             << 911   
1000           (*log_CSVec)[j] = 0.;               << 912     theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
1001         else                                  << 913     }
1002           (*log_CSVec)[j] =                   << 914   E1=E2;
1003             std::log(1. - std::exp((*log_CSVe << 915     E2*=fE; 
1004       }                                       << 916    }
1005       (*log_CSVec)[log_CSVec->size() - 1] =   << 917    
1006         (*log_CSVec)[log_CSVec->size() - 2] - << 918    
1007                                               << 919    
1008       theCSMatForScatProjToProjBackwardScatte << 920    
1009         std::log(E1), log_adjointCS, log_ESec << 921    
1010     }                                         << 922    
1011     E1 = E2;                                  << 923    
1012     E2 *= dE;                                 << 
1013   }                                           << 
1014                                               << 
1015   std::vector<G4AdjointCSMatrix*> res;           924   std::vector<G4AdjointCSMatrix*> res;
                                                   >> 925   res.clear();
                                                   >> 926   
1016   res.push_back(theCSMatForProdToProjBackward    927   res.push_back(theCSMatForProdToProjBackwardScattering);
1017   res.push_back(theCSMatForScatProjToProjBack << 928   res.push_back(theCSMatForScatProjToProjBackwardScattering); 
                                                   >> 929   
                                                   >> 930  /*
                                                   >> 931   theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
                                                   >> 932   theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
                                                   >> 933 */
                                                   >> 934 
1018                                                  935 
1019   return res;                                    936   return res;
                                                   >> 937   
                                                   >> 938   
1020 }                                                939 }
1021                                                  940 
1022 /////////////////////////////////////////////    941 ///////////////////////////////////////////////////////
1023 G4ParticleDefinition* G4AdjointCSManager::Get << 942 //
1024   G4ParticleDefinition* theFwdPartDef)        << 943 G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
1025 {                                                944 {
1026   if(theFwdPartDef->GetParticleName() == "e-" << 945  if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
1027     return G4AdjointElectron::AdjointElectron << 946  else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
1028   else if(theFwdPartDef->GetParticleName() == << 947  else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
1029     return G4AdjointGamma::AdjointGamma();    << 948  else if (theFwdPartDef ==theFwdIon) return theAdjIon;
1030   else if(theFwdPartDef->GetParticleName() == << 949  
1031     return G4AdjointProton::AdjointProton();  << 950  return 0;  
1032   else if(theFwdPartDef == fFwdIon)           << 
1033     return fAdjIon;                           << 
1034   return nullptr;                             << 
1035 }                                                951 }
1036                                               << 
1037 /////////////////////////////////////////////    952 ///////////////////////////////////////////////////////
1038 G4ParticleDefinition* G4AdjointCSManager::Get << 953 //
1039   G4ParticleDefinition* theAdjPartDef)        << 954 G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
1040 {                                                955 {
1041   if(theAdjPartDef->GetParticleName() == "adj << 956  if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
1042     return G4Electron::Electron();            << 957  else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
1043   else if(theAdjPartDef->GetParticleName() == << 958  else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
1044     return G4Gamma::Gamma();                  << 959  else if (theAdjPartDef == theAdjIon) return theFwdIon;
1045   else if(theAdjPartDef->GetParticleName() == << 960  return 0;  
1046     return G4Proton::Proton();                << 
1047   else if(theAdjPartDef == fAdjIon)           << 
1048     return fFwdIon;                           << 
1049   return nullptr;                             << 
1050 }                                                961 }
1051                                               << 
1052 /////////////////////////////////////////////    962 ///////////////////////////////////////////////////////
1053 void G4AdjointCSManager::DefineCurrentMateria << 963 //
1054   const G4MaterialCutsCouple* couple)         << 964 void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
1055 {                                                965 {
1056   if(couple != fCurrentCouple)                << 966   if(couple != currentCouple) {
1057   {                                           << 967     currentCouple   = const_cast<G4MaterialCutsCouple*> (couple);
1058     fCurrentCouple          = const_cast<G4Ma << 968     currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
1059     fCurrentMaterial        = const_cast<G4Ma << 969     currentMatIndex = couple->GetIndex();
1060     fCurrentMatIndex        = couple->GetInde << 970     lastPartDefForCS =0;
1061     fLastCSCorrectionFactor = 1.;             << 971     LastEkinForCS =0;
1062   }                                           << 972     LastCSCorrectionFactor =1.;
                                                   >> 973   }  
1063 }                                                974 }
1064                                                  975 
1065 /////////////////////////////////////////////    976 ///////////////////////////////////////////////////////
1066 void G4AdjointCSManager::DefineCurrentParticl << 977 //
1067   const G4ParticleDefinition* aPartDef)       << 978 void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
1068 {                                                979 {
1069                                               << 980   if(aPartDef != currentParticleDef) {
1070   if(aPartDef != fCurrentParticleDef)         << 981   
1071   {                                           << 982     currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
1072     fCurrentParticleDef = const_cast<G4Partic << 983   massRatio=1;
1073     fMassRatio         = 1.;                  << 984   if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
1074     if(aPartDef == fAdjIon)                   << 985     currentParticleIndex=1000000;
1075       fMassRatio = proton_mass_c2 / aPartDef- << 986   for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
1076     fCurrentParticleIndex = 1000000;          << 987       if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
1077     for(std::size_t i = 0; i < fAdjointPartic << 988     } 
1078     {                                         << 989    
1079       if(aPartDef == fAdjointParticlesInActio << 990   }  
1080         fCurrentParticleIndex = i;            << 
1081     }                                         << 
1082   }                                           << 
1083 }                                                991 }
1084                                                  992 
1085 ///////////////////////////////////////////// << 
1086 G4double G4AdjointCSManager::ComputeAdjointCS << 
1087   G4double aPrimEnergy, G4AdjointCSMatrix* an << 
1088 {                                             << 
1089   std::vector<double>* theLogPrimEnergyVector << 
1090     anAdjointCSMatrix->GetLogPrimEnergyVector << 
1091   if(theLogPrimEnergyVector->empty())         << 
1092   {                                           << 
1093     G4cout << "No data are contained in the g << 
1094     return 0.;                                << 
1095   }                                           << 
1096   G4double log_Tcut = std::log(Tcut);         << 
1097   G4double log_E    = std::log(aPrimEnergy);  << 
1098                                               << 
1099   if(aPrimEnergy <= Tcut || log_E > theLogPri << 
1100     return 0.;                                << 
1101                                                  993 
1102   G4AdjointInterpolator* theInterpolator = G4 << 
1103                                                  994 
1104   std::size_t ind =                           << 995 /////////////////////////////////////////////////////////////////////////////////////////////////
1105     theInterpolator->FindPositionForLogVector << 996 //
1106   G4double aLogPrimEnergy1, aLogPrimEnergy2;  << 997 G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
1107   G4double aLogCS1, aLogCS2;                  << 998         anAdjointCSMatrix,G4double Tcut)
1108   G4double log01, log02;                      << 999 { 
1109   std::vector<G4double>* aLogSecondEnergyVect << 1000   std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
1110   std::vector<G4double>* aLogSecondEnergyVect << 1001   if (theLogPrimEnergyVector->size() ==0){
1111   std::vector<G4double>* aLogProbVector1      << 1002   G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
1112   std::vector<G4double>* aLogProbVector2      << 1003   G4cout<<"The s"<<G4endl;
1113   std::vector<std::size_t>* aLogProbVectorInd << 1004   return 0.;
1114   std::vector<std::size_t>* aLogProbVectorInd << 1005   
1115                                               << 
1116   anAdjointCSMatrix->GetData((G4int)ind, aLog << 
1117                              aLogSecondEnergy << 
1118                              aLogProbVectorIn << 
1119   anAdjointCSMatrix->GetData(G4int(ind + 1),  << 
1120                              aLogSecondEnergy << 
1121                              aLogProbVectorIn << 
1122   if (! (aLogProbVector1 && aLogProbVector2 & << 
1123              aLogSecondEnergyVector1 && aLogS << 
1124      return  0.;                              << 
1125   }                                              1006   }
                                                   >> 1007   G4double log_Tcut = std::log(Tcut);
                                                   >> 1008   G4double log_E =std::log(aPrimEnergy);
                                                   >> 1009   
                                                   >> 1010   if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back()) return 0.;
                                                   >> 1011   
                                                   >> 1012   
1126                                                  1013 
1127   if(anAdjointCSMatrix->IsScatProjToProj())   << 1014   G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
1128   {  // case where the Tcut plays a role      << 1015  
1129     G4double log_minimum_prob1, log_minimum_p << 1016   size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
1130     log_minimum_prob1 = theInterpolator->Inte << 1017   G4double aLogPrimEnergy1,aLogPrimEnergy2;
1131       log_Tcut, *aLogSecondEnergyVector1, *aL << 1018   G4double aLogCS1,aLogCS2;
1132     log_minimum_prob2 = theInterpolator->Inte << 1019   G4double log01,log02;
1133       log_Tcut, *aLogSecondEnergyVector2, *aL << 1020   std::vector< double>* aLogSecondEnergyVector1 =0;
1134     aLogCS1 += log_minimum_prob1;             << 1021   std::vector< double>* aLogSecondEnergyVector2  =0;
1135     aLogCS2 += log_minimum_prob2;             << 1022   std::vector< double>* aLogProbVector1=0;
                                                   >> 1023   std::vector< double>* aLogProbVector2=0;
                                                   >> 1024   std::vector< size_t>* aLogProbVectorIndex1=0;
                                                   >> 1025   std::vector< size_t>* aLogProbVectorIndex2=0;
                                                   >> 1026   
                                                   >> 1027                      
                                                   >> 1028   anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
                                                   >> 1029   anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
                                                   >> 1030   if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
                                                   >> 1031   G4double log_minimum_prob1, log_minimum_prob2;
                                                   >> 1032   log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
                                                   >> 1033   log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
                                                   >> 1034   aLogCS1+= log_minimum_prob1;
                                                   >> 1035   aLogCS2+= log_minimum_prob2;
1136   }                                              1036   }
1137                                               << 1037  
1138   G4double log_adjointCS = theInterpolator->L << 1038   G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2);
1139     log_E, aLogPrimEnergy1, aLogPrimEnergy2,  << 1039   return std::exp(log_adjointCS); 
1140   return std::exp(log_adjointCS);             << 1040   
1141 }                                             << 1041   
                                                   >> 1042 }  
1142                                                  1043