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

update

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r961 r991  
    2424// ********************************************************************
    2525//
    26 // $Id: G4VEnergyLossProcess.cc,v 1.146 2009/02/19 11:25:50 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02-ref-02 $
     26// $Id: G4VEnergyLossProcess.cc,v 1.143 2008/10/17 14:46:16 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    150150  secondaryParticle(0),
    151151  nSCoffRegions(0),
    152   nDERegions(0),
    153152  idxSCoffRegions(0),
    154   idxDERegions(0),
    155153  nProcesses(0),
    156154  theDEDXTable(0),
    … …  
    178176  isIonisation(true),
    179177  useSubCutoff(false),
    180   useDeexcitation(false),
    181178  particle(0),
    182179  currentCouple(0),
    … …  
    240237           << G4endl;
    241238  delete vstrag;
    242   Clean();
     239  Clear();
    243240
    244241  if ( !baseParticle ) {
    … …  
    294291//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    295292
    296 void G4VEnergyLossProcess::Clean()
     293void G4VEnergyLossProcess::Clear()
    297294{
    298295  if(1 < verboseLevel) {
    … …  
    305302  delete [] theCrossSectionMax;
    306303  delete [] idxSCoffRegions;
    307   delete [] idxDERegions;
    308304
    309305  theDEDXAtMaxEnergy = 0;
    … …  
    316312  scProcesses.clear();
    317313  nProcesses = 0;
    318 }
    319 
    320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    321 
    322 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
    323                                                 const G4Material*,
    324                                                 G4double cut)
    325 {
    326   return cut;
    327314}
    328315
    … …  
    377364  }
    378365
    379   Clean();
     366  Clear();
    380367
    381368  // Base particle and set of models can be defined here
    … …  
    420407                                     minSubRange, verboseLevel);
    421408
    422   // Sub Cutoff and Deexcitation
    423   if (nSCoffRegions>0 || nDERegions>0) {
     409  // Sub Cutoff Regime
     410  if (nSCoffRegions>0) {
    424411    theSubCuts = modelManager->SubCutoff();
    425412
    … …  
    427414          G4ProductionCutsTable::GetProductionCutsTable();
    428415    size_t numOfCouples = theCoupleTable->GetTableSize();
    429 
    430     if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
    431     if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
     416    idxSCoffRegions = new G4int[numOfCouples];
    432417 
    433418    for (size_t j=0; j<numOfCouples; j++) {
    … …  
    436421        theCoupleTable->GetMaterialCutsCouple(j);
    437422      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
    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       }
     423      G4int reg = 0;
     424      for(G4int i=0; i<nSCoffRegions; i++) {
     425        if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1;
     426      }
     427      idxSCoffRegions[j] = reg;
    453428    }
    454429  }
    … …  
    470445      }
    471446    }
    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     }
    479447  }
    480448}
    … …  
    511479    if(isIonisation) G4cout << "  isIonisation  flag = 1";
    512480    G4cout << G4endl;
     481  }
     482}
     483
     484//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     485
     486void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
     487{
     488  G4RegionStore* regionStore = G4RegionStore::GetInstance();
     489  if(val) {
     490    useSubCutoff = true;
     491    if (!r) {r = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
     492    if (nSCoffRegions) {
     493      for (G4int i=0; i<nSCoffRegions; i++) {
     494        if (r == scoffRegions[i]) return;
     495      }
     496    }
     497    scoffRegions.push_back(r);
     498    nSCoffRegions++;
     499  } else {
     500    useSubCutoff = false;
    513501  }
    514502}
    … …  
    652640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    653641
    654 void 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 
    729 void 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 
    757 void 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   }
     642G4double G4VEnergyLossProcess::GetContinuousStepLimit(
     643                const G4Track&,
     644                G4double, G4double, G4double&)
     645{
     646  return DBL_MAX;
    781647}
    782648
    … …  
    808674  //  <<" stepLimit= "<<x<<G4endl;
    809675  return x;
     676}
     677
     678//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     679
     680G4double G4VEnergyLossProcess::GetMeanFreePath(
     681                             const G4Track& track,
     682                             G4double,
     683                             G4ForceCondition* condition)
     684
     685{
     686  *condition = NotForced;
     687  return MeanFreePath(track);
    810688}
    811689
    … …  
    928806  if (length >= fRange) {
    929807    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     }
     808    currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
     809                                       eloss, esecdep, length);
    936810    fParticleChange.SetProposedKineticEnergy(0.0);
    937     fParticleChange.ProposeLocalEnergyDeposit(eloss);
     811    fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
    938812    return &fParticleChange;
    939813  }
    … …  
    1060934      G4double tmax =
    1061935        std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
    1062       G4double emean = eloss;
    1063936      eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
    1064                                        tmax,length,emean);
     937                                       tmax,length,eloss);
    1065938      /*                           
    1066939      if(-1 < verboseLevel)
    … …  
    1077950  eloss += esecdep;
    1078951  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   }
    1087952
    1088953  // Energy balanse
    … …  
    12411106
    12421107  SelectModel(postStepScaledEnergy);
    1243   if(useDeexcitation) {
    1244     currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
    1245   }
    1246 
    12471108  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
    12481109  G4double tcut = (*theCuts)[currentMaterialIndex];
    … …  
    12791140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    12801141
    1281 G4bool 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 
    1340 G4bool 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
     1142void G4VEnergyLossProcess::PrintInfoDefinition()
     1143{
     1144  if(0 < verboseLevel) {
     1145    G4cout << G4endl << GetProcessName() << ":   for  "
     1146           << particle->GetParticleName()
     1147           << "    SubType= " << GetProcessSubType()
     1148           << G4endl
     1149           << "      dE/dx and range tables from "
     1150           << G4BestUnit(minKinEnergy,"Energy")
     1151           << " to " << G4BestUnit(maxKinEnergy,"Energy")
     1152           << " in " << nBins << " bins" << G4endl
     1153           << "      Lambda tables from threshold to "
     1154           << G4BestUnit(maxKinEnergy,"Energy")
     1155           << " in " << nBins << " bins, spline: "
     1156           << (G4LossTableManager::Instance())->SplineFlag()
    13511157           << G4endl;
    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 
    1400 G4bool 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 
    1415 G4bool 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 << ">"
     1158    if(theRangeTableForLoss && isIonisation) {
     1159      G4cout << "      finalRange(mm)= " << finalRange/mm
     1160             << ", dRoverRange= " << dRoverRange
     1161             << ", integral: " << integral
     1162             << ", fluct: " << lossFluctuationFlag
     1163             << ", linLossLimit= " << linLossLimit
     1164             << G4endl;
     1165    }
     1166    PrintInfo();
     1167    modelManager->DumpModelList(verboseLevel);
     1168    if(theCSDARangeTable && isIonisation) {
     1169      G4cout << "      CSDA range table up"
     1170             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
     1171             << " in " << nBinsCSDA << " bins" << G4endl;
     1172    }
     1173    if(nSCoffRegions>0 && isIonisation) {
     1174      G4cout << "      Subcutoff sampling in " << nSCoffRegions
     1175             << " regions" << G4endl;
     1176    }
     1177    if(2 < verboseLevel) {
     1178      G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
     1179      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
     1180      G4cout << "non restricted DEDXTable address= "
     1181             << theDEDXunRestrictedTable << G4endl;
     1182      if(theDEDXunRestrictedTable && isIonisation) {
     1183           G4cout << (*theDEDXunRestrictedTable) << G4endl;
     1184      }
     1185      if(theDEDXSubTable && isIonisation) {
     1186        G4cout << (*theDEDXSubTable) << G4endl;
     1187      }
     1188      G4cout << "      CSDARangeTable address= " << theCSDARangeTable
    14351189             << G4endl;
    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"
     1190      if(theCSDARangeTable && isIonisation) {
     1191        G4cout << (*theCSDARangeTable) << G4endl;
     1192      }
     1193      G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss
    14431194             << G4endl;
    1444     }
    1445   }
    1446   return res;
    1447 }
    1448 
    1449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1450 
    1451 G4double 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 
    1469 G4double 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 
    1492 G4double 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 
    1503 G4double 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 
    1513 G4double 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 
    1525 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
    1526                 const G4Track&,
    1527                 G4double, G4double, G4double&)
    1528 {
    1529   return DBL_MAX;
    1530 }
    1531 
    1532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    1533 
    1534 G4PhysicsVector* 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  
    1551 void 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;
     1195      if(theRangeTableForLoss && isIonisation) {
     1196             G4cout << (*theRangeTableForLoss) << G4endl;
     1197      }
     1198      G4cout << "      InverseRangeTable address= " << theInverseRangeTable
     1199             << G4endl;
     1200      if(theInverseRangeTable && isIonisation) {
     1201             G4cout << (*theInverseRangeTable) << G4endl;
     1202      }
     1203      G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
     1204      if(theLambdaTable && isIonisation) {
     1205        G4cout << (*theLambdaTable) << G4endl;
     1206      }
     1207      G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
     1208      if(theSubLambdaTable && isIonisation) {
     1209        G4cout << (*theSubLambdaTable) << G4endl;
     1210      }
    15711211    }
    15721212  }
    … …  
    17371377//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    17381378
     1379G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
     1380                 const G4MaterialCutsCouple* couple, G4double cut)
     1381{
     1382  //  G4double cut  = (*theCuts)[couple->GetIndex()];
     1383  //  G4int nbins = nLambdaBins;
     1384  G4double tmin =
     1385    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
     1386             minKinEnergy);
     1387  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
     1388  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
     1389  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
     1390
     1391  return v;
     1392}
     1393
     1394//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1395
     1396G4double G4VEnergyLossProcess::CrossSectionPerVolume(
     1397         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
     1398{
     1399  // Cross section per volume is calculated
     1400  DefineMaterial(couple);
     1401  G4double cross = 0.0;
     1402  G4bool b;
     1403  if(theLambdaTable) {
     1404    cross =
     1405      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
     1406  } else {
     1407    SelectModel(kineticEnergy);
     1408    cross =
     1409      currentModel->CrossSectionPerVolume(currentMaterial,
     1410                                          particle, kineticEnergy,
     1411                                          (*theCuts)[currentMaterialIndex]);
     1412  }
     1413
     1414  return cross;
     1415}
     1416
     1417//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1418
     1419G4bool G4VEnergyLossProcess::StorePhysicsTable(
     1420       const G4ParticleDefinition* part, const G4String& directory,
     1421       G4bool ascii)
     1422{
     1423  G4bool res = true;
     1424  if ( baseParticle || part != particle ) return res;
     1425
     1426  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
     1427    {res = false;}
     1428
     1429  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
     1430    {res = false;}
     1431
     1432  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
     1433    {res = false;}
     1434
     1435  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
     1436    {res = false;}
     1437
     1438  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
     1439    {res = false;}
     1440
     1441  if(isIonisation &&
     1442     !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
     1443    {res = false;}
     1444
     1445  if(isIonisation &&
     1446     !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
     1447    {res = false;}
     1448 
     1449  if(isIonisation &&
     1450     !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
     1451    {res = false;}
     1452 
     1453  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
     1454    {res = false;}
     1455
     1456  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
     1457    {res = false;}
     1458
     1459  if ( res ) {
     1460    if(0 < verboseLevel) {
     1461      G4cout << "Physics tables are stored for " << particle->GetParticleName()
     1462             << " and process " << GetProcessName()
     1463             << " in the directory <" << directory
     1464             << "> " << G4endl;
     1465    }
     1466  } else {
     1467    G4cout << "Fail to store Physics Tables for "
     1468           << particle->GetParticleName()
     1469           << " and process " << GetProcessName()
     1470           << " in the directory <" << directory
     1471           << "> " << G4endl;
     1472  }
     1473  return res;
     1474}
     1475
     1476//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1477
     1478G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
     1479                                        G4PhysicsTable* aTable, G4bool ascii,
     1480                                        const G4String& directory,
     1481                                        const G4String& tname)
     1482{
     1483  G4bool res = true;
     1484  if ( aTable ) {
     1485    const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
     1486    if( !aTable->StorePhysicsTable(name,ascii)) res = false;
     1487  }
     1488  return res;
     1489}
     1490
     1491//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1492
     1493G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
     1494       const G4ParticleDefinition* part, const G4String& directory,
     1495       G4bool ascii)
     1496{
     1497  G4bool res = true;
     1498  const G4String particleName = part->GetParticleName();
     1499
     1500  if(1 < verboseLevel) {
     1501    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
     1502           << particleName << " and process " << GetProcessName()
     1503           << "; tables_are_built= " << tablesAreBuilt
     1504           << G4endl;
     1505  }
     1506  if(particle == part) {
     1507
     1508    //    G4bool yes = true;
     1509    if ( !baseParticle ) {
     1510
     1511      G4bool fpi = true;
     1512      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
     1513        {fpi = false;}
     1514
     1515      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
     1516        {fpi = false;}
     1517
     1518      if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
     1519        {res = false;}
     1520
     1521      if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
     1522        {res = false;}
     1523
     1524      if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
     1525        {res = false;}
     1526
     1527      if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
     1528        {res = false;}
     1529
     1530      if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
     1531        {res = false;}
     1532
     1533      G4bool yes = false;
     1534      if(nSCoffRegions > 0) {yes = true;}
     1535
     1536      if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
     1537        {res = false;}
     1538
     1539      if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
     1540        {res = false;}
     1541
     1542      if(!fpi) yes = false;
     1543      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
     1544        {res = false;}
     1545    }
     1546  }
     1547
     1548  return res;
     1549}
     1550
     1551//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
     1552
     1553G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
     1554                                           G4PhysicsTable* aTable, G4bool ascii,
     1555                                           const G4String& directory,
     1556                                           const G4String& tname,
     1557                                           G4bool mandatory)
     1558{
     1559  G4bool res = true;
     1560  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
     1561  G4bool yes = aTable->ExistPhysicsTable(filename);
     1562  if(yes) {
     1563    yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
     1564    if((G4LossTableManager::Instance())->SplineFlag()) {
     1565      size_t n = aTable->length();
     1566      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
     1567    }
     1568  }
     1569  if(yes) {
     1570    if (0 < verboseLevel) {
     1571      G4cout << tname << " table for " << part->GetParticleName()
     1572             << " is Retrieved from <" << filename << ">"
     1573             << G4endl;
     1574    }
     1575  } else {
     1576    if(mandatory) res = false;
     1577    if(mandatory || 1 < verboseLevel) {
     1578      G4cout << tname << " table for " << part->GetParticleName()
     1579             << " from file <"
     1580             << filename << "> is not Retrieved"
     1581             << G4endl;
     1582    }
     1583  }
     1584  return res;
     1585}
     1586
     1587//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1588 
     1589void G4VEnergyLossProcess::AddCollaborativeProcess(
     1590            G4VEnergyLossProcess* p)
     1591{
     1592  G4bool add = true;
     1593  if(p->GetProcessName() != "eBrem") add = false;
     1594  if(add && nProcesses > 0) {
     1595    for(G4int i=0; i<nProcesses; i++) {
     1596      if(p == scProcesses[i]) {
     1597        add = false;
     1598        break;
     1599      }
     1600    }
     1601  }
     1602  if(add) {
     1603    scProcesses.push_back(p);
     1604    nProcesses++;
     1605    if (1 < verboseLevel) {
     1606      G4cout << "### The process " << p->GetProcessName()
     1607             << " is added to the list of collaborative processes of "
     1608             << GetProcessName() << G4endl;
     1609    }
     1610  }
     1611}
     1612
     1613//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1614
     1615G4double G4VEnergyLossProcess::GetDEDXDispersion(
     1616                                  const G4MaterialCutsCouple *couple,
     1617                                  const G4DynamicParticle* dp,
     1618                                        G4double length)
     1619{
     1620  DefineMaterial(couple);
     1621  G4double ekin = dp->GetKineticEnergy();
     1622  SelectModel(ekin*massRatio);
     1623  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
     1624  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
     1625  G4double d = 0.0;
     1626  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
     1627  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
     1628  return d;
     1629}
     1630
     1631//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1632
     1633void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*)
     1634{}
     1635
     1636//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     1637
Note: See TracChangeset for help on using the changeset viewer.