Geant4  10.01.p02
G4VPartonStringModel.cc
Go to the documentation of this file.
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 // $Id: G4VPartonStringModel.cc 83684 2014-09-09 12:37:39Z gcosmo $
28 //
30 // GEANT 4 class implementation file
31 //
32 // ---------------- G4VPartonStringModel ----------------
33 // by Gunter Folger, May 1998.
34 // abstract class for all Parton String Models
35 // ------------------------------------------------------------
36 // debug switch
37 //#define debug_PartonStringModel
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 partyicles are constructed.
52  G4ShortLivedConstructor ShortLived;
53  ShortLived.ConstructParticle();
54 }
55 
57 {
58 }
59 
61  const G4DynamicParticle &aPrimary)
62 {
63  G4ExcitedStringVector * strings = NULL;
64 
65  G4DynamicParticle thePrimary=aPrimary;
66 
67 #ifdef debug_PartonStringModel
68  G4cout<<G4endl;
69  G4cout<<"-----------------------Parton-String model is runnung ------------"<<G4endl;
70  G4cout<<"Projectile Name Mass "<<thePrimary.GetDefinition()->GetParticleName()<<" "
71  <<thePrimary.GetMass()<<G4endl;
72  G4cout<<" Momentum "<<thePrimary.Get4Momentum()<<G4endl;
73  G4cout<<"Target nucleus A Z "<<theNucleus.GetA_asInt()<<" "
74  <<theNucleus.GetZ_asInt()<<G4endl<<G4endl;
75 #endif
76 
78  G4LorentzVector Ptmp=thePrimary.Get4Momentum();
79  toZ.rotateZ(-1*Ptmp.phi());
80  toZ.rotateY(-1*Ptmp.theta());
81  thePrimary.Set4Momentum(toZ*Ptmp);
82  G4LorentzRotation toLab(toZ.inverse());
83 
84  G4int attempts = 0, maxAttempts=20;
85  while ( strings == NULL )
86  {
87  if (attempts++ > maxAttempts )
88  {
89  throw G4HadronicException(__FILE__, __LINE__,
90  "G4VPartonStringModel::Scatter(): fails to generate strings");
91  }
92  theThis->Init(theNucleus,thePrimary);
93 
94  strings = GetStrings();
95  }
96 
97  G4double stringEnergy(0);
98  G4LorentzVector SumStringMom(0.,0.,0.,0.);
99 
100 #ifdef debug_PartonStringModel
101  G4cout<<"Parton-String model: Number of produced strings "<<strings->size()<<G4endl;
102 #endif
103 
104  for ( unsigned int astring=0; astring < strings->size(); astring++)
105  {
106 // rotate string to lab frame, models have it aligned to z
107  if((*strings)[astring]->IsExcited())
108  {
109  stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
110  stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
111  (*strings)[astring]->LorentzRotate(toLab);
112  SumStringMom+=(*strings)[astring]->Get4Momentum();
113 #ifdef debug_PartonStringModel
114 G4cout<<"String No "<<astring+1<<" "<<(*strings)[astring]->Get4Momentum()<<" "
115  <<(*strings)[astring]->Get4Momentum().mag()<<G4endl;
116 #endif
117  }
118  else
119  {
120  stringEnergy += (*strings)[astring]->GetKineticTrack()->Get4Momentum().t();
121  (*strings)[astring]->LorentzRotate(toLab);
122  SumStringMom+=(*strings)[astring]->GetKineticTrack()->Get4Momentum();
123 #ifdef debug_PartonStringModel
124 G4cout<<"A track No "<<astring+1<<" "
125  <<(*strings)[astring]->GetKineticTrack()->Get4Momentum()<<" "
126  <<(*strings)[astring]->GetKineticTrack()->Get4Momentum().mag()<<G4endl;
127 #endif
128  }
129  }
130 
131  G4double InvMass=SumStringMom.mag();
132  G4V3DNucleus * ResNucleus=theThis->GetWoundedNucleus();
133 
134  // loop over wounded nucleus
135  G4Nucleon * theNuclNucleon = ResNucleus->StartLoop() ?
136  ResNucleus->GetNextNucleon() : NULL;
137  while( theNuclNucleon )
138  {
139  if(theNuclNucleon->AreYouHit())
140  {
141  G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
142  theNuclNucleon->SetMomentum(tmp);
143  }
144  theNuclNucleon = ResNucleus->GetNextNucleon();
145  }
146 
147  G4V3DNucleus * ProjResNucleus=theThis->GetProjectileNucleus();
148 
149 #ifdef debug_PartonStringModel
150  G4ThreeVector hitNucleonMomentum(0.,0.,0.);
151 #endif
152 
153  if(ProjResNucleus != 0)
154  {
155  theNuclNucleon = ProjResNucleus->StartLoop() ?
156  ProjResNucleus->GetNextNucleon() : NULL;
157  while( theNuclNucleon )
158  {
159  if(theNuclNucleon->AreYouHit())
160  {
161  G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
162  #ifdef debug_PartonStringModel
163  hitNucleonMomentum += tmp.vect();
164  #endif
165  theNuclNucleon->SetMomentum(tmp);
166  }
167  theNuclNucleon = ProjResNucleus->GetNextNucleon();
168  }
169  }
170 
171 #ifdef debug_PartonStringModel
172  G4cout<<"Parton-String model: SumStringMom "<<SumStringMom<<G4endl;
173 
174  G4V3DNucleus * fancynucleus=theThis->GetWoundedNucleus();
175 
176  // loop over wounded nucleus
177  G4int hits(0), charged_hits(0);
178 // G4ThreeVector hitNucleonMomentum(0.,0.,0.); // Uzhi Feb. 2014
179  G4Nucleon * theCurrentNucleon = fancynucleus->StartLoop() ?
180  fancynucleus->GetNextNucleon() : NULL;
181  while( theCurrentNucleon )
182  {
183  if(theCurrentNucleon->AreYouHit())
184  {
185  hits++;
186  hitNucleonMomentum += theCurrentNucleon->Get4Momentum().vect();
187  if ( theCurrentNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hits;
188  }
189  theCurrentNucleon = fancynucleus->GetNextNucleon();
190  }
191 
192  G4int initialZ=fancynucleus->GetCharge();
193  G4int initialA=fancynucleus->GetMassNumber();
194  G4double initial_mass=
196  G4double final_mass(0.);
197  if(initialA-hits != 0) final_mass =
199  initialZ-charged_hits, initialA-hits);
200  G4cout << "G4VPSM: " <<G4endl
201  << "strE "<<stringEnergy <<G4endl
202  << "Hit targeet nucleons "<<hits <<G4endl
203  << "Primary "<<Ptmp.e() <<G4endl
204  << "SumStringE "<<SumStringMom.e() <<G4endl
205  << "Target nucleus intial "<<initial_mass <<G4endl
206  << "Target nucleus final "<<final_mass <<G4endl
207  << "Excitation estimate "
208  <<Ptmp.e() + initial_mass - final_mass - stringEnergy << G4endl;
209  G4cout << "momentum balance "
210  << thePrimary.GetMomentum() + hitNucleonMomentum -
211  SumStringMom.vect()<< G4endl;
212 #endif
213 
214 // Fragment strings
215 
216  G4KineticTrackVector * theResult = 0;
217  G4double SumMass(0.);
218  attempts = 0;
219  maxAttempts=100;
220  do
221  {
222  attempts++;
223  if(theResult != 0)
224  {
225  std::for_each(theResult->begin(), theResult->end(), DeleteKineticTrack());
226  delete theResult;
227  }
228 
229  theResult = stringFragmentationModel->FragmentStrings(strings);
230 
231  if(attempts > maxAttempts ) break;
232 
233 #ifdef debug_PartonStringModel
234  G4cout<<"Parton-String model: Number of produced particles "<<theResult->size()<<G4endl;
235  G4LorentzVector SumPsecondr(0.,0.,0.,0.);
236 #endif
237 
238  SumMass=0.;
239 
240  for ( unsigned int i=0; i < theResult->size(); i++)
241  {
242  SumMass+=(*theResult)[i]->GetDefinition()->GetPDGMass();
243  //SumP+=(*theResult)[i]->Get4Momentum();
244 #ifdef debug_PartonStringModel
245  G4cout<<i<<" : "<<(*theResult)[i]->GetDefinition()->GetParticleName()<<" "
246  <<(*theResult)[i]->Get4Momentum()<<" "
247  <<(*theResult)[i]->Get4Momentum().mag()<<G4endl;
248  SumPsecondr+=(*theResult)[i]->Get4Momentum();
249 #endif
250  }
251 #ifdef debug_PartonStringModel
252  G4cout<<"SumP secondaries "<<SumPsecondr<<G4endl;
253 #endif
254  } while(SumMass > InvMass);
255 
256  std::for_each(strings->begin(), strings->end(), DeleteString() );
257  delete strings;
258 
259 #ifdef debug_PartonStringModel
260  G4cout<<"End of string model work ------------"<<G4endl<<G4endl;
261 #endif
262  return theResult;
263 }
264 
265 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const
266 {
267  outFile << GetModelName() << " has no description yet.\n";
268 }
269 
271 { return 0;}
272 
G4VPartonStringModel(const G4String &modelName="Parton String Model")
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
std::vector< G4ExcitedString * > G4ExcitedStringVector
virtual G4KineticTrackVector * FragmentStrings(const G4ExcitedStringVector *theStrings)=0
G4VPartonStringModel * theThis
virtual G4int GetCharge()=0
void SetMomentum(G4LorentzVector &aMomentum)
Definition: G4Nucleon.hh:70
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4bool StartLoop()=0
virtual G4V3DNucleus * GetWoundedNucleus() const =0
virtual G4int GetMassNumber()=0
virtual G4V3DNucleus * GetProjectileNucleus() const
virtual const G4LorentzVector & Get4Momentum() const
Definition: G4Nucleon.hh:72
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4KineticTrackVector * Scatter(const G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
virtual void ModelDescription(std::ostream &outFile) const
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
G4VStringFragmentation * stringFragmentationModel
G4double GetMass() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1249
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
virtual G4String GetModelName() const
static G4ParticleTable * GetParticleTable()
virtual void Init(const G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)=0
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
virtual G4ExcitedStringVector * GetStrings()=0
#define G4endl
Definition: G4ios.hh:61
virtual G4Nucleon * GetNextNucleon()=0
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetMomentum() const
CLHEP::HepLorentzVector G4LorentzVector