ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VPartonStringModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VPartonStringModel.cc
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 //
29 // GEANT 4 class implementation file
30 //
31 // ---------------- G4VPartonStringModel ----------------
32 // by Gunter Folger, May 1998.
33 // abstract class for all Parton String Models
34 // ------------------------------------------------------------
35 // debug switch
36 //#define debug_PartonStringModel
37 // ------------------------------------------------------------
38 
39 #include "G4VPartonStringModel.hh"
40 #include "G4ios.hh"
42 
43 #include "G4ParticleTable.hh"
44 #include "G4IonTable.hh"
45 
47  : G4VHighEnergyGenerator(modelName),
48  stringFragmentationModel(0),
49  theThis(0)
50 {
51 // Make shure Shotrylived particles are constructed.
52  G4ShortLivedConstructor ShortLived;
53  ShortLived.ConstructParticle();
54 }
55 
57 {
58 }
59 
61  const G4DynamicParticle &aPrimary)
62 {
63  G4ExcitedStringVector * strings = NULL;
64  G4DynamicParticle thePrimary=aPrimary;
65  G4LorentzVector SumStringMom(0.,0.,0.,0.);
66  G4KineticTrackVector * theResult = 0;
67  G4Nucleon * theNuclNucleon(0);
68 
69  #ifdef debug_PartonStringModel
70  G4cout<<G4endl;
71  G4cout<<"-----------------------Parton-String model is runnung ------------"<<G4endl;
72  G4cout<<"Projectile Name Mass "<<thePrimary.GetDefinition()->GetParticleName()<<" "
73  <<thePrimary.GetMass()<<G4endl;
74  G4cout<<" Momentum "<<thePrimary.Get4Momentum()<<G4endl;
75  G4cout<<"Target nucleus A Z "<<theNucleus.GetA_asInt()<<" "
76  <<theNucleus.GetZ_asInt()<<G4endl<<G4endl;
77  G4int Bsum=thePrimary.GetDefinition()->GetBaryonNumber() + theNucleus.GetA_asInt();
78  G4int Qsum=thePrimary.GetDefinition()->GetPDGCharge() + theNucleus.GetZ_asInt();
79  G4cout<<"Initial baryon number "<<Bsum<<G4endl;
80  G4cout<<"Initial charge "<<Qsum<<G4endl;
81  G4cout<<"-------------- Parton-String model: Generation of strings -------"<<G4endl<<G4endl;
82  Bsum -= theNucleus.GetA_asInt(); Qsum -= theNucleus.GetZ_asInt();
84  Bsum -= thePrimary.GetDefinition()->GetBaryonNumber();
85  Qsum -= thePrimary.GetDefinition()->GetPDGCharge();
86  }
87  G4int QsumSec(0), BsumSec(0);
88  G4LorentzVector SumPsecondr(0.,0.,0.,0.);
89  #endif
90 
92  G4LorentzVector Ptmp=thePrimary.Get4Momentum();
93  toZ.rotateZ(-1*Ptmp.phi());
94  toZ.rotateY(-1*Ptmp.theta());
95  thePrimary.Set4Momentum(toZ*Ptmp);
96  G4LorentzRotation toLab(toZ.inverse());
97 
98  G4bool Success=true;
99  G4int attempts = 0, maxAttempts=1000;
100  do
101  {
102  if (attempts++ > maxAttempts )
103  {
104  theThis->Init(theNucleus,thePrimary); // To put a nucleus into ground state
105  // But marks of hitted nucleons are left. They must be erased.
106  G4V3DNucleus * ResNucleus=theThis->GetWoundedNucleus();
107  theNuclNucleon = ResNucleus->StartLoop() ? ResNucleus->GetNextNucleon() : NULL;
108  while( theNuclNucleon )
109  {
110  if(theNuclNucleon->AreYouHit()) theNuclNucleon->Hit(nullptr);
111  theNuclNucleon = ResNucleus->GetNextNucleon();
112  }
113 
114  G4V3DNucleus * ProjResNucleus=theThis->GetProjectileNucleus();
115  if(ProjResNucleus != 0)
116  {
117  theNuclNucleon = ProjResNucleus->StartLoop() ? ProjResNucleus->GetNextNucleon() : NULL;
118  while( theNuclNucleon )
119  {
120  if(theNuclNucleon->AreYouHit()) theNuclNucleon->Hit(nullptr);
121  theNuclNucleon = ProjResNucleus->GetNextNucleon();
122  }
123  }
124 
126  ed << "Projectile Name Mass " <<thePrimary.GetDefinition()->GetParticleName()
127  << " " << thePrimary.GetMass()<< G4endl;
128  ed << "           Momentum  " << thePrimary.Get4Momentum() <<G4endl;
129  ed << "Target nucleus   A Z " << theNucleus.GetA_asInt() << " "
130  << theNucleus.GetZ_asInt() <<G4endl;
131  ed << "Initial states of projectile and target nucleus will be returned!"<<G4endl;
132  G4Exception( "G4VPartonStringModel::Scatter(): fails to generate or fragment strings ",
133  "HAD_PARTON_STRING_001", JustWarning, ed );
134 
135  G4ThreeVector Position(0.,0.,2*ResNucleus->GetOuterRadius());
136  G4KineticTrack* Hadron = new G4KineticTrack(aPrimary.GetParticleDefinition(), 0.,
137  Position, aPrimary.Get4Momentum());
138  if(theResult == 0) theResult = new G4KineticTrackVector();
139  theResult->push_back(Hadron);
140  return theResult;
141  }
142 
143  Success=true;
144 
145  theThis->Init(theNucleus,thePrimary);
146 
147  strings = GetStrings();
148 
149  if (strings->size() == 0) { Success=false; continue; }
150 
151  G4double stringEnergy(0);
152  SumStringMom=G4LorentzVector(0.,0.,0.,0.);
153 
154  #ifdef debug_PartonStringModel
155  G4cout<<"------------ Parton-String model: Number of produced strings ---- "<<strings->size()<<G4endl;
156  #endif
157 
158  for ( unsigned int astring=0; astring < strings->size(); astring++)
159  {
160  // rotate string to lab frame, models have it aligned to z
161  if((*strings)[astring]->IsExcited())
162  {
163  stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
164  stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
165  (*strings)[astring]->LorentzRotate(toLab);
166  SumStringMom+=(*strings)[astring]->Get4Momentum();
167  #ifdef debug_PartonStringModel
168  G4cout<<"String No "<<astring+1<<" "<<(*strings)[astring]->Get4Momentum()<<" "
169  <<(*strings)[astring]->Get4Momentum().mag()
170  <<" Partons "<<(*strings)[astring]->GetLeftParton()->GetDefinition()->GetPDGEncoding()
171  <<" "<<(*strings)[astring]->GetRightParton()->GetDefinition()->GetPDGEncoding()<<G4endl;
172  #endif
173  }
174  else
175  {
176  stringEnergy += (*strings)[astring]->GetKineticTrack()->Get4Momentum().t();
177  (*strings)[astring]->LorentzRotate(toLab);
178  SumStringMom+=(*strings)[astring]->GetKineticTrack()->Get4Momentum();
179  #ifdef debug_PartonStringModel
180  G4cout<<"A track No "<<astring+1<<" "
181  <<(*strings)[astring]->GetKineticTrack()->Get4Momentum()<<" "
182  <<(*strings)[astring]->GetKineticTrack()->Get4Momentum().mag()<<" "
183  <<(*strings)[astring]->GetKineticTrack()->GetDefinition()->GetParticleName()<<G4endl;
184  #endif
185  }
186  }
187 
188  #ifdef debug_PartonStringModel
189  G4cout<<G4endl<<"SumString4Mom "<<SumStringMom<<G4endl;
190  G4LorentzVector TargetResidual4Momentum(0.,0.,0.,0.);
191  G4LorentzVector ProjectileResidual4Momentum(0.,0.,0.,0.);
192  G4int hitsT(0), charged_hitsT(0);
193  G4int hitsP(0), charged_hitsP(0);
194  G4double ExcitationEt(0.), ExcitationEp(0.);
195  #endif
196 
197  G4V3DNucleus * ProjResNucleus=theThis->GetProjectileNucleus();
198 
199  G4int numberProtonProjectileResidual( 0 ), numberNeutronProjectileResidual( 0 );
200  if(ProjResNucleus != 0)
201  {
202  theNuclNucleon = ProjResNucleus->StartLoop() ? ProjResNucleus->GetNextNucleon() : NULL;
203  G4int numberProtonProjectileHits( 0 ), numberNeutronProjectileHits( 0 );
204  while( theNuclNucleon )
205  {
206  if(theNuclNucleon->AreYouHit())
207  {
208  G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
209  #ifdef debug_PartonStringModel
210  ProjectileResidual4Momentum += tmp;
211  hitsP++;
212  if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hitsP;
213  ExcitationEp +=theNuclNucleon->GetBindingEnergy();
214  #endif
215  theNuclNucleon->SetMomentum(tmp);
216  if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++numberProtonProjectileHits;
217  if ( theNuclNucleon->GetDefinition() == G4Neutron::Neutron() ) ++numberNeutronProjectileHits;
218  }
219  theNuclNucleon = ProjResNucleus->GetNextNucleon();
220  }
221  #ifdef debug_PartonStringModel
222  G4cout<<"Projectile residual A, Z and E* "
223  <<thePrimary.GetDefinition()->GetBaryonNumber() - hitsP<<" "
224  <<thePrimary.GetDefinition()->GetPDGCharge() - charged_hitsP<<" "
225  <<ExcitationEp<<G4endl;
226  G4cout<<"Projectile residual 4 momentum "<<ProjectileResidual4Momentum<<G4endl;
227  #endif
228  numberProtonProjectileResidual = thePrimary.GetDefinition()->GetPDGCharge() - numberProtonProjectileHits;
229  numberNeutronProjectileResidual = thePrimary.GetDefinition()->GetBaryonNumber()
230  - thePrimary.GetDefinition()->GetPDGCharge() - numberNeutronProjectileHits;
231  }
232 
233  G4V3DNucleus * ResNucleus=theThis->GetWoundedNucleus();
234 
235  // loop over wounded nucleus
236  theNuclNucleon = ResNucleus->StartLoop() ? ResNucleus->GetNextNucleon() : NULL;
237  G4int numberProtonTargetHits( 0 ), numberNeutronTargetHits( 0 );
238  while( theNuclNucleon )
239  {
240  if(theNuclNucleon->AreYouHit())
241  {
242  G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
243  #ifdef debug_PartonStringModel
244  TargetResidual4Momentum += tmp;
245  hitsT++;
246  if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hitsT;
247  ExcitationEt +=theNuclNucleon->GetBindingEnergy();
248  #endif
249  theNuclNucleon->SetMomentum(tmp);
250  if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++numberProtonTargetHits;
251  if ( theNuclNucleon->GetDefinition() == G4Neutron::Neutron() ) ++numberNeutronTargetHits;
252  }
253  theNuclNucleon = ResNucleus->GetNextNucleon();
254  }
255 
256  #ifdef debug_PartonStringModel
257  G4cout<<"Target residual A, Z and E* "
258  <<theNucleus.GetA_asInt() - hitsT<<" "
259  <<theNucleus.GetZ_asInt() - charged_hitsT<<" "
260  <<ExcitationEt<<G4endl;
261  G4cout<<"Target residual 4 momentum "<<TargetResidual4Momentum<<G4endl;
262  Bsum+=( hitsT + hitsP);
263  Qsum+=(charged_hitsT + charged_hitsP);
264  G4cout<<"Hitted # of nucleons of projectile and target "<<hitsP<<" "<<hitsT<<G4endl;
265  G4cout<<"Hitted # of protons of projectile and target "
266  <<charged_hitsP<<" "<<charged_hitsT<<G4endl<<G4endl;
267  G4cout<<"Bsum Qsum "<<Bsum<<" "<<Qsum<<G4endl<<G4endl;
268  #endif
269 
270  // Re-sample in the case of unphysical nuclear residual:
271  // 1 (H), 2 (2He), and 3 (3Li) protons alone without neutrons can exist, but not more;
272  // no bound states of 2 or more neutrons without protons can exist.
273  G4int numberProtonTargetResidual = theNucleus.GetZ_asInt() - numberProtonTargetHits;
274  G4int numberNeutronTargetResidual = theNucleus.GetA_asInt() - theNucleus.GetZ_asInt() - numberNeutronTargetHits;
275  G4bool unphysicalResidual = false;
276  if ( ( numberProtonTargetResidual > 3 && numberNeutronTargetResidual == 0 ) ||
277  ( numberProtonTargetResidual == 0 && numberNeutronTargetResidual > 1 ) ) {
278  unphysicalResidual = true;
279  //G4cout << "***UNPHYSICAL TARGET RESIDUAL*** Z=" << numberProtonTargetResidual
280  // << " ; N=" << numberNeutronTargetResidual;
281  }
282  if ( ( numberProtonProjectileResidual > 3 && numberNeutronProjectileResidual == 0 ) ||
283  ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual > 1 ) ) {
284  unphysicalResidual = true;
285  //G4cout << "***UNPHYSICAL PROJECTILE RESIDUAL*** Z=" << numberProtonProjectileResidual
286  // << " ; N=" << numberNeutronProjectileResidual;
287  }
288  if ( unphysicalResidual ) {
289  //G4cout << " -> REJECTING COLLISION because of unphysical residual !" << G4endl;
290  Success = false;
291  continue;
292  }
293 
294  //=========================================================================================
295  // Fragment strings
296  #ifdef debug_PartonStringModel
297  G4cout<<"---------------- Attempt to fragment strings ------------- "<<attempts<<G4endl;
298  #endif
299 
300  G4double InvMass=SumStringMom.mag();
301  G4double SumMass(0.);
302 
303  #ifdef debug_PartonStringModel
304  QsumSec=0; BsumSec=0;
305  SumPsecondr=G4LorentzVector(0.,0.,0.,0.);
306  #endif
307 
308  if(theResult != 0)
309  {
310  std::for_each(theResult->begin(), theResult->end(), DeleteKineticTrack());
311  delete theResult;
312  }
313 
314  theResult = stringFragmentationModel->FragmentStrings(strings);
315 
316  #ifdef debug_PartonStringModel
317  G4cout<<"String fragmentation success (OK - #0, Not OK - 0) : "<<theResult<<G4endl;
318  #endif
319 
320  if(theResult == 0) {Success=false; continue;}
321 
322  #ifdef debug_PartonStringModel
323  G4cout<<"Parton-String model: Number of produced particles "<<theResult->size()<<G4endl;
324  SumPsecondr = G4LorentzVector(0.,0.,0.,0.);
325  QsumSec = 0; BsumSec = 0;
326  #endif
327 
328  SumMass=0.;
329 
330  for ( unsigned int i=0; i < theResult->size(); i++)
331  {
332  SumMass+=(*theResult)[i]->Get4Momentum().mag();
333  #ifdef debug_PartonStringModel
334  G4cout<<i<<" : "<<(*theResult)[i]->GetDefinition()->GetParticleName()<<" "
335  <<(*theResult)[i]->Get4Momentum()<<" "
336  <<(*theResult)[i]->Get4Momentum().mag()<<" "
337  <<(*theResult)[i]->GetDefinition()->GetPDGMass()<<G4endl;
338  SumPsecondr+=(*theResult)[i]->Get4Momentum();
339  BsumSec += (*theResult)[i]->GetDefinition()->GetBaryonNumber();
340  QsumSec += (*theResult)[i]->GetDefinition()->GetPDGCharge();
341  #endif
342  }
343 
344  #ifdef debug_PartonStringModel
345  G4cout<<G4endl<<"-----------------------Parton-String model: balances -------------"<<G4endl;
346  if(Qsum != QsumSec) {
347  G4cout<<"Charge is not conserved!!! ----"<<G4endl;
348  G4cout<<" Qsum != QsumSec "<<Qsum<<" "<<QsumSec<<G4endl;
349  }
350  if(Bsum != BsumSec) {
351  G4cout<<"Baryon number is not conserved!!!"<<G4endl;
352  G4cout<<" Bsum != BsumSec "<<Bsum<<" "<<BsumSec<<G4endl;
353  }
354  #endif
355 
356  if((SumMass > InvMass)||(SumMass == 0.)) {Success=false;}
357  std::for_each(strings->begin(), strings->end(), DeleteString() );
358  delete strings;
359 
360  } while(!Success);
361 
362  #ifdef debug_PartonStringModel
363  G4cout <<"Baryon number balance "<<Bsum-BsumSec<<G4endl;
364  G4cout <<"Charge balance "<<Qsum-QsumSec<<G4endl;
365  G4cout <<"4 momentum balance "<<SumStringMom-SumPsecondr<<G4endl;
366  G4cout<<"---------------------End of Parton-String model work -------------"<<G4endl<<G4endl;
367  #endif
368 
369  return theResult;
370 }
371 
372 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const
373 {
374  outFile << GetModelName() << " has no description yet.\n";
375 }
376 
378 { return 0; }
379