Geant4_10
G4UrQMD1_3Model.cc
Go to the documentation of this file.
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 // * *
21 // * Parts of this code which have been developed by Abdel-Waged *
22 // * et al under contract (31-465) to the King Abdul-Aziz City for *
23 // * Science and Technology (KACST), the National Centre of *
24 // * Mathematics and Physics (NCMP), Saudi Arabia. *
25 // * *
26 // * By using, copying, modifying or distributing the software (or *
27 // * any work based on the software) you agree to acknowledge its *
28 // * use in resulting scientific publications, and indicate your *
29 // * acceptance of all terms of the Geant4 Software license. *
30 // ********************************************************************
31 //
34 //
35 // $Id: G4UrQMD1_3Model.cc 77519 2013-11-25 10:54:57Z gcosmo $
36 //
38 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39 //
40 // MODULE: G4UrQMD1_3Model.hh
41 //
42 // Version: 0.B
43 // Date: 25/01/12
44 // Authors: Kh. Abdel-Waged and Nuha Felemban
45 // Revised by: V.V. Uzhinskii
46 // SPONSERED BY
47 // Customer: KAUST/NCMP
48 // Contract: 31-465
49 //
50 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 //
52 #ifdef G4_USE_URQMD
53 
54 #include "G4UrQMD1_3Model.hh"
55 #include "G4UrQMD1_3Interface.hh"
56 //-------------------------------
57 #include "globals.hh"
58 #include "G4DynamicParticle.hh"
59 #include "G4IonTable.hh"
60 #include "G4CollisionOutput.hh"
61 #include "G4V3DNucleus.hh"
62 #include "G4Track.hh"
63 #include "G4Nucleus.hh"
64 #include "G4LorentzRotation.hh"
65 
66 #include "G4ParticleDefinition.hh"
67 #include "G4ParticleTable.hh"
68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh"
70 
71 //AND->
72 #include "G4Version.hh"
73 //AND<-
74 //----------------new_anti
75 #include "G4AntiHe3.hh"
76 #include "G4AntiDeuteron.hh"
77 #include "G4AntiTriton.hh"
78 #include "G4AntiAlpha.hh"
79 //---------------------------
80 #include <fstream>
81 #include <string>
82 
84 
85 //
87  :G4VIntraNuclearTransportModel(nam), verbose(0)
88 {
89 
90 
91  if (verbose > 3) {
92  G4cout << " >>> G4UrQMD1_3Model default constructor" << G4endl;
93  }
94 
95 
96 
97 //
98 // Set the minimum and maximum range for the UrQMD model
99 
100 // SetMinEnergy(100.0*MeV);
101 // SetMaxEnergy(200.0*GeV);
102 
103 //
104 
105 //
106  WelcomeMessage();
107 //
108  CurrentEvent=0;
109 //
110 
111 InitialiseDataTables();
112 
113 
114 //
115 }
117 //
118 // Destructor
119 //
122 
124  G4V3DNucleus* ) {
125  return 0;
126 }
127 
129 //
130 // ApplyYourself
131 //
132 // Member function to process an event, and get information about the products.
133 
134 
136  const G4HadProjectile &theTrack, G4Nucleus &theTarget)
137 {
138 //anti_new
139 // ------------------define anti_light_nucleus
140 const G4ParticleDefinition* anti_deu =
142 
143 const G4ParticleDefinition* anti_he3=
145 
146 const G4ParticleDefinition* anti_tri=
148 
149 const G4ParticleDefinition* anti_alp=
151 //---------------------------------------------------
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 //
157  theResult.Clear();
158  theResult.SetStatusChange(stopAndKill);
159 
160  G4DynamicParticle* cascadeParticle=0;
161 //
162 //
163 // Get relevant information about the projectile and target (A, Z, energy/nuc,
164 // momentum, etc).
165 //
166 
167 
168  const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
169  const G4double AP = definitionP->GetBaryonNumber();
170  const G4double ZP = definitionP->GetPDGCharge();
171  G4double AT = theTarget.GetN();
172  G4double ZT = theTarget.GetZ();
173 // -----------------------------------------------
174  G4int id=definitionP->GetPDGEncoding(); //get particle encoding
175 // ------------------------------------------------
176  G4int AP1 = G4lrint(AP);
177  G4int ZP1 = G4lrint(ZP);
178  G4int AT1 = G4lrint(AT);
179  G4int ZT1 = G4lrint(ZT);
180 // G4cout<<"------ap1--=="<<AP1<<"---zp1---=="<<ZP1<<"---id-=="<<id<< G4endl;
181 //
182 // ****************************************************************************
183 // The following is the parameters necessary to initiate Uinit() and UrQMD()
184 // ----------------------------------------------------------------------------
186  urqmdparams_.u_spproj=1; // projectile is a special particle
187 
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 
198  } else if (id==2212) {
199  urqmdparams_.u_ap=1;
200  urqmdparams_.u_zp=1;
201 
202 
203  } else if(id==-2212){
204  urqmdparams_.u_ap=-1;
205  urqmdparams_.u_zp=-1;
206  } else if(id==2112){
207  urqmdparams_.u_ap=1;
208  urqmdparams_.u_zp=-1;
209 
210  } else if(id==-2112){
211  urqmdparams_.u_ap=-1;
212  urqmdparams_.u_zp=1;
213 
214  } else if(id==211) {
215  urqmdparams_.u_ap=101;
216  urqmdparams_.u_zp=2;
217  } else if(id==-211) {
218  urqmdparams_.u_ap=101;
219  urqmdparams_.u_zp=-2;
220  } else if(id==321) {
221  urqmdparams_.u_ap=106;
222  urqmdparams_.u_zp=1;
223  } else if(id==-321) {
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 //------------------------------------------------------------
282 
283 if (CurrentEvent==0)
284 {
285 G4cout << "\n creation of table, wait-------"<<G4endl;
286 
287 G4cout << "\n"<<G4endl;
288 
289 G4int io=0;
290 
291 uinit_ (&io);
292 
293 
294 G4cout << "\n end to create table "<<G4endl;
295 
296 CurrentEvent=1;
297 }
299 
300 //#ifdef debug_G4UrQMD1_3Model
301 
302 G4cout <<"UrQMDModel running-------------" <<G4endl;
303 
304 urqmd_ ();
305 
306 //#endif
307 
308 //G4cout <<"Number of produced particles: " <<sys_.npart<<G4endl;
309 
310 G4int n = sys_.npart; //no of produced particles
311 if (n<2)
312 {
313 G4cout <<"===============Warning================"<<G4endl;
314 G4cout <<"======================================"<<G4endl;
315 
316 G4cout <<"Number of produced particles is very low: " <<sys_.npart<<G4endl;
317 G4cout <<"============================================"<<G4endl;
318 
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 
333 
334 G4int pid=pdgid_ (&isys_.ityp[i], &isys_.iso3[i]);
335 
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 //
340 
341 
344 
345  if (pd)
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;
351 
352  G4double et = (coor_.p0[i]) *GeV;
353 
354 
355 // ------------------------------Use only "Lorentz vector"----------
356 
357  G4LorentzVector lorenzvec = G4LorentzVector(px,py,pz,et);
358 
359  cascadeParticle = new G4DynamicParticle(pd, lorenzvec); //
360 
361  theResult.AddSecondary(cascadeParticle);
362 
363 //======================================================================
364 
365 
366 
367  } //if
368 } //for
369 
370 } //if warning
371 
372 //=======================================================================
373 if (verbose >= 3) {
374 
375 //
376  G4double TotalEafter = 0.0;
377  G4ThreeVector TotalPafter;
378  G4double charge = 0.0;
379  G4int baryon = 0;
380  G4int nSecondaries = theResult.GetNumberOfSecondaries();
381 
382  for (G4int j=0; j<nSecondaries; j++) {
383  TotalEafter += theResult.GetSecondary(j)->
384  GetParticle()->GetTotalEnergy();
385 
386  TotalPafter += theResult.GetSecondary(j)->
387  GetParticle()->GetMomentum();
388 
389  G4ParticleDefinition *pd = theResult.GetSecondary(j)->
390  GetParticle()->GetDefinition();
391 
392  charge += pd->GetPDGCharge();
393  baryon += pd->GetBaryonNumber();
394 
395  } //for secondaries
396 
397  G4cout <<"----------------------------------------"
398  <<"----------------------------------------"
399  <<G4endl;
400  G4cout <<"Total energy before collision = " <<TotalEbefore
401  <<" MeV" <<G4endl;
402  G4cout <<"Total energy after collision = " <<TotalEafter //MeV
403  <<" MeV" <<G4endl;
404 
405  G4cout <<"----------------------------------------"<<G4endl;
406 
407  G4cout <<"Total momentum before collision = " <<Pbefore //MeV
408  <<" MeV/c" <<G4endl;
409  G4cout <<"Total momentum after collision = " <<TotalPafter //MeV
410  <<" MeV/c" <<G4endl;
411  G4cout <<"----------------------------------------"<<G4endl;
412 
413  if (verbose >= 4) {
414  G4cout <<"Total charge before collision = " <<(ZP+ZT)*eplus
415  <<G4endl;
416  G4cout <<"Total charge after collision = " <<charge
417  <<G4endl;
418 
419  G4cout <<"----------------------------------------"<<G4endl;
420 
421  G4cout <<"Total baryon number before collision = "<<AP+AT
422  <<G4endl;
423  G4cout <<"Total baryon number after collision = "<<baryon
424  <<G4endl;
425  G4cout <<"----------------------------------------"<<G4endl;
426 
427  } //if verbose4
428 
429  G4cout <<"----------------------------------------"
430  <<"----------------------------------------"
431  <<G4endl;
432 
433  } //if verbose3
434 
435 
436 return &theResult;
437 } //G4hadfinal
438 
439 
440 //---------------------------------------------------------------------
441 
442 //---------------------------------------------------------------------
443 //
444 // WelcomeMessage
445 //
446 void G4UrQMD1_3Model::WelcomeMessage () const
447 {
448  G4cout <<G4endl;
449  G4cout <<" *****************************************************************"
450  <<G4endl;
451  G4cout <<" Interface to G4UrQMD_1.3 activated"
452  <<G4endl;
453  G4cout <<" Version number : 00.00.0B File date : 25/01/12" <<G4endl;
454  G4cout <<" (Interface written by Kh. Abdel-Waged et al. for the KACST/NCMP)"
455  <<G4endl;
456  G4cout <<G4endl;
457  G4cout <<" *****************************************************************"
458  <<G4endl;
459  G4cout << G4endl;
460 
461  return;
462 }
463 
464 void G4UrQMD1_3Model::InitialiseDataTables ()
465 {
466 //
467 //
468 // The next line is to make sure the block data statements are
469 // executed.
470 //
471 
473 
474 
477 //G4int ranseed=-time_ ();
478 // Fixed seed ///////////////////////////
479 
480 G4int ranseed=1097569630;
481 
482 G4cout <<"\n seed: "<<ranseed<<G4endl;
483 
484 sseed_ (&ranseed);
485 
486 loginit_();
487 
488 }
489 
490 #endif //G4_USE_URQMD
static G4AntiHe3 * AntiHe3()
Definition: G4AntiHe3.cc:94
G4double AtomicMass(const G4double A, const G4double Z) const
Definition: G4Nucleus.cc:240
G4double py[nmax]
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4HadSecondary * GetSecondary(size_t i)
struct ccurqmd13sys sys_
void loginit_()
static G4AntiDeuteron * AntiDeuteron()
int pdgid_(int *, int *)
void uinit_(int *)
void sseed_(int *)
static G4AntiAlpha * AntiAlpha()
Definition: G4AntiAlpha.cc:89
Definition of the G4UrQMD1_3Model class.
int G4int
Definition: G4Types.hh:78
G4double p0[nmax]
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
Char_t n[5]
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
struct ccurqmd13ffermi ffermi_
struct ccurqmd13isys isys_
G4double GetKineticEnergy() const
G4ErrorTarget * theTarget
Definition: errprop.cc:59
G4double GetEnergyDeposit()
Definition: G4Nucleus.hh:184
struct ccurqmd13coor coor_
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theTarget)
const G4LorentzVector & Get4Momentum() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double px[nmax]
void g4urqmdblockdata_()
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &, G4Nucleus &)
virtual ~G4UrQMD1_3Model()
void urqmd_()
G4UrQMD1_3Model(const G4String &name="UrQMD1_3")
G4double ffermpz[nmax]
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
struct ccurqmd13urqmdparams urqmdparams_
static G4AntiTriton * AntiTriton()
Definition: G4AntiTriton.cc:94
G4double GetPDGCharge() const
Definition of the G4UrQMD1_3Interface class.
G4int GetNumberOfSecondaries() const
void AddSecondary(G4DynamicParticle *aP)
G4double ffermpx[nmax]
G4double pz[nmax]
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector
G4double ffermpy[nmax]