Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr02/src/G4UrQMD1_3Model.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //
 32 /// \file hadronic/Hadr02/src/G4UrQMD1_3Model.cc
 33 /// \brief Implementation of the G4UrQMD1_3Model class
 34 //
 35 //
 36 ///////////////////////////////////////////////////////////////////////////////
 37 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 38 //
 39 // MODULE:          G4UrQMD1_3Model.hh
 40 //
 41 // Version:          0.B
 42 // Date:           25/01/12
 43 // Authors:        Kh. Abdel-Waged and Nuha Felemban
 44 // Revised by:      V.V. Uzhinskii
 45 //                  SPONSERED BY
 46 // Customer:        KAUST/NCMP
 47 // Contract:        31-465
 48 //
 49 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 50 //
 51 #ifdef G4_USE_URQMD
 52 
 53 #  include "G4UrQMD1_3Model.hh"
 54 
 55 #  include "G4UrQMD1_3Interface.hh"
 56 //-------------------------------
 57 #  include "G4CollisionOutput.hh"
 58 #  include "G4DynamicParticle.hh"
 59 #  include "G4IonTable.hh"
 60 #  include "G4LorentzRotation.hh"
 61 #  include "G4Nucleus.hh"
 62 #  include "G4ParticleDefinition.hh"
 63 #  include "G4ParticleTable.hh"
 64 #  include "G4PhysicalConstants.hh"
 65 #  include "G4SystemOfUnits.hh"
 66 #  include "G4Track.hh"
 67 #  include "G4V3DNucleus.hh"
 68 #  include "globals.hh"
 69 
 70 // AND->
 71 #  include "G4Version.hh"
 72 // AND<-
 73 //----------------new_anti
 74 #  include "G4AntiAlpha.hh"
 75 #  include "G4AntiDeuteron.hh"
 76 #  include "G4AntiHe3.hh"
 77 #  include "G4AntiTriton.hh"
 78 //---------------------------
 79 #  include <fstream>
 80 #  include <string>
 81 
 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 83 //
 84 G4UrQMD1_3Model::G4UrQMD1_3Model(const G4String& nam)
 85   : G4VIntraNuclearTransportModel(nam), verbose(0)
 86 {
 87   if (verbose > 3) {
 88     G4cout << " >>> G4UrQMD1_3Model default constructor" << G4endl;
 89   }
 90 
 91   //
 92   // Set the minimum and maximum range for the UrQMD model
 93 
 94   //  SetMinEnergy(100.0*MeV);
 95   //  SetMaxEnergy(200.0*GeV);
 96 
 97   //
 98 
 99   //
100   WelcomeMessage();
101   //
102   CurrentEvent = 0;
103   //
104 
105   InitialiseDataTables();
106 
107   //
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 // Destructor
112 //
113 G4UrQMD1_3Model::~G4UrQMD1_3Model() {}
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
117 G4ReactionProductVector* G4UrQMD1_3Model::Propagate(G4KineticTrackVector*, G4V3DNucleus*)
118 {
119   return 0;
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123 //
124 // ApplyYourself
125 //
126 // Member function to process an event, and get information about the products.
127 
128 G4HadFinalState* G4UrQMD1_3Model::ApplyYourself(const G4HadProjectile& theTrack,
129                                                 G4Nucleus& theTarget)
130 {
131   // anti_new
132   //   ------------------define anti_light_nucleus
133   const G4ParticleDefinition* anti_deu = G4AntiDeuteron::AntiDeuteron();
134 
135   const G4ParticleDefinition* anti_he3 = G4AntiHe3::AntiHe3();
136 
137   const G4ParticleDefinition* anti_tri = G4AntiTriton::AntiTriton();
138 
139   const G4ParticleDefinition* anti_alp = G4AntiAlpha::AntiAlpha();
140   //---------------------------------------------------
141   //
142   // The secondaries will be returned in G4HadFinalState &theResult -
143   // initialise this.  The original track will always be discontinued and
144   // secondaries followed.
145   //
146   theResult.Clear();
147   theResult.SetStatusChange(stopAndKill);
148 
149   G4DynamicParticle* cascadeParticle = 0;
150   //
151   //
152   // Get relevant information about the projectile and target (A, Z, energy/nuc,
153   // momentum, etc).
154   //
155 
156   const G4ParticleDefinition* definitionP = theTrack.GetDefinition();
157   const G4double AP = definitionP->GetBaryonNumber();
158   const G4double ZP = definitionP->GetPDGCharge();
159   G4double AT = theTarget.GetN();
160   G4double ZT = theTarget.GetZ();
161   //  -----------------------------------------------
162   G4int id = definitionP->GetPDGEncoding();  // get particle encoding
163   // ------------------------------------------------
164   G4int AP1 = G4lrint(AP);
165   G4int ZP1 = G4lrint(ZP);
166   G4int AT1 = G4lrint(AT);
167   G4int ZT1 = G4lrint(ZT);
168   //  G4cout<<"------ap1--=="<<AP1<<"---zp1---=="<<ZP1<<"---id-=="<<id<< G4endl;
169   //
170   // ****************************************************************************
171   // The following is the parameters necessary to initiate Uinit() and UrQMD()
172   // ----------------------------------------------------------------------------
173   urqmdparams_.u_sptar = 0;  //! 0= normal proj/target, 1=special proj/tar
174   urqmdparams_.u_spproj = 1;  // projectile is a special particle
175 
176   // new_anti
177 
178   if (AP1 > 1 || definitionP == anti_deu || definitionP == anti_he3 || definitionP == anti_tri
179       || definitionP == anti_alp)
180   {
181     urqmdparams_.u_ap = AP1;
182     urqmdparams_.u_zp = ZP1;
183 
184     urqmdparams_.u_spproj = 0;
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) {  // ! K0bar
223     urqmdparams_.u_ap = -106;
224     urqmdparams_.u_zp = 1;
225   }
226   else {
227     G4cout << " Sorry, No definition for particle for UrQMD::" << id << "found" << G4endl;
228 
229     // AND->
230 #  if G4VERSION_NUMBER >= 950
231     // New signature (9.5) for G4Exception
232     // Using G4HadronicException
233     throw G4HadronicException(__FILE__, __LINE__, "Sorry, no definition for particle for UrQMD");
234 #  else
235     G4Exception(" ");
236 #  endif
237     // AND<-
238   }  // end if id
239   //-------------------------------------------------------
240 
241   urqmdparams_.u_at = AT1;  // Target identified
242   urqmdparams_.u_zt = ZT1;
243   //----------------------------------------------------
244   //  identify Energy
245   //
246   G4ThreeVector Pbefore = theTrack.Get4Momentum().vect();
247   G4double T = theTrack.GetKineticEnergy();
248   G4double E = theTrack.GetTotalEnergy();
249   G4double TotalEbefore = E * AP1 + theTarget.AtomicMass(AT1, ZT1) + theTarget.GetEnergyDeposit();
250   //    -----------------------------------------------------------------
251 
252   if (AP1 > 1) {
253     urqmdparams_.u_elab = T / (AP1 * GeV);  // Units are GeV/nuc for UrQMD
254 
255     E = E / AP1;  // Units are GeV/nuc
256   }
257   else {
258     urqmdparams_.u_elab = T / GeV;  // units are GeV
259 
260     TotalEbefore = E + theTarget.AtomicMass(AT1, ZT1) + theTarget.GetEnergyDeposit();
261   }
262 
263   //------------------------------------------------------------
264   // identify impact parameter
265   urqmdparams_.u_imp = -(1.1 * std::pow(G4double(AT1), (1. / 3.)));
266   // units are in fm for UrQMD;
267   //------------------------------------------------------------
268   ///////////////////////// initialise/////////////////////
269 
270   if (CurrentEvent == 0) {
271     G4cout << "\n creation of table, wait-------" << G4endl;
272 
273     G4cout << "\n" << G4endl;
274 
275     G4int io = 0;
276 
277     uinit_(&io);
278 
279     G4cout << "\n end to create  table " << G4endl;
280 
281     CurrentEvent = 1;
282   }
283   ////////////////////////////////////////////////////////
284 
285   // #ifdef debug_G4UrQMD1_3Model
286 
287   G4cout << "UrQMDModel running-------------" << G4endl;
288 
289   urqmd_();
290 
291   // #endif
292 
293   // G4cout <<"Number of produced particles:  " <<sys_.npart<<G4endl;
294 
295   G4int n = sys_.npart;  // no of produced particles
296   if (n < 2) {
297     G4cout << "===============Warning================" << G4endl;
298     G4cout << "======================================" << G4endl;
299 
300     G4cout << "Number of produced particles is very low:  " << sys_.npart << G4endl;
301     G4cout << "============================================" << G4endl;
302 
303 // AND->
304 #  if G4VERSION_NUMBER >= 950
305     // New signature (9.5) for G4Exception
306     // Using G4HadronicException instead of base class
307     throw G4HadronicException(__FILE__, __LINE__, "Number of produced particle is very low");
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_.iso3[i]);
316 
317       // Particle is a final state secondary and not a nucleus.
318       // Determine what this secondary particle is, and if valid, load dynamic
319       // parameters.
320       //
321 
322       G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle(pid);
323 
324       if (pd) {
325         G4double px = (coor_.px[i] + ffermi_.ffermpx[i]) * GeV;
326         // units are in MeV/c for G4
327         G4double py = (coor_.py[i] + ffermi_.ffermpy[i]) * GeV;
328         G4double pz = (coor_.pz[i] + ffermi_.ffermpz[i]) * GeV;
329 
330         G4double et = (coor_.p0[i]) * GeV;
331 
332         //    ------------------------------Use only "Lorentz vector"----------
333 
334         G4LorentzVector lorenzvec = G4LorentzVector(px, py, pz, et);
335 
336         cascadeParticle = new G4DynamicParticle(pd, lorenzvec);  //
337 
338         theResult.AddSecondary(cascadeParticle);
339 
340         //======================================================================
341 
342       }  // if
343     }  // for
344 
345   }  // if warning
346 
347   //=======================================================================
348   if (verbose >= 3) {
349     //
350     G4double TotalEafter = 0.0;
351     G4ThreeVector TotalPafter;
352     G4double charge = 0.0;
353     G4int baryon = 0;
354     G4int nSecondaries = theResult.GetNumberOfSecondaries();
355 
356     for (G4int j = 0; j < nSecondaries; j++) {
357       TotalEafter += theResult.GetSecondary(j)->GetParticle()->GetTotalEnergy();
358 
359       TotalPafter += theResult.GetSecondary(j)->GetParticle()->GetMomentum();
360 
361       G4ParticleDefinition* pd = theResult.GetSecondary(j)->GetParticle()->GetDefinition();
362 
363       charge += pd->GetPDGCharge();
364       baryon += pd->GetBaryonNumber();
365 
366     }  // for secondaries
367 
368     G4cout << "----------------------------------------"
369            << "----------------------------------------" << G4endl;
370     G4cout << "Total energy before collision  = " << TotalEbefore  /// MeV
371            << " MeV" << G4endl;
372     G4cout << "Total energy after collision    = " << TotalEafter  // MeV
373            << " MeV" << G4endl;
374 
375     G4cout << "----------------------------------------" << G4endl;
376 
377     G4cout << "Total momentum before collision = " << Pbefore  // MeV
378            << " MeV/c" << G4endl;
379     G4cout << "Total momentum after collision  = " << TotalPafter  // MeV
380            << " MeV/c" << G4endl;
381     G4cout << "----------------------------------------" << G4endl;
382 
383     if (verbose >= 4) {
384       G4cout << "Total charge before collision  = " << (ZP + ZT) * eplus << G4endl;
385       G4cout << "Total charge after collision    = " << charge << G4endl;
386 
387       G4cout << "----------------------------------------" << G4endl;
388 
389       G4cout << "Total baryon number before collision = " << AP + AT << G4endl;
390       G4cout << "Total baryon number after collision  = " << baryon << G4endl;
391       G4cout << "----------------------------------------" << G4endl;
392 
393     }  // if verbose4
394 
395     G4cout << "----------------------------------------"
396            << "----------------------------------------" << G4endl;
397 
398   }  // if verbose3
399 
400   return &theResult;
401 }  // G4hadfinal
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 //
405 // WelcomeMessage
406 //
407 void G4UrQMD1_3Model::WelcomeMessage() const
408 {
409   G4cout << G4endl;
410   G4cout << " *****************************************************************" << G4endl;
411   G4cout << " Interface to        G4UrQMD_1.3                      activated" << G4endl;
412   G4cout << " Version number : 00.00.0B          File date : 25/01/12" << G4endl;
413   G4cout << " (Interface written by Kh. Abdel-Waged et al. for the KACST/NCMP)" << G4endl;
414   G4cout << G4endl;
415   G4cout << " *****************************************************************" << G4endl;
416   G4cout << G4endl;
417 
418   return;
419 }
420 
421 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
422 
423 void G4UrQMD1_3Model::InitialiseDataTables()
424 {
425   //
426   //
427   // The next line is to make sure the block data statements are
428   // executed.
429   //
430 
431   g4urqmdblockdata_();
432 
433   ///////////////////////////////////////////////////
434   /////// Dynamic seed //////////////////////////////
435   // G4int ranseed=-time_ ();
436   //     Fixed seed  ///////////////////////////
437 
438   G4int ranseed = 1097569630;
439 
440   G4cout << "\n seed:  " << ranseed << G4endl;
441 
442   sseed_(&ranseed);
443 
444   loginit_();
445 }
446 
447 #endif  // G4_USE_URQMD
448