| 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
29 // -------------------------------------------
30 // GEANT 4 class implemetation file
31 //
32 // ---------------- G4DiffractiveExcitati
33 // by Gunter Folger, October 1998.
34 // diffractive Excitation used by strin
35 // Take a projectile and a target
36 // excite the projectile and targe
37 // Essential changed by V. Uzhinsky in Novemb
38 // in order to put it in a correspondence wit
39 // model. Variant of FRITIOF with nucleon de-
40 // Other changes by V.Uzhinsky in May 2007 we
41 // meson-nucleon interactions. Additional cha
42 // were introduced in December 2006. They tre
43 // processes more exactly.
44 // Correct treatment of the diffraction disso
45 // Mass distributions for resonances and uu-d
46 // and dd-diquarks suppression in neutrons we
47 // -------------------------------------------
48
49 #include "globals.hh"
50 #include "Randomize.hh"
51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh"
53
54 #include "G4DiffractiveExcitation.hh"
55 #include "G4FTFParameters.hh"
56 #include "G4ElasticHNScattering.hh"
57
58 #include "G4RotationMatrix.hh"
59 #include "G4ParticleDefinition.hh"
60 #include "G4ParticleTable.hh"
61 #include "G4SampleResonance.hh"
62 #include "G4VSplitableHadron.hh"
63 #include "G4ExcitedString.hh"
64 #include "G4Neutron.hh"
65
66 #include "G4Exp.hh"
67 #include "G4Log.hh"
68 #include "G4Pow.hh"
69
70 //#include "G4ios.hh"
71
72 //============================================
73
74 //#define debugFTFexcitation
75 //#define debug_heavyHadrons
76
77 //============================================
78
79 G4DiffractiveExcitation::G4DiffractiveExcitati
80
81
82 //============================================
83
84 G4DiffractiveExcitation::~G4DiffractiveExcitat
85
86
87 //============================================
88
89 G4bool G4DiffractiveExcitation::ExciteParticip
90
91
92
93
94 #ifdef debugFTFexcitation
95 G4cout << G4endl << "FTF ExciteParticipants
96 #endif
97
98 CommonVariables common;
99
100 // Projectile parameters
101 common.Pprojectile = projectile->Get4Momentu
102 if ( common.Pprojectile.z() < 0.0 ) return f
103 common.ProjectilePDGcode = projectile->GetDe
104 common.absProjectilePDGcode = std::abs( comm
105 common.M0projectile = projectile->GetDefinit
106 G4double ProjectileRapidity = common.Pprojec
107
108 // Target parameters
109 common.Ptarget = target->Get4Momentum();
110 common.TargetPDGcode = target->GetDefinition
111 common.absTargetPDGcode = std::abs( common.T
112 common.M0target = target->GetDefinition()->G
113 G4double TargetRapidity = common.Ptarget.rap
114
115 // Kinematical properties of the interaction
116 G4LorentzVector Psum = common.Pprojectile +
117 common.S = Psum.mag2();
118 common.SqrtS = std::sqrt( common.S );
119
120 // Check off-shellness of the participants
121 G4bool toBePutOnMassShell = true; //Uzhi Au
122 common.MminProjectile = common.BrW.GetMinimu
123 /* Uzhi Aug.2019
124 if ( common.M0projectile < common.MminProjec
125 toBePutOnMassShell = true;
126 common.M0projectile = common.BrW.SampleMas
127
128
129 }
130 */
131 common.M0projectile2 = common.M0projectile *
132 common.ProjectileDiffStateMinMass = thePa
133 common.ProjectileNonDiffStateMinMass = thePa
134 if ( common.M0projectile > common.Projectile
135 common.ProjectileDiffStateMinMass = com
136 common.ProjectileNonDiffStateMinMass = com
137 if ( common.absProjectilePDGcode > 3000 )
138 common.ProjectileDiffStateMinMass +=
139 common.ProjectileNonDiffStateMinMass +=
140 }
141 }
142 common.MminTarget = common.BrW.GetMinimumMas
143 /* Uzhi Aug.2019
144 if ( common.M0target < common.MminTarget ) {
145 toBePutOnMassShell = true;
146 common.M0target = common.BrW.SampleMass( t
147 t
148 +
149 }
150 */
151 common.M0target2 = common.M0target * common.
152 common.TargetDiffStateMinMass = theParame
153 common.TargetNonDiffStateMinMass = theParame
154 if ( common.M0target > common.TargetDiffStat
155 common.TargetDiffStateMinMass = common.
156 common.TargetNonDiffStateMinMass = common.
157 if ( common.absTargetPDGcode > 3000 ) { /
158 common.TargetDiffStateMinMass += 140.
159 common.TargetNonDiffStateMinMass += 140.
160 }
161 };
162 #ifdef debugFTFexcitation
163 G4cout << "Proj Targ PDGcodes " << common.Pr
164 << "Mprojectile Y " << common.Pproje
165 << "M0projectile Y " << common.M0proj
166 G4cout << "Mtarget Y " << common.Ptarge
167 << "M0target Y " << common.M0targ
168 G4cout << "Pproj " << common.Pprojectile <<
169 #endif
170
171 // Transform momenta to cms and then rotate
172 common.toCms = G4LorentzRotation( -1 * Psum.
173 G4LorentzVector Ptmp = common.toCms * common
174 if ( Ptmp.pz() <= 0.0 ) return false; // "S
175 common.toCms.rotateZ( -1*Ptmp.phi() );
176 common.toCms.rotateY( -1*Ptmp.theta() );
177 common.toLab = common.toCms.inverse();
178 common.Pprojectile.transform( common.toCms )
179 common.Ptarget.transform( common.toCms );
180
181 G4double SumMasses = common.M0projectile + c
182 #ifdef debugFTFexcitation
183 G4cout << "SqrtS " << common.SqrtS << G4
184 << " " << common.M0target << " " << S
185 #endif
186 if ( common.SqrtS < SumMasses ) return false
187
188 common.PZcms2 = ( sqr( common.S ) + sqr( com
189 - 2.0 * ( common.S * ( com
190 + common.M0proje
191 #ifdef debugFTFexcitation
192 G4cout << "PZcms2 after toBePutOnMassShell "
193 #endif
194 if ( common.PZcms2 < 0.0 ) return false; //
195
196 common.PZcms = std::sqrt( common.PZcms2 );
197 if ( toBePutOnMassShell ) {
198 if ( common.Pprojectile.z() > 0.0 ) {
199 common.Pprojectile.setPz( common.PZcms
200 common.Ptarget.setPz( -common.PZcms
201 } else {
202 common.Pprojectile.setPz( -common.PZcms
203 common.Ptarget.setPz( common.PZcms
204 }
205 common.Pprojectile.setE( std::sqrt( common
206 + comm
207 + comm
208 + comm
209 common.Ptarget.setE( std::sqrt( common.M0t
210 + common.P
211 + common.P
212 + common.P
213 }
214 #ifdef debugFTFexcitation
215 G4cout << "Start --------------------" << G4
216 << " " << common.ProjectileDiffStateM
217 << G4endl
218 << "Targ M0 Mdif Mndif " << common.M0
219 << " " << common.TargetNonDiffStateMi
220 << "Proj CMS " << common.Pprojectile
221 #endif
222
223 // Check for possible quark exchange
224 ProjectileRapidity = common.Pprojectile.rapi
225 TargetRapidity = common.Ptarget.rapidity();
226 G4double QeNoExc = theParameters->GetProcPro
227 G4double QeExc = theParameters->GetProcPro
228 theParameters->GetProcPro
229 common.ProbProjectileDiffraction =
230 theParameters->GetProcProb( 2, ProjectileR
231 common.ProbTargetDiffraction =
232 theParameters->GetProcProb( 3, ProjectileR
233 common.ProbOfDiffraction = common.ProbProjec
234 #ifdef debugFTFexcitation
235 G4cout << "Proc Probs " << QeNoExc << " " <<
236 << common.ProbProjectileDiffraction <
237 << "ProjectileRapidity " << Projectil
238 #endif
239
240 if ( QeNoExc + QeExc + common.ProbProjectile
241 QeNoExc = 1.0 - QeExc - common.ProbProject
242 }
243 if ( QeExc + QeNoExc != 0.0 ) {
244 common.ProbExc = QeExc / ( QeExc + QeNoExc
245 }
246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
247 common.ProbProjectileDiffraction /= ( 1.0
248 common.ProbTargetDiffraction /= ( 1.0
249 }
250 #ifdef debugFTFexcitation
251 G4cout << "Proc Probs " << QeNoExc << " " <<
252 << common.ProbProjectileDiffraction <
253 << "ProjectileRapidity " << Projectil
254 #endif
255
256 // Try out charge exchange
257 G4int returnCode = 1;
258 if ( G4UniformRand() < QeExc + QeNoExc ) {
259 returnCode = ExciteParticipants_doChargeEx
260
261 }
262
263 G4bool returnResult = false;
264 if ( returnCode == 0 ) {
265 returnResult = true; // Successfully ende
266 } else if ( returnCode == 1 ) {
267
268 common.ProbOfDiffraction = common.ProbProj
269 #ifdef debugFTFexcitation
270 G4cout << "Excitation --------------------
271 << "Proj M0 MdMin MndMin " << commo
272 << common.ProjectileDiffStateMinMas
273 << G4endl
274 << "Targ M0 MdMin MndMin " << commo
275 << " " << common.TargetNonDiffState
276 << G4endl
277 << "Prob: ProjDiff TargDiff + Sum "
278 << common.ProbTargetDiffraction <<
279 #endif
280 if ( common.ProbOfDiffraction != 0.0 ) {
281 common.ProbProjectileDiffraction /= comm
282 } else {
283 common.ProbProjectileDiffraction = 0.0;
284 }
285 #ifdef debugFTFexcitation
286 G4cout << "Prob: ProjDiff TargDiff + Sum "
287 << common.ProbTargetDiffraction <<
288 #endif
289 common.ProjectileDiffStateMinMass2 = sq
290 common.ProjectileNonDiffStateMinMass2 = sq
291 common.TargetDiffStateMinMass2 = sq
292 common.TargetNonDiffStateMinMass2 = sq
293 // Choose between diffraction and non-diff
294 if ( G4UniformRand() < common.ProbOfDiffra
295
296 returnResult = ExciteParticipants_doDiff
297
298 } else {
299
300 returnResult = ExciteParticipants_doNonD
301
302 }
303 if ( returnResult ) {
304 common.Pprojectile += common.Qmomentum;
305 common.Ptarget -= common.Qmomentum;
306 // Transform back and update SplitableHa
307 common.Pprojectile.transform( common.toL
308 common.Ptarget.transform( common.toLab )
309 #ifdef debugFTFexcitation
310 G4cout << "Mproj " << common.Pprojectile
311 << G4endl;
312 #endif
313 projectile->Set4Momentum( common.Pprojec
314 target->Set4Momentum( common.Ptarget );
315 projectile->IncrementCollisionCount( 1 )
316 target->IncrementCollisionCount( 1 );
317 }
318 }
319
320 return returnResult;
321 }
322
323 //--------------------------------------------
324
325 G4int G4DiffractiveExcitation::
326 ExciteParticipants_doChargeExchange( G4VSplita
327 G4VSplita
328 G4FTFPara
329 G4Elastic
330 G4Diffrac
331 // First of the three utility methods used o
332 // it does the sampling for the charge-excha
333 // This method returns an integer code - ins
334 // "0" : successfully ended and nothing el
335 // "1" : successfully completed, but the w
336 // "99" : unsuccessfully ended, nothing els
337 G4int returnCode = 99;
338
339 G4double DeltaProbAtQuarkExchange = theParam
340 G4ParticleDefinition* TestParticle = 0;
341 G4double MtestPr = 0.0, MtestTr = 0.0;
342
343 #ifdef debugFTFexcitation
344 G4cout << "Q exchange ----------------------
345 #endif
346
347 // The target hadron is always a nucleon (i.
348 // never an antinucleon), therefore only a q
349 // exchanged between the projectile hadron a
350 // we could get a quark-quark-antiquark syst
351 // This implies that any projectile meson or
352 // a constituent quark in all cases - can ha
353 // hadron. Instead, any projectile anti-bary
354 // with a target hadron (because it has only
355 // projectile baryons, instead can have char
356
357 G4int NewProjCode = 0, NewTargCode = 0, Proj
358 // Projectile unpacking
359 if ( common.absProjectilePDGcode < 1000 ) {
360 UnpackMeson( common.ProjectilePDGcode, Pr
361 } else {
362 UnpackBaryon( common.ProjectilePDGcode, Pr
363 }
364 // Target unpacking
365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;
366 UnpackBaryon( common.TargetPDGcode, TargQ1,
367 #ifdef debugFTFexcitation
368 G4cout << "Proj Quarks " << ProjQ1 << " " <<
369 << "Targ Quarks " << TargQ1 << " " <<
370 #endif
371
372 // Sampling of exchanged quarks
373 G4int ProjExchangeQ = 0, TargExchangeQ = 0;
374 if ( common.absProjectilePDGcode < 1000 ) {
375
376 G4bool isProjQ1Quark = false;
377 ProjExchangeQ = ProjQ2;
378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark
379 isProjQ1Quark = true;
380 ProjExchangeQ = ProjQ1;
381 }
382 // Exchange of non-identical quarks is all
383 G4int NpossibleStates = 0;
384 if ( ProjExchangeQ != TargQ1 ) NpossibleSt
385 if ( ProjExchangeQ != TargQ2 ) NpossibleSt
386 if ( ProjExchangeQ != TargQ3 ) NpossibleSt
387 G4int Nsampled = (G4int)G4RandFlat::shootI
388 NpossibleStates = 0;
389 if ( ProjExchangeQ != TargQ1 ) {
390 if ( ++NpossibleStates == Nsampled ) {
391 TargExchangeQ = TargQ1; TargQ1 = ProjE
392 isProjQ1Quark ? ProjQ1 = TargExchangeQ
393 }
394 }
395 if ( ProjExchangeQ != TargQ2 ) {
396 if ( ++NpossibleStates == Nsampled ) {
397 TargExchangeQ = TargQ2; TargQ2 = ProjE
398 isProjQ1Quark ? ProjQ1 = TargExchangeQ
399 }
400 }
401 if ( ProjExchangeQ != TargQ3 ) {
402 if ( ++NpossibleStates == Nsampled ) {
403 TargExchangeQ = TargQ3; TargQ3 = ProjE
404 isProjQ1Quark ? ProjQ1 = TargExchangeQ
405 }
406 }
407 #ifdef debugFTFexcitation
408 G4cout << "Exchanged Qs in Pr Tr " << Proj
409 #endif
410
411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ
412 G4bool ProjExcited = false;
413 const G4int maxNumberOfAttempts = 50;
414 G4int attempts = 0;
415 while ( attempts++ < maxNumberOfAttempts )
416
417 // Determination of a new projectile ID
418 G4double ProbSpin0 = 0.5;
419 G4double Ksi = G4UniformRand();
420 if ( aProjQ1 == aProjQ2 ) {
421 if ( G4UniformRand() < ProbSpin0 ) {
422 if ( aProjQ1 < 3 ) {
423 NewProjCode = 111;
424 if ( Ksi < 0.5 ) {
425 NewProjCode = 221;
426 if ( Ksi < 0.25 ) {
427 NewProjCode = 331;
428 }
429 }
430 } else if ( aProjQ1 == 3 ) {
431 NewProjCode = 221;
432 if ( Ksi < 0.5 ) {
433 NewProjCode = 331;
434 }
435 } else if ( aProjQ1 == 4 ) {
436 NewProjCode = 441; // eta
437 } else if ( aProjQ1 == 5 ) {
438 NewProjCode = 551; // eta
439 }
440 } else {
441 if ( aProjQ1 < 3 ) {
442 NewProjCode = 113;
443 if ( Ksi < 0.5 ) {
444 NewProjCode = 223;
445 }
446 } else if ( aProjQ1 == 3 ) {
447 NewProjCode = 333;
448 } else if ( aProjQ1 == 4 ) {
449 NewProjCode = 443; // J/p
450 } else if ( aProjQ1 == 5 ) {
451 NewProjCode = 553; // Ups
452 }
453 }
454 } else {
455 if ( aProjQ1 > aProjQ2 ) {
456 NewProjCode = aProjQ1*100 + aProjQ2*
457 } else {
458 NewProjCode = aProjQ2*100 + aProjQ1*
459 }
460 }
461 #ifdef debugFTFexcitation
462 G4cout << "NewProjCode " << NewProjCode
463 #endif
464 // Decide (with 50% probability) whether
465 // but not in the case of charmed and bo
466 // there are no excited charmed and bott
467 ProjExcited = false;
468 if ( aProjQ1 <= 3 && aProjQ2 <= 3 &&
469 NewProjCode += 2; // Excited meson (J
470 ProjExcited = true;
471 }
472
473 // So far we have used the absolute valu
474 // now look at their signed values to se
475 G4int value = ProjQ1, absValue = aProjQ1
476 for ( G4int iQ = 0; iQ < 2; ++iQ ) { //
477 if ( iQ == 1 ) {
478 value = ProjQ2; absValue = aProjQ2;
479 }
480 if ( absValue == 2 || absValue == 4
481 Qquarks += 2*value/absValue; // u,
482 } else {
483 Qquarks -= value/absValue; // d,
484 }
485 }
486 // If Qquarks is positive, the sign of N
487 // If Qquarks is negative, then the sign
488 // If Qquarks is zero, then we need to d
489 // 1. If aProjQ1 and aProjQ2 are the sam
490 // (because the antiparticle is the s
491 // 2. If aProjQ1 and aProjQ2 are not the
492 // we have only two possibilities:
493 // a. aProjQ1 and aProjQ2 are two dif
494 // (s,d) or (b,d), or (b,s). In th
495 // is fine (because the heaviest o
496 // to be anti-quark belonging to t
497 // implies a meson with positive P
498 // b. aProjQ1 and aProjQ2 are two dif
499 // The heaviest of the two (c) has
500 // in the projectile, therefore th
501 // reverse: 421 -> -421 anti_D0 (c
502 if ( Qquarks < 0 || ( Qquarks == 0 &&
503 NewProjCode *= -1;
504 }
505 #ifdef debugFTFexcitation
506 G4cout << "NewProjCode +2 or 0 " << NewP
507 G4cout<<"+++++++++++++++++++++++++++++++
508 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquark
509 G4cout<<"+++++++++++++++++++++++++++++++
510 #endif
511
512 // Proj
513 TestParticle = G4ParticleTable::GetParti
514 if ( ! TestParticle ) continue;
515 common.MminProjectile = common.BrW.GetMi
516 if ( common.SqrtS - common.M0target < co
517 MtestPr = common.BrW.SampleMass( TestPar
518
519 #ifdef debugFTFexcitation
520 G4cout << "TestParticle Name " << NewPro
521 << G4endl
522 << "MtestPart MtestPart0 "<<Mtest
523 << "M0projectile projectile PDGMa
524 << projectile->GetDefinition()->G
525 #endif
526
527 // Targ
528 NewTargCode = NewNucleonId( TargQ1, Targ
529 #ifdef debugFTFexcitation
530 G4cout << "New TrQ " << TargQ1 << " " <<
531 << "NewTargCode " << NewTargCode
532 #endif
533
534 // If the target has not heavy (charm or
535 // see whether a Delta isobar can be cre
536 if ( TargQ1 <= 3 && TargQ2 <= 3 && T
537 if ( TargQ1 != TargQ2 && TargQ1 != T
538 if ( G4UniformRand() < 0.5 ) {
539 NewTargCode += 2;
540 } else if ( G4UniformRand() < 0.75 )
541 NewTargCode = 3122; // Lambda
542 }
543 } else if ( TargQ1 == TargQ2 && Targ
544 NewTargCode += 2; ProjExcited = true
545 } else if ( target->GetDefinition()->G
546 if ( G4UniformRand() > DeltaProbAtQu
547 NewTargCode += 2; ProjExcited = tr
548 }
549 } else if ( ! ProjExcited &&
550 G4UniformRand() < DeltaPro
551 common.SqrtS > common.M0pr
552 G4ParticleTable::GetP
553 NewTargCode += 2; // Create Delta i
554 }
555 }
556
557 // Protection against:
558 // - Lambda* (i.e. excited Lambda state)
559 // - Sigma* (i.e. excited Sigma hyperon
560 // - Xi* (i.e. excited Xi hyperon st
561 if ( NewTargCode == 3124 || // Lambda
562 NewTargCode == 3224 || // Sigma*+ NOT
563 NewTargCode == 3214 || // Sigma*0 NOT
564 NewTargCode == 3114 || // Sigma*- NOT
565 NewTargCode == 3324 || // Xi*0 NOT
566 NewTargCode == 3314 ) { // Xi*- NOT
567 //G4cout << "G4DiffractiveExcitation::Excite
568 // << NewTargCode << G4endl;
569 NewTargCode -= 2; // Corresponding ground-s
570 }
571
572 // Special treatment for charmed and bot
573 // so we need to transform them by hand
574 #ifdef debug_heavyHadrons
575 G4int initialNewTargCode = NewTargCode;
576 #endif
577 if ( NewTargCode == 4322 ) NewTargC
578 else if ( NewTargCode == 4312 ) NewTargC
579 else if ( NewTargCode == 5312 ) NewTargC
580 else if ( NewTargCode == 5322 ) NewTargC
581 #ifdef debug_heavyHadrons
582 if ( NewTargCode != initialNewTargCode )
583 G4cout << "G4DiffractiveExcitation::ExcitePa
584 << "\t target heavy baryon with pdgCo
585 << " into pdgCode=" << NewTargCode <<
586 }
587 #endif
588
589 TestParticle = G4ParticleTable::GetParti
590 if ( ! TestParticle ) continue;
591 #ifdef debugFTFexcitation
592 G4cout << "New targ " << NewTargCode <<
593 #endif
594 common.MminTarget = common.BrW.GetMinimu
595 if ( common.SqrtS - MtestPr < common.Mmi
596 MtestTr = common.BrW.SampleMass( TestPar
597
598 if ( common.SqrtS > MtestPr + MtestTr )
599
600 } // End of while loop
601 if ( attempts >= maxNumberOfAttempts ) ret
602
603 if ( MtestPr >= common.Pprojectile.mag()
604 common.M0projectile = MtestPr;
605 }
606 #ifdef debugFTFexcitation
607 G4cout << "M0projectile After check " << c
608 #endif
609 common.M0projectile2 = common.M0projectile
610 common.ProjectileDiffStateMinMass = com
611 common.ProjectileNonDiffStateMinMass = com
612 if ( MtestTr >= common.Ptarget.mag() ||
613 common.M0target = MtestTr;
614 }
615 common.M0target2 = common.M0target * commo
616 #ifdef debugFTFexcitation
617 G4cout << "New targ M0 M0^2 " << common.M0
618 #endif
619 common.TargetDiffStateMinMass = common.
620 common.TargetNonDiffStateMinMass = common.
621
622 } else { // of the if ( common.absProjectil
623
624 // If it is a projectile anti-baryon, no q
625 // therefore returns immediately 1 (which
626 // needs to be continued").
627 if ( common.ProjectilePDGcode < 0 ) return
628
629 // Choose randomly, with equal probability
630 // projectile or target hadron for selecti
631 G4bool isProjectileExchangedQ = false;
632 G4int firstQ = TargQ1, secondQ =
633 G4int otherFirstQ = ProjQ1, otherSecondQ =
634 if ( G4UniformRand() < 0.5 ) {
635 isProjectileExchangedQ = true;
636 firstQ = ProjQ1; secondQ = Pro
637 otherFirstQ = TargQ1; otherSecondQ = Tar
638 }
639 // Choose randomly, with equal probability
640 // selected (projectile or target) hadron
641 G4int exchangedQ = 0;
642 G4double Ksi = G4UniformRand();
643 if ( Ksi < 0.333333 ) {
644 exchangedQ = firstQ;
645 } else if ( 0.333333 <= Ksi && Ksi < 0.6
646 exchangedQ = secondQ;
647 } else {
648 exchangedQ = thirdQ;
649 }
650 #ifdef debugFTFexcitation
651 G4cout << "Exchange Qs isProjectile Q " <<
652 #endif
653 // The exchanged quarks (one of the projec
654 // are always accepted if they have differ
655 // same flavour) they are accepted only wi
656 G4double probSame = theParameters->GetProb
657 const G4int MaxCount = 100;
658 G4int count = 0, otherExchangedQ = 0;
659 do {
660 if ( exchangedQ != otherFirstQ || G4Un
661 otherExchangedQ = otherFirstQ; otherFi
662 } else {
663 if ( exchangedQ != otherSecondQ || G
664 otherExchangedQ = otherSecondQ; othe
665 } else {
666 if ( exchangedQ != otherThirdQ ||
667 otherExchangedQ = otherThirdQ; oth
668 }
669 }
670 }
671 } while ( otherExchangedQ == 0 && ++coun
672 if ( count >= MaxCount ) return returnCode
673 // Swap (between projectile and target had
674 // as "exchanged" quarks.
675 if ( Ksi < 0.333333 ) {
676 firstQ = exchangedQ;
677 } else if ( 0.333333 <= Ksi && Ksi < 0.6
678 secondQ = exchangedQ;
679 } else {
680 thirdQ = exchangedQ;
681 }
682 if ( isProjectileExchangedQ ) {
683 ProjQ1 = firstQ; ProjQ2 = secondQ;
684 TargQ1 = otherFirstQ; TargQ2 = otherSeco
685 } else {
686 TargQ1 = firstQ; TargQ2 = secondQ;
687 ProjQ1 = otherFirstQ; ProjQ2 = otherSeco
688 }
689 #ifdef debugFTFexcitation
690 G4cout << "Exchange Qs Pr Tr " << ( isPro
691 << " " << ( isProjectileExchangedQ
692 #endif
693
694 NewProjCode = NewNucleonId( ProjQ1, ProjQ2
695 NewTargCode = NewNucleonId( TargQ1, TargQ2
696 // Decide whether the new projectile hadro
697 // then decide whether the new target hadr
698 // Notice that a Delta particle has the la
699 // whereas a nucleon has "2" (because its
700 // If the new projectile hadron or the new
701 // constituent quark, then skip this part
702 // excited charm and bottom hadrons).
703 for ( G4int iHadron = 0; iHadron < 2; iHad
704 // First projectile hadron, then target
705 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2,
706 G4double massConstraint = common.M0targe
707 G4bool isHadronADelta = ( projectile->Ge
708 if ( iHadron == 1 ) { // Target hadron
709 codeQ1 = TargQ1, codeQ2 = TargQ2, code
710 massConstraint = common.M0projectile;
711 isHadronADelta = ( target->GetDefiniti
712 }
713 if ( codeQ1 > 3 || codeQ2 > 3 || cod
714 if ( codeQ1 == codeQ2 && codeQ1 == cod
715 newHadCode += 2; // Delta++ (uuu) or
716 } else if ( isHadronADelta ) { // Hadro
717 if ( G4UniformRand() > DeltaProbAtQuar
718 newHadCode += 2; // Delta+ (uud) or
719 } else {
720 newHadCode += 0; // No delta (so th
721 }
722 } else { // Hadron (projectile or targe
723 if ( G4UniformRand() < DeltaProbAtQuar
724 common.SqrtS > G4ParticleTable::G
725 + massConstraint )
726 newHadCode += 2; // Delta+ (uud) or
727 } else {
728 newHadCode += 0; // No delta (so th
729 }
730 }
731 if ( iHadron == 0 ) { // Projectile had
732 NewProjCode = newHadCode;
733 } else { // Target hadron
734 NewTargCode = newHadCode;
735 }
736 }
737 #ifdef debugFTFexcitation
738 G4cout << "NewProjCode NewTargCode " << Ne
739 #endif
740
741 // Protection against:
742 // - Lambda* (i.e. excited Lambda state)
743 // - Sigma* (i.e. excited Sigma hyperon s
744 // - Xi* (i.e. excited Xi hyperon stat
745 if ( NewProjCode == 3124 || // Lambda*
746 NewProjCode == 3224 || // Sigma*+ NOT ex
747 NewProjCode == 3214 || // Sigma*0 NOT ex
748 NewProjCode == 3114 || // Sigma*- NOT ex
749 NewProjCode == 3324 || // Xi*0 NOT ex
750 NewProjCode == 3314 ) { // Xi*- NOT ex
751 //G4cout << "G4DiffractiveExcitation::Ex
752 // << NewProjCode << G4endl;
753 NewProjCode -= 2; // Corresponding grou
754 }
755 if ( NewTargCode == 3124 || // Lambda*
756 NewTargCode == 3224 || // Sigma*+ NOT ex
757 NewTargCode == 3214 || // Sigma*0 NOT ex
758 NewTargCode == 3114 || // Sigma*- NOT ex
759 NewTargCode == 3324 || // Xi*0 NOT ex
760 NewTargCode == 3314 ) { // Xi*- NOT ex
761 //G4cout << "G4DiffractiveExcitation::Ex
762 // << NewTargCode << G4endl;
763 NewTargCode -= 2; // Corresponding grou
764 }
765
766 // Special treatment for charmed and botto
767 // so we need to transform them by hand to
768 #ifdef debug_heavyHadrons
769 G4int initialNewProjCode = NewProjCode, in
770 #endif
771 if ( NewProjCode == 4322 ) NewProjCod
772 else if ( NewProjCode == 4312 ) NewProjCod
773 else if ( NewProjCode == 5312 ) NewProjCod
774 else if ( NewProjCode == 5322 ) NewProjCod
775 if ( NewTargCode == 4322 ) NewTargCod
776 else if ( NewTargCode == 4312 ) NewTargCod
777 else if ( NewTargCode == 5312 ) NewTargCod
778 else if ( NewTargCode == 5322 ) NewTargCod
779 #ifdef debug_heavyHadrons
780 if ( NewProjCode != initialNewProjCode ||
781 G4cout << "G4DiffractiveExcitation::Exci
782 << "\t heavy baryon into an exist
783 if ( NewProjCode != initialNewProjCode )
784 G4cout << "\t \t projectile baryon with pdgC
785 << " into pdgCode=" << NewProjCode <<
786 }
787 if ( NewTargCode != initialNewTargCode )
788 G4cout << "\t \t target baryon with pd
789 << " into pdgCode=" << NewTargCode <<
790 }
791 }
792 #endif
793
794 // Sampling of the masses of the projectil
795 // Because of energy conservation, the ord
796 // randomly, half of the time we sample fi
797 // then the projectile nucleon mass, and t
798 // sample first the projectile nucleon mas
799 G4VSplitableHadron* firstHadron = target;
800 G4VSplitableHadron* secondHadron = project
801 G4int firstHadronCode = NewTargCode, secon
802 G4double massConstraint = common.M0project
803 G4bool isFirstTarget = true;
804 if ( G4UniformRand() < 0.5 ) {
805 // Sample first the projectile nucleon m
806 firstHadron = projectile; secondHa
807 firstHadronCode = NewProjCode; secondHa
808 massConstraint = common.M0target;
809 isFirstTarget = false;
810 }
811 G4double Mtest1st = 0.0, Mtest2nd = 0.0, M
812 for ( int iSamplingCase = 0; iSamplingCase
813 G4VSplitableHadron* aHadron = firstHadro
814 G4int aHadronCode = firstHadronCode;
815 if ( iSamplingCase == 1 ) { // Second n
816 aHadron = secondHadron; aHadronCode =
817 }
818 G4double MtestHadron = 0.0, MminHadron =
819 if ( aHadron->GetStatus() == 1 || aHad
820 TestParticle = G4ParticleTable::GetPar
821 if ( ! TestParticle ) return returnCod
822 MminHadron = common.BrW.GetMinimumMass
823 if ( common.SqrtS - massConstraint < M
824 if ( TestParticle->GetPDGWidth() == 0.
825 MtestHadron = common.BrW.SampleMass(
826 } else {
827 const G4int maxNumberOfAttempts = 50
828 G4int attempts = 0;
829 while ( attempts < maxNumberOfAttemp
830 attempts++;
831 MtestHadron = common.BrW.SampleMas
832
833 if ( common.SqrtS < MtestHadron +
834 continue; // Kinematically unac
835 } else {
836 break; // Kinematically acce
837 }
838 }
839 if ( attempts >= maxNumberOfAttempts
840 }
841 }
842 if ( iSamplingCase == 0 ) {
843 Mtest1st = MtestHadron; Mmin1st = Mmi
844 } else {
845 Mtest2nd = MtestHadron; Mmin2nd = Mmi
846 }
847 } // End for loop on the two sampling cas
848 if ( isFirstTarget ) {
849 MtestTr = Mtest1st; MtestPr = Mtest2n
850 common.MminTarget = Mmin1st; common.Mmi
851 } else {
852 MtestTr = Mtest2nd; MtestPr = Mtest1s
853 common.MminTarget = Mmin2nd; common.Mmi
854 }
855
856 if ( MtestPr != 0.0 ) {
857 common.M0projectile = MtestPr;
858 common.M0projectile2 = common.M0projecti
859 common.ProjectileDiffStateMinMass = c
860 common.ProjectileNonDiffStateMinMass = c
861 }
862 if ( MtestTr != 0.0 ) {
863 common.M0target = MtestTr;
864 common.M0target2 = common.M0target * com
865 common.TargetDiffStateMinMass = commo
866 common.TargetNonDiffStateMinMass = commo
867 }
868
869 } // End of if ( common.absProjectilePDGcod
870
871 // If we assume that final state hadrons aft
872 // in the ground states, we have to put
873 if ( common.SqrtS < common.M0projectile + co
874
875 common.PZcms2 = ( sqr( common.S ) + sqr( com
876 - 2.0 * ( common.S * ( com
877 + common.M0proje
878 #ifdef debugFTFexcitation
879 G4cout << "At the end// NewProjCode " << New
880 << "At the end// NewTargCode " << New
881 << "M0pr M0tr SqS " << common.M0pro
882 << common.SqrtS << G4endl
883 << "M0pr2 M0tr2 SqS " << common.M0pro
884 << common.SqrtS << G4endl
885 << "PZcms2 after the change " << comm
886 #endif
887 if ( common.PZcms2 < 0.0 ) return returnCode
888
889 projectile->SetDefinition( G4ParticleTable::
890 target->SetDefinition( G4ParticleTable::GetP
891 common.PZcms = std::sqrt( common.PZcms2 );
892 common.Pprojectile.setPz( common.PZcms );
893 common.Pprojectile.setE( std::sqrt( common.M
894 common.Ptarget.setPz( -common.PZcms );
895 common.Ptarget.setE( std::sqrt( common.M0tar
896 #ifdef debugFTFexcitation
897 G4cout << "Proj Targ and Proj+Targ in CMS" <
898 << common.Ptarget << G4endl << common
899 #endif
900
901 if ( projectile->GetStatus() != 0 ) projecti
902 if ( target->GetStatus() != 0 ) target->
903
904 // Check for possible excitation of the part
905 if ( common.SqrtS < common.M0projectile + co
906 common.SqrtS < common.ProjectileDiffSta
907 common.ProbOfDiffraction == 0.0 ) commo
908
909 if ( G4UniformRand() > common.ProbExc ) { /
910 #ifdef debugFTFexcitation
911 G4cout << "Make elastic scattering of new
912 #endif
913 common.Pprojectile.transform( common.toLab
914 common.Ptarget.transform( common.toLab );
915 projectile->Set4Momentum( common.Pprojecti
916 target->Set4Momentum( common.Ptarget );
917 G4bool Result = theElastic->ElasticScatter
918 #ifdef debugFTFexcitation
919 G4cout << "Result of el. scatt " << Result
920 << G4endl << projectile->Get4Moment
921 << projectile->Get4Momentum() + tar
922 << (projectile->Get4Momentum() + ta
923 #endif
924 if ( Result ) returnCode = 0; // successf
925 return returnCode;
926 }
927
928 #ifdef debugFTFexcitation
929 G4cout << "Make excitation of new hadrons" <
930 #endif
931
932 // Redefinition of ProbOfDiffraction because
933 common.ProbOfDiffraction = common.ProbProjec
934 if ( common.ProbOfDiffraction != 0.0 ) {
935 common.ProbProjectileDiffraction /= common
936 common.ProbTargetDiffraction /= common
937 }
938
939 return returnCode = 1; // successfully comp
940 }
941
942 //--------------------------------------------
943
944 G4bool G4DiffractiveExcitation::
945 ExciteParticipants_doDiffraction( G4VSplitable
946 G4FTFParamet
947 G4Diffractiv
948 // Second of the three utility methods used
949 // it does the sampling for the diffraction
950
951 G4bool isProjectileDiffraction = false;
952 if ( G4UniformRand() < common.ProbProjectile
953 isProjectileDiffraction = true;
954 #ifdef debugFTFexcitation
955 G4cout << "projectile diffraction" << G4en
956 #endif
957 common.ProjMassT2 = common.ProjectileDiffS
958 common.ProjMassT = common.ProjectileDiffS
959 common.TargMassT2 = common.M0target2;
960 common.TargMassT = common.M0target;
961 } else {
962 #ifdef debugFTFexcitation
963 G4cout << "Target diffraction" << G4endl;
964 #endif
965 common.ProjMassT2 = common.M0projectile2;
966 common.ProjMassT = common.M0projectile;
967 common.TargMassT2 = common.TargetDiffState
968 common.TargMassT = common.TargetDiffState
969 }
970
971 // Check whether the interaction is possible
972 if ( common.SqrtS < common.ProjMassT + commo
973
974 common.PZcms2 = ( sqr( common.S ) + sqr( com
975 - 2.0 * ( common.S * ( com
976 + common.ProjMassT2
977 if ( common.PZcms2 < 0.0 ) return false;
978 common.maxPtSquare = common.PZcms2;
979
980 G4double DiffrAveragePt2 = theParameters->Ge
981 G4bool loopCondition = true;
982 G4int whilecount = 0;
983 do { // Generate pt and mass of projectile
984
985 whilecount++;
986 if ( whilecount > 1000 ) {
987 common.Qmomentum = G4LorentzVector( 0.0,
988 return false; // Ignore this interacti
989 };
990
991 common.Qmomentum = G4LorentzVector( Gaussi
992 common.Pt2 = G4ThreeVector( common.Qmoment
993 if ( isProjectileDiffraction ) { // proje
994 common.ProjMassT2 = common.ProjectileDif
995 common.TargMassT2 = common.M0target2 + c
996 } else { // targe
997 common.ProjMassT2 = common.M0projectile2
998 common.TargMassT2 = common.TargetDiffSta
999 }
1000 common.ProjMassT = std::sqrt( common.Proj
1001 common.TargMassT = std::sqrt( common.Targ
1002 if ( common.SqrtS < common.ProjMassT + co
1003
1004 common.PZcms2 = ( sqr( common.S ) + sqr(
1005 - 2.0 * ( common.S * (
1006 + common.Proj
1007 if ( common.PZcms2 < 0.0 ) continue;
1008
1009 common.PZcms = std::sqrt( common.PZcms2 )
1010 if ( isProjectileDiffraction ) { // proj
1011 common.PMinusMin = std::sqrt( common.Pr
1012 common.PMinusMax = common.SqrtS - commo
1013 common.PMinusNew = ChooseP( common.PMin
1014 common.TMinusNew = common.SqrtS - commo
1015 common.Qminus = common.Ptarget.minus()
1016 common.TPlusNew = common.TargMassT2 / c
1017 common.Qplus = common.Ptarget.plus() -
1018 common.Qmomentum.setPz( (common.Qplus -
1019 common.Qmomentum.setE( (common.Qplus +
1020 loopCondition = ( ( common.Pprojectile
1021 common.ProjectileDiff
1022 } else { // targ
1023 common.TPlusMin = std::sqrt( common.Tar
1024 common.TPlusMax = common.SqrtS - common
1025 common.TPlusNew = ChooseP( common.TPlus
1026 common.PPlusNew = common.SqrtS - common
1027 common.Qplus = common.PPlusNew - common
1028 common.PMinusNew = common.ProjMassT2 /
1029 common.Qminus = common.PMinusNew - comm
1030 common.Qmomentum.setPz( (common.Qplus -
1031 common.Qmomentum.setE( (common.Qplus +
1032 loopCondition = ( ( common.Ptarget - c
1033 common.TargetDiffSta
1034 }
1035
1036 } while ( loopCondition ); /* Loop checkin
1037 // Repeat the sampling because ther
1038
1039 if ( isProjectileDiffraction ) { // projec
1040 projectile->SetStatus( 0 );
1041 if ( projectile->GetStatus() == 2 ) proje
1042 if ( target->GetStatus() == 1 && target
1043 } else { // target
1044 target->SetStatus( 0 );
1045 }
1046
1047 return true;
1048 }
1049
1050 //-------------------------------------------
1051
1052 G4bool G4DiffractiveExcitation::
1053 ExciteParticipants_doNonDiffraction( G4VSplit
1054 G4VSplit
1055 G4FTFPar
1056 G4Diffra
1057 // Third of the three utility methods used
1058 // it does the sampling for the non-diffrac
1059
1060 #ifdef debugFTFexcitation
1061 G4cout << "Non-diffraction process" << G4en
1062 #endif
1063
1064 // Check whether the interaction is possibl
1065 common.ProjMassT2 = common.ProjectileNonDif
1066 common.ProjMassT = common.ProjectileNonDif
1067 common.TargMassT2 = common.TargetNonDiffSta
1068 common.TargMassT = common.TargetNonDiffSta
1069 if ( common.SqrtS < common.ProjMassT + comm
1070
1071 common.PZcms2 = ( sqr( common.S ) + sqr( co
1072 - 2.0 * ( common.S * ( comm
1073 + common.ProjMassT2 *
1074 common.maxPtSquare = common.PZcms2;
1075
1076 G4int whilecount = 0;
1077 do { // Generate pt and masses
1078
1079 whilecount++;
1080 if ( whilecount > 1000 ) {
1081 common.Qmomentum = G4LorentzVector( 0.0
1082 return false; // Ignore this interacti
1083 };
1084
1085 common.Qmomentum = G4LorentzVector( Gauss
1086
1087 common.Pt2 = G4ThreeVector( common.Qmomen
1088 common.ProjMassT2 = common.ProjectileNonD
1089 common.ProjMassT = std::sqrt( common.Pro
1090 common.TargMassT2 = common.TargetNonDiffS
1091 common.TargMassT = std::sqrt( common.Tar
1092 if ( common.SqrtS < common.ProjMassT + co
1093
1094 common.PZcms2 =( sqr( common.S ) + sqr( c
1095 - 2.0 * ( common.S * ( c
1096 + common.ProjM
1097 if ( common.PZcms2 < 0.0 ) continue;
1098
1099 common.PZcms = std::sqrt( common.PZcms2 )
1100 common.PMinusMin = std::sqrt( common.Proj
1101 common.PMinusMax = common.SqrtS - common.
1102 common.TPlusMin = std::sqrt( common.TargM
1103 common.TPlusMax = common.SqrtS - common.P
1104
1105 if ( G4UniformRand() <= 0.5 ) { // Rando
1106 if ( G4UniformRand() < theParameters->G
1107 common.PMinusNew = ChooseP( common.PM
1108 } else {
1109 common.PMinusNew = ( common.PMinusMax
1110 }
1111 if ( G4UniformRand() < theParameters->G
1112 common.TPlusNew = ChooseP( common.TPl
1113 } else {
1114 common.TPlusNew = ( common.TPlusMax -
1115 }
1116 } else {
1117 if ( G4UniformRand() < theParameters->G
1118 common.TPlusNew = ChooseP( common.TPl
1119 } else {
1120 common.TPlusNew = ( common.TPlusMax -
1121 }
1122 if ( G4UniformRand() < theParameters->G
1123 common.PMinusNew = ChooseP( common.PM
1124 } else {
1125 common.PMinusNew = ( common.PMinusMax
1126 }
1127 }
1128
1129 common.Qminus = common.PMinusNew - common
1130 common.Qplus = -( common.TPlusNew - commo
1131 common.Qmomentum.setPz( (common.Qplus - c
1132 common.Qmomentum.setE( (common.Qplus + c
1133 #ifdef debugFTFexcitation
1134 G4cout <<"Sampled: Mpr, MdifPr, Mtr, Mdif
1135 << ( common.Pprojectile + common.Q
1136 << common.ProjectileNonDiffStateMi
1137 << ( common.Ptarget - common.Qmome
1138 << common.TargetNonDiffStateMinMas
1139 #endif
1140
1141 } while ( ( common.Pprojectile + common.Qmo
1142 common.ProjectileNonDiffStateMinM
1143 ( common.Ptarget - common.Qmoment
1144 common.TargetNonDiffStateMinMass2
1145
1146 projectile->SetStatus( 0 );
1147 target->SetStatus( 0 );
1148 return true;
1149 }
1150
1151
1152 //===========================================
1153
1154 void G4DiffractiveExcitation::CreateStrings(
1155
1156
1157
1158
1159
1160 //G4cout << "Create Strings SplitUp " << ha
1161 // << "Defin " << hadron->GetDefiniti
1162 // << "Defin " << hadron->GetDefiniti
1163
1164 G4bool HadronIsString = hadron->IsSplit();
1165 if( ! HadronIsString ) hadron->SplitUp();
1166
1167 G4Parton* start = hadron->GetNextParton();
1168 if ( start == nullptr ) {
1169 G4cout << " G4FTFModel::String() Error: N
1170 FirstString = 0; SecondString = 0;
1171 return;
1172 }
1173
1174 G4Parton* end = hadron->GetNextParton();
1175 if ( end == nullptr ) {
1176 G4cout << " G4FTFModel::String() Error: N
1177 FirstString = 0; SecondString = 0;
1178 return;
1179 }
1180
1181 //G4cout << start << " " << start->GetPDGco
1182 // << G4endl
1183 // << "Create string " << start->GetP
1184
1185 if ( HadronIsString ) {
1186 if ( isProjectile ) {
1187 FirstString = new G4ExcitedString( end,
1188 } else {
1189 FirstString = new G4ExcitedString( end,
1190 }
1191 FirstString->SetTimeOfCreation( hadron->G
1192 FirstString->SetPosition( hadron->GetPosi
1193 SecondString = 0;
1194 return;
1195 }
1196
1197 G4LorentzVector Phadron = hadron->Get4Momen
1198 //G4cout << "String mom " << Phadron << G4e
1199 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0
1200 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1201 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 )
1202 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0
1203 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0
1204
1205 G4int PDGcode_startQ = std::abs( start->Get
1206 G4int PDGcode_endQ = std::abs( end->Get
1207 //G4cout << "PDGcode_startQ " << PDGcode_st
1208
1209 G4double Wmin( 0.0 );
1210 if ( isProjectile ) {
1211 Wmin = theParameters->GetProjMinDiffMass(
1212 } else {
1213 Wmin = theParameters->GetTarMinDiffMass()
1214 }
1215
1216 G4double W = hadron->Get4Momentum().mag();
1217 G4double W2 = W*W;
1218 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 );
1219 G4bool Kink = false;
1220
1221 if ( ! ( ( start->GetDefinition()->GetParti
1222 end->GetDefinition()->GetParti
1223 ( start->GetDefinition()->GetParti
1224 end->GetDefinition()->GetParti
1225 // Kinky strings are allowed only for qq-
1226 // Kinky strings are impossible for other
1227 // according to the analysis of Pbar P in
1228
1229 if ( W > Wmin ) { // Kink is possible
1230 if ( hadron->GetStatus() == 0 ) {
1231 G4double Pt2kink = theParameters->Get
1232 if ( Pt2kink ) {
1233 Pt = std::sqrt( Pt2kink * ( G4Pow::
1234 } else {
1235 Pt = 0.0;
1236 }
1237 } else {
1238 Pt = 0.0;
1239 }
1240
1241 if ( Pt > 500.0*MeV ) {
1242 G4double Ymax = G4Log( W/2.0/Pt + std
1243 G4double Y = Ymax*( 1.0 - 2.0*G4Unifo
1244 x1 = 1.0 - Pt/W * G4Exp( Y );
1245 x3 = 1.0 - Pt/W * G4Exp(-Y );
1246 //x2 = 2.0 - x1 - x3;
1247
1248 G4double Mass_startQ = 650.0*MeV;
1249 if ( PDGcode_startQ < 3 ) Mass_start
1250 if ( PDGcode_startQ == 3 ) Mass_start
1251 if ( PDGcode_startQ == 4 ) Mass_start
1252 if ( PDGcode_startQ == 5 ) Mass_start
1253 G4double Mass_endQ = 650.0*MeV;
1254 if ( PDGcode_endQ < 3 ) Mass_endQ =
1255 if ( PDGcode_endQ == 3 ) Mass_endQ =
1256 if ( PDGcode_endQ == 4 ) Mass_endQ =
1257 if ( PDGcode_endQ == 5 ) Mass_endQ =
1258
1259 G4double P2_1 = W2*x1*x1/4.0 - Mass_e
1260 G4double P2_3 = W2*x3*x3/4.0 - Mass_s
1261 G4double P2_2 = sqr( (2.0 - x1 - x3)*
1262 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1263 Kink = false;
1264 } else {
1265 G4double P_1 = std::sqrt( P2_1 );
1266 G4double P_2 = std::sqrt( P2_2 );
1267 G4double P_3 = std::sqrt( P2_3 );
1268 G4double CosT12 = ( P2_3 - P2_1 - P
1269 G4double CosT13 = ( P2_2 - P2_1 - P
1270
1271 if ( std::abs( CosT12 ) > 1.0 ||
1272 Kink = false;
1273 } else {
1274 Kink = true;
1275 Pt = P_2 * std::sqrt( 1.0 - CosT1
1276 Pstart.setPx( -Pt ); Pstart.setPy
1277 Pend.setPx( 0.0 ); Pend.setPy(
1278 Pkink.setPx( Pt ); Pkink.setPy(
1279 Pstart.setE( x3*W/2.0 );
1280 Pkink.setE( Pkink.vect().mag() );
1281 Pend.setE( x1*W/2.0 );
1282
1283 G4double XkQ = GetQuarkFractionOf
1284 if ( Pkink.getZ() > 0.0 ) {
1285 if ( XkQ > 0.5 ) {
1286 PkinkQ1 = XkQ*Pkink;
1287 } else {
1288 PkinkQ1 = (1.0 - XkQ)*Pkink;
1289 }
1290 } else {
1291 if ( XkQ > 0.5 ) {
1292 PkinkQ1 = (1.0 - XkQ)*Pkink;
1293 } else {
1294 PkinkQ1 = XkQ*Pkink;
1295 }
1296 }
1297
1298 PkinkQ2 = Pkink - PkinkQ1;
1299 // Minimizing Pt1^2+Pt3^2
1300 G4double Cos2Psi = ( sqr(x1) - sq
1301 std::sqrt( sqr
1302 G4double Psi = std::acos( Cos2Psi
1303
1304 G4LorentzRotation Rotate;
1305 if ( isProjectile ) {
1306 Rotate.rotateY( Psi );
1307 } else {
1308 Rotate.rotateY( pi + Psi );
1309 }
1310 Rotate.rotateZ( twopi * G4Uniform
1311 Pstart *= Rotate;
1312 Pkink *= Rotate;
1313 PkinkQ1 *= Rotate;
1314 PkinkQ2 *= Rotate;
1315 Pend *= Rotate;
1316 }
1317 } // End of if ( P2_1 <= 0.0 || P2
1318 } // End of if ( Pt > 500.0*MeV )
1319 } // End of if ( W > Wmin ) : check for a
1320 } // end of qq-q string selection
1321
1322 if ( Kink ) { // Kink is possible
1323
1324 //G4cout << "Kink is sampled!" << G4endl;
1325 std::vector< G4double > QuarkProbabilitie
1326 theParameters->GetQuarkProbabilitiesA
1327
1328 G4int QuarkInGluon( 1 ); G4double Ksi = G
1329 for ( unsigned int Iq = 0; Iq < 3; Iq++ )
1330 //G4cout << "Iq " << Iq << G4endl;
1331 if ( Ksi > QuarkProbabilitiesAtGluonSpl
1332 }
1333 //G4cout << "Last Iq " << QuarkInGluon <<
1334 G4Parton* Gquark = new G4Parton( QuarkInG
1335 G4Parton* Ganti_quark = new G4Parton( -Qu
1336 //G4cout << "Lorentz " << G4endl;
1337
1338 G4LorentzRotation toCMS( -1 * Phadron.boo
1339 G4LorentzRotation toLab( toCMS.inverse()
1340 //G4cout << "Pstart " << Pstart << G4endl
1341 //G4cout << "Pend " << Pend << G4endl;
1342 //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1343 //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1344 //G4cout << "Pstart " << Pstart << G4endl
1345
1346 Pstart.transform( toLab ); start->Set4Mo
1347 PkinkQ1.transform( toLab );
1348 PkinkQ2.transform( toLab );
1349 Pend.transform( toLab ); end->Set4Mome
1350 //G4cout << "Pstart " << Pstart << G4endl
1351 //G4cout << "Pend " << Pend << G4endl;
1352 //G4cout << "Defin " << hadron->GetDefini
1353 //G4cout << "Defin " << hadron->GetDefini
1354
1355 //G4int absPDGcode = std::abs( hadron->Ge
1356 G4int absPDGcode = 1500;
1357 if ( start->GetDefinition()->GetParticleS
1358 end->GetDefinition()->GetParticleSub
1359 absPDGcode = 110;
1360 }
1361 //G4cout << "absPDGcode " << absPDGcode <
1362
1363 if ( absPDGcode < 1000 ) { // meson
1364 if ( isProjectile ) { // Projectile
1365 if ( end->GetDefinition()->GetPDGEnco
1366 FirstString = new G4ExcitedString(
1367 SecondString = new G4ExcitedString(
1368 Ganti_quark->Set4Momentum( PkinkQ1
1369 Gquark->Set4Momentum( PkinkQ2 );
1370 } else { // Anti_Quark on the end
1371 FirstString = new G4ExcitedString(
1372 SecondString = new G4ExcitedString(
1373 Gquark->Set4Momentum( PkinkQ1 );
1374 Ganti_quark->Set4Momentum( PkinkQ2
1375 }
1376 } else { // Target
1377 if ( end->GetDefinition()->GetPDGEnco
1378 FirstString = new G4ExcitedString(
1379 SecondString = new G4ExcitedString(
1380 Ganti_quark->Set4Momentum( PkinkQ2
1381 Gquark->Set4Momentum( PkinkQ1 );
1382 } else { // Anti_Quark on the end
1383 FirstString = new G4ExcitedString(
1384 SecondString = new G4ExcitedString(
1385 Gquark->Set4Momentum( PkinkQ2 );
1386 Ganti_quark->Set4Momentum( PkinkQ1
1387 }
1388 }
1389 } else { // Baryon/AntiBaryon
1390 if ( isProjectile ) { // Projectile
1391 if ( end->GetDefinition()->GetParticl
1392 end->GetDefinition()->GetPDGEnco
1393 FirstString = new G4ExcitedString(
1394 SecondString = new G4ExcitedString(
1395 Gquark->Set4Momentum( PkinkQ1 );
1396 Ganti_quark->Set4Momentum( PkinkQ2
1397 } else { /
1398 FirstString = new G4ExcitedString(
1399 SecondString = new G4ExcitedString(
1400 Ganti_quark->Set4Momentum( PkinkQ1
1401 Gquark->Set4Momentum( PkinkQ2 );
1402 }
1403 } else { // Target
1404 if ( end->GetDefinition()->GetParticl
1405 end->GetDefinition()->GetPDGEnco
1406 Gquark->Set4Momentum( PkinkQ1 );
1407 Ganti_quark->Set4Momentum( PkinkQ2
1408 FirstString = new G4ExcitedString(
1409 SecondString = new G4ExcitedString(
1410 } else { // Anti_DiQuark on the end
1411 FirstString = new G4ExcitedString(
1412 SecondString = new G4ExcitedString(
1413 Gquark->Set4Momentum( PkinkQ2 );
1414 Ganti_quark->Set4Momentum( PkinkQ1
1415 }
1416 }
1417 }
1418
1419 FirstString->SetTimeOfCreation( hadron->G
1420 FirstString->SetPosition( hadron->GetPosi
1421 SecondString->SetTimeOfCreation( hadron->
1422 SecondString->SetPosition( hadron->GetPos
1423
1424 } else { // Kink is impossible
1425
1426 if ( isProjectile ) {
1427 FirstString = new G4ExcitedString( end,
1428 } else {
1429 FirstString = new G4ExcitedString( end,
1430 }
1431 FirstString->SetTimeOfCreation( hadron->G
1432 FirstString->SetPosition( hadron->GetPosi
1433 SecondString = 0;
1434 // momenta of string ends
1435 G4LorentzVector HadronMom = hadron->Get4M
1436 G4LorentzVector Pstart1 = G4LorentzVector
1437 G4LorentzVector Pend1 = G4LorentzVector
1438 G4double Pz = HadronMom.pz();
1439 G4double Eh = HadronMom.e();
1440 G4double Pt2 = sqr( HadronMom.px() ) + sq
1441 G4double Mt2 = HadronMom.mt2();
1442 G4double Exp = std::sqrt( sqr(Pz) + ( sqr
1443 G4double Pzq = Pz/2.0 - Exp;
1444 G4double Eq = std::sqrt( sqr(Pzq) + Pt2
1445 G4double Pzqq = Pz/2.0 + Exp;
1446 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt
1447 start->Set4Momentum( Pstart1 );
1448 end->Set4Momentum( Pend1 );
1449 Pstart = Pstart1; Pend = Pend1;
1450
1451 } // End of "if (Kink)"
1452
1453 //G4cout << "Quarks in the string at creati
1454 // << " " << FirstString->GetLeftPart
1455 // << FirstString << " " << SecondStr
1456
1457 #ifdef G4_FTFDEBUG
1458 G4cout << " generated string flavors
1459 << end->GetPDGcode() << G4endl << "
1460 << start->Get4Momentum() << "mass :
1461 << " generated string momenta: Diqua
1462 << end->Get4Momentum().mag() << G4en
1463 << Pstart + Pend << G4endl << " Orig
1464 << hadron->Get4Momentum() << " "<<ha
1465 #endif
1466
1467 return;
1468 }
1469
1470
1471 //===========================================
1472
1473 G4double G4DiffractiveExcitation::ChooseP( G4
1474 // Choose an x between Xmin and Xmax with P
1475 // To be improved...
1476 G4double range = Pmax - Pmin;
1477 if ( Pmin <= 0.0 || range <= 0.0 ) {
1478 G4cout << " Pmin, range : " << Pmin << "
1479 throw G4HadronicException( __FILE__, __LI
1480 "G4Diffractive
1481 }
1482 G4double P = Pmin * G4Pow::GetInstance()->p
1483 //G4double P = (Pmax - Pmin) * G4UniformRan
1484 return P;
1485 }
1486
1487
1488 //===========================================
1489
1490 G4ThreeVector G4DiffractiveExcitation::Gaussi
1491 // @@ this method is used in FTFModel as w
1492 G4double Pt2( 0.0 );
1493 if ( AveragePt2 <= 0.0 ) {
1494 Pt2 = 0.0;
1495 } else {
1496 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor
1497 ( G4Ex
1498 }
1499 G4double Pt = std::sqrt( Pt2 );
1500 G4double phi = G4UniformRand() * twopi;
1501 return G4ThreeVector( Pt * std::cos( phi ),
1502 }
1503
1504
1505 //===========================================
1506
1507 G4double G4DiffractiveExcitation::GetQuarkFra
1508 G4double z, yf;
1509 const G4int maxNumberOfLoops = 10000;
1510 G4int loopCounter = 0;
1511 do {
1512 z = zmin + G4UniformRand() * (zmax - zmin
1513 yf = z*z + sqr(1.0 - z);
1514 } while ( ( G4UniformRand() > yf ) &&
1515 ++loopCounter < maxNumberOfLoops
1516 if ( loopCounter >= maxNumberOfLoops ) {
1517 z = 0.5*(zmin + zmax); // Just something
1518 }
1519 return z;
1520 }
1521
1522
1523 //===========================================
1524
1525 void G4DiffractiveExcitation::UnpackMeson( co
1526 G4int absIdPDG = std::abs( IdPDG );
1527 if ( ! ( absIdPDG == 111 || absIdPDG == 221
1528 absIdPDG == 441 || absIdPDG == 443
1529 // All other projectile mesons (including
1530 Q1 = absIdPDG / 100;
1531 Q2 = (absIdPDG % 100) / 10;
1532 G4int anti = 1 - 2 * ( std::max( Q1, Q2 )
1533 if ( IdPDG < 0 ) anti *= -1;
1534 Q1 *= anti;
1535 Q2 *= -1 * anti;
1536 } else {
1537 if ( absIdPDG == 441 || absIdPDG == 443 )
1538 Q1 = 4; Q2 = -4;
1539 } else if ( absIdPDG == 553 ) {
1540 Q1 = 5; Q2 = -5;
1541 } else {
1542 if ( G4UniformRand() < 0.5 ) {
1543 Q1 = 1; Q2 = -1;
1544 } else {
1545 Q1 = 2; Q2 = -2;
1546 }
1547 }
1548 }
1549 return;
1550 }
1551
1552
1553 //===========================================
1554
1555 void G4DiffractiveExcitation::UnpackBaryon( G
1556 G
1557 Q1 = IdPDG / 1000;
1558 Q2 = (IdPDG % 1000) / 100;
1559 Q3 = (IdPDG % 100) / 10;
1560 return;
1561 }
1562
1563
1564 //===========================================
1565
1566 G4int G4DiffractiveExcitation::NewNucleonId(
1567 // Order the three integers in such a way t
1568 G4int TmpQ( 0 );
1569 if ( Q3 > Q2 ) {
1570 TmpQ = Q2;
1571 Q2 = Q3;
1572 Q3 = TmpQ;
1573 } else if ( Q3 > Q1 ) {
1574 TmpQ = Q1;
1575 Q1 = Q3;
1576 Q3 = TmpQ;
1577 }
1578 if ( Q2 > Q1 ) {
1579 TmpQ = Q1;
1580 Q1 = Q2;
1581 Q2 = TmpQ;
1582 }
1583 // By now Q1 >= Q2 >= Q3
1584 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 +
1585 return NewCode;
1586 }
1587
1588
1589 //===========================================
1590
1591 G4DiffractiveExcitation::G4DiffractiveExcitat
1592 throw G4HadronicException( __FILE__, __LINE
1593 "G4DiffractiveEx
1594 }
1595
1596
1597 //===========================================
1598
1599 const G4DiffractiveExcitation & G4Diffractive
1600 throw G4HadronicException( __FILE__, __LINE
1601 "G4DiffractiveEx
1602 return *this;
1603 }
1604
1605
1606 //===========================================
1607
1608 G4bool G4DiffractiveExcitation::operator==( c
1609 throw G4HadronicException( __FILE__, __LINE
1610 "G4DiffractiveEx
1611 }
1612
1613
1614 //===========================================
1615
1616 G4bool G4DiffractiveExcitation::operator!= (
1617 throw G4HadronicException( __FILE__, __LINE
1618 "G4DiffractiveEx
1619 }
1620
1621