| 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 // $Id: G4ElasticHNScattering.cc,v 1.14 2009-12-16 17:51:13 gunter Exp $
28 <<
29 // ------------------------------------------- 28 // ------------------------------------------------------------
30 // GEANT 4 class implemetation file 29 // GEANT 4 class implemetation file
31 // 30 //
32 // ---------------- G4ElasticHNScattering 31 // ---------------- G4ElasticHNScattering --------------
33 // by V. Uzhinsky, March 200 32 // by V. Uzhinsky, March 2008.
34 // elastic scattering used by Frit 33 // elastic scattering used by Fritiof model
35 // Take a projectile and a tar << 34 // Take a projectile and a target
36 // scatter the projectile and << 35 // scatter the projectile and target
37 // ------------------------------------------- 36 // ---------------------------------------------------------------------
38 37
>> 38
39 #include "globals.hh" 39 #include "globals.hh"
40 #include "Randomize.hh" 40 #include "Randomize.hh"
41 #include "G4PhysicalConstants.hh" <<
42 #include "G4SystemOfUnits.hh" <<
43 41
44 #include "G4ElasticHNScattering.hh" 42 #include "G4ElasticHNScattering.hh"
45 #include "G4LorentzRotation.hh" 43 #include "G4LorentzRotation.hh"
46 #include "G4ThreeVector.hh" 44 #include "G4ThreeVector.hh"
47 #include "G4ParticleDefinition.hh" 45 #include "G4ParticleDefinition.hh"
48 #include "G4VSplitableHadron.hh" 46 #include "G4VSplitableHadron.hh"
49 #include "G4ExcitedString.hh" 47 #include "G4ExcitedString.hh"
50 #include "G4FTFParameters.hh" 48 #include "G4FTFParameters.hh"
>> 49 //#include "G4ios.hh"
51 50
52 #include "G4SampleResonance.hh" << 51 G4ElasticHNScattering::G4ElasticHNScattering()
53 << 52 {
54 #include "G4Exp.hh" <<
55 #include "G4Log.hh" <<
56 <<
57 //============================================ <<
58 <<
59 G4ElasticHNScattering::G4ElasticHNScattering() <<
60 <<
61 <<
62 //============================================ <<
63 <<
64 G4bool G4ElasticHNScattering::ElasticScatterin <<
65 <<
66 <<
67 projectile->IncrementCollisionCount( 1 ); <<
68 target->IncrementCollisionCount( 1 ); <<
69 <<
70 if ( projectile->Get4Momentum().z() < 0.0 ) <<
71 <<
72 // Projectile parameters <<
73 G4LorentzVector Pprojectile = projectile->Ge <<
74 G4double M0projectile = Pprojectile.mag(); <<
75 G4double M0projectile2 = M0projectile * M0pr <<
76 <<
77 // Target parameters <<
78 G4LorentzVector Ptarget = target->Get4Moment <<
79 G4double M0target = Ptarget.mag(); <<
80 G4double M0target2 = M0target * M0target; <<
81 <<
82 G4double AveragePt2 = theParameters->GetAvar <<
83 <<
84 // Transform momenta to cms and then rotate <<
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 c <<
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 ) retur <<
102 <<
103 PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( <<
104 - 2*S*M0projectile2 - 2*S*M0targe <<
105 <<
106 PZcms = ( PZcms2 > 0.0 ? std::sqrt( PZcms2 ) <<
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( A <<
120 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma <<
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 ) <<
127 if ( loopCounter >= maxNumberOfLoops ) { <<
128 return false; <<
129 } <<
130 <<
131 PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( Ta <<
132 - 2.0*S*ProjMassT2 - 2.0*S*TargMa <<
133 <<
134 if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // t <<
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 <<
142 Pprojectile.transform( toLab ); <<
143 Ptarget.transform( toLab ); <<
144 <<
145 // Calculation of the creation time <<
146 projectile->SetTimeOfCreation( target->GetTi <<
147 projectile->SetPosition( target->GetPosition <<
148 <<
149 // Creation time and position of target nucl <<
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 } 53 }
160 54
161 << 55 G4bool G4ElasticHNScattering::
162 //============================================ << 56 ElasticScattering (G4VSplitableHadron *projectile,
163 << 57 G4VSplitableHadron *target,
164 G4ThreeVector G4ElasticHNScattering::GaussianP << 58 G4FTFParameters *theParameters) const
165 << 59 {
166 // @@ this method is used in FTFModel as wel << 60 // -------------------- Projectile parameters -----------------------------------
167 G4double Pt2( 0.0 ); << 61 G4LorentzVector Pprojectile=projectile->Get4Momentum();
168 if ( AveragePt2 <= 0.0 ) { << 62
169 Pt2 = 0.0; << 63 if(Pprojectile.z() < 0.)
170 } else { << 64 {
171 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Uniform << 65 target->SetStatus(2);
172 } << 66 return false;
173 G4double Pt = ( Pt2 > 0.0 ? std::sqrt( Pt2 ) << 67 }
174 G4double phi = G4UniformRand() * twopi; << 68
175 return G4ThreeVector( Pt * std::cos( phi ), << 69 G4bool PutOnMassShell(false);
>> 70
>> 71 G4double M0projectile = Pprojectile.mag();
>> 72 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
>> 73 {
>> 74 PutOnMassShell=true;
>> 75 M0projectile=projectile->GetDefinition()->GetPDGMass();
>> 76 }
>> 77
>> 78 G4double Mprojectile2 = M0projectile * M0projectile;
>> 79
>> 80 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
>> 81
>> 82 // -------------------- Target parameters ----------------------------------------------
>> 83
>> 84 G4LorentzVector Ptarget=target->Get4Momentum();
>> 85
>> 86 G4double M0target = Ptarget.mag();
>> 87
>> 88 if(M0target < target->GetDefinition()->GetPDGMass())
>> 89 {
>> 90 PutOnMassShell=true;
>> 91 M0target=target->GetDefinition()->GetPDGMass();
>> 92 }
>> 93
>> 94 G4double Mtarget2 = M0target * M0target;
>> 95
>> 96 // Transform momenta to cms and then rotate parallel to z axis;
>> 97
>> 98 G4LorentzVector Psum;
>> 99 Psum=Pprojectile+Ptarget;
>> 100
>> 101 G4LorentzRotation toCms(-1*Psum.boostVector());
>> 102
>> 103 G4LorentzVector Ptmp=toCms*Pprojectile;
>> 104
>> 105 if ( Ptmp.pz() <= 0. )
>> 106 {
>> 107 // "String" moving backwards in CMS, abort collision !!
>> 108 //G4cout << " abort Collision!! " << G4endl;
>> 109 target->SetStatus(2);
>> 110 return false;
>> 111 }
>> 112
>> 113 toCms.rotateZ(-1*Ptmp.phi());
>> 114 toCms.rotateY(-1*Ptmp.theta());
>> 115
>> 116 G4LorentzRotation toLab(toCms.inverse());
>> 117
>> 118 Pprojectile.transform(toCms);
>> 119 Ptarget.transform(toCms);
>> 120
>> 121 // ---------------------- Putting on mass-on-shell, if needed ------------------------
>> 122 G4double PZcms2, PZcms;
>> 123
>> 124 G4double S=Psum.mag2();
>> 125 // G4double SqrtS=std::sqrt(S);
>> 126
>> 127 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
>> 128 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
>> 129
>> 130 if(PZcms2 < 0.)
>> 131 { // It can be in an interaction with off-shell nuclear nucleon
>> 132 if(M0projectile > projectile->GetDefinition()->GetPDGMass())
>> 133 { // An attempt to de-excite the projectile
>> 134 // It is assumed that the target is in the ground state
>> 135 M0projectile = projectile->GetDefinition()->GetPDGMass();
>> 136 Mprojectile2=M0projectile*M0projectile;
>> 137 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
>> 138 2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2)
>> 139 /4./S;
>> 140 if(PZcms2 < 0.){ return false;} // Non succesful attempt after the de-excitation
>> 141 }
>> 142 else // if(M0projectile > projectile->GetDefinition()->GetPDGMass())
>> 143 {
>> 144 target->SetStatus(2);
>> 145 return false; // The projectile was not excited,
>> 146 // but the energy was too low to put
>> 147 // the target nucleon on mass-shell
>> 148 } // end of if(M0projectile > projectile->GetDefinition()->GetPDGMass())
>> 149 } // end of if(PZcms2 < 0.)
>> 150
>> 151 PZcms = std::sqrt(PZcms2);
>> 152
>> 153 if(PutOnMassShell)
>> 154 {
>> 155 if(Pprojectile.z() > 0.)
>> 156 {
>> 157 Pprojectile.setPz( PZcms);
>> 158 Ptarget.setPz( -PZcms);
>> 159 }
>> 160 else // if(Pprojectile.z() > 0.)
>> 161 {
>> 162 Pprojectile.setPz(-PZcms);
>> 163 Ptarget.setPz( PZcms);
>> 164 };
>> 165
>> 166 Pprojectile.setE(std::sqrt(Mprojectile2+
>> 167 Pprojectile.x()*Pprojectile.x()+
>> 168 Pprojectile.y()*Pprojectile.y()+
>> 169 PZcms2));
>> 170 Ptarget.setE(std::sqrt( Mtarget2 +
>> 171 Ptarget.x()*Ptarget.x()+
>> 172 Ptarget.y()*Ptarget.y()+
>> 173 PZcms2));
>> 174 } // end of if(PutOnMassShell)
>> 175
>> 176 G4double maxPtSquare = PZcms2;
>> 177
>> 178 // ------ Now we can calculate the transfered Pt --------------------------
>> 179 G4double Pt2;
>> 180 G4double ProjMassT2, ProjMassT;
>> 181 G4double TargMassT2, TargMassT;
>> 182
>> 183 G4LorentzVector Qmomentum;
>> 184
>> 185 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
>> 186
>> 187 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
>> 188
>> 189 ProjMassT2=Mprojectile2+Pt2;
>> 190 ProjMassT =std::sqrt(ProjMassT2);
>> 191
>> 192 TargMassT2=Mtarget2+Pt2;
>> 193 TargMassT =std::sqrt(TargMassT2);
>> 194
>> 195 PZcms2=(S*S+ProjMassT2*ProjMassT2+
>> 196 TargMassT2*TargMassT2-
>> 197 2.*S*ProjMassT2-2.*S*TargMassT2-
>> 198 2.*ProjMassT2*TargMassT2)/4./S;
>> 199 if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem
>> 200 PZcms =std::sqrt(PZcms2);
>> 201
>> 202 Pprojectile.setPz( PZcms);
>> 203 Ptarget.setPz( -PZcms);
>> 204
>> 205 Pprojectile += Qmomentum;
>> 206 Ptarget -= Qmomentum;
>> 207
>> 208 // Transform back and update SplitableHadron Participant.
>> 209 Pprojectile.transform(toLab);
>> 210 Ptarget.transform(toLab);
>> 211 /* // Maybe it will be needed for an exact calculations--------------------
>> 212 G4double TargetMomentum=std::sqrt(Ptarget.x()*Ptarget.x()+
>> 213 Ptarget.y()*Ptarget.y()+
>> 214 Ptarget.z()*Ptarget.z());
>> 215 */
>> 216
>> 217 // Calculation of the creation time ---------------------
>> 218 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
>> 219 projectile->SetPosition(target->GetPosition());
>> 220 // Creation time and position of target nucleon were determined at
>> 221 // ReggeonCascade() of G4FTFModel
>> 222 // ------------------------------------------------------
>> 223
>> 224 projectile->Set4Momentum(Pprojectile);
>> 225 target->Set4Momentum(Ptarget);
>> 226
>> 227 projectile->IncrementCollisionCount(1);
>> 228 target->IncrementCollisionCount(1);
>> 229
>> 230 return true;
>> 231 }
>> 232
>> 233
>> 234 // --------- private methods ----------------------
>> 235
>> 236 G4ThreeVector G4ElasticHNScattering::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
>> 237 { // @@ this method is used in FTFModel as well. Should go somewhere common!
>> 238
>> 239 G4double Pt2(0.);
>> 240 if(AveragePt2 <= 0.) {Pt2=0.;}
>> 241 else
>> 242 {
>> 243 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
>> 244 (std::exp(-maxPtSquare/AveragePt2)-1.));
>> 245 }
>> 246 G4double Pt=std::sqrt(Pt2);
>> 247 G4double phi=G4UniformRand() * twopi;
>> 248
>> 249 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
>> 250 }
>> 251
>> 252 G4ElasticHNScattering::G4ElasticHNScattering(const G4ElasticHNScattering &)
>> 253 {
>> 254 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering copy contructor not meant to be called");
>> 255 }
>> 256
>> 257
>> 258 G4ElasticHNScattering::~G4ElasticHNScattering()
>> 259 {
>> 260 }
>> 261
>> 262
>> 263 const G4ElasticHNScattering & G4ElasticHNScattering::operator=(const G4ElasticHNScattering &)
>> 264 {
>> 265 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering = operator meant to be called");
>> 266 return *this;
>> 267 }
>> 268
>> 269
>> 270 int G4ElasticHNScattering::operator==(const G4ElasticHNScattering &) const
>> 271 {
>> 272 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering == operator meant to be called");
>> 273 return false;
>> 274 }
>> 275
>> 276 int G4ElasticHNScattering::operator!=(const G4ElasticHNScattering &) const
>> 277 {
>> 278 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering != operator meant to be called");
>> 279 return true;
176 } 280 }
177 <<
178 <<
179 //============================================ <<
180 <<
181 G4ElasticHNScattering::G4ElasticHNScattering( <<
182 throw G4HadronicException( __FILE__, __LINE_ <<
183 "G4ElasticHNScatt <<
184 } <<
185 <<
186 <<
187 //============================================ <<
188 <<
189 G4ElasticHNScattering::~G4ElasticHNScattering( <<
190 <<
191 <<
192 //============================================ <<
193 <<
194 const G4ElasticHNScattering & G4ElasticHNScatt <<
195 throw G4HadronicException( __FILE__, __LINE_ <<
196 "G4ElasticHNScatt <<
197 } <<
198 <<
199 <<
200 //============================================ <<
201 <<
202 G4bool G4ElasticHNScattering::operator==( cons <<
203 throw G4HadronicException( __FILE__, __LINE__ <<
204 "G4ElasticHNScatte <<
205 } <<
206 <<
207 <<
208 //============================================ <<
209 <<
210 G4bool G4ElasticHNScattering::operator!=( cons <<
211 throw G4HadronicException( __FILE__, __LINE_ <<
212 "G4ElasticHNScatte <<
213 } <<
214 <<
215 281