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