| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // particle_hp -- source file 26 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 080721 To be "ClearHistories" effective, th << 30 // 080721 To be "ClearHistories" effective, the selection scheme of angular distribution is modified by T. Koi
31 // by T. Koi <<
32 // 31 //
33 // P. Arce, Dec-2014 Conversion neutron_hp to 32 // P. Arce, Dec-2014 Conversion neutron_hp to particle_hp
34 // V. Ivanchenko, July-2023 Basic revision of <<
35 // 33 //
36 <<
37 #include "G4ParticleHPContEnergyAngular.hh" 34 #include "G4ParticleHPContEnergyAngular.hh"
38 #include "G4HadronicException.hh" <<
39 #include "G4Neutron.hh" <<
40 <<
41 G4ParticleHPContEnergyAngular::G4ParticleHPCon <<
42 { <<
43 theProjectile = (nullptr == proj) ? G4Neutro <<
44 currentMeanEnergy.Put(-2); <<
45 fCacheAngular.Put(nullptr); <<
46 } <<
47 <<
48 G4ParticleHPContEnergyAngular::~G4ParticleHPCo <<
49 { <<
50 delete[] theAngular; <<
51 delete fCacheAngular.Get(); <<
52 } <<
53 <<
54 void G4ParticleHPContEnergyAngular::Init(std:: <<
55 { <<
56 aDataFile >> theTargetCode >> theAngularRep <<
57 const std::size_t esize = nEnergy > 0 ? nEne <<
58 theAngular = new G4ParticleHPContAngularPar[ <<
59 theManager.Init(aDataFile); <<
60 for (G4int i = 0; i < nEnergy; ++i) { <<
61 theAngular[i].Init(aDataFile, theProjectil <<
62 theAngular[i].SetInterpolation(theInterpol <<
63 #ifndef PHP_AS_HP <<
64 theAngular[i].PrepareTableInterpolation(); <<
65 #endif <<
66 } <<
67 } <<
68 35
69 G4ReactionProduct* G4ParticleHPContEnergyAngul << 36 G4ReactionProduct * G4ParticleHPContEnergyAngular::Sample(G4double anEnergy, G4double massCode, G4double /*mass*/)
70 <<
71 { 37 {
72 G4ReactionProduct* result; << 38 G4ReactionProduct * result;
73 G4int i(0); << 39 G4int i(0);
74 G4int it(0); << 40 G4int it(0);
75 for (i = 0; i < nEnergy; ++i) { << 41 for(i=0;i<nEnergy;i++)
76 it = i; << 42 {
77 #ifdef PHP_AS_HP << 43 it = i;
78 if (theAngular[i].GetEnergy() > anEnergy) << 44 #ifdef PHP_AS_HP
>> 45 if(theAngular[i].GetEnergy()>anEnergy) break;
79 #else 46 #else
80 if (theAngular[i].GetEnergy() >= anEnergy) << 47 if(theAngular[i].GetEnergy()>=anEnergy) break;
81 #endif 48 #endif
82 } << 49 }
83 G4double targetMass = GetTarget()->GetMass() << 50 if( std::getenv("G4PHPTEST") ) G4cout << i << " G4ParticleHPContEnergyAngular dataE " << theAngular[i].GetEnergy() << " > " << anEnergy << " it_theAngular " << it << " interpolation " << theInterpolation << G4endl; //GDEB
84 /* << 51 G4double targetMass = GetTarget()->GetMass();
85 G4cout << "G4ParticleHPContEnergyAngular::Sa << 52 if(it==0)
86 << " M=" << targetMass << G4endl; << 53 {
87 */ << 54 theAngular[0].SetTarget(GetTarget());
88 if (it == 0) { << 55 theAngular[0].SetTargetCode(theTargetCode);
89 theAngular[0].SetTarget(GetTarget()); << 56 theAngular[0].SetPrimary(GetProjectileRP());
90 theAngular[0].SetTargetCode(theTargetCode) << 57 result = theAngular[0].Sample(anEnergy, massCode, targetMass,
91 theAngular[0].SetPrimary(GetProjectileRP() << 58 theAngularRep, theInterpolation);
92 result = theAngular[0].Sample(anEnergy, ma << 59 currentMeanEnergy.Put( theAngular[0].MeanEnergyOfThisInteraction() );
93 currentMeanEnergy.Put(theAngular[0].MeanEn << 60 }
94 } << 61 else
95 else { << 62 {
96 // interpolation through alternating sampl << 63 // interpolation through alternating sampling. This needs improvement @@@
97 // This is the cause of the He3 problem !! << 64 // This is the cause of the He3 problem !!!!!!!!
98 // See to it, if you can improve this. << 65 // See to it, if you can improve this.
99 // G4double random = G4UniformRand(); << 66 //080714 TK commnet Randomizing use angular distribution
100 // G4double deltaE = theAngular[it].GetEne << 67 //080714 TK Always use the upper side distribution. enabling ClearHistories method.
101 // G4double offset = theAngular[it].GetEne << 68 //G4double random = G4UniformRand();
102 // if(random<offset/deltaE) it--; << 69 //G4double deltaE = theAngular[it].GetEnergy()-theAngular[it-1].GetEnergy();
103 //--- create new << 70 //G4double offset = theAngular[it].GetEnergy()-anEnergy;
104 // if (theManager.GetScheme(0) != LINLIN) << 71 //if(random<offset/deltaE) it--;
105 // // asserted in G4ParticleHPContEnerg << 72 //--- create new
>> 73 // if( theManager.GetScheme(0) != LINLIN ) { // asserted in G4ParticleHPContEnergyAngular::init there is only one range
106 #ifdef PHP_AS_HP 74 #ifdef PHP_AS_HP
107 theAngular[it].SetTarget(GetTarget()); << 75 theAngular[it].SetTarget(GetTarget());
108 theAngular[it].SetTargetCode(theTargetCode << 76 theAngular[it].SetTargetCode(theTargetCode);
109 theAngular[it].SetPrimary(GetProjectileRP( << 77 theAngular[it].SetPrimary(GetProjectileRP());
110 result = theAngular[it].Sample(anEnergy, m << 78 result = theAngular[it].Sample(anEnergy, massCode, targetMass,
111 currentMeanEnergy.Put(theAngular[it].MeanE << 79 theAngularRep, theInterpolation);
112 << 80 currentMeanEnergy.Put( theAngular[it].MeanEnergyOfThisInteraction() );
113 #else 81 #else
114 if (fCacheAngular.Get() == nullptr) { << 82 if( std::getenv("G4PHPTEST") ) G4cout << i << " G4ParticleHPContEnergyAngular To BUILDBYINTERPOLATION " << it << " : " << theAngular[it].GetEnergy()<< " , " << theAngular[it].GetNEnergies() << " " << it-1 << " : " << theAngular[it-1].GetEnergy()<< " : " << theAngular[it-1].GetNEnergies() << G4endl; //GDEB
115 auto angpar = new G4ParticleHPContAngula << 83
116 fCacheAngular.Put(angpar); << 84 // G4ParticleHPContAngularPar * fAngular = new G4ParticleHPContAngularPar(theProjectile ); //fix start
117 } << 85 if (fCacheAngular.Get() == NULL) {
118 fCacheAngular.Get()->SetInterpolation(theI << 86 G4ParticleHPContAngularPar* angpar = new G4ParticleHPContAngularPar(theProjectile);
119 fCacheAngular.Get()->BuildByInterpolation( << 87 fCacheAngular.Put(angpar);
120 << 88 }
121 fCacheAngular.Get()->SetTarget(GetTarget() << 89 fCacheAngular.Get()->SetInterpolation(theInterpolation);
122 fCacheAngular.Get()->SetTargetCode(theTarg << 90 fCacheAngular.Get()->BuildByInterpolation( anEnergy, theManager.GetScheme(0), (theAngular[it-1]), (theAngular[it]) );
123 fCacheAngular.Get()->SetPrimary(GetProject << 91 fCacheAngular.Get()->SetTarget(GetTarget());
124 << 92 fCacheAngular.Get()->SetTargetCode(theTargetCode);
125 result = << 93 fCacheAngular.Get()->SetPrimary(GetProjectileRP());
126 fCacheAngular.Get()->Sample(anEnergy, ma << 94 result = fCacheAngular.Get()->Sample(anEnergy, massCode, targetMass,
127 currentMeanEnergy.Put(fCacheAngular.Get()- << 95 theAngularRep, theInterpolation);
>> 96 currentMeanEnergy.Put( fCacheAngular.Get()->MeanEnergyOfThisInteraction() );
>> 97 fCacheAngular.Get()->ClearHistories();
>> 98
>> 99 // delete fAngular; //fix end
>> 100
128 #endif 101 #endif
129 } // end (it != 0) branch << 102 }
130 //G4cout << "+!+ Emin=" << currentMeanEnergy << 103
131 return result; << 104 // G4cout << " 0 0 @@@ G4ParticleHPContEnergyAngular::Sample " << result->GetDefinition()->GetParticleName() << " E= " << result->GetKineticEnergy() << G4endl;//GDEB
>> 105 return result;
132 } 106 }
133 107
134 G4double G4ParticleHPContEnergyAngular::MeanEn << 108 G4double G4ParticleHPContEnergyAngular::
>> 109 MeanEnergyOfThisInteraction()
135 { 110 {
136 //G4cout << "+++ Emin=" << currentMeanEnergy << 111 G4double result(0);
137 G4double result(0); << 112 if(currentMeanEnergy.Get()<-1)
138 if (currentMeanEnergy.Get() < -1.0) { << 113 {
139 throw G4HadronicException(__FILE__, __LINE << 114 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContEnergyAngular: Logical error in Product class");
140 "G4ParticleHPCon << 115 }
141 } << 116 else
142 result = currentMeanEnergy.Get(); << 117 {
143 << 118 result = currentMeanEnergy.Get();
144 currentMeanEnergy.Put(-2.0); << 119 }
145 return result; << 120 currentMeanEnergy.Put( -2 );
>> 121 return result;
146 } 122 }
147 123
>> 124
>> 125
148 void G4ParticleHPContEnergyAngular::ClearHisto 126 void G4ParticleHPContEnergyAngular::ClearHistories()
149 { << 127 {
150 if (theAngular != nullptr) { << 128 if ( theAngular!= NULL )
151 for (G4int i = 0; i < nEnergy; ++i) << 129 {
152 theAngular[i].ClearHistories(); << 130 for ( G4int i = 0 ; i< nEnergy ; i++ )
153 } << 131 theAngular[i].ClearHistories();
154 << 132 }
155 // Added fCacheAngular ClearHistories() - th <<
156 // Maybe theAngular does not even need Clear <<
157 if (fCacheAngular.Get() != nullptr) fCacheAn <<
158 } 133 }
159 134