| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 28
29 // ------------------------------------------- 29 // ------------------------------------------------------------
30 // GEANT 4 class implemetation file 30 // GEANT 4 class implemetation file
31 // 31 //
32 // ---------------- G4ElasticHNScattering 32 // ---------------- G4ElasticHNScattering --------------
33 // by V. Uzhinsky, March 200 33 // by V. Uzhinsky, March 2008.
34 // elastic scattering used by Frit 34 // elastic scattering used by Fritiof model
35 // Take a projectile and a tar 35 // Take a projectile and a target
36 // scatter the projectile and 36 // scatter the projectile and target
37 // ------------------------------------------- 37 // ---------------------------------------------------------------------
38 38
39 #include "globals.hh" 39 #include "globals.hh"
40 #include "Randomize.hh" 40 #include "Randomize.hh"
41 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 43
44 #include "G4ElasticHNScattering.hh" 44 #include "G4ElasticHNScattering.hh"
45 #include "G4LorentzRotation.hh" 45 #include "G4LorentzRotation.hh"
46 #include "G4ThreeVector.hh" 46 #include "G4ThreeVector.hh"
47 #include "G4ParticleDefinition.hh" 47 #include "G4ParticleDefinition.hh"
48 #include "G4VSplitableHadron.hh" 48 #include "G4VSplitableHadron.hh"
49 #include "G4ExcitedString.hh" 49 #include "G4ExcitedString.hh"
50 #include "G4FTFParameters.hh" 50 #include "G4FTFParameters.hh"
51 51
52 #include "G4SampleResonance.hh" 52 #include "G4SampleResonance.hh"
53 53
54 #include "G4Exp.hh" 54 #include "G4Exp.hh"
55 #include "G4Log.hh" 55 #include "G4Log.hh"
56 56
57 //============================================ 57 //============================================================================
58 58
59 G4ElasticHNScattering::G4ElasticHNScattering() 59 G4ElasticHNScattering::G4ElasticHNScattering() {}
60 60
61 61
62 //============================================ 62 //============================================================================
63 63
64 G4bool G4ElasticHNScattering::ElasticScatterin 64 G4bool G4ElasticHNScattering::ElasticScattering( G4VSplitableHadron* projectile,
65 65 G4VSplitableHadron* target,
66 66 G4FTFParameters* theParameters ) const {
67 projectile->IncrementCollisionCount( 1 ); 67 projectile->IncrementCollisionCount( 1 );
68 target->IncrementCollisionCount( 1 ); 68 target->IncrementCollisionCount( 1 );
69 69
70 if ( projectile->Get4Momentum().z() < 0.0 ) 70 if ( projectile->Get4Momentum().z() < 0.0 ) return false; //Uzhi Aug.2019
71 71
72 // Projectile parameters 72 // Projectile parameters
73 G4LorentzVector Pprojectile = projectile->Ge 73 G4LorentzVector Pprojectile = projectile->Get4Momentum();
74 G4double M0projectile = Pprojectile.mag(); 74 G4double M0projectile = Pprojectile.mag();
75 G4double M0projectile2 = M0projectile * M0pr 75 G4double M0projectile2 = M0projectile * M0projectile;
76 76
77 // Target parameters 77 // Target parameters
78 G4LorentzVector Ptarget = target->Get4Moment 78 G4LorentzVector Ptarget = target->Get4Momentum();
79 G4double M0target = Ptarget.mag(); 79 G4double M0target = Ptarget.mag();
80 G4double M0target2 = M0target * M0target; 80 G4double M0target2 = M0target * M0target;
81 81
82 G4double AveragePt2 = theParameters->GetAvar 82 G4double AveragePt2 = theParameters->GetAvaragePt2ofElasticScattering();
83 83
84 // Transform momenta to cms and then rotate 84 // Transform momenta to cms and then rotate parallel to z axis;
85 G4LorentzVector Psum; 85 G4LorentzVector Psum;
86 Psum = Pprojectile + Ptarget; 86 Psum = Pprojectile + Ptarget;
87 G4LorentzRotation toCms( -1*Psum.boostVector 87 G4LorentzRotation toCms( -1*Psum.boostVector() );
88 G4LorentzVector Ptmp = toCms*Pprojectile; 88 G4LorentzVector Ptmp = toCms*Pprojectile;
89 if ( Ptmp.pz() <= 0.0 ) return false; 89 if ( Ptmp.pz() <= 0.0 ) return false;
90 //"String" moving backwards in CMS, abort c 90 //"String" moving backwards in CMS, abort collision !
91 //G4cout << " abort Collision! " << G4endl; 91 //G4cout << " abort Collision! " << G4endl;
92 toCms.rotateZ( -1*Ptmp.phi() ); 92 toCms.rotateZ( -1*Ptmp.phi() );
93 toCms.rotateY( -1*Ptmp.theta() ); 93 toCms.rotateY( -1*Ptmp.theta() );
94 G4LorentzRotation toLab( toCms.inverse() ); 94 G4LorentzRotation toLab( toCms.inverse() );
95 Pprojectile.transform( toCms ); 95 Pprojectile.transform( toCms );
96 Ptarget.transform( toCms ); 96 Ptarget.transform( toCms );
97 97
98 G4double PZcms2, PZcms; 98 G4double PZcms2, PZcms;
99 G4double S = Psum.mag2(); 99 G4double S = Psum.mag2();
100 G4double SqrtS = std::sqrt( S ); 100 G4double SqrtS = std::sqrt( S );
101 if ( SqrtS < M0projectile + M0target ) retur 101 if ( SqrtS < M0projectile + M0target ) return false;
102 102
103 PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( 103 PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( M0target2 )
104 - 2*S*M0projectile2 - 2*S*M0targe 104 - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
105 105
106 PZcms = ( PZcms2 > 0.0 ? std::sqrt( PZcms2 ) 106 PZcms = ( PZcms2 > 0.0 ? std::sqrt( PZcms2 ) : 0.0 );
107 107
108 G4double maxPtSquare = PZcms2; 108 G4double maxPtSquare = PZcms2;
109 109
110 // Now we can calculate the transferred Pt 110 // Now we can calculate the transferred Pt
111 G4double Pt2; 111 G4double Pt2;
112 G4double ProjMassT2, ProjMassT; 112 G4double ProjMassT2, ProjMassT;
113 G4double TargMassT2, TargMassT; 113 G4double TargMassT2, TargMassT;
114 G4LorentzVector Qmomentum; 114 G4LorentzVector Qmomentum;
115 115
116 const G4int maxNumberOfLoops = 1000; 116 const G4int maxNumberOfLoops = 1000;
117 G4int loopCounter = 0; 117 G4int loopCounter = 0;
118 do { 118 do {
119 Qmomentum = G4LorentzVector( GaussianPt( A 119 Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0.0 );
120 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma 120 Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
121 ProjMassT2 = M0projectile2 + Pt2; 121 ProjMassT2 = M0projectile2 + Pt2;
122 ProjMassT = std::sqrt( ProjMassT2 ); 122 ProjMassT = std::sqrt( ProjMassT2 );
123 TargMassT2 = M0target2 + Pt2; 123 TargMassT2 = M0target2 + Pt2;
124 TargMassT = std::sqrt( TargMassT2 ); 124 TargMassT = std::sqrt( TargMassT2 );
125 } while ( ( SqrtS < ProjMassT + TargMassT ) 125 } while ( ( SqrtS < ProjMassT + TargMassT ) &&
126 ++loopCounter < maxNumberOfLoops ) 126 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
127 if ( loopCounter >= maxNumberOfLoops ) { 127 if ( loopCounter >= maxNumberOfLoops ) {
128 return false; 128 return false;
129 } 129 }
130 130
131 PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( Ta 131 PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( TargMassT2 )
132 - 2.0*S*ProjMassT2 - 2.0*S*TargMa 132 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
133 133
134 if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // t 134 if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // to avoid the exactness problem
135 PZcms = std::sqrt( PZcms2 ); 135 PZcms = std::sqrt( PZcms2 );
136 Pprojectile.setPz( PZcms ); 136 Pprojectile.setPz( PZcms );
137 Ptarget.setPz( -PZcms ); 137 Ptarget.setPz( -PZcms );
138 Pprojectile += Qmomentum; 138 Pprojectile += Qmomentum;
139 Ptarget -= Qmomentum; 139 Ptarget -= Qmomentum;
140 140
141 // Transform back and update SplitableHadron 141 // Transform back and update SplitableHadron Participant.
142 Pprojectile.transform( toLab ); 142 Pprojectile.transform( toLab );
143 Ptarget.transform( toLab ); 143 Ptarget.transform( toLab );
144 144
145 // Calculation of the creation time 145 // Calculation of the creation time
146 projectile->SetTimeOfCreation( target->GetTi 146 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
147 projectile->SetPosition( target->GetPosition 147 projectile->SetPosition( target->GetPosition() );
148 148
149 // Creation time and position of target nucl 149 // Creation time and position of target nucleon were determined at
150 // ReggeonCascade() of G4FTFModel 150 // ReggeonCascade() of G4FTFModel
151 151
152 projectile->Set4Momentum( Pprojectile ); 152 projectile->Set4Momentum( Pprojectile );
153 target->Set4Momentum( Ptarget ); 153 target->Set4Momentum( Ptarget );
154 154
155 //projectile->IncrementCollisionCount( 1 ); 155 //projectile->IncrementCollisionCount( 1 );
156 //target->IncrementCollisionCount( 1 ); 156 //target->IncrementCollisionCount( 1 );
157 157
158 return true; 158 return true;
159 } 159 }
160 160
161 161
162 //============================================ 162 //============================================================================
163 163
164 G4ThreeVector G4ElasticHNScattering::GaussianP 164 G4ThreeVector G4ElasticHNScattering::GaussianPt( G4double AveragePt2,
165 165 G4double maxPtSquare ) const {
166 // @@ this method is used in FTFModel as wel 166 // @@ this method is used in FTFModel as well. Should go somewhere common!
167 G4double Pt2( 0.0 ); 167 G4double Pt2( 0.0 );
168 if ( AveragePt2 <= 0.0 ) { 168 if ( AveragePt2 <= 0.0 ) {
169 Pt2 = 0.0; 169 Pt2 = 0.0;
170 } else { 170 } else {
171 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Uniform 171 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() * ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
172 } 172 }
173 G4double Pt = ( Pt2 > 0.0 ? std::sqrt( Pt2 ) 173 G4double Pt = ( Pt2 > 0.0 ? std::sqrt( Pt2 ) : 0.0 );
174 G4double phi = G4UniformRand() * twopi; 174 G4double phi = G4UniformRand() * twopi;
175 return G4ThreeVector( Pt * std::cos( phi ), 175 return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
176 } 176 }
177 177
178 178
179 //============================================ 179 //============================================================================
180 180
181 G4ElasticHNScattering::G4ElasticHNScattering( 181 G4ElasticHNScattering::G4ElasticHNScattering( const G4ElasticHNScattering& ) {
182 throw G4HadronicException( __FILE__, __LINE_ 182 throw G4HadronicException( __FILE__, __LINE__,
183 "G4ElasticHNScatt 183 "G4ElasticHNScattering copy constructor not meant to be called" );
184 } 184 }
185 185
186 186
187 //============================================ 187 //============================================================================
188 188
189 G4ElasticHNScattering::~G4ElasticHNScattering( 189 G4ElasticHNScattering::~G4ElasticHNScattering() {}
190 190
191 191
192 //============================================ 192 //============================================================================
193 193
194 const G4ElasticHNScattering & G4ElasticHNScatt 194 const G4ElasticHNScattering & G4ElasticHNScattering::operator=( const G4ElasticHNScattering& ) {
195 throw G4HadronicException( __FILE__, __LINE_ 195 throw G4HadronicException( __FILE__, __LINE__,
196 "G4ElasticHNScatt 196 "G4ElasticHNScattering = operator not meant to be called" );
197 } 197 }
198 198
199 199
200 //============================================ 200 //============================================================================
201 201
202 G4bool G4ElasticHNScattering::operator==( cons 202 G4bool G4ElasticHNScattering::operator==( const G4ElasticHNScattering& ) const {
203 throw G4HadronicException( __FILE__, __LINE__ 203 throw G4HadronicException( __FILE__, __LINE__,
204 "G4ElasticHNScatte 204 "G4ElasticHNScattering == operator not meant to be called" );
205 } 205 }
206 206
207 207
208 //============================================ 208 //============================================================================
209 209
210 G4bool G4ElasticHNScattering::operator!=( cons 210 G4bool G4ElasticHNScattering::operator!=( const G4ElasticHNScattering& ) const {
211 throw G4HadronicException( __FILE__, __LINE_ 211 throw G4HadronicException( __FILE__, __LINE__,
212 "G4ElasticHNScatte 212 "G4ElasticHNScattering != operator not meant to be called" );
213 } 213 }
214 214
215 215