Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4FTFAnnihilation.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//
28
29// ------------------------------------------------------------
30// GEANT 4 class implemetation file
31//
32// ---------------- G4FTFAnnihilation --------------
33// by V. Uzhinsky, Spring 2011.
34// Take a projectile and a target
35// make annihilation or re-orangement of quarks and anti-quarks.
36// Ideas of Quark-Gluon-String model my A. Capella and A.B. Kaidalov
37// are implemented.
38// ---------------------------------------------------------------------
39
40#include "globals.hh"
41#include "Randomize.hh"
42#include "G4PhysicalConstants.hh"
43#include "G4SystemOfUnits.hh"
44
45#include "G4DiffractiveSplitableHadron.hh"
46#include "G4DiffractiveExcitation.hh"
47#include "G4FTFParameters.hh"
48#include "G4ElasticHNScattering.hh"
49#include "G4FTFAnnihilation.hh"
50
51#include "G4LorentzRotation.hh"
52#include "G4RotationMatrix.hh"
53#include "G4ThreeVector.hh"
54#include "G4ParticleDefinition.hh"
55#include "G4VSplitableHadron.hh"
56#include "G4ExcitedString.hh"
57#include "G4ParticleTable.hh"
58#include "G4Neutron.hh"
59#include "G4ParticleDefinition.hh"
60
61#include "G4Exp.hh"
62#include "G4Log.hh"
63#include "G4Pow.hh"
64
65//#include "UZHI_diffraction.hh"
66
67#include "G4ParticleTable.hh"
68
69//============================================================================
70
71//#define debugFTFannih
72
73
74//============================================================================
75
76G4FTFAnnihilation::G4FTFAnnihilation() {}
77
78
79//============================================================================
80
81G4FTFAnnihilation::~G4FTFAnnihilation() {}
82
83
84//============================================================================
85
86G4bool G4FTFAnnihilation::Annihilate( G4VSplitableHadron* projectile,
87 G4VSplitableHadron* target,
88 G4VSplitableHadron*& AdditionalString,
89 G4FTFParameters* theParameters ) const {
90
91 #ifdef debugFTFannih
92 G4cout << "---------------------------- Annihilation----------------" << G4endl;
93 #endif
94
95 CommonVariables common;
96
97 // Projectile parameters
98 common.Pprojectile = projectile->Get4Momentum();
99 G4int ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
100 if ( ProjectilePDGcode > 0 ) {
101 target->SetStatus( 3 );
102 return false;
103 }
104 G4double M0projectile2 = common.Pprojectile.mag2();
105
106 // Target parameters
107 G4int TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
108 common.Ptarget = target->Get4Momentum();
109 G4double M0target2 = common.Ptarget.mag2();
110
111 #ifdef debugFTFannih
112 G4cout << "PDG codes " << ProjectilePDGcode << " " << TargetPDGcode << G4endl
113 << "Pprojec " << common.Pprojectile << " " << common.Pprojectile.mag() << G4endl
114 << "Ptarget " << common.Ptarget << " " << common.Ptarget.mag() << G4endl
115 << "M0 proj target " << std::sqrt( M0projectile2 )
116 << " " << std::sqrt( M0target2 ) << G4endl;
117 #endif
118
119 // Kinematical properties of the interactions
120 G4LorentzVector Psum = common.Pprojectile + common.Ptarget; // 4-momentum in CMS
121 common.S = Psum.mag2();
122 common.SqrtS = std::sqrt( common.S );
123 #ifdef debugFTFannih
124 G4cout << "Psum SqrtS S " << Psum << " " << common.SqrtS << " " << common.S << G4endl;
125 #endif
126
127 // Transform momenta to cms and then rotate parallel to z axis
128 G4LorentzRotation toCms( -1*Psum.boostVector() );
129 G4LorentzVector Ptmp( toCms*common.Pprojectile );
130 toCms.rotateZ( -1*Ptmp.phi() );
131 toCms.rotateY( -1*Ptmp.theta() );
132 common.toLab = toCms.inverse();
133
134 if ( G4UniformRand() <= G4Pow::GetInstance()->powA( 1880.0/common.SqrtS, 4.0 ) ) {
135 common.RotateStrings = true;
136 common.RandomRotation.rotateZ( 2.0*pi*G4UniformRand() );
137 common.RandomRotation.rotateY( std::acos( 2.0*G4UniformRand() - 1.0 ) );
138 common.RandomRotation.rotateZ( 2.0*pi*G4UniformRand() );
139 }
140
141 G4double MesonProdThreshold = projectile->GetDefinition()->GetPDGMass() +
142 target->GetDefinition()->GetPDGMass() +
143 ( 2.0*140.0 + 16.0 )*MeV; // 2 Mpi + DeltaE
144 G4double Prel2 = sqr(common.S) + sqr(M0projectile2) + sqr(M0target2)
145 - 2.0*( common.S*(M0projectile2 + M0target2) + M0projectile2*M0target2 );
146 Prel2 /= common.S;
147 G4double X_a = 0.0, X_b = 0.0, X_c = 0.0, X_d = 0.0;
148 if ( Prel2 <= 0.0 ) {
149 // Annihilation at rest! Values are copied from Parameters
150 X_a = 625.1; // mb // 3-shirt diagram
151 X_b = 0.0; // mb // anti-quark-quark annihilation
152 X_c = 49.989; // mb // 2 Q-Qbar string creation
153 X_d = 6.614; // mb // One Q-Qbar string
154 #ifdef debugFTFannih
155 G4cout << "Annih at Rest X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d
156 << G4endl;
157 #endif
158 } else { // Annihilation in flight!
159 G4double FlowF = 1.0 / std::sqrt( Prel2 )*GeV;
160 // Process cross sections
161 X_a = 25.0*FlowF; // mb 3-shirt diagram
162 if ( common.SqrtS < MesonProdThreshold ) {
163 X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( ( MesonProdThreshold - common.SqrtS )/GeV, 2.5 );
164 } else {
165 X_b = 6.8*GeV / common.SqrtS; // mb anti-quark-quark annihilation
166 }
167 if ( projectile->GetDefinition()->GetPDGMass() + target->GetDefinition()->GetPDGMass()
168 > common.SqrtS ) {
169 X_b = 0.0;
170 }
171 // This can be in an interaction of low energy anti-baryon with off-shell nuclear nucleon
172 X_c = 2.0 * FlowF * sqr( projectile->GetDefinition()->GetPDGMass() +
173 target->GetDefinition()->GetPDGMass() ) / common.S; // mb re-arrangement of
174 // 2 quarks and 2 anti-quarks
175 X_d = 23.3*GeV*GeV / common.S; // mb anti-quark-quark string creation
176 #ifdef debugFTFannih
177 G4cout << "Annih in Flight X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d
178 << G4endl << "SqrtS MesonProdThreshold " << common.SqrtS << " " << MesonProdThreshold
179 << G4endl;
180 #endif
181 }
182
183 G4bool isUnknown = false;
184 if ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) { // Target proton or Delta+
185 if ( ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214 ) { // anti_proton or anti_Delta+
186 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // Pbar P
187 } else if ( ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114 ) { // anti_neutron or anti_Delta0
188 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // NeutrBar P
189 } else if ( ProjectilePDGcode == -3122 ) { // anti_Lambda (no anti_Lambda* in PDG)
190 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar P
191 } else if ( ProjectilePDGcode == -3112 ) { // anti_Sigma- (no anti_Sigma*- in G4)
192 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma-Bar P
193 } else if ( ProjectilePDGcode == -3212 ) { // anti_Sigma0 (no anti_Sigma*0 in G4)
194 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar P
195 } else if ( ProjectilePDGcode == -3222 ) { // anti_Sigma+ (no anti_Sigma*+ in G4)
196 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma+Bar P
197 } else if ( ProjectilePDGcode == -3312 ) { // anti_Xi- (no anti_Xi*- in G4)
198 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi-Bar P
199 } else if ( ProjectilePDGcode == -3322 ) { // anti_Xi0 (no anti_Xi*0 in G4)
200 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi0Bar P
201 } else if ( ProjectilePDGcode == -3334 ) { // anti_Omega- (no anti_Omega*- in PDG)
202 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar P
203 } else {
204 isUnknown = true;
205 }
206 } else if ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) { // Target neutron or Delta0
207 if ( ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214 ) { // anti_proton or anti_Delta+
208 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // Pbar N
209 } else if ( ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114 ) { // anti_neutron or anti_Delta0
210 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // NeutrBar N
211 } else if ( ProjectilePDGcode == -3122 ) { // anti_Lambda (no anti_Lambda* in PDG)
212 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar N
213 } else if ( ProjectilePDGcode == -3112 ) { // anti_Sigma- (no anti_Sigma*- in G4)
214 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma-Bar N
215 } else if ( ProjectilePDGcode == -3212 ) { // anti_Sigma0 (no anti_Sigma*0 in G4)
216 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar N
217 } else if ( ProjectilePDGcode == -3222 ) { // anti_Sigma+ (no anti_Sigma*+ in G4)
218 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma+Bar N
219 } else if ( ProjectilePDGcode == -3312 ) { // anti_Xi- (no anti_Xi*- in G4)
220 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi-Bar N
221 } else if ( ProjectilePDGcode == -3322 ) { // anti_Xi0 (no anti_Xi*0 in G4)
222 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi0Bar N
223 } else if ( ProjectilePDGcode == -3334 ) { // anti_Omega- (no anti_Omega*- in PDG)
224 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar N
225 } else {
226 isUnknown = true;
227 }
228 } else {
229 isUnknown = true;
230 }
231 if ( isUnknown ) {
232 G4cout << "Unknown anti-baryon for FTF annihilation: PDGcodes - "
233 << ProjectilePDGcode << " " << TargetPDGcode << G4endl;
234 }
235
236 #ifdef debugFTFannih
237 G4cout << "Annih Actual X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl;
238 #endif
239
240 G4double Xannihilation = X_a + X_b + X_c + X_d;
241
242 // Projectile unpacking
243 UnpackBaryon( ProjectilePDGcode, common.AQ[0], common.AQ[1], common.AQ[2] );
244
245 // Target unpacking
246 UnpackBaryon( TargetPDGcode, common.Q[0], common.Q[1], common.Q[2] );
247
248 G4double Ksi = G4UniformRand();
249
250 if ( Ksi < X_a / Xannihilation ) {
251 return Create3QuarkAntiQuarkStrings( projectile, target, AdditionalString, theParameters, common );
252 }
253
254 G4int resultCode = 99;
255 if ( Ksi < (X_a + X_b) / Xannihilation ) {
256 resultCode = Create1DiquarkAntiDiquarkString( projectile, target, common );
257 if ( resultCode == 0 ) {
258 return true;
259 } else if ( resultCode == 99 ) {
260 return false;
261 }
262 }
263
264 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation ) {
265 resultCode = Create2QuarkAntiQuarkStrings( projectile, target, theParameters, common );
266 if ( resultCode == 0 ) {
267 return true;
268 } else if ( resultCode == 99 ) {
269 return false;
270 }
271 }
272
273 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation ) {
274 return Create1QuarkAntiQuarkString( projectile, target, theParameters, common );
275 }
276
277 return true;
278}
279
280
281//-----------------------------------------------------------------------
282
283G4bool G4FTFAnnihilation::
284Create3QuarkAntiQuarkStrings( G4VSplitableHadron* projectile,
285 G4VSplitableHadron* target,
286 G4VSplitableHadron*& AdditionalString,
287 G4FTFParameters* theParameters,
288 G4FTFAnnihilation::CommonVariables& common ) const {
289 // Simulation of 3 anti-quark - quark strings creation
290
291 #ifdef debugFTFannih
292 G4cout << "Process a, 3 shirt diagram" << G4endl;
293 #endif
294
295 // Sampling kinematical properties of quark. It can be done before string's creation
296
297 const G4int maxNumberOfLoops = 1000;
298 G4double MassQ2 = 0.0; // Simplest case is considered with Mass_Q = 0.0
299 // In principle, this must work with Mass_Q != 0.0
300 G4double Quark_Xs[6];
301 G4ThreeVector Quark_Mom[6];
302
303 G4double Alfa_R = 0.5;
304 G4double AveragePt2 = 200.0*200.0, maxPtSquare = common.S;
305 G4double ScaleFactor = 1.0;
306 G4double Alfa = 0.0, Beta = 0.0;
307
308 G4int NumberOfTries = 0, loopCounter = 0;
309
310 do {
311 // Sampling X's of anti-baryon and baryon
312 G4double x1 = 0.0, x2 = 0.0, x3 = 0.0;
313 G4double Product = 1.0;
314 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // anti-baryon (1st case), baryon (2nd case)
315
316 G4double r1 = G4UniformRand(), r2 = G4UniformRand();
317 if ( Alfa_R == 1.0 ) {
318 x1 = 1.0 - std::sqrt( r1 );
319 x2 = (1.0 - x1) * r2;
320 } else {
321 x1 = sqr( r1 );
322 x2 = (1.0 - x1) * sqr( std::sin( pi/2.0*r2 ) );
323 }
324 x3 = 1.0 - x1 - x2;
325
326 G4int index = iCase*3; // 0 for anti-baryon, 3 for baryon
327 Quark_Xs[index] = x1; Quark_Xs[index+1] = x2; Quark_Xs[index+2] = x3;
328 Product *= (x1*x2*x3);
329 }
330
331 if ( Product == 0.0 ) continue;
332
333 ++NumberOfTries;
334 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
335 // After a large number of tries, it is better to reduce the values of <Pt^2>
336 ScaleFactor /= 2.0;
337 AveragePt2 *= ScaleFactor;
338 }
339
340 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
341 for ( G4int i = 0; i < 6; ++i ) {
342 Quark_Mom [i] = GaussianPt( AveragePt2, maxPtSquare );
343 PtSum += Quark_Mom[i];
344 }
345
346 PtSum /= 6.0;
347 Alfa = 0.0; Beta = 0.0;
348
349 for ( G4int i = 0; i < 6; ++i ) { // Loop over the quarks and (anti-)quarks
350 Quark_Mom[i] -= PtSum;
351
352 G4double val = ( Quark_Mom[i].mag2() + MassQ2 ) / Quark_Xs[i];
353 if ( i < 3 ) { // anti-baryon
354 Alfa += val;
355 } else { // baryon (iCase == 1)
356 Beta += val;
357 }
358 }
359
360 } while ( ( std::sqrt( Alfa ) + std::sqrt( Beta ) > common.SqrtS ) &&
361 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
362
363 if ( loopCounter >= maxNumberOfLoops ) {
364 return false;
365 }
366
367 G4double DecayMomentum2 = sqr(common.S) + sqr(Alfa) + sqr(Beta)
368 - 2.0*( common.S*(Alfa + Beta) + Alfa*Beta );
369
370 G4double WminusTarget = 0.0, WplusProjectile = 0.0;
371 WminusTarget = ( common.S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.SqrtS;
372 WplusProjectile = common.SqrtS - Beta/WminusTarget;
373
374 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // anti-baryon (1st case), baryon (2nd case)
375 G4int index = iCase*3; // 0 for anti-baryon, 3 for baryon
376 G4double w = WplusProjectile; // for anti-baryon
377 if ( iCase == 1 ) w = - WminusTarget; // for baryon
378 for ( G4int i = 0; i < 3; ++i ) {
379 G4double Pz = w * Quark_Xs[index+i] / 2.0 -
380 ( Quark_Mom[index+i].mag2() + MassQ2 ) /
381 ( 2.0 * w * Quark_Xs[index+i] );
382 Quark_Mom[index+i].setZ( Pz );
383 }
384 }
385
386 // Sampling of anti-quark order in projectile
387 G4int SampledCase = (G4int)G4RandFlat::shootInt( 6 );
388 G4int Tmp1 = 0, Tmp2 = 0;
389 switch ( SampledCase ) {
390 case 1 : Tmp1 = common.AQ[1]; common.AQ[1] = common.AQ[2]; common.AQ[2] = Tmp1; break;
391 case 2 : Tmp1 = common.AQ[0]; common.AQ[0] = common.AQ[1]; common.AQ[1] = Tmp1; break;
392 case 3 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = common.AQ[2];
393 common.AQ[1] = Tmp1; common.AQ[2] = Tmp2; break;
394 case 4 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = Tmp2;
395 common.AQ[1] = common.AQ[2]; common.AQ[2] = Tmp1; break;
396 case 5 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = common.AQ[2];
397 common.AQ[1] = Tmp2; common.AQ[2] = Tmp1; break;
398 }
399
400 // Set the string properties
401 // An anti quark - quark pair can have the quantum number of either a scalar meson
402 // or a vector meson: the last digit of the PDG code is, respectively, 1 and 3.
403 // For simplicity only scalar is considered here.
404 G4int NewCode = 0, antiQuark = 0, quark = 0;
405 G4ParticleDefinition* TestParticle = nullptr;
406 for ( G4int iString = 0; iString < 3; ++iString ) { // Loop over the 3 string cases
407 if ( iString == 0 ) {
408 antiQuark = common.AQ[0]; quark = common.Q[0];
409 projectile->SetFirstParton( antiQuark );
410 projectile->SetSecondParton( quark );
411 projectile->SetStatus( 0 );
412 } else if ( iString == 1 ) {
413 quark = common.Q[1]; antiQuark = common.AQ[1];
414 target->SetFirstParton( quark );
415 target->SetSecondParton( antiQuark );
416 target->SetStatus( 0 );
417 } else { // iString == 2
418 antiQuark = common.AQ[2]; quark = common.Q[2];
419 }
420 G4int absAntiQuark = std::abs( antiQuark ), absQuark = std::abs( quark );
421 G4double aKsi = G4UniformRand();
422 if ( absAntiQuark == absQuark ) {
423 if ( absAntiQuark != 3 ) { // Not yet considered the case absAntiQuark 4 (charm) and 5 (bottom)
424 NewCode = 111; // Pi0-meson
425 if ( aKsi < 0.5 ) {
426 NewCode = 221; // Eta -meson
427 if ( aKsi < 0.25 ) {
428 NewCode = 331; // Eta'-meson
429 }
430 }
431 } else {
432 NewCode = 221; // Eta -meson
433 if ( aKsi < 0.5 ) {
434 NewCode = 331; // Eta'-meson
435 }
436 }
437 } else { // Vector mesons - rho, omega, phi (not yet considered the analogous cases for charm and bottom)
438 if ( absAntiQuark > absQuark ) {
439 NewCode = absAntiQuark*100 + absQuark*10 + 1; NewCode *= absAntiQuark/antiQuark;
440 } else {
441 NewCode = absQuark*100 + absAntiQuark*10 + 1; NewCode *= absQuark/quark;
442 }
443 }
444 if ( iString == 2 ) AdditionalString = new G4DiffractiveSplitableHadron();
445 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
446 if ( ! TestParticle ) return false;
447 if ( iString == 0 ) {
448 projectile->SetDefinition( TestParticle );
449 theParameters->SetProjMinDiffMass( 0.5 ); // 0.5 GeV : Min diffractive mass of pi-meson
450 theParameters->SetProjMinNonDiffMass( 0.5 ); // It must be self-consistent with Parameters
451 } else if ( iString == 1 ) {
452 target->SetDefinition( TestParticle );
453 theParameters->SetTarMinDiffMass( 0.5 );
454 theParameters->SetTarMinNonDiffMass( 0.5 );
455 } else { // iString == 2
456 AdditionalString->SetDefinition( TestParticle );
457 AdditionalString->SetFirstParton( common.AQ[2] );
458 AdditionalString->SetSecondParton( common.Q[2] );
459 AdditionalString->SetStatus( 0 );
460 }
461 } // End of the for loop over the 3 string cases
462
463 // 1st string AQ[0]-Q[0], 2nd string AQ[1]-Q[1], 3rd string AQ[2]-Q[2]
464
465 G4LorentzVector Pstring1, Pstring2, Pstring3;
466 G4int QuarkOrder[3] = { 0 };
467 G4double YstringMax = 0.0, YstringMin = 0.0;
468 for ( G4int i = 0; i < 3; ++i ) {
469 G4ThreeVector tmp = Quark_Mom[i] + Quark_Mom[i+3];
470 G4LorentzVector Pstring( tmp, std::sqrt( Quark_Mom[i].mag2() + MassQ2 ) +
471 std::sqrt( Quark_Mom[i+3].mag2() + MassQ2 ) );
472 // Add protection for rapidity = 0.5*ln( (E+Pz)/(E-Pz) )
473 G4double Ystring = 0.0;
474 if ( Pstring.e() > 1.0e-30 ) {
475 if ( Pstring.e() + Pstring.pz() < 1.0e-30 ) { // Very small numerator in the logarithm
476 Ystring = -1.0e30; // A very large negative value (E ~= -Pz)
477 if ( Pstring.e() - Pstring.pz() < 1.0e-30 ) { // Very small denominator in the logarithm
478 Ystring = 1.0e30; // A very large positive value (E ~= Pz)
479 } else { // Normal case
480 Ystring = Pstring.rapidity();
481 }
482 }
483 }
484 // Keep ordering in rapidity: "1" highest, "2" middle, "3" smallest
485 if ( i == 0 ) {
486 Pstring1 = Pstring; YstringMax = Ystring;
487 QuarkOrder[0] = 0;
488 } else if ( i == 1 ) {
489 if ( Ystring > YstringMax ) {
490 Pstring2 = Pstring1; YstringMin = YstringMax;
491 Pstring1 = Pstring; YstringMax = Ystring;
492 QuarkOrder[0] = 1; QuarkOrder[1] = 0;
493 } else {
494 Pstring2 = Pstring; YstringMin = Ystring;
495 QuarkOrder[1] = 1;
496 }
497 } else { // i == 2
498 if ( Ystring > YstringMax ) {
499 Pstring3 = Pstring2;
500 Pstring2 = Pstring1;
501 Pstring1 = Pstring;
502 QuarkOrder[1] = QuarkOrder[0];
503 QuarkOrder[2] = QuarkOrder[1];
504 QuarkOrder[0] = 2;
505 } else if ( Ystring > YstringMin ) {
506 Pstring3 = Pstring2;
507 Pstring2 = Pstring;
508 } else {
509 Pstring3 = Pstring;
510 QuarkOrder[2] = 2;
511 }
512 }
513 }
514
515 G4LorentzVector Quark_4Mom[6];
516 for ( G4int i = 0; i < 6; ++i ) {
517 Quark_4Mom[i] = G4LorentzVector( Quark_Mom[i], std::sqrt( Quark_Mom[i].mag2() + MassQ2 ) );
518 if ( common.RotateStrings ) Quark_4Mom[i] *= common.RandomRotation;
519 Quark_4Mom[i].transform( common.toLab );
520 }
521
522 projectile->Splitting();
523 projectile->GetNextAntiParton()->Set4Momentum( Quark_4Mom[QuarkOrder[0]] );
524 projectile->GetNextParton()->Set4Momentum( Quark_4Mom[QuarkOrder[0]+3] );
525
526 target->Splitting();
527 target->GetNextParton()->Set4Momentum( Quark_4Mom[QuarkOrder[2]] );
528 target->GetNextAntiParton()->Set4Momentum( Quark_4Mom[QuarkOrder[2]+3] );
529
530 AdditionalString->Splitting();
531 AdditionalString->GetNextAntiParton()->Set4Momentum( Quark_4Mom[QuarkOrder[1]] );
532 AdditionalString->GetNextParton()->Set4Momentum( Quark_4Mom[QuarkOrder[1]+3] );
533
534 common.Pprojectile = Pstring1; // Highest rapidity
535 common.Ptarget = Pstring3; // Lowest rapidity
536 G4LorentzVector LeftString( Pstring2 ); // Middle rapidity
537
538 if ( common.RotateStrings ) {
539 common.Pprojectile *= common.RandomRotation;
540 common.Ptarget *= common.RandomRotation;
541 LeftString *= common.RandomRotation;
542 }
543
544 common.Pprojectile.transform( common.toLab );
545 common.Ptarget.transform( common.toLab );
546 LeftString.transform( common.toLab );
547
548 // Calculation of the creation time
549 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
550 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
551 projectile->SetPosition( target->GetPosition() );
552 AdditionalString->SetTimeOfCreation( target->GetTimeOfCreation() );
553 AdditionalString->SetPosition( target->GetPosition() );
554
555 projectile->Set4Momentum( common.Pprojectile );
556 AdditionalString->Set4Momentum( LeftString );
557 target->Set4Momentum( common.Ptarget );
558
559 projectile->IncrementCollisionCount( 1 );
560 AdditionalString->IncrementCollisionCount( 1 );
561 target->IncrementCollisionCount( 1 );
562
563 return true;
564}
565
566
567//-----------------------------------------------------------------------
568
569G4int G4FTFAnnihilation::
570Create1DiquarkAntiDiquarkString( G4VSplitableHadron* projectile,
571 G4VSplitableHadron* target,
572 G4FTFAnnihilation::CommonVariables& common ) const {
573 // Simulation of anti-diquark-diquark string creation.
574 // This method returns an integer code - instead of a boolean, with the following meaning:
575 // "0" : successfully ended and nothing else needs to be done;
576 // "1" : successfully completed, but the work needs to be continued;
577 // "99" : unsuccessfully ended, nothing else can be done.
578
579 #ifdef debugFTFannih
580 G4cout << "Process b, quark - anti-quark annihilation, di-q - anti-di-q string" << G4endl;
581 #endif
582
583 G4int CandidatsN = 0, CandAQ[9][2] = {}, CandQ[9][2] = {};
584 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // index of the 3 constituent anti-quarks of the antibaryon projectile
585 for ( G4int iQ = 0; iQ < 3; ++iQ ) { // index of the 3 constituent quarks of the target nucleon
586 if ( -common.AQ[iAQ] == common.Q[iQ] ) { // antiquark - quark that can annihilate
587 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
588 // of the (anti-baryon) projectile or (nucleon) target.
589 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
590 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
591 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
592 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
593 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
594 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
595 ++CandidatsN;
596 }
597 }
598 }
599
600 // Remaining two (anti-)quarks that form the (anti-)diquark
601 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, LeftQ2 = 0;
602 if ( CandidatsN != 0 ) {
603 G4int SampledCase = (G4int)G4RandFlat::shootInt( CandidatsN );
604 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0] ];
605 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1] ];
606 LeftQ1 = common.Q[ CandQ[SampledCase][0] ];
607 LeftQ2 = common.Q[ CandQ[SampledCase][1] ];
608
609 // Build anti-diquark and diquark : the last digit can be either 3 - for all combinations
610 // of anti-quark - anti-quark and quark - quark - or 1 - only when the two anti-quarks
611 // or quarks are different. For simplicity, only 3 is considered.
612 G4int Anti_DQ = 0, DQ = 0;
613 if ( std::abs( LeftAQ1 ) > std::abs( LeftAQ2 ) ) {
614 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3;
615 } else {
616 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3;
617 }
618 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 ) ) {
619 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3;
620 } else {
621 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3;
622 }
623
624 // Set the string properties
625 projectile->SetFirstParton( DQ );
626 projectile->SetSecondParton( Anti_DQ );
627
628 // It is assumed that quark and di-quark masses are 0.
629 G4LorentzVector Pquark = G4LorentzVector( 0.0, 0.0, -common.SqrtS/2.0, common.SqrtS/2.0 );
630 G4LorentzVector Paquark = G4LorentzVector( 0.0, 0.0, common.SqrtS/2.0, common.SqrtS/2.0 );
631
632 if ( common.RotateStrings ) {
633 Pquark *= common.RandomRotation;
634 Paquark *= common.RandomRotation;
635 }
636
637 Pquark.transform( common.toLab );
638 Paquark.transform( common.toLab );
639
640 projectile->GetNextParton()->Set4Momentum( Pquark );
641 projectile->GetNextAntiParton()->Set4Momentum( Paquark );
642
643 projectile->Splitting();
644
645 projectile->SetStatus( 0 );
646 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4
647 common.Pprojectile.setPx( 0.0 );
648 common.Pprojectile.setPy( 0.0 );
649 common.Pprojectile.setPz( 0.0 );
650 common.Pprojectile.setE( common.SqrtS );
651 common.Pprojectile.transform( common.toLab );
652
653 // Calculation of the creation time
654 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
655 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
656 projectile->SetPosition( target->GetPosition() );
657 projectile->Set4Momentum( common.Pprojectile );
658
659 projectile->IncrementCollisionCount( 1 );
660 target->IncrementCollisionCount( 1 );
661
662 return 0; // Completed successfully: nothing else to be done
663 } // End of if ( CandidatsN != 0 )
664
665 // If we allow the string to interact with other nuclear nucleons, we have to
666 // set up MinDiffrMass in Parameters, and ascribe a PDGEncoding. To be done yet!
667
668 return 1; // Successfully ended, but the work is not over
669}
670
671
672//-----------------------------------------------------------------------
673
674G4int G4FTFAnnihilation::
675Create2QuarkAntiQuarkStrings( G4VSplitableHadron* projectile,
676 G4VSplitableHadron* target,
677 G4FTFParameters* theParameters,
678 G4FTFAnnihilation::CommonVariables& common ) const {
679 // Simulation of 2 anti-quark-quark strings creation.
680 // This method returns an integer code - instead of a boolean, with the following meaning:
681 // "0" : successfully ended and nothing else needs to be done;
682 // "1" : successfully completed, but the work needs to be continued;
683 // "99" : unsuccessfully ended, nothing else can be done.
684
685 #ifdef debugFTFannih
686 G4cout << "Process c, quark - anti-quark and string junctions annihilation, 2 strings left."
687 << G4endl;
688 #endif
689
690 // Sampling kinematical properties: 1st string LeftAQ1-LeftQ1, 2nd string LeftAQ2-LeftQ2
691 G4ThreeVector Quark_Mom[4];
692 G4double Quark_Xs[4];
693 G4double AveragePt2 = 200.0*200.0, maxPtSquare = common.S, MassQ2 = 0.0, ScaleFactor = 1.0;
694 G4int NumberOfTries = 0, loopCounter = 0;
695 const G4int maxNumberOfLoops = 1000;
696 G4double Alfa = 0.0, Beta = 0.0;
697 G4double WminusTarget = 0.0, WplusProjectile = 0.0, Alfa_R = 0.5;
698 do {
699 // Sampling X's of the 2 quarks and 2 anti-quarks
700
701 G4double Product = 1.0;
702 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // Loop over the two strings
703 G4double x = 0.0, r = G4UniformRand();
704 if ( Alfa_R == 1.0 ) {
705 if ( iCase == 0 ) { // first string
706 x = std::sqrt( r );
707 } else { // second string
708 x = 1.0 - std::sqrt( r );
709 }
710 } else {
711 x = sqr( std::sin( pi/2.0*r ) );
712 }
713 G4int index = iCase*2; // 0 for the first string, 2 for the second string
714 Quark_Xs[index] = x ; Quark_Xs[index+1] = 1.0 - x ;
715 Product *= x*(1.0-x);
716 }
717
718 if ( Product == 0.0 ) continue;
719
720 ++NumberOfTries;
721 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
722 // After a large number of tries, it is better to reduce the values of <Pt^2>
723 ScaleFactor /= 2.0;
724 AveragePt2 *= ScaleFactor;
725 }
726
727 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
728 for( G4int i = 0; i < 4; ++i ) {
729 Quark_Mom[i] = GaussianPt( AveragePt2, maxPtSquare );
730 PtSum += Quark_Mom[i];
731 }
732
733 PtSum /= 4.0;
734 for ( G4int i = 0; i < 4; ++i ) {
735 Quark_Mom[i] -= PtSum;
736 }
737
738 Alfa = 0.0; Beta = 0.0;
739 for ( G4int iCase = 0; iCase < 2; ++iCase ) {
740 G4int index = iCase * 2;
741 for ( G4int i = 0; i < 2; ++i ) {
742 G4double val = ( Quark_Mom[index+i].mag2() + MassQ2 ) / Quark_Xs[index+i];
743 if ( iCase == 0 ) { // first string
744 Alfa += val;
745 } else { // second string
746 Beta += val;
747 }
748 }
749 }
750
751 } while ( ( std::sqrt( Alfa ) + std::sqrt( Beta ) > common.SqrtS ) &&
752 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
753
754 if ( loopCounter >= maxNumberOfLoops ) {
755 return 99; // unsuccessfully ended, nothing else can be done
756 }
757
758 G4double DecayMomentum2 = sqr(common.S) + sqr(Alfa) + sqr(Beta)
759 - 2.0*( common.S*(Alfa + Beta) + Alfa*Beta );
760 WminusTarget = ( common.S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.SqrtS;
761 WplusProjectile = common.SqrtS - Beta/WminusTarget;
762
763 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // Loop over the two strings
764 G4int index = iCase*2; // 0 for the first string, 2 for the second string
765 for ( G4int i = 0; i < 2; ++i ) {
766 G4double w = WplusProjectile; // For the first string
767 if ( iCase == 1 ) w = - WminusTarget; // For the second string
768 G4double Pz = w * Quark_Xs[index+i] / 2.0
769 - ( Quark_Mom[index+i].mag2() + MassQ2 ) /
770 ( 2.0 * w * Quark_Xs[index+i] );
771 Quark_Mom[index+i].setZ( Pz );
772 }
773 }
774
775 G4int CandidatsN = 0, CandAQ[9][2] = {}, CandQ[9][2] = {};
776 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, LeftQ2 = 0;
777 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // index of the 3 constituent anti-quarks of the antibaryon projectile
778 for ( G4int iQ = 0; iQ < 3; ++iQ ) { // index of the 3 constituent quarks of the nucleon target
779 if ( -common.AQ[iAQ] == common.Q[iQ] ) { // antiquark - quark that can annihilate
780 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
781 // of the (anti-baryon) projectile or (nucleon) target.
782 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
783 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
784 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
785 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
786 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
787 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
788 ++CandidatsN;
789 }
790 }
791 }
792
793 if ( CandidatsN != 0 ) {
794 G4int SampledCase = (G4int)G4RandFlat::shootInt( CandidatsN );
795 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0] ];
796 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1] ];
797 if ( G4UniformRand() < 0.5 ) {
798 LeftQ1 = common.Q[ CandQ[SampledCase][0] ];
799 LeftQ2 = common.Q[ CandQ[SampledCase][1] ];
800 } else {
801 LeftQ2 = common.Q[ CandQ[SampledCase][0] ];
802 LeftQ1 = common.Q[ CandQ[SampledCase][1] ];
803 }
804
805 // Set the string properties
806 // An anti quark - quark pair can have the quantum number of either a scalar meson
807 // or a vector meson: the last digit of the PDG code is, respectively, 1 and 3.
808 // For simplicity only scalar is considered here.
809 G4int NewCode = 0, antiQuark = 0, quark = 0;
810 G4ParticleDefinition* TestParticle = nullptr;
811 for ( G4int iString = 0; iString < 2; ++iString ) { // Loop over the 2 string cases
812 if ( iString == 0 ) {
813 antiQuark = LeftAQ1; quark = LeftQ1;
814 projectile->SetFirstParton( antiQuark );
815 projectile->SetSecondParton( quark );
816 projectile->SetStatus( 0 );
817 } else { // iString == 1
818 quark = LeftQ2; antiQuark = LeftAQ2;
819 target->SetFirstParton( quark );
820 target->SetSecondParton( antiQuark );
821 target->SetStatus( 0 );
822 }
823 G4int absAntiQuark = std::abs( antiQuark ), absQuark = std::abs( quark );
824 G4double aKsi = G4UniformRand();
825 if ( absAntiQuark == absQuark ) {
826 if ( absAntiQuark != 3 ) {
827 NewCode = 111; // Pi0-meson
828 if ( aKsi < 0.5 ) {
829 NewCode = 221; // Eta -meson
830 if ( aKsi < 0.25 ) {
831 NewCode = 331; // Eta'-meson
832 }
833 }
834 } else {
835 NewCode = 221; // Eta -meson
836 if ( aKsi < 0.5 ) {
837 NewCode = 331; // Eta'-meson
838 }
839 }
840 } else {
841 if ( absAntiQuark > absQuark ) {
842 NewCode = absAntiQuark*100 + absQuark*10 + 1; NewCode *= absAntiQuark/antiQuark;
843 } else {
844 NewCode = absQuark*100 + absAntiQuark*10 + 1; NewCode *= absQuark/quark;
845 }
846 }
847 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
848 if ( ! TestParticle ) return 99; // unsuccessfully ended, nothing else can be done
849 if ( iString == 0 ) {
850 projectile->SetDefinition( TestParticle );
851 theParameters->SetProjMinDiffMass( 0.5 );
852 theParameters->SetProjMinNonDiffMass( 0.5 );
853 } else { // iString == 1
854 target->SetDefinition( TestParticle );
855 theParameters->SetTarMinDiffMass( 0.5 );
856 theParameters->SetTarMinNonDiffMass( 0.5 );
857 }
858 } // End of loop over the 2 string cases
859
860 G4int QuarkOrder[2];
861 G4LorentzVector Pstring1, Pstring2;
862 G4double Ystring1 = 0.0, Ystring2 = 0.0;
863
864 for ( G4int iCase = 0; iCase < 2; ++iCase ) { // Loop over the two strings
865 G4ThreeVector tmp = Quark_Mom[iCase] + Quark_Mom[iCase+2];
866 G4LorentzVector Pstring( tmp, std::sqrt( Quark_Mom[iCase].mag2() + MassQ2 ) +
867 std::sqrt( Quark_Mom[iCase+2].mag2() + MassQ2 ) );
868 // Add protection for rapidity = 0.5*ln( (E+Pz)/(E-Pz) )
869 G4double Ystring = 0.0;
870 if ( Pstring.e() > 1.0e-30 ) {
871 if ( Pstring.e() + Pstring.pz() < 1.0e-30 ) { // Very small numerator in the logarithm
872 Ystring = -1.0e30; // A very large negative value (E ~= -Pz)
873 if ( Pstring.e() - Pstring.pz() < 1.0e-30 ) { // Very small denominator in the logarithm
874 Ystring = 1.0e30; // A very large positive value (E ~= Pz)
875 } else { // Normal case
876 Ystring = Pstring.rapidity();
877 }
878 }
879 }
880 if ( iCase == 0 ) { // For the first string
881 Pstring1 = Pstring; Ystring1 = Ystring;
882 } else { // For the second string
883 Pstring2 = Pstring; Ystring2 = Ystring;
884 }
885 }
886 if ( Ystring1 > Ystring2 ) {
887 common.Pprojectile = Pstring1; common.Ptarget = Pstring2;
888 QuarkOrder[0] = 0; QuarkOrder[1] = 1;
889 } else {
890 common.Pprojectile = Pstring2; common.Ptarget = Pstring1;
891 QuarkOrder[0] = 1; QuarkOrder[1] = 0;
892 }
893
894 if ( common.RotateStrings ) {
895 common.Pprojectile *= common.RandomRotation;
896 common.Ptarget *= common.RandomRotation;
897 }
898
899 common.Pprojectile.transform( common.toLab );
900 common.Ptarget.transform( common.toLab );
901
902 G4LorentzVector Quark_4Mom[4];
903 for ( G4int i = 0; i < 4; ++i ) {
904 Quark_4Mom[i] = G4LorentzVector( Quark_Mom[i], std::sqrt( Quark_Mom[i].mag2() + MassQ2 ) );
905 if ( common.RotateStrings ) Quark_4Mom[i] *= common.RandomRotation;
906 Quark_4Mom[i].transform( common.toLab );
907 }
908
909 projectile->Splitting();
910 projectile->GetNextAntiParton()->Set4Momentum( Quark_4Mom[QuarkOrder[0]] );
911 projectile->GetNextParton()->Set4Momentum( Quark_4Mom[QuarkOrder[0]+2] );
912
913 target->Splitting();
914 target->GetNextParton()->Set4Momentum( Quark_4Mom[QuarkOrder[1]] );
915 target->GetNextAntiParton()->Set4Momentum( Quark_4Mom[QuarkOrder[1]+2] );
916
917 // Calculation of the creation time
918 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
919 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
920 projectile->SetPosition( target->GetPosition() );
921 projectile->Set4Momentum( common.Pprojectile );
922 target->Set4Momentum( common.Ptarget );
923
924 projectile->IncrementCollisionCount( 1 );
925 target->IncrementCollisionCount( 1 );
926
927 return 0; // Completed successfully: nothing else to be done
928 } // End of if ( CandidatsN != 0 )
929
930 return 1; // Successfully ended, but the work is not over
931}
932
933
934//-----------------------------------------------------------------------
935
936G4bool G4FTFAnnihilation::
937Create1QuarkAntiQuarkString( G4VSplitableHadron* projectile,
938 G4VSplitableHadron* target,
939 G4FTFParameters* theParameters,
940 G4FTFAnnihilation::CommonVariables& common ) const {
941 // Simulation of anti-quark - quark string creation
942
943 #ifdef debugFTFannih
944 G4cout << "Process d, only 1 quark - anti-quark string" << G4endl;
945 #endif
946
947 // Determine the set of candidates anti-quark - quark pairs that do not annihilate.
948 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
949 // of the (anti-baryon) projectile or (nucleon) target.
950 G4int CandidatsN = 0, CandAQ[36], CandQ[36];
951 G4int LeftAQ = 0, LeftQ = 0;
952 for ( G4int iAQ1 = 0; iAQ1 < 3; ++iAQ1 ) {
953 for ( G4int iAQ2 = 0; iAQ2 < 3; ++iAQ2 ) {
954 if ( iAQ1 != iAQ2 ) {
955 for ( G4int iQ1 = 0; iQ1 < 3; ++iQ1 ) {
956 for ( G4int iQ2 = 0; iQ2 < 3; ++iQ2 ) {
957 if ( iQ1 != iQ2 ) {
958 if ( -common.AQ[iAQ1] == common.Q[iQ1] && -common.AQ[iAQ2] == common.Q[iQ2] ) {
959 if ( ( iAQ1 == 0 && iAQ2 == 1 ) || ( iAQ1 == 1 && iAQ2 == 0 ) ) {
960 CandAQ[CandidatsN] = 2;
961 } else if ( ( iAQ1 == 0 && iAQ2 == 2 ) || ( iAQ1 == 2 && iAQ2 == 0 ) ) {
962 CandAQ[CandidatsN] = 1;
963 } else if ( ( iAQ1 == 1 && iAQ2 == 2 ) || ( iAQ1 == 2 && iAQ2 == 1 ) ) {
964 CandAQ[CandidatsN] = 0;
965 }
966 if ( ( iQ1 == 0 && iQ2 == 1 ) || ( iQ1 == 1 && iQ2 == 0 ) ) {
967 CandQ[CandidatsN] = 2;
968 } else if ( ( iQ1 == 0 && iQ2 == 2 ) || ( iQ1 == 2 && iQ2 == 0 ) ) {
969 CandQ[CandidatsN] = 1;
970 } else if ( ( iQ1 == 1 && iQ2 == 2 ) || ( iQ1 == 2 && iQ2 == 1 ) ) {
971 CandQ[CandidatsN] = 0;
972 }
973 ++CandidatsN;
974 }
975 }
976 }
977 }
978 }
979 }
980 }
981
982 if ( CandidatsN != 0 ) {
983 G4int SampledCase = (G4int)G4RandFlat::shootInt( CandidatsN );
984 LeftAQ = common.AQ[ CandAQ[SampledCase] ];
985 LeftQ = common.Q[ CandQ[SampledCase] ];
986
987 // Set the string properties
988 projectile->SetFirstParton( LeftQ );
989 projectile->SetSecondParton( LeftAQ );
990 projectile->SetStatus( 0 );
991 G4int aAQ = std::abs( LeftAQ ), aQ = std::abs( LeftQ );
992 G4int NewCode = 0;
993 G4double aKsi = G4UniformRand();
994 // The string can have the quantum number of either a scalar or a vector (whose last digit
995 // of the PDG code is, respectively, 1 and 3). For simplicity only scalar is considered here.
996 if ( aAQ == aQ ) {
997 if ( aAQ != 3 ) {
998 NewCode = 111; // Pi0-meson
999 if ( aKsi < 0.5 ) {
1000 NewCode = 221; // Eta -meson
1001 if ( aKsi < 0.25 ) {
1002 NewCode = 331; // Eta'-meson
1003 }
1004 }
1005 } else {
1006 NewCode = 221; // Eta -meson
1007 if ( aKsi < 0.5 ) {
1008 NewCode = 331; // Eta'-meson
1009 }
1010 }
1011 } else {
1012 if ( aAQ > aQ ) {
1013 NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ;
1014 } else {
1015 NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ;
1016 }
1017 }
1018
1019 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
1020 if ( ! TestParticle ) return false;
1021 projectile->SetDefinition( TestParticle );
1022 theParameters->SetProjMinDiffMass( 0.5 );
1023 theParameters->SetProjMinNonDiffMass( 0.5 );
1024
1025 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4
1026 common.Pprojectile.setPx( 0.0 );
1027 common.Pprojectile.setPy( 0.0 );
1028 common.Pprojectile.setPz( 0.0 );
1029 common.Pprojectile.setE( common.SqrtS );
1030
1031 common.Pprojectile.transform( common.toLab );
1032
1033 G4LorentzVector Pquark = G4LorentzVector( 0.0, 0.0, -common.SqrtS/2.0, common.SqrtS/2.0 );
1034 G4LorentzVector Paquark = G4LorentzVector( 0.0, 0.0, +common.SqrtS/2.0, common.SqrtS/2.0 );
1035
1036 if ( common.RotateStrings ) {
1037 Pquark *= common.RandomRotation; Paquark *= common.RandomRotation;
1038 }
1039 Pquark.transform(common.toLab); projectile->GetNextParton()->Set4Momentum(Pquark);
1040 Paquark.transform(common.toLab); projectile->GetNextAntiParton()->Set4Momentum(Paquark);
1041
1042 projectile->Splitting();
1043
1044 // Calculation of the creation time
1045 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
1046 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1047 projectile->SetPosition( target->GetPosition() );
1048 projectile->Set4Momentum( common.Pprojectile );
1049
1050 projectile->IncrementCollisionCount( 1 );
1051 target->IncrementCollisionCount( 1 );
1052
1053 return true;
1054 } // End of if ( CandidatsN != 0 )
1055
1056 return true;
1057}
1058
1059
1060//============================================================================
1061
1062G4double G4FTFAnnihilation::ChooseX( G4double /* Alpha */, G4double /* Beta */ ) const {
1063 // If for sampling Xs other values of Alfa and Beta instead of 0.5 will be
1064 // chosen the method will be implemented
1065 //G4double tmp = Alpha*Beta;
1066 //tmp *= 1.0;
1067 return 0.5;
1068}
1069
1070
1071//============================================================================
1072
1073G4ThreeVector G4FTFAnnihilation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1074 // @@ this method is used in FTFModel as well. Should go somewhere common!
1075 G4double Pt2 = 0.0;
1076 if ( AveragePt2 <= 0.0 ) {
1077 Pt2 = 0.0;
1078 } else {
1079 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1080 ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
1081 }
1082 G4double Pt = std::sqrt( Pt2 );
1083 G4double phi = G4UniformRand() * twopi;
1084 return G4ThreeVector ( Pt*std::cos( phi ), Pt*std::sin( phi ), 0.0 );
1085}
1086
1087
1088//============================================================================
1089
1090void G4FTFAnnihilation::UnpackBaryon( G4int IdPDG, G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1091 G4int AbsId = std::abs( IdPDG );
1092 Q1 = AbsId / 1000;
1093 Q2 = ( AbsId % 1000 ) / 100;
1094 Q3 = ( AbsId % 100 ) / 10;
1095 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = -Q3; } // Anti-baryon
1096 return;
1097}
1098
1099
1100//============================================================================
1101
1102G4FTFAnnihilation::G4FTFAnnihilation( const G4FTFAnnihilation& ) {
1103 throw G4HadronicException( __FILE__, __LINE__,
1104 "G4FTFAnnihilation copy constructor not meant to be called" );
1105}
1106
1107
1108//============================================================================
1109
1110const G4FTFAnnihilation & G4FTFAnnihilation::operator=( const G4FTFAnnihilation& ) {
1111 throw G4HadronicException( __FILE__, __LINE__,
1112 "G4FTFAnnihilation = operator not meant to be called" );
1113}
1114
1115
1116//============================================================================
1117
1118G4bool G4FTFAnnihilation::operator==( const G4FTFAnnihilation& ) const {
1119 throw G4HadronicException( __FILE__, __LINE__,
1120 "G4FTFAnnihilation == operator not meant to be called" );
1121}
1122
1123
1124//============================================================================
1125
1126G4bool G4FTFAnnihilation::operator!=( const G4FTFAnnihilation& ) const {
1127 throw G4HadronicException( __FILE__, __LINE__,
1128 "G4DiffractiveExcitation != operator not meant to be called" );
1129}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:33
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:158
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:170
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:173
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4FTFAnnihilation::~G4FTFAnnihilation
virtual ~G4FTFAnnihilation()
Definition: G4FTFAnnihilation.cc:81
G4FTFAnnihilation::Annihilate
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
Definition: G4FTFAnnihilation.cc:86
G4FTFParameters::SetTarMinNonDiffMass
void SetTarMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:319
G4FTFParameters::SetProjMinDiffMass
void SetProjMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:307
G4FTFParameters::SetTarMinDiffMass
void SetTarMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:315
G4FTFParameters::SetProjMinNonDiffMass
void SetProjMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:311
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:586
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4Parton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
G4VSplitableHadron::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition: G4VSplitableHadron.hh:145
G4VSplitableHadron::SetStatus
void SetStatus(const G4int aStatus)
Definition: G4VSplitableHadron.hh:165
G4VSplitableHadron::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4VSplitableHadron.hh:135
G4VSplitableHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4VSplitableHadron.hh:120
G4VSplitableHadron::GetNextParton
virtual G4Parton * GetNextParton()=0
G4VSplitableHadron::SetSecondParton
virtual void SetSecondParton(G4int PDGcode)=0
G4VSplitableHadron::GetNextAntiParton
virtual G4Parton * GetNextAntiParton()=0
G4VSplitableHadron::SetDefinition
void SetDefinition(const G4ParticleDefinition *aDefinition)
Definition: G4VSplitableHadron.hh:130
G4VSplitableHadron::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4VSplitableHadron.hh:125
G4VSplitableHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4VSplitableHadron.hh:160
G4VSplitableHadron::IncrementCollisionCount
void IncrementCollisionCount(G4int aCount)
Definition: G4VSplitableHadron.hh:140
G4VSplitableHadron::SetFirstParton
virtual void SetFirstParton(G4int PDGcode)=0
G4VSplitableHadron::SetPosition
void SetPosition(const G4ThreeVector &aPosition)
Definition: G4VSplitableHadron.hh:155
sqr
T sqr(const T &x)
Definition: templates.hh:128