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