Geant4  10.00.p03
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 74627 2013-10-17 07:04:38Z 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  do { // while ( theInteractions.size() == 0 )
115 
116  std::pair< G4double, G4double > theImpactParameter;
117  theImpactParameter = theNucleus->ChooseImpactXandY( xyradius );
118  impactX = theImpactParameter.first;
119  impactY = theImpactParameter.second;
120 
121  #ifdef debugFTFparticipant
122  G4cout << "New interaction list," << " b= "
123  << std::sqrt( sqr(impactX ) + sqr( impactY ) )/fermi << G4endl;
124  #endif
125 
126  G4ThreeVector thePosition( impactX, impactY, 0.0 );
127  primarySplitable->SetPosition( thePosition );
128 
131 
132  #ifdef debugFTFparticipant
133  G4int TrN( 0 );
134  #endif
135 
136  while ( ( nucleon = theNucleus->GetNextNucleon() ) ) {
137 
138  G4double impact2 = sqr( impactX - nucleon->GetPosition().x() ) +
139  sqr( impactY - nucleon->GetPosition().y() );
140 
141  if ( theParameters->GetProbabilityOfInteraction( impact2/fermi/fermi ) >
142  G4UniformRand() ) {
143  primarySplitable->SetStatus( 1 ); // It takes part in the interaction
144  G4VSplitableHadron* targetSplitable = 0;
145  if ( ! nucleon->AreYouHit() ) {
146  targetSplitable = new G4DiffractiveSplitableHadron( *nucleon );
147  nucleon->Hit( targetSplitable );
148  targetSplitable->SetStatus( 1 ); // It takes part in the interaction
149 
150  #ifdef debugFTFparticipant
151  G4cout << "Participated nucleons #, " << TrN << " " << "Splitable Pr* Tr* "
152  << primarySplitable << " " << targetSplitable << G4endl;
153  #endif
154 
155  }
156  G4InteractionContent* aInteraction = new G4InteractionContent( primarySplitable );
157  G4Nucleon* PrNucleon = 0;
158  aInteraction->SetProjectileNucleon( PrNucleon );
159  aInteraction->SetTarget( targetSplitable );
160  aInteraction->SetTargetNucleon( nucleon );
161  aInteraction->SetStatus( 1 );
162  aInteraction->SetInteractionTime( ( primarySplitable->GetPosition().z() +
163  nucleon->GetPosition().z() ) / betta_z );
164  theInteractions.push_back( aInteraction );
165  }
166 
167  #ifdef debugFTFparticipant
168  TrN++;
169  #endif
170 
171  }
172 
173  } while ( theInteractions.size() == 0 );
174 
175  #ifdef debugFTFparticipant
176  G4cout << "Number of Hit nucleons " << theInteractions.size() << "\t Bx " << impactX/fermi
177  << "\t By " << impactY/fermi << "\t B "
178  << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl << G4endl;
179  #endif
180 
181  //SortInteractionsIncT(); // Not need because nucleons are sorted in increasing z-coordinates.
182  ShiftInteractionTime(); // To put correct times and z-coordinates
183  return;
184 
185  } // end of if ( theProjectileNucleus == 0 )
186 
187  // Projectile and target are nuclei
188 
189  #ifdef debugFTFparticipant
190  G4cout << "Projectile and target are nuclei" << G4endl;
191  #endif
192 
193  G4double xyradius;
194  xyradius = theProjectileNucleus->GetOuterRadius() + // Range of impact parameter sampling
195  theNucleus->GetOuterRadius() + deltaxy;
196 
197  G4double impactX( 0.0 ), impactY( 0.0 );
198 
199  do { // while ( theInteractions.size() == 0 )
200 
201  std::pair< G4double, G4double > theImpactParameter;
202  theImpactParameter = theNucleus->ChooseImpactXandY( xyradius );
203  impactX = theImpactParameter.first;
204  impactY = theImpactParameter.second;
205 
206  #ifdef debugFTFparticipant
207  G4cout << "New interaction list, " << "b "
208  << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl;
209  #endif
210 
211  G4ThreeVector theBeamPosition( impactX, impactY, 0.0 );
212 
214  G4Nucleon* ProjectileNucleon;
215 
216  #ifdef debugFTFparticipant
217  G4int PrNuclN( 0 );
218  #endif
219 
220  while ( ( ProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) ) {
221 
222  G4VSplitableHadron* ProjectileSplitable = 0;
224  G4Nucleon* TargetNucleon = 0;
225 
226  #ifdef debugFTFparticipant
227  G4int TrNuclN( 0 );
228  #endif
229 
230  while ( ( TargetNucleon = theNucleus->GetNextNucleon() ) ) {
231 
232  G4double impact2 = sqr( impactX + ProjectileNucleon->GetPosition().x() -
233  TargetNucleon->GetPosition().x() ) +
234  sqr( impactY + ProjectileNucleon->GetPosition().y() -
235  TargetNucleon->GetPosition().y() );
236  G4VSplitableHadron* TargetSplitable = 0;
237  if ( theParameters->GetProbabilityOfInteraction( impact2/fermi/fermi ) >
238  G4UniformRand() ) { // An interaction has happend!
239 
240  #ifdef debugFTFparticipant
241  G4cout << G4endl << "An Interaction has happend" << G4endl << "Proj N mom " << PrNuclN
242  << " " << ProjectileNucleon->Get4Momentum() << "-------------" << G4endl
243  << "Targ N mom " << TrNuclN << " " << TargetNucleon->Get4Momentum() << G4endl
244  << "PrN TrN Z coords " << ProjectileNucleon->GetPosition().z()/fermi
245  << " " << TargetNucleon->GetPosition().z()/fermi
246  << " " << ProjectileNucleon->GetPosition().z()/fermi +
247  TargetNucleon->GetPosition().z()/fermi << G4endl;
248  #endif
249 
250  if ( ! ProjectileNucleon->AreYouHit() ) {
251  // Projectile nucleon was not involved until now.
252  ProjectileSplitable = new G4DiffractiveSplitableHadron( *ProjectileNucleon );
253  ProjectileNucleon->Hit( ProjectileSplitable );
254  ProjectileSplitable->SetStatus( 1 ); // It takes part in the interaction
255  } else { // Projectile nucleon was involved before.
256  ProjectileSplitable = ProjectileNucleon->GetSplitableHadron();
257  }
258 
259  if ( ! TargetNucleon->AreYouHit() ) { // Target nucleon was not involved until now
260  TargetSplitable = new G4DiffractiveSplitableHadron( *TargetNucleon );
261  TargetNucleon->Hit( TargetSplitable );
262  TargetSplitable->SetStatus( 1 ); // It takes part in the interaction
263  } else { // Target nucleon was involved before.
264  TargetSplitable = TargetNucleon->GetSplitableHadron();
265  }
266 
267  G4InteractionContent* anInteraction = new G4InteractionContent( ProjectileSplitable );
268  anInteraction->SetTarget( TargetSplitable );
269  anInteraction->SetProjectileNucleon( ProjectileNucleon );
270  anInteraction->SetTargetNucleon( TargetNucleon );
271  anInteraction->SetInteractionTime( ( ProjectileNucleon->GetPosition().z() +
272  TargetNucleon->GetPosition().z() ) / betta_z );
273  anInteraction->SetStatus( 1 );
274 
275  #ifdef debugFTFparticipant
276  G4cout << "Part anInteraction->GetInteractionTime() "
277  << anInteraction->GetInteractionTime()/fermi << G4endl
278  << "Splitable Pr* Tr* " << ProjectileSplitable << " "
279  << TargetSplitable << G4endl;
280  #endif
281 
282  theInteractions.push_back( anInteraction );
283 
284  } // End of an Interaction has happend!
285 
286  #ifdef debugFTFparticipant
287  TrNuclN++;
288  #endif
289 
290  } // End of while ( ( TargetNucleon = theNucleus->GetNextNucleon() ) )
291 
292  #ifdef debugFTFparticipant
293  PrNuclN++;
294  #endif
295 
296  } // End of while ( ( ProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) )
297 
298  if ( theInteractions.size() != 0 ) theProjectileNucleus->DoTranslation( theBeamPosition );
299 
300  } while ( theInteractions.size() == 0 );
301 
304 
305  #ifdef debugFTFparticipant
306  G4cout << G4endl << "Number of primary collisions " << theInteractions.size()
307  << "\t Bx " << impactX/fermi << "\t By " << impactY/fermi
308  << "\t B " << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl
309  << "FTF participant End. #######################" << G4endl << G4endl;
310  #endif
311 
312  return;
313 }
314 
315 
316 //============================================================================
317 
319  const G4InteractionContent* Int2 ) {
320  return Int1->GetInteractionTime() < Int2->GetInteractionTime();
321 }
322 
323 
324 //============================================================================
325 
326 void G4FTFParticipants::SortInteractionsIncT() { // on increased T
327  if ( theInteractions.size() < 2 ) return; // Avoid unnecesary work
328  std::sort( theInteractions.begin(), theInteractions.end(), G4FTFPartHelperForSortInT );
329 }
330 
331 
332 //============================================================================
333 
335  G4double InitialTime = theInteractions[0]->GetInteractionTime();
336  for ( unsigned int i = 1; i < theInteractions.size(); i++ ) {
337  G4double InterTime = theInteractions[i]->GetInteractionTime() - InitialTime;
338  theInteractions[i]->SetInteractionTime( InterTime );
339  G4InteractionContent* aCollision = theInteractions[i];
340  G4VSplitableHadron* projectile = aCollision->GetProjectile();
341  G4VSplitableHadron* target = aCollision->GetTarget();
342  G4ThreeVector prPosition = projectile->GetPosition();
343  prPosition.setZ( target->GetPosition().z() );
344  projectile->SetPosition( prPosition );
345  projectile->SetTimeOfCreation( InterTime );
346  target->SetTimeOfCreation( InterTime );
347  }
348  return;
349 }
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:87
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:93
void SetTargetNucleon(G4Nucleon *aNucleon)
void SetTarget(G4VSplitableHadron *aTarget)