ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ElasticHNScattering.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ElasticHNScattering.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 implemetation file
31 //
32 // ---------------- G4ElasticHNScattering --------------
33 // by V. Uzhinsky, March 2008.
34 // elastic scattering used by Fritiof model
35 // Take a projectile and a target
36 // scatter the projectile and target
37 // ---------------------------------------------------------------------
38 
39 #include "globals.hh"
40 #include "Randomize.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 #include "G4ElasticHNScattering.hh"
45 #include "G4LorentzRotation.hh"
46 #include "G4ThreeVector.hh"
47 #include "G4ParticleDefinition.hh"
48 #include "G4VSplitableHadron.hh"
49 #include "G4ExcitedString.hh"
50 #include "G4FTFParameters.hh"
51 
52 #include "G4SampleResonance.hh"
53 
54 #include "G4Exp.hh"
55 #include "G4Log.hh"
56 
57 //============================================================================
58 
60 
61 
62 //============================================================================
63 
66  G4FTFParameters* theParameters ) const {
67  projectile->IncrementCollisionCount( 1 );
68  target->IncrementCollisionCount( 1 );
69 
70  if ( projectile->Get4Momentum().z() < 0.0 ) return false; //Uzhi Aug.2019
71 
72  // Projectile parameters
73  G4LorentzVector Pprojectile = projectile->Get4Momentum();
74  G4double M0projectile = Pprojectile.mag();
75  G4double M0projectile2 = M0projectile * M0projectile;
76 
77  // Target parameters
78  G4LorentzVector Ptarget = target->Get4Momentum();
79  G4double M0target = Ptarget.mag();
80  G4double M0target2 = M0target * M0target;
81 
82  G4double AveragePt2 = theParameters->GetAvaragePt2ofElasticScattering();
83 
84  // Transform momenta to cms and then rotate parallel to z axis;
85  G4LorentzVector Psum;
86  Psum = Pprojectile + Ptarget;
87  G4LorentzRotation toCms( -1*Psum.boostVector() );
88  G4LorentzVector Ptmp = toCms*Pprojectile;
89  if ( Ptmp.pz() <= 0.0 ) return false;
90  //"String" moving backwards in CMS, abort collision !
91  //G4cout << " abort Collision! " << G4endl;
92  toCms.rotateZ( -1*Ptmp.phi() );
93  toCms.rotateY( -1*Ptmp.theta() );
94  G4LorentzRotation toLab( toCms.inverse() );
95  Pprojectile.transform( toCms );
96  Ptarget.transform( toCms );
97 
98  G4double PZcms2, PZcms;
99  G4double S = Psum.mag2();
100  G4double SqrtS = std::sqrt( S );
101  if ( SqrtS < M0projectile + M0target ) return false;
102 
103  PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( M0target2 )
104  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
105 
106  PZcms = ( PZcms2 > 0.0 ? std::sqrt( PZcms2 ) : 0.0 );
107 
108  G4double maxPtSquare = PZcms2;
109 
110  // Now we can calculate the transferred Pt
111  G4double Pt2;
112  G4double ProjMassT2, ProjMassT;
113  G4double TargMassT2, TargMassT;
114  G4LorentzVector Qmomentum;
115 
116  const G4int maxNumberOfLoops = 1000;
117  G4int loopCounter = 0;
118  do {
119  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0.0 );
120  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
121  ProjMassT2 = M0projectile2 + Pt2;
122  ProjMassT = std::sqrt( ProjMassT2 );
123  TargMassT2 = M0target2 + Pt2;
124  TargMassT = std::sqrt( TargMassT2 );
125  } while ( ( SqrtS < ProjMassT + TargMassT ) &&
126  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
127  if ( loopCounter >= maxNumberOfLoops ) {
128  return false;
129  }
130 
131  PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( TargMassT2 )
132  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
133 
134  if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // to avoid the exactness problem
135  PZcms = std::sqrt( PZcms2 );
136  Pprojectile.setPz( PZcms );
137  Ptarget.setPz( -PZcms );
138  Pprojectile += Qmomentum;
139  Ptarget -= Qmomentum;
140 
141  // Transform back and update SplitableHadron Participant.
142  Pprojectile.transform( toLab );
143  Ptarget.transform( toLab );
144 
145  // Calculation of the creation time
146  projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
147  projectile->SetPosition( target->GetPosition() );
148 
149  // Creation time and position of target nucleon were determined at
150  // ReggeonCascade() of G4FTFModel
151 
152  projectile->Set4Momentum( Pprojectile );
153  target->Set4Momentum( Ptarget );
154 
155  //projectile->IncrementCollisionCount( 1 );
156  //target->IncrementCollisionCount( 1 );
157 
158  return true;
159 }
160 
161 
162 //============================================================================
163 
165  G4double maxPtSquare ) const {
166  // @@ this method is used in FTFModel as well. Should go somewhere common!
167  G4double Pt2( 0.0 );
168  if ( AveragePt2 <= 0.0 ) {
169  Pt2 = 0.0;
170  } else {
171  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() * ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
172  }
173  G4double Pt = ( Pt2 > 0.0 ? std::sqrt( Pt2 ) : 0.0 );
175  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
176 }
177 
178 
179 //============================================================================
180 
182  throw G4HadronicException( __FILE__, __LINE__,
183  "G4ElasticHNScattering copy constructor not meant to be called" );
184 }
185 
186 
187 //============================================================================
188 
190 
191 
192 //============================================================================
193 
195  throw G4HadronicException( __FILE__, __LINE__,
196  "G4ElasticHNScattering = operator not meant to be called" );
197 }
198 
199 
200 //============================================================================
201 
203  throw G4HadronicException( __FILE__, __LINE__,
204  "G4ElasticHNScattering == operator not meant to be called" );
205 }
206 
207 
208 //============================================================================
209 
211  throw G4HadronicException( __FILE__, __LINE__,
212  "G4ElasticHNScattering != operator not meant to be called" );
213 }
214