| Geant4 Cross Reference |
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 /// \file hadronic/Hadr02/src/HadronicInelasticModelCRMC.cc
27 /// \brief Implementation of the HadronicInelasticModelCRMC class methods
28 //
29 //
30 //---------------------------------------------------------------------------
31 //
32 #ifdef G4_USE_CRMC
33
34 # include "HadronicInelasticModelCRMC.hh"
35
36 # include "G4IonTable.hh"
37 # include "G4NucleiProperties.hh"
38 # include "G4ParticleDefinition.hh"
39 # include "G4ParticleTable.hh"
40 # include "G4SystemOfUnits.hh"
41 # include "G4ThreeVector.hh"
42 # include "Randomize.hh"
43
44 # include <cmath>
45 # include <cstdlib>
46 # include <iostream>
47 # include <string>
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 # define MAX_ENERGY_LAB_GEV 10000000.
52 # define MAX_ENERGY_CMS_GEV \
53 30000. // assuming that the target is <=100 times heavier than the projectile
54
55 # define IGNORE_PARTICLE_UNKNOWN_PDGID false
56 # define USE_ENERGY_CORR false
57 # define ENERGY_NON_CONSERVATION_RESAMPLE false
58 # define ENERGY_NON_CONSERVATION_EMAX_GEV 0.999
59 # define ENERGY_NON_CONSERVATION_FRACTION_MAX 0.00001
60 # define ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT 10
61 # define ENERGY_NON_CONSERVATION_FRACTION_MAX_ENERGYTRY ENERGY_NON_CONSERVATION_EMAX_GEV
62
63 # define SPLIT_MULTI_NEUTRONS_MAXN 10
64 # define PARTICLE_MULTI_NEUTRONS_ERRORCODE -1
65
66 # define ERROR_REPORT_EMAIL "andrii.tykhonov@SPAMNOTcern.ch"
67 # define CRMC_CONFIG_FILE_ENV_VARIABLE "CRMC_CONFIG_FILE"
68
69 //***********************************
70 // CRMC ION DEFINITION
71 // ID = CRMC_ION_COEF_0 +
72 // CRMC_ION_COEF_Z * Z +
73 // CRMC_ION_COEF_A * A
74 //
75 # define CRMC_ION_COEF_0 1000000000
76 # define CRMC_ION_COEF_Z 10000
77 # define CRMC_ION_COEF_A 10
78 //***********************************
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 HadronicInelasticModelCRMC::HadronicInelasticModelCRMC(int model, const G4String& modelName)
83 : G4HadronicInteraction(modelName), fPrintDebug(false)
84 {
85 SetMaxEnergy(MAX_ENERGY_LAB_GEV * GeV);
86
87 // int model = 1; // Epos (use temporary), it is faster
88 // int model = 12; // Dpmjet
89 int seed = 123456789;
90 // int seed = CLHEP::HepRandom::getTheSeed(); // Returns 0 which is invalid
91 int produce_tables = 0; // CRMC default, see CRMCoptions.cc in the CRMC package
92 fTypeOutput = 0; // CRMC default, see CRMCoptions.cc in the CRMC package
93 static std::string crmc_param =
94 GetCrmcParamPath(); //"crmc.param"; // CRMC default, see CRMCoptions.cc in the CRMC package
95
96 fInterface = new CRMCinterface();
97 fInterface->init(model);
98
99 // open FORTRAN IO at first call
100 fInterface->crmc_init(MAX_ENERGY_CMS_GEV, seed, model, produce_tables, fTypeOutput,
101 crmc_param.c_str(), "", 0);
102
103 // final state
104 finalState = new G4HadFinalState();
105
106 // geant4 particle helpers (tables)
107 fParticleTable = G4ParticleTable::GetParticleTable();
108 fIonTable = fParticleTable->GetIonTable();
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112
113 std::string HadronicInelasticModelCRMC::GetCrmcParamPath()
114 {
115 std::string crmcParamPath = std::getenv(CRMC_CONFIG_FILE_ENV_VARIABLE);
116 if (crmcParamPath == "") {
117 std::ostringstream errorstr;
118 errorstr << "CRMC ERROR: could not find crmc param file, please check "
119 << CRMC_CONFIG_FILE_ENV_VARIABLE << " envornoment variable!";
120 std::string error(errorstr.str());
121 std::cout << error << std::endl;
122 throw error;
123 }
124 std::cout << "Using CRMC parameter file: " << crmcParamPath << std::endl;
125 return crmcParamPath;
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129
130 HadronicInelasticModelCRMC::~HadronicInelasticModelCRMC() {}
131
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133
134 G4HadFinalState* HadronicInelasticModelCRMC::ApplyYourself(const G4HadProjectile& aTrack,
135 G4Nucleus& targetNucleus)
136 {
137 //* leanup data vectors
138 gCRMC_data.Clean();
139
140 //* cleanup geant4 final state vector
141 finalState->Clear();
142 finalState->SetStatusChange(
143 G4HadFinalStateStatus::stopAndKill); // TODO: check: inelastic collisions kills previos
144 // particles?
145
146 //* git input particles parameters
147 int id_proj = aTrack.GetDefinition()->GetPDGEncoding();
148 int id_targ = targetNucleus.GetZ_asInt() * 10000 + targetNucleus.GetA_asInt() * 10;
149 double p_proj = aTrack.Get4Momentum().pz() / GeV;
150 double e_proj = aTrack.Get4Momentum().e() / GeV;
151 double p_targ = 0.;
152 double e_targ =
153 G4NucleiProperties::GetNuclearMass(targetNucleus.GetA_asInt(), targetNucleus.GetZ_asInt())
154 / GeV;
155 double e_initial = e_proj + e_targ;
156 // ... bug fix (March 2, 2020 - momentum per nucleon!)
157 double a_proj = (double)(aTrack.GetDefinition()->GetAtomicMass()); // GetAtomicNumber());
158 if (a_proj < 1.0)
159 a_proj = 1.0; // explanation: if particle is not an ion/proton, the GetAtomicMass returns 0
160 double a_targ = (double)(targetNucleus.GetA_asInt());
161
162 //* DEBUG messages
163 if (fPrintDebug) {
164 std::cout << "\n\n\n\n\n\n\n==============================================" << std::endl;
165 std::cout << "Start interaction" << std::endl;
166 std::cout << "id_proj=" << id_proj << std::endl;
167 std::cout << "id_targ=" << id_targ << std::endl;
168 std::cout << "p_proj=" << p_proj << std::endl;
169 std::cout << "p_targ=" << p_targ << std::endl;
170 }
171
172 // set up input particle type and energy
173 fInterface->crmc_set(1, // fNCollision,
174 p_proj / a_proj, // fCfg.fProjectileMomentum (per nucleon!!!),
175 p_targ / a_targ, // fCfg.fTargetMomentum (per nucleon!!!),
176 id_proj, // fCfg.fProjectileId,
177 id_targ); // fCfg.fTargetId);
178
179 //=================================================
180 // sample 1 interaction until the energy
181 // conservation is fulfilled
182 int resample_attampts = 1;
183 double max_energy_diff = ENERGY_NON_CONSERVATION_EMAX_GEV;
184 double energy_diff_coef = 1.;
185 double forbid_energy_corr = false;
186 while (true) {
187 // run one interaction
188 fInterface->crmc_generate(fTypeOutput, // fCfg.fTypoaut,
189 1, // iColl+1,
190 gCRMC_data.fNParticles, gCRMC_data.fImpactParameter,
191 gCRMC_data.fPartId[0], gCRMC_data.fPartPx[0], gCRMC_data.fPartPy[0],
192 gCRMC_data.fPartPz[0], gCRMC_data.fPartEnergy[0],
193 gCRMC_data.fPartMass[0], gCRMC_data.fPartStatus[0]);
194
195 // split Z=0 A>1 "particles" into multiple neutrons
196 SplitMultiNeutrons(gCRMC_data);
197
198 // energy check
199 double e_final = 0;
200 for (int i = 0; i < gCRMC_data.fNParticles; i++) {
201 if (gCRMC_data.fPartStatus[i] != 1) continue; // only final state particles
202 G4ParticleDefinition* pdef;
203 int Z_test = (gCRMC_data.fPartId[i] - CRMC_ION_COEF_0) / CRMC_ION_COEF_Z;
204 int A_test =
205 (gCRMC_data.fPartId[i] - CRMC_ION_COEF_0 - CRMC_ION_COEF_Z * Z_test) / CRMC_ION_COEF_A;
206 if (fPrintDebug) {
207 std::cout << std::endl;
208 std::cout << "**********************************************************************"
209 << std::endl;
210 std::cout << "PDG test: " << gCRMC_data.fPartId[i] << std::endl;
211 std::cout << "fIonTable->GetIon(Z_test, A_test) = "
212 << fIonTable->GetIon(Z_test, A_test) << std::endl;
213 std::cout << "ParticleTable->FindParticle(gCRMC_data.fPartId[i]) = "
214 << fParticleTable->FindParticle(gCRMC_data.fPartId[i]) << std::endl;
215 std::cout << "**********************************************************************"
216 << std::endl;
217 }
218
219 // pdef = fParticleTable->FindParticle(gCRMC_data.fPartId[i]);
220 int pdef_errorcode;
221 pdef = GetParticleDefinition(gCRMC_data.fPartId[i], pdef_errorcode);
222 if (!pdef && IGNORE_PARTICLE_UNKNOWN_PDGID) {
223 continue;
224 }
225
226 double p2 = std::pow(gCRMC_data.fPartPx[i], 2) + std::pow(gCRMC_data.fPartPy[i], 2)
227 + std::pow(gCRMC_data.fPartPz[i], 2);
228 double mass = pdef->GetPDGMass() / GeV;
229 e_final += std::sqrt(mass * mass + p2);
230 }
231
232 // Check if we need to resample again...
233 double diff = std::fabs(e_final - e_initial);
234 if (e_final != 0. && e_initial != 0. && USE_ENERGY_CORR) energy_diff_coef = e_final / e_initial;
235 if (fPrintDebug) {
236 std::cout << "# e_initial = " << e_initial << " GeV" << std::endl;
237 std::cout << "# e_final = " << e_final << " GeV" << std::endl;
238 std::cout << "# energy_diff_coef = " << energy_diff_coef << std::endl;
239 }
240
241 // energy conservation check, if yes
242 if (!ENERGY_NON_CONSERVATION_RESAMPLE) {
243 // ===== NOCHECK ========== NOCHECK ============== NOCHECK ========
244 break;
245 // ===== NOCHECK ========== NOCHECK ============== NOCHECK ========
246 }
247 else if (diff < max_energy_diff || diff / e_initial < ENERGY_NON_CONSERVATION_FRACTION_MAX) {
248 // ===== OK ========== OK ============== OK ========
249 forbid_energy_corr = true;
250 break; // everything is fine, no need to resample, break the re-sampling loop
251 // ===== OK ========== OK ============== OK ========
252 }
253 else if (resample_attampts < ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT) {
254 resample_attampts++;
255 std::cout << std::endl;
256 std::cout
257 << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
258 << std::endl;
259 std::cout
260 << "# #"
261 << std::endl;
262 std::cout
263 << "# [HadronicInelasticModelCRMC::ApplyYourself]: Energy non conservation detected: #"
264 << std::endl;
265 std::cout << "# e_initial = " << e_initial << " GeV" << std::endl;
266 std::cout << "# e_final = " << e_final << " GeV" << std::endl;
267 std::cout << "# diff = " << diff << " GeV" << std::endl;
268 std::cout << "# Running attempt #" << resample_attampts << std::endl;
269 std::cout
270 << "# #"
271 << std::endl;
272 std::cout
273 << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
274 << std::endl;
275 std::cout << std::endl;
276 }
277 else if (max_energy_diff < ENERGY_NON_CONSERVATION_FRACTION_MAX_ENERGYTRY) {
278 std::cout << std::endl;
279 std::cout
280 << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
281 << std::endl;
282 std::cout << "# reached maximum number of attempts = "
283 << ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT
284 << " ==> increasing twice the energy threshold!" << std::endl;
285 max_energy_diff *= 2.;
286 std::cout << "# max_energy_diff = " << max_energy_diff << std::endl;
287 std::cout
288 << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
289 << std::endl;
290 std::cout << std::endl;
291 resample_attampts = 1;
292 }
293 else {
294 std::cout << std::endl;
295 std::cout
296 << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
297 << std::endl;
298 std::cout << "# reached maximum number of attempts = "
299 << ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT << "not resampling any more!"
300 << std::endl;
301 std::cout
302 << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
303 << std::endl;
304 std::cout << std::endl;
305 // ===== FAIL ========== FAIL ============== FAIl ========
306 break;
307 // ===== FAIL ========== FAIL ============== FAIl ========
308 }
309 }
310 // ... finished sampling one interaction
311 //=================================================
312
313 // ... for DEBUG messages
314 double totalenergy = 0;
315 double totalz = 0;
316 double eaftertest0 = 0.;
317 double eaftertest1 = 0.;
318 double eaftertest2 = 0.;
319
320 // save secondary particles for outputa
321 for (int i = 0; i < gCRMC_data.fNParticles; i++) {
322 //* Keep only final state particles
323 // .. (-9 is the beam, 2 is a particle which decayed and 1 is final)
324 if (gCRMC_data.fPartStatus[i] != 1) continue;
325
326 if (fPrintDebug) {
327 std::cout << "\n\nSecondary:" << std::endl
328 << gCRMC_data.fPartId[i] << std::endl
329 << gCRMC_data.fPartPx[i] << std::endl
330 << gCRMC_data.fPartPy[i] << std::endl
331 << gCRMC_data.fPartPz[i] << std::endl
332 << gCRMC_data.fPartEnergy[i] << " ENERGY " << std::endl;
333 }
334
335 // G4ParticleDefinition* pdef = fParticleTable->FindParticle(gCRMC_data.fPartId[i]);
336 int pdef_errorcode;
337 G4ParticleDefinition* pdef = GetParticleDefinition(gCRMC_data.fPartId[i], pdef_errorcode);
338 if (!pdef) {
339 if (IGNORE_PARTICLE_UNKNOWN_PDGID) {
340 std::cout << std::endl;
341 std::cout << "*****************************************************************************"
342 "***************************"
343 << std::endl;
344 std::cout << " -- WARNING "
345 "-----------------------------------------------------------------------------"
346 "------ WARNING -- "
347 << std::endl;
348 std::cout
349 << " [HadronicInelasticModelCRMC] Geant4 could not find particle definition for PDG ID = "
350 << gCRMC_data.fPartId[i] << std::endl;
351 std::cout << " [HadronicInelasticModelCRMC] Ignoring this particle. This might cause "
352 "energy non-conservation!"
353 << std::endl;
354 std::cout << " -- WARNING "
355 "-----------------------------------------------------------------------------"
356 "------ WARNING -- "
357 << std::endl;
358 std::cout << "*****************************************************************************"
359 "***************************"
360 << std::endl;
361 continue;
362 }
363 else {
364 std::cout << std::endl;
365 std::cout << "*****************************************************************************"
366 "***************************"
367 << std::endl;
368 std::cout << " -- ERROR "
369 "-----------------------------------------------------------------------------"
370 "------ ERROR -- "
371 << std::endl;
372 std::cout
373 << " [HadronicInelasticModelCRMC] Geant4 could not find particle definition for PDG ID = "
374 << gCRMC_data.fPartId[i] << std::endl;
375 std::cout << " [HadronicInelasticModelCRMC] Throwing exception! Please report to: "
376 << ERROR_REPORT_EMAIL << std::endl;
377 std::cout << " -- ERROR "
378 "-----------------------------------------------------------------------------"
379 "------ ERROR -- "
380 << std::endl;
381 std::cout << "*****************************************************************************"
382 "***************************"
383 << std::endl;
384 throw;
385 }
386 }
387
388 double part_e_corr = 1.;
389 double part_p_corr = 1.;
390 if (USE_ENERGY_CORR && !forbid_energy_corr && energy_diff_coef != 0) {
391 part_e_corr = 1. / energy_diff_coef;
392 double pbefore2 = std::pow(gCRMC_data.fPartPx[i], 2) + std::pow(gCRMC_data.fPartPy[i], 2)
393 + std::pow(gCRMC_data.fPartPz[i], 2);
394 double mass2 =
395 std::pow(pdef->GetPDGMass() / GeV, 2); // std::pow(gCRMC_data.fPartEnergy[i],2) - pbefore2;
396 double ebefore2 = pbefore2 + mass2;
397 double pafter2 = ebefore2 * part_e_corr * part_e_corr - mass2;
398 if (pbefore2) part_p_corr = std::sqrt(std::fabs(pafter2 / pbefore2));
399 if (fPrintDebug) std::cout << "part_p_corr=" << part_p_corr << std::endl;
400 eaftertest0 += std::sqrt(mass2 + pbefore2);
401 eaftertest1 += std::sqrt(mass2 + pafter2);
402 }
403
404 G4DynamicParticle* part =
405 new G4DynamicParticle(pdef, G4ThreeVector(gCRMC_data.fPartPx[i] * GeV * part_p_corr,
406 gCRMC_data.fPartPy[i] * GeV * part_p_corr,
407 gCRMC_data.fPartPz[i] * GeV * part_p_corr));
408 eaftertest2 += part->GetTotalEnergy();
409 finalState->AddSecondary(part);
410 totalenergy += gCRMC_data.fPartEnergy[i];
411 totalz += gCRMC_data.fPartPz[i];
412 }
413
414 if (fPrintDebug) {
415 std::cout << "totalenergy (GeV) = " << totalenergy << std::endl;
416 std::cout << "totalz (GeV) = " << totalz << std::endl;
417 std::cout << "initialz (GeV) = " << p_proj + p_targ << std::endl;
418 std::cout << "eaftertest0 = " << eaftertest0 << std::endl;
419 std::cout << "eaftertest1 = " << eaftertest1 << std::endl;
420 std::cout << "eaftertest2 = " << eaftertest2 << std::endl;
421 std::cout << "Finishing interaction: " << std::endl;
422 const G4LorentzVector& p1 = aTrack.Get4Momentum();
423 std::cout << "e=" << p1.e() << " px=" << p1.px() << " py=" << p1.py() << " pz=" << p1.pz()
424 << std::endl;
425 std::cout << aTrack.GetDefinition()->GetAtomicNumber() << std::endl;
426 std::cout << aTrack.GetDefinition()->GetPDGCharge() << std::endl;
427 std::cout << targetNucleus.GetA_asInt() << std::endl;
428 std::cout << targetNucleus.GetZ_asInt() << std::endl;
429 std::cout << "Stop interaction" << std::endl;
430 std::cout << "==============================================\n\n\n\n\n\n" << std::endl;
431 }
432 // std::cout<<"finalState->GetNumberOfSecondaries()="<<finalState->GetNumberOfSecondaries()<<
433 // std::endl; // Debugging info
434 return finalState;
435 }
436
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438
439 G4bool HadronicInelasticModelCRMC::IsApplicable(const G4HadProjectile&, G4Nucleus&)
440 {
441 return true;
442 }
443
444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
445
446 G4ParticleDefinition* HadronicInelasticModelCRMC::GetParticleDefinition(long particle_id,
447 int& error_code)
448 {
449 G4ParticleDefinition* pdef = fParticleTable->FindParticle(particle_id);
450 if (!pdef && particle_id > CRMC_ION_COEF_0) {
451 int Z = (particle_id - CRMC_ION_COEF_0) / CRMC_ION_COEF_Z;
452 int A = (particle_id - CRMC_ION_COEF_0 - CRMC_ION_COEF_Z * Z) / CRMC_ION_COEF_A;
453 if (IsMultiNeutron(Z, A)) {
454 error_code = PARTICLE_MULTI_NEUTRONS_ERRORCODE;
455 pdef = NULL;
456 }
457 else {
458 pdef = fIonTable->GetIon(Z, A);
459 }
460 }
461 return pdef;
462 }
463
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
465
466 bool HadronicInelasticModelCRMC::IsMultiNeutron(int Z, int A)
467 {
468 bool result = false;
469 if (!Z && A > 1) {
470 if (A <= SPLIT_MULTI_NEUTRONS_MAXN) {
471 result = true;
472 }
473 else {
474 std::cout << " [HadronicInelasticModelCRMC::IsMultiNeutron] ERROR A=" << A
475 << " is higher than " << SPLIT_MULTI_NEUTRONS_MAXN << " throwing exception!"
476 << std::endl;
477 throw;
478 }
479 }
480 return result;
481 }
482
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
484
485 void HadronicInelasticModelCRMC::SplitMultiNeutrons(CRMCdata& CRMC_data)
486 {
487 for (int i = 0; i < CRMC_data.fNParticles; i++) {
488 // check if it is a final-state secondary particle
489 if (CRMC_data.fPartStatus[i] != 1) continue;
490
491 int pdef_errorcode;
492 GetParticleDefinition(CRMC_data.fPartId[i], pdef_errorcode);
493 if (pdef_errorcode != PARTICLE_MULTI_NEUTRONS_ERRORCODE) continue;
494
495 //
496 int particle_id = gCRMC_data.fPartId[i];
497 int Z = (particle_id - CRMC_ION_COEF_0) / CRMC_ION_COEF_Z;
498 int A = (particle_id - CRMC_ION_COEF_0 - CRMC_ION_COEF_Z * Z) / CRMC_ION_COEF_A;
499 if (Z != 0 || A < 2) {
500 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] ERROR consistency check "
501 "failed! Throwing exception! "
502 << std::endl;
503 throw;
504 }
505
506 //
507 std::cout << std::endl;
508 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO splitting the floowing "
509 "particle into neutrons: "
510 << std::endl;
511 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO Z = " << Z << std::endl;
512 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO A = " << A << std::endl;
513 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartId = "
514 << CRMC_data.fPartId[i] << std::endl;
515 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartPx = "
516 << CRMC_data.fPartPx[i] << std::endl;
517 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartPy = "
518 << CRMC_data.fPartPy[i] << std::endl;
519 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartPz = "
520 << CRMC_data.fPartPz[i] << std::endl;
521 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartEnergy = "
522 << CRMC_data.fPartEnergy[i] << std::endl;
523 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartMass = "
524 << CRMC_data.fPartMass[i] << std::endl;
525 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO fPartStatus = "
526 << CRMC_data.fPartStatus[i] << std::endl;
527
528 //
529 int NEUTRON_PDG_ID = 2112;
530 G4ParticleDefinition* p_n_def = fParticleTable->FindParticle(NEUTRON_PDG_ID);
531 double m_n = p_n_def->GetPDGMass() / GeV;
532 double e_n = CRMC_data.fPartEnergy[i] / A;
533 int status_n = CRMC_data.fPartStatus[i];
534 if (e_n < m_n) {
535 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] WARNING neutron energy "
536 << e_n << " lower than neutron mass " << m_n << " assigning e_n = m_n "
537 << std::endl;
538 e_n = m_n;
539 }
540 double p_tot_before = std::sqrt(CRMC_data.fPartPx[i] * CRMC_data.fPartPx[i]
541 + CRMC_data.fPartPy[i] * CRMC_data.fPartPy[i]
542 + CRMC_data.fPartPz[i] * CRMC_data.fPartPz[i]);
543 double p_tot_after = std::sqrt(e_n * e_n - m_n * m_n);
544 double px_n = 0;
545 double py_n = 0;
546 double pz_n = 0;
547 if (p_tot_before > 0. && p_tot_after > 0.) {
548 px_n = CRMC_data.fPartPx[i] * p_tot_after / p_tot_before;
549 py_n = CRMC_data.fPartPy[i] * p_tot_after / p_tot_before;
550 pz_n = CRMC_data.fPartPz[i] * p_tot_after / p_tot_before;
551 }
552 for (int j = 0; j < A; j++) {
553 int i_neutron = j ? CRMC_data.fNParticles + j : i;
554 CRMC_data.fPartId[i_neutron] = NEUTRON_PDG_ID;
555 CRMC_data.fPartPx[i_neutron] = px_n;
556 CRMC_data.fPartPy[i_neutron] = py_n;
557 CRMC_data.fPartPz[i_neutron] = pz_n;
558 CRMC_data.fPartEnergy[i_neutron] = e_n;
559 CRMC_data.fPartMass[i_neutron] = m_n;
560 CRMC_data.fPartStatus[i_neutron] = status_n;
561 }
562 CRMC_data.fNParticles += A - 1;
563
564 //
565 std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] done for a particle. "
566 << std::endl;
567 std::cout << std::endl;
568 }
569 }
570
571 #endif // G4_USE_CRMC
572