Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QGSParticipants.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 #include <utility>
27 
28 #include "G4QGSParticipants.hh"
29 #include "globals.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4LorentzVector.hh"
32 
33 // Class G4QGSParticipants
34 
35 // HPW Feb 1999
36 // Promoting model parameters from local variables class properties
38 
39 G4QGSParticipants::G4QGSParticipants() : theDiffExcitaton(), //0.7*GeV, 250*MeV, 250*MeV),
40  ModelMode(SOFT),
41  //nCutMax(7),ThresholdParameter(0.45*GeV),
42  nCutMax(7),ThresholdParameter(0.000*GeV),
43  QGSMThreshold(3*GeV),theNucleonRadius(1.5*fermi)
44 {
45 }
46 
48 : G4VParticipants(),ModelMode(right.ModelMode), nCutMax(right.nCutMax),
49  ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
50  theNucleonRadius(right.theNucleonRadius)
51 {
52 }
53 
54 
56 {
57 }
58 
60 {
61  // Find the collisions and collition conditions
62  G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary);
63 
64  // now build the parton pairs. HPW
65  SplitHadrons();
66 
67  // soft collisions first HPW, ordering is vital
69 
70  // the rest is diffractive HPW
72 
73  // clean-up, if necessary
74  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
75  theInteractions.clear();
76  std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
77  theTargets.clear();
78  delete aProjectile;
79 }
80 
82 {
83  G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE); // @@@ check the TRUE
84  G4PomeronCrossSection theProbability(thePrimary.GetDefinition()); // @@@ should be data member
85  G4double outerRadius = theNucleus->GetOuterRadius();
86  // Check reaction threshold
87 
89  G4Nucleon * pNucleon = theNucleus->GetNextNucleon();
90  G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
91  //--DEBUG--G4cout << " qgspart- " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
92  //--DEBUG-- << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
93  G4double s_nucleus = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
94  G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass();
95  ModelMode = SOFT;
96  if (sqr(ThresholdMass + ThresholdParameter) > s_nucleus)
97  {
99  //throw G4HadronicException(__FILE__, __LINE__,
100  // "Initial energy is too low. The 4-vectors of the input are inconsistant with the particle masses.");
101  }
102  if (sqr(ThresholdMass + QGSMThreshold) > s_nucleus) // thus only diffractive in cascade!
103  {
105  }
106 
107  // first find the collisions HPW
108  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
109  theInteractions.clear();
110  G4int totalCuts = 0;
111 
112  #ifdef debug_QGS
113  G4double eK = thePrimary.GetKineticEnergy()/GeV;
114  #endif
115  #ifdef debug_G4QGSParticipants
116  G4double impactUsed = 0;
117  G4LorentzVector intNuclMom;
118  #endif
119 
120  const G4int maxNumberOfLoops = 1000;
121  G4int loopCounter = -1;
122  while ( (theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 26.10.2015, A.Ribon */
123  {
124  // choose random impact parameter HPW
125  std::pair<G4double, G4double> theImpactParameter;
126  theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
127  G4double impactX = theImpactParameter.first;
128  G4double impactY = theImpactParameter.second;
129 
130  // loop over nucleons to find collisions
132  G4int nucleonCount = 0; // debug
134  #ifdef debug_G4QGSParticipants
135  intNuclMom=aPrimaryMomentum;
136  #endif
137  while( (pNucleon = theNucleus->GetNextNucleon()) ) /* Loop checking, 26.10.2015, A.Ribon */
138  {
139  if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV)
140  {
141  break;
142  }
143  nucleonCount++; // debug
144  // Needs to be moved to Probability class @@@
145  G4LorentzVector nucleonMomentum=pNucleon->Get4Momentum();
146  nucleonMomentum.setE(nucleonMomentum.e()-pNucleon->GetBindingEnergy());
147  G4double s_nucleon = (aPrimaryMomentum + nucleonMomentum).mag2();
148  G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
149  sqr(impactY - pNucleon->GetPosition().y());
150  G4double Probability = theProbability.GetInelasticProbability(s_nucleon, Distance2);
151  // test for inelastic collision
152  G4double rndNumber = G4UniformRand();
153  // ModelMode = DIFFRACTIVE;
154  if (Probability > rndNumber)
155  {
156  #ifdef debug_G4QGSParticipants
157  G4cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl;
158  G4cout << " qgspart+ " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
159  << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
160  intNuclMom += nucleonMomentum;
161  #endif
162  pNucleon->SetMomentum(nucleonMomentum);
163  G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
165  theTargets.push_back(aTarget);
166  pNucleon->Hit(aTarget);
167  if ((theProbability.GetDiffractiveProbability(s_nucleon, Distance2)/Probability > G4UniformRand()
168  &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ))
169  {
170  // diffractive interaction occurs
171  if(IsSingleDiffractive())
172  {
173  theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget);
174  } else {
175  theDiffExcitaton.ExciteParticipants(aProjectile, aTarget);
176  }
177  G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
178  aInteraction->SetTarget(aTarget);
179  theInteractions.push_back(aInteraction);
180  aInteraction->SetNumberOfDiffractiveCollisions(1);
181  totalCuts += 1;
182  } else {
183  // nondiffractive soft interaction occurs
184  // sample nCut+1 (cut Pomerons) pairs of strings can be produced
185  G4int nCut;
186  G4double * running = new G4double[nCutMax];
187  running[0] = 0;
188  for(nCut = 0; nCut < nCutMax; nCut++)
189  {
190  running[nCut] = theProbability.GetCutPomeronProbability(s_nucleon, Distance2, nCut + 1);
191  if(nCut!=0) running[nCut] += running[nCut-1];
192  }
193  G4double random = running[nCutMax-1]*G4UniformRand();
194  for(nCut = 0; nCut < nCutMax; nCut++)
195  {
196  if(running[nCut] > random) break;
197  }
198  delete [] running;
199  nCut = 0;
200  aTarget->IncrementCollisionCount(nCut+1);
201  aProjectile->IncrementCollisionCount(nCut+1);
202  G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
203  aInteraction->SetTarget(aTarget);
204  aInteraction->SetNumberOfSoftCollisions(nCut+1);
205  theInteractions.push_back(aInteraction);
206  totalCuts += nCut+1;
207  #ifdef debug_G4QGSParticipants
208  impactUsed=Distance2;
209  #endif
210  }
211  }
212  }
213 
214  #ifdef debug_G4QGSParticipants
215  G4cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl;
216  G4cout << " Interact 4-Vect " << intNuclMom << G4endl;
217  #endif
218 
219  }
220 
221  if ( loopCounter >= maxNumberOfLoops ) {
223  ed << " loopCounter exceeds maxNumberOfLoops : forced exit! " << G4endl;
224  G4Exception( "G4QGSParticipants::SelectInteractions ", "HAD_QGS_001", JustWarning, ed );
225  }
226 
227  #ifdef debug_G4QGSParticipants
228  G4cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl;
229  G4cout << "Impact Parameter used = "<<impactUsed<<G4endl;
230  #endif
231 
232  return aProjectile;
233 }
234 
236 {
237  // remove the "G4PartonPair::PROJECTILE", etc., which are not necessary. @@@
238  unsigned int i;
239  for(i = 0; i < theInteractions.size(); i++)
240  {
241  G4InteractionContent* anIniteraction = theInteractions[i];
242  G4VSplitableHadron* aProjectile = anIniteraction->GetProjectile();
243  G4Parton* aParton = aProjectile->GetNextParton();
244  G4PartonPair * aPartonPair;
245  // projectile first HPW
246  if (aParton)
247  {
248  aPartonPair = new G4PartonPair(aParton, aProjectile->GetNextAntiParton(),
250 
251  #ifdef debug_G4QGSPart_PDiffColl
252  G4cout << "DiffPair Pro " << aPartonPair->GetParton1()->GetPDGcode() << " "
253  << aPartonPair->GetParton1()->Get4Momentum() << " "
254  << aPartonPair->GetParton1()->GetX() << " " << G4endl;
255  G4cout << " " << aPartonPair->GetParton2()->GetPDGcode() << " "
256  << aPartonPair->GetParton2()->Get4Momentum() << " "
257  << aPartonPair->GetParton2()->GetX() << " " << G4endl;
258  #endif
259 
260  thePartonPairs.push_back(aPartonPair);
261  }
262  // then target HPW
263  G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
264  aParton = aTarget->GetNextParton();
265  if (aParton)
266  {
267  aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(),
269 
270  #ifdef debug_G4QGSPart_PDiffColl
271  G4cout << "DiffPair Tgt " << aPartonPair->GetParton1()->GetPDGcode() << " "
272  << aPartonPair->GetParton1()->Get4Momentum() << " "
273  << aPartonPair->GetParton1()->GetX() << " " << G4endl;
274  G4cout << " " << aPartonPair->GetParton2()->GetPDGcode() << " "
275  << aPartonPair->GetParton2()->Get4Momentum() << " "
276  << aPartonPair->GetParton2()->GetX() << " " << G4endl;
277  #endif
278 
279  thePartonPairs.push_back(aPartonPair);
280  }
281  }
282 }
283 
285 {
286  std::vector<G4InteractionContent*>::iterator i;
287  G4LorentzVector str4Mom;
288  i = theInteractions.begin();
289  while ( i != theInteractions.end() ) /* Loop checking, 10.08.2015, A.Ribon */
290  {
291  G4InteractionContent* anIniteraction = *i;
292  G4PartonPair * aPair = NULL;
293  if (anIniteraction->GetNumberOfSoftCollisions())
294  {
295  G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
296  G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
297  for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
298  {
299  aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
301 
302  #ifdef debug_G4QGSPart_PSoftColl
303  G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
304  << aPair->GetParton1()->Get4Momentum() << " "
305  << aPair->GetParton1()->GetX() << " " << G4endl;
306  G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
307  << aPair->GetParton2()->Get4Momentum() << " "
308  << aPair->GetParton2()->GetX() << " " << G4endl;
309  #endif
310  #ifdef debug_G4QGSParticipants
311  str4Mom += aPair->GetParton1()->Get4Momentum();
312  str4Mom += aPair->GetParton2()->Get4Momentum();
313  #endif
314 
315  thePartonPairs.push_back(aPair);
316  aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
318 
319  #ifdef debug_G4QGSPart_PSoftColl
320  G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
321  << aPair->GetParton1()->Get4Momentum() << " "
322  << aPair->GetParton1()->GetX() << " " << G4endl;
323  G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
324  << aPair->GetParton2()->Get4Momentum() << " "
325  << aPair->GetParton2()->GetX() << " " << G4endl;
326  #endif
327  #ifdef debug_G4QGSParticipants
328  str4Mom += aPair->GetParton1()->Get4Momentum();
329  str4Mom += aPair->GetParton2()->Get4Momentum();
330  #endif
331 
332  thePartonPairs.push_back(aPair);
333  }
334  delete *i;
335  i=theInteractions.erase(i); // i now points to the next interaction
336  } else {
337  i++;
338  }
339  }
340 
341  #ifdef debug_G4QGSPart_PSoftColl
342  G4cout << " string 4 mom " << str4Mom << G4endl;
343  #endif
344 }
345 
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
G4Parton * GetParton2(void)
Definition: G4PartonPair.hh:76
G4int GetPDGcode() const
Definition: G4Parton.hh:127
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void SetMomentum(G4LorentzVector &aMomentum)
Definition: G4Nucleon.hh:70
double x() const
virtual G4bool StartLoop()=0
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
const G4double QGSMThreshold
virtual const G4LorentzVector & Get4Momentum() const
Definition: G4Nucleon.hh:72
virtual G4Parton * GetNextAntiParton()=0
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
G4Parton * GetParton1(void)
Definition: G4PartonPair.hh:71
virtual G4double GetOuterRadius()=0
virtual ~G4QGSParticipants()
const G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
virtual G4VSplitableHadron * SelectInteractions(const G4ReactionProduct &thePrimary)
virtual G4Parton * GetNextParton()=0
#define TRUE
Definition: globals.hh:55
void IncrementCollisionCount(G4int aCount)
G4QGSDiffractiveExcitation theDiffExcitaton
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
G4V3DNucleus * theNucleus
G4double GetX()
Definition: G4Parton.hh:87
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetKineticEnergy() const
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
Definition: G4V3DNucleus.hh:87
G4double GetTotalEnergy() const
void SetNumberOfDiffractiveCollisions(int)
std::vector< G4PartonPair * > thePartonPairs
G4double GetPDGMass() const
std::vector< G4VSplitableHadron * > theTargets
double y() const
static constexpr double GeV
Definition: G4SIunits.hh:217
G4VSplitableHadron * GetTarget() const
G4VSplitableHadron * GetProjectile() const
void BuildInteractions(const G4ReactionProduct &thePrimary)
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
virtual G4Nucleon * GetNextNucleon()=0
void Hit(G4VSplitableHadron *aHit)
Definition: G4Nucleon.hh:90
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103
G4SingleDiffractiveExcitation theSingleDiffExcitation
G4double GetBindingEnergy() const
Definition: G4Nucleon.hh:75
std::vector< G4InteractionContent * > theInteractions
const G4double theNucleonRadius
G4double GetMass() const
const G4double ThresholdParameter
G4ThreadLocal G4int G4QGSParticipants_NPart
void SetTarget(G4VSplitableHadron *aTarget)