68{
71
72 G4int iZ = aTarg.GetZ_asInt();
73 G4int iA = aTarg.GetA_asInt();
75 if (aTarg.GetIsotope() != nullptr) iM = aTarg.GetIsotope()->Getm();
76
78
81
84 if (aGIDITarget == nullptr) {
85
91 }
92
93
94
97 std::vector<G4GIDI_Product>* products;
98 do {
100 loop++;
101 } while (products == nullptr && loop < loopMax);
102
103
104
105
106
107
108 if (loop > loopMax - 1) {
109
110
116
117
118
119
120
121
122
123
124 } else {
127
128
133 std::vector<G4int> lightProductIndex;
134 std::vector<G4int> heavyProductIndex;
135 for (
G4int i = 0; i < int( products->size() ); i++ ) {
136 productA = (*products)[i].A;
137 if (productA < 5) {
138 lightProductIndex.push_back(i);
139 GZtot += (*products)[i].Z;
140 GAtot += productA;
141 } else {
142 heavyProductIndex.push_back(i);
143 }
144 }
145
146
147
148
149
150
151
152 auto rng = std::default_random_engine {};
153 std::shuffle(heavyProductIndex.begin(), heavyProductIndex.end(), rng);
154
155
156
157
158
159 std::vector<G4int> savedHeavyIndex;
161 for (
G4int i = 0; i < int(heavyProductIndex.size() ); i++) {
162 itest = heavyProductIndex[i];
163 productA = (*products)[itest].A;
164 productZ = (*products)[itest].Z;
165 if ((GAtot + productA <= iTotA) && (GZtot + productZ <= iTotZ) ) {
166 savedHeavyIndex.push_back(itest);
167 GZtot += productZ;
168 GAtot += productA;
169 }
170 }
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
193 for (
G4int i = 0; i < int(lightProductIndex.size() ); i++) {
194 itest = lightProductIndex[i];
195 productZ = (*products)[itest].
Z;
196 productA = (*products)[itest].A;
198 if (productA == 1 && productZ == 0) {
200 } else if (productA == 1 && productZ == 1) {
202 } else if (productA == 2 && productZ == 1) {
204 } else if (productA == 3 && productZ == 1) {
206 } else if (productA == 4 && productZ == 2) {
208 } else {
210 }
211
213 (*products)[itest].py*MeV,
214 (*products)[itest].pz*MeV );
215 Psum += momentum;
217
219 }
220
222 for (
G4int i = 0; i < int(savedHeavyIndex.size() ); i++) {
223 itest = savedHeavyIndex[i];
224 productZ = (*products)[itest].Z;
225 productA = (*products)[itest].A;
226 productM = (*products)[itest].m;
229 productA,
230 productM) );
232 (*products)[itest].py*MeV,
233 (*products)[itest].pz*MeV );
234 Psum += momentum;
236
238 }
239
240
241
242
243
244
245 if (iTotA - GAtot > 1) {
247 if (iTotZ == GZtot) {
248
249
250
252 } else {
254 }
256
258 }
259 }
260
261 delete products;
262
264
266}
static G4Deuteron * Deuteron()
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetMomentum(const G4ThreeVector &momentum)
std::vector< G4GIDI_Product > * getOthersFinalState(double e_in, double temperature, double(*rng)(void *), void *rngState)
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void SetEnergyChange(G4double anEnergy)
void SetMomentumChange(const G4ThreeVector &aV)
const G4Material * GetMaterial() const
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4HadFinalState theParticleChange
static G4IonTable * GetIonTable()
G4int GetNucleusEncoding(G4int iZ, G4int iA, G4int iM)
G4LENDManager * lend_manager
G4GIDI_target * get_target_from_map(G4int nuclear_code)
G4double GetTemperature() const
static G4Neutron * Neutron()
G4int GetAtomicNumber() const
G4int GetAtomicMass() const
static G4Proton * Proton()
static G4Triton * Triton()