| 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 // ---------------- G4QGSDiffractiveExcit 30 // ---------------- G4QGSDiffractiveExcitation --------------
31 // by Gunter Folger, October 1998. 31 // by Gunter Folger, October 1998.
32 // diffractive Excitation used by stri 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 // Essential changed by V. Uzhinsky in Novemb 35 // Essential changed by V. Uzhinsky in November - December 2006
36 // in order to put it in a correspondence wit 36 // in order to put it in a correspondence with original FRITIOF
37 // model. Variant of FRITIOF with nucleon de- 37 // model. Variant of FRITIOF with nucleon de-excitation is implemented.
38 // ------------------------------------------- 38 // ---------------------------------------------------------------------
39 39
40 // Modified: 40 // Modified:
41 // 25-05-07 : G.Folger 41 // 25-05-07 : G.Folger
42 // move from management/G4DiffractiveExc 42 // move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation
43 // 43 //
44 44
45 #include "globals.hh" 45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh" 48 #include "Randomize.hh"
49 49
50 #include "G4QGSDiffractiveExcitation.hh" 50 #include "G4QGSDiffractiveExcitation.hh"
51 #include "G4LorentzRotation.hh" 51 #include "G4LorentzRotation.hh"
52 #include "G4ThreeVector.hh" 52 #include "G4ThreeVector.hh"
53 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh"
54 #include "G4VSplitableHadron.hh" 54 #include "G4VSplitableHadron.hh"
55 #include "G4ExcitedString.hh" 55 #include "G4ExcitedString.hh"
56 //#include "G4ios.hh" 56 //#include "G4ios.hh"
57 57
58 #include "G4Exp.hh" 58 #include "G4Exp.hh"
59 #include "G4Log.hh" 59 #include "G4Log.hh"
60 #include "G4Pow.hh" 60 #include "G4Pow.hh"
61 61
62 //============================================ 62 //============================================================================
63 63
64 //#define debugDoubleDiffraction 64 //#define debugDoubleDiffraction
65 65
66 //============================================ 66 //============================================================================
67 67
68 G4QGSDiffractiveExcitation::G4QGSDiffractiveEx 68 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation()
69 { 69 {
70 } 70 }
71 71
72 G4QGSDiffractiveExcitation::~G4QGSDiffractiveE 72 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation()
73 { 73 {
74 } 74 }
75 75
76 76
77 G4bool G4QGSDiffractiveExcitation:: 77 G4bool G4QGSDiffractiveExcitation::
78 ExciteParticipants(G4VSplitableHadron *project 78 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4bool ) const
79 { 79 {
80 #ifdef debugDoubleDiffraction 80 #ifdef debugDoubleDiffraction
81 G4cout<<G4endl<<"G4QGSDiffractiveExcitatio 81 G4cout<<G4endl<<"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<G4endl;
82 G4cout<<"Proj Targ "<<projectile->GetDefin 82 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetParticleName()<<" "<<target->GetDefinition()->GetParticleName()<<G4endl;
83 G4cout<<"Proj 4 Mom "<<projectile->Get4Mom 83 G4cout<<"Proj 4 Mom "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
84 G4cout<<"Targ 4 Mom "<<target->Get4Momentu 84 G4cout<<"Targ 4 Mom "<<target->Get4Momentum() <<" "<<target->Get4Momentum().mag() <<G4endl;
85 #endif 85 #endif
86 86
87 G4LorentzVector Pprojectile=projectile->Get4 87 G4LorentzVector Pprojectile=projectile->Get4Momentum();
88 88
89 // -------------------- Projectile parameter 89 // -------------------- Projectile parameters -----------------------------------
90 G4bool PutOnMassShell=0; 90 G4bool PutOnMassShell=0;
91 91
92 G4double M0projectile = Pprojectile.mag(); 92 G4double M0projectile = Pprojectile.mag(); // Without de-excitation
93 93
94 if(M0projectile < projectile->GetDefinition( 94 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
95 { 95 {
96 PutOnMassShell=1; 96 PutOnMassShell=1;
97 M0projectile=projectile->GetDefinition()-> 97 M0projectile=projectile->GetDefinition()->GetPDGMass();
98 } 98 }
99 99
100 // -------------------- Target parameters -- 100 // -------------------- Target parameters ----------------------------------------------
101 G4LorentzVector Ptarget=target->Get4Momentum 101 G4LorentzVector Ptarget=target->Get4Momentum();
102 102
103 G4double M0target = Ptarget.mag(); 103 G4double M0target = Ptarget.mag();
104 104
105 if(M0target < target->GetDefinition()->GetPD 105 if(M0target < target->GetDefinition()->GetPDGMass())
106 { 106 {
107 PutOnMassShell=1; 107 PutOnMassShell=1;
108 M0target=target->GetDefinition()->GetPDGMa 108 M0target=target->GetDefinition()->GetPDGMass();
109 } 109 }
110 110
111 G4LorentzVector Psum=Pprojectile+Ptarget; 111 G4LorentzVector Psum=Pprojectile+Ptarget;
112 G4double S=Psum.mag2(); 112 G4double S=Psum.mag2();
113 G4double SqrtS=std::sqrt(S); 113 G4double SqrtS=std::sqrt(S);
114 114
115 if (SqrtS < M0projectile + M0target) {return 115 if (SqrtS < M0projectile + M0target) {return false;} // The model cannot work for pp-interactions
116 116 // at Plab < 1.3 GeV/c.
117 117
118 G4double Mprojectile2 = M0projectile * M0pro 118 G4double Mprojectile2 = M0projectile * M0projectile;
119 G4double Mtarget2 = M0target * M0tar 119 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
120 120
121 // Transform momenta to cms and then rotate 121 // Transform momenta to cms and then rotate parallel to z axis;
122 122
123 G4LorentzRotation toCms(-1*Psum.boostVector( 123 G4LorentzRotation toCms(-1*Psum.boostVector());
124 124
125 G4LorentzVector Ptmp=toCms*Pprojectile; 125 G4LorentzVector Ptmp=toCms*Pprojectile;
126 126
127 if ( Ptmp.pz() <= 0. ) 127 if ( Ptmp.pz() <= 0. )
128 { 128 {
129 // "String" moving backwards in CMS, abor 129 // "String" moving backwards in CMS, abort collision !!
130 //G4cout << " abort Collision!! " << G4end 130 //G4cout << " abort Collision!! " << G4endl;
131 return false; 131 return false;
132 } 132 }
133 133
134 toCms.rotateZ(-1*Ptmp.phi()); 134 toCms.rotateZ(-1*Ptmp.phi());
135 toCms.rotateY(-1*Ptmp.theta()); 135 toCms.rotateY(-1*Ptmp.theta());
136 136
137 G4LorentzRotation toLab(toCms.inverse()); 137 G4LorentzRotation toLab(toCms.inverse());
138 138
139 Pprojectile.transform(toCms); 139 Pprojectile.transform(toCms);
140 Ptarget.transform(toCms); 140 Ptarget.transform(toCms);
141 141
142 G4double PZcms2=(S*S+Mprojectile2*Mprojectil 142 G4double PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
143 2*S*Mprojectile2-2*S*Mtarget2-2*Mproj 143 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
144 144
145 if (PZcms2 < 0) {return false;} // It can 145 if (PZcms2 < 0) {return false;} // It can be in an interaction with off-shell nuclear nucleon
146 146
147 G4double PZcms = std::sqrt(PZcms2); 147 G4double PZcms = std::sqrt(PZcms2);
148 148
149 if (PutOnMassShell) 149 if (PutOnMassShell)
150 { 150 {
151 if (Pprojectile.z() > 0.) 151 if (Pprojectile.z() > 0.)
152 { 152 {
153 Pprojectile.setPz( PZcms); 153 Pprojectile.setPz( PZcms);
154 Ptarget.setPz( -PZcms); 154 Ptarget.setPz( -PZcms);
155 } 155 }
156 else 156 else
157 { 157 {
158 Pprojectile.setPz(-PZcms); 158 Pprojectile.setPz(-PZcms);
159 Ptarget.setPz( PZcms); 159 Ptarget.setPz( PZcms);
160 }; 160 };
161 161
162 Pprojectile.setE(std::sqrt(Mprojectile2+sq 162 Pprojectile.setE(std::sqrt(Mprojectile2+sqr(Pprojectile.x())+sqr(Pprojectile.y())+PZcms2));
163 Ptarget.setE( std::sqrt( Mtarget2+sq 163 Ptarget.setE( std::sqrt( Mtarget2+sqr( Ptarget.x())+sqr( Ptarget.y())+PZcms2));
164 } 164 }
165 165
166 G4double maxPtSquare = PZcms2; 166 G4double maxPtSquare = PZcms2;
167 167
168 #ifdef debugDoubleDiffraction 168 #ifdef debugDoubleDiffraction
169 G4cout << "Pprojectile after boost to CMS: 169 G4cout << "Pprojectile after boost to CMS: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
170 G4cout << "Ptarget after boost to CMS: 170 G4cout << "Ptarget after boost to CMS: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
171 #endif 171 #endif
172 172
173 G4int PrPDGcode=projectile->GetDefinition 173 G4int PrPDGcode=projectile->GetDefinition()->GetPDGEncoding();
174 G4int absPrPDGcode=std::abs(PrPDGcode); 174 G4int absPrPDGcode=std::abs(PrPDGcode);
175 G4double MinPrDiffMass(0.); 175 G4double MinPrDiffMass(0.);
176 G4double AveragePt2(0.); 176 G4double AveragePt2(0.);
177 177
178 if (M0projectile <= projectile->GetDefinitio 178 if (M0projectile <= projectile->GetDefinition()->GetPDGMass())
179 { // Normal projectile 179 { // Normal projectile
180 if( absPrPDGcode > 1000 ) 180 if( absPrPDGcode > 1000 ) //------Projectile is baryon --------
181 { 181 {
182 if ( absPrPDGcode > 4000 && absPrPDGcode 182 if ( absPrPDGcode > 4000 && absPrPDGcode < 6000 ) // Projectile is a charm or bottom baryon
183 { 183 {
184 MinPrDiffMass = projectile->GetDefinit 184 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV
185 AveragePt2 = 0.3; 185 AveragePt2 = 0.3; // GeV^2
186 } 186 }
187 else 187 else
188 { 188 {
189 MinPrDiffMass = 1.16; // 189 MinPrDiffMass = 1.16; // GeV
190 AveragePt2 = 0.3; 190 AveragePt2 = 0.3; // GeV^2
191 } 191 }
192 } 192 }
193 else if( absPrPDGcode == 211 || PrPDGcode 193 else if( absPrPDGcode == 211 || PrPDGcode == 111) //------Projectile is Pion -----------
194 { 194 {
195 MinPrDiffMass = 1.0; 195 MinPrDiffMass = 1.0; // GeV
196 AveragePt2 = 0.3; 196 AveragePt2 = 0.3; // GeV^2
197 } 197 }
198 else if( absPrPDGcode == 321 || absPrPDGco 198 else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310) //-Projectile is Kaon-
199 { 199 {
200 MinPrDiffMass = 1.1; 200 MinPrDiffMass = 1.1; // GeV
201 AveragePt2 = 0.3; 201 AveragePt2 = 0.3; // GeV^2
202 } 202 }
203 else if( absPrPDGcode > 400 && absPrPDGcod 203 else if( absPrPDGcode > 400 && absPrPDGcode < 600) // Projectile is a charm or bottom meson
204 { 204 {
205 MinPrDiffMass = projectile->GetDefinitio 205 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV
206 AveragePt2 = 0.3; 206 AveragePt2 = 0.3; // GeV^2
207 } 207 }
208 else 208 else //------Projectile is undefined, Nucleon assumed
209 { 209 {
210 MinPrDiffMass = 1.16; 210 MinPrDiffMass = 1.16; // GeV
211 AveragePt2 = 0.3; 211 AveragePt2 = 0.3; // GeV^2
212 } 212 }
213 } 213 }
214 else 214 else
215 { // Excited projectile 215 { // Excited projectile
216 MinPrDiffMass = M0projectile + 220.0*MeV; 216 MinPrDiffMass = M0projectile + 220.0*MeV;
217 AveragePt2 = 0.3; 217 AveragePt2 = 0.3;
218 } 218 }
219 219
220 MinPrDiffMass = MinPrDiffMass * GeV; 220 MinPrDiffMass = MinPrDiffMass * GeV;
221 AveragePt2 = AveragePt2 * GeV*GeV; 221 AveragePt2 = AveragePt2 * GeV*GeV;
222 //------------------------------------------ 222 //---------------------------------------------
223 G4double MinTrDiffMass = 1.16*GeV; 223 G4double MinTrDiffMass = 1.16*GeV;
224 224
225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) { 225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) {return false;} // The model cannot work at low energy
226 226
227 G4double MinPrDiffMass2 = MinPrDiffMass * Mi 227 G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
228 G4double MinTrDiffMass2 = MinTrDiffMass * Mi 228 G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
229 229
230 G4double Pt2; 230 G4double Pt2;
231 G4double ProjMassT2, ProjMassT; 231 G4double ProjMassT2, ProjMassT;
232 G4double TargMassT2, TargMassT; 232 G4double TargMassT2, TargMassT;
233 G4double PMinusNew, TPlusNew; 233 G4double PMinusNew, TPlusNew;
234 234
235 G4LorentzVector Qmomentum; 235 G4LorentzVector Qmomentum;
236 G4double Qminus, Qplus; 236 G4double Qminus, Qplus;
237 237
238 G4int whilecount=0; 238 G4int whilecount=0;
239 do { 239 do {
240 if (whilecount++ >= 500 && (whilecount%100 240 if (whilecount++ >= 500 && (whilecount%100)==0)
241 if (whilecount > 1000 ) { 241 if (whilecount > 1000 ) {
242 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 242 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
243 return false; // Ignore this interactio 243 return false; // Ignore this interaction
244 } 244 }
245 245
246 // Generate pt 246 // Generate pt
247 Qmomentum=G4LorentzVector(GaussianPt(Avera 247 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
248 248
249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2() 249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
250 ProjMassT2=MinPrDiffMass2+Pt2; 250 ProjMassT2=MinPrDiffMass2+Pt2;
251 ProjMassT =std::sqrt(ProjMassT2); 251 ProjMassT =std::sqrt(ProjMassT2);
252 252
253 TargMassT2=MinTrDiffMass2+Pt2; 253 TargMassT2=MinTrDiffMass2+Pt2;
254 TargMassT =std::sqrt(TargMassT2); 254 TargMassT =std::sqrt(TargMassT2);
255 255
256 if (SqrtS < ProjMassT + TargMassT) continu 256 if (SqrtS < ProjMassT + TargMassT) continue;
257 257
258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMass 258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMassT2*TargMassT2-
259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjM 259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)/4./S;
260 if (PZcms2 < 0 ) {PZcms2=0;}; 260 if (PZcms2 < 0 ) {PZcms2=0;};
261 PZcms =std::sqrt(PZcms2); 261 PZcms =std::sqrt(PZcms2);
262 262
263 G4double PMinusMin=std::sqrt(ProjMassT2+PZ 263 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
264 G4double PMinusMax=SqrtS-TargMassT; 264 G4double PMinusMax=SqrtS-TargMassT;
265 265
266 PMinusNew=ChooseP(PMinusMin,PMinusMax); 266 PMinusNew=ChooseP(PMinusMin,PMinusMax);
267 Qminus=PMinusNew-Pprojectile.minus(); 267 Qminus=PMinusNew-Pprojectile.minus();
268 268
269 G4double TPlusMin=std::sqrt(TargMassT2+PZc 269 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
270 G4double TPlusMax=SqrtS-ProjMassT; 270 G4double TPlusMax=SqrtS-ProjMassT;
271 271
272 TPlusNew=ChooseP(TPlusMin, TPlusMax); 272 TPlusNew=ChooseP(TPlusMin, TPlusMax);
273 Qplus=-(TPlusNew-Ptarget.plus()); 273 Qplus=-(TPlusNew-Ptarget.plus());
274 274
275 Qmomentum.setPz( (Qplus-Qminus)/2 ); 275 Qmomentum.setPz( (Qplus-Qminus)/2 );
276 Qmomentum.setE( (Qplus+Qminus)/2 ); 276 Qmomentum.setE( (Qplus+Qminus)/2 );
277 277
278 } while ( (Pprojectile+Qmomentum).mag2() < 278 } while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 ||
279 (Ptarget -Qmomentum).mag2() < MinTrD 279 (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
280 280
281 Pprojectile += Qmomentum; 281 Pprojectile += Qmomentum;
282 Ptarget -= Qmomentum; 282 Ptarget -= Qmomentum;
283 283
284 // Transform back and update SplitableHadron 284 // Transform back and update SplitableHadron Participant.
285 Pprojectile.transform(toLab); 285 Pprojectile.transform(toLab);
286 Ptarget.transform(toLab); 286 Ptarget.transform(toLab);
287 287
288 #ifdef debugDoubleDiffraction 288 #ifdef debugDoubleDiffraction
289 G4cout << "Pprojectile after boost to Lab: 289 G4cout << "Pprojectile after boost to Lab: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
290 G4cout << "Ptarget after boost to Lab: 290 G4cout << "Ptarget after boost to Lab: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
291 #endif 291 #endif
292 292
293 target->Set4Momentum(Ptarget); 293 target->Set4Momentum(Ptarget);
294 projectile->Set4Momentum(Pprojectile); 294 projectile->Set4Momentum(Pprojectile);
295 295
296 return true; 296 return true;
297 } 297 }
298 298
299 299
300 G4ExcitedString * G4QGSDiffractiveExcitation:: 300 G4ExcitedString * G4QGSDiffractiveExcitation::
301 String(G4VSplitableHadron * hadron, G4bool isP 301 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
302 { 302 {
303 hadron->SplitUp(); 303 hadron->SplitUp();
304 G4Parton *start= hadron->GetNextParton(); 304 G4Parton *start= hadron->GetNextParton();
305 if ( start==NULL) 305 if ( start==NULL)
306 { G4cout << " G4QGSDiffractiveExcitation::St 306 { G4cout << " G4QGSDiffractiveExcitation::String() Error:No start parton found"<< G4endl;
307 return NULL; 307 return NULL;
308 } 308 }
309 G4Parton *end = hadron->GetNextParton(); 309 G4Parton *end = hadron->GetNextParton();
310 if ( end==NULL) 310 if ( end==NULL)
311 { G4cout << " G4QGSDiffractiveExcitation::St 311 { G4cout << " G4QGSDiffractiveExcitation::String() Error:No end parton found"<< G4endl;
312 return NULL; 312 return NULL;
313 } 313 }
314 314
315 G4ExcitedString * string; 315 G4ExcitedString * string;
316 if ( isProjectile ) 316 if ( isProjectile )
317 { 317 {
318 string= new G4ExcitedString(end,start, +1) 318 string= new G4ExcitedString(end,start, +1);
319 } else { 319 } else {
320 string= new G4ExcitedString(start,end, -1) 320 string= new G4ExcitedString(start,end, -1);
321 } 321 }
322 322
323 string->SetPosition(hadron->GetPosition()); 323 string->SetPosition(hadron->GetPosition());
324 324
325 // momenta of string ends 325 // momenta of string ends
326 326
327 G4double maxAvailMomentumSquared=sqr(hadron- 327 G4double maxAvailMomentumSquared=sqr(hadron->Get4Momentum().mag()/2.);
328 328
329 G4double widthOfPtSquare = 0.5*sqr(GeV); 329 G4double widthOfPtSquare = 0.5*sqr(GeV);
330 G4ThreeVector pt=GaussianPt(widthOfPtSquare, 330 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
331 331
332 G4LorentzVector Pstart(G4LorentzVector(pt,0. 332 G4LorentzVector Pstart(G4LorentzVector(pt,0.));
333 G4LorentzVector Pend; 333 G4LorentzVector Pend;
334 Pend.setPx(hadron->Get4Momentum().px() - pt. 334 Pend.setPx(hadron->Get4Momentum().px() - pt.x());
335 Pend.setPy(hadron->Get4Momentum().py() - pt. 335 Pend.setPy(hadron->Get4Momentum().py() - pt.y());
336 336
337 G4double tm1=hadron->Get4Momentum().minus() 337 G4double tm1=hadron->Get4Momentum().minus() +
338 ( Pend.perp2()-Pstart.perp2() ) / had 338 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
339 339
340 G4double tm2= std::sqrt( std::max(0., sqr(tm 340 G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
341 4. * Pend.perp2() * hadron->Get4Momen 341 4. * Pend.perp2() * hadron->Get4Momentum().minus()
342 / hadron->Get4Momentum().plus() )); 342 / hadron->Get4Momentum().plus() ));
343 343
344 G4int Sign= isProjectile ? -1 : 1; 344 G4int Sign= isProjectile ? -1 : 1;
345 345
346 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 346 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
347 G4double startMinus= hadron->Get4Momentum(). 347 G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
348 348
349 G4double startPlus= Pstart.perp2() / startM 349 G4double startPlus= Pstart.perp2() / startMinus;
350 G4double endPlus = hadron->Get4Momentum().p 350 G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
351 351
352 Pstart.setPz(0.5*(startPlus - startMinus)); 352 Pstart.setPz(0.5*(startPlus - startMinus));
353 Pstart.setE(0.5*(startPlus + startMinus)); 353 Pstart.setE(0.5*(startPlus + startMinus));
354 354
355 Pend.setPz(0.5*(endPlus - endMinus)); 355 Pend.setPz(0.5*(endPlus - endMinus));
356 Pend.setE(0.5*(endPlus + endMinus)); 356 Pend.setE(0.5*(endPlus + endMinus));
357 357
358 start->Set4Momentum(Pstart); 358 start->Set4Momentum(Pstart);
359 end->Set4Momentum(Pend); 359 end->Set4Momentum(Pend);
360 360
361 #ifdef debugQGSdiffExictation 361 #ifdef debugQGSdiffExictation
362 G4cout << " generated string flavors 362 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
363 G4cout << " generated string momenta: qu 363 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
364 G4cout << " generated string momenta: Diqu 364 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
365 G4cout << " sum of ends 365 G4cout << " sum of ends " << Pstart+Pend << G4endl;
366 G4cout << " Original 366 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
367 #endif 367 #endif
368 368
369 return string; 369 return string;
370 } 370 }
371 371
372 372
373 // --------- private methods ----------------- 373 // --------- private methods ----------------------
374 374
375 G4double G4QGSDiffractiveExcitation::ChooseP(G 375 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
376 { 376 {
377 // choose an x between Xmin and Xmax with P( 377 // choose an x between Xmin and Xmax with P(x) ~ 1/x
378 // to be improved... 378 // to be improved...
379 379
380 G4double range=Pmax-Pmin; 380 G4double range=Pmax-Pmin;
381 381
382 if ( Pmin <= 0. || range <=0. ) 382 if ( Pmin <= 0. || range <=0. )
383 { 383 {
384 G4cout << " Pmin, range : " << Pmin << " , 384 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
385 throw G4HadronicException(__FILE__, __LINE 385 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
386 } 386 }
387 387
388 G4double P; 388 G4double P;
389 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmi 389 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmin,G4UniformRand());
390 //debug-hpw cout << "DiffractiveX "<<x<<G4en 390 //debug-hpw cout << "DiffractiveX "<<x<<G4endl;
391 return P; 391 return P;
392 } 392 }
393 393
394 394
395 G4ThreeVector G4QGSDiffractiveExcitation::Gaus 395 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
396 { // @@ this method is used in FTF 396 { // @@ this method is used in FTFModel as well. Should go somewhere common!
397 G4double Pt2; 397 G4double Pt2;
398 398
399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand 399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.));
400 400
401 G4double Pt=std::sqrt(Pt2); 401 G4double Pt=std::sqrt(Pt2);
402 402
403 G4double phi=G4UniformRand() * twopi; 403 G4double phi=G4UniformRand() * twopi;
404 404
405 return G4ThreeVector (Pt*std::cos(phi), Pt*s 405 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
406 } 406 }
407 407