Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNAModelInterface.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/dna/models/src/G4DNAModelInterface.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNAModelInterface.cc (Version 10.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 //                                                 26 //
 27 // Contact authors: S. Meylan, C. Villagrasa       27 // Contact authors: S. Meylan, C. Villagrasa
                                                   >>  28 //
 28 // email: sylvain.meylan@symalgo-tech.com, car     29 // email: sylvain.meylan@symalgo-tech.com, carmen.villagrasa@irsn.fr
 29 // updated : Hoang Tran : 6/1/2023 clean code  << 
 30                                                    30 
 31 #include "G4DNAModelInterface.hh"                  31 #include "G4DNAModelInterface.hh"
 32                                                <<  32 #include "G4PhysicalConstants.hh"
                                                   >>  33 #include "G4SystemOfUnits.hh"
 33 #include "G4DNAMolecularMaterial.hh"               34 #include "G4DNAMolecularMaterial.hh"
 34 #include "G4LossTableManager.hh"               <<  35 
 35 #include "G4ParticleChangeForGamma.hh"         <<  36 
 36 #include "G4VDNAModel.hh"                      <<  37 G4DNAModelInterface::G4DNAModelInterface(const G4String &nam)
 37 #include "G4VEmModel.hh"                       <<  38     : G4VEmModel(nam), fName(nam), fpParticleChangeForGamma(0), fSampledMat("")
 38 G4DNAModelInterface::G4DNAModelInterface(const <<  39 {
                                                   >>  40 
                                                   >>  41 }
 39                                                    42 
 40 //....oooOO0OOooo........oooOO0OOooo........oo     43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 41                                                    44 
 42 void G4DNAModelInterface::Initialise(const G4P <<  45 G4DNAModelInterface::~G4DNAModelInterface()
 43 {                                                  46 {
 44   // Those two statements are necessary to ove <<  47     // Loop on all the registered models to properly delete them (free the memory)
 45   // (ionisation, elastic, etc...). Indeed, wi <<  48     for(unsigned int i=0, ie = fRegisteredModels.size(); i<ie; ++i)
 46   // themselves their energy limits per materi <<  49     {
 47   // in the G4DNAProcess classes.              <<  50             if(fRegisteredModels.at(i) != nullptr) delete fRegisteredModels.at(i);
 48   //                                           <<  51     }
                                                   >>  52 }
 49                                                    53 
 50   fpG4_WATER = G4Material::GetMaterial("G4_WAT <<  54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 51                                                    55 
 52   SetLowEnergyLimit(0.);                       <<  56 void G4DNAModelInterface::Initialise(const G4ParticleDefinition* particle,
 53   SetHighEnergyLimit(DBL_MAX);                 <<  57                                   const G4DataVector& cuts)
                                                   >>  58 {
                                                   >>  59     // Those two statements are necessary to override the energy limits set in the G4DNAProcesses (ionisation, elastic, etc...).
                                                   >>  60     // Indeed, with the ModelInterface system, the model define themselves their energy limits per material and particle.
                                                   >>  61     // Therefore, such a limit should not be in the G4DNAProcess classes.
                                                   >>  62     //
                                                   >>  63     SetLowEnergyLimit(0.);
                                                   >>  64     SetHighEnergyLimit(DBL_MAX);
 54                                                    65 
 55   fpParticleChangeForGamma = GetParticleChange <<  66     fpParticleChangeForGamma = GetParticleChangeForGamma();
 56                                                    67 
 57   // Loop on all the registered models to init <<  68     // Loop on all the registered models to initialise them
 58   for (auto & fRegisteredModel : fRegisteredMo <<  69     for(unsigned int i=0, ie = fRegisteredModels.size(); i<ie; ++i)
 59     fRegisteredModel->SetParticleChange(fpPart <<  70     {
 60     fRegisteredModel->Initialise(particle, cut <<  71         fRegisteredModels.at(i)->Initialise(particle, cuts, fpParticleChangeForGamma);
 61   }                                            <<  72     }
 62   // used to retrieve the model corresponding  << 
 63   BuildMaterialParticleModelTable(particle);   << 
 64                                                    73 
 65   BuildMaterialMolPerVolTable();               << 
 66                                                    74 
 67   StreamInfo(G4cout);                          <<  75     // Build the [material][particle]=Models table
                                                   >>  76     // used to retrieve the model corresponding to the current material/particle couple
                                                   >>  77     BuildMaterialParticleModelTable(particle);
 68                                                    78 
                                                   >>  79     BuildMaterialMolPerVolTable();
 69 }                                                  80 }
 70                                                    81 
 71 //....oooOO0OOooo........oooOO0OOooo........oo     82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 72                                                    83 
 73 G4double G4DNAModelInterface::CrossSectionPerV     84 G4double G4DNAModelInterface::CrossSectionPerVolume(const G4Material* material,
 74                                                <<  85                                                  const G4ParticleDefinition* p,
 75                                                <<  86                                                  G4double ekin,
 76 {                                              <<  87                                                  G4double emin,
 77   // Method to return the crossSection * nbMol <<  88                                                  G4double emax)
 78   // Process class then calculates the path.   <<  89 {
 79   // The cross section is calculated in the re <<  90     // Method to return the crossSection * nbMoleculePerUnitVolume to the process class.
 80   // Two cases are handled here: normal materi <<  91     // Process class then calculates the path.
 81   //                                           <<  92     // The cross section is calculated in the registered model(s) and this class just call the method
 82   // Idea:                                     <<  93     // Two cases are handled here: normal material and composite material.
 83   // *** Simple material ***                   <<  94     //
 84   // Ask for the cross section of the chosen m <<  95     // Idea:
 85   // Multiply it by the number of medium molec <<  96     // *** Simple material ***
 86   // Return the value.                         <<  97     // Ask for the cross section of the chosen model.
 87   // *** Composite material ***                <<  98     // Multiply it by the number of medium molecules per volume unit.
 88   // Ask for the cross section of the chosen m <<  99     // Return the value.
 89   // Apply a factor to each cross section and  << 100     // *** Composite material ***
 90   // component per composite volume unit. The  << 101     // Ask for the cross section of the chosen model for each component.
 91                                                << 102     // Apply a factor to each cross section and sum the results. The factor is the molecule number of component per composite volume unit.
 92   // To reset the sampledMat variable.         << 103     // The total cross section is returned.
 93   // Can be used by user to retrieve current c << 104 
 94   fSampledMat = 0;                             << 105     // To reset the sampledMat variable.
 95                                                << 106     // Can be used by user to retrieve current component
 96   // This is the value to be sum up and to be  << 107     fSampledMat = "";
 97   G4double crossSectionTimesNbMolPerVol(0.);   << 108 
 98                                                << 109     // This is the value to be sum up and to be returned at then end
 99   // Reset the map saving the material and the << 110     G4double crossSectionTimesNbMolPerVol (0);
100   // Used in SampleSecondaries if the interact << 111 
101   // composite                                 << 112     // Reset the map saving the material and the cumulated corresponding cross section
102   fMaterialCS.clear();                         << 113     // Used in SampleSecondaries if the interaction is selected for the step and if the material is a composite
103                                                << 114     fMaterialCS.clear();
104   // This is the value to be used by SampleSec << 115 
105   fCSsumTot = 0;                               << 116      // This is the value to be used by SampleSecondaries
106                                                << 117     fCSsumTot = 0;
107   // *****************************             << 118 
108   // Material is not a composite               << 119     // *****************************
109   // *****************************             << 120     // Material is not a composite
110   //                                           << 121     // *****************************
111   if (material->GetMatComponents().empty()) {  << 122     //
112     // Get the material name                   << 123     if(material->GetMatComponents().empty())
113     const size_t & materialID = material->GetI << 124     {
114                                                << 125         // Get the material name
115     // Use the table to get the  model         << 126         const G4String& materialName = material->GetName();
116     auto model = SelectModel(materialID, p, ek << 127 
117                                                << 128         // Use the table to get the  model
118     // Get the nunber of molecules per volume  << 129         G4VDNAModel* model = GetDNAModel(materialName, p->GetParticleName(), ekin);
119                                                << 130 
120     // Calculate the cross section times the n << 131         // Get the nunber of molecules per volume unit for that material
121     if (model != nullptr) {                    << 132         G4double nbOfMoleculePerVolumeUnit = GetNumMoleculePerVolumeUnitForMaterial(material);
122       if (dynamic_cast<G4VDNAModel*>(model) == << 133 
123         // water material models only          << 134         // Calculate the cross section times the number of molecules
124         crossSectionTimesNbMolPerVol = model-> << 135         if(model != 0)
125       }                                        << 136             crossSectionTimesNbMolPerVol = nbOfMoleculePerVolumeUnit * model->CrossSectionPerVolume(material, materialName, p, ekin, emin, emax);
126       else {                                   << 137         else // no model was selected, we are out of the energy ranges
127         crossSectionTimesNbMolPerVol = model-> << 138             crossSectionTimesNbMolPerVol = 0.;
128       }                                        << 139     }
129     }                                          << 140 
130     else  // no model was selected, we are out << 141     // ********************************
131       crossSectionTimesNbMolPerVol = 0.;       << 142     // Material is a composite
132   }                                            << 143     // ********************************
133                                                << 144     //
134   // ********************************          << 145     else
135   // Material is a composite                   << 146     {
136   // ********************************          << 147         // Copy the map in a local variable
137   //                                           << 148         // Otherwise we get segmentation fault and iterator pointing to nowhere: do not know why...
138   else {                                       << 149         // Maybe MatComponents map is overrided by something somewhere ?
139     // Copy the map in a local variable        << 150         std::map<G4Material*, G4double> componentsMap = material->GetMatComponents();
140     // Otherwise we get segmentation fault and << 151 
141     // Maybe MatComponents map is overrided by << 152         // Retrieve the iterator
142     auto componentsMap = material->GetMatCompo << 153         std::map<G4Material*, G4double>::const_iterator it = componentsMap.begin();
143                                                << 154 
144     G4cout << material->GetName() << G4endl;   << 155         // Get the size
145                                                << 156         unsigned int componentNumber = componentsMap.size();
146     // Loop on all the components              << 157 
147     for (const auto& it : componentsMap) {     << 158         // Loop on all the components
148       // Get the current component             << 159         //for(it = material->GetMatComponents().begin(); it!=material->GetMatComponents().end();++it)
149       auto component = it.first;               << 160         for(unsigned int i=0; i<componentNumber; ++i)
150       // Get the current component mass fracti << 161         {
151       // G4double massFraction = it->second;   << 162             // Get the current component
152                                                << 163             G4Material* component = it->first;
153       // Get the number of component molecules << 164 
154       G4double nbMoleculeOfComponentInComposit << 165             // Get the current component mass fraction
155         GetNumMolPerVolUnitForComponentInCompo << 166             //G4double massFraction = it->second;
156       G4cout << " ==========>component : " <<  << 167 
157              << " nbMoleculeOfComponentInCompo << 168             // Get the number of component molecules in a volume unit of composite material
158              << G4endl;                        << 169             G4double nbMoleculeOfComponentInCompositeMat = GetNumMolPerVolUnitForComponentInComposite(component, material);
159                                                << 170 
160       // Get the current component name        << 171             // Get the current component name
161       const std::size_t & componentID = compon << 172             const G4String componentName = component->GetName();
162                                                << 173 
163       // Retrieve the model corresponding to t << 174             // Retrieve the model corresponding to the current component (ie material)
164       auto model = SelectModel(componentID, p, << 175             G4VDNAModel* model = GetDNAModel(componentName, p->GetParticleName(), ekin);
165                                                << 176 
166       // Add the component part of the cross s << 177             // Add the component part of the cross section to the cross section variable.
167       // The component cross section is multip << 178             // The component cross section is multiplied by the total molecule number in the composite scaled by the mass fraction.
168       // scaled by the mass fraction.          << 179             if(model != 0)
169       G4double crossSection;                   << 180                 crossSectionTimesNbMolPerVol =
170       if (model != nullptr) {                  << 181                         nbMoleculeOfComponentInCompositeMat * model->CrossSectionPerVolume(component, componentName, p, ekin, emin, emax);
171         if (dynamic_cast<G4VDNAModel*>(model)  << 182             else // no model was selected, we are out of the energy ranges
172           // water models                      << 183                 crossSectionTimesNbMolPerVol = 0.;
173           crossSection =                       << 184 
174             model->CrossSectionPerVolume(compo << 185             // Save the component name and its calculated crossSectionTimesNbMolPerVol
175             / GetNumMoleculePerVolumeUnitForMa << 186             // To be used by sampling secondaries if the interaction is selected for the step
176         }                                      << 187             fMaterialCS[componentName] = crossSectionTimesNbMolPerVol;
177         else {                                 << 188 
178           crossSection = model->CrossSectionPe << 189             // Save the component name and its calculated crossSectionTimesNbMolPerVol
179                          / GetNumMoleculePerVo << 190             // To be used by sampling secondaries if the interaction is selected for the step
                                                   >> 191             fCSsumTot += crossSectionTimesNbMolPerVol;
                                                   >> 192 
                                                   >> 193             // Move forward the iterator
                                                   >> 194             ++it;
180         }                                         195         }
181         crossSectionTimesNbMolPerVol = nbMolec << 196 
182       }                                        << 197         crossSectionTimesNbMolPerVol = fCSsumTot;
183       else  // no model was selected, we are o << 198 
184       {                                        << 199     }
185         crossSectionTimesNbMolPerVol = 0.;     << 200 
186       }                                        << 201     // return the cross section times the number of molecules
187                                                << 202     // the path of the interaction will be calculated using that value
188       // Save the component name and its calcu << 203     return crossSectionTimesNbMolPerVol;
189       // To be used by sampling secondaries if << 
190       fMaterialCS[componentID] = crossSectionT << 
191                                                << 
192       // Save the component name and its calcu << 
193       // To be used by sampling secondaries if << 
194       fCSsumTot += crossSectionTimesNbMolPerVo << 
195     }                                          << 
196     crossSectionTimesNbMolPerVol = fCSsumTot;  << 
197   }                                            << 
198                                                << 
199   // return the cross section times the number << 
200   // the path of the interaction will be calcu << 
201   return crossSectionTimesNbMolPerVol;         << 
202 }                                                 204 }
203                                                   205 
204 //....oooOO0OOooo........oooOO0OOooo........oo    206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
205                                                   207 
206 void G4DNAModelInterface::SampleSecondaries(st    208 void G4DNAModelInterface::SampleSecondaries(std::vector<G4DynamicParticle*>* fVect,
207                                             co << 209                                          const G4MaterialCutsCouple* couple,
208                                             co << 210                                          const G4DynamicParticle* aDynamicParticle,
209                                             G4 << 211                                          G4double tmin,
210 {                                              << 212                                          G4double tmax)
211   // To call the sampleSecondaries method of t << 213 {
212   // In the case of composite material, we nee << 214     // To call the sampleSecondaries method of the registered model(s)
213   // To do so we use a random sampling on the  << 215     // In the case of composite material, we need to choose a component to call the method from.
214   // CrossSectionPerVolume method. If we enter << 216     // To do so we use a random sampling on the crossSectionTimesNbMolPerVol used in CrossSectionPerVolume method.
215   // (and process) has been chosen for the cur << 217     // If we enter that method it means the corresponding interaction (and process) has been chosen for the current step.
216                                                << 218 
217   std::size_t materialID;                      << 219     G4String materialName;
218                                                << 220 
219   // *******************************           << 221     // *******************************
220   // Material is not a composite               << 222     // Material is not a composite
221   // *******************************           << 223     // *******************************
222   //                                           << 224     //
223   if (couple->GetMaterial()->GetMatComponents( << 225     if(couple->GetMaterial()->GetMatComponents().empty())
224     materialID = couple->GetMaterial()->GetInd << 226     {
225   }                                            << 227         materialName = couple->GetMaterial()->GetName();
226                                                << 228     }
227   // ****************************              << 229 
228   // Material is a composite                   << 230     // ****************************
229   // ****************************              << 
230   //                                           << 
231   else {                                       << 
232     // Material is a composite                    231     // Material is a composite
233     // We need to select a component           << 232     // ****************************
                                                   >> 233     //
                                                   >> 234     else
                                                   >> 235     {
                                                   >> 236         // Material is a composite
                                                   >> 237         // We need to select a component
                                                   >> 238 
                                                   >> 239         // We select a random number between 0 and fCSSumTot
                                                   >> 240         G4double rand = G4UniformRand()*fCSsumTot;
                                                   >> 241 
                                                   >> 242         G4double cumulCS (0);
                                                   >> 243 
                                                   >> 244         G4bool result = false;
                                                   >> 245 
                                                   >> 246         // We loop on each component cumulated cross section
                                                   >> 247         //
                                                   >> 248         // Retrieve the iterators
                                                   >> 249         std::map<const G4String , G4double>::const_iterator it = fMaterialCS.begin();
                                                   >> 250         std::map<const G4String , G4double>::const_iterator ite = fMaterialCS.end();
                                                   >> 251         // While this is true we do not have found our component.
                                                   >> 252         while(rand>cumulCS)
                                                   >> 253         {
                                                   >> 254             // Check if the sampling is ok
                                                   >> 255             if(it==ite)
                                                   >> 256             {
                                                   >> 257                 G4Exception("G4DNAModelManager::SampleSecondaries","em0006",
                                                   >> 258                             FatalException,
                                                   >> 259                             "The random component selection has failed: we ran into the end of the map without having a selected component");
                                                   >> 260                 return; // to make some compilers happy
                                                   >> 261             }
                                                   >> 262 
                                                   >> 263             // Set the cumulated value for the iteration
                                                   >> 264             cumulCS += it->second;
                                                   >> 265 
                                                   >> 266             // Check if we have reach the material to be selected
                                                   >> 267             // The DBL_MAX is here to take into account a return DBL_MAX in CSPerVol for the elastic model
                                                   >> 268             // to force elastic sampleSecondaries where the particle can be killed.
                                                   >> 269             // Used when paticle energy is lower than limit.
                                                   >> 270             if(rand<cumulCS || cumulCS >= DBL_MAX)
                                                   >> 271             {
                                                   >> 272                 // we have our selected material
                                                   >> 273                 materialName = it->first;
                                                   >> 274                 result = true;
                                                   >> 275                 break;
                                                   >> 276             }
                                                   >> 277 
                                                   >> 278             // make the iterator move forward
                                                   >> 279             ++it;
                                                   >> 280         }
234                                                   281 
235     // We select a random number between 0 and << 282         // Check that we get a result
236     G4double rand = G4UniformRand() * fCSsumTo << 283         if(!result)
                                                   >> 284         {
                                                   >> 285             // it is possible to end up here if the return DBL_MAX of CSPerVol in the elastic model is not taken into account
                                                   >> 286 
                                                   >> 287             G4Exception("G4DNAModelManager::SampleSecondaries","em0006",
                                                   >> 288                         FatalException,
                                                   >> 289                         "The random component selection has failed: while loop ended without a selected component.");
                                                   >> 290             return; // to make some compilers happy
                                                   >> 291         }
237                                                   292 
238     G4double cumulCS(0);                       << 293     }
239                                                   294 
240     G4bool result = false;                     << 295     // **************************************
                                                   >> 296     // Call the SampleSecondaries method
                                                   >> 297     // **************************************
241                                                   298 
242     // We loop on each component cumulated cro << 299     // Rename material if modified NIST material
243     //                                         << 300     // This is needed when material is obtained from G4MaterialCutsCouple
244     // Retrieve the iterators                  << 301     if(materialName.find("_MODIFIED")!=G4String::npos)
245     auto it = fMaterialCS.begin();             << 302     {
246     auto ite = fMaterialCS.end();              << 303         materialName = materialName.substr(0,materialName.size()-9);
247     // While this is true we do not have found << 304     }
248     while (rand > cumulCS) {                   << 
249       // Check if the sampling is ok           << 
250       if (it == ite) {                         << 
251         G4Exception(                           << 
252           "G4DNAModelManager::SampleSecondarie << 
253           "The random component selection has  << 
254           "having a selected component");      << 
255         return;  // to make some compilers hap << 
256       }                                        << 
257                                                << 
258       // Set the cumulated value for the itera << 
259       cumulCS += it->second;                   << 
260                                                << 
261       // Check if we have reach the material t << 
262       // The DBL_MAX is here to take into acco << 
263       // to force elastic sampleSecondaries wh << 
264       // Used when paticle energy is lower tha << 
265       if (rand < cumulCS || cumulCS >= DBL_MAX << 
266         // we have our selected material       << 
267         materialID = it->first;                << 
268         result = true;                         << 
269         break;                                 << 
270       }                                        << 
271                                                << 
272       // make the iterator move forward        << 
273       ++it;                                    << 
274     }                                          << 
275                                                << 
276     // Check that we get a result              << 
277     if (!result) {                             << 
278       // it is possible to end up here if the  << 
279       // taken into account                    << 
280                                                << 
281       G4Exception("G4DNAModelManager::SampleSe << 
282                   "The random component select << 
283                   "component.");               << 
284       return;  // to make some compilers happy << 
285     }                                          << 
286   }                                            << 
287                                                << 
288   // **************************************    << 
289   // Call the SampleSecondaries method         << 
290   // **************************************    << 
291                                                << 
292   // Rename material if modified NIST material << 
293   // This is needed when material is obtained  << 
294   //  if (materialName.find("_MODIFIED") != G4 << 
295   //    materialName = materialName.substr(0,  << 
296   //  }                                        << 
297                                                << 
298   fSampledMat = materialID;                    << 
299                                                << 
300   auto model = SelectModel(materialID, aDynami << 
301                            aDynamicParticle->G << 
302                                                   305 
303   model->SampleSecondaries(fVect, couple, aDyn << 306     fSampledMat = materialName;
                                                   >> 307 
                                                   >> 308     G4VDNAModel* model = GetDNAModel(materialName,
                                                   >> 309                                      aDynamicParticle->GetParticleDefinition()->GetParticleName(),
                                                   >> 310                                      aDynamicParticle->GetKineticEnergy() );
                                                   >> 311             //fMaterialParticleModelTable[materialName][aDynamicParticle->GetDefinition()->GetParticleName()][0];
                                                   >> 312 
                                                   >> 313     model->SampleSecondaries(fVect, couple, materialName, aDynamicParticle, fpParticleChangeForGamma, tmin, tmax);
304 }                                                 314 }
305                                                   315 
306 void G4DNAModelInterface::RegisterModel(G4VEmM << 316 void G4DNAModelInterface::RegisterModel(G4VDNAModel* model)
307 {                                                 317 {
308   fRegisteredModels.push_back(model);          << 318     fRegisteredModels.push_back(model);
                                                   >> 319 }
                                                   >> 320 
                                                   >> 321 void G4DNAModelInterface::RegisterModel(G4VEmModel* model, const G4ParticleDefinition* particle)
                                                   >> 322 {
                                                   >> 323     G4DNADummyModel* dummyWrapper = new G4DNADummyModel("G4_WATER", particle, model->GetName(), model);
                                                   >> 324 
                                                   >> 325     RegisterModel(dummyWrapper);
309 }                                                 326 }
310                                                   327 
311 //....oooOO0OOooo........oooOO0OOooo........oo    328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
312                                                   329 
313 void G4DNAModelInterface::BuildMaterialParticl    330 void G4DNAModelInterface::BuildMaterialParticleModelTable(const G4ParticleDefinition* p)
314 {                                                 331 {
315   // Method to build a map: [material][particl << 332     // Method to build a map: [material][particle] = Model*.
316   // The map is used to retrieve the correct m << 333     // The map is used to retrieve the correct model for the current particle/material couple.
                                                   >> 334 
                                                   >> 335     // Get the current particle name
                                                   >> 336     const G4String& pName = p->GetParticleName();
                                                   >> 337 
                                                   >> 338     // Retrieve the iterator
                                                   >> 339     G4MaterialTable::iterator it;
317                                                   340 
318   // Loop on all materials registered in the s << 341     // Loop on all materials registered in the simulation
319   for (auto it : *G4Material::GetMaterialTable << 342     for(it = G4Material::GetMaterialTable()->begin(); it!=G4Material::GetMaterialTable()->end(); ++it)
320     // Get the material pointer                << 343     {
321     G4Material* mat = it;                      << 344         // Get the material pointer
322     // Get the map                             << 345         G4Material* mat = *it;
323     // Check that the material is not a compos << 346 
324     auto componentMap = mat->GetMatComponents( << 347         // Get the map
325     if (componentMap.empty()) {                << 348         std::map<G4Material*, G4double> componentMap = mat->GetMatComponents();
326       // Get the material name                 << 349 
327       const std::size_t & matID = mat->GetInde << 350         // Get the number of component within the composite
328       InsertModelInTable(matID, p);            << 351         unsigned int compositeSize = componentMap.size();
329     }                                          << 352 
330     // if the material is a composite material << 353         // Check that the material is not a composite material
331     // register them                           << 354         if(componentMap.empty())
332     else {                                     << 355         {
333       // Loop on all the components of the mat << 356             // Get the material name
334       for (const auto& itComp : componentMap)  << 357             const G4String& matName = mat->GetName();
335         G4Material* component = itComp.first;  << 358 
336         // Check that the component is not its << 359             // Insert the model in the table.
337         if (!component->GetMatComponents().emp << 360             InsertModelInTable(matName, pName);
338           std::ostringstream oss;              << 361         }
339           oss << "Material " << mat->GetName() << 362         // if the material is a composite material then we need to loop on all its components to register them
340           oss << " " << component->GetName();  << 363         else
341           G4Exception("G4DNAModelManager::Buil << 364         {
342                       FatalException, oss.str( << 365             // Retrieve the component map begin iterator
343           return;  // to make some compilers h << 366             std::map<G4Material*, G4double>::const_iterator itComp = componentMap.begin();
                                                   >> 367 
                                                   >> 368             // Loop on all the components of the material
                                                   >> 369             //for(itComp = mat->GetMatComponents().begin(); itComp != eitComp; ++itComp)
                                                   >> 370             for(unsigned int k=0; k<compositeSize; ++k)
                                                   >> 371             {
                                                   >> 372                 G4Material* component = itComp->first;
                                                   >> 373 
                                                   >> 374 //                // Check that the component is not itself a composite
                                                   >> 375 //                if(component->GetMatComponents().size()!=0)
                                                   >> 376 //                {
                                                   >> 377 //                    std::ostringstream oss;
                                                   >> 378 //                    oss<<"Material "<<mat->GetName()<<" is a composite and its component ";
                                                   >> 379 //                    oss<<component->GetName()<<" is also a composite material. Building composite with other composites is not implemented yet";
                                                   >> 380 //                    oss<<G4endl;
                                                   >> 381 //                    G4Exception("G4DNAModelManager::BuildMaterialParticleModelTable","em0006",
                                                   >> 382 //                                FatalException, oss.str().c_str());
                                                   >> 383 //                    return; // to make some compilers happy
                                                   >> 384 //                }
                                                   >> 385 
                                                   >> 386                 // Get the current component name
                                                   >> 387                 const G4String compName = component->GetName();
                                                   >> 388 
                                                   >> 389                 // If there is a model then insert the model corresponding to the component in the table
                                                   >> 390                 // contains a if statement to check we have not registered the material as a component or a normal material before.
                                                   >> 391                 InsertModelInTable(compName, pName);
                                                   >> 392 
                                                   >> 393                 // move forward the iterator
                                                   >> 394                 ++itComp;
                                                   >> 395             }
344         }                                         396         }
345         // Get the current component name      << 
346         const std::size_t & compID = component << 
347         // If there is a model then insert the << 
348         // contains a if statement to check we << 
349         // normal material before.             << 
350         InsertModelInTable(compID, p);         << 
351         // move forward the iterator           << 
352       }                                        << 
353     }                                             397     }
354   }                                            << 
355 }                                                 398 }
356                                                   399 
357 //....oooOO0OOooo........oooOO0OOooo........oo    400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358                                                   401 
359 void G4DNAModelInterface::BuildMaterialMolPerV    402 void G4DNAModelInterface::BuildMaterialMolPerVolTable()
360 {                                                 403 {
361   // To be sure the G4DNAMolecularMaterial is  << 404     // To be sure the G4DNAMolecularMaterial is initialized
362   G4DNAMolecularMaterial::Instance()->Initiali << 405     G4DNAMolecularMaterial::Instance()->Initialize();
363                                                   406 
364   G4MaterialTable* materialTable = G4Material: << 407     G4MaterialTable* materialTable = G4Material::GetMaterialTable();
365                                                   408 
366   // Loop on all the materials inside the "mat << 409     // Loop on all the materials inside the "materialTable"
367   for (auto currentMaterial : *materialTable)  << 410     for(size_t i=0, ie=materialTable->size(); i<ie; i++)
368     // Current material                        << 411     {
369     // Current material name                   << 412         // Current material
370     const std::size_t & currentMatID = current << 413         G4Material* currentMaterial = materialTable->at(i);
371                                                << 414 
372     // Will the material be used in this inter << 415         // Current material name
373     // Loop on all the materials that can be d << 416         const G4String& currentMatName = currentMaterial->GetName();
374     auto it = fMaterialParticleModelTable.begi << 417 
375     auto ite = fMaterialParticleModelTable.end << 418         // Will the material be used in this interface instance ?
376     for (; it != ite; it++) {                  << 419         // Loop on all the materials that can be dealt with in this class
377       const std::size_t & materialID = it->fir << 420         MaterialParticleModelTable::iterator it = fMaterialParticleModelTable.begin();
378                                                << 421         MaterialParticleModelTable::iterator ite = fMaterialParticleModelTable.end();
379       if (materialID == currentMatID) {        << 422         for(; it != ite; it++)
380         const std::vector<G4double>* numMolPer << 423         {
381           G4DNAMolecularMaterial::Instance()-> << 424             const G4String& materialName = it->first;
382         fMaterialMolPerVol[materialID] = numMo << 425 
383       }                                        << 426             if(materialName == currentMatName)
                                                   >> 427             {
                                                   >> 428                 const std::vector<double>* numMolPerVolForMat = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(currentMaterial);
                                                   >> 429                 fMaterialMolPerVol[materialName] = numMolPerVolForMat;
                                                   >> 430             }
                                                   >> 431         }
384     }                                             432     }
385   }                                            << 
386 }                                                 433 }
387                                                   434 
388 //....oooOO0OOooo........oooOO0OOooo........oo    435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389                                                   436 
390 void G4DNAModelInterface::InsertModelInTable(c << 437 void G4DNAModelInterface::InsertModelInTable(const G4String& matName, const G4String& pName)
391 {                                                 438 {
392   // To insert the model(s) in the table Mater << 439     // To insert the model(s) in the table Material Particule -> Model(s)
393                                                   440 
394   // First, we need to check if the current ma << 441     // First, we need to check if the current material has already been inserted in the table.
395   // This is possible because of the composite << 442     // This is possible because of the composite material. We could add a component M1 and then try to add the independant M1 material.
396   // add the independant M1 material. This cas << 443     // This case must be avoided. Checking if M1 is already in the table is the way to avoid it.
397   // table is the way to avoid it.             << 444     //
398   //                                           << 445     // Chech if the current material and particle are already in the table.
399   // Check if the current material and particl << 446     // If they are: do nothing.
400   // If they are: do nothing.                  << 447     // If they are not: add the model(s)
401   // If they are not: add the model(s)         << 448     //
402   //                                           << 449     // Check for the material
403   // Check for the material                    << 450     if(fMaterialParticleModelTable.find(matName) == fMaterialParticleModelTable.end())
404   if (fMaterialParticleModelTable.find(matID)  << 451     {
405     // Check for the particle                  << 452         // Check for the particle
406     if (fMaterialParticleModelTable[matID].fin << 453         if(fMaterialParticleModelTable[matName].find(pName) == fMaterialParticleModelTable[matName].end())
407       G4int modelNbForMaterial = 0;            << 454         {
408       for (const auto& it : fRegisteredModels) << 455             G4int modelNbForMaterial (0);
409         auto model = dynamic_cast<G4VDNAModel* << 456 
410         if (model != nullptr) {                << 457             // Loop on all models registered in the simulation to check:
411           if (model->IsParticleExistingInModel << 458             // 1- if they can be applied to the current material
412             fMaterialParticleModelTable[matID] << 459             // 2- if they can be applied to the current particle
413             // and add one to the "there is a  << 460             for(unsigned int i=0, ie=fRegisteredModels.size(); i<ie; ++i)
414             ++modelNbForMaterial;              << 461             {
415           }                                    << 462                 // check if the model is correct for material and particle (previous 1 and 2)
                                                   >> 463                 if(fRegisteredModels[i]->IsParticleExistingInModelForMaterial(pName, matName))
                                                   >> 464                 {
                                                   >> 465                     // if yes then add the model in the map
                                                   >> 466                     fMaterialParticleModelTable[matName][pName].push_back(fRegisteredModels[i]);
                                                   >> 467 
                                                   >> 468                     // and add one to the "there is a model" material flag
                                                   >> 469                     ++modelNbForMaterial;
                                                   >> 470                 }
                                                   >> 471             }
                                                   >> 472 
                                                   >> 473             // The model(s) applicable to the currently selected material should be in the map.
                                                   >> 474             // We check if there are several models for the material.
                                                   >> 475             if(modelNbForMaterial>1)
                                                   >> 476             {
                                                   >> 477                 // If there are several models for a given material and particle couple it could be
                                                   >> 478                 // because of the energy ranges. We will check if the energy ranges are coherent.
                                                   >> 479 
                                                   >> 480                 // Get the models (vector)
                                                   >> 481                 std::vector<G4VDNAModel*>& models = fMaterialParticleModelTable[matName][pName];
                                                   >> 482 
                                                   >> 483                 // Declare a map to sort the limits (G4double) and a model "id" (G4int).
                                                   >> 484                 // The model id is created on the fly here.
                                                   >> 485                 // The idea is to fill a map with [limit] = counter. This map will be auto-sorted
                                                   >> 486                 // and we will check by iterating on it that the counter order is maintained.
                                                   >> 487 
                                                   >> 488                 // Delcare the map
                                                   >> 489                 std::map<G4double, G4int, std::less<G4double> > sortMap;
                                                   >> 490 
                                                   >> 491                 G4double smallDiff = 0.01 *eV;
                                                   >> 492 
                                                   >> 493                 // Loop on all the model for the current couple
                                                   >> 494                 // and fill a map with [lim] = modelNumber
                                                   >> 495                 for(unsigned int ii=0, em=models.size(); ii<em; ++ii)
                                                   >> 496                 {
                                                   >> 497                     G4double lowLim = models[ii]->GetLowELimit(matName, pName);
                                                   >> 498                     G4double highLim = models[ii]->GetHighELimit(matName, pName);
                                                   >> 499 
                                                   >> 500                     if(sortMap.find(lowLim) != sortMap.end() )
                                                   >> 501                     {
                                                   >> 502                         lowLim += smallDiff;
                                                   >> 503                     }
                                                   >> 504 
                                                   >> 505                     sortMap[lowLim] = ii;
                                                   >> 506 
                                                   >> 507                     if(sortMap.find(highLim) != sortMap.end() )
                                                   >> 508                     {
                                                   >> 509                         highLim -= smallDiff;
                                                   >> 510                     }
                                                   >> 511 
                                                   >> 512                     sortMap[highLim] = ii;
                                                   >> 513                 }
                                                   >> 514 
                                                   >> 515                 // The map has been created and ordered at this point.
                                                   >> 516                 // We will check the map order.
                                                   >> 517 
                                                   >> 518                 // Loop on the sortMap with iterator and check the order is correct.
                                                   >> 519                 std::map<G4double, G4int>::iterator it = sortMap.begin();
                                                   >> 520 
                                                   >> 521                 // First energy limit value
                                                   >> 522                 G4double dummyLim = it->first - smallDiff;
                                                   >> 523 
                                                   >> 524                 // Loop on all the models again.
                                                   >> 525                 // The goal is to check if for each limit pairs we have the same model number
                                                   >> 526                 // and that the upper and lower limit are consistent.
                                                   >> 527                 for(unsigned int ii=0, eii=models.size(); ii<eii; ++ii)
                                                   >> 528                 {
                                                   >> 529                     G4double lim1 = it->first - smallDiff;
                                                   >> 530                     G4int count1 = it->second;
                                                   >> 531 
                                                   >> 532                     // Iterate
                                                   >> 533                     ++it;
                                                   >> 534 
                                                   >> 535                     G4double lim2 = it->first + smallDiff;
                                                   >> 536                     G4int count2 = it->second;
                                                   >> 537 
                                                   >> 538                     // Iterate
                                                   >> 539                     ++it;
                                                   >> 540 
                                                   >> 541                     // Check model number and energy limit consistency
                                                   >> 542                     // std::abs(dummyLim - lim1) > 1.*eV because we cannot do (dummyLim != lim1)
                                                   >> 543                     // without experimenting precision loss. Therefore, the std::abs(...) > tolerance is the usual way of avoiding
                                                   >> 544                     // the issue.
                                                   >> 545                     if( (count1 != count2) || ( std::abs(dummyLim - lim1) > 1.*eV ) )
                                                   >> 546                     {
                                                   >> 547                         // Error
                                                   >> 548 
                                                   >> 549                         std::ostringstream oss;
                                                   >> 550                         oss<<"The material "<<matName<<" and the particle "<<pName;
                                                   >> 551                         oss<<" have several models registered for the "<<fName<<" interaction and their energy ranges ";
                                                   >> 552                         oss<<"do not match. \nEnergy ranges: \n";
                                                   >> 553 
                                                   >> 554                         for(int iii=0, eiii=models.size(); iii<eiii; ++iii)
                                                   >> 555                         {
                                                   >> 556                             oss<<models[iii]->GetName()<<"\n";
                                                   >> 557                             oss<<"low: "<<models[iii]->GetLowELimit(matName, pName)/eV<<" eV \n";
                                                   >> 558                             oss<<"high: "<<models[iii]->GetHighELimit(matName, pName)/eV<<" eV \n";
                                                   >> 559                         }
                                                   >> 560 
                                                   >> 561                         G4Exception("G4DNAModelManager::InsertModelInTable","em0006",
                                                   >> 562                                     FatalException, oss.str().c_str());
                                                   >> 563                         return; // to make some compilers happy
                                                   >> 564                     }
                                                   >> 565 
                                                   >> 566                     dummyLim = lim2;
                                                   >> 567                 }
                                                   >> 568 
                                                   >> 569                 // If we are here then everything was ok.
                                                   >> 570             }
                                                   >> 571             // no model for the material case
                                                   >> 572             else if(modelNbForMaterial==0)
                                                   >> 573             {
                                                   >> 574 //                std::ostringstream oss;
                                                   >> 575 //                oss<<"The material "<<matName<<" and the particle "<<pName;
                                                   >> 576 //                oss<<" does not have any model registered for the "<<fName<<" interaction. ";
                                                   >> 577 
                                                   >> 578 //                G4Exception("G4DNAModelManager::InsertModelInTable","em0006",
                                                   >> 579 //                            FatalException, oss.str().c_str());
                                                   >> 580 //                return; // to make some compilers happy
                                                   >> 581             }
416         }                                         582         }
417         else {                                 << 583     }
418           auto index = fpG4_WATER->GetIndex(); << 584 }
419           fMaterialParticleModelTable[index][p << 585 
420           ++modelNbForMaterial;                << 586 G4VDNAModel *G4DNAModelInterface::GetDNAModel(const G4String &material, const G4String &particle, G4double ekin)
                                                   >> 587 {
                                                   >> 588     // Output pointer
                                                   >> 589     G4VDNAModel* model = 0;
                                                   >> 590 
                                                   >> 591     // Get a reference to all the models for the couple (material and particle)
                                                   >> 592     std::vector<G4VDNAModel*>& models = fMaterialParticleModelTable[material][particle];
                                                   >> 593 
                                                   >> 594     // We must choose one of the model(s) accordingly to the particle energy and the model energy range(s)
                                                   >> 595 
                                                   >> 596     //G4bool isOneModelSelected = false;
                                                   >> 597 
                                                   >> 598     // Loop on all the models within the models vector and check if ekin is within the energy range.
                                                   >> 599     for(int i=0, ie=models.size(); i<ie; ++i)
                                                   >> 600     {
                                                   >> 601         // ekin is in the energy range: we select the model and stop the loop.
                                                   >> 602         if( ekin >= models[i]->GetLowELimit(material, particle)
                                                   >> 603                 && ekin < models[i]->GetHighELimit(material, particle) )
                                                   >> 604         {
                                                   >> 605             // Select the model
                                                   >> 606             model = models[i];
                                                   >> 607 
                                                   >> 608             // Boolean flag
                                                   >> 609             //isOneModelSelected = true;
                                                   >> 610 
                                                   >> 611             // Quit the for loop
                                                   >> 612             break;
421         }                                         613         }
422       }                                        << 614 
423       if (modelNbForMaterial == 0) {           << 615         // ekin is not in the energy range: we continue the loop.
424         std::ostringstream oss;                << 616     }
425         oss << "The material " << (*G4Material << 617 
426             << " and the particle " << p->GetP << 618 //    // If no model was selected then fatal error
427         oss << " does not have any model regis << 619 //    if(!isOneModelSelected)
428         G4Exception("G4DNAModelInterface::Inse << 620 //    {
429                     oss.str().c_str());        << 621 //        G4String msg = "No model has ";
430         return;  // to make some compilers hap << 622 //        msg += ekin/eV;
431       }                                        << 623 //        msg += " eV in its energy range. Therefore nothing was selected.";
432     }                                          << 624 
433   }                                            << 625 //        G4Exception("G4DNAModelManager::GetDNAModel","em0006",
434 }                                              << 626 //                    FatalException,
435                                                << 627 //                    msg);
436 G4VEmModel* G4DNAModelInterface::SelectModel(c << 628 //    }
437                                              c << 629 
438                                              c << 630     // Return a pointer to the selected model
439 {                                              << 631     return model;
440   // Output pointer                            << 
441   G4VEmModel* model = nullptr;                 << 
442                                                << 
443   // Get a reference to all the models for the << 
444   auto modelData = fMaterialParticleModelTable << 
445                                                << 
446   // We must choose one of the model(s) accord << 
447   // range(s)                                  << 
448                                                << 
449   // Loop on all the models within the models  << 
450   auto DNAModel = dynamic_cast<G4VDNAModel*>(m << 
451   G4double lowL, highL;                        << 
452   if (DNAModel == nullptr) {                   << 
453     // ekin is in the energy range: we select  << 
454     lowL = modelData->LowEnergyLimit();        << 
455     highL = modelData->HighEnergyLimit();      << 
456     if (ekin >= lowL && ekin < highL) {        << 
457       // Select the model                      << 
458                                                << 
459       model = modelData;                       << 
460       // return model;                         << 
461       //  Quit the for loop                    << 
462       // break;                                << 
463     }                                          << 
464     // ekin is not in the energy range: we con << 
465   }                                            << 
466   else {                                       << 
467     // ekin is in the energy range: we select  << 
468     lowL = DNAModel->GetLowELimit(materialID,  << 
469     highL = DNAModel->GetHighELimit(materialID << 
470     if (ekin >= lowL && ekin < highL) {        << 
471       // Select the model                      << 
472       model = modelData;                       << 
473       // return model;                         << 
474       //  Quit the for loop                    << 
475       // break;                                << 
476     }                                          << 
477     // ekin is not in the energy range: we con << 
478   }                                            << 
479   //}                                          << 
480   return model;                                << 
481 }                                                 632 }
482                                                   633 
483 //....oooOO0OOooo........oooOO0OOooo........oo    634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
484                                                   635 
485 G4double G4DNAModelInterface::GetNumMoleculePe    636 G4double G4DNAModelInterface::GetNumMoleculePerVolumeUnitForMaterial(const G4Material* mat)
486 {                                                 637 {
487   return fMaterialMolPerVol[mat->GetIndex()]-> << 638     return fMaterialMolPerVol[mat->GetName()]->at(mat->GetIndex() );
488 }                                                 639 }
489                                                   640 
490 //....oooOO0OOooo........oooOO0OOooo........oo    641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
491                                                   642 
492 G4double                                       << 643 G4double G4DNAModelInterface::GetNumMolPerVolUnitForComponentInComposite(const G4Material* component, const G4Material* composite)
493 G4DNAModelInterface::GetNumMolPerVolUnitForCom << 644 {
494                                                << 645     return fMaterialMolPerVol[component->GetName() ]->at(composite->GetIndex() );
495 {                                              << 
496   return fMaterialMolPerVol[component->GetInde << 
497 }                                              << 
498                                                << 
499 void G4DNAModelInterface::StreamInfo(std::ostr << 
500 {                                              << 
501   G4long prec = os.precision(5);               << 
502   os << "===================================== << 
503      << this->GetName() << "      ============ << 
504      << "\n";                                  << 
505   os << std::setw(15) << "Material#" << std::s << 
506      << std::setw(17) << "LowLimit(MeV)" << st << 
507      << "Fast" << std::setw(13) << "Stationary << 
508   for (const auto& it1 : fMaterialParticleMode << 
509     os << std::setw(15) << (*G4Material::GetMa << 
510     for (const auto& it2 : it1.second) {       << 
511       os << std::setw(13) << it2.first->GetPar << 
512       os << std::setw(35) << it2.second->GetNa << 
513       auto DNAModel = dynamic_cast<G4VDNAModel << 
514       if (DNAModel == nullptr) {               << 
515         os << std::setw(17) << it2.second->Low << 
516         os << std::setw(17) << it2.second->Hig << 
517       }                                        << 
518       else {                                   << 
519         auto lowL = DNAModel->GetLowELimit(it1 << 
520         auto highL = DNAModel->GetHighELimit(i << 
521         os << std::setw(17) << lowL;           << 
522         os << std::setw(17) << highL;          << 
523       }                                        << 
524       os << std::setw(13) << "no";             << 
525       os << std::setw(13) << "no";             << 
526       os << std::setw(13) << "no" << G4endl;   << 
527     }                                          << 
528   }                                            << 
529                                                << 
530   os << "===================================== << 
531         "===================================== << 
532      << G4endl;                                << 
533   os.precision(prec);                          << 
534 }                                                 646 }
535                                                   647