| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file hadronic/Hadr02/src/HadronicInelasti
27 /// \brief Implementation of the HadronicInela
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........oo
50
51 # define MAX_ENERGY_LAB_GEV 10000000.
52 # define MAX_ENERGY_CMS_GEV \
53 30000. // assuming that the target is <=1
54
55 # define IGNORE_PARTICLE_UNKNOWN_PDGID false
56 # define USE_ENERGY_CORR false
57 # define ENERGY_NON_CONSERVATION_RESAMPLE fal
58 # define ENERGY_NON_CONSERVATION_EMAX_GEV 0.9
59 # define ENERGY_NON_CONSERVATION_FRACTION_MAX
60 # define ENERGY_NON_CONSERVATION_FRACTION_MAX
61 # define ENERGY_NON_CONSERVATION_FRACTION_MAX
62
63 # define SPLIT_MULTI_NEUTRONS_MAXN 10
64 # define PARTICLE_MULTI_NEUTRONS_ERRORCODE -1
65
66 # define ERROR_REPORT_EMAIL "andrii.tykhonov@
67 # define CRMC_CONFIG_FILE_ENV_VARIABLE "CRMC_
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........oo
81
82 HadronicInelasticModelCRMC::HadronicInelasticM
83 : G4HadronicInteraction(modelName), fPrintDe
84 {
85 SetMaxEnergy(MAX_ENERGY_LAB_GEV * GeV);
86
87 // int model = 1; // Epos (use temporary), i
88 // int model = 12; // Dpmjet
89 int seed = 123456789;
90 // int seed = CLHEP::HepRandom::getTheSeed()
91 int produce_tables = 0; // CRMC default, se
92 fTypeOutput = 0; // CRMC default, see CRMCo
93 static std::string crmc_param =
94 GetCrmcParamPath(); //"crmc.param"; // CR
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, se
101 crmc_param.c_str(), ""
102
103 // final state
104 finalState = new G4HadFinalState();
105
106 // geant4 particle helpers (tables)
107 fParticleTable = G4ParticleTable::GetParticl
108 fIonTable = fParticleTable->GetIonTable();
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oo
112
113 std::string HadronicInelasticModelCRMC::GetCrm
114 {
115 std::string crmcParamPath = std::getenv(CRMC
116 if (crmcParamPath == "") {
117 std::ostringstream errorstr;
118 errorstr << "CRMC ERROR: could not find cr
119 << CRMC_CONFIG_FILE_ENV_VARIABLE
120 std::string error(errorstr.str());
121 std::cout << error << std::endl;
122 throw error;
123 }
124 std::cout << "Using CRMC parameter file: " <
125 return crmcParamPath;
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oo
129
130 HadronicInelasticModelCRMC::~HadronicInelastic
131
132 //....oooOO0OOooo........oooOO0OOooo........oo
133
134 G4HadFinalState* HadronicInelasticModelCRMC::A
135
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); // T
144 // p
145
146 //* git input particles parameters
147 int id_proj = aTrack.GetDefinition()->GetPDG
148 int id_targ = targetNucleus.GetZ_asInt() * 1
149 double p_proj = aTrack.Get4Momentum().pz() /
150 double e_proj = aTrack.Get4Momentum().e() /
151 double p_targ = 0.;
152 double e_targ =
153 G4NucleiProperties::GetNuclearMass(targetN
154 / GeV;
155 double e_initial = e_proj + e_targ;
156 // ... bug fix (March 2, 2020 - momentum per
157 double a_proj = (double)(aTrack.GetDefinitio
158 if (a_proj < 1.0)
159 a_proj = 1.0; // explanation: if particle
160 double a_targ = (double)(targetNucleus.GetA_
161
162 //* DEBUG messages
163 if (fPrintDebug) {
164 std::cout << "\n\n\n\n\n\n\n==============
165 std::cout << "Start interaction" << std::e
166 std::cout << "id_proj=" << id_proj << std:
167 std::cout << "id_targ=" << id_targ << std:
168 std::cout << "p_proj=" << p_proj << std::e
169 std::cout << "p_targ=" << p_targ << std::e
170 }
171
172 // set up input particle type and energy
173 fInterface->crmc_set(1, // fNCollision,
174 p_proj / a_proj, // fC
175 p_targ / a_targ, // fC
176 id_proj, // fCfg.fProj
177 id_targ); // fCfg.fTar
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_CONSERVA
184 double energy_diff_coef = 1.;
185 double forbid_energy_corr = false;
186 while (true) {
187 // run one interaction
188 fInterface->crmc_generate(fTypeOutput, //
189 1, // iColl+1,
190 gCRMC_data.fNPar
191 gCRMC_data.fPart
192 gCRMC_data.fPart
193 gCRMC_data.fPart
194
195 // split Z=0 A>1 "particles" into multiple
196 SplitMultiNeutrons(gCRMC_data);
197
198 // energy check
199 double e_final = 0;
200 for (int i = 0; i < gCRMC_data.fNParticles
201 if (gCRMC_data.fPartStatus[i] != 1) cont
202 G4ParticleDefinition* pdef;
203 int Z_test = (gCRMC_data.fPartId[i] - CR
204 int A_test =
205 (gCRMC_data.fPartId[i] - CRMC_ION_COEF
206 if (fPrintDebug) {
207 std::cout << std::endl;
208 std::cout << "************************
209 << std::endl;
210 std::cout << "PDG test: " << gCRMC_dat
211 std::cout << "fIonTable->GetIon(Z_test
212 << fIonTable->GetIon(Z_test,
213 std::cout << "ParticleTable->FindParti
214 << fParticleTable->FindParti
215 std::cout << "************************
216 << std::endl;
217 }
218
219 // pdef = fParticleTable->FindParticle(g
220 int pdef_errorcode;
221 pdef = GetParticleDefinition(gCRMC_data.
222 if (!pdef && IGNORE_PARTICLE_UNKNOWN_PDG
223 continue;
224 }
225
226 double p2 = std::pow(gCRMC_data.fPartPx[
227 + std::pow(gCRMC_data.fPartP
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_initia
234 if (e_final != 0. && e_initial != 0. && US
235 if (fPrintDebug) {
236 std::cout << "# e_initial = " << e_initi
237 std::cout << "# e_final = " << e_final
238 std::cout << "# energy_diff_coef = " <<
239 }
240
241 // energy conservation check, if yes
242 if (!ENERGY_NON_CONSERVATION_RESAMPLE) {
243 // ===== NOCHECK ========== NOCHECK ====
244 break;
245 // ===== NOCHECK ========== NOCHECK ====
246 }
247 else if (diff < max_energy_diff || diff /
248 // ===== OK ========== OK ==============
249 forbid_energy_corr = true;
250 break; // everything is fine, no need t
251 // ===== OK ========== OK ==============
252 }
253 else if (resample_attampts < ENERGY_NON_CO
254 resample_attampts++;
255 std::cout << std::endl;
256 std::cout
257 << "#==== WARNING ==== WARNING ==== WA
258 << std::endl;
259 std::cout
260 << "#
261 << std::endl;
262 std::cout
263 << "# [HadronicInelasticModelCRMC::App
264 << std::endl;
265 std::cout << "# e_initial = " << e_initi
266 std::cout << "# e_final = " << e_final
267 std::cout << "# diff = " << diff <<
268 std::cout << "# Running attempt #" << re
269 std::cout
270 << "#
271 << std::endl;
272 std::cout
273 << "#==== WARNING ==== WARNING ==== WA
274 << std::endl;
275 std::cout << std::endl;
276 }
277 else if (max_energy_diff < ENERGY_NON_CONS
278 std::cout << std::endl;
279 std::cout
280 << "#==== WARNING ==== WARNING ==== WA
281 << std::endl;
282 std::cout << "# reached maximum number o
283 << ENERGY_NON_CONSERVATION_FRA
284 << " ==> increasing twice the
285 max_energy_diff *= 2.;
286 std::cout << "# max_energy_diff = " << m
287 std::cout
288 << "#==== WARNING ==== WARNING ==== WA
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 ==== WA
297 << std::endl;
298 std::cout << "# reached maximum number o
299 << ENERGY_NON_CONSERVATION_FRA
300 << std::endl;
301 std::cout
302 << "#==== WARNING ==== WARNING ==== WA
303 << std::endl;
304 std::cout << std::endl;
305 // ===== FAIL ========== FAIL ==========
306 break;
307 // ===== 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;
322 //* Keep only final state particles
323 // .. (-9 is the beam, 2 is a particle whi
324 if (gCRMC_data.fPartStatus[i] != 1) contin
325
326 if (fPrintDebug) {
327 std::cout << "\n\nSecondary:" << std::en
328 << gCRMC_data.fPartId[i] << st
329 << gCRMC_data.fPartPx[i] << st
330 << gCRMC_data.fPartPy[i] << st
331 << gCRMC_data.fPartPz[i] << st
332 << gCRMC_data.fPartEnergy[i] <
333 }
334
335 // G4ParticleDefinition* pdef = fParticleT
336 int pdef_errorcode;
337 G4ParticleDefinition* pdef = GetParticleDe
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] Ge
350 << gCRMC_data.fPartId[i] << std::end
351 std::cout << " [HadronicInelasticModel
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] Ge
374 << gCRMC_data.fPartId[i] << std::end
375 std::cout << " [HadronicInelasticModel
376 << ERROR_REPORT_EMAIL << std
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
391 part_e_corr = 1. / energy_diff_coef;
392 double pbefore2 = std::pow(gCRMC_data.fP
393 + std::pow(gCRMC_data.
394 double mass2 =
395 std::pow(pdef->GetPDGMass() / GeV, 2);
396 double ebefore2 = pbefore2 + mass2;
397 double pafter2 = ebefore2 * part_e_corr
398 if (pbefore2) part_p_corr = std::sqrt(st
399 if (fPrintDebug) std::cout << "part_p_co
400 eaftertest0 += std::sqrt(mass2 + pbefore
401 eaftertest1 += std::sqrt(mass2 + pafter2
402 }
403
404 G4DynamicParticle* part =
405 new G4DynamicParticle(pdef, G4ThreeVecto
406
407
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) = " << tot
416 std::cout << "totalz (GeV) = " << tot
417 std::cout << "initialz (GeV) = " << p_p
418 std::cout << "eaftertest0 = " << eaf
419 std::cout << "eaftertest1 = " << eaf
420 std::cout << "eaftertest2 = " << eaf
421 std::cout << "Finishing interaction: " <<
422 const G4LorentzVector& p1 = aTrack.Get4Mom
423 std::cout << "e=" << p1.e() << " px=" << p
424 << std::endl;
425 std::cout << aTrack.GetDefinition()->GetAt
426 std::cout << aTrack.GetDefinition()->GetPD
427 std::cout << targetNucleus.GetA_asInt() <<
428 std::cout << targetNucleus.GetZ_asInt() <<
429 std::cout << "Stop interaction" << std::en
430 std::cout << "============================
431 }
432 // std::cout<<"finalState->GetNumberOfSecond
433 // std::endl; // Debugging info
434 return finalState;
435 }
436
437 //....oooOO0OOooo........oooOO0OOooo........oo
438
439 G4bool HadronicInelasticModelCRMC::IsApplicabl
440 {
441 return true;
442 }
443
444 //....oooOO0OOooo........oooOO0OOooo........oo
445
446 G4ParticleDefinition* HadronicInelasticModelCR
447
448 {
449 G4ParticleDefinition* pdef = fParticleTable-
450 if (!pdef && particle_id > CRMC_ION_COEF_0)
451 int Z = (particle_id - CRMC_ION_COEF_0) /
452 int A = (particle_id - CRMC_ION_COEF_0 - C
453 if (IsMultiNeutron(Z, A)) {
454 error_code = PARTICLE_MULTI_NEUTRONS_ERR
455 pdef = NULL;
456 }
457 else {
458 pdef = fIonTable->GetIon(Z, A);
459 }
460 }
461 return pdef;
462 }
463
464 //....oooOO0OOooo........oooOO0OOooo........oo
465
466 bool HadronicInelasticModelCRMC::IsMultiNeutro
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 << " [HadronicInelasticModelCR
475 << " is higher than " << SPLIT
476 << std::endl;
477 throw;
478 }
479 }
480 return result;
481 }
482
483 //....oooOO0OOooo........oooOO0OOooo........oo
484
485 void HadronicInelasticModelCRMC::SplitMultiNeu
486 {
487 for (int i = 0; i < CRMC_data.fNParticles; i
488 // check if it is a final-state secondary
489 if (CRMC_data.fPartStatus[i] != 1) continu
490
491 int pdef_errorcode;
492 GetParticleDefinition(CRMC_data.fPartId[i]
493 if (pdef_errorcode != PARTICLE_MULTI_NEUTR
494
495 //
496 int particle_id = gCRMC_data.fPartId[i];
497 int Z = (particle_id - CRMC_ION_COEF_0) /
498 int A = (particle_id - CRMC_ION_COEF_0 - C
499 if (Z != 0 || A < 2) {
500 std::cout << " [HadronicInelasticModelCR
501 "failed! Throwing exception
502 << std::endl;
503 throw;
504 }
505
506 //
507 std::cout << std::endl;
508 std::cout << " [HadronicInelasticModelCRMC
509 "particle into neutrons: "
510 << std::endl;
511 std::cout << " [HadronicInelasticModelCRMC
512 std::cout << " [HadronicInelasticModelCRMC
513 std::cout << " [HadronicInelasticModelCRMC
514 << CRMC_data.fPartId[i] << std::
515 std::cout << " [HadronicInelasticModelCRMC
516 << CRMC_data.fPartPx[i] << std::
517 std::cout << " [HadronicInelasticModelCRMC
518 << CRMC_data.fPartPy[i] << std::
519 std::cout << " [HadronicInelasticModelCRMC
520 << CRMC_data.fPartPz[i] << std::
521 std::cout << " [HadronicInelasticModelCRMC
522 << CRMC_data.fPartEnergy[i] << s
523 std::cout << " [HadronicInelasticModelCRMC
524 << CRMC_data.fPartMass[i] << std
525 std::cout << " [HadronicInelasticModelCRMC
526 << CRMC_data.fPartStatus[i] << s
527
528 //
529 int NEUTRON_PDG_ID = 2112;
530 G4ParticleDefinition* p_n_def = fParticleT
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 << " [HadronicInelasticModelCR
536 << e_n << " lower than neutron
537 << std::endl;
538 e_n = m_n;
539 }
540 double p_tot_before = std::sqrt(CRMC_data.
541 + CRMC_dat
542 + CRMC_dat
543 double p_tot_after = std::sqrt(e_n * e_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_afte
549 py_n = CRMC_data.fPartPy[i] * p_tot_afte
550 pz_n = CRMC_data.fPartPz[i] * p_tot_afte
551 }
552 for (int j = 0; j < A; j++) {
553 int i_neutron = j ? CRMC_data.fNParticle
554 CRMC_data.fPartId[i_neutron] = NEUTRON_P
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] = statu
561 }
562 CRMC_data.fNParticles += A - 1;
563
564 //
565 std::cout << " [HadronicInelasticModelCRMC
566 << std::endl;
567 std::cout << std::endl;
568 }
569 }
570
571 #endif // G4_USE_CRMC
572