Geant4  10.02.p01
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 }
47 
49 : G4VParticipants(),ModelMode(right.ModelMode), nCutMax(right.nCutMax),
50  ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
51  theNucleonRadius(right.theNucleonRadius)
52 {
53 }
54 
55 
57 {
58 }
59 
61 {
62 
63  // Find the collisions and collition conditions
64  G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary);
65 
66  // now build the parton pairs. HPW
67  SplitHadrons();
68 
69  // soft collisions first HPW, ordering is vital
71 
72  // the rest is diffractive HPW
74 
75  // clean-up, if necessary
76  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
77  theInteractions.clear();
78  std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
79  theTargets.clear();
80  delete aProjectile;
81 }
82 
84 {
85  G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE); // @@@ check the TRUE
86  G4PomeronCrossSection theProbability(thePrimary.GetDefinition()); // @@@ should be data member
87  G4double outerRadius = theNucleus->GetOuterRadius();
88  // Check reaction threshold
89 
91  G4Nucleon * pNucleon = theNucleus->GetNextNucleon();
92  G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
93  //--DEBUG--G4cout << " qgspart- " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
94  //--DEBUG-- << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
95  G4double s_nucleus = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
96  G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass();
97  ModelMode = SOFT;
98  if (sqr(ThresholdMass + ThresholdParameter) > s_nucleus)
99  {
101  //throw G4HadronicException(__FILE__, __LINE__, "Initial energy is too low. The 4-vectors of the input are inconsistant with the particle masses.");
102  }
103  if (sqr(ThresholdMass + QGSMThreshold) > s_nucleus) // thus only diffractive in cascade!
104  {
106  }
107 
108  // first find the collisions HPW
109  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
110  theInteractions.clear();
111  G4int totalCuts = 0;
112 
113  #ifdef debug_QGS
114  G4double eK = thePrimary.GetKineticEnergy()/GeV;
115  #endif
116  #ifdef debug_G4QGSParticipants
117  G4double impactUsed = 0;
118  G4LorentzVector intNuclMom;
119  #endif
120 
121  const G4int maxNumberOfLoops = 1000;
122  G4int loopCounter = -1;
123  while ( (theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 26.10.2015, A.Ribon */
124  {
125  // choose random impact parameter HPW
126  std::pair<G4double, G4double> theImpactParameter;
127  theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
128  G4double impactX = theImpactParameter.first;
129  G4double impactY = theImpactParameter.second;
130 
131  // loop over nucleons to find collisions
133  G4int nucleonCount = 0; // debug
135  #ifdef debug_G4QGSParticipants
136  intNuclMom=aPrimaryMomentum;
137  #endif
138  while( (pNucleon = theNucleus->GetNextNucleon()) ) /* Loop checking, 26.10.2015, A.Ribon */
139  {
140  if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV)
141  {
142  break;
143  }
144  nucleonCount++; // debug
145  // Needs to be moved to Probability class @@@
146  G4LorentzVector nucleonMomentum=pNucleon->Get4Momentum();
147  nucleonMomentum.setE(nucleonMomentum.e()-pNucleon->GetBindingEnergy());
148  G4double s_nucleon = (aPrimaryMomentum + nucleonMomentum).mag2();
149  G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
150  sqr(impactY - pNucleon->GetPosition().y());
151  G4double Probability = theProbability.GetInelasticProbability(s_nucleon, Distance2);
152  // test for inelastic collision
153  G4double rndNumber = G4UniformRand();
154  // ModelMode = DIFFRACTIVE;
155  if (Probability > rndNumber)
156  {
157  #ifdef debug_G4QGSParticipants
158  G4cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl;
159  G4cout << " qgspart+ " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
160  << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
161  intNuclMom += nucleonMomentum;
162  #endif
163  pNucleon->SetMomentum(nucleonMomentum);
164  G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
166  theTargets.push_back(aTarget);
167  pNucleon->Hit(aTarget);
168  if ((theProbability.GetDiffractiveProbability(s_nucleon, Distance2)/Probability > G4UniformRand()
169  &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ))
170  {
171  // diffractive interaction occurs
172  if(IsSingleDiffractive())
173  {
174  theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget);
175  }
176  else
177  {
178  theDiffExcitaton.ExciteParticipants(aProjectile, aTarget);
179  }
180  G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
181  aInteraction->SetTarget(aTarget);
182  theInteractions.push_back(aInteraction);
183  aInteraction->SetNumberOfDiffractiveCollisions(1);
184  totalCuts += 1;
185  }
186  else
187  {
188  // nondiffractive soft interaction occurs
189  // sample nCut+1 (cut Pomerons) pairs of strings can be produced
190  G4int nCut;
191  G4double * running = new G4double[nCutMax];
192  running[0] = 0;
193  for(nCut = 0; nCut < nCutMax; nCut++)
194  {
195  running[nCut] = theProbability.GetCutPomeronProbability(s_nucleon, Distance2, nCut + 1);
196  if(nCut!=0) running[nCut] += running[nCut-1];
197  }
198  G4double random = running[nCutMax-1]*G4UniformRand();
199  for(nCut = 0; nCut < nCutMax; nCut++)
200  {
201  if(running[nCut] > random) break;
202  }
203  delete [] running;
204  nCut = 0;
205  aTarget->IncrementCollisionCount(nCut+1);
206  aProjectile->IncrementCollisionCount(nCut+1);
207  G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
208  aInteraction->SetTarget(aTarget);
209  aInteraction->SetNumberOfSoftCollisions(nCut+1);
210  theInteractions.push_back(aInteraction);
211  totalCuts += nCut+1;
212  #ifdef debug_G4QGSParticipants
213  impactUsed=Distance2;
214  #endif
215  }
216  }
217  }
218  #ifdef debug_G4QGSParticipants
219  G4cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl;
220  G4cout << " Interact 4-Vect " << intNuclMom << G4endl;
221  #endif
222  }
223  if ( loopCounter >= maxNumberOfLoops ) {
225  ed << " loopCounter exceeds maxNumberOfLoops : forced exit! " << G4endl;
226  G4Exception( "G4QGSParticipants::SelectInteractions ", "HAD_QGS_001", JustWarning, ed );
227  }
228  #ifdef debug_G4QGSParticipants
229  G4cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl;
230  G4cout << "Impact Parameter used = "<<impactUsed<<G4endl;
231  #endif
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(),
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  thePartonPairs.push_back(aPartonPair);
260  }
261  // then target HPW
262  G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
263  aParton = aTarget->GetNextParton();
264  if (aParton)
265  {
266  aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(),
269  #ifdef debug_G4QGSPart_PDiffColl
270  G4cout << "DiffPair Tgt " << aPartonPair->GetParton1()->GetPDGcode() << " "
271  << aPartonPair->GetParton1()->Get4Momentum() << " "
272  << aPartonPair->GetParton1()->GetX() << " " << G4endl;
273  G4cout << " " << aPartonPair->GetParton2()->GetPDGcode() << " "
274  << aPartonPair->GetParton2()->Get4Momentum() << " "
275  << aPartonPair->GetParton2()->GetX() << " " << G4endl;
276  #endif
277  thePartonPairs.push_back(aPartonPair);
278  }
279  }
280 }
281 
283 {
284  std::vector<G4InteractionContent*>::iterator i;
285  G4LorentzVector str4Mom;
286  i = theInteractions.begin();
287  while ( i != theInteractions.end() ) /* Loop checking, 10.08.2015, A.Ribon */
288  {
289  G4InteractionContent* anIniteraction = *i;
290  G4PartonPair * aPair = NULL;
291  if (anIniteraction->GetNumberOfSoftCollisions())
292  {
293  G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
294  G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
295  for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
296  {
297  aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
299  #ifdef debug_G4QGSPart_PSoftColl
300  G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
301  << aPair->GetParton1()->Get4Momentum() << " "
302  << aPair->GetParton1()->GetX() << " " << G4endl;
303  G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
304  << aPair->GetParton2()->Get4Momentum() << " "
305  << aPair->GetParton2()->GetX() << " " << G4endl;
306  #endif
307  #ifdef debug_G4QGSParticipants
308  str4Mom += aPair->GetParton1()->Get4Momentum();
309  str4Mom += aPair->GetParton2()->Get4Momentum();
310  #endif
311  thePartonPairs.push_back(aPair);
312  aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
314  #ifdef debug_G4QGSPart_PSoftColl
315  G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
316  << aPair->GetParton1()->Get4Momentum() << " "
317  << aPair->GetParton1()->GetX() << " " << G4endl;
318  G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
319  << aPair->GetParton2()->Get4Momentum() << " "
320  << aPair->GetParton2()->GetX() << " " << G4endl;
321  #endif
322  #ifdef debug_G4QGSParticipants
323  str4Mom += aPair->GetParton1()->Get4Momentum();
324  str4Mom += aPair->GetParton2()->Get4Momentum();
325  #endif
326  thePartonPairs.push_back(aPair);
327  }
328  delete *i;
329  i=theInteractions.erase(i); // i now points to the next interaction
330  } else {
331  i++;
332  }
333  }
334  #ifdef debug_G4QGSPart_PSoftColl
335  G4cout << " string 4 mom " << str4Mom << G4endl;
336  #endif
337 }
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
G4Parton * GetParton2(void)
Definition: G4PartonPair.hh:79
G4int GetPDGcode() const
Definition: G4Parton.hh:124
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void SetMomentum(G4LorentzVector &aMomentum)
Definition: G4Nucleon.hh:70
virtual G4bool StartLoop()=0
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:140
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:74
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)
static const double GeV
Definition: G4SIunits.hh:214
G4QGSDiffractiveExcitation theDiffExcitaton
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
G4V3DNucleus * theNucleus
G4double GetX()
Definition: G4Parton.hh:86
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
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
G4SingleDiffractiveExcitation theSingleDiffExcitation
G4double GetBindingEnergy() const
Definition: G4Nucleon.hh:75
std::vector< G4InteractionContent * > theInteractions
const G4double theNucleonRadius
G4double GetMass() const
const G4double ThresholdParameter
static const double fermi
Definition: G4SIunits.hh:102
G4ThreadLocal G4int G4QGSParticipants_NPart
void SetTarget(G4VSplitableHadron *aTarget)
CLHEP::HepLorentzVector G4LorentzVector