Geant4  10.02.p02
G4FTFParticipants.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: G4FTFParticipants.cc 91914 2015-08-11 07:00:39Z gcosmo $
28 // GEANT4 tag $Name: $
29 //
30 
31 // ------------------------------------------------------------
32 // GEANT 4 class implementation file
33 //
34 // ---------------- G4FTFParticipants----------------
35 // by Gunter Folger, June 1998.
36 // class finding colliding particles in FTFPartonStringModel
37 // Changed in a part by V. Uzhinsky in oder to put in correcpondence
38 // with original FRITIOF mode. November - December 2006.
39 // Ajusted for (anti) nucleus - nucleus interactions by V. Uzhinsky.
40 // (February 2011)
41 // ------------------------------------------------------------
42 
43 #include <utility>
44 #include <vector>
45 #include <algorithm>
46 
47 #include "G4FTFParticipants.hh"
48 #include "G4ios.hh"
49 #include "Randomize.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4FTFParameters.hh"
53 #include "G4VSplitableHadron.hh"
54 
55 
56 //============================================================================
57 
58 //#define debugFTFparticipant
59 
60 
61 //============================================================================
62 
63 G4FTFParticipants::G4FTFParticipants() : currentInteraction( -1 ) {}
64 
65 
66 //============================================================================
67 
69  G4VParticipants(), currentInteraction( -1 )
70 {
71  G4Exception( "G4FTFParticipants::G4FTFParticipants()", "HAD_FTF_001",
72  FatalException, " Must not use copy ctor()" );
73 }
74 
75 
76 //============================================================================
77 
79 
80 
81 //============================================================================
82 
84  G4FTFParameters* theParameters ) {
85 
86  #ifdef debugFTFparticipant
87  G4cout << "Participants::GetList" << G4endl
88  << "thePrimary " << thePrimary.GetMomentum() << G4endl << G4endl;
89  #endif
90 
91  G4double betta_z = thePrimary.GetMomentum().z() / thePrimary.GetTotalEnergy();
92  if ( betta_z < 1.0e-10 ) betta_z = 1.0e-10;
93 
94  StartLoop(); // reset Loop over Interactions
95 
96  for ( unsigned int i = 0; i < theInteractions.size(); i++ ) delete theInteractions[i];
97  theInteractions.clear();
98 
99  G4double deltaxy = 2.0 * fermi; // Extra nuclear radius
100 
101  if ( theProjectileNucleus == 0 ) { // Hadron-nucleus or anti-baryon-nucleus interactions
102 
103  G4double impactX( 0.0 ), impactY( 0.0 );
104 
105  G4VSplitableHadron* primarySplitable = new G4DiffractiveSplitableHadron( thePrimary );
106 
107  #ifdef debugFTFparticipant
108  G4cout << "Hadron-nucleus or anti-baryon-nucleus interactions" << G4endl;
109  #endif
110 
111  G4double xyradius;
112  xyradius = theNucleus->GetOuterRadius() + deltaxy; // Range of impact parameter sampling
113 
114  const G4int maxNumberOfLoops = 1000;
115  G4int loopCounter = 0;
116  do {
117 
118  std::pair< G4double, G4double > theImpactParameter;
119  theImpactParameter = theNucleus->ChooseImpactXandY( xyradius );
120  impactX = theImpactParameter.first;
121  impactY = theImpactParameter.second;
122 
123  #ifdef debugFTFparticipant
124  G4cout << "New interaction list," << " b= "
125  << std::sqrt( sqr(impactX ) + sqr( impactY ) )/fermi << G4endl;
126  #endif
127 
128  G4ThreeVector thePosition( impactX, impactY, 0.0 );
129  primarySplitable->SetPosition( thePosition );
130 
133 
134  #ifdef debugFTFparticipant
135  G4int TrN( 0 );
136  #endif
137 
138  while ( ( nucleon = theNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
139 
140  G4double impact2 = sqr( impactX - nucleon->GetPosition().x() ) +
141  sqr( impactY - nucleon->GetPosition().y() );
142 
143  if ( theParameters->GetProbabilityOfInteraction( impact2/fermi/fermi ) >
144  G4UniformRand() ) {
145  primarySplitable->SetStatus( 1 ); // It takes part in the interaction
146  G4VSplitableHadron* targetSplitable = 0;
147  if ( ! nucleon->AreYouHit() ) {
148  targetSplitable = new G4DiffractiveSplitableHadron( *nucleon );
149  nucleon->Hit( targetSplitable );
150  targetSplitable->SetStatus( 1 ); // It takes part in the interaction
151 
152  #ifdef debugFTFparticipant
153  G4cout << "Participated nucleons #, " << TrN << " " << "Splitable Pr* Tr* "
154  << primarySplitable << " " << targetSplitable << G4endl;
155  #endif
156 
157  }
158  G4InteractionContent* aInteraction = new G4InteractionContent( primarySplitable );
159  G4Nucleon* PrNucleon = 0;
160  aInteraction->SetProjectileNucleon( PrNucleon );
161  aInteraction->SetTarget( targetSplitable );
162  aInteraction->SetTargetNucleon( nucleon );
163  aInteraction->SetStatus( 1 );
164  aInteraction->SetInteractionTime( ( primarySplitable->GetPosition().z() +
165  nucleon->GetPosition().z() ) / betta_z );
166  theInteractions.push_back( aInteraction );
167  }
168 
169  #ifdef debugFTFparticipant
170  TrN++;
171  #endif
172 
173  }
174 
175  } while ( ( theInteractions.size() == 0 ) &&
176  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
177  if ( loopCounter >= maxNumberOfLoops ) {
178  #ifdef debugFTFparticipant
179  G4cout << "BAD situation: forced exit from the while loop!" << G4endl;
180  #endif
181  return;
182  }
183 
184  #ifdef debugFTFparticipant
185  G4cout << "Number of Hit nucleons " << theInteractions.size() << "\t Bx " << impactX/fermi
186  << "\t By " << impactY/fermi << "\t B "
187  << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl << G4endl;
188  #endif
189 
190  //SortInteractionsIncT(); // Not need because nucleons are sorted in increasing z-coordinates.
191  ShiftInteractionTime(); // To put correct times and z-coordinates
192  return;
193 
194  } // end of if ( theProjectileNucleus == 0 )
195 
196  // Projectile and target are nuclei
197 
198  #ifdef debugFTFparticipant
199  G4cout << "Projectile and target are nuclei" << G4endl;
200  #endif
201 
202 //G4cout<<theProjectileNucleus->GetOuterRadius()/fermi<<" "<<theNucleus->GetOuterRadius()/fermi<<" "<<deltaxy/fermi<<G4endl;
203 
204  G4double xyradius;
205  xyradius = theProjectileNucleus->GetOuterRadius() + // Range of impact parameter sampling
206  theNucleus->GetOuterRadius() + deltaxy;
207 
208  G4double impactX( 0.0 ), impactY( 0.0 );
209 
210  const G4int maxNumberOfLoops = 1000;
211  G4int loopCounter = 0;
212  do {
213 
214  std::pair< G4double, G4double > theImpactParameter;
215  theImpactParameter = theNucleus->ChooseImpactXandY( xyradius );
216  impactX = theImpactParameter.first;
217  impactY = theImpactParameter.second;
218 
219  #ifdef debugFTFparticipant
220  G4cout << "New interaction list, " << "b "
221  << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl;
222  #endif
223 
224  G4ThreeVector theBeamPosition( impactX, impactY, 0.0 );
225 
227  G4Nucleon* ProjectileNucleon;
228 
229  #ifdef debugFTFparticipant
230  G4int PrNuclN( 0 );
231  #endif
232 
233  while ( ( ProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
234 
235  G4VSplitableHadron* ProjectileSplitable = 0;
237  G4Nucleon* TargetNucleon = 0;
238 
239  #ifdef debugFTFparticipant
240  G4int TrNuclN( 0 );
241  #endif
242 
243  while ( ( TargetNucleon = theNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
244 
245  G4double impact2 = sqr( impactX + ProjectileNucleon->GetPosition().x() -
246  TargetNucleon->GetPosition().x() ) +
247  sqr( impactY + ProjectileNucleon->GetPosition().y() -
248  TargetNucleon->GetPosition().y() );
249  G4VSplitableHadron* TargetSplitable = 0;
250  if ( theParameters->GetProbabilityOfInteraction( impact2/fermi/fermi ) >
251  G4UniformRand() ) { // An interaction has happend!
252 
253  #ifdef debugFTFparticipant
254  G4cout << G4endl << "An Interaction has happend" << G4endl << "Proj N mom " << PrNuclN
255  << " " << ProjectileNucleon->Get4Momentum() << "-------------" << G4endl
256  << "Targ N mom " << TrNuclN << " " << TargetNucleon->Get4Momentum() << G4endl
257  << "PrN TrN Z coords " << ProjectileNucleon->GetPosition().z()/fermi
258  << " " << TargetNucleon->GetPosition().z()/fermi
259  << " " << ProjectileNucleon->GetPosition().z()/fermi +
260  TargetNucleon->GetPosition().z()/fermi << G4endl;
261  #endif
262 
263  if ( ! ProjectileNucleon->AreYouHit() ) {
264  // Projectile nucleon was not involved until now.
265  ProjectileSplitable = new G4DiffractiveSplitableHadron( *ProjectileNucleon );
266  ProjectileNucleon->Hit( ProjectileSplitable );
267  ProjectileSplitable->SetStatus( 1 ); // It takes part in the interaction
268  } else { // Projectile nucleon was involved before.
269  ProjectileSplitable = ProjectileNucleon->GetSplitableHadron();
270  }
271 
272  if ( ! TargetNucleon->AreYouHit() ) { // Target nucleon was not involved until now
273  TargetSplitable = new G4DiffractiveSplitableHadron( *TargetNucleon );
274  TargetNucleon->Hit( TargetSplitable );
275  TargetSplitable->SetStatus( 1 ); // It takes part in the interaction
276  } else { // Target nucleon was involved before.
277  TargetSplitable = TargetNucleon->GetSplitableHadron();
278  }
279 
280  G4InteractionContent* anInteraction = new G4InteractionContent( ProjectileSplitable );
281  anInteraction->SetTarget( TargetSplitable );
282  anInteraction->SetProjectileNucleon( ProjectileNucleon );
283  anInteraction->SetTargetNucleon( TargetNucleon );
284  anInteraction->SetInteractionTime( ( ProjectileNucleon->GetPosition().z() +
285  TargetNucleon->GetPosition().z() ) / betta_z );
286  anInteraction->SetStatus( 1 );
287 
288  #ifdef debugFTFparticipant
289  G4cout << "Part anInteraction->GetInteractionTime() "
290  << anInteraction->GetInteractionTime()/fermi << G4endl
291  << "Splitable Pr* Tr* " << ProjectileSplitable << " "
292  << TargetSplitable << G4endl;
293  #endif
294 
295  theInteractions.push_back( anInteraction );
296 
297  } // End of an Interaction has happend!
298 
299  #ifdef debugFTFparticipant
300  TrNuclN++;
301  #endif
302 
303  } // End of while ( ( TargetNucleon = theNucleus->GetNextNucleon() ) )
304 
305  #ifdef debugFTFparticipant
306  PrNuclN++;
307  #endif
308 
309  } // End of while ( ( ProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) )
310 
311  if ( theInteractions.size() != 0 ) theProjectileNucleus->DoTranslation( theBeamPosition );
312 
313  } while ( ( theInteractions.size() == 0 ) &&
314  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
315  if ( loopCounter >= maxNumberOfLoops ) {
316  #ifdef debugFTFparticipant
317  G4cout << "BAD situation: forced exit from the while loop!" << G4endl;
318  #endif
319  return;
320  }
321 
324 
325  #ifdef debugFTFparticipant
326  G4cout << G4endl << "Number of primary collisions " << theInteractions.size()
327  << "\t Bx " << impactX/fermi << "\t By " << impactY/fermi
328  << "\t B " << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl
329  << "FTF participant End. #######################" << G4endl << G4endl;
330  #endif
331  return;
332 }
333 
334 
335 //============================================================================
336 
338  const G4InteractionContent* Int2 ) {
339  return Int1->GetInteractionTime() < Int2->GetInteractionTime();
340 }
341 
342 
343 //============================================================================
344 
345 void G4FTFParticipants::SortInteractionsIncT() { // on increased T
346  if ( theInteractions.size() < 2 ) return; // Avoid unnecesary work
347  std::sort( theInteractions.begin(), theInteractions.end(), G4FTFPartHelperForSortInT );
348 }
349 
350 
351 //============================================================================
352 
354  G4double InitialTime = theInteractions[0]->GetInteractionTime();
355  for ( unsigned int i = 1; i < theInteractions.size(); i++ ) {
356  G4double InterTime = theInteractions[i]->GetInteractionTime() - InitialTime;
357  theInteractions[i]->SetInteractionTime( InterTime );
358  G4InteractionContent* aCollision = theInteractions[i];
359  G4VSplitableHadron* projectile = aCollision->GetProjectile();
360  G4VSplitableHadron* target = aCollision->GetTarget();
361  G4ThreeVector prPosition = projectile->GetPosition();
362  prPosition.setZ( target->GetPosition().z() );
363  projectile->SetPosition( prPosition );
364  projectile->SetTimeOfCreation( InterTime );
365  target->SetTimeOfCreation( InterTime );
366  }
367  return;
368 }
369 
370 
371 //============================================================================
372 
374  for ( size_t i = 0; i < theInteractions.size(); i++ ) {
375  if ( theInteractions[ i ] ) {
376  delete theInteractions[ i ];
377  theInteractions[ i ] = 0;
378  }
379  }
380  theInteractions.clear();
381  currentInteraction = -1;
382 }
383 
void GetList(const G4ReactionProduct &thePrimary, G4FTFParameters *theParameters)
CLHEP::Hep3Vector G4ThreeVector
std::vector< G4InteractionContent * > theInteractions
virtual G4bool StartLoop()=0
void SetInteractionTime(G4double aValue)
void SetTimeOfCreation(G4double aTime)
virtual const G4LorentzVector & Get4Momentum() const
Definition: G4Nucleon.hh:72
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
int G4int
Definition: G4Types.hh:78
void SetStatus(const G4int aStatus)
G4VSplitableHadron * GetSplitableHadron() const
Definition: G4Nucleon.hh:96
virtual G4double GetOuterRadius()=0
void SetProjectileNucleon(G4Nucleon *aNucleon)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4bool nucleon(G4int ityp)
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
G4V3DNucleus * theProjectileNucleus
G4V3DNucleus * theNucleus
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
Definition: G4V3DNucleus.hh:87
G4double GetTotalEnergy() const
void SetStatus(G4int aValue)
void SetPosition(const G4ThreeVector &aPosition)
G4VSplitableHadron * GetTarget() const
G4VSplitableHadron * GetProjectile() const
G4ThreeVector GetMomentum() const
const G4ThreeVector & GetPosition() const
#define G4endl
Definition: G4ios.hh:61
G4double GetProbabilityOfInteraction(const G4double impactsquare)
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
G4double GetInteractionTime() const
virtual void DoTranslation(const G4ThreeVector &theShift)=0
bool G4FTFPartHelperForSortInT(const G4InteractionContent *Int1, const G4InteractionContent *Int2)
static const double fermi
Definition: G4SIunits.hh:102
void SetTargetNucleon(G4Nucleon *aNucleon)
void SetTarget(G4VSplitableHadron *aTarget)