Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.123 2008/01/11 19:55:29 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.146 2009/02/19 11:25:50 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    146146
    147147G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name,
    148   G4ProcessType type): G4VContinuousDiscreteProcess(name, type),
     148                                           G4ProcessType type):
     149  G4VContinuousDiscreteProcess(name, type),
     150  secondaryParticle(0),
    149151  nSCoffRegions(0),
     152  nDERegions(0),
    150153  idxSCoffRegions(0),
     154  idxDERegions(0),
    151155  nProcesses(0),
    152156  theDEDXTable(0),
    … …  
    165169  theEnergyOfCrossSectionMax(0),
    166170  theCrossSectionMax(0),
    167   particle(0),
    168171  baseParticle(0),
    169   secondaryParticle(0),
    170   currentCouple(0),
    171   nBins(120),
    172   nBinsCSDA(70),
    173   nWarnings(0),
    174   linLossLimit(0.05),
    175172  minSubRange(0.1),
    176   lambdaFactor(0.8),
    177   mfpKinEnergy(0.0),
    178173  lossFluctuationFlag(true),
    179174  rndmStepFlag(false),
    180175  tablesAreBuilt(false),
    181176  integral(true),
     177  isIon(false),
    182178  isIonisation(true),
    183   useSubCutoff(false)
     179  useSubCutoff(false),
     180  useDeexcitation(false),
     181  particle(0),
     182  currentCouple(0),
     183  nWarnings(0),
     184  mfpKinEnergy(0.0)
    184185{
    185186  SetVerboseLevel(1);
    186187
    187   // Size of tables
     188  // low energy limit
    188189  lowestKinEnergy  = 1.*eV;
     190
     191  // Size of tables assuming spline
    189192  minKinEnergy     = 0.1*keV;
    190193  maxKinEnergy     = 100.0*TeV;
     194  nBins            = 84;
    191195  maxKinEnergyCSDA = 1.0*GeV;
     196  nBinsCSDA        = 35;
     197
     198  // default linear loss limit for spline
     199  linLossLimit  = 0.01;
    192200
    193201  // default dRoverRange and finalRange
    194202  SetStepFunction(0.2, 1.0*mm);
    195203
    196   theElectron = G4Electron::Electron();
    197   thePositron = G4Positron::Positron();
     204  // default lambda factor
     205  lambdaFactor  = 0.8;
     206
     207  // particle types
     208  theElectron   = G4Electron::Electron();
     209  thePositron   = G4Positron::Positron();
     210  theGenericIon = 0;
    198211
    199212  // run time objects
    … …  
    208221  fluctModel = 0;
    209222
    210   scoffRegions.clear();
    211   scProcesses.clear();
    212223  scTracks.reserve(5);
    213224  secParticles.reserve(5);
    … …  
    215226  // Data for stragling of ranges from ICRU'37 report
    216227  const G4int nrbins = 7;
    217   vstrag = new G4PhysicsLogVector(keV, GeV, nrbins);
     228  vstrag = new G4PhysicsLogVector(keV, GeV, nrbins-1);
     229  vstrag->SetSpline(true);
    218230  G4double s[nrbins] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9};
    219231  for(G4int i=0; i<nrbins; i++) {vstrag->PutValue(i, s[i]);}
    … …  
    228240           << G4endl;
    229241  delete vstrag;
    230   Clear();
     242  Clean();
    231243
    232244  if ( !baseParticle ) {
    … …  
    282294//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    283295
    284 void G4VEnergyLossProcess::Clear()
    285 {
    286   if(1 < verboseLevel)
     296void G4VEnergyLossProcess::Clean()
     297{
     298  if(1 < verboseLevel) {
    287299    G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName()
    288     << G4endl;
    289 
     300           << G4endl;
     301  }
    290302  delete [] theDEDXAtMaxEnergy;
    291303  delete [] theRangeAtMaxEnergy;
    … …  
    293305  delete [] theCrossSectionMax;
    294306  delete [] idxSCoffRegions;
     307  delete [] idxDERegions;
    295308
    296309  theDEDXAtMaxEnergy = 0;
    … …  
    300313  tablesAreBuilt = false;
    301314
    302   scTracks.clear();
     315  //scTracks.clear();
    303316  scProcesses.clear();
    304 }
    305 
    306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    307 
    308 void G4VEnergyLossProcess::PreparePhysicsTable(
    309      const G4ParticleDefinition& part)
    310 {
    311 
    312   // Are particle defined?
    313   if( !particle ) {
    314     if(part.GetParticleType() == "nucleus" &&
    315        part.GetParticleSubType() == "generic")
    316          particle = G4GenericIon::GenericIon();
    317     else particle = &part;
    318   }
    319 
     317  nProcesses = 0;
     318}
     319
     320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     321
     322G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
     323                                                const G4Material*,
     324                                                G4double cut)
     325{
     326  return cut;
     327}
     328
     329//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     330
     331void
     332G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
     333{
    320334  if(1 < verboseLevel) {
    321335    G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
    322336           << GetProcessName()
    323337           << " for " << part.GetParticleName()
    324            << " local: " << particle->GetParticleName()
    325338           << G4endl;
    326339  }
    327 
    328   G4LossTableManager* lManager = G4LossTableManager::Instance();
    329 
    330   if(&part != particle) {
    331     if(part.GetParticleType() == "nucleus") lManager->RegisterIon(&part, this);
    332     else                          lManager->RegisterExtraParticle(&part, this);
    333     return;
    334   }
    335 
    336   Clear();
    337340
    338341  currentCouple = 0;
    … …  
    341344  fRange        = DBL_MAX;
    342345  preStepKinEnergy = 0.0;
     346  chargeSqRatio = 1.0;
     347  massRatio = 1.0;
     348  reduceFactor = 1.0;
     349
     350  G4LossTableManager* lManager = G4LossTableManager::Instance();
     351
     352  // Are particle defined?
     353  if( !particle ) {
     354    particle = &part;
     355    if(part.GetParticleType() == "nucleus") {
     356      if(!theGenericIon) theGenericIon = G4GenericIon::GenericIon();
     357      if(particle == theGenericIon) { isIon = true; }
     358      else if(part.GetPDGCharge() > eplus) {
     359        isIon = true;
     360
     361        // generic ions created on-fly
     362        if(part.GetPDGCharge() > 2.5*eplus) {
     363          particle = theGenericIon;
     364        }
     365      }
     366    }
     367  }
     368
     369  if( particle != &part) {
     370    if(part.GetParticleType() == "nucleus") {
     371      isIon = true;
     372      lManager->RegisterIon(&part, this);
     373    } else {
     374      lManager->RegisterExtraParticle(&part, this);
     375    }
     376    return;
     377  }
     378
     379  Clean();
    343380
    344381  // Base particle and set of models can be defined here
    … …  
    372409  G4double initialCharge = particle->GetPDGCharge();
    373410  G4double initialMass   = particle->GetPDGMass();
    374   chargeSquare = initialCharge*initialCharge/(eplus*eplus);
    375   chargeSqRatio = 1.0;
    376   massRatio = 1.0;
    377   reduceFactor = 1.0;
    378411
    379412  if (baseParticle) {
    … …  
    387420                                     minSubRange, verboseLevel);
    388421
    389   // Sub Cutoff Regime
    390   scProcesses.clear();
    391   nProcesses = 0;
    392  
    393   if (nSCoffRegions>0) {
     422  // Sub Cutoff and Deexcitation
     423  if (nSCoffRegions>0 || nDERegions>0) {
    394424    theSubCuts = modelManager->SubCutoff();
    395425
    … …  
    397427          G4ProductionCutsTable::GetProductionCutsTable();
    398428    size_t numOfCouples = theCoupleTable->GetTableSize();
    399     idxSCoffRegions = new G4int[numOfCouples];
     429
     430    if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
     431    if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
    400432 
    401433    for (size_t j=0; j<numOfCouples; j++) {
    … …  
    404436        theCoupleTable->GetMaterialCutsCouple(j);
    405437      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
    406       G4int reg = 0;
    407       for(G4int i=0; i<nSCoffRegions; i++) {
    408         if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1;
    409       }
    410       idxSCoffRegions[j] = reg;
     438     
     439      if(nSCoffRegions>0) {
     440        G4bool reg = false;
     441        for(G4int i=0; i<nSCoffRegions; i++) {
     442          if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
     443        }
     444        idxSCoffRegions[j] = reg;
     445      }
     446      if(nDERegions>0) {
     447        G4bool reg = false;
     448        for(G4int i=0; i<nDERegions; i++) {
     449          if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
     450        }
     451        idxDERegions[j] = reg;
     452      }
    411453    }
    412454  }
    … …  
    416458  if (1 < verboseLevel) {
    417459    G4cout << "G4VEnergyLossProcess::Initialise() is done "
     460           << " for local " << particle->GetParticleName()
     461           << " isIon= " << isIon
    418462           << " chargeSqRatio= " << chargeSqRatio
    419463           << " massRatio= " << massRatio
    … …  
    426470      }
    427471    }
     472    if (nDERegions) {
     473      G4cout << " Deexcitation is ON for regions: " << G4endl;
     474      for (G4int i=0; i<nDERegions; i++) {
     475        const G4Region* r = deRegions[i];
     476        G4cout << "           " << r->GetName() << G4endl;
     477      }
     478    }
    428479  }
    429480}
    … …  
    442493  }
    443494
    444   if(!tablesAreBuilt && &part == particle)
    445     G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
    446 
    447   if(0 < verboseLevel && (&part == particle) && !baseParticle) {
    448     PrintInfoDefinition();
    449     safetyHelper->InitialiseHelper();
     495  if(&part == particle) {
     496    if(!tablesAreBuilt) {
     497      G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
     498    }
     499    if(!baseParticle) {
     500      if(0 < verboseLevel) PrintInfoDefinition();
     501   
     502      // needs to be done only once
     503      safetyHelper->InitialiseHelper();
     504    }
    450505  }
    451506
    … …  
    456511    if(isIonisation) G4cout << "  isIonisation  flag = 1";
    457512    G4cout << G4endl;
    458   }
    459 }
    460 
    461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    462 
    463 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
    464 {
    465   G4RegionStore* regionStore = G4RegionStore::GetInstance();
    466   if(val) {
    467     useSubCutoff = true;
    468     if (!r) r = regionStore->GetRegion("DefaultRegionForTheWorld", false);
    469     if (nSCoffRegions) {
    470       for (G4int i=0; i<nSCoffRegions; i++) {
    471         if (r == scoffRegions[i]) return;
    472       }
    473     }
    474     scoffRegions.push_back(r);
    475     nSCoffRegions++;
    476   } else {
    477     useSubCutoff = false;
    478513  }
    479514}
    … …  
    528563  for(size_t i=0; i<numOfCouples; i++) {
    529564
    530     if(1 < verboseLevel)
     565    if(1 < verboseLevel) {
    531566      G4cout << "G4VEnergyLossProcess::BuildDEDXVector flag=  "
    532567             << table->GetFlag(i) << G4endl;
    533 
     568    }
    534569    if (table->GetFlag(i)) {
    535570
    … …  
    538573        theCoupleTable->GetMaterialCutsCouple(i);
    539574      G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
     575      aVector->SetSpline((G4LossTableManager::Instance())->SplineFlag());
     576
    540577      modelManager->FillDEDXVector(aVector, couple, tType);
    541578
    … …  
    580617           << G4endl;
    581618  }
    582   if(!table) return table;
     619  if(!table) {return table;}
    583620
    584621  // Access to materials
    … …  
    615652//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    616653
    617 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
    618                 const G4Track&,
    619                 G4double, G4double, G4double&)
    620 {
    621   return DBL_MAX;
     654void G4VEnergyLossProcess::PrintInfoDefinition()
     655{
     656  if(0 < verboseLevel) {
     657    G4cout << G4endl << GetProcessName() << ":   for  "
     658           << particle->GetParticleName()
     659           << "    SubType= " << GetProcessSubType()
     660           << G4endl
     661           << "      dE/dx and range tables from "
     662           << G4BestUnit(minKinEnergy,"Energy")
     663           << " to " << G4BestUnit(maxKinEnergy,"Energy")
     664           << " in " << nBins << " bins" << G4endl
     665           << "      Lambda tables from threshold to "
     666           << G4BestUnit(maxKinEnergy,"Energy")
     667           << " in " << nBins << " bins, spline: "
     668           << (G4LossTableManager::Instance())->SplineFlag()
     669           << G4endl;
     670    if(theRangeTableForLoss && isIonisation) {
     671      G4cout << "      finalRange(mm)= " << finalRange/mm
     672             << ", dRoverRange= " << dRoverRange
     673             << ", integral: " << integral
     674             << ", fluct: " << lossFluctuationFlag
     675             << ", linLossLimit= " << linLossLimit
     676             << G4endl;
     677    }
     678    PrintInfo();
     679    modelManager->DumpModelList(verboseLevel);
     680    if(theCSDARangeTable && isIonisation) {
     681      G4cout << "      CSDA range table up"
     682             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
     683             << " in " << nBinsCSDA << " bins" << G4endl;
     684    }
     685    if(nSCoffRegions>0 && isIonisation) {
     686      G4cout << "      Subcutoff sampling in " << nSCoffRegions
     687             << " regions" << G4endl;
     688    }
     689    if(2 < verboseLevel) {
     690      G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
     691      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
     692      G4cout << "non restricted DEDXTable address= "
     693             << theDEDXunRestrictedTable << G4endl;
     694      if(theDEDXunRestrictedTable && isIonisation) {
     695           G4cout << (*theDEDXunRestrictedTable) << G4endl;
     696      }
     697      if(theDEDXSubTable && isIonisation) {
     698        G4cout << (*theDEDXSubTable) << G4endl;
     699      }
     700      G4cout << "      CSDARangeTable address= " << theCSDARangeTable
     701             << G4endl;
     702      if(theCSDARangeTable && isIonisation) {
     703        G4cout << (*theCSDARangeTable) << G4endl;
     704      }
     705      G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
     706             << G4endl;
     707      if(theRangeTableForLoss && isIonisation) {
     708             G4cout << (*theRangeTableForLoss) << G4endl;
     709      }
     710      G4cout << "      InverseRangeTable address= " << theInverseRangeTable
     711             << G4endl;
     712      if(theInverseRangeTable && isIonisation) {
     713             G4cout << (*theInverseRangeTable) << G4endl;
     714      }
     715      G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     716      if(theLambdaTable && isIonisation) {
     717        G4cout << (*theLambdaTable) << G4endl;
     718      }
     719      G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
     720      if(theSubLambdaTable && isIonisation) {
     721        G4cout << (*theSubLambdaTable) << G4endl;
     722      }
     723    }
     724  }
     725}
     726
     727//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     728
     729void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
     730{
     731  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     732  const G4Region* reg = r;
     733  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     734
     735  // the region is in the list
     736  if (nSCoffRegions) {
     737    for (G4int i=0; i<nSCoffRegions; i++) {
     738      if (reg == scoffRegions[i]) {
     739        if(!val) deRegions[i] = 0;
     740        return;
     741      }
     742    }
     743  }
     744
     745  // new region
     746  if(val) {
     747    useSubCutoff = true;
     748    scoffRegions.push_back(reg);
     749    nSCoffRegions++;
     750  } else {
     751    useSubCutoff = false;
     752  }
     753}
     754
     755//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     756
     757void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
     758{
     759  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     760  const G4Region* reg = r;
     761  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     762
     763  // the region is in the list
     764  if (nDERegions) {
     765    for (G4int i=0; i<nDERegions; i++) {
     766      if (reg == deRegions[i]) {
     767        if(!val) deRegions[i] = 0;
     768        return;
     769      }
     770    }
     771  }
     772
     773  // new region
     774  if(val) {
     775    useDeexcitation = true;
     776    deRegions.push_back(reg);
     777    nDERegions++;
     778  } else {
     779    useDeexcitation = false;
     780  }
    622781}
    623782
    … …  
    641800    if(x > finalRange && y < currentMinStep) {
    642801      x = y + finalRange*(1.0 - dRoverRange)*(2.0 - finalRange/fRange);
    643     } else if (rndmStepFlag) x = SampleRange();
    644     //    G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
    645     //    <<" range= "<<fRange <<" cMinSt="<<currentMinStep
    646     //    <<" safety= " << safety<< " limit= " << x <<G4endl;
    647   }
    648   //  G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
     802    } else if (rndmStepFlag) {x = SampleRange();}
     803    //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
     804    //  <<" range= "<<fRange <<" cMinSt="<<currentMinStep
     805    //  << " limit= " << x <<G4endl;
     806  }
     807  //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
    649808  //  <<" stepLimit= "<<x<<G4endl;
    650809  return x;
    651 }
    652 
    653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    654 
    655 G4double G4VEnergyLossProcess::GetMeanFreePath(
    656                              const G4Track& track,
    657                              G4double,
    658                              G4ForceCondition* condition)
    659 
    660 {
    661   *condition = NotForced;
    662   return MeanFreePath(track);
    663810}
    664811
    … …  
    673820  *condition = NotForced;
    674821  G4double x = DBL_MAX;
     822
     823  // initialisation of material, mass, charge, model at the beginning of the step
     824  DefineMaterial(track.GetMaterialCutsCouple());
     825
     826  const G4ParticleDefinition* currPart = track.GetDefinition();
     827  if(theGenericIon == particle) {
     828    massRatio = proton_mass_c2/currPart->GetPDGMass();
     829  } 
     830  preStepKinEnergy    = track.GetKineticEnergy();
     831  preStepScaledEnergy = preStepKinEnergy*massRatio;
     832  SelectModel(preStepScaledEnergy);
     833
     834  if(isIon) {
     835    chargeSqRatio =
     836      currentModel->GetChargeSquareRatio(currPart,currentMaterial,preStepKinEnergy);
     837    reduceFactor  = 1.0/(chargeSqRatio*massRatio);
     838  }
     839  //G4cout << "q2= " << chargeSqRatio << " massRatio= " << massRatio << G4endl;
     840  // initialisation for sampling of the interaction length
    675841  if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
    676   InitialiseStep(track);
    677 
     842  if(theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
     843
     844  // compute mean free path
    678845  if(preStepScaledEnergy < mfpKinEnergy) {
    679846    if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy);
    … …  
    686853    if (theNumberOfInteractionLengthLeft < 0.0) {
    687854      // beggining of tracking (or just after DoIt of this process)
     855      //G4cout<<"G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength Reset"<<G4endl;
    688856      ResetNumberOfInteractionLengthLeft();
    689857    } else if(currentInteractionLength < DBL_MAX) {
    … …  
    702870      G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
    703871      G4cout << "[ " << GetProcessName() << "]" << G4endl;
    704       G4cout << " for " << particle->GetParticleName()
     872      G4cout << " for " << currPart->GetParticleName()
    705873             << " in Material  " <<  currentMaterial->GetName()
    706874             << " Ekin(MeV)= " << preStepKinEnergy/MeV
    … …  
    710878    }
    711879#endif
    712 
    713880    // zero cross section case
    714881  } else {
    … …  
    738905  if(length <= DBL_MIN) return &fParticleChange;
    739906  G4double eloss  = 0.0;
    740 
    741   /*
     907  G4double esecdep = 0.0;
     908 
     909  /* 
    742910  if(-1 < verboseLevel) {
    743911    const G4ParticleDefinition* d = track.GetDefinition();
    … …  
    755923  */
    756924
     925  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
     926
    757927  // stopping
    758928  if (length >= fRange) {
     929    eloss = preStepKinEnergy;
     930    if (useDeexcitation) {
     931      if(idxDERegions[currentMaterialIndex]) {
     932        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     933        if(eloss < 0.0) eloss = 0.0;
     934      }
     935    }
    759936    fParticleChange.SetProposedKineticEnergy(0.0);
    760     fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
     937    fParticleChange.ProposeLocalEnergyDeposit(eloss);
    761938    return &fParticleChange;
    762939  }
    … …  
    772949      GetScaledRangeForScaledEnergy(preStepScaledEnergy) - length/reduceFactor;
    773950    eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
     951   
    774952    /*
    775953    if(-1 < verboseLevel)
    … …  
    781959             << " eloss0(MeV)= "
    782960             << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
     961             << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
    783962             << G4endl;
    784963    */
    785964  }
    786965
    787   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
    788   G4VEmModel* currentModel = SelectModel(preStepScaledEnergy);
    789   /*   
     966  /*   
    790967  G4double eloss0 = eloss;
    791968  if(-1 < verboseLevel ) {
    … …  
    801978  G4double cut  = (*theCuts)[currentMaterialIndex];
    802979  G4double esec = 0.0;
    803   G4double esecdep = 0.0;
    804980
    805981  // SubCutOff
    … …  
    807983    if(idxSCoffRegions[currentMaterialIndex]) {
    808984
    809       G4double currentMinSafety = 0.0;
     985      G4bool yes = false;
    810986      G4StepPoint* prePoint  = step.GetPreStepPoint();
    811       G4StepPoint* postPoint = step.GetPostStepPoint();
    812       G4double preSafety  = prePoint->GetSafety();
    813       G4double postSafety = preSafety - length;
    814       G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
    815 
    816       // recompute safety
    817       if(prePoint->GetStepStatus() != fGeomBoundary &&
    818          postPoint->GetStepStatus() != fGeomBoundary) {
    819 
    820         /*
    821         //      G4bool yes = (track.GetTrackID() == 5512);
    822         G4bool yes = false;
    823         if(yes)
    824           G4cout << "G4VEnergyLoss: presafety= " << preSafety
    825                  << " rcut= " << rcut << "  length= " << length
    826                  << " dir " << track.GetMomentumDirection()
    827                  << G4endl;
    828         */
    829 
    830         if(preSafety < rcut)
     987
     988      // Check boundary
     989      if(prePoint->GetStepStatus() == fGeomBoundary) yes = true;
     990
     991      // Check PrePoint
     992      else {
     993        G4double preSafety  = prePoint->GetSafety();
     994        G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
     995
     996        // recompute presafety
     997        if(preSafety < rcut) {
    831998          preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition());
    832 
    833         //if(yes) {
    834         //  G4cout << "G4VEnergyLoss: newsafety= " << preSafety << G4endl;
    835           //       if(preSafety==0.0 && track.GetTrackID() == 5512 ) exit(1);
    836         //}
    837         if(postSafety < rcut)
    838           postSafety = safetyHelper->ComputeSafety(postPoint->GetPosition());
    839         /*     
    840           if(-1 < verboseLevel)
    841           G4cout << "Subcutoff: presafety(mm)= " << preSafety/mm
    842                  << " postsafety(mm)= " << postSafety/mm
    843                  << " rcut(mm)= " << rcut/mm
    844                  << G4endl;
    845         */
    846         currentMinSafety = std::min(preSafety,postSafety);
    847       }
    848 
     999        }
     1000
     1001        if(preSafety < rcut) yes = true;
     1002
     1003        // Check PostPoint
     1004        else {
     1005          G4double postSafety = preSafety - length;
     1006          if(postSafety < rcut) {
     1007            postSafety =
     1008              safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
     1009            if(postSafety < rcut) yes = true;
     1010          }
     1011        }
     1012      }
     1013 
    8491014      // Decide to start subcut sampling
    850       if(currentMinSafety < rcut) {
     1015      if(yes) {
    8511016
    8521017        cut = (*theSubCuts)[currentMaterialIndex];
    … …  
    8561021                                currentModel,currentMaterialIndex,
    8571022                                esecdep);
     1023        // add bremsstrahlung sampling
    8581024        /*
    8591025        if(nProcesses > 0) {
    … …  
    8761042            esec += e;
    8771043            pParticleChange->AddSecondary(t);
    878             //  mom -= t->GetMomentum();
    8791044          }     
    880           //        fParticleChange.SetProposedMomentum(mom);           
    8811045        }
    8821046      }
    … …  
    8851049
    8861050  // Corrections, which cannot be tabulated
    887   CorrectionsAlongStep(currentCouple, dynParticle, eloss, length);
     1051  currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
     1052                                     eloss, esecdep, length);
    8881053
    8891054  // Sample fluctuations
    … …  
    8951060      G4double tmax =
    8961061        std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
     1062      G4double emean = eloss;
    8971063      eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
    898                                        tmax,length,eloss);
    899       /*           
     1064                                       tmax,length,emean);
     1065      /*                            
    9001066      if(-1 < verboseLevel)
    9011067      G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
    9021068             << " fluc= " << (eloss-eloss0)/MeV
    903              << " currentChargeSquare= " << chargeSquare
     1069             << " ChargeSqRatio= " << chargeSqRatio
    9041070             << " massRatio= " << massRatio
    9051071             << " tmax= " << tmax
    … …  
    9111077  eloss += esecdep;
    9121078  if(eloss < 0.0) eloss = 0.0;
     1079
     1080  // deexcitation
     1081  else if (useDeexcitation) {
     1082    if(idxDERegions[currentMaterialIndex]) {
     1083      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
     1084      if(eloss < 0.0) eloss = 0.0;
     1085    }
     1086  }
    9131087
    9141088  // Energy balanse
    … …  
    9171091    eloss  = preStepKinEnergy - esec;
    9181092    finalT = 0.0;
     1093  } else if(isIon) {
     1094    fParticleChange.SetProposedCharge(
     1095      currentModel->GetParticleCharge(track.GetDefinition(),currentMaterial,finalT));
    9191096  }
    9201097
    … …  
    9221099  fParticleChange.ProposeLocalEnergyDeposit(eloss);
    9231100
    924   /*
     1101  /* 
    9251102  if(-1 < verboseLevel) {
    9261103    G4cout << "Final value eloss(MeV)= " << eloss/MeV
    … …  
    9311108           << G4endl;
    9321109  }
    933   */
     1110  */ 
    9341111
    9351112  return &fParticleChange;
    … …  
    9431120       G4VEmModel* model,
    9441121       G4int idx,
    945        G4double& extraEdep)
     1122       G4double& /*extraEdep*/)
    9461123{
    9471124  // Fast check weather subcutoff can work
    … …  
    9521129  const G4Track* track = step.GetTrack();
    9531130  const G4DynamicParticle* dp = track->GetDynamicParticle();
     1131  G4double e = dp->GetKineticEnergy()*massRatio;
    9541132  G4bool b;
    955   G4double cross =
    956     chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(dp->GetKineticEnergy(),b));
     1133  G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(e,b));
    9571134  G4double length = step.GetStepLength();
    9581135
    9591136  // negligible probability to get any interaction
    9601137  if(length*cross < perMillion) return;
    961   /*   
     1138  /*     
    9621139  if(-1 < verboseLevel)
    9631140    G4cout << "<<< Subcutoff for " << GetProcessName()
    … …  
    9701147  // Sample subcutoff secondaries
    9711148  G4StepPoint* preStepPoint = step.GetPreStepPoint();
     1149  G4StepPoint* postStepPoint = step.GetPostStepPoint();
    9721150  G4ThreeVector prepoint = preStepPoint->GetPosition();
    973   G4ThreeVector dr = step.GetPostStepPoint()->GetPosition() - prepoint;
     1151  G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
    9741152  G4double pretime = preStepPoint->GetGlobalTime();
    975   //  G4double dt = length/preStepPoint->GetVelocity();
     1153  G4double dt = postStepPoint->GetGlobalTime() - pretime;
     1154  //G4double dt = length/preStepPoint->GetVelocity();
    9761155  G4double fragment = 0.0;
    9771156
    … …  
    9921171
    9931172      G4bool addSec = true;
     1173      /*
    9941174      // do not track very low-energy delta-electrons
    9951175      if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
    … …  
    10041184        }
    10051185      }
     1186      */
    10061187      if(addSec) {
    1007         //      G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
    1008         G4Track* t = new G4Track((*it), pretime, r);
     1188        G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
     1189        //G4Track* t = new G4Track((*it), pretime, r);
    10091190        t->SetTouchableHandle(track->GetTouchableHandle());
    10101191        tracks.push_back(t);
    10111192
    1012         /*
    1013           if(-1 < verboseLevel)
    1014           G4cout << "New track " << p->GetDefinition()->GetParticleName()
    1015           << " e(keV)= " << p->GetKineticEnergy()/keV
    1016           << " fragment= " << fragment
    1017           << G4endl;
     1193        /*     
     1194        if(-1 < verboseLevel)
     1195          G4cout << "New track " << t->GetDefinition()->GetParticleName()
     1196                 << " e(keV)= " << t->GetKineticEnergy()/keV
     1197                << " fragment= " << fragment
     1198                << G4endl;
    10181199        */
    10191200      }
    … …  
    10591240  }
    10601241
    1061   G4VEmModel* currentModel = SelectModel(postStepScaledEnergy);
     1242  SelectModel(postStepScaledEnergy);
     1243  if(useDeexcitation) {
     1244    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
     1245  }
     1246
    10621247  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
    10631248  G4double tcut = (*theCuts)[currentMaterialIndex];
    … …  
    10941279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    10951280
    1096 void G4VEnergyLossProcess::PrintInfoDefinition()
    1097 {
    1098   if(0 < verboseLevel) {
    1099     G4cout << G4endl << GetProcessName() << ":   tables are built for  "
    1100            << particle->GetParticleName()
    1101            << G4endl
    1102            << "      dE/dx and range tables from "
    1103            << G4BestUnit(minKinEnergy,"Energy")
    1104            << " to " << G4BestUnit(maxKinEnergy,"Energy")
    1105            << " in " << nBins << " bins." << G4endl
    1106            << "      Lambda tables from threshold to "
    1107            << G4BestUnit(maxKinEnergy,"Energy")
    1108            << " in " << nBins << " bins."
     1281G4bool G4VEnergyLossProcess::StorePhysicsTable(
     1282       const G4ParticleDefinition* part, const G4String& directory,
     1283       G4bool ascii)
     1284{
     1285  G4bool res = true;
     1286  if ( baseParticle || part != particle ) return res;
     1287
     1288  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
     1289    {res = false;}
     1290
     1291  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
     1292    {res = false;}
     1293
     1294  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
     1295    {res = false;}
     1296
     1297  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
     1298    {res = false;}
     1299
     1300  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
     1301    {res = false;}
     1302
     1303  if(isIonisation &&
     1304     !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
     1305    {res = false;}
     1306
     1307  if(isIonisation &&
     1308     !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
     1309    {res = false;}
     1310 
     1311  if(isIonisation &&
     1312     !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
     1313    {res = false;}
     1314 
     1315  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
     1316    {res = false;}
     1317
     1318  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
     1319    {res = false;}
     1320
     1321  if ( res ) {
     1322    if(0 < verboseLevel) {
     1323      G4cout << "Physics tables are stored for " << particle->GetParticleName()
     1324             << " and process " << GetProcessName()
     1325             << " in the directory <" << directory
     1326             << "> " << G4endl;
     1327    }
     1328  } else {
     1329    G4cout << "Fail to store Physics Tables for "
     1330           << particle->GetParticleName()
     1331           << " and process " << GetProcessName()
     1332           << " in the directory <" << directory
     1333           << "> " << G4endl;
     1334  }
     1335  return res;
     1336}
     1337
     1338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1339
     1340G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
     1341       const G4ParticleDefinition* part, const G4String& directory,
     1342       G4bool ascii)
     1343{
     1344  G4bool res = true;
     1345  const G4String particleName = part->GetParticleName();
     1346
     1347  if(1 < verboseLevel) {
     1348    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
     1349           << particleName << " and process " << GetProcessName()
     1350           << "; tables_are_built= " << tablesAreBuilt
    11091351           << G4endl;
    1110     PrintInfo();
    1111     if(theRangeTableForLoss && isIonisation)
    1112       G4cout << "      Step function: finalRange(mm)= " << finalRange/mm
    1113              << ", dRoverRange= " << dRoverRange
    1114              << ", integral: " << integral
    1115              << ", fluct: " << lossFluctuationFlag
    1116              << G4endl;
    1117    
    1118     if(theCSDARangeTable && isIonisation)
    1119       G4cout << "      CSDA range table up"
    1120              << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
    1121              << " in " << nBinsCSDA << " bins." << G4endl;
    1122    
    1123     if(nSCoffRegions>0)
    1124       G4cout << "      Subcutoff sampling in " << nSCoffRegions
    1125              << " regions" << G4endl;
    1126 
    1127     if(2 < verboseLevel) {
    1128       G4cout << "DEDXTable address= " << theDEDXTable << G4endl;
    1129       if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
    1130       G4cout << "non restricted DEDXTable address= "
    1131              << theDEDXunRestrictedTable << G4endl;
    1132       if(theDEDXunRestrictedTable && isIonisation)
    1133            G4cout << (*theDEDXunRestrictedTable) << G4endl;
    1134       if(theDEDXSubTable && isIonisation) G4cout << (*theDEDXSubTable)
    1135                                                  << G4endl;
    1136       G4cout << "CSDARangeTable address= " << theCSDARangeTable
     1352  }
     1353  if(particle == part) {
     1354
     1355    //    G4bool yes = true;
     1356    if ( !baseParticle ) {
     1357
     1358      G4bool fpi = true;
     1359      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
     1360        {fpi = false;}
     1361
     1362      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
     1363        {fpi = false;}
     1364
     1365      if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
     1366        {res = false;}
     1367
     1368      if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
     1369        {res = false;}
     1370
     1371      if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
     1372        {res = false;}
     1373
     1374      if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
     1375        {res = false;}
     1376
     1377      if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
     1378        {res = false;}
     1379
     1380      G4bool yes = false;
     1381      if(nSCoffRegions > 0) {yes = true;}
     1382
     1383      if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
     1384        {res = false;}
     1385
     1386      if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
     1387        {res = false;}
     1388
     1389      if(!fpi) yes = false;
     1390      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
     1391        {res = false;}
     1392    }
     1393  }
     1394
     1395  return res;
     1396}
     1397
     1398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1399
     1400G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
     1401                                        G4PhysicsTable* aTable, G4bool ascii,
     1402                                        const G4String& directory,
     1403                                        const G4String& tname)
     1404{
     1405  G4bool res = true;
     1406  if ( aTable ) {
     1407    const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
     1408    if( !aTable->StorePhysicsTable(name,ascii)) res = false;
     1409  }
     1410  return res;
     1411}
     1412
     1413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1414
     1415G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
     1416                                           G4PhysicsTable* aTable, G4bool ascii,
     1417                                           const G4String& directory,
     1418                                           const G4String& tname,
     1419                                           G4bool mandatory)
     1420{
     1421  G4bool res = true;
     1422  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
     1423  G4bool yes = aTable->ExistPhysicsTable(filename);
     1424  if(yes) {
     1425    yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
     1426    if((G4LossTableManager::Instance())->SplineFlag()) {
     1427      size_t n = aTable->length();
     1428      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
     1429    }
     1430  }
     1431  if(yes) {
     1432    if (0 < verboseLevel) {
     1433      G4cout << tname << " table for " << part->GetParticleName()
     1434             << " is Retrieved from <" << filename << ">"
    11371435             << G4endl;
    1138       if(theCSDARangeTable && isIonisation) G4cout << (*theCSDARangeTable)
    1139             << G4endl;
    1140       G4cout << "RangeTableForLoss address= " << theRangeTableForLoss
     1436    }
     1437  } else {
     1438    if(mandatory) res = false;
     1439    if(mandatory || 1 < verboseLevel) {
     1440      G4cout << tname << " table for " << part->GetParticleName()
     1441             << " from file <"
     1442             << filename << "> is not Retrieved"
    11411443             << G4endl;
    1142       if(theRangeTableForLoss && isIonisation)
    1143              G4cout << (*theRangeTableForLoss) << G4endl;
    1144       G4cout << "InverseRangeTable address= " << theInverseRangeTable
    1145              << G4endl;
    1146       if(theInverseRangeTable && isIonisation)
    1147              G4cout << (*theInverseRangeTable) << G4endl;
    1148       G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
    1149       if(theLambdaTable && isIonisation) G4cout << (*theLambdaTable) << G4endl;
    1150       G4cout << "SubLambdaTable address= " << theSubLambdaTable << G4endl;
    1151       if(theSubLambdaTable && isIonisation) G4cout << (*theSubLambdaTable)
    1152              << G4endl;
     1444    }
     1445  }
     1446  return res;
     1447}
     1448
     1449//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1450
     1451G4double G4VEnergyLossProcess::GetDEDXDispersion(
     1452                                  const G4MaterialCutsCouple *couple,
     1453                                  const G4DynamicParticle* dp,
     1454                                        G4double length)
     1455{
     1456  DefineMaterial(couple);
     1457  G4double ekin = dp->GetKineticEnergy();
     1458  SelectModel(ekin*massRatio);
     1459  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
     1460  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
     1461  G4double d = 0.0;
     1462  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
     1463  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
     1464  return d;
     1465}
     1466
     1467//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1468
     1469G4double G4VEnergyLossProcess::CrossSectionPerVolume(
     1470         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
     1471{
     1472  // Cross section per volume is calculated
     1473  DefineMaterial(couple);
     1474  G4double cross = 0.0;
     1475  G4bool b;
     1476  if(theLambdaTable) {
     1477    cross =
     1478      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
     1479  } else {
     1480    SelectModel(kineticEnergy);
     1481    cross =
     1482      currentModel->CrossSectionPerVolume(currentMaterial,
     1483                                          particle, kineticEnergy,
     1484                                          (*theCuts)[currentMaterialIndex]);
     1485  }
     1486
     1487  return cross;
     1488}
     1489
     1490//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1491
     1492G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
     1493{
     1494  DefineMaterial(track.GetMaterialCutsCouple());
     1495  preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
     1496  G4double x = DBL_MAX;
     1497  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
     1498  return x;
     1499}
     1500
     1501//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1502
     1503G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
     1504                                                   G4double x, G4double y,
     1505                                                   G4double& z)
     1506{
     1507  G4GPILSelection sel;
     1508  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
     1509}
     1510
     1511//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1512
     1513G4double G4VEnergyLossProcess::GetMeanFreePath(
     1514                             const G4Track& track,
     1515                             G4double,
     1516                             G4ForceCondition* condition)
     1517
     1518{
     1519  *condition = NotForced;
     1520  return MeanFreePath(track);
     1521}
     1522
     1523//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1524
     1525G4double G4VEnergyLossProcess::GetContinuousStepLimit(
     1526                const G4Track&,
     1527                G4double, G4double, G4double&)
     1528{
     1529  return DBL_MAX;
     1530}
     1531
     1532//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1533
     1534G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
     1535                 const G4MaterialCutsCouple* couple, G4double cut)
     1536{
     1537  //  G4double cut  = (*theCuts)[couple->GetIndex()];
     1538  //  G4int nbins = nLambdaBins;
     1539  G4double tmin =
     1540    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
     1541             minKinEnergy);
     1542  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
     1543  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
     1544  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
     1545
     1546  return v;
     1547}
     1548
     1549//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1550 
     1551void G4VEnergyLossProcess::AddCollaborativeProcess(
     1552            G4VEnergyLossProcess* p)
     1553{
     1554  G4bool add = true;
     1555  if(p->GetProcessName() != "eBrem") add = false;
     1556  if(add && nProcesses > 0) {
     1557    for(G4int i=0; i<nProcesses; i++) {
     1558      if(p == scProcesses[i]) {
     1559        add = false;
     1560        break;
     1561      }
     1562    }
     1563  }
     1564  if(add) {
     1565    scProcesses.push_back(p);
     1566    nProcesses++;
     1567    if (1 < verboseLevel) {
     1568      G4cout << "### The process " << p->GetProcessName()
     1569             << " is added to the list of collaborative processes of "
     1570             << GetProcessName() << G4endl;
    11531571    }
    11541572  }
    … …  
    11821600  } else if(fSubRestricted == tType) {   
    11831601    theDEDXSubTable = p;
    1184   } else if(fIonisation == tType && theIonisationTable != p) {   
     1602  } else if(fIsIonisation == tType && theIonisationTable != p) {   
    11851603    if(theIonisationTable) theIonisationTable->clearAndDestroy();
    11861604    theIonisationTable = p;
    1187   } else if(fSubIonisation == tType && theIonisationSubTable != p) {   
     1605  } else if(fIsSubIonisation == tType && theIonisationSubTable != p) {   
    11881606    if(theIonisationSubTable) theIonisationSubTable->clearAndDestroy();
    11891607    theIonisationSubTable = p;
    … …  
    13191737//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    13201738
    1321 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
    1322                  const G4MaterialCutsCouple* couple, G4double cut)
    1323 {
    1324   //  G4double cut  = (*theCuts)[couple->GetIndex()];
    1325   //  G4int nbins = nLambdaBins;
    1326   G4double tmin =
    1327     std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
    1328              minKinEnergy);
    1329   if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
    1330   G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
    1331   return v;
    1332 }
    1333 
    1334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1335 
    1336 G4double G4VEnergyLossProcess::MicroscopicCrossSection(
    1337          G4double kineticEnergy, const G4MaterialCutsCouple* couple)
    1338 {
    1339   // Cross section per atom is calculated
    1340   DefineMaterial(couple);
    1341   G4double cross = 0.0;
    1342   G4bool b;
    1343   if(theLambdaTable)
    1344     cross =
    1345       ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b)/
    1346       currentMaterial->GetTotNbOfAtomsPerVolume();
    1347 
    1348   return cross;
    1349 }
    1350 
    1351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1352 
    1353 G4bool G4VEnergyLossProcess::StorePhysicsTable(
    1354        const G4ParticleDefinition* part, const G4String& directory,
    1355        G4bool ascii)
    1356 {
    1357   G4bool res = true;
    1358   if ( baseParticle || part != particle ) return res;
    1359 
    1360   if ( theDEDXTable ) {
    1361     const G4String name = GetPhysicsTableFileName(part,directory,"DEDX",ascii);
    1362     if( !theDEDXTable->StorePhysicsTable(name,ascii)) res = false;
    1363   }
    1364 
    1365   if ( theDEDXunRestrictedTable ) {
    1366     const G4String name =
    1367       GetPhysicsTableFileName(part,directory,"DEDXnr",ascii);
    1368     if( !theDEDXTable->StorePhysicsTable(name,ascii)) res = false;
    1369   }
    1370 
    1371   if ( theDEDXSubTable ) {
    1372     const G4String name =
    1373       GetPhysicsTableFileName(part,directory,"SubDEDX",ascii);
    1374     if( !theDEDXSubTable->StorePhysicsTable(name,ascii)) res = false;
    1375   }
    1376 
    1377   if ( theIonisationTable ) {
    1378     const G4String name =
    1379       GetPhysicsTableFileName(part,directory,"Ionisation",ascii);
    1380     if( !theIonisationTable->StorePhysicsTable(name,ascii)) res = false;
    1381   }
    1382 
    1383   if ( theIonisationSubTable ) {
    1384     const G4String name =
    1385       GetPhysicsTableFileName(part,directory,"SubIonisation",ascii);
    1386     if( !theIonisationSubTable->StorePhysicsTable(name,ascii)) res = false;
    1387   }
    1388 
    1389   if ( theCSDARangeTable && isIonisation ) {
    1390     const G4String name =
    1391       GetPhysicsTableFileName(part,directory,"CSDARange",ascii);
    1392     if( !theCSDARangeTable->StorePhysicsTable(name,ascii)) res = false;
    1393   }
    1394 
    1395   if ( theRangeTableForLoss && isIonisation ) {
    1396     const G4String name =
    1397       GetPhysicsTableFileName(part,directory,"Range",ascii);
    1398     if( !theRangeTableForLoss->StorePhysicsTable(name,ascii)) res = false;
    1399   }
    1400 
    1401   if ( theInverseRangeTable && isIonisation ) {
    1402     const G4String name =
    1403       GetPhysicsTableFileName(part,directory,"InverseRange",ascii);
    1404     if( !theInverseRangeTable->StorePhysicsTable(name,ascii)) res = false;
    1405   }
    1406 
    1407   if ( theLambdaTable  && isIonisation) {
    1408     const G4String name =
    1409       GetPhysicsTableFileName(part,directory,"Lambda",ascii);
    1410     if( !theLambdaTable->StorePhysicsTable(name,ascii)) res = false;
    1411   }
    1412 
    1413   if ( theSubLambdaTable  && isIonisation) {
    1414     const G4String name =
    1415       GetPhysicsTableFileName(part,directory,"SubLambda",ascii);
    1416     if( !theSubLambdaTable->StorePhysicsTable(name,ascii)) res = false;
    1417   }
    1418   if ( res ) {
    1419     if(0 < verboseLevel) {
    1420       G4cout << "Physics tables are stored for " << particle->GetParticleName()
    1421              << " and process " << GetProcessName()
    1422              << " in the directory <" << directory
    1423              << "> " << G4endl;
    1424     }
    1425   } else {
    1426     G4cout << "Fail to store Physics Tables for "
    1427            << particle->GetParticleName()
    1428            << " and process " << GetProcessName()
    1429            << " in the directory <" << directory
    1430            << "> " << G4endl;
    1431   }
    1432   return res;
    1433 }
    1434 
    1435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
    1436 
    1437 G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
    1438        const G4ParticleDefinition* part, const G4String& directory,
    1439        G4bool ascii)
    1440 {
    1441   G4bool res = true;
    1442   const G4String particleName = part->GetParticleName();
    1443 
    1444   if(1 < verboseLevel) {
    1445     G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
    1446            << particleName << " and process " << GetProcessName()
    1447            << "; tables_are_built= " << tablesAreBuilt
    1448            << G4endl;
    1449   }
    1450   if(particle == part) {
    1451 
    1452     G4bool yes = true;
    1453     G4bool fpi = true;
    1454     if ( !baseParticle ) {
    1455       G4String filename;
    1456 
    1457       filename = GetPhysicsTableFileName(part,directory,"DEDX",ascii);
    1458       yes = theDEDXTable->ExistPhysicsTable(filename);
    1459       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1460                     theDEDXTable,filename,ascii);
    1461       if(yes) {
    1462         if (0 < verboseLevel) {
    1463           G4cout << "DEDX table for " << particleName
    1464                  << " is Retrieved from <"
    1465                  << filename << ">"
    1466                  << G4endl;
    1467         }
    1468       } else {
    1469         fpi = false;
    1470         if (1 < verboseLevel) {
    1471           G4cout << "DEDX table for " << particleName << " from file <"
    1472                  << filename << "> is not Retrieved"
    1473                  << G4endl;
    1474         }
    1475       }
    1476 
    1477       filename = GetPhysicsTableFileName(part,directory,"Range",ascii);
    1478       yes = theRangeTableForLoss->ExistPhysicsTable(filename);
    1479       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1480                     theRangeTableForLoss,filename,ascii);
    1481       if(yes) {
    1482         if (0 < verboseLevel) {
    1483           G4cout << "Range table for loss for " << particleName
    1484                  << " is Retrieved from <"
    1485                  << filename << ">"
    1486                  << G4endl;
    1487         }
    1488       } else {
    1489         if(fpi) {
    1490           res = false;
    1491           G4cout << "Range table for loss for " << particleName
    1492                  << " from file <"
    1493                  << filename << "> is not Retrieved"
    1494                  << G4endl;
    1495         }
    1496       }
    1497 
    1498       filename = GetPhysicsTableFileName(part,directory,"DEDXnr",ascii);
    1499       yes = theDEDXunRestrictedTable->ExistPhysicsTable(filename);
    1500       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1501                     theDEDXunRestrictedTable,filename,ascii);
    1502       if(yes) {
    1503         if (0 < verboseLevel) {
    1504           G4cout << "Non-restricted DEDX table for " << particleName
    1505                  << " is Retrieved from <"
    1506                  << filename << ">"
    1507                  << G4endl;
    1508         }
    1509       } else {
    1510         if (1 < verboseLevel) {
    1511           G4cout << "Non-restricted DEDX table for " << particleName
    1512                  << " from file <"
    1513                  << filename << "> is not Retrieved"
    1514                  << G4endl;
    1515         }
    1516       }
    1517 
    1518       filename = GetPhysicsTableFileName(part,directory,"CSDARange",ascii);
    1519       yes = theCSDARangeTable->ExistPhysicsTable(filename);
    1520       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1521                     theCSDARangeTable,filename,ascii);
    1522       if(yes) {
    1523         if (0 < verboseLevel) {
    1524           G4cout << "CSDA Range table for " << particleName
    1525                  << " is Retrieved from <"
    1526                  << filename << ">"
    1527                  << G4endl;
    1528         }
    1529       } else {
    1530         G4cout << "CSDA Range table for loss for " << particleName
    1531                << " does not exist"
    1532                << G4endl;
    1533       }
    1534 
    1535       filename = GetPhysicsTableFileName(part,directory,"InverseRange",ascii);
    1536       yes = theInverseRangeTable->ExistPhysicsTable(filename);
    1537       if(yes)  yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1538                      theInverseRangeTable,filename,ascii);
    1539       if(yes) {
    1540         if (0 < verboseLevel) {
    1541           G4cout << "InverseRange table for " << particleName
    1542                  << " is Retrieved from <"
    1543                  << filename << ">"
    1544                  << G4endl;
    1545         }
    1546       } else {
    1547         if(fpi) {
    1548           res = false;
    1549           G4cout << "InverseRange table for " << particleName
    1550                  << " from file <"
    1551                  << filename << "> is not Retrieved"
    1552                  << G4endl;
    1553 
    1554         }
    1555       }
    1556 
    1557       filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
    1558       yes = theLambdaTable->ExistPhysicsTable(filename);
    1559       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1560                     theLambdaTable,filename,ascii);
    1561       if(yes) {
    1562         if (0 < verboseLevel) {
    1563           G4cout << "Lambda table for " << particleName
    1564                  << " is Retrieved from <"
    1565                  << filename << ">"
    1566                  << G4endl;
    1567         }
    1568       } else {
    1569         if(fpi) {
    1570           res = false;
    1571           G4cout << "Lambda table for " << particleName << " from file <"
    1572                  << filename << "> is not Retrieved"
    1573                  << G4endl;
    1574         }
    1575       }
    1576 
    1577       filename = GetPhysicsTableFileName(part,directory,"SubDEDX",ascii);
    1578       yes = theDEDXSubTable->ExistPhysicsTable(filename);
    1579       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1580                     theDEDXSubTable,filename,ascii);
    1581       if(yes) {
    1582         if (0 < verboseLevel) {
    1583           G4cout << "SubDEDX table for " << particleName
    1584                  << " is Retrieved from <"
    1585                  << filename << ">"
    1586                  << G4endl;
    1587         }
    1588       } else {
    1589         if(nSCoffRegions) {
    1590           res=false;
    1591           G4cout << "SubDEDX table for " << particleName << " from file <"
    1592                  << filename << "> is not Retrieved"
    1593                  << G4endl;
    1594         }
    1595       }
    1596 
    1597       filename = GetPhysicsTableFileName(part,directory,"SubLambda",ascii);
    1598       yes = theSubLambdaTable->ExistPhysicsTable(filename);
    1599       if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1600                     theSubLambdaTable,filename,ascii);
    1601       if(yes) {
    1602         if (0 < verboseLevel) {
    1603           G4cout << "SubLambda table for " << particleName
    1604                  << " is Retrieved from <"
    1605                  << filename << ">"
    1606                  << G4endl;
    1607         }
    1608       } else {
    1609         if(nSCoffRegions) {
    1610           res=false;
    1611           G4cout << "SubLambda table for " << particleName << " from file <"
    1612                  << filename << "> is not Retrieved"
    1613                  << G4endl;
    1614         }
    1615       }
    1616 
    1617       filename = GetPhysicsTableFileName(part,directory,"Ionisation",ascii);
    1618       yes = theIonisationTable->ExistPhysicsTable(filename);
    1619       if(yes) {
    1620         theIonisationTable =
    1621           G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable);
    1622        
    1623         yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1624                     theIonisationTable,filename,ascii);
    1625       }
    1626       if(yes) {
    1627         if (0 < verboseLevel) {
    1628           G4cout << "Ionisation table for " << particleName
    1629                  << " is Retrieved from <"
    1630                  << filename << ">"
    1631                  << G4endl;
    1632         }
    1633       }
    1634 
    1635       filename = GetPhysicsTableFileName(part,directory,"SubIonisation",ascii);
    1636       yes = theIonisationSubTable->ExistPhysicsTable(filename);
    1637       if(yes) {
    1638         theIonisationSubTable =
    1639           G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable);
    1640         yes = G4PhysicsTableHelper::RetrievePhysicsTable(
    1641                     theIonisationSubTable,filename,ascii);
    1642       }
    1643       if(yes) {
    1644         if (0 < verboseLevel) {
    1645           G4cout << "SubIonisation table for " << particleName
    1646                  << " is Retrieved from <"
    1647                  << filename << ">"
    1648                  << G4endl;
    1649         }
    1650       }
    1651     }
    1652   }
    1653 
    1654   return res;
    1655 }
    1656 
    1657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1658  
    1659 void G4VEnergyLossProcess::AddCollaborativeProcess(
    1660             G4VEnergyLossProcess* p)
    1661 {
    1662   G4bool add = true;
    1663   if(nProcesses > 0) {
    1664     for(G4int i=0; i<nProcesses; i++) {
    1665       if(p == scProcesses[i]) {
    1666         add = false;
    1667         break;
    1668       }
    1669     }
    1670   }
    1671   if(add) {
    1672     scProcesses.push_back(p);
    1673     nProcesses++;
    1674     if (0 < verboseLevel)
    1675       G4cout << "### The process " << p->GetProcessName()
    1676              << " is added to the list of collaborative processes of "
    1677              << GetProcessName() << G4endl;
    1678   }
    1679 }
    1680 
    1681 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1682 
    1683 G4double G4VEnergyLossProcess::GetDEDXDispersion(
    1684                                   const G4MaterialCutsCouple *couple,
    1685                                   const G4DynamicParticle* dp,
    1686                                         G4double length)
    1687 {
    1688   DefineMaterial(couple);
    1689   G4double ekin = dp->GetKineticEnergy();
    1690   G4VEmModel* currentModel = SelectModel(ekin*massRatio);
    1691   G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
    1692   tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
    1693   G4double d = 0.0;
    1694   G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
    1695   if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
    1696   return d;
    1697 }
    1698 
    1699 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1700 
    1701 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*)
    1702 {}
    1703 
    1704 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1705 
Note: See TracChangeset for help on using the changeset viewer.