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