00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <fstream>
00030 #include <iomanip>
00031
00032 #include "G4AdjointCSManager.hh"
00033
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4AdjointCSMatrix.hh"
00037 #include "G4AdjointInterpolator.hh"
00038 #include "G4AdjointCSMatrix.hh"
00039 #include "G4VEmAdjointModel.hh"
00040 #include "G4ElementTable.hh"
00041 #include "G4Element.hh"
00042 #include "G4ParticleDefinition.hh"
00043 #include "G4Element.hh"
00044 #include "G4VEmProcess.hh"
00045 #include "G4VEnergyLossProcess.hh"
00046 #include "G4PhysicsTable.hh"
00047 #include "G4PhysicsLogVector.hh"
00048 #include "G4PhysicsTableHelper.hh"
00049 #include "G4Electron.hh"
00050 #include "G4Gamma.hh"
00051 #include "G4Proton.hh"
00052 #include "G4AdjointElectron.hh"
00053 #include "G4AdjointGamma.hh"
00054 #include "G4AdjointProton.hh"
00055 #include "G4ProductionCutsTable.hh"
00056 #include "G4ProductionCutsTable.hh"
00057
00058 G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
00060
00061 G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
00062 { if(theInstance == 0) {
00063 static G4AdjointCSManager ins;
00064 theInstance = &ins;
00065 }
00066 return theInstance;
00067 }
00068
00070
00071 G4AdjointCSManager::G4AdjointCSManager()
00072 { CrossSectionMatrixesAreBuilt=false;
00073 TotalSigmaTableAreBuilt=false;
00074 theTotalForwardSigmaTableVector.clear();
00075 theTotalAdjointSigmaTableVector.clear();
00076 listOfForwardEmProcess.clear();
00077 listOfForwardEnergyLossProcess.clear();
00078 theListOfAdjointParticlesInAction.clear();
00079 EminForFwdSigmaTables.clear();
00080 EminForAdjSigmaTables.clear();
00081 EkinofFwdSigmaMax.clear();
00082 EkinofAdjSigmaMax.clear();
00083 listSigmaTableForAdjointModelScatProjToProj.clear();
00084 listSigmaTableForAdjointModelProdToProj.clear();
00085 Tmin=0.1*keV;
00086 Tmax=100.*TeV;
00087 nbins=320;
00088
00089 RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
00090 RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
00091 RegisterAdjointParticle(G4AdjointProton::AdjointProton());
00092
00093 verbose = 1;
00094
00095 lastPartDefForCS =0;
00096 LastEkinForCS =0;
00097 LastCSCorrectionFactor =1.;
00098
00099 forward_CS_mode = true;
00100
00101 currentParticleDef = 0;
00102 currentCouple =0;
00103 currentMaterial=0;
00104 lastMaterial=0;
00105
00106
00107 theAdjIon = 0;
00108 theFwdIon = 0;
00109
00110
00111
00112
00113 }
00115
00116 G4AdjointCSManager::~G4AdjointCSManager()
00117 {;
00118 }
00120
00121 size_t G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
00122 {listOfAdjointEMModel.push_back(aModel);
00123 listSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable);
00124 listSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable);
00125 return listOfAdjointEMModel.size() -1;
00126
00127 }
00129
00130 void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
00131 {
00132 G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
00133 if (anAdjPartDef && aProcess){
00134 RegisterAdjointParticle(anAdjPartDef);
00135 G4int index=-1;
00136
00137 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
00138 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
00139 }
00140 listOfForwardEmProcess[index]->push_back(aProcess);
00141 }
00142 }
00144
00145 void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
00146 {
00147 G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
00148 if (anAdjPartDef && aProcess){
00149 RegisterAdjointParticle(anAdjPartDef);
00150 G4int index=-1;
00151 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
00152 if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
00153 }
00154 listOfForwardEnergyLossProcess[index]->push_back(aProcess);
00155 }
00156 }
00158
00159 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
00160 { G4int index=-1;
00161 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
00162 if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
00163 }
00164
00165 if (index ==-1){
00166 listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>());
00167 theTotalForwardSigmaTableVector.push_back(new G4PhysicsTable);
00168 theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
00169 listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
00170 theListOfAdjointParticlesInAction.push_back(aPartDef);
00171 EminForFwdSigmaTables.push_back(std::vector<G4double> ());
00172 EminForAdjSigmaTables.push_back(std::vector<G4double> ());
00173 EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
00174 EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
00175
00176 }
00177 }
00179
00180 void G4AdjointCSManager::BuildCrossSectionMatrices()
00181 {
00182 if (CrossSectionMatrixesAreBuilt) return;
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 theAdjointCSMatricesForScatProjToProj.clear();
00203 theAdjointCSMatricesForProdToProj.clear();
00204 const G4ElementTable* theElementTable = G4Element::GetElementTable();
00205 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00206
00207 G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
00208 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
00209 G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
00210 G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
00211 if (aModel->GetUseMatrix()){
00212 std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
00213 std::vector<G4AdjointCSMatrix*>* aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
00214 aListOfMat1->clear();
00215 aListOfMat2->clear();
00216 if (aModel->GetUseMatrixPerElement()){
00217 if (aModel->GetUseOnlyOneMatrixForAllElements()){
00218 std::vector<G4AdjointCSMatrix*>
00219 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
00220 aListOfMat1->push_back(two_matrices[0]);
00221 aListOfMat2->push_back(two_matrices[1]);
00222 }
00223 else {
00224 for (size_t j=0; j<theElementTable->size();j++){
00225 G4Element* anElement=(*theElementTable)[j];
00226 G4int Z = int(anElement->GetZ());
00227 G4int A = int(anElement->GetA());
00228 std::vector<G4AdjointCSMatrix*>
00229 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
00230 aListOfMat1->push_back(two_matrices[0]);
00231 aListOfMat2->push_back(two_matrices[1]);
00232 }
00233 }
00234 }
00235 else {
00236 for (size_t j=0; j<theMaterialTable->size();j++){
00237 G4Material* aMaterial=(*theMaterialTable)[j];
00238 std::vector<G4AdjointCSMatrix*>
00239 two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
00240 aListOfMat1->push_back(two_matrices[0]);
00241 aListOfMat2->push_back(two_matrices[1]);
00242 }
00243
00244 }
00245 theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
00246 theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
00247 aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
00248 }
00249 else { G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
00250 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
00251 theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
00252 theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
00253
00254 }
00255 }
00256 G4cout<<" All adjoint cross section matrices are computed!"<<G4endl;
00257 G4cout<<"======================================================================"<<G4endl;
00258
00259 CrossSectionMatrixesAreBuilt = true;
00260
00261
00262 }
00263
00264
00266
00267 void G4AdjointCSManager::BuildTotalSigmaTables()
00268 { if (TotalSigmaTableAreBuilt) return;
00269
00270
00271 const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
00272
00273
00274
00275 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
00276 listSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy();
00277 listSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy();
00278 for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
00279 listSigmaTableForAdjointModelScatProjToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
00280 listSigmaTableForAdjointModelProdToProj[i]->push_back(new G4PhysicsLogVector(Tmin, Tmax, nbins));
00281 }
00282 }
00283
00284
00285
00286 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
00287 G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
00288 DefineCurrentParticle(thePartDef);
00289 theTotalForwardSigmaTableVector[i]->clearAndDestroy();
00290 theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
00291 EminForFwdSigmaTables[i].clear();
00292 EminForAdjSigmaTables[i].clear();
00293 EkinofFwdSigmaMax[i].clear();
00294 EkinofAdjSigmaMax[i].clear();
00295
00296
00297 for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
00298 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317 G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins);
00318 G4bool Emin_found=false;
00319 G4double sigma_max =0.;
00320 G4double e_sigma_max =0.;
00321 for(size_t l=0; l<aVector->GetVectorLength(); l++) {
00322 G4double totCS=0.;
00323 G4double e=aVector->GetLowEdgeEnergy(l);
00324 for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
00325 totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
00326 }
00327 for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
00328 if (thePartDef == theAdjIon) {
00329 size_t mat_index = couple->GetIndex();
00330 G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
00331 G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
00332 (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
00333 }
00334 G4double e1=e/massRatio;
00335 totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
00336 }
00337 aVector->PutValue(l,totCS);
00338 if (totCS>sigma_max){
00339 sigma_max=totCS;
00340 e_sigma_max = e;
00341
00342 }
00343
00344
00345 if (totCS>0 && !Emin_found) {
00346 EminForFwdSigmaTables[i].push_back(e);
00347 Emin_found=true;
00348 }
00349
00350
00351 }
00352
00353
00354 EkinofFwdSigmaMax[i].push_back(e_sigma_max);
00355
00356
00357 if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
00358
00359 theTotalForwardSigmaTableVector[i]->push_back(aVector);
00360
00361
00362 Emin_found=false;
00363 sigma_max=0;
00364 e_sigma_max =0.;
00365 G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins);
00366 for(eindex=0; eindex<aVector->GetVectorLength(); eindex++) {
00367 G4double e=aVector->GetLowEdgeEnergy(eindex);
00368 G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio);
00369 aVector1->PutValue(eindex,totCS);
00370 if (totCS>sigma_max){
00371 sigma_max=totCS;
00372 e_sigma_max = e;
00373
00374 }
00375
00376 if (totCS>0 && !Emin_found) {
00377 EminForAdjSigmaTables[i].push_back(e);
00378 Emin_found=true;
00379 }
00380
00381 }
00382
00383 EkinofAdjSigmaMax[i].push_back(e_sigma_max);
00384 if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
00385
00386 theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
00387
00388 }
00389 }
00390 TotalSigmaTableAreBuilt =true;
00391
00392 }
00394
00395 G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
00396 const G4MaterialCutsCouple* aCouple)
00397 { DefineCurrentMaterial(aCouple);
00398 DefineCurrentParticle(aPartDef);
00399 G4bool b;
00400 return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
00401
00402
00403
00404 }
00406
00407 G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
00408 const G4MaterialCutsCouple* aCouple)
00409 { DefineCurrentMaterial(aCouple);
00410 DefineCurrentParticle(aPartDef);
00411 G4bool b;
00412 return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
00413
00414
00415 }
00417
00418 G4double G4AdjointCSManager::GetAdjointSigma(G4double Ekin_nuc, size_t index_model,G4bool is_scat_proj_to_proj,
00419 const G4MaterialCutsCouple* aCouple)
00420 { DefineCurrentMaterial(aCouple);
00421 G4bool b;
00422 if (is_scat_proj_to_proj) return (((*listSigmaTableForAdjointModelScatProjToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
00423 else return (((*listSigmaTableForAdjointModelProdToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
00424 }
00426
00427 void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
00428 const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
00429 { DefineCurrentMaterial(aCouple);
00430 DefineCurrentParticle(aPartDef);
00431 emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
00432 emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
00433
00434
00435
00436 }
00438
00439 void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
00440 const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
00441 { DefineCurrentMaterial(aCouple);
00442 DefineCurrentParticle(aPartDef);
00443 e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
00444 G4bool b;
00445 sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
00446 e_sigma_max/=massRatio;
00447
00448
00449 }
00451
00452 void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
00453 const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
00454 { DefineCurrentMaterial(aCouple);
00455 DefineCurrentParticle(aPartDef);
00456 e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
00457 G4bool b;
00458 sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
00459 e_sigma_max/=massRatio;
00460
00461
00462 }
00464
00465 G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
00466 G4double& fwd_TotCS)
00467 { G4double corr_fac = 1.;
00468 if (forward_CS_mode) {
00469 fwd_TotCS=PrefwdCS;
00470 if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
00471 DefineCurrentMaterial(aCouple);
00472 PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
00473 PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
00474 LastEkinForCS = PreStepEkin;
00475 lastPartDefForCS = aPartDef;
00476 if (PrefwdCS >0. && PreadjCS >0.) {
00477 forward_CS_is_used = true;
00478 LastCSCorrectionFactor = PrefwdCS/PreadjCS;
00479 }
00480 else {
00481 forward_CS_is_used = false;
00482 LastCSCorrectionFactor = 1.;
00483
00484 }
00485
00486 }
00487 corr_fac =LastCSCorrectionFactor;
00488
00489
00490
00491 }
00492 else {
00493 forward_CS_is_used = false;
00494 LastCSCorrectionFactor = 1.;
00495 }
00496 fwd_TotCS=PrefwdCS;
00497 fwd_is_used = forward_CS_is_used;
00498 return corr_fac;
00499 }
00501
00502 G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
00503 const G4MaterialCutsCouple* aCouple, G4double step_length)
00504 { G4double corr_fac = 1.;
00505
00506
00507 G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
00508 G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
00509 if (!forward_CS_is_used || pre_adjCS ==0. || after_fwdCS==0.) {
00510 forward_CS_is_used=false;
00511 G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
00512 corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
00513 LastCSCorrectionFactor = 1.;
00514 }
00515 else {
00516 LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
00517 }
00518
00519
00520
00521 return corr_fac;
00522 }
00524
00525 G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
00526 {
00527 return 1./LastCSCorrectionFactor;
00528
00529 }
00531
00532 G4double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
00533 G4VEmAdjointModel* aModel,
00534 G4double PrimEnergy,
00535 G4double Tcut,
00536 G4bool IsScatProjToProjCase,
00537 std::vector<G4double>& CS_Vs_Element)
00538 {
00539
00540 G4double EminSec=0;
00541 G4double EmaxSec=0;
00542
00543 if (IsScatProjToProjCase){
00544 EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
00545 EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
00546 }
00547 else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
00548 EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
00549 EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
00550 }
00551 if (EminSec >= EmaxSec) return 0.;
00552
00553
00554 G4bool need_to_compute=false;
00555 if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
00556 lastMaterial =aMaterial;
00557 lastPrimaryEnergy = PrimEnergy;
00558 lastTcut=Tcut;
00559 listOfIndexOfAdjointEMModelInAction.clear();
00560 listOfIsScatProjToProjCase.clear();
00561 lastAdjointCSVsModelsAndElements.clear();
00562 need_to_compute=true;
00563
00564 }
00565 size_t ind=0;
00566 if (!need_to_compute){
00567 need_to_compute=true;
00568 for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){
00569 size_t ind1=listOfIndexOfAdjointEMModelInAction[i];
00570 if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){
00571 need_to_compute=false;
00572 CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind];
00573 }
00574 ind++;
00575 }
00576 }
00577
00578 if (need_to_compute){
00579 size_t ind_model=0;
00580 for (size_t i=0;i<listOfAdjointEMModel.size();i++){
00581 if (aModel == listOfAdjointEMModel[i]){
00582 ind_model=i;
00583 break;
00584 }
00585 }
00586 G4double Tlow=Tcut;
00587 if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit();
00588 listOfIndexOfAdjointEMModelInAction.push_back(ind_model);
00589 listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
00590 CS_Vs_Element.clear();
00591 if (!aModel->GetUseMatrix()){
00592 CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
00593
00594
00595 }
00596 else if (aModel->GetUseMatrixPerElement()){
00597 size_t n_el = aMaterial->GetNumberOfElements();
00598 if (aModel->GetUseOnlyOneMatrixForAllElements()){
00599 G4AdjointCSMatrix* theCSMatrix;
00600 if (IsScatProjToProjCase){
00601 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0];
00602 }
00603 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0];
00604 G4double CS =0.;
00605 if (PrimEnergy > Tlow)
00606 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
00607 G4double factor=0.;
00608 for (size_t i=0;i<n_el;i++){
00609
00610 factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
00611 }
00612 CS *=factor;
00613 CS_Vs_Element.push_back(CS);
00614
00615 }
00616 else {
00617 for (size_t i=0;i<n_el;i++){
00618 size_t ind_el = aMaterial->GetElement(i)->GetIndex();
00619
00620 G4AdjointCSMatrix* theCSMatrix;
00621 if (IsScatProjToProjCase){
00622 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
00623 }
00624 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
00625 G4double CS =0.;
00626 if (PrimEnergy > Tlow)
00627 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
00628
00629 CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
00630 }
00631 }
00632
00633 }
00634 else {
00635 size_t ind_mat = aMaterial->GetIndex();
00636 G4AdjointCSMatrix* theCSMatrix;
00637 if (IsScatProjToProjCase){
00638 theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
00639 }
00640 else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat];
00641 G4double CS =0.;
00642 if (PrimEnergy > Tlow)
00643 CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
00644 CS_Vs_Element.push_back(CS);
00645
00646
00647 }
00648 lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
00649
00650 }
00651
00652
00653 G4double CS=0;
00654 for (size_t i=0;i<CS_Vs_Element.size();i++){
00655 CS+=CS_Vs_Element[i];
00656
00657 }
00658 return CS;
00659 }
00661
00662 G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial,
00663 G4VEmAdjointModel* aModel,
00664 G4double PrimEnergy,
00665 G4double Tcut,
00666 G4bool IsScatProjToProjCase)
00667 { std::vector<G4double> CS_Vs_Element;
00668 G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
00669 G4double rand_var= G4UniformRand();
00670 G4double SumCS=0.;
00671 size_t ind=0;
00672 for (size_t i=0;i<CS_Vs_Element.size();i++){
00673 SumCS+=CS_Vs_Element[i];
00674 if (rand_var<=SumCS/CS){
00675 ind=i;
00676 break;
00677 }
00678 }
00679
00680 return const_cast<G4Element*>(aMaterial->GetElement(ind));
00681
00682
00683
00684 }
00686
00687 G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* aCouple,
00688 G4ParticleDefinition* aPartDef,
00689 G4double Ekin)
00690 {
00691 G4double TotalCS=0.;
00692
00693 DefineCurrentMaterial(aCouple);
00694
00695
00696 std::vector<G4double> CS_Vs_Element;
00697 G4double CS;
00698 for (size_t i=0; i<listOfAdjointEMModel.size();i++){
00699
00700 G4double Tlow=0;
00701 if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
00702 else {
00703 G4ParticleDefinition* theDirSecondPartDef =
00704 GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
00705 size_t idx=56;
00706 if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
00707 else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
00708 else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
00709 if (idx <56) {
00710 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
00711 Tlow =(*aVec)[aCouple->GetIndex()];
00712 }
00713
00714
00715 }
00716 if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
00717 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
00718 CS=ComputeAdjointCS(currentMaterial,
00719 listOfAdjointEMModel[i],
00720 Ekin, Tlow,true,CS_Vs_Element);
00721 TotalCS += CS;
00722 (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,CS);
00723 }
00724 if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
00725 CS = ComputeAdjointCS(currentMaterial,
00726 listOfAdjointEMModel[i],
00727 Ekin, Tlow,false, CS_Vs_Element);
00728 TotalCS += CS;
00729 (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,CS);
00730 }
00731
00732 }
00733 else {
00734 (*listSigmaTableForAdjointModelScatProjToProj[i])[currentMatIndex]->PutValue(eindex,0.);
00735 (*listSigmaTableForAdjointModelProdToProj[i])[currentMatIndex]->PutValue(eindex,0.);
00736
00737 }
00738 }
00739 return TotalCS;
00740
00741
00742 }
00744
00745 std::vector<G4AdjointCSMatrix*>
00746 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
00747 G4int nbin_pro_decade)
00748 {
00749 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
00750 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
00751
00752
00753
00754
00755 G4double EkinMin =aModel->GetLowEnergyLimit();
00756 G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
00757 G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
00758 if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
00759
00760
00761
00762
00763 G4double fE=std::pow(10.,1./nbin_pro_decade);
00764 G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
00765 G4double E1=EkinMin;
00766 while (E1 <EkinMaxForProd){
00767 E1=std::max(EkinMin,E2);
00768 E1=std::min(EkinMaxForProd,E1);
00769 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
00770 if (aMat.size()>=2) {
00771 std::vector< double>* log_ESecVec=aMat[0];
00772 std::vector< double>* log_CSVec=aMat[1];
00773 G4double log_adjointCS=log_CSVec->back();
00774
00775 for (size_t j=0;j<log_CSVec->size();j++) {
00776 if (j==0) (*log_CSVec)[j] = 0.;
00777 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
00778 }
00779 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
00780 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
00781 }
00782 E1=E2;
00783 E2*=fE;
00784 }
00785
00786
00787
00788
00789 E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
00790 E1=EkinMin;
00791 while (E1 <EkinMaxForScat){
00792 E1=std::max(EkinMin,E2);
00793 E1=std::min(EkinMaxForScat,E1);
00794 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
00795 if (aMat.size()>=2) {
00796 std::vector< double>* log_ESecVec=aMat[0];
00797 std::vector< double>* log_CSVec=aMat[1];
00798 G4double log_adjointCS=log_CSVec->back();
00799
00800 for (size_t j=0;j<log_CSVec->size();j++) {
00801 if (j==0) (*log_CSVec)[j] = 0.;
00802 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
00803 }
00804 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
00805 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
00806 }
00807 E1=E2;
00808 E2*=fE;
00809 }
00810
00811
00812 std::vector<G4AdjointCSMatrix*> res;
00813 res.clear();
00814 res.push_back(theCSMatForProdToProjBackwardScattering);
00815 res.push_back(theCSMatForScatProjToProjBackwardScattering);
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830 return res;
00831
00832
00833 }
00835
00836 std::vector<G4AdjointCSMatrix*>
00837 G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel* aModel,
00838 G4Material* aMaterial,
00839 G4int nbin_pro_decade)
00840 {
00841 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
00842 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
00843
00844
00845
00846
00847 G4double EkinMin =aModel->GetLowEnergyLimit();
00848 G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
00849 G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
00850 if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860 G4double fE=std::pow(10.,1./nbin_pro_decade);
00861 G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
00862 G4double E1=EkinMin;
00863 while (E1 <EkinMaxForProd){
00864 E1=std::max(EkinMin,E2);
00865 E1=std::min(EkinMaxForProd,E1);
00866 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
00867 if (aMat.size()>=2) {
00868 std::vector< double>* log_ESecVec=aMat[0];
00869 std::vector< double>* log_CSVec=aMat[1];
00870 G4double log_adjointCS=log_CSVec->back();
00871
00872
00873 for (size_t j=0;j<log_CSVec->size();j++) {
00874
00875 if (j==0) (*log_CSVec)[j] = 0.;
00876 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
00877
00878 }
00879 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
00880 theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
00881 }
00882
00883
00884
00885 E1=E2;
00886 E2*=fE;
00887 }
00888
00889
00890
00891
00892 E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
00893 E1=EkinMin;
00894 while (E1 <EkinMaxForScat){
00895 E1=std::max(EkinMin,E2);
00896 E1=std::min(EkinMaxForScat,E1);
00897 std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
00898 if (aMat.size()>=2) {
00899 std::vector< double>* log_ESecVec=aMat[0];
00900 std::vector< double>* log_CSVec=aMat[1];
00901 G4double log_adjointCS=log_CSVec->back();
00902
00903 for (size_t j=0;j<log_CSVec->size();j++) {
00904
00905 if (j==0) (*log_CSVec)[j] = 0.;
00906 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
00907
00908
00909 }
00910 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
00911
00912 theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
00913 }
00914 E1=E2;
00915 E2*=fE;
00916 }
00917
00918
00919
00920
00921
00922
00923
00924 std::vector<G4AdjointCSMatrix*> res;
00925 res.clear();
00926
00927 res.push_back(theCSMatForProdToProjBackwardScattering);
00928 res.push_back(theCSMatForScatProjToProjBackwardScattering);
00929
00930
00931
00932
00933
00934
00935
00936 return res;
00937
00938
00939 }
00940
00942
00943 G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
00944 {
00945 if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
00946 else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
00947 else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
00948 else if (theFwdPartDef ==theFwdIon) return theAdjIon;
00949
00950 return 0;
00951 }
00953
00954 G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
00955 {
00956 if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
00957 else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
00958 else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
00959 else if (theAdjPartDef == theAdjIon) return theFwdIon;
00960 return 0;
00961 }
00963
00964 void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
00965 {
00966 if(couple != currentCouple) {
00967 currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
00968 currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
00969 currentMatIndex = couple->GetIndex();
00970 lastPartDefForCS =0;
00971 LastEkinForCS =0;
00972 LastCSCorrectionFactor =1.;
00973 }
00974 }
00975
00977
00978 void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
00979 {
00980 if(aPartDef != currentParticleDef) {
00981
00982 currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
00983 massRatio=1;
00984 if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
00985 currentParticleIndex=1000000;
00986 for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
00987 if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
00988 }
00989
00990 }
00991 }
00992
00993
00994
00996
00997 G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
00998 anAdjointCSMatrix,G4double Tcut)
00999 {
01000 std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
01001 if (theLogPrimEnergyVector->size() ==0){
01002 G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
01003 G4cout<<"The s"<<G4endl;
01004 return 0.;
01005
01006 }
01007 G4double log_Tcut = std::log(Tcut);
01008 G4double log_E =std::log(aPrimEnergy);
01009
01010 if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back()) return 0.;
01011
01012
01013
01014 G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
01015
01016 size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
01017 G4double aLogPrimEnergy1,aLogPrimEnergy2;
01018 G4double aLogCS1,aLogCS2;
01019 G4double log01,log02;
01020 std::vector< double>* aLogSecondEnergyVector1 =0;
01021 std::vector< double>* aLogSecondEnergyVector2 =0;
01022 std::vector< double>* aLogProbVector1=0;
01023 std::vector< double>* aLogProbVector2=0;
01024 std::vector< size_t>* aLogProbVectorIndex1=0;
01025 std::vector< size_t>* aLogProbVectorIndex2=0;
01026
01027
01028 anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
01029 anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
01030 if (anAdjointCSMatrix->IsScatProjToProjCase()){
01031 G4double log_minimum_prob1, log_minimum_prob2;
01032 log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
01033 log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
01034 aLogCS1+= log_minimum_prob1;
01035 aLogCS2+= log_minimum_prob2;
01036 }
01037
01038 G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2);
01039 return std::exp(log_adjointCS);
01040
01041
01042 }
