| 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 // * 20 // * *
21 // * Parts of this code which have been devel 21 // * Parts of this code which have been developed by Abdel-Waged *
22 // * et al under contract (31-465) to the King 22 // * et al under contract (31-465) to the King Abdul-Aziz City for *
23 // * Science and Technology (KACST), the Natio 23 // * Science and Technology (KACST), the National Centre of *
24 // * Mathematics and Physics (NCMP), Saudi Ara 24 // * Mathematics and Physics (NCMP), Saudi Arabia. *
25 // * 25 // * *
26 // * By using, copying, modifying or distri 26 // * By using, copying, modifying or distributing the software (or *
27 // * any work based on the software) you ag 27 // * any work based on the software) you agree to acknowledge its *
28 // * use in resulting scientific publicati 28 // * use in resulting scientific publications, and indicate your *
29 // * acceptance of all terms of the Geant4 Sof 29 // * acceptance of all terms of the Geant4 Software license. *
30 // ******************************************* 30 // ********************************************************************
31 // 31 //
32 /// \file hadronic/Hadr02/src/G4UrQMD1_3Model. 32 /// \file hadronic/Hadr02/src/G4UrQMD1_3Model.cc
33 /// \brief Implementation of the G4UrQMD1_3Mod 33 /// \brief Implementation of the G4UrQMD1_3Model class
34 // 34 //
>> 35 // $Id: G4UrQMD1_3Model.cc 81932 2014-06-06 15:39:45Z gcosmo $
35 // 36 //
36 ////////////////////////////////////////////// 37 ///////////////////////////////////////////////////////////////////////////////
37 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 38 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38 // 39 //
39 // MODULE: G4UrQMD1_3Model.hh 40 // MODULE: G4UrQMD1_3Model.hh
40 // 41 //
41 // Version: 0.B 42 // Version: 0.B
42 // Date: 25/01/12 43 // Date: 25/01/12
43 // Authors: Kh. Abdel-Waged and Nuha Fe 44 // Authors: Kh. Abdel-Waged and Nuha Felemban
44 // Revised by: V.V. Uzhinskii << 45 // Revised by: V.V. Uzhinskii
45 // SPONSERED BY 46 // SPONSERED BY
46 // Customer: KAUST/NCMP 47 // Customer: KAUST/NCMP
47 // Contract: 31-465 48 // Contract: 31-465
48 // 49 //
49 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 50 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50 // << 51 //
51 #ifdef G4_USE_URQMD 52 #ifdef G4_USE_URQMD
52 53
53 # include "G4UrQMD1_3Model.hh" << 54 #include "G4UrQMD1_3Model.hh"
54 << 55 #include "G4UrQMD1_3Interface.hh"
55 # include "G4UrQMD1_3Interface.hh" <<
56 //------------------------------- 56 //-------------------------------
57 # include "G4CollisionOutput.hh" << 57 #include "globals.hh"
58 # include "G4DynamicParticle.hh" << 58 #include "G4DynamicParticle.hh"
59 # include "G4IonTable.hh" << 59 #include "G4IonTable.hh"
60 # include "G4LorentzRotation.hh" << 60 #include "G4CollisionOutput.hh"
61 # include "G4Nucleus.hh" << 61 #include "G4V3DNucleus.hh"
62 # include "G4ParticleDefinition.hh" << 62 #include "G4Track.hh"
63 # include "G4ParticleTable.hh" << 63 #include "G4Nucleus.hh"
64 # include "G4PhysicalConstants.hh" << 64 #include "G4LorentzRotation.hh"
65 # include "G4SystemOfUnits.hh" << 65
66 # include "G4Track.hh" << 66 #include "G4ParticleDefinition.hh"
67 # include "G4V3DNucleus.hh" << 67 #include "G4ParticleTable.hh"
68 # include "globals.hh" << 68 #include "G4PhysicalConstants.hh"
69 << 69 #include "G4SystemOfUnits.hh"
70 // AND-> << 70
71 # include "G4Version.hh" << 71 //AND->
72 // AND<- << 72 #include "G4Version.hh"
>> 73 //AND<-
73 //----------------new_anti 74 //----------------new_anti
74 # include "G4AntiAlpha.hh" << 75 #include "G4AntiHe3.hh"
75 # include "G4AntiDeuteron.hh" << 76 #include "G4AntiDeuteron.hh"
76 # include "G4AntiHe3.hh" << 77 #include "G4AntiTriton.hh"
77 # include "G4AntiTriton.hh" << 78 #include "G4AntiAlpha.hh"
78 //--------------------------- 79 //---------------------------
79 # include <fstream> << 80 #include <fstream>
80 # include <string> << 81 #include <string>
>> 82
>> 83 /////////////////////////////////////////////////////////////////////////////////
81 84
82 //....oooOO0OOooo........oooOO0OOooo........oo <<
83 // 85 //
84 G4UrQMD1_3Model::G4UrQMD1_3Model(const G4Strin 86 G4UrQMD1_3Model::G4UrQMD1_3Model(const G4String& nam)
85 : G4VIntraNuclearTransportModel(nam), verbos << 87 :G4VIntraNuclearTransportModel(nam), verbose(0)
86 { 88 {
>> 89
>> 90
87 if (verbose > 3) { 91 if (verbose > 3) {
88 G4cout << " >>> G4UrQMD1_3Model default co 92 G4cout << " >>> G4UrQMD1_3Model default constructor" << G4endl;
89 } 93 }
90 94
91 // <<
92 // Set the minimum and maximum range for the <<
93 95
94 // SetMinEnergy(100.0*MeV); <<
95 // SetMaxEnergy(200.0*GeV); <<
96 96
97 // << 97 //
>> 98 // Set the minimum and maximum range for the UrQMD model
>> 99
>> 100 // SetMinEnergy(100.0*MeV);
>> 101 // SetMaxEnergy(200.0*GeV);
98 102
99 // << 103 //
>> 104
>> 105 //
100 WelcomeMessage(); 106 WelcomeMessage();
101 // << 107 //
102 CurrentEvent = 0; << 108 CurrentEvent=0;
103 // << 109 //
104 110
105 InitialiseDataTables(); << 111 InitialiseDataTables();
106 112
107 // <<
108 } <<
109 113
110 //....oooOO0OOooo........oooOO0OOooo........oo << 114 //
>> 115 }
>> 116 ////////////////////////////////////////////////////////////////////////////////
>> 117 //
111 // Destructor 118 // Destructor
112 // 119 //
113 G4UrQMD1_3Model::~G4UrQMD1_3Model() {} << 120 G4UrQMD1_3Model::~G4UrQMD1_3Model (){}
>> 121 ////////////////////////////////////////////////////////////////////////////////
114 122
115 //....oooOO0OOooo........oooOO0OOooo........oo << 123 G4ReactionProductVector* G4UrQMD1_3Model::Propagate(G4KineticTrackVector* ,
116 << 124 G4V3DNucleus* ) {
117 G4ReactionProductVector* G4UrQMD1_3Model::Prop <<
118 { <<
119 return 0; 125 return 0;
120 } 126 }
121 127
122 //....oooOO0OOooo........oooOO0OOooo........oo << 128 ////////////////////////////////////////////////////////////////////////////////
123 // 129 //
124 // ApplyYourself 130 // ApplyYourself
125 // 131 //
126 // Member function to process an event, and ge 132 // Member function to process an event, and get information about the products.
127 133
128 G4HadFinalState* G4UrQMD1_3Model::ApplyYoursel << 134
129 << 135 G4HadFinalState *G4UrQMD1_3Model::ApplyYourself (
>> 136 const G4HadProjectile &theTrack, G4Nucleus &theTarget)
130 { 137 {
131 // anti_new << 138 //anti_new
132 // ------------------define anti_light_nuc << 139 // ------------------define anti_light_nucleus
133 const G4ParticleDefinition* anti_deu = G4Ant << 140 const G4ParticleDefinition* anti_deu =
134 << 141 G4AntiDeuteron::AntiDeuteron();
135 const G4ParticleDefinition* anti_he3 = G4Ant << 142
136 << 143 const G4ParticleDefinition* anti_he3=
137 const G4ParticleDefinition* anti_tri = G4Ant << 144 G4AntiHe3::AntiHe3();
138 << 145
139 const G4ParticleDefinition* anti_alp = G4Ant << 146 const G4ParticleDefinition* anti_tri=
140 //------------------------------------------ << 147 G4AntiTriton::AntiTriton();
141 // << 148
142 // The secondaries will be returned in G4Had << 149 const G4ParticleDefinition* anti_alp=
143 // initialise this. The original track will << 150 G4AntiAlpha::AntiAlpha();
144 // secondaries followed. << 151 //---------------------------------------------------
145 // << 152 //
>> 153 // The secondaries will be returned in G4HadFinalState &theResult -
>> 154 // initialise this. The original track will always be discontinued and
>> 155 // secondaries followed.
>> 156 //
146 theResult.Clear(); 157 theResult.Clear();
147 theResult.SetStatusChange(stopAndKill); 158 theResult.SetStatusChange(stopAndKill);
148 159
149 G4DynamicParticle* cascadeParticle = 0; << 160 G4DynamicParticle* cascadeParticle=0;
150 // << 161 //
151 // << 162 //
152 // Get relevant information about the projec << 163 // Get relevant information about the projectile and target (A, Z, energy/nuc,
153 // momentum, etc). << 164 // momentum, etc).
154 // << 165 //
155 << 166
156 const G4ParticleDefinition* definitionP = th << 167
157 const G4double AP = definitionP->GetBaryonNu << 168 const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
158 const G4double ZP = definitionP->GetPDGCharg << 169 const G4double AP = definitionP->GetBaryonNumber();
159 G4double AT = theTarget.GetN(); << 170 const G4double ZP = definitionP->GetPDGCharge();
160 G4double ZT = theTarget.GetZ(); << 171 G4double AT = theTarget.GetN();
161 // ---------------------------------------- << 172 G4double ZT = theTarget.GetZ();
162 G4int id = definitionP->GetPDGEncoding(); / << 173 // -----------------------------------------------
163 // ----------------------------------------- << 174 G4int id=definitionP->GetPDGEncoding(); //get particle encoding
>> 175 // ------------------------------------------------
164 G4int AP1 = G4lrint(AP); 176 G4int AP1 = G4lrint(AP);
165 G4int ZP1 = G4lrint(ZP); 177 G4int ZP1 = G4lrint(ZP);
166 G4int AT1 = G4lrint(AT); 178 G4int AT1 = G4lrint(AT);
167 G4int ZT1 = G4lrint(ZT); 179 G4int ZT1 = G4lrint(ZT);
168 // G4cout<<"------ap1--=="<<AP1<<"---zp1--- << 180 // G4cout<<"------ap1--=="<<AP1<<"---zp1---=="<<ZP1<<"---id-=="<<id<< G4endl;
169 // << 181 //
170 // ***************************************** << 182 // ****************************************************************************
171 // The following is the parameters necessary << 183 // The following is the parameters necessary to initiate Uinit() and UrQMD()
172 // ----------------------------------------- << 184 // ----------------------------------------------------------------------------
173 urqmdparams_.u_sptar = 0; //! 0= normal pro << 185 urqmdparams_.u_sptar=0; //!0= normal proj/target, 1=special proj/tar
174 urqmdparams_.u_spproj = 1; // projectile is << 186 urqmdparams_.u_spproj=1; // projectile is a special particle
175 << 187
176 // new_anti << 188 //new_anti
>> 189
>> 190 if (AP1>1 ||definitionP==anti_deu ||definitionP==anti_he3
>> 191 ||definitionP==anti_tri ||definitionP==anti_alp)
>> 192 {
>> 193
>> 194 urqmdparams_.u_ap=AP1;
>> 195 urqmdparams_.u_zp=ZP1;
>> 196
>> 197 urqmdparams_.u_spproj=0;
>> 198 } else if (id==2212) { //!proton
>> 199 urqmdparams_.u_ap=1;
>> 200 urqmdparams_.u_zp=1;
>> 201
>> 202
>> 203 } else if(id==-2212){ //! anti-proton
>> 204 urqmdparams_.u_ap=-1;
>> 205 urqmdparams_.u_zp=-1;
>> 206 } else if(id==2112){ //! neutron
>> 207 urqmdparams_.u_ap=1;
>> 208 urqmdparams_.u_zp=-1;
>> 209
>> 210 } else if(id==-2112){ //! anti-neutron
>> 211 urqmdparams_.u_ap=-1;
>> 212 urqmdparams_.u_zp=1;
>> 213
>> 214 } else if(id==211) { //! pi+
>> 215 urqmdparams_.u_ap=101;
>> 216 urqmdparams_.u_zp=2;
>> 217 } else if(id==-211) { //! pi-
>> 218 urqmdparams_.u_ap=101;
>> 219 urqmdparams_.u_zp=-2;
>> 220 } else if(id==321) { //! K+
>> 221 urqmdparams_.u_ap=106;
>> 222 urqmdparams_.u_zp=1;
>> 223 } else if(id==-321) { //! K-
>> 224 urqmdparams_.u_ap=-106;
>> 225 urqmdparams_.u_zp=-1;
>> 226 } else if(id==130 || id==310) { // ! K0
>> 227 urqmdparams_.u_ap=106;
>> 228 urqmdparams_.u_zp=-1;
>> 229 } else if(id==-130 || id==-310){ // ! K0bar
>> 230 urqmdparams_.u_ap=-106;
>> 231 urqmdparams_.u_zp=1;
>> 232 } else {
>> 233
>> 234 G4cout << " Sorry, No definition for particle for UrQMD::"
>> 235 <<id<< "found" << G4endl;
>> 236
>> 237 //AND->
>> 238 #if G4VERSION_NUMBER>=950
>> 239 //New signature (9.5) for G4Exception
>> 240 //Using G4HadronicException
>> 241 throw G4HadronicException(__FILE__,__LINE__,
>> 242 "Sorry, no definition for particle for UrQMD");
>> 243 #else
>> 244 G4Exception(" ");
>> 245 #endif
>> 246 //AND<-
>> 247 } //end if id
>> 248 //-------------------------------------------------------
>> 249
>> 250 urqmdparams_.u_at=AT1; // Target identified
>> 251 urqmdparams_.u_zt=ZT1;
>> 252 //----------------------------------------------------
>> 253 // identify Energy
>> 254 //
>> 255 G4ThreeVector Pbefore = theTrack.Get4Momentum().vect();
>> 256 G4double T = theTrack.GetKineticEnergy();
>> 257 G4double E = theTrack.GetTotalEnergy();
>> 258 G4double TotalEbefore = E*AP1 +
>> 259 theTarget.AtomicMass(AT1, ZT1) + theTarget.GetEnergyDeposit();
>> 260 // -----------------------------------------------------------------
>> 261
>> 262 if (AP1>1) {
>> 263 urqmdparams_.u_elab=T/(AP1*GeV); // Units are GeV/nuc for UrQMD
>> 264
>> 265 E = E/AP1; // Units are GeV/nuc
>> 266
>> 267 } else {
>> 268
>> 269 urqmdparams_.u_elab=T/GeV; //units are GeV
>> 270
>> 271 TotalEbefore = E +
>> 272 theTarget.AtomicMass(AT1, ZT1) + theTarget.GetEnergyDeposit();
>> 273 }
>> 274
>> 275
>> 276 //------------------------------------------------------------
>> 277 // identify impact parameter
>> 278 urqmdparams_.u_imp=-(1.1 * std::pow(G4double(AT1),(1./3.)));
>> 279 //units are in fm for UrQMD;
>> 280 //------------------------------------------------------------
>> 281 ///////////////////////// initialise/////////////////////
177 282
178 if (AP1 > 1 || definitionP == anti_deu || de << 283 if (CurrentEvent==0)
179 || definitionP == anti_alp) << 284 {
180 { << 285 G4cout << "\n creation of table, wait-------"<<G4endl;
181 urqmdparams_.u_ap = AP1; <<
182 urqmdparams_.u_zp = ZP1; <<
183 286
184 urqmdparams_.u_spproj = 0; << 287 G4cout << "\n"<<G4endl;
185 } <<
186 else if (id == 2212) { //! proton <<
187 urqmdparams_.u_ap = 1; <<
188 urqmdparams_.u_zp = 1; <<
189 } <<
190 else if (id == -2212) { //! anti-proton <<
191 urqmdparams_.u_ap = -1; <<
192 urqmdparams_.u_zp = -1; <<
193 } <<
194 else if (id == 2112) { //! neutron <<
195 urqmdparams_.u_ap = 1; <<
196 urqmdparams_.u_zp = -1; <<
197 } <<
198 else if (id == -2112) { //! anti-neutron <<
199 urqmdparams_.u_ap = -1; <<
200 urqmdparams_.u_zp = 1; <<
201 } <<
202 else if (id == 211) { //! pi+ <<
203 urqmdparams_.u_ap = 101; <<
204 urqmdparams_.u_zp = 2; <<
205 } <<
206 else if (id == -211) { //! pi- <<
207 urqmdparams_.u_ap = 101; <<
208 urqmdparams_.u_zp = -2; <<
209 } <<
210 else if (id == 321) { //! K+ <<
211 urqmdparams_.u_ap = 106; <<
212 urqmdparams_.u_zp = 1; <<
213 } <<
214 else if (id == -321) { //! K- <<
215 urqmdparams_.u_ap = -106; <<
216 urqmdparams_.u_zp = -1; <<
217 } <<
218 else if (id == 130 || id == 310) { // ! K0 <<
219 urqmdparams_.u_ap = 106; <<
220 urqmdparams_.u_zp = -1; <<
221 } <<
222 else if (id == -130 || id == -310) { // ! K <<
223 urqmdparams_.u_ap = -106; <<
224 urqmdparams_.u_zp = 1; <<
225 } <<
226 else { <<
227 G4cout << " Sorry, No definition for parti <<
228 288
229 // AND-> << 289 G4int io=0;
230 # if G4VERSION_NUMBER >= 950 <<
231 // New signature (9.5) for G4Exception <<
232 // Using G4HadronicException <<
233 throw G4HadronicException(__FILE__, __LINE <<
234 # else <<
235 G4Exception(" "); <<
236 # endif <<
237 // AND<- <<
238 } // end if id <<
239 //------------------------------------------ <<
240 <<
241 urqmdparams_.u_at = AT1; // Target identifi <<
242 urqmdparams_.u_zt = ZT1; <<
243 //------------------------------------------ <<
244 // identify Energy <<
245 // <<
246 G4ThreeVector Pbefore = theTrack.Get4Momentu <<
247 G4double T = theTrack.GetKineticEnergy(); <<
248 G4double E = theTrack.GetTotalEnergy(); <<
249 G4double TotalEbefore = E * AP1 + theTarget. <<
250 // -------------------------------------- <<
251 290
252 if (AP1 > 1) { << 291 uinit_ (&io);
253 urqmdparams_.u_elab = T / (AP1 * GeV); // <<
254 292
255 E = E / AP1; // Units are GeV/nuc <<
256 } <<
257 else { <<
258 urqmdparams_.u_elab = T / GeV; // units a <<
259 293
260 TotalEbefore = E + theTarget.AtomicMass(AT << 294 G4cout << "\n end to create table "<<G4endl;
261 } <<
262 295
263 //------------------------------------------ << 296 CurrentEvent=1;
264 // identify impact parameter << 297 }
265 urqmdparams_.u_imp = -(1.1 * std::pow(G4doub << 298 ////////////////////////////////////////////////////////
266 // units are in fm for UrQMD; <<
267 //------------------------------------------ <<
268 ///////////////////////// initialise//////// <<
269 299
270 if (CurrentEvent == 0) { << 300 //#ifdef debug_G4UrQMD1_3Model
271 G4cout << "\n creation of table, wait----- <<
272 301
273 G4cout << "\n" << G4endl; << 302 G4cout <<"UrQMDModel running-------------" <<G4endl;
274 303
275 G4int io = 0; << 304 urqmd_ ();
276 305
277 uinit_(&io); << 306 //#endif
278 307
279 G4cout << "\n end to create table " << G4 << 308 //G4cout <<"Number of produced particles: " <<sys_.npart<<G4endl;
280 309
281 CurrentEvent = 1; << 310 G4int n = sys_.npart; //no of produced particles
282 } << 311 if (n<2)
283 //////////////////////////////////////////// << 312 {
>> 313 G4cout <<"===============Warning================"<<G4endl;
>> 314 G4cout <<"======================================"<<G4endl;
284 315
285 // #ifdef debug_G4UrQMD1_3Model << 316 G4cout <<"Number of produced particles is very low: " <<sys_.npart<<G4endl;
>> 317 G4cout <<"============================================"<<G4endl;
286 318
287 G4cout << "UrQMDModel running-------------" << 319 //AND->
>> 320 #if G4VERSION_NUMBER>=950
>> 321 //New signature (9.5) for G4Exception
>> 322 //Using G4HadronicException instead of base class
>> 323 throw G4HadronicException(__FILE__,__LINE__,
>> 324 "Number of produced particle is very low");
>> 325 #else
>> 326 G4Exception(" "); //stop
>> 327 #endif
>> 328 //AND<-
>> 329 } else {
>> 330 for (G4int i=0; i<n; i++)
>> 331 {
>> 332
288 333
289 urqmd_(); << 334 G4int pid=pdgid_ (&isys_.ityp[i], &isys_.iso3[i]);
290 335
291 // #endif << 336 // Particle is a final state secondary and not a nucleus.
>> 337 // Determine what this secondary particle is, and if valid, load dynamic
>> 338 // parameters.
>> 339 //
292 340
293 // G4cout <<"Number of produced particles: <<
294 341
295 G4int n = sys_.npart; // no of produced par << 342 G4ParticleDefinition* pd=
296 if (n < 2) { << 343 G4ParticleTable::GetParticleTable()->FindParticle(pid);
297 G4cout << "===============Warning========= <<
298 G4cout << "=============================== <<
299 344
300 G4cout << "Number of produced particles is << 345 if (pd)
301 G4cout << "=============================== << 346 {
>> 347 G4double px = (coor_.px[i]+ffermi_.ffermpx[i])* GeV;
>> 348 //units are in MeV/c for G4
>> 349 G4double py = (coor_.py[i]+ffermi_.ffermpy[i])* GeV;
>> 350 G4double pz = (coor_.pz[i]+ffermi_.ffermpz[i])* GeV;
302 351
303 // AND-> << 352 G4double et = (coor_.p0[i]) *GeV;
304 # if G4VERSION_NUMBER >= 950 <<
305 // New signature (9.5) for G4Exception <<
306 // Using G4HadronicException instead of ba <<
307 throw G4HadronicException(__FILE__, __LINE <<
308 # else <<
309 G4Exception(" "); // stop <<
310 # endif <<
311 // AND<- <<
312 } <<
313 else { <<
314 for (G4int i = 0; i < n; i++) { <<
315 G4int pid = pdgid_(&isys_.ityp[i], &isys <<
316 353
317 // Particle is a final state secondary a <<
318 // Determine what this secondary particl <<
319 // parameters. <<
320 // <<
321 354
322 G4ParticleDefinition* pd = G4ParticleTab << 355 // ------------------------------Use only "Lorentz vector"----------
323 356
324 if (pd) { << 357 G4LorentzVector lorenzvec = G4LorentzVector(px,py,pz,et);
325 G4double px = (coor_.px[i] + ffermi_.f <<
326 // units are in MeV/c for G4 <<
327 G4double py = (coor_.py[i] + ffermi_.f <<
328 G4double pz = (coor_.pz[i] + ffermi_.f <<
329 358
330 G4double et = (coor_.p0[i]) * GeV; << 359 cascadeParticle = new G4DynamicParticle(pd, lorenzvec); //
331 360
332 // ------------------------------Us << 361 theResult.AddSecondary(cascadeParticle);
333 362
334 G4LorentzVector lorenzvec = G4LorentzV << 363 //======================================================================
335 364
336 cascadeParticle = new G4DynamicParticl <<
337 365
338 theResult.AddSecondary(cascadeParticle <<
339 366
340 //==================================== << 367 } //if
>> 368 } //for
341 369
342 } // if << 370 } //if warning
343 } // for <<
344 371
345 } // if warning << 372 //=======================================================================
>> 373 if (verbose >= 3) {
346 374
347 //========================================== << 375 //
348 if (verbose >= 3) { <<
349 // <<
350 G4double TotalEafter = 0.0; 376 G4double TotalEafter = 0.0;
351 G4ThreeVector TotalPafter; 377 G4ThreeVector TotalPafter;
352 G4double charge = 0.0; << 378 G4double charge = 0.0;
353 G4int baryon = 0; << 379 G4int baryon = 0;
354 G4int nSecondaries = theResult.GetNumberOf << 380 G4int nSecondaries = theResult.GetNumberOfSecondaries();
355 << 381
356 for (G4int j = 0; j < nSecondaries; j++) { << 382 for (G4int j=0; j<nSecondaries; j++) {
357 TotalEafter += theResult.GetSecondary(j) << 383 TotalEafter += theResult.GetSecondary(j)->
>> 384 GetParticle()->GetTotalEnergy();
358 385
359 TotalPafter += theResult.GetSecondary(j) << 386 TotalPafter += theResult.GetSecondary(j)->
>> 387 GetParticle()->GetMomentum();
360 388
361 G4ParticleDefinition* pd = theResult.Get << 389 G4ParticleDefinition *pd = theResult.GetSecondary(j)->
>> 390 GetParticle()->GetDefinition();
362 391
363 charge += pd->GetPDGCharge(); 392 charge += pd->GetPDGCharge();
364 baryon += pd->GetBaryonNumber(); 393 baryon += pd->GetBaryonNumber();
365 394
366 } // for secondaries << 395 } //for secondaries
367 396
368 G4cout << "------------------------------- << 397 G4cout <<"----------------------------------------"
369 << "------------------------------- << 398 <<"----------------------------------------"
370 G4cout << "Total energy before collision << 399 <<G4endl;
371 << " MeV" << G4endl; << 400 G4cout <<"Total energy before collision = " <<TotalEbefore ///MeV
372 G4cout << "Total energy after collision << 401 <<" MeV" <<G4endl;
373 << " MeV" << G4endl; << 402 G4cout <<"Total energy after collision = " <<TotalEafter //MeV
374 << 403 <<" MeV" <<G4endl;
375 G4cout << "------------------------------- << 404
376 << 405 G4cout <<"----------------------------------------"<<G4endl;
377 G4cout << "Total momentum before collision << 406
378 << " MeV/c" << G4endl; << 407 G4cout <<"Total momentum before collision = " <<Pbefore //MeV
379 G4cout << "Total momentum after collision << 408 <<" MeV/c" <<G4endl;
380 << " MeV/c" << G4endl; << 409 G4cout <<"Total momentum after collision = " <<TotalPafter //MeV
381 G4cout << "------------------------------- << 410 <<" MeV/c" <<G4endl;
>> 411 G4cout <<"----------------------------------------"<<G4endl;
382 412
383 if (verbose >= 4) { 413 if (verbose >= 4) {
384 G4cout << "Total charge before collision << 414 G4cout <<"Total charge before collision = " <<(ZP+ZT)*eplus
385 G4cout << "Total charge after collision << 415 <<G4endl;
>> 416 G4cout <<"Total charge after collision = " <<charge
>> 417 <<G4endl;
386 418
387 G4cout << "----------------------------- << 419 G4cout <<"----------------------------------------"<<G4endl;
388 420
389 G4cout << "Total baryon number before co << 421 G4cout <<"Total baryon number before collision = "<<AP+AT
390 G4cout << "Total baryon number after col << 422 <<G4endl;
391 G4cout << "----------------------------- << 423 G4cout <<"Total baryon number after collision = "<<baryon
>> 424 <<G4endl;
>> 425 G4cout <<"----------------------------------------"<<G4endl;
392 426
393 } // if verbose4 << 427 } //if verbose4
394 428
395 G4cout << "------------------------------- << 429 G4cout <<"----------------------------------------"
396 << "------------------------------- << 430 <<"----------------------------------------"
>> 431 <<G4endl;
397 432
398 } // if verbose3 << 433 } //if verbose3
399 434
400 return &theResult; <<
401 } // G4hadfinal <<
402 435
403 //....oooOO0OOooo........oooOO0OOooo........oo << 436 return &theResult;
>> 437 } //G4hadfinal
>> 438
>> 439
>> 440 //---------------------------------------------------------------------
>> 441
>> 442 //---------------------------------------------------------------------
404 // 443 //
405 // WelcomeMessage 444 // WelcomeMessage
406 // 445 //
407 void G4UrQMD1_3Model::WelcomeMessage() const << 446 void G4UrQMD1_3Model::WelcomeMessage () const
408 { 447 {
409 G4cout << G4endl; << 448 G4cout <<G4endl;
410 G4cout << " ******************************** << 449 G4cout <<" *****************************************************************"
411 G4cout << " Interface to G4UrQMD_1.3 << 450 <<G4endl;
412 G4cout << " Version number : 00.00.0B << 451 G4cout <<" Interface to G4UrQMD_1.3 activated"
413 G4cout << " (Interface written by Kh. Abdel- << 452 <<G4endl;
414 G4cout << G4endl; << 453 G4cout <<" Version number : 00.00.0B File date : 25/01/12" <<G4endl;
415 G4cout << " ******************************** << 454 G4cout <<" (Interface written by Kh. Abdel-Waged et al. for the KACST/NCMP)"
>> 455 <<G4endl;
>> 456 G4cout <<G4endl;
>> 457 G4cout <<" *****************************************************************"
>> 458 <<G4endl;
416 G4cout << G4endl; 459 G4cout << G4endl;
417 460
418 return; 461 return;
419 } 462 }
420 463
421 //....oooOO0OOooo........oooOO0OOooo........oo << 464 void G4UrQMD1_3Model::InitialiseDataTables ()
422 <<
423 void G4UrQMD1_3Model::InitialiseDataTables() <<
424 { 465 {
425 // << 466 //
426 // << 467 //
427 // The next line is to make sure the block d << 468 // The next line is to make sure the block data statements are
428 // executed. << 469 // executed.
429 // << 470 //
>> 471
>> 472 g4urqmdblockdata_ ();
>> 473
430 474
431 g4urqmdblockdata_(); << 475 ///////////////////////////////////////////////////
>> 476 /////// Dynamic seed //////////////////////////////
>> 477 //G4int ranseed=-time_ ();
>> 478 // Fixed seed ///////////////////////////
432 479
433 //////////////////////////////////////////// << 480 G4int ranseed=1097569630;
434 /////// Dynamic seed /////////////////////// <<
435 // G4int ranseed=-time_ (); <<
436 // Fixed seed ///////////////////////// <<
437 481
438 G4int ranseed = 1097569630; << 482 G4cout <<"\n seed: "<<ranseed<<G4endl;
439 483
440 G4cout << "\n seed: " << ranseed << G4endl; << 484 sseed_ (&ranseed);
441 485
442 sseed_(&ranseed); << 486 loginit_();
443 487
444 loginit_(); <<
445 } 488 }
446 489
447 #endif // G4_USE_URQMD << 490 #endif //G4_USE_URQMD
448 491