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