| 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 // GEANT 4 class implemetation file 28 // GEANT 4 class implemetation file
29 // 29 //
30 // ---------------- G4SingleDiffractiveEx 30 // ---------------- G4SingleDiffractiveExcitation --------------
31 // by Gunter Folger, October 1998. 31 // by Gunter Folger, October 1998.
32 // diffractive Excitation used by strings 32 // diffractive Excitation used by strings models
33 // Take a projectile and a target 33 // Take a projectile and a target
34 // excite the projectile and target 34 // excite the projectile and target
35 // ------------------------------------------- 35 // ------------------------------------------------------------
36 36
37 #include "G4SingleDiffractiveExcitation.hh" 37 #include "G4SingleDiffractiveExcitation.hh"
38 #include "globals.hh" 38 #include "globals.hh"
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 #include "Randomize.hh" 41 #include "Randomize.hh"
42 #include "G4LorentzRotation.hh" 42 #include "G4LorentzRotation.hh"
43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh"
44 #include "G4ParticleDefinition.hh" 44 #include "G4ParticleDefinition.hh"
45 #include "G4VSplitableHadron.hh" 45 #include "G4VSplitableHadron.hh"
46 #include "G4ExcitedString.hh" 46 #include "G4ExcitedString.hh"
47 47
48 #include "G4Log.hh" 48 #include "G4Log.hh"
49 #include "G4Pow.hh" 49 #include "G4Pow.hh"
50 50
51 //#define debugSingleDiffraction 51 //#define debugSingleDiffraction
52 52
53 G4SingleDiffractiveExcitation::G4SingleDiffrac 53 G4SingleDiffractiveExcitation::G4SingleDiffractiveExcitation(){}
54 54
55 G4SingleDiffractiveExcitation::~G4SingleDiffra 55 G4SingleDiffractiveExcitation::~G4SingleDiffractiveExcitation(){}
56 56
57 G4bool G4SingleDiffractiveExcitation:: 57 G4bool G4SingleDiffractiveExcitation::
58 ExciteParticipants( G4VSplitableHadron *projec 58 ExciteParticipants( G4VSplitableHadron *projectile, G4VSplitableHadron *target,
59 G4bool ProjectileDiffracti 59 G4bool ProjectileDiffraction ) const
60 { 60 {
61 #ifdef debugSingleDiffraction 61 #ifdef debugSingleDiffraction
62 G4cout<<G4endl<<"G4SingleDiffractiveExcita 62 G4cout<<G4endl<<"G4SingleDiffractiveExcitation::ExciteParticipants"<<G4endl;
63 #endif 63 #endif
64 64
65 G4LorentzVector Pprojectile=projectile->Get4 65 G4LorentzVector Pprojectile=projectile->Get4Momentum();
66 G4double Mprojectile = projectile->GetDef 66 G4double Mprojectile = projectile->GetDefinition()->GetPDGMass();
67 G4double Mprojectile2=sqr(projectile->GetDef 67 G4double Mprojectile2=sqr(projectile->GetDefinition()->GetPDGMass());
68 68
69 G4LorentzVector Ptarget=target->Get4Momentum 69 G4LorentzVector Ptarget=target->Get4Momentum();
70 G4double Mtarget = target->GetDefinition( 70 G4double Mtarget = target->GetDefinition()->GetPDGMass();
71 G4double Mtarget2=sqr(target->GetDefinition( 71 G4double Mtarget2=sqr(target->GetDefinition()->GetPDGMass());
72 72
73 #ifdef debugSingleDiffraction 73 #ifdef debugSingleDiffraction
74 G4cout<<"Proj Targ "<<projectile->GetDefin 74 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl;
75 G4cout<<"Pr Tr 4-Mom "<<Pprojectile<<" "<< 75 G4cout<<"Pr Tr 4-Mom "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl
76 <<" "<<Ptarget <<" "<< 76 <<" "<<Ptarget <<" "<<Ptarget.mag() <<G4endl;
77 #endif 77 #endif
78 78
79 G4LorentzVector Psum=Pprojectile+Ptarget; 79 G4LorentzVector Psum=Pprojectile+Ptarget;
80 G4double SqrtS=Psum.mag(); 80 G4double SqrtS=Psum.mag();
81 G4double S =Psum.mag2(); 81 G4double S =Psum.mag2();
82 82
83 #ifdef debugSingleDiffraction 83 #ifdef debugSingleDiffraction
84 G4cout<<"SqrtS-Mprojectile-Mtarget "<<Sqrt 84 G4cout<<"SqrtS-Mprojectile-Mtarget "<<SqrtS<<" "<<Mprojectile<<" "<<Mtarget
85 <<" "<<SqrtS-Mprojectile-Mtarget<<G4 85 <<" "<<SqrtS-Mprojectile-Mtarget<<G4endl;
86 #endif 86 #endif
87 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) 87 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) {
88 #ifdef debugSingleDiffraction 88 #ifdef debugSingleDiffraction
89 G4cerr<<"Projectile: "<<projectile->GetDef 89 G4cerr<<"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<" "
90 <<Pprojectile<<" "<<Pprojectile.mag( 90 <<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
91 G4cerr<<"Target: "<<target->GetDefinit 91 G4cerr<<"Target: "<<target->GetDefinition()->GetPDGEncoding()<<" "
92 <<Ptarget<<" "<<Ptarget.mag()<<G4end 92 <<Ptarget<<" "<<Ptarget.mag()<<G4endl;
93 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = " 93 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "<<Pprojectile.mag()+Ptarget.mag()<<G4endl;
94 #endif 94 #endif
95 return true; 95 return true;
96 } 96 }
97 97
98 G4LorentzRotation toCms(-1*Psum.boostVector( 98 G4LorentzRotation toCms(-1*Psum.boostVector());
99 99
100 G4LorentzVector Ptmp=toCms*Pprojectile; 100 G4LorentzVector Ptmp=toCms*Pprojectile;
101 101
102 if ( Ptmp.pz() <= 0. ) 102 if ( Ptmp.pz() <= 0. )
103 { 103 {
104 // "String" moving backwards in CMS, abor 104 // "String" moving backwards in CMS, abort collision !!
105 // G4cout << " abort Collision!! " 105 // G4cout << " abort Collision!! " << G4endl;
106 return false; 106 return false;
107 } 107 }
108 108
109 toCms.rotateZ(-1*Ptmp.phi()); 109 toCms.rotateZ(-1*Ptmp.phi());
110 toCms.rotateY(-1*Ptmp.theta()); 110 toCms.rotateY(-1*Ptmp.theta());
111 111
112 G4LorentzRotation toLab(toCms.inverse()); 112 G4LorentzRotation toLab(toCms.inverse());
113 113
114 Pprojectile.transform(toCms); 114 Pprojectile.transform(toCms);
115 Ptarget.transform(toCms); 115 Ptarget.transform(toCms);
116 #ifdef debugSingleDiffraction 116 #ifdef debugSingleDiffraction
117 G4cout << "Pprojectile in CMS " << Pproje 117 G4cout << "Pprojectile in CMS " << Pprojectile << G4endl;
118 G4cout << "Ptarget in CMS " << Ptarge 118 G4cout << "Ptarget in CMS " << Ptarget << G4endl;
119 #endif 119 #endif
120 G4double maxPtSquare=sqr(Ptarget.pz()); 120 G4double maxPtSquare=sqr(Ptarget.pz());
121 121
122 G4double ProjectileMinDiffrMass(0.), TargetM 122 G4double ProjectileMinDiffrMass(0.), TargetMinDiffrMass(0.);
123 G4double AveragePt2(0.); 123 G4double AveragePt2(0.);
124 G4int absPDGcode=std::abs(projectile->GetDef 124 G4int absPDGcode=std::abs(projectile->GetDefinition()->GetPDGEncoding());
125 125
126 if ( ProjectileDiffraction ) { 126 if ( ProjectileDiffraction ) {
127 if ( absPDGcode > 1000 ) 127 if ( absPDGcode > 1000 ) //------Projectile is baryon --------
128 { 128 {
129 if ( absPDGcode > 4000 && absPDGcode < 6 129 if ( absPDGcode > 4000 && absPDGcode < 6000 ) // Projectile is a charm or bottom baryon
130 { 130 {
131 ProjectileMinDiffrMass = projectile->G 131 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV
132 AveragePt2 = 0.3; 132 AveragePt2 = 0.3; // GeV^2
133 } 133 }
134 else 134 else
135 { 135 {
136 ProjectileMinDiffrMass = 1.16; 136 ProjectileMinDiffrMass = 1.16; // GeV
137 AveragePt2 = 0.3; 137 AveragePt2 = 0.3; // GeV^2
138 } 138 }
139 } 139 }
140 else if( absPDGcode == 211 || absPDGcode = 140 else if( absPDGcode == 211 || absPDGcode == 111) //------Projectile is Pion -----------
141 { 141 {
142 ProjectileMinDiffrMass = 1.0; 142 ProjectileMinDiffrMass = 1.0; // GeV
143 AveragePt2 = 0.3; 143 AveragePt2 = 0.3; // GeV^2
144 } 144 }
145 else if( absPDGcode == 321 || absPDGcode = 145 else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon
146 { 146 {
147 ProjectileMinDiffrMass = 1.1; 147 ProjectileMinDiffrMass = 1.1; // GeV
148 AveragePt2 = 0.3; 148 AveragePt2 = 0.3; // GeV^2
149 } 149 }
150 else if( absPDGcode == 22) 150 else if( absPDGcode == 22) //------Projectile is Gamma -----------
151 { 151 {
152 ProjectileMinDiffrMass = 0.25; 152 ProjectileMinDiffrMass = 0.25; // GeV
153 AveragePt2 = 0.36; 153 AveragePt2 = 0.36; // GeV^2
154 } 154 }
155 else if( absPDGcode > 400 && absPDGcode < 155 else if( absPDGcode > 400 && absPDGcode < 600) // Projectile is a charm or bottom meson
156 { 156 {
157 ProjectileMinDiffrMass = projectile->Get 157 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV
158 AveragePt2 = 0.3; 158 AveragePt2 = 0.3; // GeV^2
159 } 159 }
160 else 160 else //------Projectile is undefined, Nucleon assumed
161 { 161 {
162 ProjectileMinDiffrMass = 1.1; 162 ProjectileMinDiffrMass = 1.1; // GeV
163 AveragePt2 = 0.3; 163 AveragePt2 = 0.3; // GeV^2
164 }; 164 };
165 165
166 ProjectileMinDiffrMass = ProjectileMinDiff 166 ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV;
167 Mprojectile2=sqr(ProjectileMinDiffrMass); 167 Mprojectile2=sqr(ProjectileMinDiffrMass);
168 } 168 }
169 else 169 else
170 { 170 {
171 TargetMinDiffrMass = 1.16*GeV; 171 TargetMinDiffrMass = 1.16*GeV; // For target nucleon
172 Mtarget2 = sqr( TargetMinDiffrMass) ; 172 Mtarget2 = sqr( TargetMinDiffrMass) ;
173 AveragePt2 = 0.3; 173 AveragePt2 = 0.3; // GeV^2
174 } // end of if ( ProjectileDiffraction ) 174 } // end of if ( ProjectileDiffraction )
175 175
176 AveragePt2 = AveragePt2 * GeV*GeV; 176 AveragePt2 = AveragePt2 * GeV*GeV;
177 177
178 G4double Pt2, PZcms, PZcms2; 178 G4double Pt2, PZcms, PZcms2;
179 G4double ProjMassT2, ProjMassT; 179 G4double ProjMassT2, ProjMassT;
180 G4double TargMassT2, TargMassT; 180 G4double TargMassT2, TargMassT;
181 G4double PMinusMin, PMinusMax; 181 G4double PMinusMin, PMinusMax;
182 G4double TPlusMin, TPlusMax; 182 G4double TPlusMin, TPlusMax;
183 G4double PMinusNew, PPlusNew, TPlusNew, TMin 183 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
184 184
185 G4LorentzVector Qmomentum; 185 G4LorentzVector Qmomentum;
186 G4double Qminus, Qplus; 186 G4double Qminus, Qplus;
187 187
188 G4int whilecount=0; 188 G4int whilecount=0;
189 do { 189 do {
190 whilecount++; 190 whilecount++;
191 191
192 if (whilecount > 1000 ) 192 if (whilecount > 1000 )
193 { 193 {
194 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 194 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
195 return false; // Ignore this intera 195 return false; // Ignore this interaction
196 } 196 }
197 197
198 // Generate pt 198 // Generate pt
199 Qmomentum=G4LorentzVector(GaussianPt(Avera 199 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
200 200
201 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma 201 Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
202 202
203 ProjMassT2 = Mprojectile2 + Pt2; 203 ProjMassT2 = Mprojectile2 + Pt2;
204 ProjMassT = std::sqrt( ProjMassT2 ); 204 ProjMassT = std::sqrt( ProjMassT2 );
205 TargMassT2 = Mtarget2 + Pt2; 205 TargMassT2 = Mtarget2 + Pt2;
206 TargMassT = std::sqrt( TargMassT2 ); 206 TargMassT = std::sqrt( TargMassT2 );
207 207
208 #ifdef debugSingleDiffraction 208 #ifdef debugSingleDiffraction
209 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjM 209 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl;
210 #endif 210 #endif
211 if ( SqrtS < ProjMassT + TargMassT ) conti 211 if ( SqrtS < ProjMassT + TargMassT ) continue;
212 212
213 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + T 213 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
214 - 2.0*S*ProjMassT2 - 2.0*S*TargM 214 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
215 215
216 if ( PZcms2 < 0 ) continue; 216 if ( PZcms2 < 0 ) continue;
217 217
218 PZcms = std::sqrt( PZcms2 ); 218 PZcms = std::sqrt( PZcms2 );
219 219
220 if ( ProjectileDiffraction ) 220 if ( ProjectileDiffraction )
221 { // The projectile will fragment, t 221 { // The projectile will fragment, the target will saved.
222 PMinusMin = std::sqrt( ProjMassT2 + PZcm 222 PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
223 PMinusMax = SqrtS - TargMassT; 223 PMinusMax = SqrtS - TargMassT;
224 224
225 PMinusNew = ChooseX( PMinusMin, PMinusMa 225 PMinusNew = ChooseX( PMinusMin, PMinusMax );
226 TMinusNew = SqrtS - PMinusNew; 226 TMinusNew = SqrtS - PMinusNew;
227 227
228 Qminus = Ptarget.minus() - TMinusNew; 228 Qminus = Ptarget.minus() - TMinusNew;
229 TPlusNew = TargMassT2 / TMinusNew; 229 TPlusNew = TargMassT2 / TMinusNew;
230 Qplus = Ptarget.plus() - TPlusNew; 230 Qplus = Ptarget.plus() - TPlusNew;
231 231
232 } else { // The target will fragment, the 232 } else { // The target will fragment, the projectile will saved.
233 TPlusMin = std::sqrt( TargMassT2 + PZcms 233 TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
234 TPlusMax = SqrtS - ProjMassT; 234 TPlusMax = SqrtS - ProjMassT;
235 235
236 TPlusNew = ChooseX( TPlusMin, TPlusMax ) 236 TPlusNew = ChooseX( TPlusMin, TPlusMax );
237 PPlusNew = SqrtS - TPlusNew; 237 PPlusNew = SqrtS - TPlusNew;
238 238
239 Qplus = PPlusNew - Pprojectile.plus(); 239 Qplus = PPlusNew - Pprojectile.plus();
240 PMinusNew = ProjMassT2 / PPlusNew; 240 PMinusNew = ProjMassT2 / PPlusNew;
241 Qminus = PMinusNew - Pprojectile.minus() 241 Qminus = PMinusNew - Pprojectile.minus();
242 } 242 }
243 243
244 Qmomentum.setPz( (Qplus - Qminus)/2 ); 244 Qmomentum.setPz( (Qplus - Qminus)/2 );
245 Qmomentum.setE( (Qplus + Qminus)/2 ); 245 Qmomentum.setE( (Qplus + Qminus)/2 );
246 246
247 #ifdef debugSingleDiffraction 247 #ifdef debugSingleDiffraction
248 G4cout<<ProjectileDiffraction<<" "<<( Pp 248 G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl;
249 G4cout<<!ProjectileDiffraction<<" "<<( P 249 G4cout<<!ProjectileDiffraction<<" "<<( Ptarget - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl;
250 #endif 250 #endif
251 251
252 } while ( ( ProjectileDiffraction&&( Pprojec 252 } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) ||
253 (!ProjectileDiffraction&&( Ptarget 253 (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) );
254 // Repeat the sampling because there was n 254 // Repeat the sampling because there was not any excitation
255 255
256 Pprojectile += Qmomentum; 256 Pprojectile += Qmomentum;
257 257
258 Ptarget -= Qmomentum; 258 Ptarget -= Qmomentum;
259 259
260 // Transform back and update SplitableHadron 260 // Transform back and update SplitableHadron Participant.
261 Pprojectile.transform(toLab); 261 Pprojectile.transform(toLab);
262 Ptarget.transform(toLab); 262 Ptarget.transform(toLab);
263 263
264 #ifdef debugSingleDiffraction 264 #ifdef debugSingleDiffraction
265 G4cout << "Pprojectile in Lab. " << Pproj 265 G4cout << "Pprojectile in Lab. " << Pprojectile << G4endl;
266 G4cout << "Ptarget in Lab. " << Ptarg 266 G4cout << "Ptarget in Lab. " << Ptarget << G4endl;
267 G4cout << "G4SingleDiffractiveExcitation- 267 G4cout << "G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() << G4endl;
268 G4cout << "G4SingleDiffractiveExcitation- 268 G4cout << "G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() << G4endl;
269 #endif 269 #endif
270 270
271 target->Set4Momentum(Ptarget); 271 target->Set4Momentum(Ptarget);
272 projectile->Set4Momentum(Pprojectile); 272 projectile->Set4Momentum(Pprojectile);
273 273
274 return true; 274 return true;
275 } 275 }
276 276
277 // --------- private methods ----------------- 277 // --------- private methods ----------------------
278 278
279 G4double G4SingleDiffractiveExcitation::Choose 279 G4double G4SingleDiffractiveExcitation::ChooseX(G4double Xmin, G4double Xmax) const
280 { 280 {
281 // choose an x between Xmin and Xmax with P( 281 // choose an x between Xmin and Xmax with P(x) ~ 1/x
282 G4double range=Xmax-Xmin; 282 G4double range=Xmax-Xmin;
283 283
284 if ( Xmin <= 0. || range <=0. ) 284 if ( Xmin <= 0. || range <=0. )
285 { 285 {
286 G4cout << " Xmin, range : " << Xmin << " , 286 G4cout << " Xmin, range : " << Xmin << " , " << range << G4endl;
287 throw G4HadronicException(__FILE__, __LINE 287 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation::ChooseX : Invalid arguments ");
288 } 288 }
289 289
290 G4double x = Xmin*G4Pow::GetInstance()->powA 290 G4double x = Xmin*G4Pow::GetInstance()->powA(Xmax/Xmin, G4UniformRand() );
291 // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin) 291 // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin) - G4UniformRand() * (1.0/std::sqrt(Xmin) - 1.0/std::sqrt(Xmax)));
292 return x; 292 return x;
293 } 293 }
294 294
295 295
296 G4ThreeVector G4SingleDiffractiveExcitation::G 296 G4ThreeVector G4SingleDiffractiveExcitation::GaussianPt(G4double widthSquare, G4double maxPtSquare) const
297 { // @@ this method is used in FTF 297 { // @@ this method is used in FTFModel as well. Should go somewhere common!
298 298
299 G4double pt2; 299 G4double pt2;
300 300
301 const G4int maxNumberOfLoops = 1000; 301 const G4int maxNumberOfLoops = 1000;
302 G4int loopCounter = 0; 302 G4int loopCounter = 0;
303 do { 303 do {
304 pt2=-widthSquare * G4Log( G4UniformRand() 304 pt2=-widthSquare * G4Log( G4UniformRand() );
305 } while ( ( pt2 > maxPtSquare) && ++loopCoun 305 } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
306 if ( loopCounter >= maxNumberOfLoops ) { 306 if ( loopCounter >= maxNumberOfLoops ) {
307 pt2 = 0.99*maxPtSquare; // Just an accept 307 pt2 = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration.
308 } 308 }
309 309
310 pt2=std::sqrt(pt2); 310 pt2=std::sqrt(pt2);
311 311
312 G4double phi=G4UniformRand() * twopi; 312 G4double phi=G4UniformRand() * twopi;
313 313
314 return G4ThreeVector (pt2*std::cos(phi), pt2 314 return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);
315 } 315 }
316 316
317 317