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