Hadr02/src/G4UrQMD1_3Model.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // *                                                                  *
0021 // * Parts of this code which have been  developed by Abdel-Waged     *
0022 // * et al under contract (31-465) to the King Abdul-Aziz City for    *
0023 // * Science and Technology (KACST), the National Centre of           *
0024 // * Mathematics and Physics (NCMP), Saudi Arabia.                    *
0025 // *                                                                  *
0026 // * By using,  copying,  modifying or  distributing the software (or *
0027 // * any work based  on the software)  you  agree  to acknowledge its *
0028 // * use  in  resulting  scientific  publications,  and indicate your *
0029 // * acceptance of all terms of the Geant4 Software license.          *
0030 // ********************************************************************
0031 //
0032 /// \file hadronic/Hadr02/src/G4UrQMD1_3Model.cc
0033 /// \brief Implementation of the G4UrQMD1_3Model class
0034 //
0035 //
0036 ///////////////////////////////////////////////////////////////////////////////
0037 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0038 //
0039 // MODULE:          G4UrQMD1_3Model.hh
0040 //
0041 // Version:          0.B
0042 // Date:           25/01/12
0043 // Authors:        Kh. Abdel-Waged and Nuha Felemban
0044 // Revised by:      V.V. Uzhinskii
0045 //                  SPONSERED BY
0046 // Customer:        KAUST/NCMP
0047 // Contract:        31-465
0048 //
0049 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 //
0051 #ifdef G4_USE_URQMD
0052
0053 #  include "G4UrQMD1_3Model.hh"
0054
0055 #  include "G4UrQMD1_3Interface.hh"
0056 //-------------------------------
0057 #  include "G4CollisionOutput.hh"
0058 #  include "G4DynamicParticle.hh"
0059 #  include "G4IonTable.hh"
0060 #  include "G4LorentzRotation.hh"
0061 #  include "G4Nucleus.hh"
0062 #  include "G4ParticleDefinition.hh"
0063 #  include "G4ParticleTable.hh"
0064 #  include "G4PhysicalConstants.hh"
0065 #  include "G4SystemOfUnits.hh"
0066 #  include "G4Track.hh"
0067 #  include "G4V3DNucleus.hh"
0068 #  include "globals.hh"
0069
0070 // AND->
0071 #  include "G4Version.hh"
0072 // AND<-
0073 //----------------new_anti
0074 #  include "G4AntiAlpha.hh"
0075 #  include "G4AntiDeuteron.hh"
0076 #  include "G4AntiHe3.hh"
0077 #  include "G4AntiTriton.hh"
0078 //---------------------------
0079 #  include <fstream>
0080 #  include <string>
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083 //
0084 G4UrQMD1_3Model::G4UrQMD1_3Model(const G4String& nam)
0085   : G4VIntraNuclearTransportModel(nam), verbose(0)
0086 {
0087   if (verbose > 3) {
0088     G4cout << " >>> G4UrQMD1_3Model default constructor" << G4endl;
0089   }
0090
0091   //
0092   // Set the minimum and maximum range for the UrQMD model
0093
0094   //  SetMinEnergy(100.0*MeV);
0095   //  SetMaxEnergy(200.0*GeV);
0096
0097   //
0098
0099   //
0100   WelcomeMessage();
0101   //
0102   CurrentEvent = 0;
0103   //
0104
0105   InitialiseDataTables();
0106
0107   //
0108 }
0109
0110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0111 // Destructor
0112 //
0113 G4UrQMD1_3Model::~G4UrQMD1_3Model() {}
0114
0115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0116
0117 G4ReactionProductVector* G4UrQMD1_3Model::Propagate(G4KineticTrackVector*, G4V3DNucleus*)
0118 {
0119   return 0;
0120 }
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123 //
0124 // ApplyYourself
0125 //
0126 // Member function to process an event, and get information about the products.
0127
0128 G4HadFinalState* G4UrQMD1_3Model::ApplyYourself(const G4HadProjectile& theTrack,
0129                                                 G4Nucleus& theTarget)
0130 {
0131   // anti_new
0132   //   ------------------define anti_light_nucleus
0133   const G4ParticleDefinition* anti_deu = G4AntiDeuteron::AntiDeuteron();
0134
0135   const G4ParticleDefinition* anti_he3 = G4AntiHe3::AntiHe3();
0136
0137   const G4ParticleDefinition* anti_tri = G4AntiTriton::AntiTriton();
0138
0139   const G4ParticleDefinition* anti_alp = G4AntiAlpha::AntiAlpha();
0140   //---------------------------------------------------
0141   //
0142   // The secondaries will be returned in G4HadFinalState &theResult -
0143   // initialise this.  The original track will always be discontinued and
0144   // secondaries followed.
0145   //
0146   theResult.Clear();
0147   theResult.SetStatusChange(stopAndKill);
0148
0149   G4DynamicParticle* cascadeParticle = 0;
0150   //
0151   //
0152   // Get relevant information about the projectile and target (A, Z, energy/nuc,
0153   // momentum, etc).
0154   //
0155
0156   const G4ParticleDefinition* definitionP = theTrack.GetDefinition();
0157   const G4double AP = definitionP->GetBaryonNumber();
0158   const G4double ZP = definitionP->GetPDGCharge();
0159   G4double AT = theTarget.GetN();
0160   G4double ZT = theTarget.GetZ();
0161   //  -----------------------------------------------
0162   G4int id = definitionP->GetPDGEncoding();  // get particle encoding
0163   // ------------------------------------------------
0164   G4int AP1 = G4lrint(AP);
0165   G4int ZP1 = G4lrint(ZP);
0166   G4int AT1 = G4lrint(AT);
0167   G4int ZT1 = G4lrint(ZT);
0168   //  G4cout<<"------ap1--=="<<AP1<<"---zp1---=="<<ZP1<<"---id-=="<<id<< G4endl;
0169   //
0170   // ****************************************************************************
0171   // The following is the parameters necessary to initiate Uinit() and UrQMD()
0172   // ----------------------------------------------------------------------------
0173   urqmdparams_.u_sptar = 0;  //! 0= normal proj/target, 1=special proj/tar
0174   urqmdparams_.u_spproj = 1;  // projectile is a special particle
0175
0176   // new_anti
0177
0178   if (AP1 > 1 || definitionP == anti_deu || definitionP == anti_he3 || definitionP == anti_tri
0179       || definitionP == anti_alp)
0180   {
0181     urqmdparams_.u_ap = AP1;
0182     urqmdparams_.u_zp = ZP1;
0183
0184     urqmdparams_.u_spproj = 0;
0185   }
0186   else if (id == 2212) {  //! proton
0187     urqmdparams_.u_ap = 1;
0188     urqmdparams_.u_zp = 1;
0189   }
0190   else if (id == -2212) {  //! anti-proton
0191     urqmdparams_.u_ap = -1;
0192     urqmdparams_.u_zp = -1;
0193   }
0194   else if (id == 2112) {  //! neutron
0195     urqmdparams_.u_ap = 1;
0196     urqmdparams_.u_zp = -1;
0197   }
0198   else if (id == -2112) {  //! anti-neutron
0199     urqmdparams_.u_ap = -1;
0200     urqmdparams_.u_zp = 1;
0201   }
0202   else if (id == 211) {  //! pi+
0203     urqmdparams_.u_ap = 101;
0204     urqmdparams_.u_zp = 2;
0205   }
0206   else if (id == -211) {  //! pi-
0207     urqmdparams_.u_ap = 101;
0208     urqmdparams_.u_zp = -2;
0209   }
0210   else if (id == 321) {  //! K+
0211     urqmdparams_.u_ap = 106;
0212     urqmdparams_.u_zp = 1;
0213   }
0214   else if (id == -321) {  //! K-
0215     urqmdparams_.u_ap = -106;
0216     urqmdparams_.u_zp = -1;
0217   }
0218   else if (id == 130 || id == 310) {  //  ! K0
0219     urqmdparams_.u_ap = 106;
0220     urqmdparams_.u_zp = -1;
0221   }
0222   else if (id == -130 || id == -310) {  // ! K0bar
0223     urqmdparams_.u_ap = -106;
0224     urqmdparams_.u_zp = 1;
0225   }
0226   else {
0227     G4cout << " Sorry, No definition for particle for UrQMD::" << id << "found" << G4endl;
0228
0229     // AND->
0230 #  if G4VERSION_NUMBER >= 950
0231     // New signature (9.5) for G4Exception
0232     // Using G4HadronicException
0233     throw G4HadronicException(__FILE__, __LINE__, "Sorry, no definition for particle for UrQMD");
0234 #  else
0235     G4Exception(" ");
0236 #  endif
0237     // AND<-
0238   }  // end if id
0239   //-------------------------------------------------------
0240
0241   urqmdparams_.u_at = AT1;  // Target identified
0242   urqmdparams_.u_zt = ZT1;
0243   //----------------------------------------------------
0244   //  identify Energy
0245   //
0246   G4ThreeVector Pbefore = theTrack.Get4Momentum().vect();
0247   G4double T = theTrack.GetKineticEnergy();
0248   G4double E = theTrack.GetTotalEnergy();
0249   G4double TotalEbefore = E * AP1 + theTarget.AtomicMass(AT1, ZT1) + theTarget.GetEnergyDeposit();
0250   //    -----------------------------------------------------------------
0251
0252   if (AP1 > 1) {
0253     urqmdparams_.u_elab = T / (AP1 * GeV);  // Units are GeV/nuc for UrQMD
0254
0255     E = E / AP1;  // Units are GeV/nuc
0256   }
0257   else {
0258     urqmdparams_.u_elab = T / GeV;  // units are GeV
0259
0260     TotalEbefore = E + theTarget.AtomicMass(AT1, ZT1) + theTarget.GetEnergyDeposit();
0261   }
0262
0263   //------------------------------------------------------------
0264   // identify impact parameter
0265   urqmdparams_.u_imp = -(1.1 * std::pow(G4double(AT1), (1. / 3.)));
0266   // units are in fm for UrQMD;
0267   //------------------------------------------------------------
0268   ///////////////////////// initialise/////////////////////
0269
0270   if (CurrentEvent == 0) {
0271     G4cout << "\n creation of table, wait-------" << G4endl;
0272
0273     G4cout << "\n" << G4endl;
0274
0275     G4int io = 0;
0276
0277     uinit_(&io);
0278
0279     G4cout << "\n end to create  table " << G4endl;
0280
0281     CurrentEvent = 1;
0282   }
0283   ////////////////////////////////////////////////////////
0284
0285   // #ifdef debug_G4UrQMD1_3Model
0286
0287   G4cout << "UrQMDModel running-------------" << G4endl;
0288
0289   urqmd_();
0290
0291   // #endif
0292
0293   // G4cout <<"Number of produced particles:  " <<sys_.npart<<G4endl;
0294
0295   G4int n = sys_.npart;  // no of produced particles
0296   if (n < 2) {
0297     G4cout << "===============Warning================" << G4endl;
0298     G4cout << "======================================" << G4endl;
0299
0300     G4cout << "Number of produced particles is very low:  " << sys_.npart << G4endl;
0301     G4cout << "============================================" << G4endl;
0302
0303 // AND->
0304 #  if G4VERSION_NUMBER >= 950
0305     // New signature (9.5) for G4Exception
0306     // Using G4HadronicException instead of base class
0307     throw G4HadronicException(__FILE__, __LINE__, "Number of produced particle is very low");
0308 #  else
0309     G4Exception(" ");  // stop
0310 #  endif
0311     // AND<-
0312   }
0313   else {
0314     for (G4int i = 0; i < n; i++) {
0315       G4int pid = pdgid_(&isys_.ityp[i], &isys_.iso3[i]);
0316
0317       // Particle is a final state secondary and not a nucleus.
0318       // Determine what this secondary particle is, and if valid, load dynamic
0319       // parameters.
0320       //
0321
0322       G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle(pid);
0323
0324       if (pd) {
0325         G4double px = (coor_.px[i] + ffermi_.ffermpx[i]) * GeV;
0326         // units are in MeV/c for G4
0327         G4double py = (coor_.py[i] + ffermi_.ffermpy[i]) * GeV;
0328         G4double pz = (coor_.pz[i] + ffermi_.ffermpz[i]) * GeV;
0329
0330         G4double et = (coor_.p0[i]) * GeV;
0331
0332         //    ------------------------------Use only "Lorentz vector"----------
0333
0334         G4LorentzVector lorenzvec = G4LorentzVector(px, py, pz, et);
0335
0336         cascadeParticle = new G4DynamicParticle(pd, lorenzvec);  //
0337
0338         theResult.AddSecondary(cascadeParticle);
0339
0340         //======================================================================
0341
0342       }  // if
0343     }  // for
0344
0345   }  // if warning
0346
0347   //=======================================================================
0348   if (verbose >= 3) {
0349     //
0350     G4double TotalEafter = 0.0;
0351     G4ThreeVector TotalPafter;
0352     G4double charge = 0.0;
0353     G4int baryon = 0;
0354     G4int nSecondaries = theResult.GetNumberOfSecondaries();
0355
0356     for (G4int j = 0; j < nSecondaries; j++) {
0357       TotalEafter += theResult.GetSecondary(j)->GetParticle()->GetTotalEnergy();
0358
0359       TotalPafter += theResult.GetSecondary(j)->GetParticle()->GetMomentum();
0360
0361       G4ParticleDefinition* pd = theResult.GetSecondary(j)->GetParticle()->GetDefinition();
0362
0363       charge += pd->GetPDGCharge();
0364       baryon += pd->GetBaryonNumber();
0365
0366     }  // for secondaries
0367
0368     G4cout << "----------------------------------------"
0369            << "----------------------------------------" << G4endl;
0370     G4cout << "Total energy before collision  = " << TotalEbefore  /// MeV
0371            << " MeV" << G4endl;
0372     G4cout << "Total energy after collision    = " << TotalEafter  // MeV
0373            << " MeV" << G4endl;
0374
0375     G4cout << "----------------------------------------" << G4endl;
0376
0377     G4cout << "Total momentum before collision = " << Pbefore  // MeV
0378            << " MeV/c" << G4endl;
0379     G4cout << "Total momentum after collision  = " << TotalPafter  // MeV
0380            << " MeV/c" << G4endl;
0381     G4cout << "----------------------------------------" << G4endl;
0382
0383     if (verbose >= 4) {
0384       G4cout << "Total charge before collision  = " << (ZP + ZT) * eplus << G4endl;
0385       G4cout << "Total charge after collision    = " << charge << G4endl;
0386
0387       G4cout << "----------------------------------------" << G4endl;
0388
0389       G4cout << "Total baryon number before collision = " << AP + AT << G4endl;
0390       G4cout << "Total baryon number after collision  = " << baryon << G4endl;
0391       G4cout << "----------------------------------------" << G4endl;
0392
0393     }  // if verbose4
0394
0395     G4cout << "----------------------------------------"
0396            << "----------------------------------------" << G4endl;
0397
0398   }  // if verbose3
0399
0400   return &theResult;
0401 }  // G4hadfinal
0402
0403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0404 //
0405 // WelcomeMessage
0406 //
0407 void G4UrQMD1_3Model::WelcomeMessage() const
0408 {
0409   G4cout << G4endl;
0410   G4cout << " *****************************************************************" << G4endl;
0411   G4cout << " Interface to        G4UrQMD_1.3                      activated" << G4endl;
0412   G4cout << " Version number : 00.00.0B          File date : 25/01/12" << G4endl;
0413   G4cout << " (Interface written by Kh. Abdel-Waged et al. for the KACST/NCMP)" << G4endl;
0414   G4cout << G4endl;
0415   G4cout << " *****************************************************************" << G4endl;
0416   G4cout << G4endl;
0417
0418   return;
0419 }
0420
0421 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0422
0423 void G4UrQMD1_3Model::InitialiseDataTables()
0424 {
0425   //
0426   //
0427   // The next line is to make sure the block data statements are
0428   // executed.
0429   //
0430
0431   g4urqmdblockdata_();
0432
0433   ///////////////////////////////////////////////////
0434   /////// Dynamic seed //////////////////////////////
0435   // G4int ranseed=-time_ ();
0436   //     Fixed seed  ///////////////////////////
0437
0438   G4int ranseed = 1097569630;
0439
0440   G4cout << "\n seed:  " << ranseed << G4endl;
0441
0442   sseed_(&ranseed);
0443
0444   loginit_();
0445 }
0446
0447 #endif  // G4_USE_URQMD