Geant4  10.01.p03
G4ElasticHNScattering.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: G4ElasticHNScattering.cc 86646 2014-11-14 13:29:39Z gcosmo $
28 //
29 
30 // ------------------------------------------------------------
31 // GEANT 4 class implemetation file
32 //
33 // ---------------- G4ElasticHNScattering --------------
34 // by V. Uzhinsky, March 2008.
35 // elastic scattering used by Fritiof model
36 // Take a projectile and a target
37 // scatter the projectile and target
38 // ---------------------------------------------------------------------
39 
40 #include "globals.hh"
41 #include "Randomize.hh"
42 #include "G4PhysicalConstants.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" // Uzhi Oct 2014
53 
54 //============================================================================
55 
57 
58 
59 //============================================================================
60 
62  G4VSplitableHadron* target,
63  G4FTFParameters* theParameters ) const {
64  projectile->IncrementCollisionCount( 1 );
65  target->IncrementCollisionCount( 1 );
66 
67  G4SampleResonance BrW; // Uzhi Oct 2014
68 
69  // Projectile parameters
70  G4LorentzVector Pprojectile = projectile->Get4Momentum();
71  if ( Pprojectile.z() < 0.0 ) return false;
72  G4bool PutOnMassShell( false );
73  G4double M0projectile = Pprojectile.mag();
74 // if ( M0projectile < projectile->GetDefinition()->GetPDGMass() ) { // Uzhi Oct 2014
75 
76  G4double MminProjectile=BrW.GetMinimumMass(projectile->GetDefinition()); // Uzhi Oct 2014
77 
78  if ( M0projectile < MminProjectile ) { // Uzhi Oct 2014
79  PutOnMassShell = true;
80  M0projectile = projectile->GetDefinition()->GetPDGMass();
81  }
82  G4double M0projectile2 = M0projectile * M0projectile;
83  G4double AveragePt2 = theParameters->GetAvaragePt2ofElasticScattering();
84 
85  // Target parameters
86  G4LorentzVector Ptarget = target->Get4Momentum();
87  G4double M0target = Ptarget.mag();
88 // if ( M0target < target->GetDefinition()->GetPDGMass() ) { // Uzhi Oct 2014
89 
90  G4double MminTarget=BrW.GetMinimumMass(target->GetDefinition()); // Uzhi Oct 2014
91 
92  if ( M0target < MminTarget ) { // Uzhi Oct 2014
93  PutOnMassShell = true;
94  M0target = target->GetDefinition()->GetPDGMass();
95  }
96  G4double M0target2 = M0target * M0target;
97 
98  // Transform momenta to cms and then rotate parallel to z axis;
99  G4LorentzVector Psum;
100  Psum = Pprojectile + Ptarget;
101  G4LorentzRotation toCms( -1*Psum.boostVector() );
102  G4LorentzVector Ptmp = toCms*Pprojectile;
103  if ( Ptmp.pz() <= 0.0 ) return false;
104  // "String" moving backwards in CMS, abort collision !
105  //G4cout << " abort Collision! " << G4endl;
106  toCms.rotateZ( -1*Ptmp.phi() );
107  toCms.rotateY( -1*Ptmp.theta() );
108  G4LorentzRotation toLab( toCms.inverse() );
109  Pprojectile.transform( toCms );
110  Ptarget.transform( toCms );
111 
112  // Putting on mass-on-shell, if needed
113  G4double PZcms2, PZcms;
114  G4double S = Psum.mag2();
115  G4double SqrtS = std::sqrt( S );
116  if ( SqrtS < M0projectile + M0target ) return false;
117 
118  PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( M0target2 )
119  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
120 
121  if ( PZcms2 < 0.0 ) { // It can be in an interaction with off-shell nuclear nucleon
122  if ( M0projectile > projectile->GetDefinition()->GetPDGMass() ) {
123  // An attempt to de-excite the projectile
124  // It is assumed that the target is in the ground state
125  M0projectile = projectile->GetDefinition()->GetPDGMass();
126  M0projectile2 = M0projectile * M0projectile;
127  PZcms2= ( S*S + sqr( M0projectile2 ) + sqr( M0target2 )
128  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
129  if ( PZcms2 < 0.0 ) { return false; } // Nonsuccesful attempt to de-excitate the projectile
130  } else {
131  return false; // The projectile was not excited, but the energy was too low to put
132  // the target nucleon on mass-shell
133  }
134  }
135 
136  PZcms = std::sqrt( PZcms2 );
137 
138  if ( PutOnMassShell ) {
139  if ( Pprojectile.z() > 0.0 ) {
140  Pprojectile.setPz( PZcms );
141  Ptarget.setPz( -PZcms );
142  } else {
143  Pprojectile.setPz( -PZcms );
144  Ptarget.setPz( PZcms );
145  };
146  Pprojectile.setE( std::sqrt( M0projectile2 + Pprojectile.x() * Pprojectile.x() +
147  Pprojectile.y() * Pprojectile.y() + PZcms2 ) );
148  Ptarget.setE( std::sqrt( M0target2 + Ptarget.x() * Ptarget.x() + Ptarget.y() * Ptarget.y() +
149  PZcms2 ) );
150  }
151 
152  G4double maxPtSquare = PZcms2;
153 
154  // Now we can calculate the transferred Pt
155  G4double Pt2;
156  G4double ProjMassT2, ProjMassT;
157  G4double TargMassT2, TargMassT;
158  G4LorentzVector Qmomentum;
159 
160  do {
161  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0.0 );
162  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
163  ProjMassT2 = M0projectile2 + Pt2;
164  ProjMassT = std::sqrt( ProjMassT2 );
165  TargMassT2 = M0target2 + Pt2;
166  TargMassT = std::sqrt( TargMassT2 );
167  } while ( SqrtS < ProjMassT + TargMassT );
168 
169  PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( TargMassT2 )
170  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
171 
172  if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // to avoid the exactness problem
173  PZcms = std::sqrt( PZcms2 );
174  Pprojectile.setPz( PZcms );
175  Ptarget.setPz( -PZcms );
176  Pprojectile += Qmomentum;
177  Ptarget -= Qmomentum;
178 
179  // Transform back and update SplitableHadron Participant.
180  Pprojectile.transform( toLab );
181  Ptarget.transform( toLab );
182 
183  // Calculation of the creation time
184  projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
185  projectile->SetPosition( target->GetPosition() );
186 
187  // Creation time and position of target nucleon were determined at
188  // ReggeonCascade() of G4FTFModel
189 
190  projectile->Set4Momentum( Pprojectile );
191  target->Set4Momentum( Ptarget );
192 
193  //projectile->IncrementCollisionCount( 1 );
194  //target->IncrementCollisionCount( 1 );
195 
196  return true;
197 }
198 
199 
200 //============================================================================
201 
203  G4double maxPtSquare ) const {
204  // @@ this method is used in FTFModel as well. Should go somewhere common!
205  G4double Pt2( 0.0 );
206  if ( AveragePt2 <= 0.0 ) {
207  Pt2 = 0.0;
208  } else {
209  Pt2 = -AveragePt2 * std::log( 1.0 + G4UniformRand() *
210  ( std::exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
211  }
212  G4double Pt = std::sqrt( Pt2 );
213  G4double phi = G4UniformRand() * twopi;
214  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
215 }
216 
217 
218 //============================================================================
219 
221  throw G4HadronicException( __FILE__, __LINE__,
222  "G4ElasticHNScattering copy contructor not meant to be called" );
223 }
224 
225 
226 //============================================================================
227 
229 
230 
231 //============================================================================
232 
234  throw G4HadronicException( __FILE__, __LINE__,
235  "G4ElasticHNScattering = operator not meant to be called" );
236 }
237 
238 
239 //============================================================================
240 
242  throw G4HadronicException( __FILE__, __LINE__,
243  "G4ElasticHNScattering == operator not meant to be called" );
244 }
245 
246 
247 //============================================================================
248 
250  throw G4HadronicException( __FILE__, __LINE__,
251  "G4ElasticHNScattering != operator not meant to be called" );
252 }
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
CLHEP::Hep3Vector G4ThreeVector
int operator==(const G4ElasticHNScattering &right) const
CLHEP::HepLorentzRotation G4LorentzRotation
void SetTimeOfCreation(G4double aTime)
G4double GetAvaragePt2ofElasticScattering()
const G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:93
const G4ElasticHNScattering & operator=(const G4ElasticHNScattering &right)
bool G4bool
Definition: G4Types.hh:79
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
void IncrementCollisionCount(G4int aCount)
const G4LorentzVector & Get4Momentum() const
G4double GetMinimumMass(const G4ParticleDefinition *p) const
G4double GetPDGMass() const
void SetPosition(const G4ThreeVector &aPosition)
const G4ThreeVector & GetPosition() const
int operator!=(const G4ElasticHNScattering &right) const
void Set4Momentum(const G4LorentzVector &a4Momentum)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector