source: trunk/source/processes/electromagnetic/utils/src/G4EmModelManager.cc @ 869

Last change on this file since 869 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 21.7 KB
Line 
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// $Id: G4EmModelManager.cc,v 1.40 2007/11/09 11:35:54 vnivanch Exp $
27// GEANT4 tag $Name:  $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4EmModelManager
35//
36// Author:        Vladimir Ivanchenko
37//
38// Creation date: 07.05.2002
39//
40// Modifications:
41//
42// 23-12-02 V.Ivanchenko change interface in order to move
43//                           to cut per region
44// 20-01-03 Migrade to cut per region (V.Ivanchenko)
45// 24-01-03 Make models region aware (V.Ivanchenko)
46// 13-02-03 The set of models is defined for region (V.Ivanchenko)
47// 06-03-03 Fix in energy intervals for models (V.Ivanchenko)
48// 13-04-03 Add startFromNull (V.Ivanchenko)
49// 13-05-03 Add calculation of precise range (V.Ivanchenko)
50// 16-07-03 Replace G4Material by G4MaterialCutCouple in dE/dx and CrossSection
51//          calculation (V.Ivanchenko)
52// 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
53// 03-11-03 Substitute STL vector for G4RegionModels (V.Ivanchenko)
54// 26-01-04 Fix in energy range conditions (V.Ivanchenko)
55// 24-03-05 Remove check or IsInCharge (V.Ivanchenko)
56// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
57// 18-08-05 Fix cut for e+e- pair production (V.Ivanchenko)
58// 29-11-05 Add protection for arithmetic operations with cut=DBL_MAX (V.Ivanchenko)
59// 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
60// 13-05-06 Add GetModel by index method (VI)
61// 15-03-07 Add maxCutInRange (V.Ivanchenko)
62// 12-04-07 Add verbosity at destruction (V.Ivanchenko)
63//
64// Class Description:
65//
66// It is the unified energy loss process it calculates the continuous
67// energy loss for charged particles using a set of Energy Loss
68// models valid for different energy regions. There are a possibility
69// to create and access to dE/dx and range tables, or to calculate
70// that information on fly.
71// -------------------------------------------------------------------
72//
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76#include "G4EmModelManager.hh"
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& list, G4DataVector& lowE)
81{
82  nModelsForRegion      = nMod;
83  theListOfModelIndexes = new G4int [nModelsForRegion];
84  lowKineticEnergy      = new G4double [nModelsForRegion];
85  for (G4int i=0; i<nModelsForRegion; i++) {
86    theListOfModelIndexes[i] = list[i];
87    lowKineticEnergy[i] = lowE[i];
88  }
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
93G4RegionModels::~G4RegionModels()
94{
95  delete [] theListOfModelIndexes;
96  delete [] lowKineticEnergy;
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
101#include "G4LossTableManager.hh"
102#include "G4Step.hh"
103#include "G4ParticleDefinition.hh"
104#include "G4PhysicsVector.hh"
105#include "G4Gamma.hh"
106#include "G4Positron.hh"
107#include "G4MaterialCutsCouple.hh"
108#include "G4ProductionCutsTable.hh"
109#include "G4Region.hh"
110#include "G4RegionStore.hh"
111#include "G4Gamma.hh"
112#include "G4Positron.hh"
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115
116G4EmModelManager::G4EmModelManager():
117  nEmModels(0),
118  nRegions(0),
119  nCouples(0),
120  idxOfRegionModels(0),
121  setOfRegionModels(0),
122  minSubRange(0.1),
123  particle(0),
124  verboseLevel(0)
125{
126  models.clear();
127  flucModels.clear();
128  regions.clear();
129  orderOfModels.clear();
130  upperEkin.clear();
131  maxCutInRange    = 12.*cm;
132  maxSubCutInRange = 0.7*mm;
133  theGamma = G4Gamma::Gamma();
134  thePositron = G4Positron::Positron();
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139G4EmModelManager::~G4EmModelManager()
140{
141  G4int i,j;
142  Clear();
143  if(1 < verboseLevel) {
144    G4cout << "G4EmModelManager:: delete models ";
145    if(particle) G4cout << " for " << particle->GetParticleName();
146    G4cout << " nModels=" << nEmModels <<G4endl;
147  }
148
149  for(i = 0; i<nEmModels; i++) {
150    orderOfModels[i] = 1;
151  }
152  for(i = 0; i<nEmModels; i++) {
153    if (orderOfModels[i]) {
154      orderOfModels[i] = 0;
155      for(j = i+1; j<nEmModels; j++) {
156        if(models[i] == models[j]) orderOfModels[j] = 0;
157      }
158      G4String nam = models[i]->GetName();
159      if(nam != "PAI" && nam != "PAIModel" ) delete models[i];
160    }
161  }
162  for(i = 0; i<nEmModels; i++) {
163    orderOfModels[i] = 1;
164  }
165  for(i = 0; i<nEmModels; i++) {
166    if (orderOfModels[i]) {
167      orderOfModels[i] = 0;
168      for(j = i+1; j<nEmModels; j++) {
169        if(flucModels[i] == flucModels[j]) orderOfModels[j] = 0;
170      }
171      delete flucModels[i];
172    }
173  }
174  if(1 < verboseLevel) 
175    G4cout << "G4EmModelManager:: models are deleted!" << G4endl;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4EmModelManager::Clear()
181{
182  if(1 < verboseLevel) {
183    G4cout << "G4EmModelManager::Clear()";
184    if(particle) G4cout << " for " << particle->GetParticleName();
185    G4cout << G4endl;
186  }
187
188  theCuts.clear();
189  theSubCuts.clear();
190  upperEkin.clear();
191  if(idxOfRegionModels) delete [] idxOfRegionModels;
192  if(setOfRegionModels && nRegions) {
193    for(G4int i=0; i<nRegions; i++) {
194      delete (setOfRegionModels[i]);
195    }
196    delete [] setOfRegionModels;
197  }
198  idxOfRegionModels = 0;
199  setOfRegionModels = 0;
200}
201
202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203
204void G4EmModelManager::AddEmModel(G4int num, G4VEmModel* p,
205                                  G4VEmFluctuationModel* fm, const G4Region* r)
206{
207  if(!p) {
208    G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined." 
209           << G4endl;
210    return;
211  }
212  models.push_back(p);
213  flucModels.push_back(fm);
214  regions.push_back(r);
215  orderOfModels.push_back(num);
216  p->DefineForRegion(r);
217  if (nEmModels>0) {
218    G4int idx = nEmModels;
219    do {idx--;} while (idx && num < orderOfModels[idx]);
220    if (num >= orderOfModels[idx] && num <= orderOfModels[idx+1]) idx++;
221    if (idx < nEmModels) {
222      models[nEmModels] = models[idx];
223      flucModels[nEmModels] = flucModels[idx];
224      regions[nEmModels] = regions[idx];
225      orderOfModels[nEmModels] = orderOfModels[idx];
226      models[idx] = p;
227      flucModels[idx] = fm;
228      regions[idx] = r;
229      orderOfModels[idx] = num;
230    }
231  }
232  nEmModels++;
233}
234
235//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237void G4EmModelManager::UpdateEmModel(const G4String& nam, 
238                                     G4double emin, G4double emax)
239{
240  if (nEmModels) {
241    for(G4int i=0; i<nEmModels; i++) {
242      if(nam == models[i]->GetName()) {
243        models[i]->SetLowEnergyLimit(emin);
244        models[i]->SetHighEnergyLimit(emax);
245        break;
246      }
247    }
248  }
249  G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <"
250         << nam << "> is found out"
251         << G4endl;
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255
256G4VEmModel* G4EmModelManager::GetModel(G4int i)
257{
258  G4VEmModel* m = 0;
259  if(i >= 0 && i < nEmModels) m = models[i];
260  else if(verboseLevel > 0) 
261    G4cout << "G4EmModelManager::GetModel WARNING: "
262           << "index " << i << " is wrong Nmodels= "
263           << nEmModels << G4endl;
264  return m;
265}
266
267//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
268
269const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p,
270                                                 const G4ParticleDefinition* sp,
271                                                       G4double theMinSubRange,
272                                                       G4int val)
273{
274  verboseLevel = val;
275  if(1 < verboseLevel) {
276    G4cout << "G4EmModelManager::Initialise() for "
277           << p->GetParticleName()
278           << G4endl;
279  }
280  // Are models defined?
281  if(!nEmModels) {
282    G4Exception("G4EmModelManager::Initialise without any model defined");
283  }
284  particle = p;
285  secondaryParticle = sp;
286  minSubRange = theMinSubRange;
287
288  Clear();
289  G4RegionStore* regionStore = G4RegionStore::GetInstance();
290  const G4Region* world = 
291    regionStore->GetRegion("DefaultRegionForTheWorld", false);
292
293  // Identify the list of regions with different set of models
294  nRegions = 1;
295  std::vector<const G4Region*> set;
296  set.push_back(world);
297  G4bool isWorld = false;
298
299  for (G4int ii=0; ii<nEmModels; ii++) {
300    const G4Region* r = regions[ii];
301    if ( r == 0 || r == world) {
302      isWorld = true;
303      regions[ii] = world;
304    } else {
305      G4bool newRegion = true;
306      if (nRegions>1) {
307        for (G4int j=1; j<nRegions; j++) {
308          if ( r == set[j] ) newRegion = false;
309        }
310      }
311      if (newRegion) {
312        set.push_back(r);
313        nRegions++;
314      }
315    }
316  }
317
318  G4ProductionCutsTable* theCoupleTable=
319    G4ProductionCutsTable::GetProductionCutsTable();
320  G4int numOfCouples = theCoupleTable->GetTableSize();
321  idxOfRegionModels = new G4int[numOfCouples+1];
322  idxOfRegionModels[numOfCouples] = 0;
323  setOfRegionModels = new G4RegionModels*[nRegions];
324  upperEkin.resize(nEmModels);
325
326  // Order models for regions
327  for (G4int reg=0; reg<nRegions; reg++) {
328    const G4Region* region = set[reg];
329
330    G4int n = 0;
331
332    std::vector<G4int>    modelAtRegion;
333    G4DataVector            eLow;
334    G4DataVector            eHigh;
335    modelAtRegion.clear();
336    eLow.clear();
337    eHigh.clear();
338
339    if(isWorld || 0 < reg) {
340
341      for (G4int ii=0; ii<nEmModels; ii++) {
342
343        G4VEmModel* model = models[ii];
344        if ( region == regions[ii] ) {
345
346          G4double tmin = model->LowEnergyLimit();
347          G4double tmax = model->HighEnergyLimit();
348          if (n>0) tmin = std::max(tmin, eHigh[n-1]);
349
350          if(1 < verboseLevel) {
351            G4cout << "Model #" << ii
352                   << " <" << model->GetName() << "> for region <";
353            if (region) G4cout << region->GetName();
354            G4cout << ">  "
355                 << " tmin(MeV)= " << tmin/MeV
356                 << "; tmax(MeV)= " << tmax/MeV
357                 << "; order= " << orderOfModels[ii]
358                 << G4endl;
359          }
360
361          if (tmin < tmax) {
362            modelAtRegion.push_back(ii);
363            eLow.push_back(tmin);
364            eHigh.push_back(tmax);
365            upperEkin[ii] = tmax;
366            n++;
367          }
368        }
369      }
370    } else {
371      n = 1;
372      models.push_back(0);
373      modelAtRegion.push_back(nEmModels);
374      eLow.push_back(0.0);
375      eHigh.push_back(DBL_MAX);
376      upperEkin.push_back(DBL_MAX);
377    }
378    eLow[0] = 0.0;
379
380    if(1 < verboseLevel) {
381      G4cout << "New G4RegionModels set with " << n << " models for region <";
382      if (region) G4cout << region->GetName();
383      G4cout << ">  Elow(MeV)= ";
384      for(G4int ii=0; ii<n; ii++) {G4cout << eLow[ii]/MeV << " ";}
385      G4cout << G4endl;
386    }
387    G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow);
388    setOfRegionModels[reg] = rm;
389  }
390
391  // Access to materials and build cuts
392
393  for(G4int i=0; i<numOfCouples; i++) {
394
395    const G4MaterialCutsCouple* couple = 
396      theCoupleTable->GetMaterialCutsCouple(i);
397    const G4Material* material = couple->GetMaterial();
398    const G4ProductionCuts* pcuts = couple->GetProductionCuts();
399 
400    G4int reg = nRegions;
401    do {reg--;} while (reg>0 && pcuts != (set[reg]->GetProductionCuts()));
402    idxOfRegionModels[i] = reg;
403
404    if(1 < verboseLevel) {
405      G4cout << "G4EmModelManager::Initialise() for "
406             << material->GetName()
407             << " indexOfCouple= " << i
408             << " indexOfRegion= " << reg
409             << G4endl;
410    }
411
412    G4double cut = DBL_MAX;
413    G4double subcut = DBL_MAX;
414    if(secondaryParticle) {
415      size_t idx = 1;
416      if( secondaryParticle == theGamma ) idx = 0;
417      cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i];
418      if( secondaryParticle == thePositron && cut < DBL_MAX )
419        cut += (*theCoupleTable->GetEnergyCutsVector(2))[i] + 
420               2.0*electron_mass_c2;
421
422      // compute subcut
423      if( cut < DBL_MAX ) subcut = minSubRange*cut;
424      if(pcuts->GetProductionCut(idx) < maxCutInRange) {
425        G4double tcutmax = 
426          theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
427                                               material,maxSubCutInRange);
428        if(tcutmax < subcut) subcut = tcutmax;
429      }
430    }
431
432    G4int nm = setOfRegionModels[reg]->NumberOfModels();
433    for(G4int j=0; j<nm; j++) {
434
435      G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
436
437      G4double tcutmin = model->MinEnergyCut(particle, couple);
438
439      if(cut < tcutmin) cut = tcutmin;
440      if(subcut < tcutmin) subcut = tcutmin;
441      if(1 < verboseLevel) {
442            G4cout << "The model # " << j
443                   << "; tcutmin(MeV)= " << tcutmin/MeV
444                   << "; tcut(MeV)= " << cut/MeV
445                   << "; tsubcut(MeV)= " << subcut/MeV
446                   << " for " << particle->GetParticleName()
447                   << G4endl;
448      }
449    }
450    theCuts.push_back(cut);
451    theSubCuts.push_back(subcut);
452  }
453
454  for(G4int jj=0; jj<nEmModels; jj++) {
455    models[jj]->Initialise(particle, theCuts);
456    if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle);
457  }
458
459
460  if(1 < verboseLevel) {
461    G4cout << "G4EmModelManager for " << particle->GetParticleName() 
462           << " is initialised; nRegions=  " << nRegions
463           << G4endl;
464  }
465
466  return &theCuts;
467}
468
469//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
470
471void G4EmModelManager::FillDEDXVector(G4PhysicsVector* aVector,
472                                      const G4MaterialCutsCouple* couple,
473                                      G4EmTableType tType)
474{
475
476  // vectors to provide continues dE/dx
477  G4DataVector factor;
478  G4DataVector dedxLow;
479  G4DataVector dedxHigh;
480
481  G4double e;
482
483  size_t i = couple->GetIndex();
484  G4double cut = theCuts[i];
485  G4double subcut = 0.0;
486
487  if(fTotal == tType) cut = DBL_MAX;
488  else if(fSubRestricted == tType) subcut = theSubCuts[i];
489
490  if(1 < verboseLevel) {
491    G4cout << "G4EmModelManager::FillDEDXVector() for "
492           << couple->GetMaterial()->GetName()
493           << "  Ecut(MeV)= " << cut
494           << "  Esubcut(MeV)= " << subcut
495           << "  Type " << tType
496           << "  for " << particle->GetParticleName()
497           << G4endl;
498  }
499
500  G4int reg  = idxOfRegionModels[i];
501  const G4RegionModels* regModels = setOfRegionModels[reg];
502  G4int nmod = regModels->NumberOfModels();
503  factor.resize(nmod);
504  dedxLow.resize(nmod);
505  dedxHigh.resize(nmod);
506
507  if(1 < verboseLevel) {
508      G4cout << "There are " << nmod << " models for "
509             << couple->GetMaterial()->GetName()
510             << " at the region #" << reg
511             << G4endl;
512  }
513
514
515  // calculate factors to provide continuity of energy loss
516  factor[0] = 1.0;
517  G4int j;
518
519  G4int totBinsLoss = aVector->GetVectorLength();
520
521  dedxLow[0]  = 0.0;
522
523  e = upperEkin[regModels->ModelIndex(0)];
524  G4VEmModel* model = models[regModels->ModelIndex(0)]; 
525  dedxHigh[0] = 0.0;
526  if(model && cut > subcut) {
527    dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);
528    if(subcut > 0.0) 
529      dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);
530  }
531  if(nmod > 1) {
532    for(j=1; j<nmod; j++) {
533
534      e = upperEkin[regModels->ModelIndex(j-1)];
535      G4int idx = regModels->ModelIndex(j); 
536
537      dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
538      if(subcut > 0.0) 
539        dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
540      if(subcut == cut) dedxLow[j] = 0.0;
541
542      e = upperEkin[idx];
543      dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
544      if(subcut > 0.0) 
545        dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
546      if(subcut == cut) dedxHigh[j] = 0.0;
547    }
548
549    for(j=1; j<nmod; j++) {
550      if(dedxLow[j] > 0.0) factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0);
551      else                 factor[j] = 0.0;
552    }
553
554    if(2 < verboseLevel) {
555      G4cout << "Loop over " << totBinsLoss << " bins start " << G4endl;
556    }
557  }
558
559  // Calculate energy losses vector
560  for(j=0; j<totBinsLoss; j++) {
561
562    G4double e = aVector->GetLowEdgeEnergy(j);
563    G4double fac = 1.0;
564
565      // Choose a model of energy losses
566    G4int k = 0;
567    if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
568      do {
569        k++;
570        fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
571      } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
572    }
573
574    model = models[regModels->ModelIndex(k)];
575    G4double dedx = 0.0;
576    G4double dedx0 = 0.0;
577    if(model && cut > subcut) {
578      dedx = model->ComputeDEDX(couple,particle,e,cut); 
579      dedx0 = dedx;
580      if(subcut > 0.0) dedx -= model->ComputeDEDX(couple,particle,e,subcut); 
581      dedx *= fac;
582    }
583
584    if(dedx < 0.0) dedx = 0.0;
585    if(2 < verboseLevel) {
586        G4cout << "Material= " << couple->GetMaterial()->GetName()
587               << "   E(MeV)= " << e/MeV
588               << "  dEdx(MeV/mm)= " << dedx*mm/MeV
589               << "  dEdx0(MeV/mm)= " << dedx0*mm/MeV
590               << "  fac= " << fac
591               << G4endl;
592    }
593    aVector->PutValue(j, dedx);
594  }
595}
596
597//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
598
599void G4EmModelManager::FillLambdaVector(G4PhysicsVector* aVector,
600                                        const G4MaterialCutsCouple* couple,
601                                        G4bool startFromNull,
602                                        G4EmTableType tType)
603{
604  // vectors to provide continues cross section
605  G4DataVector factor;
606  G4DataVector sigmaLow;
607  G4DataVector sigmaHigh;
608
609  G4double e;
610
611  size_t i = couple->GetIndex();
612  G4double cut  = theCuts[i];
613  G4double tmax = DBL_MAX;
614  if (fSubRestricted == tType) {
615    tmax = cut;
616    cut  = theSubCuts[i];
617  }
618
619  if(1 < verboseLevel) {
620    G4cout << "G4EmModelManager::FillLambdaVector() for particle "
621           << particle->GetParticleName()
622           << " in " << couple->GetMaterial()->GetName()
623           << "  Ecut(MeV)= " << cut
624           << "  Emax(MeV)= " << tmax
625           << "  Type " << tType   
626           << G4endl;
627  }
628
629  G4int reg  = idxOfRegionModels[i];
630  const G4RegionModels* regModels = setOfRegionModels[reg];
631  G4int nmod = regModels->NumberOfModels();
632  factor.resize(nmod);
633  sigmaLow.resize(nmod);
634  sigmaHigh.resize(nmod);
635
636  if(2 < verboseLevel) {
637      G4cout << "There are " << nmod << " models for "
638             << couple->GetMaterial()->GetName() << G4endl;
639  }
640
641  // calculate factors to provide continuity of energy loss
642  factor[0] = 1.0;
643  G4int j;
644  G4int totBinsLambda = aVector->GetVectorLength();
645
646  sigmaLow[0]  = 0.0;
647
648  e = upperEkin[regModels->ModelIndex(0)];
649  G4VEmModel* model = models[regModels->ModelIndex(0)];
650  sigmaHigh[0] = 0.0;
651  if(model) sigmaHigh[0] = model->CrossSection(couple,particle,e,cut,tmax);
652
653  if(2 < verboseLevel) {
654      G4cout << "### For material " << couple->GetMaterial()->GetName()
655             << "  " << nmod
656             << " models"
657             << " Ecut(MeV)= " << cut/MeV
658             << " Emax(MeV)= " << e/MeV
659             << " nbins= "  << totBinsLambda
660             << G4endl;
661      G4cout << " model #0   eUp= " << e
662             << "  sigmaUp= " << sigmaHigh[0] << G4endl;
663  }
664
665
666  if(nmod > 1) {
667
668    for(j=1; j<nmod; j++) {
669
670      e  = upperEkin[regModels->ModelIndex(j-1)];
671      sigmaLow[j] = 
672        models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
673      e  = upperEkin[regModels->ModelIndex(j)];
674      sigmaHigh[j] = 
675        models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
676      if(1 < verboseLevel) {
677        G4cout << " model #" << j << "   eUp= " << e
678               << "  sigmaUp= " << sigmaHigh[j]
679               << "  sigmaDown= " << sigmaLow[j]
680               << G4endl;
681      }
682    }
683    for(j=1; j<nmod; j++) {
684      if(sigmaLow[j] > 0.0) factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0);
685      else                  factor[j] = 0.0;
686    }
687  }
688
689  // Calculate lambda vector
690  for(j=0; j<totBinsLambda; j++) {
691
692    e = aVector->GetLowEdgeEnergy(j);
693
694    // Choose a model of energy losses
695    G4int k = 0;
696    G4double fac = 1.0;
697    if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
698      do {
699        k++;
700        fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
701      } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
702    }
703
704    model = models[regModels->ModelIndex(k)];
705    G4double cross = 0.0;
706    if(model) cross = model->CrossSection(couple,particle,e,cut,tmax)*fac;
707    if(j==0 && startFromNull) cross = 0.0;
708
709    if(2 < verboseLevel) {
710      G4cout << "FillLambdaVector: " << j << ".   e(MeV)= " << e/MeV
711             << "  cross(1/mm)= " << cross*mm
712             << " fac= " << fac << " k= " << k
713             << " model= " << regModels->ModelIndex(k)
714             << G4endl;
715    }
716    if(cross < 0.0) cross = 0.0;
717
718    aVector->PutValue(j, cross);
719  }
720}
721
722//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.