Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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