Hadr02/src/G4HIJING_Model.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 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 //
0028 // MODULE:          G4HIJING_Model.hh
0029 //
0030 // Version:        1.B
0031 // Date:           10/09/2013
0032 // Authors:        Khaled Abdel-Waged
0033 // Institute:      Umm Al-Qura University
0034 // Country:        Saudi Arabia
0035 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 //
0037 #include "G4HIJING_Model.hh"
0038 #ifdef G4_USE_HIJING
0039 #  include "G4HIJING_Interface.hh"
0040 //-------------------------------
0041 #  include "G4CollisionOutput.hh"
0042 #  include "G4DynamicParticle.hh"
0043 #  include "G4IonTable.hh"
0044 #  include "G4LorentzRotation.hh"
0045 #  include "G4Nucleus.hh"
0046 #  include "G4ParticleDefinition.hh"
0047 #  include "G4ParticleTable.hh"
0048 #  include "G4Track.hh"
0049 #  include "G4V3DNucleus.hh"
0050 #  include "globals.hh"
0051
0052 // AND->
0053 #  include "G4Version.hh"
0054 // AND<-
0055 //----------------new_anti
0056 #  include "G4AntiAlpha.hh"
0057 #  include "G4AntiDeuteron.hh"
0058 #  include "G4AntiHe3.hh"
0059 #  include "G4AntiTriton.hh"
0060 //---------------------------
0061 #  include "HistoManager.hh"  //newkhaled
0062
0063 #  include "G4SystemOfUnits.hh"
0064
0065 #  include <fstream>
0066 #  include <string>
0067 ///////////////////////////////////////////////////////////////////////////
0068
0069 //
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071 G4HIJING_Model::G4HIJING_Model(const G4String& nam) : G4VIntraNuclearTransportModel(nam), verbose(0)
0072 {
0073   if (verbose > 3) {
0074     G4cout << " >>> G4HIJING_Model default constructor" << G4endl;
0075   }
0076
0077 #  ifdef G4ANALYSIS_USE
0078   fHistoManager = HistoManager::GetPointer();  // new_khaled
0079 #  endif
0080
0081   //
0082   // Set the minimum and maximum range for the HIJING model
0083
0084   SetMinEnergy(4.0 * GeV);
0085   //  SetMaxEnergy(2000.0*TeV);
0086
0087   //
0088
0089   //
0090   WelcomeMessage();
0091   //
0092   CurrentEvent = 0;
0093
0094   //
0095
0096   InitialiseDataTables();
0097
0098   //
0099 }
0100 ////////////////////////////////////////////////////////////////////////////////
0101 //
0102 // Destructor
0103 //
0104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0105 G4HIJING_Model::~G4HIJING_Model() {}
0106 ////////////////////////////////////////////////////////////////////////////////
0107
0108 G4ReactionProductVector* G4HIJING_Model::Propagate(G4KineticTrackVector*, G4V3DNucleus*)
0109 {
0110   return 0;
0111 }
0112
0113 ////////////////////////////////////////////////////////////////////////////////
0114 //
0115 // ApplyYourself
0116 //
0117 // Member function to process an event, and get information about the products.
0118
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120 G4HadFinalState* G4HIJING_Model::ApplyYourself(const G4HadProjectile& theTrack,
0121                                                G4Nucleus& theTarget)
0122 {
0123   G4cout << "HERE I AM" << G4endl;
0124   // anti_new
0125   //   ------------------define anti_light_nucleus
0126   const G4ParticleDefinition* anti_deu = G4AntiDeuteron::AntiDeuteron();
0127
0128   const G4ParticleDefinition* anti_he3 = G4AntiHe3::AntiHe3();
0129
0130   const G4ParticleDefinition* anti_tri = G4AntiTriton::AntiTriton();
0131
0132   const G4ParticleDefinition* anti_alp = G4AntiAlpha::AntiAlpha();
0133
0134   //---------------------------------------------------
0135   //
0136   // The secondaries will be returned in G4HadFinalState &theResult -
0137   // initialise this.  The original track will always be discontinued and
0138   // secondaries followed.
0139   //
0140   theResult.Clear();
0141   theResult.SetStatusChange(stopAndKill);
0142
0143   G4DynamicParticle* cascadeParticle = 0;
0144   //
0145   //
0146   // Get relevant information about the projectile and target (A, Z, energy/nuc,
0147   // momentum, etc).
0148   //
0149
0150   const G4ParticleDefinition* definitionP = theTrack.GetDefinition();
0151   const G4double AP = definitionP->GetBaryonNumber();
0152   const G4double ZP = definitionP->GetPDGCharge();
0153   G4int AT = theTarget.GetN_asInt();
0154   G4int ZT = theTarget.GetZ_asInt();
0155   //  -----------------------------------------------
0156   G4int id = definitionP->GetPDGEncoding();  // get particle encoding
0157
0158   //      G4cout<<"particle id=========       "<<id<<G4endl;
0159   // ------------------------------------------------
0160   G4int AP1 = G4lrint(AP);
0161   G4int ZP1 = G4lrint(ZP);
0162   G4int AT1 = AT;
0163   G4int ZT1 = ZT;
0164
0165   // ****************************************************************************
0166   // The following is the parameters necessary to initiate HIJSET() and HIJING()
0167   // ----------------------------------------------------------------------------
0168   //           hiparnt_.ihpr2[3]=0;     //switch off(=0) /  on(=1) jet quenching
0169   //           hiparnt_.ihpr2[2]=1;     //switch on triggered Jet production
0170   // ---------------------------------------------------------------------------
0171   //        hiparnt_.ihnt2[0]=AP1;  //Projectile
0172   hiparnt_.ihnt2[1] = ZP1;
0173   hiparnt_.ihnt2[2] = AT1;  // Target
0174   hiparnt_.ihnt2[3] = ZT1;
0175   hiparnt_.ihnt2[5] = 0;  // Special Target
0176
0177   if (AP1 > 1 || definitionP == anti_deu || definitionP == anti_he3 || definitionP == anti_tri
0178       || definitionP == anti_alp)
0179
0180   {
0181     hiparnt_.ihnt2[0] = AP1;
0182     hiparnt_.ihnt2[4] = 0;  // Special Projectile
0183   }
0184   else if (id == 2212) {  //! proton
0185
0186     hiparnt_.ihnt2[0] = 1;
0187     hiparnt_.ihnt2[4] = 2212;
0188   }
0189   else if (id == -2212) {  //! anti-proton
0190
0191     hiparnt_.ihnt2[0] = 1;
0192     hiparnt_.ihnt2[4] = -2212;
0193   }
0194   else if (id == 2112) {  //! neutron
0195
0196     hiparnt_.ihnt2[0] = 1;
0197     hiparnt_.ihnt2[4] = 2112;
0198   }
0199   else if (id == -2112) {  //! anti-neutron
0200
0201     hiparnt_.ihnt2[0] = 1;
0202     hiparnt_.ihnt2[4] = -2112;
0203   }
0204   else if (id == 211) {  //! pi+
0205     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
0206     hiparnt_.ihnt2[4] = 211;
0207   }
0208   else if (id == -211) {  //! pi-
0209
0210     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
0211     hiparnt_.ihnt2[4] = -211;
0212   }
0213   else if (id == 321) {  //! K+
0214
0215     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
0216     hiparnt_.ihnt2[4] = 321;
0217   }
0218   else if (id == -321) {  //! K-
0219
0220     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
0221     hiparnt_.ihnt2[4] = -321;
0222   }
0223   else {
0224     G4cout << " Sorry, No definition for PROJECTLE for HIJING::" << id << "found" << G4endl;
0225
0226     // AND->
0227 #  if G4VERSION_NUMBER >= 950
0228     // New signature (9.5) for G4Exception
0229     // Using G4HadronicException
0230     throw G4HadronicException(__FILE__, __LINE__, "Sorry, no definition for PROJECTILE for HIJING");
0231 #  else
0232     G4Exception(" ");
0233 #  endif
0234     // AND<-
0235   }  // end if id
0236
0237   //-------------------------------------------------------
0238   //  -------------identify mass -------------------------
0239
0240   G4int id_n = 2112;
0241   G4int id_p = 2212;
0242
0243   hiparnt_.hint1[7] = std::max(ulmass_(&id_n), ulmass_(&id_p));
0244
0245   hiparnt_.hint1[8] = hiparnt_.hint1[7];
0246
0247   if (hiparnt_.ihnt2[4] != 0) hiparnt_.hint1[7] = ulmass_(&hiparnt_.ihnt2[4]);
0248   // rest mass of the projectile HIJING
0249
0250   //----------------------------------------------------
0251   //  identify Energy
0252   //
0253
0254   G4double m = hiparnt_.hint1[7];  // mass in GeV
0255
0256   G4ThreeVector P3 = theTrack.Get4Momentum().vect() / GeV;
0257   // momentum in GeV
0258
0259   G4double Pbeam = P3.z();
0260   // momentum in z-direction
0261
0262   G4double Ebeam = Eplab(m, Pbeam);
0263   // calculate Energy of beam
0264
0265   // G4cout<<"mass= "<<m<<"  P3= "<<P3<<endl;
0266
0267   //---------------------------Beam ---------------------------------------
0268
0269   // Lab frame: beam moves in negative z-direction
0270
0271   G4LorentzVector lab = G4LorentzVector(0.0, 0.0, -1.0 * Pbeam, Ebeam + m);
0272
0273   G4double TotalPbefore = -1.0 * lab.z();
0274   // Calculate z-Momentum before collision
0275   //
0276   G4double TotalEbefore = lab.e();
0277   // Calculate Energy before collision
0278
0279   //   --------------------------------------------------------
0280   //                     Turn to CM frame:
0281   //   ---------------------------------------------------------
0282
0283   G4LorentzVector cms = G4LorentzVector(0.0, 0.0, 0.0, lab.mag());
0284
0285   // ----------------------Get relative speed between frames---------
0286   // ----------------------------------------------------------------
0287   G4LorentzVector Psum = (lab + cms);  // 4-Momentum sum
0288   G4double beta_rel = Psum.beta();
0289
0290   //---------------------Transform to equal frame--------------------
0291   //-----------------------------------------------------------------
0292
0293   Psum.boost(0.0, 0.0, -1.0 * beta_rel);
0294
0295   //-----------------Get equal speed velocity between frames--------
0296   G4double betann = Psum.beta();
0297   // G4double gama= Psum.gamma();
0298
0299   // ----------Colliding CM Energy per nucleon-nucleon for HIJING-
0300   // ----------------------------------------------------
0301
0302   G4double Ecms = lab.mag();  // CM energy for HIJING
0303   efrm = Ecms;  // units are in GeV for HIJING
0304
0305   ///////////////////////// initialise/////////////////////
0306
0307   if (CurrentEvent == 0) {
0308     G4cout << "\n initialise HIJING, wait-------" << G4endl;
0309
0310     G4cout << "\n" << G4endl;
0311
0312     // hijset_ (&efrm,&AP1,&ZP1,&AT1,&ZT1);
0313
0314     hijset_(&efrm);
0315
0316     G4cout << "\n end initialize " << G4endl;
0317
0318     CurrentEvent = 1;
0319   }
0320   ////////////////////////////////////////////////////////
0321   //------------------------------------------------------------
0322   // identify impact parameter
0323   bmin = 0.0;
0324   //   bmax=0.5;
0325
0326   bmax = hiparnt_.hipr1[33] + hiparnt_.hipr1[34];
0327
0328   //----------------------------------------------
0329
0330   do {
0331     G4cout << "HIJING_Model running-------------" << G4endl;
0332
0333     hijing_(&bmin, &bmax);
0334
0335     Nproduce = himain1_.natt;  // no of produced particles
0336
0337     if (Nproduce < 2) {
0338       G4cout << "===============Warning=====================================" << G4endl;
0339       G4cout << "-----------------------------------------------------------" << G4endl;
0340       G4cout << "Number of produced particles is very low:  " << himain1_.natt << G4endl;
0341       G4cout << "------------------------------------------------------------" << G4endl;
0342       G4cout << "============================================================" << G4endl;
0343     }
0344   } while (Nproduce < 2);
0345   // =============================================================================
0346
0347   G4double BB = hiparnt_.hint1[18];  // impact parameter HINT1(19)
0348   //     cout<<"HIJING=====impact parameter= "<<BB<<endl;
0349
0350   for (G4int i = 0; i < Nproduce; i++) {
0351     G4int pid = himain2_.katt[0][i];
0352
0353     // Particle is a final state secondary and not a nucleus.
0354     // Determine what this secondary particle is, and if valid, load dynamic
0355     // parameters.
0356     //
0357     //   G4cout<<"pid================"<<pid<<G4endl;
0358
0359     G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle(pid);
0360     ///////////////////////////////////////////////////////////////
0361     //  exclude beam nucleons as produced particles
0362     //     cout<<" himain2_.katt[1][i]== "<<himain2_.katt[1][i]<<endl;
0363     //    if(himain2_.katt[1][i]==0 || himain2_.katt[1][i]==10) continue;
0364     //  -----------------------------------------------------------
0365     //      --------------reject neutral particles by calling luchge <new>
0366     //         G4int chg_HIJ=luchge_ (&pid);
0367     //        if (chg_HIJ==0) continue;
0368
0369     if (pd) {
0370       // units are in MeV/c for G4
0371
0372       G4double px = himain2_.patt[0][i] * GeV;
0373       G4double py = himain2_.patt[1][i] * GeV;
0374       G4double pz = himain2_.patt[2][i] * GeV;
0375       G4double et = himain2_.patt[3][i] * GeV;
0376
0377       //    ------------------------------Use  "Lorentz vector"----------
0378       G4LorentzVector lorenzCM = G4LorentzVector(px, py, pz, et);
0379       //    Move to the lab frame
0380       lorenzCM.boost(0.0, 0.0, -1.0 * betann);
0381       G4LorentzVector lorenzLab =
0382         G4LorentzVector(lorenzCM.px(), lorenzCM.py(), -1.0 * lorenzCM.pz(), lorenzCM.e());
0383       //-------------------------------------------------------------------
0384       cascadeParticle = new G4DynamicParticle(pd, lorenzLab);
0385
0386       theResult.AddSecondary(cascadeParticle);
0387
0388     }  // if pd
0389
0390   }  // for
0391
0392 #  ifdef G4ANALYSIS_USE  // khaled new
0393   fHistoManager->StoreSecondaries(BB, theResult);
0394 #  endif
0395   //} //if warning
0396
0397   //
0398
0399   //=======================================================================
0400   if (verbose >= 3) {
0401     //
0402     G4double TotalEafter = 0.0;
0403     G4ThreeVector TotalPafter;
0404     G4double charge = 0.0;
0405     G4int baryon = 0;
0406     G4int nSecondaries = theResult.GetNumberOfSecondaries();
0407
0408     for (G4int j = 0; j < nSecondaries; j++) {
0409       TotalEafter += theResult.GetSecondary(j)->GetParticle()->GetTotalEnergy() / GeV;
0410
0411       TotalPafter += theResult.GetSecondary(j)->GetParticle()->GetMomentum() / GeV;
0412
0413       G4ParticleDefinition* pd = theResult.GetSecondary(j)->GetParticle()->GetDefinition();
0414
0415       charge += pd->GetPDGCharge();
0416       baryon += pd->GetBaryonNumber();
0417
0418     }  // for secondaries
0419
0420     G4cout << "----------------------------------------"
0421            << "----------------------------------------" << G4endl;
0422     G4cout << "Total energy before collision  = " << TotalEbefore  /// GeV
0423            << " GeV" << G4endl;
0424     G4cout << "Total energy after collision    = "
0425            << TotalEafter / nSecondaries
0426            // GeV
0427            << " GeV" << G4endl;
0428
0429     G4cout << "----------------------------------------" << G4endl;
0430
0431     G4cout << "Total momentum before collision = "
0432            << TotalPbefore
0433            // GeV
0434            << " GeV/c" << G4endl;
0435     G4cout << "Total momentum after collision  = "
0436            << TotalPafter.z() / nSecondaries
0437            // GeV
0438            << " GeV/c" << G4endl;
0439     G4cout << "----------------------------------------" << G4endl;
0440
0441     if (verbose >= 4) {
0442       G4cout << "Total charge before collision  = " << (ZP + ZT)  //
0443              << G4endl;
0444       G4cout << "Total charge after collision    = " << charge << G4endl;
0445
0446       G4cout << "----------------------------------------" << G4endl;
0447
0448       G4cout << "Total baryon number before collision = " << AP + AT << G4endl;
0449       G4cout << "Total baryon number after collision  = " << baryon << G4endl;
0450       G4cout << "----------------------------------------" << G4endl;
0451
0452     }  // if verbose4
0453
0454     G4cout << "----------------------------------------"
0455            << "----------------------------------------" << G4endl;
0456
0457   }  // if verbose3
0458
0459   return &theResult;
0460 }  // G4hadfinal
0461
0462 //---------------------------------------------------------------------
0463
0464 //---------------------------------------------------------------------
0465 //
0466 // WelcomeMessage
0467 //
0468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0469 void G4HIJING_Model::WelcomeMessage() const
0470 {
0471   G4cout << G4endl;
0472   G4cout << " *****************************************************************" << G4endl;
0473   G4cout << " Interface to        G4HIJING_Model                      activated" << G4endl;
0474   G4cout << " Version number : 01.00.0B          File date : 10/09/2013" << G4endl;
0475   G4cout << "  Interface written by    Khaled Abdel-Waged              " << G4endl;
0476   G4cout << "                       Umm Al-Qura University             " << G4endl;
0477   G4cout << "                         SAUDI ARABIA                     " << G4endl;
0478   G4cout << G4endl;
0479   G4cout << " *****************************************************************" << G4endl;
0480   G4cout << G4endl;
0481   return;
0482 }
0483
0484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0485 void G4HIJING_Model::InitialiseDataTables()
0486 {
0487   //
0488   //
0489   // The next line is to make sure the block data statements are
0490   // executed.
0491   //
0492
0493   g4hijingblockdata_();
0494 }
0495
0496 G4double G4HIJING_Model::Eplab(G4double m, G4double P)
0497 {
0498   G4double Eb = std::sqrt(P * P + m * m);
0499   return Eb;
0500 }
0501 #endif