Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr02/src/G4HIJING_Model.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 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 27 //
 28 // MODULE:          G4HIJING_Model.hh
 29 //
 30 // Version:        1.B
 31 // Date:           10/09/2013
 32 // Authors:        Khaled Abdel-Waged
 33 // Institute:      Umm Al-Qura University
 34 // Country:        Saudi Arabia
 35 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 36 //
 37 #include "G4HIJING_Model.hh"
 38 #ifdef G4_USE_HIJING
 39 #  include "G4HIJING_Interface.hh"
 40 //-------------------------------
 41 #  include "G4CollisionOutput.hh"
 42 #  include "G4DynamicParticle.hh"
 43 #  include "G4IonTable.hh"
 44 #  include "G4LorentzRotation.hh"
 45 #  include "G4Nucleus.hh"
 46 #  include "G4ParticleDefinition.hh"
 47 #  include "G4ParticleTable.hh"
 48 #  include "G4Track.hh"
 49 #  include "G4V3DNucleus.hh"
 50 #  include "globals.hh"
 51 
 52 // AND->
 53 #  include "G4Version.hh"
 54 // AND<-
 55 //----------------new_anti
 56 #  include "G4AntiAlpha.hh"
 57 #  include "G4AntiDeuteron.hh"
 58 #  include "G4AntiHe3.hh"
 59 #  include "G4AntiTriton.hh"
 60 //---------------------------
 61 #  include "HistoManager.hh"  //newkhaled
 62 
 63 #  include "G4SystemOfUnits.hh"
 64 
 65 #  include <fstream>
 66 #  include <string>
 67 ///////////////////////////////////////////////////////////////////////////
 68 
 69 //
 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 71 G4HIJING_Model::G4HIJING_Model(const G4String& nam) : G4VIntraNuclearTransportModel(nam), verbose(0)
 72 {
 73   if (verbose > 3) {
 74     G4cout << " >>> G4HIJING_Model default constructor" << G4endl;
 75   }
 76 
 77 #  ifdef G4ANALYSIS_USE
 78   fHistoManager = HistoManager::GetPointer();  // new_khaled
 79 #  endif
 80 
 81   //
 82   // Set the minimum and maximum range for the HIJING model
 83 
 84   SetMinEnergy(4.0 * GeV);
 85   //  SetMaxEnergy(2000.0*TeV);
 86 
 87   //
 88 
 89   //
 90   WelcomeMessage();
 91   //
 92   CurrentEvent = 0;
 93 
 94   //
 95 
 96   InitialiseDataTables();
 97 
 98   //
 99 }
100 ////////////////////////////////////////////////////////////////////////////////
101 //
102 // Destructor
103 //
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 G4HIJING_Model::~G4HIJING_Model() {}
106 ////////////////////////////////////////////////////////////////////////////////
107 
108 G4ReactionProductVector* G4HIJING_Model::Propagate(G4KineticTrackVector*, G4V3DNucleus*)
109 {
110   return 0;
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 //
115 // ApplyYourself
116 //
117 // Member function to process an event, and get information about the products.
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 G4HadFinalState* G4HIJING_Model::ApplyYourself(const G4HadProjectile& theTrack,
121                                                G4Nucleus& theTarget)
122 {
123   G4cout << "HERE I AM" << G4endl;
124   // anti_new
125   //   ------------------define anti_light_nucleus
126   const G4ParticleDefinition* anti_deu = G4AntiDeuteron::AntiDeuteron();
127 
128   const G4ParticleDefinition* anti_he3 = G4AntiHe3::AntiHe3();
129 
130   const G4ParticleDefinition* anti_tri = G4AntiTriton::AntiTriton();
131 
132   const G4ParticleDefinition* anti_alp = G4AntiAlpha::AntiAlpha();
133 
134   //---------------------------------------------------
135   //
136   // The secondaries will be returned in G4HadFinalState &theResult -
137   // initialise this.  The original track will always be discontinued and
138   // secondaries followed.
139   //
140   theResult.Clear();
141   theResult.SetStatusChange(stopAndKill);
142 
143   G4DynamicParticle* cascadeParticle = 0;
144   //
145   //
146   // Get relevant information about the projectile and target (A, Z, energy/nuc,
147   // momentum, etc).
148   //
149 
150   const G4ParticleDefinition* definitionP = theTrack.GetDefinition();
151   const G4double AP = definitionP->GetBaryonNumber();
152   const G4double ZP = definitionP->GetPDGCharge();
153   G4int AT = theTarget.GetN_asInt();
154   G4int ZT = theTarget.GetZ_asInt();
155   //  -----------------------------------------------
156   G4int id = definitionP->GetPDGEncoding();  // get particle encoding
157 
158   //      G4cout<<"particle id=========       "<<id<<G4endl;
159   // ------------------------------------------------
160   G4int AP1 = G4lrint(AP);
161   G4int ZP1 = G4lrint(ZP);
162   G4int AT1 = AT;
163   G4int ZT1 = ZT;
164 
165   // ****************************************************************************
166   // The following is the parameters necessary to initiate HIJSET() and HIJING()
167   // ----------------------------------------------------------------------------
168   //           hiparnt_.ihpr2[3]=0;     //switch off(=0) /  on(=1) jet quenching
169   //           hiparnt_.ihpr2[2]=1;     //switch on triggered Jet production
170   // ---------------------------------------------------------------------------
171   //        hiparnt_.ihnt2[0]=AP1;  //Projectile
172   hiparnt_.ihnt2[1] = ZP1;
173   hiparnt_.ihnt2[2] = AT1;  // Target
174   hiparnt_.ihnt2[3] = ZT1;
175   hiparnt_.ihnt2[5] = 0;  // Special Target
176 
177   if (AP1 > 1 || definitionP == anti_deu || definitionP == anti_he3 || definitionP == anti_tri
178       || definitionP == anti_alp)
179 
180   {
181     hiparnt_.ihnt2[0] = AP1;
182     hiparnt_.ihnt2[4] = 0;  // Special Projectile
183   }
184   else if (id == 2212) {  //! proton
185 
186     hiparnt_.ihnt2[0] = 1;
187     hiparnt_.ihnt2[4] = 2212;
188   }
189   else if (id == -2212) {  //! anti-proton
190 
191     hiparnt_.ihnt2[0] = 1;
192     hiparnt_.ihnt2[4] = -2212;
193   }
194   else if (id == 2112) {  //! neutron
195 
196     hiparnt_.ihnt2[0] = 1;
197     hiparnt_.ihnt2[4] = 2112;
198   }
199   else if (id == -2112) {  //! anti-neutron
200 
201     hiparnt_.ihnt2[0] = 1;
202     hiparnt_.ihnt2[4] = -2112;
203   }
204   else if (id == 211) {  //! pi+
205     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
206     hiparnt_.ihnt2[4] = 211;
207   }
208   else if (id == -211) {  //! pi-
209 
210     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
211     hiparnt_.ihnt2[4] = -211;
212   }
213   else if (id == 321) {  //! K+
214 
215     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
216     hiparnt_.ihnt2[4] = 321;
217   }
218   else if (id == -321) {  //! K-
219 
220     hiparnt_.ihnt2[0] = 1;  // needed by HIJING
221     hiparnt_.ihnt2[4] = -321;
222   }
223   else {
224     G4cout << " Sorry, No definition for PROJECTLE for HIJING::" << id << "found" << G4endl;
225 
226     // AND->
227 #  if G4VERSION_NUMBER >= 950
228     // New signature (9.5) for G4Exception
229     // Using G4HadronicException
230     throw G4HadronicException(__FILE__, __LINE__, "Sorry, no definition for PROJECTILE for HIJING");
231 #  else
232     G4Exception(" ");
233 #  endif
234     // AND<-
235   }  // end if id
236 
237   //-------------------------------------------------------
238   //  -------------identify mass -------------------------
239 
240   G4int id_n = 2112;
241   G4int id_p = 2212;
242 
243   hiparnt_.hint1[7] = std::max(ulmass_(&id_n), ulmass_(&id_p));
244 
245   hiparnt_.hint1[8] = hiparnt_.hint1[7];
246 
247   if (hiparnt_.ihnt2[4] != 0) hiparnt_.hint1[7] = ulmass_(&hiparnt_.ihnt2[4]);
248   // rest mass of the projectile HIJING
249 
250   //----------------------------------------------------
251   //  identify Energy
252   //
253 
254   G4double m = hiparnt_.hint1[7];  // mass in GeV
255 
256   G4ThreeVector P3 = theTrack.Get4Momentum().vect() / GeV;
257   // momentum in GeV
258 
259   G4double Pbeam = P3.z();
260   // momentum in z-direction
261 
262   G4double Ebeam = Eplab(m, Pbeam);
263   // calculate Energy of beam
264 
265   // G4cout<<"mass= "<<m<<"  P3= "<<P3<<endl;
266 
267   //---------------------------Beam ---------------------------------------
268 
269   // Lab frame: beam moves in negative z-direction
270 
271   G4LorentzVector lab = G4LorentzVector(0.0, 0.0, -1.0 * Pbeam, Ebeam + m);
272 
273   G4double TotalPbefore = -1.0 * lab.z();
274   // Calculate z-Momentum before collision
275   //
276   G4double TotalEbefore = lab.e();
277   // Calculate Energy before collision
278 
279   //   --------------------------------------------------------
280   //                     Turn to CM frame:
281   //   ---------------------------------------------------------
282 
283   G4LorentzVector cms = G4LorentzVector(0.0, 0.0, 0.0, lab.mag());
284 
285   // ----------------------Get relative speed between frames---------
286   // ----------------------------------------------------------------
287   G4LorentzVector Psum = (lab + cms);  // 4-Momentum sum
288   G4double beta_rel = Psum.beta();
289 
290   //---------------------Transform to equal frame--------------------
291   //-----------------------------------------------------------------
292 
293   Psum.boost(0.0, 0.0, -1.0 * beta_rel);
294 
295   //-----------------Get equal speed velocity between frames--------
296   G4double betann = Psum.beta();
297   // G4double gama= Psum.gamma();
298 
299   // ----------Colliding CM Energy per nucleon-nucleon for HIJING-
300   // ----------------------------------------------------
301 
302   G4double Ecms = lab.mag();  // CM energy for HIJING
303   efrm = Ecms;  // units are in GeV for HIJING
304 
305   ///////////////////////// initialise/////////////////////
306 
307   if (CurrentEvent == 0) {
308     G4cout << "\n initialise HIJING, wait-------" << G4endl;
309 
310     G4cout << "\n" << G4endl;
311 
312     // hijset_ (&efrm,&AP1,&ZP1,&AT1,&ZT1);
313 
314     hijset_(&efrm);
315 
316     G4cout << "\n end initialize " << G4endl;
317 
318     CurrentEvent = 1;
319   }
320   ////////////////////////////////////////////////////////
321   //------------------------------------------------------------
322   // identify impact parameter
323   bmin = 0.0;
324   //   bmax=0.5;
325 
326   bmax = hiparnt_.hipr1[33] + hiparnt_.hipr1[34];
327 
328   //----------------------------------------------
329 
330   do {
331     G4cout << "HIJING_Model running-------------" << G4endl;
332 
333     hijing_(&bmin, &bmax);
334 
335     Nproduce = himain1_.natt;  // no of produced particles
336 
337     if (Nproduce < 2) {
338       G4cout << "===============Warning=====================================" << G4endl;
339       G4cout << "-----------------------------------------------------------" << G4endl;
340       G4cout << "Number of produced particles is very low:  " << himain1_.natt << G4endl;
341       G4cout << "------------------------------------------------------------" << G4endl;
342       G4cout << "============================================================" << G4endl;
343     }
344   } while (Nproduce < 2);
345   // =============================================================================
346 
347   G4double BB = hiparnt_.hint1[18];  // impact parameter HINT1(19)
348   //     cout<<"HIJING=====impact parameter= "<<BB<<endl;
349 
350   for (G4int i = 0; i < Nproduce; i++) {
351     G4int pid = himain2_.katt[0][i];
352 
353     // Particle is a final state secondary and not a nucleus.
354     // Determine what this secondary particle is, and if valid, load dynamic
355     // parameters.
356     //
357     //   G4cout<<"pid================"<<pid<<G4endl;
358 
359     G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->FindParticle(pid);
360     ///////////////////////////////////////////////////////////////
361     //  exclude beam nucleons as produced particles
362     //     cout<<" himain2_.katt[1][i]== "<<himain2_.katt[1][i]<<endl;
363     //    if(himain2_.katt[1][i]==0 || himain2_.katt[1][i]==10) continue;
364     //  -----------------------------------------------------------
365     //      --------------reject neutral particles by calling luchge <new>
366     //         G4int chg_HIJ=luchge_ (&pid);
367     //        if (chg_HIJ==0) continue;
368 
369     if (pd) {
370       // units are in MeV/c for G4
371 
372       G4double px = himain2_.patt[0][i] * GeV;
373       G4double py = himain2_.patt[1][i] * GeV;
374       G4double pz = himain2_.patt[2][i] * GeV;
375       G4double et = himain2_.patt[3][i] * GeV;
376 
377       //    ------------------------------Use  "Lorentz vector"----------
378       G4LorentzVector lorenzCM = G4LorentzVector(px, py, pz, et);
379       //    Move to the lab frame
380       lorenzCM.boost(0.0, 0.0, -1.0 * betann);
381       G4LorentzVector lorenzLab =
382         G4LorentzVector(lorenzCM.px(), lorenzCM.py(), -1.0 * lorenzCM.pz(), lorenzCM.e());
383       //-------------------------------------------------------------------
384       cascadeParticle = new G4DynamicParticle(pd, lorenzLab);
385 
386       theResult.AddSecondary(cascadeParticle);
387 
388     }  // if pd
389 
390   }  // for
391 
392 #  ifdef G4ANALYSIS_USE  // khaled new
393   fHistoManager->StoreSecondaries(BB, theResult);
394 #  endif
395   //} //if warning
396 
397   //
398 
399   //=======================================================================
400   if (verbose >= 3) {
401     //
402     G4double TotalEafter = 0.0;
403     G4ThreeVector TotalPafter;
404     G4double charge = 0.0;
405     G4int baryon = 0;
406     G4int nSecondaries = theResult.GetNumberOfSecondaries();
407 
408     for (G4int j = 0; j < nSecondaries; j++) {
409       TotalEafter += theResult.GetSecondary(j)->GetParticle()->GetTotalEnergy() / GeV;
410 
411       TotalPafter += theResult.GetSecondary(j)->GetParticle()->GetMomentum() / GeV;
412 
413       G4ParticleDefinition* pd = theResult.GetSecondary(j)->GetParticle()->GetDefinition();
414 
415       charge += pd->GetPDGCharge();
416       baryon += pd->GetBaryonNumber();
417 
418     }  // for secondaries
419 
420     G4cout << "----------------------------------------"
421            << "----------------------------------------" << G4endl;
422     G4cout << "Total energy before collision  = " << TotalEbefore  /// GeV
423            << " GeV" << G4endl;
424     G4cout << "Total energy after collision    = "
425            << TotalEafter / nSecondaries
426            // GeV
427            << " GeV" << G4endl;
428 
429     G4cout << "----------------------------------------" << G4endl;
430 
431     G4cout << "Total momentum before collision = "
432            << TotalPbefore
433            // GeV
434            << " GeV/c" << G4endl;
435     G4cout << "Total momentum after collision  = "
436            << TotalPafter.z() / nSecondaries
437            // GeV
438            << " GeV/c" << G4endl;
439     G4cout << "----------------------------------------" << G4endl;
440 
441     if (verbose >= 4) {
442       G4cout << "Total charge before collision  = " << (ZP + ZT)  //
443              << G4endl;
444       G4cout << "Total charge after collision    = " << charge << G4endl;
445 
446       G4cout << "----------------------------------------" << G4endl;
447 
448       G4cout << "Total baryon number before collision = " << AP + AT << G4endl;
449       G4cout << "Total baryon number after collision  = " << baryon << G4endl;
450       G4cout << "----------------------------------------" << G4endl;
451 
452     }  // if verbose4
453 
454     G4cout << "----------------------------------------"
455            << "----------------------------------------" << G4endl;
456 
457   }  // if verbose3
458 
459   return &theResult;
460 }  // G4hadfinal
461 
462 //---------------------------------------------------------------------
463 
464 //---------------------------------------------------------------------
465 //
466 // WelcomeMessage
467 //
468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
469 void G4HIJING_Model::WelcomeMessage() const
470 {
471   G4cout << G4endl;
472   G4cout << " *****************************************************************" << G4endl;
473   G4cout << " Interface to        G4HIJING_Model                      activated" << G4endl;
474   G4cout << " Version number : 01.00.0B          File date : 10/09/2013" << G4endl;
475   G4cout << "  Interface written by    Khaled Abdel-Waged              " << G4endl;
476   G4cout << "                       Umm Al-Qura University             " << G4endl;
477   G4cout << "                         SAUDI ARABIA                     " << G4endl;
478   G4cout << G4endl;
479   G4cout << " *****************************************************************" << G4endl;
480   G4cout << G4endl;
481   return;
482 }
483 
484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
485 void G4HIJING_Model::InitialiseDataTables()
486 {
487   //
488   //
489   // The next line is to make sure the block data statements are
490   // executed.
491   //
492 
493   g4hijingblockdata_();
494 }
495 
496 G4double G4HIJING_Model::Eplab(G4double m, G4double P)
497 {
498   G4double Eb = std::sqrt(P * P + m * m);
499   return Eb;
500 }
501 #endif
502