Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 #include "G4PhysicalConstants.hh"
 53 
 54 
 55 //============================================================================
 56 
 57 //#define debugFTFparticipant
 58 
 59 
 60 //============================================================================
 61 
 62 G4FTFParticipants::G4FTFParticipants() : Bimpact( 0.0 ), BinInterval( false ),
 63                                          Bmin2( -1.0 ), Bmax2( -1.0 ),
 64                                          currentInteraction( -1 )
 65 {}
 66 
 67 
 68 //============================================================================
 69 
 70 G4FTFParticipants::~G4FTFParticipants() {}
 71 
 72 //============================================================================
 73 
 74 void G4FTFParticipants::GetList( const G4ReactionProduct& thePrimary, 
 75                                  G4FTFParameters* theParameters ) { 
 76 
 77   #ifdef debugFTFparticipant
 78   G4cout << "Participants::GetList" << G4endl 
 79          << "thePrimary " << thePrimary.GetMomentum() << G4endl << G4endl;
 80   #endif
 81 
 82   G4double betta_z = thePrimary.GetMomentum().z() / thePrimary.GetTotalEnergy();
 83   if ( betta_z < 1.0e-10 ) betta_z = 1.0e-10;
 84 
 85   StartLoop(); // reset Loop over Interactions
 86 
 87   for ( unsigned int i = 0; i < theInteractions.size(); i++ ) delete theInteractions[i];
 88   theInteractions.clear();
 89 
 90   G4double deltaxy = 2.0 * fermi;  // Extra nuclear radius
 91 
 92   if ( theProjectileNucleus == nullptr ) {  // Hadron-nucleus or anti-baryon-nucleus interactions
 93 
 94     G4double impactX( 0.0 ), impactY( 0.0 );
 95     G4double B( 0.0 ), B2( 0.0 );
 96 
 97     G4VSplitableHadron* primarySplitable = new G4DiffractiveSplitableHadron( thePrimary );
 98 
 99     #ifdef debugFTFparticipant
100     G4cout << "Hadron-nucleus or anti-baryon-nucleus interactions" << G4endl;
101     #endif
102 
103     G4double xyradius = theNucleus->GetOuterRadius() + deltaxy;  // Range of impact parameter sampling
104 
105     const G4int maxNumberOfLoops = 1000;
106     G4int loopCounter = 0;                                          
107     do {
108 
109       std::pair< G4double, G4double > theImpactParameter;
110       if ( SampleBinInterval() ) {
111         B2 = GetBmin2() + G4UniformRand() * ( GetBmax2() - GetBmin2() );
112         B  = B2 > 0.0 ? std::sqrt( B2 ) : 0.0;
113         G4double Phi = twopi * G4UniformRand();
114         impactX = B * std::cos( Phi );
115         impactY = B * std::sin( Phi );
116         SetImpactParameter( B );
117       } else {
118         theImpactParameter = theNucleus->ChooseImpactXandY( xyradius );
119         impactX = theImpactParameter.first; 
120         impactY = theImpactParameter.second;
121         SetImpactParameter( std::sqrt( sqr(impactX) + sqr(impactY) ) );
122       }
123 
124       #ifdef debugFTFparticipant
125       G4cout << "New interaction list," << " b[fm]= " 
126              << std::sqrt( sqr(impactX ) + sqr( impactY ) )/fermi << G4endl;
127       #endif
128 
129       G4ThreeVector thePosition( impactX, impactY, 0.0 );
130       primarySplitable->SetPosition( thePosition );
131 
132       theNucleus->StartLoop();
133       G4Nucleon* nucleon;
134 
135       #ifdef debugFTFparticipant
136       G4int TrN( 0 );
137       #endif
138 
139       while ( ( nucleon = theNucleus->GetNextNucleon() ) ) {  /* Loop checking, 10.08.2015, A.Ribon */
140 
141         G4double impact2 = sqr( impactX - nucleon->GetPosition().x() ) +
142                            sqr( impactY - nucleon->GetPosition().y() );
143 
144         if ( theParameters->GetProbabilityOfInteraction( impact2/fermi/fermi ) > 
145              G4UniformRand() ) {
146           primarySplitable->SetStatus( 1 );  // It takes part in the interaction
147           G4VSplitableHadron* targetSplitable = 0;
148           if ( ! nucleon->AreYouHit() ) {
149             targetSplitable = new G4DiffractiveSplitableHadron( *nucleon );
150             nucleon->Hit( targetSplitable );
151             targetSplitable->SetStatus( 1 ); // It takes part in the interaction
152 
153             #ifdef debugFTFparticipant
154             G4cout << "Participated nucleons #, " << TrN << " " << "Splitable Pr* Tr* "
155                    << primarySplitable << " " << targetSplitable << G4endl;
156             #endif
157 
158           }
159           G4InteractionContent* aInteraction = new G4InteractionContent( primarySplitable );
160           G4Nucleon* PrNucleon = 0;
161           aInteraction->SetProjectileNucleon( PrNucleon );
162           aInteraction->SetTarget( targetSplitable );
163           aInteraction->SetTargetNucleon( nucleon );     
164           aInteraction->SetStatus( 1 );                  
165           aInteraction->SetInteractionTime( ( primarySplitable->GetPosition().z() + 
166                                               nucleon->GetPosition().z() ) / betta_z );
167           theInteractions.push_back( aInteraction );
168         }
169 
170         #ifdef debugFTFparticipant
171         TrN++;
172         #endif
173 
174       } 
175 
176     } while ( ( theInteractions.size() == 0 ) && 
177               ++loopCounter < maxNumberOfLoops );  /* Loop checking, 10.08.2015, A.Ribon */
178     if ( loopCounter >= maxNumberOfLoops ) {
179       #ifdef debugFTFparticipant
180       G4cout << "BAD situation: forced exit from the while loop!" << G4endl;
181       #endif
182       return;
183     }
184 
185     #ifdef debugFTFparticipant
186     G4cout << "Number of Hit nucleons " << theInteractions.size() << "\t Bx[fm] " << impactX/fermi
187            << "\t By[fm] " << impactY/fermi << "\t B[fm] " 
188            << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl << G4endl;
189     #endif
190 
191     //SortInteractionsIncT(); // Not need because nucleons are sorted in increasing z-coordinates.
192     ShiftInteractionTime();  // To put correct times and z-coordinates
193     return;
194 
195   }  // end of if ( theProjectileNucleus == 0 )
196 
197   // Projectile and target are nuclei
198 
199   #ifdef debugFTFparticipant
200   G4cout << "Projectile and target are nuclei" << G4endl;
201   #endif
202 
203   //G4cout<<theProjectileNucleus->GetOuterRadius()/fermi<<" "<<theNucleus->GetOuterRadius()/fermi<<" "<<deltaxy/fermi<<G4endl;
204 
205   // Range of impact parameter sampling
206   G4double xyradius = theProjectileNucleus->GetOuterRadius() + theNucleus->GetOuterRadius() + deltaxy;
207 
208   G4double impactX( 0.0 ), impactY( 0.0 );
209   G4double B( 0.0 ), B2( 0.0 );
210 
211   const G4int maxNumberOfLoops = 1000;
212   G4int loopCounter = 0;
213   do {
214 
215     std::pair< G4double, G4double > theImpactParameter;
216     if ( SampleBinInterval() ) {
217       B2 = GetBmin2() + G4UniformRand() * ( GetBmax2() - GetBmin2() );
218       B  = B2 > 0.0 ? std::sqrt( B2 ) : 0.0;  // In G4 internal units (mm)
219       G4double Phi = twopi * G4UniformRand();
220       impactX = B * std::cos( Phi ); 
221       impactY = B * std::sin( Phi );
222       SetImpactParameter( B );
223     } else {
224       theImpactParameter = theNucleus->ChooseImpactXandY( xyradius );
225       impactX = theImpactParameter.first; 
226       impactY = theImpactParameter.second;
227       SetImpactParameter( std::sqrt( sqr(impactX) + sqr(impactY) ) );
228     }
229 
230     #ifdef debugFTFparticipant
231     G4cout << "New interaction list, " << "b[fm] "
232            << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl;
233     #endif
234 
235     G4ThreeVector theBeamPosition( impactX, impactY, 0.0 );                    
236 
237     theProjectileNucleus->StartLoop();
238     G4Nucleon* ProjectileNucleon;
239 
240     #ifdef debugFTFparticipant
241     G4int PrNuclN( 0 );
242     #endif
243 
244     while ( ( ProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) ) {  /* Loop checking, 10.08.2015, A.Ribon */
245  
246       G4VSplitableHadron* ProjectileSplitable = 0;
247       theNucleus->StartLoop();
248       G4Nucleon* TargetNucleon = 0;
249 
250       #ifdef debugFTFparticipant
251       G4int TrNuclN( 0 );
252       #endif
253 
254       while ( ( TargetNucleon = theNucleus->GetNextNucleon() ) ) {  /* Loop checking, 10.08.2015, A.Ribon */
255 
256         G4double impact2 = sqr( impactX + ProjectileNucleon->GetPosition().x() - 
257                                 TargetNucleon->GetPosition().x() ) +
258                            sqr( impactY + ProjectileNucleon->GetPosition().y() -
259                                 TargetNucleon->GetPosition().y() );
260         G4VSplitableHadron* TargetSplitable = 0;
261         if ( theParameters->GetProbabilityOfInteraction( impact2/fermi/fermi ) >
262              G4UniformRand() ) {  // An interaction has happend!
263 
264           #ifdef debugFTFparticipant
265           G4cout << G4endl << "An Interaction has happend" << G4endl << "Proj N mom " << PrNuclN
266                  << " " << ProjectileNucleon->Get4Momentum() << "-------------" << G4endl
267                  << "Targ N mom " << TrNuclN << " " << TargetNucleon->Get4Momentum() << G4endl
268                  << "PrN TrN Z coords [fm]" << ProjectileNucleon->GetPosition().z()/fermi 
269                  << " " << TargetNucleon->GetPosition().z()/fermi 
270                  << " " << ProjectileNucleon->GetPosition().z()/fermi + 
271                            TargetNucleon->GetPosition().z()/fermi << G4endl;
272           #endif
273 
274           if ( ! ProjectileNucleon->AreYouHit() ) { 
275             // Projectile nucleon was not involved until now.
276             ProjectileSplitable = new G4DiffractiveSplitableHadron( *ProjectileNucleon );
277             ProjectileNucleon->Hit( ProjectileSplitable );
278             ProjectileSplitable->SetStatus( 1 );  // It takes part in the interaction
279           } else {  // Projectile nucleon was involved before.
280             ProjectileSplitable = ProjectileNucleon->GetSplitableHadron();
281           }
282 
283           if ( ! TargetNucleon->AreYouHit() ) {  // Target nucleon was not involved until now
284             TargetSplitable = new G4DiffractiveSplitableHadron( *TargetNucleon );
285             TargetNucleon->Hit( TargetSplitable );
286             TargetSplitable->SetStatus( 1 );   // It takes part in the interaction
287           } else {  // Target nucleon was involved before.
288             TargetSplitable = TargetNucleon->GetSplitableHadron();
289           }
290 
291           G4InteractionContent* anInteraction = new G4InteractionContent( ProjectileSplitable );
292           anInteraction->SetTarget( TargetSplitable );
293           anInteraction->SetProjectileNucleon( ProjectileNucleon );
294           anInteraction->SetTargetNucleon( TargetNucleon );
295           anInteraction->SetInteractionTime( ( ProjectileNucleon->GetPosition().z() + 
296                                                TargetNucleon->GetPosition().z() ) / betta_z );
297           anInteraction->SetStatus( 1 );                     
298 
299           #ifdef debugFTFparticipant
300           G4cout << "Part anInteraction->GetInteractionTime() [fm] " 
301                  << anInteraction->GetInteractionTime()/fermi << G4endl
302                  << "Splitable Pr* Tr* " << ProjectileSplitable << " " 
303                  << TargetSplitable << G4endl;
304           #endif
305 
306           theInteractions.push_back( anInteraction );
307 
308         } // End of an Interaction has happend!
309 
310         #ifdef debugFTFparticipant
311         TrNuclN++;
312         #endif
313 
314       }  // End of while ( ( TargetNucleon = theNucleus->GetNextNucleon() ) )
315 
316       #ifdef debugFTFparticipant
317       PrNuclN++;
318       #endif
319 
320     }  // End of while ( ( ProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) )
321 
322     if ( theInteractions.size() != 0 ) theProjectileNucleus->DoTranslation( theBeamPosition );
323 
324   } while ( ( theInteractions.size() == 0 ) &&
325             ++loopCounter < maxNumberOfLoops );  /* Loop checking, 10.08.2015, A.Ribon */
326   if ( loopCounter >= maxNumberOfLoops ) {
327     #ifdef debugFTFparticipant
328     G4cout << "BAD situation: forced exit from the while loop!" << G4endl;
329     #endif
330     return;
331   }
332 
333   SortInteractionsIncT();
334   ShiftInteractionTime();
335 
336   #ifdef debugFTFparticipant
337   G4cout << G4endl << "Number of primary collisions " << theInteractions.size() 
338          << "\t Bx[fm] " << impactX/fermi << "\t By[fm] " << impactY/fermi
339          << "\t B[fm] " << std::sqrt( sqr( impactX ) + sqr( impactY ) )/fermi << G4endl
340          << "FTF participant End. #######################" << G4endl << G4endl;
341   #endif
342   return;
343 }
344 
345 
346 //============================================================================
347 
348 bool G4FTFPartHelperForSortInT( const G4InteractionContent* Int1, 
349                                 const G4InteractionContent* Int2 ) {
350   return Int1->GetInteractionTime() < Int2->GetInteractionTime();
351 }
352 
353 
354 //============================================================================
355 
356 void G4FTFParticipants::SortInteractionsIncT() {  // on increased T 
357   if ( theInteractions.size() < 2 ) return;  // Avoid unnecesary work
358   std::sort( theInteractions.begin(), theInteractions.end(), G4FTFPartHelperForSortInT ); 
359 }
360 
361 
362 //============================================================================
363 
364 void G4FTFParticipants::ShiftInteractionTime() {
365   G4double InitialTime = theInteractions[0]->GetInteractionTime();
366   for ( unsigned int i = 1; i < theInteractions.size(); i++ ) {
367     G4double InterTime = theInteractions[i]->GetInteractionTime() - InitialTime;
368     theInteractions[i]->SetInteractionTime( InterTime );
369     G4InteractionContent* aCollision = theInteractions[i];
370     G4VSplitableHadron* projectile = aCollision->GetProjectile();
371     G4VSplitableHadron* target = aCollision->GetTarget();
372     G4ThreeVector prPosition = projectile->GetPosition();
373     prPosition.setZ( target->GetPosition().z() );
374     projectile->SetPosition( prPosition );
375     projectile->SetTimeOfCreation( InterTime );
376     target->SetTimeOfCreation( InterTime );
377   }
378   return;
379 }
380 
381 
382 //============================================================================
383 
384 void G4FTFParticipants::Clean() {
385   for ( size_t i = 0; i < theInteractions.size(); i++ ) {
386     if ( theInteractions[ i ] ) {
387       delete theInteractions[ i ];
388       theInteractions[ i ] = 0;
389     }
390   }
391   theInteractions.clear();
392   currentInteraction = -1;
393 }
394 
395