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