Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNADingfelderChargeDecreaseModel.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 
 28 #include "G4DNADingfelderChargeDecreaseModel.hh"
 29 #include "G4PhysicalConstants.hh"
 30 #include "G4SystemOfUnits.hh"
 31 #include "G4DNAMolecularMaterial.hh"
 32 #include "G4DNAChemistryManager.hh"
 33 #include "G4Log.hh"
 34 #include "G4Pow.hh"
 35 #include "G4Alpha.hh"
 36 
 37 
 38 static G4Pow * gpow = G4Pow::GetInstance();
 39 
 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 41 
 42 using namespace std;
 43 
 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45 
 46 G4DNADingfelderChargeDecreaseModel::G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition*,
 47                                                                        const G4String& nam) :
 48 G4VEmModel(nam) 
 49 {
 50   numberOfPartialCrossSections[0] = 0;
 51   numberOfPartialCrossSections[1] = 0;
 52   numberOfPartialCrossSections[2] = 0;
 53 
 54   verboseLevel = 0;
 55   // Verbosity scale:
 56   // 0 = nothing
 57   // 1 = warning for energy non-conservation
 58   // 2 = details of energy budget
 59   // 3 = calculation of cross sections, file openings, sampling of atoms
 60   // 4 = entering in methods
 61 
 62   if (verboseLevel > 0)
 63   {
 64     G4cout << "Dingfelder charge decrease model is constructed " << G4endl;
 65   }
 66   // Selection of stationary mode
 67 
 68   statCode = false;
 69 }
 70 
 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 72 
 73 void G4DNADingfelderChargeDecreaseModel::Initialise(const G4ParticleDefinition* particle,
 74                                                     const G4DataVector& /*cuts*/)
 75 {
 76 
 77   if (verboseLevel > 3)
 78   {
 79     G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()"
 80         << G4endl;
 81   }
 82 
 83   // Energy limits
 84 
 85   G4DNAGenericIonsManager *instance;
 86   instance = G4DNAGenericIonsManager::Instance();
 87   protonDef = G4Proton::ProtonDefinition();
 88   alphaPlusPlusDef = G4Alpha::Alpha();
 89   alphaPlusDef = instance->GetIon("alpha+");
 90   hydrogenDef = instance->GetIon("hydrogen");
 91   heliumDef = instance->GetIon("helium");
 92 
 93   G4String proton;
 94   G4String alphaPlusPlus;
 95   G4String alphaPlus;
 96 
 97   // Limits
 98 
 99   proton = protonDef->GetParticleName();
100   lowEnergyLimit[proton] = 100. * eV;
101   highEnergyLimit[proton] = 100. * MeV;
102 
103   alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
104   lowEnergyLimit[alphaPlusPlus] = 1. * keV;
105   highEnergyLimit[alphaPlusPlus] = 400. * MeV;
106 
107   alphaPlus = alphaPlusDef->GetParticleName();
108   lowEnergyLimit[alphaPlus] = 1. * keV;
109   highEnergyLimit[alphaPlus] = 400. * MeV;
110 
111   //
112 
113   if (particle==protonDef)
114   {
115     SetLowEnergyLimit(lowEnergyLimit[proton]);
116     SetHighEnergyLimit(highEnergyLimit[proton]);
117   }
118 
119   if (particle==alphaPlusPlusDef)
120   {
121     SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
122     SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
123   }
124 
125   if (particle==alphaPlusDef)
126   {
127     SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
128     SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
129   }
130 
131   // Final state
132 
133   //PROTON
134   f0[0][0]=1.;
135   a0[0][0]=-0.180;
136   a1[0][0]=-3.600;
137   b0[0][0]=-18.22;
138   b1[0][0]=-1.997;
139   c0[0][0]=0.215;
140   d0[0][0]=3.550;
141   x0[0][0]=3.450;
142   x1[0][0]=5.251;
143 
144   numberOfPartialCrossSections[0] = 1;
145 
146   //ALPHA++
147   f0[0][1]=1.; a0[0][1]=0.95;
148   a1[0][1]=-2.75;
149   b0[0][1]=-23.00;
150   c0[0][1]=0.215;
151   d0[0][1]=2.95;
152   x0[0][1]=3.50;
153 
154   f0[1][1]=1.;
155   a0[1][1]=0.95;
156   a1[1][1]=-2.75;
157   b0[1][1]=-23.73;
158   c0[1][1]=0.250;
159   d0[1][1]=3.55;
160   x0[1][1]=3.72;
161 
162   x1[0][1]=-1.;
163   b1[0][1]=-1.;
164 
165   x1[1][1]=-1.;
166   b1[1][1]=-1.;
167 
168   numberOfPartialCrossSections[1] = 2;
169 
170   // ALPHA+
171   f0[0][2]=1.;
172   a0[0][2]=0.65;
173   a1[0][2]=-2.75;
174   b0[0][2]=-21.81;
175   c0[0][2]=0.232;
176   d0[0][2]=2.95;
177   x0[0][2]=3.53;
178 
179   x1[0][2]=-1.;
180   b1[0][2]=-1.;
181 
182   numberOfPartialCrossSections[2] = 1;
183 
184   //
185 
186   if( verboseLevel>0 )
187   {
188     G4cout << "Dingfelder charge decrease model is initialized " << G4endl
189     << "Energy range: "
190     << LowEnergyLimit() / keV << " keV - "
191     << HighEnergyLimit() / MeV << " MeV for "
192     << particle->GetParticleName()
193     << G4endl;
194   }
195 
196   // Initialize water density pointer
197   fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
198 
199   if (isInitialised)
200   { return;}
201   fParticleChangeForGamma = GetParticleChangeForGamma();
202   isInitialised = true;
203 
204 }
205 
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207 
208 G4double G4DNADingfelderChargeDecreaseModel::CrossSectionPerVolume(const G4Material* material,
209                                                                    const G4ParticleDefinition* particleDefinition,
210                                                                    G4double k,
211                                                                    G4double,
212                                                                    G4double)
213 {
214   if (verboseLevel > 3)
215   {
216     G4cout
217         << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel"
218         << G4endl;
219   }
220 
221   // Calculate total cross section for model
222   if (
223       particleDefinition != protonDef
224       &&
225       particleDefinition != alphaPlusPlusDef
226       &&
227       particleDefinition != alphaPlusDef   
228   )
229 
230   return 0;
231 
232   G4double lowLim = 0;
233   G4double highLim = 0;
234   G4double crossSection = 0.;
235 
236   G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
237 
238   const G4String& particleName = particleDefinition->GetParticleName();
239 
240   std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
241   pos1 = lowEnergyLimit.find(particleName);
242 
243   if (pos1 != lowEnergyLimit.end())
244   {
245     lowLim = pos1->second;
246   }
247 
248   std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
249   pos2 = highEnergyLimit.find(particleName);
250 
251   if (pos2 != highEnergyLimit.end())
252   {
253     highLim = pos2->second;
254   }
255 
256   if (k >= lowLim && k <= highLim)
257   {
258     crossSection = Sum(k,particleDefinition);
259   }
260 
261   if (verboseLevel > 2)
262   {
263     G4cout << "_______________________________________" << G4endl;
264     G4cout << "G4DNADingfelderChargeDecreaeModel" << G4endl;
265     G4cout << "Kinetic energy(eV)=" << k/eV << "particle :" << particleName << G4endl;
266     G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
267     G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*
268     waterDensity/(1./cm) << G4endl;
269     // material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
270   }
271 
272   return crossSection*waterDensity;
273   //return crossSection*material->GetAtomicNumDensityVector()[1];
274 
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 
279 void G4DNADingfelderChargeDecreaseModel::SampleSecondaries(std::vector<
280                                                                G4DynamicParticle*>* fvect,
281                                                            const G4MaterialCutsCouple* /*couple*/,
282                                                            const G4DynamicParticle* aDynamicParticle,
283                                                            G4double,
284                                                            G4double)
285 {
286   if (verboseLevel > 3)
287   {
288     G4cout
289         << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel"
290         << G4endl;
291   }
292 
293   G4double inK = aDynamicParticle->GetKineticEnergy();
294 
295   G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
296 
297   G4double particleMass = definition->GetPDGMass();
298 
299   G4int finalStateIndex = RandomSelect(inK,definition);
300 
301   G4int n = NumberOfFinalStates(definition, finalStateIndex);
302   G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
303   G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
304 
305   G4double outK = 0.;
306   
307   if (definition==G4Proton::Proton())
308   {
309    if (!statCode) outK = inK - n*(inK*electron_mass_c2/proton_mass_c2) 
310                              - waterBindingEnergy + outgoingParticleBindingEnergy;
311    else outK = inK;
312   }
313   
314   else
315   {
316    if (!statCode) outK = inK - n*(inK*electron_mass_c2/particleMass) 
317                              - waterBindingEnergy + outgoingParticleBindingEnergy;
318    else outK = inK;
319   }
320   
321   if (outK<0)
322   {
323     G4Exception("G4DNADingfelderChargeDecreaseModel::SampleSecondaries","em0004",
324         FatalException,"Final kinetic energy is negative.");
325   }
326 
327   fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
328 
329   if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(waterBindingEnergy);
330   
331   else
332   
333   {
334    if (definition==G4Proton::Proton()) 
335      fParticleChangeForGamma->ProposeLocalEnergyDeposit(n*(inK*electron_mass_c2/proton_mass_c2) 
336      + waterBindingEnergy - outgoingParticleBindingEnergy);
337    else
338      fParticleChangeForGamma->ProposeLocalEnergyDeposit(n*(inK*electron_mass_c2/particleMass) 
339      + waterBindingEnergy - outgoingParticleBindingEnergy);
340   }
341 
342   auto  dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
343       aDynamicParticle->GetMomentumDirection(),
344       outK);
345   fvect->push_back(dp);
346 
347   const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
348       G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
349           1,
350           theIncomingTrack);
351 
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355 
356 G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
357                                                               G4int finalStateIndex)
358 
359 {
360   if (particleDefinition == G4Proton::Proton())
361     return 1;
362 
363   if (particleDefinition == alphaPlusPlusDef)
364   {
365     if (finalStateIndex == 0)
366       return 1;
367     return 2;
368   }
369 
370   if (particleDefinition == alphaPlusDef)
371     return 1;
372 
373   return 0;
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 
378 G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition,
379                                                                                      G4int finalStateIndex)
380 {
381   if (particleDefinition == G4Proton::Proton())
382     return hydrogenDef;
383 
384   if (particleDefinition == alphaPlusPlusDef)
385   {
386     if (finalStateIndex == 0)
387       return alphaPlusDef;
388     return heliumDef;
389   }
390 
391   if (particleDefinition == alphaPlusDef)
392     return heliumDef;
393 
394   return nullptr;
395 }
396 
397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
398 
399 G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
400                                                                         G4int finalStateIndex)
401 {
402   // Ionization energy of first water shell
403   // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
404   // W + 10.79 eV -> W+ + e-
405 
406   if (particleDefinition == G4Proton::Proton())
407     return 10.79 * eV;
408 
409   if (particleDefinition == alphaPlusPlusDef)
410   {
411     // Binding energy for    W+ -> W++ + e-    10.79 eV
412     // Binding energy for    W  -> W+  + e-    10.79 eV
413     //
414     // Ionization energy of first water shell
415     // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
416 
417     if (finalStateIndex == 0)
418       return 10.79 * eV;
419 
420     return 10.79 * 2 * eV;
421   }
422 
423   if (particleDefinition == alphaPlusDef)
424   {
425     // Binding energy for    W+ -> W++ + e-    10.79 eV
426     // Binding energy for    W  -> W+  + e-    10.79 eV
427     //
428     // Ionization energy of first water shell
429     // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
430 
431     return 10.79 * eV;
432   }
433 
434   return 0;
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438 
439 G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
440                                                                                    G4int finalStateIndex)
441 {
442   if (particleDefinition == G4Proton::Proton())
443     return 13.6 * eV;
444 
445   if (particleDefinition == alphaPlusPlusDef)
446   {
447     // Binding energy for    He+ -> He++ + e-    54.509 eV
448     // Binding energy for    He  -> He+  + e-    24.587 eV
449 
450     if (finalStateIndex == 0)
451       return 54.509 * eV;
452 
453     return (54.509 + 24.587) * eV;
454   }
455 
456   if (particleDefinition == alphaPlusDef)
457   {
458     // Binding energy for    He+ -> He++ + e-    54.509 eV
459     // Binding energy for    He  -> He+  + e-    24.587 eV
460 
461     return 24.587 * eV;
462   }
463 
464   return 0;
465 }
466 
467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
468 
469 G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k,
470                                                                  G4int index,
471                                                                  const G4ParticleDefinition* particleDefinition)
472 {
473   G4int particleTypeIndex = 0;
474 
475   if (particleDefinition == protonDef)
476     particleTypeIndex = 0;
477 
478   if (particleDefinition == alphaPlusPlusDef)
479     particleTypeIndex = 1;
480 
481   if (particleDefinition == alphaPlusDef)
482     particleTypeIndex = 2;
483 
484   //
485   // sigma(T) = f0 10 ^ y(log10(T/eV))
486   //
487   //         /  a0 x + b0                    x < x0
488   //         |
489   // y(x) = <   a0 x + b0 - c0 (x - x0)^d0   x0 <= x < x1
490   //         |
491   //         \  a1 x + b1                    x >= x1
492   //
493   //
494   // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
495   //
496   // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
497   //
498   // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
499   // Inelastic-collision cross sections of liquid water for interactions of energetic proton
500   //
501 
502   if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
503   {
504     //
505     // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
506     //
507     // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
508     //
509     // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
510     //
511 
512     x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
513         + gpow->powA((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
514                        / (c0[index][particleTypeIndex]
515                            * d0[index][particleTypeIndex]),
516                    1. / (d0[index][particleTypeIndex] - 1.));
517     b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
518         - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
519         + b0[index][particleTypeIndex]
520         - c0[index][particleTypeIndex]
521             * gpow->powA(x1[index][particleTypeIndex]
522                            - x0[index][particleTypeIndex],
523                        d0[index][particleTypeIndex]);
524   }
525 
526   G4double x(G4Log(k / eV)/gpow->logZ(10));
527   G4double y;
528 
529   if (x < x0[index][particleTypeIndex])
530     y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
531   else if (x < x1[index][particleTypeIndex])
532     y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
533         - c0[index][particleTypeIndex]
534             * gpow->powA(x - x0[index][particleTypeIndex],
535                        d0[index][particleTypeIndex]);
536   else
537     y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
538 
539   return f0[index][particleTypeIndex] * gpow->powA(10., y) * m * m;
540 
541 }
542 
543 G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
544                                                        const G4ParticleDefinition* particleDefinition)
545 {
546   G4int particleTypeIndex = 0;
547 
548   if (particleDefinition == protonDef)
549     particleTypeIndex = 0;
550 
551   if (particleDefinition == alphaPlusPlusDef)
552     particleTypeIndex = 1;
553 
554   if (particleDefinition == alphaPlusDef)
555     particleTypeIndex = 2;
556 
557   const G4int n = numberOfPartialCrossSections[particleTypeIndex];
558   auto  values(new G4double[n]);
559   G4double value(0);
560   G4int i = n;
561 
562   while (i > 0)
563   {
564     i--;
565     values[i] = PartialCrossSection(k, i, particleDefinition);
566     value += values[i];
567   }
568 
569   value *= G4UniformRand();
570 
571   i = n;
572   while (i > 0)
573   {
574     i--;
575 
576     if (values[i] > value)
577       break;
578 
579     value -= values[i];
580   }
581 
582   delete[] values;
583 
584   return i;
585 }
586 
587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
588 
589 G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k,
590                                                  const G4ParticleDefinition* particleDefinition)
591 {
592   G4int particleTypeIndex = 0;
593 
594   if (particleDefinition == protonDef)
595     particleTypeIndex = 0;
596 
597   if (particleDefinition == alphaPlusPlusDef)
598     particleTypeIndex = 1;
599 
600   if (particleDefinition == alphaPlusDef)
601     particleTypeIndex = 2;
602 
603   G4double totalCrossSection = 0.;
604 
605   for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
606   {
607     totalCrossSection += PartialCrossSection(k, i, particleDefinition);
608   }
609 
610   return totalCrossSection;
611 }
612 
613