ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointCSManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AdjointCSManager.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include <fstream>
29 #include <iomanip>
30 
31 #include "G4AdjointCSManager.hh"
32 
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4AdjointCSMatrix.hh"
36 #include "G4AdjointInterpolator.hh"
37 #include "G4AdjointCSMatrix.hh"
38 #include "G4VEmAdjointModel.hh"
39 #include "G4ElementTable.hh"
40 #include "G4Element.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4Element.hh"
43 #include "G4VEmProcess.hh"
44 #include "G4VEnergyLossProcess.hh"
45 #include "G4PhysicsTable.hh"
46 #include "G4PhysicsLogVector.hh"
47 #include "G4PhysicsTableHelper.hh"
48 #include "G4Electron.hh"
49 #include "G4Gamma.hh"
50 #include "G4Proton.hh"
51 #include "G4AdjointElectron.hh"
52 #include "G4AdjointGamma.hh"
53 #include "G4AdjointProton.hh"
54 #include "G4ProductionCutsTable.hh"
55 #include "G4ProductionCutsTable.hh"
56 
59 //
61 {
62  if(theInstance == nullptr) {
64  theInstance = inst.Instance();
65  }
66  return theInstance;
67 }
68 
70 //
76  listOfForwardEmProcess.clear();
79  EminForFwdSigmaTables.clear();
80  EminForAdjSigmaTables.clear();
81  EkinofFwdSigmaMax.clear();
82  EkinofAdjSigmaMax.clear();
85  Tmin=0.1*keV;
86  Tmax=100.*TeV;
87  nbins=320; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
88 
92 
93  verbose = 1;
95  currentMatIndex = 0;
96  eindex = 0;
97 
98  lastPartDefForCS = nullptr;
101 
102  forward_CS_is_used = true;
103  forward_CS_mode = true;
104 
105  currentParticleDef = nullptr;
106  currentCouple =nullptr;
107  currentMaterial=nullptr;
108  lastMaterial=nullptr;
109 
110  theAdjIon = nullptr;
111  theFwdIon = nullptr;
112 
113  PreadjCS = PostadjCS = PrefwdCS = PostfwdCS = 0.0;
114 }
116 //
118 {;
119 }
121 //
123 {listOfAdjointEMModel.push_back(aModel);
126  return listOfAdjointEMModel.size() -1;
127 
128 }
130 //
132 {
133  G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
134  if (anAdjPartDef && aProcess){
135  RegisterAdjointParticle(anAdjPartDef);
136  G4int index=-1;
137 
138  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
139  if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
140  }
141  listOfForwardEmProcess[index]->push_back(aProcess);
142  }
143 }
145 //
147 {
148  G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
149  if (anAdjPartDef && aProcess){
150  RegisterAdjointParticle(anAdjPartDef);
151  G4int index=-1;
152  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
153  if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
154  }
155  listOfForwardEnergyLossProcess[index]->push_back(aProcess);
156  }
157 }
159 //
161 { G4int index=-1;
162  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
163  if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
164  }
165 
166  if (index ==-1){
167  listOfForwardEnergyLossProcess.push_back(new std::vector<G4VEnergyLossProcess*>());
170  listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
171  theListOfAdjointParticlesInAction.push_back(aPartDef);
172  EminForFwdSigmaTables.push_back(std::vector<G4double> ());
173  EminForAdjSigmaTables.push_back(std::vector<G4double> ());
174  EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
175  EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
176 
177  }
178 }
180 //
182 {
183  if (CrossSectionMatrixesAreBuilt) return;
184  //Tcut, Tmax
185  //The matrices will be computed probably just once
186  //When Tcut will change some PhysicsTable will be recomputed
187  // for each MaterialCutCouple but not all the matrices
188  //The Tcut defines a lower limit in the energy of the Projectile before the scattering
189  //In the Projectile to Scattered Projectile case we have
190  // E_ScatProj<E_Proj-Tcut
191  //Therefore in the adjoint case we have
192  // Eproj> E_ScatProj+Tcut
193  //This implies that when computing the adjoint CS we should integrate over Epro
194  // from E_ScatProj+Tcut to Emax
195  //In the Projectile to Secondary case Tcut plays a role only in the fact that
196  // Esecond should be greater than Tcut to have the possibility to have any adjoint
197  //process
198  //To avoid to recompute the matrices for all changes of MaterialCutCouple
199  //We propose to compute the matrices only once for the minimum possible Tcut and then
200  //to interpolate the probability for a new Tcut (implemented in G4VAdjointEmModel)
201 
202 
205  const G4ElementTable* theElementTable = G4Element::GetElementTable();
206  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
207 
208  G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
209  for (size_t i=0; i<listOfAdjointEMModel.size();i++){
211  G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
212  if (aModel->GetUseMatrix()){
213  std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
214  std::vector<G4AdjointCSMatrix*>* aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
215  aListOfMat1->clear();
216  aListOfMat2->clear();
217  if (aModel->GetUseMatrixPerElement()){
218  if (aModel->GetUseOnlyOneMatrixForAllElements()){
219  std::vector<G4AdjointCSMatrix*>
220  two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
221  aListOfMat1->push_back(two_matrices[0]);
222  aListOfMat2->push_back(two_matrices[1]);
223  }
224  else {
225  for (size_t j=0; j<theElementTable->size();j++){
226  G4Element* anElement=(*theElementTable)[j];
227  G4int Z = G4lrint(anElement->GetZ());
228  G4int A = G4lrint(anElement->GetN());
229  std::vector<G4AdjointCSMatrix*>
230  two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
231  aListOfMat1->push_back(two_matrices[0]);
232  aListOfMat2->push_back(two_matrices[1]);
233  }
234  }
235  }
236  else { //Per material case
237  for (size_t j=0; j<theMaterialTable->size();j++){
238  G4Material* aMaterial=(*theMaterialTable)[j];
239  std::vector<G4AdjointCSMatrix*>
240  two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
241  aListOfMat1->push_back(two_matrices[0]);
242  aListOfMat2->push_back(two_matrices[1]);
243  }
244 
245  }
246  theAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
247  theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
248  aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
249  }
250  else { G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
251  std::vector<G4AdjointCSMatrix*> two_empty_matrices;
252  theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
253  theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
254 
255  }
256  }
257  G4cout<<" All adjoint cross section matrices are computed!"<<G4endl;
258  G4cout<<"======================================================================"<<G4endl;
259 
261 
262 
263 }
264 
265 
267 //
269 { if (TotalSigmaTableAreBuilt) return;
270 
271 
273 
274 
275  //Prepare the Sigma table for all AdjointEMModel, will be filled later on
276  for (size_t i=0; i<listOfAdjointEMModel.size();i++){
277  listSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy();
278  listSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy();
279  for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
282  }
283  }
284 
285 
286 
287  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
289  DefineCurrentParticle(thePartDef);
290  theTotalForwardSigmaTableVector[i]->clearAndDestroy();
291  theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
292  EminForFwdSigmaTables[i].clear();
293  EminForAdjSigmaTables[i].clear();
294  EkinofFwdSigmaMax[i].clear();
295  EkinofAdjSigmaMax[i].clear();
296  //G4cout<<thePartDef->GetParticleName();
297 
298  for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
299  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
300 
301  /*
302  G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
303  G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
304 
305  std::fstream FileOutputAdjCS(file_name1, std::ios::out);
306  std::fstream FileOutputFwdCS(file_name2, std::ios::out);
307 
308 
309 
310  FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
311  FileOutputAdjCS<<std::setprecision(6);
312  FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
313  FileOutputFwdCS<<std::setprecision(6);
314  */
315 
316 
317  //make first the total fwd CS table for FwdProcess
319  G4bool Emin_found=false;
320  G4double sigma_max =0.;
321  G4double e_sigma_max =0.;
322  for(size_t l=0; l<aVector->GetVectorLength(); l++) {
323  G4double totCS=0.;
324  G4double e=aVector->GetLowEdgeEnergy(l);
325  for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
326  totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
327  }
328  for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
329  if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
330  size_t mat_index = couple->GetIndex();
331  G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
332  G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
333  (*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
334  }
336  totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
337  }
338  aVector->PutValue(l,totCS);
339  if (totCS>sigma_max){
340  sigma_max=totCS;
341  e_sigma_max = e;
342 
343  }
344  //FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
345 
346  if (totCS>0 && !Emin_found) {
347  EminForFwdSigmaTables[i].push_back(e);
348  Emin_found=true;
349  }
350 
351 
352  }
353  //FileOutputFwdCS.close();
354 
355  EkinofFwdSigmaMax[i].push_back(e_sigma_max);
356 
357 
358  if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
359 
360  theTotalForwardSigmaTableVector[i]->push_back(aVector);
361 
362 
363  Emin_found=false;
364  sigma_max=0;
365  e_sigma_max =0.;
366  G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins);
367  for(eindex=0; eindex<aVector->GetVectorLength(); eindex++) {
368  G4double e=aVector->GetLowEdgeEnergy(eindex);
369  G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
370  aVector1->PutValue(eindex,totCS);
371  if (totCS>sigma_max){
372  sigma_max=totCS;
373  e_sigma_max = e;
374 
375  }
376  //FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
377  if (totCS>0 && !Emin_found) {
378  EminForAdjSigmaTables[i].push_back(e);
379  Emin_found=true;
380  }
381 
382  }
383  //FileOutputAdjCS.close();
384  EkinofAdjSigmaMax[i].push_back(e_sigma_max);
385  if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
386 
387  theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
388 
389  }
390  }
392 
393 }
395 //
397  const G4MaterialCutsCouple* aCouple)
398 { DefineCurrentMaterial(aCouple);
399  DefineCurrentParticle(aPartDef);
400  G4bool b;
402 
403 
404 
405 }
407 //
409  const G4MaterialCutsCouple* aCouple)
410 { DefineCurrentMaterial(aCouple);
411  DefineCurrentParticle(aPartDef);
412  G4bool b;
414 
415 
416 }
418 //
419 G4double G4AdjointCSManager::GetAdjointSigma(G4double Ekin_nuc, size_t index_model,G4bool is_scat_proj_to_proj,
420  const G4MaterialCutsCouple* aCouple)
421 { DefineCurrentMaterial(aCouple);
422  G4bool b;
423  if (is_scat_proj_to_proj) return (((*listSigmaTableForAdjointModelScatProjToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
424  else return (((*listSigmaTableForAdjointModelProdToProj[index_model])[currentMatIndex])->GetValue(Ekin_nuc, b));
425 }
427 //
429  const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
430 { DefineCurrentMaterial(aCouple);
431  DefineCurrentParticle(aPartDef);
434 
435 
436 
437 }
439 //
441  const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
442 { DefineCurrentMaterial(aCouple);
443  DefineCurrentParticle(aPartDef);
445  G4bool b;
446  sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
447  e_sigma_max/=massRatio;
448 
449 
450 }
452 //
454  const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
455 { DefineCurrentMaterial(aCouple);
456  DefineCurrentParticle(aPartDef);
458  G4bool b;
459  sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
460  e_sigma_max/=massRatio;
461 
462 
463 }
465 //
467  G4double& fwd_TotCS)
468 { G4double corr_fac = 1.;
469  if (forward_CS_mode && aPartDef ) {
470  fwd_TotCS=PrefwdCS;
471  if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
472  DefineCurrentMaterial(aCouple);
473  PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
474  PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
475  LastEkinForCS = PreStepEkin;
476  lastPartDefForCS = aPartDef;
477  if (PrefwdCS >0. && PreadjCS >0.) {
478  forward_CS_is_used = true;
480  }
481  else {
482  forward_CS_is_used = false;
484 
485  }
486 
487  }
488  corr_fac =LastCSCorrectionFactor;
489 
490 
491 
492  }
493  else {
494  forward_CS_is_used = false;
496  }
497  fwd_TotCS=PrefwdCS;
498  fwd_is_used = forward_CS_is_used;
499  return corr_fac;
500 }
502 //
504  const G4MaterialCutsCouple* aCouple, G4double step_length)
505 { G4double corr_fac = 1.;
506  //return corr_fac;
507  //G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
508  G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
509  G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
510  if (!forward_CS_is_used || pre_adjCS ==0. || after_fwdCS==0.) {
511  forward_CS_is_used=false;
512  G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
513  corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
515  }
516  else {
517  LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
518  }
519 
520 
521 
522  return corr_fac;
523 }
525 //
527 {//return 1.;
528  return 1./LastCSCorrectionFactor;
529 
530 }
532 //
534  G4VEmAdjointModel* aModel,
535  G4double PrimEnergy,
536  G4double Tcut,
537  G4bool IsScatProjToProjCase,
538  std::vector<G4double>& CS_Vs_Element)
539 {
540 
541  G4double EminSec=0;
542  G4double EmaxSec=0;
543 
544  if (IsScatProjToProjCase){
545  EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
546  EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
547  }
548  else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
549  EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
550  EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
551  }
552  if (EminSec >= EmaxSec) return 0.;
553 
554 
555  G4bool need_to_compute=false;
556  if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
557  lastMaterial =aMaterial;
558  lastPrimaryEnergy = PrimEnergy;
559  lastTcut=Tcut;
563  need_to_compute=true;
564 
565  }
566  size_t ind=0;
567  if (!need_to_compute){
568  need_to_compute=true;
569  for (size_t i=0;i<listOfIndexOfAdjointEMModelInAction.size();i++){
571  if (aModel == listOfAdjointEMModel[ind1] && IsScatProjToProjCase == listOfIsScatProjToProjCase[i]){
572  need_to_compute=false;
573  CS_Vs_Element = lastAdjointCSVsModelsAndElements[ind];
574  }
575  ind++;
576  }
577  }
578 
579  if (need_to_compute){
580  size_t ind_model=0;
581  for (size_t i=0;i<listOfAdjointEMModel.size();i++){
582  if (aModel == listOfAdjointEMModel[i]){
583  ind_model=i;
584  break;
585  }
586  }
587  G4double Tlow=Tcut;
588  if (!listOfAdjointEMModel[ind_model]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[ind_model]->GetLowEnergyLimit();
589  listOfIndexOfAdjointEMModelInAction.push_back(ind_model);
590  listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
591  CS_Vs_Element.clear();
592  if (!aModel->GetUseMatrix()){
593  CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
594 
595 
596  }
597  else if (aModel->GetUseMatrixPerElement()){
598  size_t n_el = aMaterial->GetNumberOfElements();
599  if (aModel->GetUseOnlyOneMatrixForAllElements()){
600  G4AdjointCSMatrix* theCSMatrix;
601  if (IsScatProjToProjCase){
602  theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][0];
603  }
604  else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][0];
605  G4double CS =0.;
606  if (PrimEnergy > Tlow)
607  CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
608  G4double factor=0.;
609  for (size_t i=0;i<n_el;i++){ //this could be computed only once
610  //size_t ind_el = aMaterial->GetElement(i)->GetIndex();
611  factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
612  }
613  CS *=factor;
614  CS_Vs_Element.push_back(CS);
615 
616  }
617  else {
618  for (size_t i=0;i<n_el;i++){
619  size_t ind_el = aMaterial->GetElement(i)->GetIndex();
620  //G4cout<<aMaterial->GetName()<<G4endl;
621  G4AdjointCSMatrix* theCSMatrix;
622  if (IsScatProjToProjCase){
623  theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
624  }
625  else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
626  G4double CS =0.;
627  if (PrimEnergy > Tlow)
628  CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
629  //G4cout<<CS<<G4endl;
630  CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
631  }
632  }
633 
634  }
635  else {
636  size_t ind_mat = aMaterial->GetIndex();
637  G4AdjointCSMatrix* theCSMatrix;
638  if (IsScatProjToProjCase){
639  theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
640  }
641  else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_mat];
642  G4double CS =0.;
643  if (PrimEnergy > Tlow)
644  CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
645  CS_Vs_Element.push_back(CS);
646 
647 
648  }
649  lastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
650 
651  }
652 
653 
654  G4double CS=0;
655  for (size_t i=0;i<CS_Vs_Element.size();i++){
656  CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
657 
658  }
659  return CS;
660 }
662 //
664  G4VEmAdjointModel* aModel,
665  G4double PrimEnergy,
666  G4double Tcut,
667  G4bool IsScatProjToProjCase)
668 { std::vector<G4double> CS_Vs_Element;
669  G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
670  G4double rand_var= G4UniformRand();
671  G4double SumCS=0.;
672  size_t ind=0;
673  for (size_t i=0;i<CS_Vs_Element.size();i++){
674  SumCS+=CS_Vs_Element[i];
675  if (rand_var<=SumCS/CS){
676  ind=i;
677  break;
678  }
679  }
680 
681  return const_cast<G4Element*>(aMaterial->GetElement(ind));
682 
683 
684 
685 }
687 //
689  G4ParticleDefinition* aPartDef,
690  G4double Ekin)
691 {
692  G4double TotalCS=0.;
693 
694  DefineCurrentMaterial(aCouple);
695 
696 
697  std::vector<G4double> CS_Vs_Element;
698  G4double CS;
699  for (size_t i=0; i<listOfAdjointEMModel.size();i++){
700 
701  G4double Tlow=0;
702  if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
703  else {
704  G4ParticleDefinition* theDirSecondPartDef =
705  GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
706  size_t idx=56;
707  if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
708  else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
709  else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
710  if (idx <56) {
711  const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
712  Tlow =(*aVec)[aCouple->GetIndex()];
713  }
714 
715 
716  }
717  if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
718  if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
721  Ekin, Tlow,true,CS_Vs_Element);
722  TotalCS += CS;
724  }
725  if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
728  Ekin, Tlow,false, CS_Vs_Element);
729  TotalCS += CS;
731  }
732 
733  }
734  else {
737 
738  }
739  }
740  return TotalCS;
741 
742 
743 }
745 //
746 std::vector<G4AdjointCSMatrix*>
748  G4int nbin_pro_decade)
749 {
750  G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
751  G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
752 
753 
754  //make the vector of primary energy of the adjoint particle, could try to make this just once ?
755 
756  G4double EkinMin =aModel->GetLowEnergyLimit();
757  G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
758  G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
759  if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
760 
761 
762  //Product to projectile backward scattering
763  //-----------------------------------------
764  G4double fE=std::pow(10.,1./nbin_pro_decade);
765  G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
766  G4double E1=EkinMin;
767  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
768  while (E1 <EkinMaxForProd){
769  E1=std::max(EkinMin,E2);
770  E1=std::min(EkinMaxForProd,E1);
771  std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
772  if (aMat.size()>=2) {
773  std::vector< double>* log_ESecVec=aMat[0];
774  std::vector< double>* log_CSVec=aMat[1];
775  G4double log_adjointCS=log_CSVec->back();
776  //normalise CSVec such that it becomes a probability vector
777  for (size_t j=0;j<log_CSVec->size();j++) {
778  if (j==0) (*log_CSVec)[j] = 0.;
779  else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
780  }
781  (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
782  theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
783  }
784  E1=E2;
785  E2*=fE;
786  }
787 
788  //Scattered projectile to projectile backward scattering
789  //-----------------------------------------
790 
791  E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
792  E1=EkinMin;
793  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
794  while (E1 <EkinMaxForScat){
795  E1=std::max(EkinMin,E2);
796  E1=std::min(EkinMaxForScat,E1);
797  std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
798  if (aMat.size()>=2) {
799  std::vector< double>* log_ESecVec=aMat[0];
800  std::vector< double>* log_CSVec=aMat[1];
801  G4double log_adjointCS=log_CSVec->back();
802  //normalise CSVec such that it becomes a probability vector
803  for (size_t j=0;j<log_CSVec->size();j++) {
804  if (j==0) (*log_CSVec)[j] = 0.;
805  else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
806  }
807  (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
808  theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
809  }
810  E1=E2;
811  E2*=fE;
812  }
813 
814 
815  std::vector<G4AdjointCSMatrix*> res;
816  res.clear();
817  res.push_back(theCSMatForProdToProjBackwardScattering);
818  res.push_back(theCSMatForScatProjToProjBackwardScattering);
819 
820 
821 /*
822  G4String file_name;
823  std::stringstream astream;
824  G4String str_Z;
825  astream<<Z;
826  astream>>str_Z;
827  theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt");
828  theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
829 
830 */
831 
832 
833  return res;
834 
835 
836 }
838 //
839 std::vector<G4AdjointCSMatrix*>
841  G4Material* aMaterial,
842  G4int nbin_pro_decade)
843 {
844  G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
845  G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
846 
847 
848  //make the vector of primary energy of the adjoint particle, could try to make this just once ?
849 
850  G4double EkinMin =aModel->GetLowEnergyLimit();
851  G4double EkinMaxForScat =aModel->GetHighEnergyLimit()*0.999;
852  G4double EkinMaxForProd =aModel->GetHighEnergyLimit()*0.999;
853  if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
854 
855 
856 
857 
858 
859 
860 
861  //Product to projectile backward scattering
862  //-----------------------------------------
863  G4double fE=std::pow(10.,1./nbin_pro_decade);
864  G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
865  G4double E1=EkinMin;
866  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
867  while (E1 <EkinMaxForProd){
868  E1=std::max(EkinMin,E2);
869  E1=std::min(EkinMaxForProd,E1);
870  std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
871  if (aMat.size()>=2) {
872  std::vector< double>* log_ESecVec=aMat[0];
873  std::vector< double>* log_CSVec=aMat[1];
874  G4double log_adjointCS=log_CSVec->back();
875 
876  //normalise CSVec such that it becomes a probability vector
877  for (size_t j=0;j<log_CSVec->size();j++) {
878  //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
879  if (j==0) (*log_CSVec)[j] = 0.;
880  else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
881  //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
882  }
883  (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
884  theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
885  }
886 
887 
888 
889  E1=E2;
890  E2*=fE;
891  }
892 
893  //Scattered projectile to projectile backward scattering
894  //-----------------------------------------
895 
896  E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
897  E1=EkinMin;
898  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
899  while (E1 <EkinMaxForScat){
900  E1=std::max(EkinMin,E2);
901  E1=std::min(EkinMaxForScat,E1);
902  std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
903  if (aMat.size()>=2) {
904  std::vector< double>* log_ESecVec=aMat[0];
905  std::vector< double>* log_CSVec=aMat[1];
906  G4double log_adjointCS=log_CSVec->back();
907 
908  for (size_t j=0;j<log_CSVec->size();j++) {
909  //G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
910  if (j==0) (*log_CSVec)[j] = 0.;
911  else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
912  //G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
913 
914  }
915  (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
916 
917  theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
918  }
919  E1=E2;
920  E2*=fE;
921  }
922 
923 
924 
925 
926 
927 
928 
929  std::vector<G4AdjointCSMatrix*> res;
930  res.clear();
931 
932  res.push_back(theCSMatForProdToProjBackwardScattering);
933  res.push_back(theCSMatForScatProjToProjBackwardScattering);
934 
935  /*
936  theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
937  theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
938 */
939 
940 
941  return res;
942 
943 
944 }
945 
947 //
949 {
950  if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
951  else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
952  else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
953  else if (theFwdPartDef ==theFwdIon) return theAdjIon;
954 
955  return 0;
956 }
958 //
960 {
961  if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
962  else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
963  else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
964  else if (theAdjPartDef == theAdjIon) return theFwdIon;
965  return 0;
966 }
968 //
970 {
971  if(couple != currentCouple) {
972  currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
973  currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
974  currentMatIndex = couple->GetIndex();
975  lastPartDefForCS =0;
976  LastEkinForCS =0;
978  }
979 }
980 
982 //
984 {
985  if(aPartDef != currentParticleDef) {
986 
987  currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
988  massRatio=1;
989  if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
990  currentParticleIndex=1000000;
991  for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
993  }
994 
995  }
996 }
997 
998 
999 
1001 //
1003  anAdjointCSMatrix,G4double Tcut)
1004 {
1005  std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
1006  if (theLogPrimEnergyVector->size() ==0){
1007  G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
1008  G4cout<<"The s"<<G4endl;
1009  return 0.;
1010 
1011  }
1012  G4double log_Tcut = std::log(Tcut);
1013  G4double log_E =std::log(aPrimEnergy);
1014 
1015  if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back()) return 0.;
1016 
1017 
1018 
1020 
1021  size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
1022  G4double aLogPrimEnergy1,aLogPrimEnergy2;
1023  G4double aLogCS1,aLogCS2;
1024  G4double log01,log02;
1025  std::vector< double>* aLogSecondEnergyVector1 =0;
1026  std::vector< double>* aLogSecondEnergyVector2 =0;
1027  std::vector< double>* aLogProbVector1=0;
1028  std::vector< double>* aLogProbVector2=0;
1029  std::vector< size_t>* aLogProbVectorIndex1=0;
1030  std::vector< size_t>* aLogProbVectorIndex2=0;
1031 
1032 
1033  anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
1034  anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
1035  if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
1036  G4double log_minimum_prob1, log_minimum_prob2;
1037  log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
1038  log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
1039  aLogCS1+= log_minimum_prob1;
1040  aLogCS2+= log_minimum_prob2;
1041  }
1042 
1043  G4double log_adjointCS = theInterpolator->LinearInterpolation(log_E,aLogPrimEnergy1,aLogPrimEnergy2,aLogCS1,aLogCS2);
1044  return std::exp(log_adjointCS);
1045 
1046 
1047 }