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