| 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() ); //AR-Jun2018
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 ) {
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 ) {
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 || ProjectilePDGcode == -3124 ) {
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 || ProjectilePDGcode == -3114 ) {
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 || ProjectilePDGcode == -3214 ) {
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 || ProjectilePDGcode == -3224 ) {
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 || ProjectilePDGcode == -3314 ) {
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 || ProjectilePDGcode == -3324 ) {
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 ) {
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 ) {
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 ) {
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 || ProjectilePDGcode == -3124 ) {
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 || ProjectilePDGcode == -3114 ) {
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 || ProjectilePDGcode == -3214 ) {
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 || ProjectilePDGcode == -3224 ) {
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 || ProjectilePDGcode == -3314 ) {
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 || ProjectilePDGcode == -3324 ) {
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 ) {
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 <<
236 #ifdef debugFTFannih 235 #ifdef debugFTFannih
237 G4cout << "Annih Actual X a b c d " << X_a < 236 G4cout << "Annih Actual X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl;
238 #endif 237 #endif
239 238
240 G4double Xannihilation = X_a + X_b + X_c + X 239 G4double Xannihilation = X_a + X_b + X_c + X_d;
241 240
242 // Projectile unpacking 241 // Projectile unpacking
243 UnpackBaryon( ProjectilePDGcode, common.AQ[0 242 UnpackBaryon( ProjectilePDGcode, common.AQ[0], common.AQ[1], common.AQ[2] );
244 243
245 // Target unpacking 244 // Target unpacking
246 UnpackBaryon( TargetPDGcode, common.Q[0], co 245 UnpackBaryon( TargetPDGcode, common.Q[0], common.Q[1], common.Q[2] );
247 246
248 G4double Ksi = G4UniformRand(); 247 G4double Ksi = G4UniformRand();
249 248
250 if ( Ksi < X_a / Xannihilation ) { 249 if ( Ksi < X_a / Xannihilation ) {
251 return Create3QuarkAntiQuarkStrings( proje 250 return Create3QuarkAntiQuarkStrings( projectile, target, AdditionalString, theParameters, common );
252 } 251 }
253 252
254 G4int resultCode = 99; 253 G4int resultCode = 99;
255 if ( Ksi < (X_a + X_b) / Xannihilation ) { 254 if ( Ksi < (X_a + X_b) / Xannihilation ) {
256 resultCode = Create1DiquarkAntiDiquarkStri 255 resultCode = Create1DiquarkAntiDiquarkString( projectile, target, common );
257 if ( resultCode == 0 ) { 256 if ( resultCode == 0 ) {
258 return true; 257 return true;
259 } else if ( resultCode == 99 ) { 258 } else if ( resultCode == 99 ) {
260 return false; 259 return false;
261 } 260 }
262 } 261 }
263 262
264 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilati 263 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation ) {
265 resultCode = Create2QuarkAntiQuarkStrings( 264 resultCode = Create2QuarkAntiQuarkStrings( projectile, target, theParameters, common );
266 if ( resultCode == 0 ) { 265 if ( resultCode == 0 ) {
267 return true; 266 return true;
268 } else if ( resultCode == 99 ) { 267 } else if ( resultCode == 99 ) {
269 return false; 268 return false;
270 } 269 }
271 } 270 }
272 271
273 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xanni 272 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation ) {
274 return Create1QuarkAntiQuarkString( projec 273 return Create1QuarkAntiQuarkString( projectile, target, theParameters, common );
275 } 274 }
276 275
277 return true; 276 return true;
278 } 277 }
279 278
280 <<
281 //-------------------------------------------- 279 //-----------------------------------------------------------------------
282 280
283 G4bool G4FTFAnnihilation:: 281 G4bool G4FTFAnnihilation::
284 Create3QuarkAntiQuarkStrings( G4VSplitableHadr 282 Create3QuarkAntiQuarkStrings( G4VSplitableHadron* projectile,
285 G4VSplitableHadr 283 G4VSplitableHadron* target,
286 G4VSplitableHadr 284 G4VSplitableHadron*& AdditionalString,
287 G4FTFParameters* 285 G4FTFParameters* theParameters,
288 G4FTFAnnihilatio 286 G4FTFAnnihilation::CommonVariables& common ) const {
289 // Simulation of 3 anti-quark - quark string 287 // Simulation of 3 anti-quark - quark strings creation
290 288
291 #ifdef debugFTFannih 289 #ifdef debugFTFannih
292 G4cout << "Process a, 3 shirt diagram" << G4 290 G4cout << "Process a, 3 shirt diagram" << G4endl;
293 #endif 291 #endif
294 292
295 // Sampling kinematical properties of quark. <<
296 <<
297 const G4int maxNumberOfLoops = 1000; <<
298 G4double MassQ2 = 0.0; // Simp <<
299 // In p <<
300 G4double Quark_Xs[6]; <<
301 G4ThreeVector Quark_Mom[6]; <<
302 <<
303 G4double Alfa_R = 0.5; <<
304 G4double AveragePt2 = 200.0*200.0, maxPtSqua <<
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 <<
315 <<
316 G4double r1 = G4UniformRand(), r2 = G4Un <<
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. <<
323 } <<
324 x3 = 1.0 - x1 - x2; <<
325 <<
326 G4int index = iCase*3; // 0 for anti-ba <<
327 Quark_Xs[index] = x1; Quark_Xs[index+1] <<
328 Product *= (x1*x2*x3); <<
329 } <<
330 <<
331 if ( Product == 0.0 ) continue; <<
332 <<
333 ++NumberOfTries; <<
334 if ( NumberOfTries == 100*(NumberOfTries/1 <<
335 // After a large number of tries, it is <<
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, <<
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 <<
350 Quark_Mom[i] -= PtSum; <<
351 <<
352 G4double val = ( Quark_Mom[i].mag2() + M <<
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( B <<
361 ++loopCounter < maxNumberOfLoops ) <<
362 <<
363 if ( loopCounter >= maxNumberOfLoops ) { <<
364 return false; <<
365 } <<
366 <<
367 G4double DecayMomentum2 = sqr(common.S) + sq <<
368 - 2.0*( common.S*( <<
369 <<
370 G4double WminusTarget = 0.0, WplusProjectile <<
371 WminusTarget = ( common.S - Alfa + Beta + st <<
372 WplusProjectile = common.SqrtS - Beta/Wminus <<
373 <<
374 for ( G4int iCase = 0; iCase < 2; ++iCase ) <<
375 G4int index = iCase*3; // <<
376 G4double w = WplusProjectile; // <<
377 if ( iCase == 1 ) w = - WminusTarget; // <<
378 for ( G4int i = 0; i < 3; ++i ) { <<
379 G4double Pz = w * Quark_Xs[index+i] / 2. <<
380 ( Quark_Mom[index+i].mag2( <<
381 ( 2.0 * w * Quark_Xs[index <<
382 Quark_Mom[index+i].setZ( Pz ); <<
383 } <<
384 } <<
385 <<
386 // Sampling of anti-quark order in projectil 293 // Sampling of anti-quark order in projectile
387 G4int SampledCase = (G4int)G4RandFlat::shoot << 294 G4int SampledCase = G4RandFlat::shootInt( G4long( 6 ) );
388 G4int Tmp1 = 0, Tmp2 = 0; 295 G4int Tmp1 = 0, Tmp2 = 0;
389 switch ( SampledCase ) { 296 switch ( SampledCase ) {
390 case 1 : Tmp1 = common.AQ[1]; common.AQ[1] 297 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] 298 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 299 case 3 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = common.AQ[2];
393 common.AQ[1] = Tmp1; comm 300 common.AQ[1] = Tmp1; common.AQ[2] = Tmp2; break;
394 case 4 : Tmp1 = common.AQ[0]; Tmp2 301 case 4 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = Tmp2;
395 common.AQ[1] = common.AQ[2]; comm 302 common.AQ[1] = common.AQ[2]; common.AQ[2] = Tmp1; break;
396 case 5 : Tmp1 = common.AQ[0]; Tmp2 303 case 5 : Tmp1 = common.AQ[0]; Tmp2 = common.AQ[1]; common.AQ[0] = common.AQ[2];
397 common.AQ[1] = Tmp2; comm 304 common.AQ[1] = Tmp2; common.AQ[2] = Tmp1; break;
398 } 305 }
399 306
400 // Set the string properties 307 // Set the string properties
401 // An anti quark - quark pair can have the q 308 // 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 309 // or a vector meson: the last digit of the PDG code is, respectively, 1 and 3.
403 // For simplicity only scalar is considered 310 // For simplicity only scalar is considered here.
404 G4int NewCode = 0, antiQuark = 0, quark = 0; 311 G4int NewCode = 0, antiQuark = 0, quark = 0;
405 G4ParticleDefinition* TestParticle = nullptr 312 G4ParticleDefinition* TestParticle = nullptr;
406 for ( G4int iString = 0; iString < 3; ++iStr << 313 for ( G4int iString = 0; iString < 3; iString++ ) { // Loop over the 3 string cases
407 if ( iString == 0 ) { 314 if ( iString == 0 ) {
408 antiQuark = common.AQ[0]; quark = commo 315 antiQuark = common.AQ[0]; quark = common.Q[0];
>> 316 projectile->SplitUp();
409 projectile->SetFirstParton( antiQuark ); 317 projectile->SetFirstParton( antiQuark );
410 projectile->SetSecondParton( quark ); 318 projectile->SetSecondParton( quark );
411 projectile->SetStatus( 0 ); 319 projectile->SetStatus( 0 );
412 } else if ( iString == 1 ) { 320 } else if ( iString == 1 ) {
413 quark = common.Q[1]; antiQuark = common 321 quark = common.Q[1]; antiQuark = common.AQ[1];
>> 322 target->SplitUp();
414 target->SetFirstParton( quark ); 323 target->SetFirstParton( quark );
415 target->SetSecondParton( antiQuark ); 324 target->SetSecondParton( antiQuark );
416 target->SetStatus( 0 ); 325 target->SetStatus( 0 );
417 } else { // iString == 2 326 } else { // iString == 2
418 antiQuark = common.AQ[2]; quark = commo 327 antiQuark = common.AQ[2]; quark = common.Q[2];
419 } 328 }
420 G4int absAntiQuark = std::abs( antiQuark ) 329 G4int absAntiQuark = std::abs( antiQuark ), absQuark = std::abs( quark );
421 G4double aKsi = G4UniformRand(); 330 G4double aKsi = G4UniformRand();
422 if ( absAntiQuark == absQuark ) { 331 if ( absAntiQuark == absQuark ) {
423 if ( absAntiQuark != 3 ) { // Not yet c << 332 if ( absAntiQuark != 3 ) {
424 NewCode = 111; // Pi0-meson 333 NewCode = 111; // Pi0-meson
425 if ( aKsi < 0.5 ) { 334 if ( aKsi < 0.5 ) {
426 NewCode = 221; // Eta -meso 335 NewCode = 221; // Eta -meson
427 if ( aKsi < 0.25 ) { 336 if ( aKsi < 0.25 ) {
428 NewCode = 331; // Eta'-meso 337 NewCode = 331; // Eta'-meson
429 } 338 }
430 } 339 }
431 } else { 340 } else {
432 NewCode = 221; // Eta -meso 341 NewCode = 221; // Eta -meson
433 if ( aKsi < 0.5 ) { 342 if ( aKsi < 0.5 ) {
434 NewCode = 331; // Eta'-meso 343 NewCode = 331; // Eta'-meson
435 } 344 }
436 } 345 }
437 } else { // Vector mesons - rho, omega, p << 346 } else {
438 if ( absAntiQuark > absQuark ) { 347 if ( absAntiQuark > absQuark ) {
439 NewCode = absAntiQuark*100 + absQuark* 348 NewCode = absAntiQuark*100 + absQuark*10 + 1; NewCode *= absAntiQuark/antiQuark;
440 } else { 349 } else {
441 NewCode = absQuark*100 + absAntiQuark* 350 NewCode = absQuark*100 + absAntiQuark*10 + 1; NewCode *= absQuark/quark;
442 } 351 }
443 } 352 }
444 if ( iString == 2 ) AdditionalString = new 353 if ( iString == 2 ) AdditionalString = new G4DiffractiveSplitableHadron();
445 TestParticle = G4ParticleTable::GetParticl 354 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
446 if ( ! TestParticle ) return false; 355 if ( ! TestParticle ) return false;
447 if ( iString == 0 ) { 356 if ( iString == 0 ) {
448 projectile->SetDefinition( TestParticle 357 projectile->SetDefinition( TestParticle );
449 theParameters->SetProjMinDiffMass( 0.5 ) << 358 theParameters->SetProjMinDiffMass( 0.5 );
450 theParameters->SetProjMinNonDiffMass( 0. << 359 theParameters->SetProjMinNonDiffMass( 0.5 );
451 } else if ( iString == 1 ) { 360 } else if ( iString == 1 ) {
452 target->SetDefinition( TestParticle ); 361 target->SetDefinition( TestParticle );
453 theParameters->SetTarMinDiffMass( 0.5 ); 362 theParameters->SetTarMinDiffMass( 0.5 );
454 theParameters->SetTarMinNonDiffMass( 0.5 363 theParameters->SetTarMinNonDiffMass( 0.5 );
455 } else { // iString == 2 364 } else { // iString == 2
456 AdditionalString->SetDefinition( TestPar 365 AdditionalString->SetDefinition( TestParticle );
>> 366 AdditionalString->SplitUp();
457 AdditionalString->SetFirstParton( common 367 AdditionalString->SetFirstParton( common.AQ[2] );
458 AdditionalString->SetSecondParton( commo 368 AdditionalString->SetSecondParton( common.Q[2] );
459 AdditionalString->SetStatus( 0 ); 369 AdditionalString->SetStatus( 0 );
460 } 370 }
461 } // End of the for loop over the 3 string 371 } // End of the for loop over the 3 string cases
462 372
>> 373 // Sampling kinematical properties:
463 // 1st string AQ[0]-Q[0], 2nd string AQ[1]-Q 374 // 1st string AQ[0]-Q[0], 2nd string AQ[1]-Q[1], 3rd string AQ[2]-Q[2]
464 375
>> 376 G4ThreeVector Quark_Mom[6];
>> 377 G4double ModMom2[6];
>> 378 G4double AveragePt2 = 200.0*200.0, maxPtSquare = common.S, SumMt = 0.0, MassQ2 = 0.0,
>> 379 ScaleFactor = 1.0;
>> 380 G4int NumberOfTries = 0, loopCounter = 0;
>> 381 const G4int maxNumberOfLoops = 1000;
>> 382 do {
>> 383 NumberOfTries++;
>> 384 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
>> 385 // At large number of tries it would be better to reduce the values of <Pt^2>
>> 386 ScaleFactor /= 2.0;
>> 387 AveragePt2 *= ScaleFactor;
>> 388 }
>> 389 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
>> 390 for ( G4int i = 0; i < 6; i++ ) {
>> 391 Quark_Mom [i] = GaussianPt( AveragePt2, maxPtSquare );
>> 392 PtSum += Quark_Mom[i];
>> 393 }
>> 394 PtSum /= 6.0;
>> 395 SumMt = 0.0;
>> 396 for ( G4int i = 0; i < 6; i++ ) {
>> 397 Quark_Mom[i] -= PtSum;
>> 398 ModMom2[i] = Quark_Mom[i].mag2();
>> 399 SumMt += std::sqrt( ModMom2[i] + MassQ2 );
>> 400 }
>> 401 } while ( ( SumMt > common.SqrtS ) &&
>> 402 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
>> 403 if ( loopCounter >= maxNumberOfLoops ) {
>> 404 return false;
>> 405 }
>> 406
>> 407 // Sampling X's of anti-baryon and baryon
>> 408 G4double WminusTarget = 0.0, WplusProjectile = 0.0;
>> 409 G4double Alfa_R = 0.5; ScaleFactor = 1.0;
>> 410 G4bool Success = true;
>> 411 NumberOfTries = 0; loopCounter = 0;
>> 412 do {
>> 413 Success = true;
>> 414 NumberOfTries++;
>> 415 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
>> 416 // At large number of tries it would be better to reduce the values of Pt's
>> 417 ScaleFactor /= 2.0;
>> 418 }
>> 419 G4double Alfa = 0.0, Beta = 0.0;
>> 420 for ( G4int iCase = 0; iCase < 2; iCase++ ) { // anti-baryon (1st case), baryon (2nd case)
>> 421 G4double x1 = 0.0, x2 = 0.0;
>> 422 G4double r1 = G4UniformRand(), r2 = G4UniformRand();
>> 423 if ( Alfa_R == 1.0 ) {
>> 424 x1 = 1.0 - std::sqrt( r1 );
>> 425 x2 = (1.0 - x1) * r2;
>> 426 } else {
>> 427 x1 = sqr( r1 );
>> 428 x2 = (1.0 - x1) * sqr( std::sin( pi/2.0*r2 ) );
>> 429 }
>> 430 G4double x3 = 1.0 - x1 - x2;
>> 431 G4int index = iCase*3; // 0 for anti-baryon, 3 for baryon
>> 432 Quark_Mom[index].setZ( x1 ); Quark_Mom[index+1].setZ( x2 ); Quark_Mom[index+2].setZ( x3 );
>> 433 for ( G4int i = 0; i < 3; i++ ) { // Loop over the 3 (anti-)quarks
>> 434 if ( Quark_Mom[index+i].getZ() != 0.0 ) {
>> 435 G4double val = ( ScaleFactor * ModMom2[index+i] + MassQ2 ) / Quark_Mom[index+i].getZ();
>> 436 if ( iCase == 0 ) { // anti-baryon
>> 437 Alfa += val;
>> 438 } else { // baryon (iCase == 1)
>> 439 Beta += val;
>> 440 }
>> 441 } else {
>> 442 Success = false;
>> 443 }
>> 444 }
>> 445 }
>> 446 if ( ! Success ) continue;
>> 447 if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > common.SqrtS ) {
>> 448 Success = false;
>> 449 continue;
>> 450 }
>> 451 G4double DecayMomentum2 = sqr(common.S) + sqr(Alfa) + sqr(Beta)
>> 452 - 2.0*( common.S*(Alfa + Beta) + Alfa*Beta );
>> 453 WminusTarget = ( common.S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.SqrtS;
>> 454 WplusProjectile = common.SqrtS - Beta/WminusTarget;
>> 455 } while ( ( ! Success ) &&
>> 456 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
>> 457 if ( loopCounter >= maxNumberOfLoops ) {
>> 458 return false;
>> 459 }
>> 460
>> 461 G4double SqrtScaleF = std::sqrt( ScaleFactor );
>> 462 for ( G4int iCase = 0; iCase < 2; iCase++ ) { // anti-baryon (1st case), baryon (2nd case)
>> 463 G4int index = iCase*3; // 0 for anti-baryon, 3 for baryon
>> 464 G4double w = WplusProjectile; // for anti-baryon
>> 465 if ( iCase == 1 ) w = - WminusTarget; // for baryon
>> 466 for ( G4int i = 0; i < 3; i++ ) {
>> 467 G4double Pz = w * Quark_Mom[index+i].getZ() / 2.0 -
>> 468 ( ScaleFactor * ModMom2[index+i] + MassQ2 ) /
>> 469 ( 2.0 * w * Quark_Mom[index+i].getZ() );
>> 470 Quark_Mom[index+i].setZ( Pz );
>> 471 if ( ScaleFactor != 1.0 ) {
>> 472 Quark_Mom[index+i].setX( SqrtScaleF * Quark_Mom[index+i].getX() );
>> 473 Quark_Mom[index+i].setY( SqrtScaleF * Quark_Mom[index+i].getY() );
>> 474 }
>> 475 }
>> 476 }
465 G4LorentzVector Pstring1, Pstring2, Pstring3 477 G4LorentzVector Pstring1, Pstring2, Pstring3;
466 G4int QuarkOrder[3] = { 0 }; <<
467 G4double YstringMax = 0.0, YstringMin = 0.0; 478 G4double YstringMax = 0.0, YstringMin = 0.0;
468 for ( G4int i = 0; i < 3; ++i ) { << 479 for ( G4int i = 0; i < 3; i++ ) {
469 G4ThreeVector tmp = Quark_Mom[i] + Quark_M 480 G4ThreeVector tmp = Quark_Mom[i] + Quark_Mom[i+3];
470 G4LorentzVector Pstring( tmp, std::sqrt( Q 481 G4LorentzVector Pstring( tmp, std::sqrt( Quark_Mom[i].mag2() + MassQ2 ) +
471 std::sqrt( Q 482 std::sqrt( Quark_Mom[i+3].mag2() + MassQ2 ) );
>> 483 //AR-Jun2018 if ( common.RotateStrings ) Pstring *= common.RandomRotation;
472 // Add protection for rapidity = 0.5*ln( 484 // Add protection for rapidity = 0.5*ln( (E+Pz)/(E-Pz) )
473 G4double Ystring = 0.0; 485 G4double Ystring = 0.0;
474 if ( Pstring.e() > 1.0e-30 ) { 486 if ( Pstring.e() > 1.0e-30 ) {
475 if ( Pstring.e() + Pstring.pz() < 1.0e-3 487 if ( Pstring.e() + Pstring.pz() < 1.0e-30 ) { // Very small numerator in the logarithm
476 Ystring = -1.0e30; // A very large n 488 Ystring = -1.0e30; // A very large negative value (E ~= -Pz)
477 if ( Pstring.e() - Pstring.pz() < 1.0e 489 if ( Pstring.e() - Pstring.pz() < 1.0e-30 ) { // Very small denominator in the logarithm
478 Ystring = 1.0e30; // A very large p 490 Ystring = 1.0e30; // A very large positive value (E ~= Pz)
479 } else { // Normal case 491 } else { // Normal case
480 Ystring = Pstring.rapidity(); 492 Ystring = Pstring.rapidity();
481 } 493 }
482 } 494 }
483 } 495 }
484 // Keep ordering in rapidity: "1" highest, 496 // Keep ordering in rapidity: "1" highest, "2" middle, "3" smallest
485 if ( i == 0 ) { 497 if ( i == 0 ) {
486 Pstring1 = Pstring; YstringMax = Yst << 498 Pstring1 = Pstring; YstringMax = Ystring;
487 QuarkOrder[0] = 0; <<
488 } else if ( i == 1 ) { 499 } else if ( i == 1 ) {
489 if ( Ystring > YstringMax ) { 500 if ( Ystring > YstringMax ) {
490 Pstring2 = Pstring1; YstringMin = Yst << 501 Pstring2 = Pstring1; YstringMin = YstringMax;
491 Pstring1 = Pstring; YstringMax = Yst << 502 Pstring1 = Pstring; YstringMax = Ystring;
492 QuarkOrder[0] = 1; QuarkOrder[1] = 0; <<
493 } else { 503 } else {
494 Pstring2 = Pstring; YstringMin = Yst << 504 Pstring2 = Pstring; YstringMin = Ystring;
495 QuarkOrder[1] = 1; <<
496 } 505 }
497 } else { // i == 2 506 } else { // i == 2
498 if ( Ystring > YstringMax ) { 507 if ( Ystring > YstringMax ) {
499 Pstring3 = Pstring2; 508 Pstring3 = Pstring2;
500 Pstring2 = Pstring1; 509 Pstring2 = Pstring1;
501 Pstring1 = Pstring; 510 Pstring1 = Pstring;
502 QuarkOrder[1] = QuarkOrder[0]; <<
503 QuarkOrder[2] = QuarkOrder[1]; <<
504 QuarkOrder[0] = 2; <<
505 } else if ( Ystring > YstringMin ) { 511 } else if ( Ystring > YstringMin ) {
506 Pstring3 = Pstring2; 512 Pstring3 = Pstring2;
507 Pstring2 = Pstring; 513 Pstring2 = Pstring;
508 } else { 514 } else {
509 Pstring3 = Pstring; << 515 Pstring3 = Pstring;
510 QuarkOrder[2] = 2; <<
511 } 516 }
512 } 517 }
513 } 518 }
514 <<
515 G4LorentzVector Quark_4Mom[6]; <<
516 for ( G4int i = 0; i < 6; ++i ) { <<
517 Quark_4Mom[i] = G4LorentzVector( Quark_Mom <<
518 if ( common.RotateStrings ) Quark_4Mom[i] <<
519 Quark_4Mom[i].transform( common.toLab ); <<
520 } <<
521 <<
522 projectile->Splitting(); <<
523 projectile->GetNextAntiParton()->Set4Momentu <<
524 projectile->GetNextParton()->Set4Momentum( Q <<
525 <<
526 target->Splitting(); <<
527 target->GetNextParton()->Set4Momentum( Quark <<
528 target->GetNextAntiParton()->Set4Momentum( Q <<
529 <<
530 AdditionalString->Splitting(); <<
531 AdditionalString->GetNextAntiParton()->Set4M <<
532 AdditionalString->GetNextParton()->Set4Momen <<
533 <<
534 common.Pprojectile = Pstring1; // 519 common.Pprojectile = Pstring1; // Highest rapidity
535 common.Ptarget = Pstring3; // 520 common.Ptarget = Pstring3; // Lowest rapidity
536 G4LorentzVector LeftString( Pstring2 ); // 521 G4LorentzVector LeftString( Pstring2 ); // Middle rapidity
537 522
538 if ( common.RotateStrings ) { << 523 if ( common.RotateStrings ) { //AR-Jun2018
539 common.Pprojectile *= common.RandomRotatio 524 common.Pprojectile *= common.RandomRotation;
540 common.Ptarget *= common.RandomRotatio <<
541 LeftString *= common.RandomRotatio 525 LeftString *= common.RandomRotation;
>> 526 common.Ptarget *= common.RandomRotation;
542 } 527 }
543 528
544 common.Pprojectile.transform( common.toLab ) 529 common.Pprojectile.transform( common.toLab );
545 common.Ptarget.transform( common.toLab ); <<
546 LeftString.transform( common.toLab ); 530 LeftString.transform( common.toLab );
547 << 531 common.Ptarget.transform( common.toLab );
>> 532
548 // Calculation of the creation time 533 // Calculation of the creation time
549 // Creation time and position of target nucl 534 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
550 projectile->SetTimeOfCreation( target->GetTi 535 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
551 projectile->SetPosition( target->GetPosition 536 projectile->SetPosition( target->GetPosition() );
552 AdditionalString->SetTimeOfCreation( target- 537 AdditionalString->SetTimeOfCreation( target->GetTimeOfCreation() );
553 AdditionalString->SetPosition( target->GetPo 538 AdditionalString->SetPosition( target->GetPosition() );
554 539
555 projectile->Set4Momentum( common.Pprojectile 540 projectile->Set4Momentum( common.Pprojectile );
556 AdditionalString->Set4Momentum( LeftString ) 541 AdditionalString->Set4Momentum( LeftString );
557 target->Set4Momentum( common.Ptarget ); 542 target->Set4Momentum( common.Ptarget );
558 <<
559 projectile->IncrementCollisionCount( 1 ); 543 projectile->IncrementCollisionCount( 1 );
560 AdditionalString->IncrementCollisionCount( 1 544 AdditionalString->IncrementCollisionCount( 1 );
561 target->IncrementCollisionCount( 1 ); 545 target->IncrementCollisionCount( 1 );
562 546
563 return true; 547 return true;
564 } 548 }
565 549
566 <<
567 //-------------------------------------------- 550 //-----------------------------------------------------------------------
568 551
569 G4int G4FTFAnnihilation:: 552 G4int G4FTFAnnihilation::
570 Create1DiquarkAntiDiquarkString( G4VSplitableH 553 Create1DiquarkAntiDiquarkString( G4VSplitableHadron* projectile,
571 G4VSplitableH 554 G4VSplitableHadron* target,
572 G4FTFAnnihila 555 G4FTFAnnihilation::CommonVariables& common ) const {
573 // Simulation of anti-diquark-diquark string 556 // Simulation of anti-diquark-diquark string creation.
574 // This method returns an integer code - ins 557 // This method returns an integer code - instead of a boolean, with the following meaning:
575 // "0" : successfully ended and nothing el 558 // "0" : successfully ended and nothing else needs to be done;
576 // "1" : successfully completed, but the w 559 // "1" : successfully completed, but the work needs to be continued;
577 // "99" : unsuccessfully ended, nothing els 560 // "99" : unsuccessfully ended, nothing else can be done.
578 561
579 #ifdef debugFTFannih 562 #ifdef debugFTFannih
580 G4cout << "Process b, quark - anti-quark ann 563 G4cout << "Process b, quark - anti-quark annihilation, di-q - anti-di-q string" << G4endl;
581 #endif 564 #endif
582 565
583 G4int CandidatsN = 0, CandAQ[9][2] = {}, Can << 566 G4int CandidatsN = 0, CandAQ[9][2], CandQ[9][2];
584 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // << 567 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 ) { // << 568 for ( G4int iQ = 0; iQ < 3; iQ++ ) { // index of the 3 constituent quarks of the nucleon target
586 if ( -common.AQ[iAQ] == common.Q[iQ] ) { 569 if ( -common.AQ[iAQ] == common.Q[iQ] ) { // antiquark - quark that can annihilate
587 // Here "0", "1", "2" means, respectiv 570 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
588 // of the (anti-baryon) projectile or 571 // of the (anti-baryon) projectile or (nucleon) target.
589 if ( iAQ == 0 ) { CandAQ[CandidatsN][0 572 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
590 if ( iAQ == 1 ) { CandAQ[CandidatsN][0 573 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
591 if ( iAQ == 2 ) { CandAQ[CandidatsN][0 574 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
592 if ( iQ == 0 ) { CandQ[CandidatsN][0] 575 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
593 if ( iQ == 1 ) { CandQ[CandidatsN][0] 576 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
594 if ( iQ == 2 ) { CandQ[CandidatsN][0] 577 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
595 ++CandidatsN; << 578 CandidatsN++;
596 } 579 }
597 } 580 }
598 } 581 }
599 582
600 // Remaining two (anti-)quarks that form the 583 // Remaining two (anti-)quarks that form the (anti-)diquark
601 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, 584 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, LeftQ2 = 0;
602 if ( CandidatsN != 0 ) { 585 if ( CandidatsN != 0 ) {
603 G4int SampledCase = (G4int)G4RandFlat::sho << 586 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
604 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0 587 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0] ];
605 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1 588 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1] ];
606 LeftQ1 = common.Q[ CandQ[SampledCase][0] 589 LeftQ1 = common.Q[ CandQ[SampledCase][0] ];
607 LeftQ2 = common.Q[ CandQ[SampledCase][1] 590 LeftQ2 = common.Q[ CandQ[SampledCase][1] ];
608 591
609 // Build anti-diquark and diquark : the la 592 // Build anti-diquark and diquark : the last digit can be either 3 - for all combinations
610 // of anti-quark - anti-quark and quark - 593 // of anti-quark - anti-quark and quark - quark - or 1 - only when the two anti-quarks
611 // or quarks are different. For simplicity 594 // or quarks are different. For simplicity, only 3 is considered.
612 G4int Anti_DQ = 0, DQ = 0; 595 G4int Anti_DQ = 0, DQ = 0;
613 if ( std::abs( LeftAQ1 ) > std::abs( LeftA 596 if ( std::abs( LeftAQ1 ) > std::abs( LeftAQ2 ) ) {
614 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3 597 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3;
615 } else { 598 } else {
616 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3 599 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3;
617 } 600 }
618 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 601 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 ) ) {
619 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3; 602 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3;
620 } else { 603 } else {
621 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3; 604 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3;
622 } 605 }
623 606
624 // Set the string properties 607 // Set the string properties
>> 608 projectile->SplitUp();
625 projectile->SetFirstParton( DQ ); 609 projectile->SetFirstParton( DQ );
626 projectile->SetSecondParton( Anti_DQ ); 610 projectile->SetSecondParton( Anti_DQ );
627 611
628 // It is assumed that quark and di-quark m <<
629 G4LorentzVector Pquark = G4LorentzVector( <<
630 G4LorentzVector Paquark = G4LorentzVector( <<
631 <<
632 if ( common.RotateStrings ) { 612 if ( common.RotateStrings ) {
>> 613 G4LorentzVector Pquark = G4LorentzVector( 0.0, 0.0, common.SqrtS/2.0, common.SqrtS/2.0 );
633 Pquark *= common.RandomRotation; 614 Pquark *= common.RandomRotation;
>> 615 G4LorentzVector Paquark = G4LorentzVector( 0.0, 0.0, -common.SqrtS/2.0, common.SqrtS/2.0 );
634 Paquark *= common.RandomRotation; 616 Paquark *= common.RandomRotation;
>> 617 Pquark.transform( common.toLab ); projectile->GetNextParton()->Set4Momentum( Pquark );
>> 618 Paquark.transform( common.toLab ); projectile->GetNextAntiParton()->Set4Momentum( Paquark );
635 } 619 }
636 620
637 Pquark.transform( common.toLab ); <<
638 Paquark.transform( common.toLab ); <<
639 <<
640 projectile->GetNextParton()->Set4Momentum( <<
641 projectile->GetNextAntiParton()->Set4Momen <<
642 <<
643 projectile->Splitting(); <<
644 <<
645 projectile->SetStatus( 0 ); 621 projectile->SetStatus( 0 );
646 target->SetStatus( 4 ); // The target nuc 622 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4
647 common.Pprojectile.setPx( 0.0 ); 623 common.Pprojectile.setPx( 0.0 );
648 common.Pprojectile.setPy( 0.0 ); 624 common.Pprojectile.setPy( 0.0 );
649 common.Pprojectile.setPz( 0.0 ); 625 common.Pprojectile.setPz( 0.0 );
650 common.Pprojectile.setE( common.SqrtS ); 626 common.Pprojectile.setE( common.SqrtS );
651 common.Pprojectile.transform( common.toLab 627 common.Pprojectile.transform( common.toLab );
652 628
653 // Calculation of the creation time 629 // Calculation of the creation time
654 // Creation time and position of target nu 630 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
655 projectile->SetTimeOfCreation( target->Get 631 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
656 projectile->SetPosition( target->GetPositi 632 projectile->SetPosition( target->GetPosition() );
657 projectile->Set4Momentum( common.Pprojecti 633 projectile->Set4Momentum( common.Pprojectile );
658 <<
659 projectile->IncrementCollisionCount( 1 ); 634 projectile->IncrementCollisionCount( 1 );
660 target->IncrementCollisionCount( 1 ); 635 target->IncrementCollisionCount( 1 );
661 636
662 return 0; // Completed successfully: noth 637 return 0; // Completed successfully: nothing else to be done
663 } // End of if ( CandidatsN != 0 ) 638 } // End of if ( CandidatsN != 0 )
664 639
665 // If we allow the string to interact with o <<
666 // set up MinDiffrMass in Parameters, and as <<
667 <<
668 return 1; // Successfully ended, but the wo 640 return 1; // Successfully ended, but the work is not over
669 } 641 }
670 642
671 <<
672 //-------------------------------------------- 643 //-----------------------------------------------------------------------
673 644
674 G4int G4FTFAnnihilation:: 645 G4int G4FTFAnnihilation::
675 Create2QuarkAntiQuarkStrings( G4VSplitableHadr 646 Create2QuarkAntiQuarkStrings( G4VSplitableHadron* projectile,
676 G4VSplitableHadr 647 G4VSplitableHadron* target,
677 G4FTFParameters* 648 G4FTFParameters* theParameters,
678 G4FTFAnnihilatio 649 G4FTFAnnihilation::CommonVariables& common ) const {
679 // Simulation of 2 anti-quark-quark strings 650 // Simulation of 2 anti-quark-quark strings creation.
680 // This method returns an integer code - ins 651 // This method returns an integer code - instead of a boolean, with the following meaning:
681 // "0" : successfully ended and nothing el 652 // "0" : successfully ended and nothing else needs to be done;
682 // "1" : successfully completed, but the w 653 // "1" : successfully completed, but the work needs to be continued;
683 // "99" : unsuccessfully ended, nothing els 654 // "99" : unsuccessfully ended, nothing else can be done.
684 655
685 #ifdef debugFTFannih 656 #ifdef debugFTFannih
686 G4cout << "Process c, quark - anti-quark and 657 G4cout << "Process c, quark - anti-quark and string junctions annihilation, 2 strings left."
687 << G4endl; 658 << G4endl;
688 #endif 659 #endif
689 660
690 // Sampling kinematical properties: 1st stri << 661 G4int CandidatsN = 0, CandAQ[9][2], CandQ[9][2];
691 G4ThreeVector Quark_Mom[4]; <<
692 G4double Quark_Xs[4]; <<
693 G4double AveragePt2 = 200.0*200.0, maxPtSqua <<
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 <<
698 do { <<
699 // Sampling X's of the 2 quarks and 2 anti <<
700 <<
701 G4double Product = 1.0; <<
702 for ( G4int iCase = 0; iCase < 2; ++iCase <<
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 fir <<
714 Quark_Xs[index] = x ; Quark_Xs[index+1] <<
715 Product *= x*(1.0-x); <<
716 } <<
717 <<
718 if ( Product == 0.0 ) continue; <<
719 <<
720 ++NumberOfTries; <<
721 if ( NumberOfTries == 100*(NumberOfTries/1 <<
722 // After a large number of tries, it is <<
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, m <<
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]. <<
743 if ( iCase == 0 ) { // first string <<
744 Alfa += val; <<
745 } else { // second strin <<
746 Beta += val; <<
747 } <<
748 } <<
749 } <<
750 <<
751 } while ( ( std::sqrt( Alfa ) + std::sqrt( B <<
752 ++loopCounter < maxNumberOfLoops ) <<
753 <<
754 if ( loopCounter >= maxNumberOfLoops ) { <<
755 return 99; // unsuccessfully ended, nothi <<
756 } <<
757 <<
758 G4double DecayMomentum2 = sqr(common.S) + sq <<
759 - 2.0*( common.S*( <<
760 WminusTarget = ( common.S - Alfa + Beta + st <<
761 WplusProjectile = common.SqrtS - Beta/Wminus <<
762 <<
763 for ( G4int iCase = 0; iCase < 2; ++iCase ) <<
764 G4int index = iCase*2; // 0 for the first <<
765 for ( G4int i = 0; i < 2; ++i ) { <<
766 G4double w = WplusProjectile; / <<
767 if ( iCase == 1 ) w = - WminusTarget; / <<
768 G4double Pz = w * Quark_Xs[index+i] / 2. <<
769 - ( Quark_Mom[index+i].mag <<
770 ( 2.0 * w * Quark_Xs[ind <<
771 Quark_Mom[index+i].setZ( Pz ); <<
772 } <<
773 } <<
774 <<
775 G4int CandidatsN = 0, CandAQ[9][2] = {}, Can <<
776 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, 662 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, LeftQ2 = 0;
777 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // << 663 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 ) { // << 664 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] ) { 665 if ( -common.AQ[iAQ] == common.Q[iQ] ) { // antiquark - quark that can annihilate
780 // Here "0", "1", "2" means, respectiv 666 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
781 // of the (anti-baryon) projectile or 667 // of the (anti-baryon) projectile or (nucleon) target.
782 if ( iAQ == 0 ) { CandAQ[CandidatsN][0 668 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; }
783 if ( iAQ == 1 ) { CandAQ[CandidatsN][0 669 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; }
784 if ( iAQ == 2 ) { CandAQ[CandidatsN][0 670 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; }
785 if ( iQ == 0 ) { CandQ[CandidatsN][0] 671 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; }
786 if ( iQ == 1 ) { CandQ[CandidatsN][0] 672 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; }
787 if ( iQ == 2 ) { CandQ[CandidatsN][0] 673 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; }
788 ++CandidatsN; << 674 CandidatsN++;
789 } 675 }
790 } 676 }
791 } 677 }
792 678
793 if ( CandidatsN != 0 ) { 679 if ( CandidatsN != 0 ) {
794 G4int SampledCase = (G4int)G4RandFlat::sho << 680 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
795 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0 681 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0] ];
796 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1 682 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1] ];
797 if ( G4UniformRand() < 0.5 ) { 683 if ( G4UniformRand() < 0.5 ) {
798 LeftQ1 = common.Q[ CandQ[SampledCase][0] 684 LeftQ1 = common.Q[ CandQ[SampledCase][0] ];
799 LeftQ2 = common.Q[ CandQ[SampledCase][1] 685 LeftQ2 = common.Q[ CandQ[SampledCase][1] ];
800 } else { 686 } else {
801 LeftQ2 = common.Q[ CandQ[SampledCase][0] 687 LeftQ2 = common.Q[ CandQ[SampledCase][0] ];
802 LeftQ1 = common.Q[ CandQ[SampledCase][1] 688 LeftQ1 = common.Q[ CandQ[SampledCase][1] ];
803 } 689 }
804 690
805 // Set the string properties 691 // Set the string properties
806 // An anti quark - quark pair can have the 692 // 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 693 // or a vector meson: the last digit of the PDG code is, respectively, 1 and 3.
808 // For simplicity only scalar is considere 694 // For simplicity only scalar is considered here.
809 G4int NewCode = 0, antiQuark = 0, quark = 695 G4int NewCode = 0, antiQuark = 0, quark = 0;
810 G4ParticleDefinition* TestParticle = nullp 696 G4ParticleDefinition* TestParticle = nullptr;
811 for ( G4int iString = 0; iString < 2; ++iS << 697 for ( G4int iString = 0; iString < 2; iString++ ) { // Loop over the 2 string cases
812 if ( iString == 0 ) { 698 if ( iString == 0 ) {
813 antiQuark = LeftAQ1; quark = LeftQ1; 699 antiQuark = LeftAQ1; quark = LeftQ1;
>> 700 projectile->SplitUp();
814 projectile->SetFirstParton( antiQuark 701 projectile->SetFirstParton( antiQuark );
815 projectile->SetSecondParton( quark ); 702 projectile->SetSecondParton( quark );
816 projectile->SetStatus( 0 ); 703 projectile->SetStatus( 0 );
817 } else { // iString == 1 704 } else { // iString == 1
818 quark = LeftQ2; antiQuark = LeftAQ2; 705 quark = LeftQ2; antiQuark = LeftAQ2;
>> 706 target->SplitUp();
819 target->SetFirstParton( quark ); 707 target->SetFirstParton( quark );
820 target->SetSecondParton( antiQuark ); 708 target->SetSecondParton( antiQuark );
821 target->SetStatus( 0 ); 709 target->SetStatus( 0 );
822 } 710 }
823 G4int absAntiQuark = std::abs( antiQuark 711 G4int absAntiQuark = std::abs( antiQuark ), absQuark = std::abs( quark );
824 G4double aKsi = G4UniformRand(); 712 G4double aKsi = G4UniformRand();
825 if ( absAntiQuark == absQuark ) { 713 if ( absAntiQuark == absQuark ) {
826 if ( absAntiQuark != 3 ) { 714 if ( absAntiQuark != 3 ) {
827 NewCode = 111; // Pi0-meson 715 NewCode = 111; // Pi0-meson
828 if ( aKsi < 0.5 ) { 716 if ( aKsi < 0.5 ) {
829 NewCode = 221; // Eta -meso 717 NewCode = 221; // Eta -meson
830 if ( aKsi < 0.25 ) { 718 if ( aKsi < 0.25 ) {
831 NewCode = 331; // Eta'-meso 719 NewCode = 331; // Eta'-meson
832 } 720 }
833 } 721 }
834 } else { 722 } else {
835 NewCode = 221; // Eta -meso 723 NewCode = 221; // Eta -meson
836 if ( aKsi < 0.5 ) { 724 if ( aKsi < 0.5 ) {
837 NewCode = 331; // Eta'-meso 725 NewCode = 331; // Eta'-meson
838 } 726 }
839 } 727 }
840 } else { 728 } else {
841 if ( absAntiQuark > absQuark ) { 729 if ( absAntiQuark > absQuark ) {
842 NewCode = absAntiQuark*100 + absQuar 730 NewCode = absAntiQuark*100 + absQuark*10 + 1; NewCode *= absAntiQuark/antiQuark;
843 } else { 731 } else {
844 NewCode = absQuark*100 + absAntiQuar 732 NewCode = absQuark*100 + absAntiQuark*10 + 1; NewCode *= absQuark/quark;
845 } 733 }
846 } 734 }
847 TestParticle = G4ParticleTable::GetParti 735 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
848 if ( ! TestParticle ) return 99; // uns 736 if ( ! TestParticle ) return 99; // unsuccessfully ended, nothing else can be done
849 if ( iString == 0 ) { 737 if ( iString == 0 ) {
850 projectile->SetDefinition( TestParticl 738 projectile->SetDefinition( TestParticle );
851 theParameters->SetProjMinDiffMass( 0.5 739 theParameters->SetProjMinDiffMass( 0.5 );
852 theParameters->SetProjMinNonDiffMass( 740 theParameters->SetProjMinNonDiffMass( 0.5 );
853 } else { // iString == 1 741 } else { // iString == 1
854 target->SetDefinition( TestParticle ); 742 target->SetDefinition( TestParticle );
855 theParameters->SetTarMinDiffMass( 0.5 743 theParameters->SetTarMinDiffMass( 0.5 );
856 theParameters->SetTarMinNonDiffMass( 0 744 theParameters->SetTarMinNonDiffMass( 0.5 );
857 } 745 }
858 } // End of loop over the 2 string cases 746 } // End of loop over the 2 string cases
859 747
860 G4int QuarkOrder[2]; << 748 // Sampling kinematical properties: 1st string LeftAQ1-LeftQ1, 2nd string LeftAQ2-LeftQ2
>> 749 G4ThreeVector Quark_Mom[4];
>> 750 G4double ModMom2[4];
>> 751 G4double AveragePt2 = 200.0*200.0, maxPtSquare = common.S, SumMt = 0.0, MassQ2 = 0.0,
>> 752 ScaleFactor = 1.0;
>> 753 G4int NumberOfTries = 0, loopCounter = 0;
>> 754 const G4int maxNumberOfLoops = 1000;
>> 755 do {
>> 756 NumberOfTries++;
>> 757 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
>> 758 // At large number of tries it would be better to reduce the values of <Pt^2>
>> 759 ScaleFactor /= 2.0;
>> 760 AveragePt2 *= ScaleFactor;
>> 761 }
>> 762 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
>> 763 for( G4int i = 0; i < 4; i++ ) {
>> 764 Quark_Mom[i] = GaussianPt( AveragePt2, maxPtSquare );
>> 765 PtSum += Quark_Mom[i];
>> 766 }
>> 767 PtSum /= 4.0;
>> 768 SumMt = 0.0;
>> 769 for ( G4int i = 0; i < 4; i++ ) {
>> 770 Quark_Mom[i] -= PtSum;
>> 771 ModMom2[i] = Quark_Mom[i].mag2();
>> 772 SumMt += std::sqrt( ModMom2[i] + MassQ2 );
>> 773 }
>> 774 } while ( ( SumMt > common.SqrtS ) &&
>> 775 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
>> 776 if ( loopCounter >= maxNumberOfLoops ) {
>> 777 return 99; // unsuccessfully ended, nothing else can be done
>> 778 }
>> 779
>> 780 // Sampling X's of the two strings
>> 781 G4double WminusTarget = 0.0, WplusProjectile = 0.0, Alfa_R = 0.5; ScaleFactor = 1.0;
>> 782 G4bool Success = true;
>> 783 NumberOfTries = 0, loopCounter = 0;
>> 784 do {
>> 785 Success = true;
>> 786 NumberOfTries++;
>> 787 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
>> 788 // At large number of tries it would be better to reduce the values of Pt's
>> 789 ScaleFactor /= 2.0;
>> 790 }
>> 791 G4double Alfa = 0.0, Beta = 0.0;
>> 792 for ( G4int iCase = 0; iCase < 2; iCase++ ) { // Loop over the two strings
>> 793 G4double x = 0.0, r = G4UniformRand();
>> 794 if ( Alfa_R == 1.0 ) {
>> 795 if ( iCase == 0 ) { // first string
>> 796 x = std::sqrt( r );
>> 797 } else { // second string
>> 798 x = 1.0 - std::sqrt( r);
>> 799 }
>> 800 } else {
>> 801 x = sqr( std::sin( pi/2.0*r ) );
>> 802 }
>> 803 G4int index = iCase*2; // 0 for the first string, 2 for the second string
>> 804 Quark_Mom[index].setZ( x ); Quark_Mom[index+1].setZ( 1.0 - x );
>> 805 for ( G4int i = 0; i < 2; i++ ) {
>> 806 if ( Quark_Mom[i].getZ() != 0.0 ) {
>> 807 G4double val = ( ScaleFactor * ModMom2[index+i] + MassQ2 ) / Quark_Mom[index+i].getZ();
>> 808 if ( iCase == 0 ) { // first string
>> 809 Alfa += val;
>> 810 } else { // second string
>> 811 Beta += val;
>> 812 }
>> 813 } else {
>> 814 Success = false;
>> 815 }
>> 816 }
>> 817 }
>> 818 if ( ! Success ) continue;
>> 819 if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > common.SqrtS ) {
>> 820 Success = false;
>> 821 continue;
>> 822 }
>> 823 G4double DecayMomentum2 = sqr(common.S) + sqr(Alfa) + sqr(Beta)
>> 824 - 2.0*( common.S*(Alfa + Beta) + Alfa*Beta );
>> 825 WminusTarget = ( common.S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / common.SqrtS;
>> 826 WplusProjectile = common.SqrtS - Beta/WminusTarget;
>> 827 } while ( ( ! Success ) &&
>> 828 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
>> 829 if ( loopCounter >= maxNumberOfLoops ) {
>> 830 return 99; // unsuccessfully ended, nothing else can be done
>> 831 }
>> 832
>> 833 G4double SqrtScaleF = std::sqrt( ScaleFactor );
861 G4LorentzVector Pstring1, Pstring2; 834 G4LorentzVector Pstring1, Pstring2;
862 G4double Ystring1 = 0.0, Ystring2 = 0.0; 835 G4double Ystring1 = 0.0, Ystring2 = 0.0;
863 << 836 for ( G4int iCase = 0; iCase < 2; iCase++ ) { // Loop over the two strings
864 for ( G4int iCase = 0; iCase < 2; ++iCase << 837 G4int index = iCase*2; // 0 for the first string, 2 for the second string
>> 838 for ( G4int i = 0; i < 2; i++ ) {
>> 839 G4double w = WplusProjectile; // For the first string
>> 840 if ( iCase == 1 ) w = - WminusTarget; // For the second string
>> 841 G4double Pz = w * Quark_Mom[index+i].getZ() / 2.0
>> 842 - ( ScaleFactor * ModMom2[index+i] + MassQ2 ) /
>> 843 ( 2.0 * w * Quark_Mom[index+i].getZ() );
>> 844 Quark_Mom[index+i].setZ( Pz );
>> 845 if ( ScaleFactor != 1.0 ) {
>> 846 Quark_Mom[index+i].setX( SqrtScaleF * Quark_Mom[index+i].getX() );
>> 847 Quark_Mom[index+i].setY( SqrtScaleF * Quark_Mom[index+i].getY() );
>> 848 }
>> 849 }
>> 850 }
>> 851 for ( G4int iCase = 0; iCase < 2; iCase++ ) { // Loop over the two strings
865 G4ThreeVector tmp = Quark_Mom[iCase] + Q 852 G4ThreeVector tmp = Quark_Mom[iCase] + Quark_Mom[iCase+2];
866 G4LorentzVector Pstring( tmp, std::sqrt( 853 G4LorentzVector Pstring( tmp, std::sqrt( Quark_Mom[iCase].mag2() + MassQ2 ) +
867 std::sqrt( 854 std::sqrt( Quark_Mom[iCase+2].mag2() + MassQ2 ) );
>> 855 //AR-Jun2018 if ( common.RotateStrings ) Pstring *= common.RandomRotation;
868 // Add protection for rapidity = 0.5*ln 856 // Add protection for rapidity = 0.5*ln( (E+Pz)/(E-Pz) )
869 G4double Ystring = 0.0; 857 G4double Ystring = 0.0;
870 if ( Pstring.e() > 1.0e-30 ) { 858 if ( Pstring.e() > 1.0e-30 ) {
871 if ( Pstring.e() + Pstring.pz() < 1.0e 859 if ( Pstring.e() + Pstring.pz() < 1.0e-30 ) { // Very small numerator in the logarithm
872 Ystring = -1.0e30; // A very large 860 Ystring = -1.0e30; // A very large negative value (E ~= -Pz)
873 if ( Pstring.e() - Pstring.pz() < 1. 861 if ( Pstring.e() - Pstring.pz() < 1.0e-30 ) { // Very small denominator in the logarithm
874 Ystring = 1.0e30; // A very large 862 Ystring = 1.0e30; // A very large positive value (E ~= Pz)
875 } else { // Normal case 863 } else { // Normal case
876 Ystring = Pstring.rapidity(); 864 Ystring = Pstring.rapidity();
877 } 865 }
878 } 866 }
879 } 867 }
880 if ( iCase == 0 ) { // For the first st 868 if ( iCase == 0 ) { // For the first string
881 Pstring1 = Pstring; Ystring1 = Ystring 869 Pstring1 = Pstring; Ystring1 = Ystring;
882 } else { // For the second s 870 } else { // For the second string
883 Pstring2 = Pstring; Ystring2 = Ystring 871 Pstring2 = Pstring; Ystring2 = Ystring;
884 } 872 }
885 } 873 }
886 if ( Ystring1 > Ystring2 ) { 874 if ( Ystring1 > Ystring2 ) {
887 common.Pprojectile = Pstring1; common.P 875 common.Pprojectile = Pstring1; common.Ptarget = Pstring2;
888 QuarkOrder[0] = 0; QuarkOrder[1] = 1; <<
889 } else { 876 } else {
890 common.Pprojectile = Pstring2; common.P 877 common.Pprojectile = Pstring2; common.Ptarget = Pstring1;
891 QuarkOrder[0] = 1; QuarkOrder[1] = 0; <<
892 } 878 }
893 879
894 if ( common.RotateStrings ) { << 880 if ( common.RotateStrings ) { //AR-Jun2018
895 common.Pprojectile *= common.RandomRotat 881 common.Pprojectile *= common.RandomRotation;
896 common.Ptarget *= common.RandomRotat 882 common.Ptarget *= common.RandomRotation;
897 } 883 }
898 884
899 common.Pprojectile.transform( common.toLab 885 common.Pprojectile.transform( common.toLab );
900 common.Ptarget.transform( common.toLab ); 886 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_M <<
905 if ( common.RotateStrings ) Quark_4Mom[i <<
906 Quark_4Mom[i].transform( common.toLab ); <<
907 } <<
908 <<
909 projectile->Splitting(); <<
910 projectile->GetNextAntiParton()->Set4Momen <<
911 projectile->GetNextParton()->Set4Momentum( <<
912 <<
913 target->Splitting(); <<
914 target->GetNextParton()->Set4Momentum( Qua <<
915 target->GetNextAntiParton()->Set4Momentum( <<
916 887
917 // Calculation of the creation time 888 // Calculation of the creation time
918 // Creation time and position of target nu 889 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
919 projectile->SetTimeOfCreation( target->Get 890 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
920 projectile->SetPosition( target->GetPositi 891 projectile->SetPosition( target->GetPosition() );
921 projectile->Set4Momentum( common.Pprojecti 892 projectile->Set4Momentum( common.Pprojectile );
922 target->Set4Momentum( common.Ptarget ); 893 target->Set4Momentum( common.Ptarget );
923 <<
924 projectile->IncrementCollisionCount( 1 ); 894 projectile->IncrementCollisionCount( 1 );
925 target->IncrementCollisionCount( 1 ); 895 target->IncrementCollisionCount( 1 );
926 896
927 return 0; // Completed successfully: noth 897 return 0; // Completed successfully: nothing else to be done
928 } // End of if ( CandidatsN != 0 ) 898 } // End of if ( CandidatsN != 0 )
929 899
930 return 1; // Successfully ended, but the wo 900 return 1; // Successfully ended, but the work is not over
931 } 901 }
932 902
933 <<
934 //-------------------------------------------- 903 //-----------------------------------------------------------------------
935 904
936 G4bool G4FTFAnnihilation:: 905 G4bool G4FTFAnnihilation::
937 Create1QuarkAntiQuarkString( G4VSplitableHadro 906 Create1QuarkAntiQuarkString( G4VSplitableHadron* projectile,
938 G4VSplitableHadro 907 G4VSplitableHadron* target,
939 G4FTFParameters* 908 G4FTFParameters* theParameters,
940 G4FTFAnnihilation 909 G4FTFAnnihilation::CommonVariables& common ) const {
941 // Simulation of anti-quark - quark string c 910 // Simulation of anti-quark - quark string creation
942 911
943 #ifdef debugFTFannih 912 #ifdef debugFTFannih
944 G4cout << "Process d, only 1 quark - anti-qu 913 G4cout << "Process d, only 1 quark - anti-quark string" << G4endl;
945 #endif 914 #endif
946 915
947 // Determine the set of candidates anti-quar 916 // Determine the set of candidates anti-quark - quark pairs that do not annihilate.
948 // Here "0", "1", "2" means, respectively, " 917 // Here "0", "1", "2" means, respectively, "first", "second" and "third" constituent
949 // of the (anti-baryon) projectile or (nucle 918 // of the (anti-baryon) projectile or (nucleon) target.
950 G4int CandidatsN = 0, CandAQ[36], CandQ[36]; 919 G4int CandidatsN = 0, CandAQ[36], CandQ[36];
951 G4int LeftAQ = 0, LeftQ = 0; 920 G4int LeftAQ = 0, LeftQ = 0;
952 for ( G4int iAQ1 = 0; iAQ1 < 3; ++iAQ1 ) { << 921 for ( G4int iAQ1 = 0; iAQ1 < 3; iAQ1++ ) {
953 for ( G4int iAQ2 = 0; iAQ2 < 3; ++iAQ2 ) { << 922 for ( G4int iAQ2 = 0; iAQ2 < 3; iAQ2++ ) {
954 if ( iAQ1 != iAQ2 ) { 923 if ( iAQ1 != iAQ2 ) {
955 for ( G4int iQ1 = 0; iQ1 < 3; ++iQ1 ) << 924 for ( G4int iQ1 = 0; iQ1 < 3; iQ1++ ) {
956 for ( G4int iQ2 = 0; iQ2 < 3; ++iQ2 << 925 for ( G4int iQ2 = 0; iQ2 < 3; iQ2++ ) {
957 if ( iQ1 != iQ2 ) { 926 if ( iQ1 != iQ2 ) {
958 if ( -common.AQ[iAQ1] == common. 927 if ( -common.AQ[iAQ1] == common.Q[iQ1] && -common.AQ[iAQ2] == common.Q[iQ2] ) {
959 if ( ( iAQ1 == 0 && i 928 if ( ( iAQ1 == 0 && iAQ2 == 1 ) || ( iAQ1 == 1 && iAQ2 == 0 ) ) {
960 CandAQ[CandidatsN] = 2; 929 CandAQ[CandidatsN] = 2;
961 } else if ( ( iAQ1 == 0 && i 930 } else if ( ( iAQ1 == 0 && iAQ2 == 2 ) || ( iAQ1 == 2 && iAQ2 == 0 ) ) {
962 CandAQ[CandidatsN] = 1; 931 CandAQ[CandidatsN] = 1;
963 } else if ( ( iAQ1 == 1 && i 932 } else if ( ( iAQ1 == 1 && iAQ2 == 2 ) || ( iAQ1 == 2 && iAQ2 == 1 ) ) {
964 CandAQ[CandidatsN] = 0; 933 CandAQ[CandidatsN] = 0;
965 } 934 }
966 if ( ( iQ1 == 0 && 935 if ( ( iQ1 == 0 && iQ2 == 1 ) || ( iQ1 == 1 && iQ2 == 0 ) ) {
967 CandQ[CandidatsN] = 2; 936 CandQ[CandidatsN] = 2;
968 } else if ( ( iQ1 == 0 && 937 } else if ( ( iQ1 == 0 && iQ2 == 2 ) || ( iQ1 == 2 && iQ2 == 0 ) ) {
969 CandQ[CandidatsN] = 1; 938 CandQ[CandidatsN] = 1;
970 } else if ( ( iQ1 == 1 && 939 } else if ( ( iQ1 == 1 && iQ2 == 2 ) || ( iQ1 == 2 && iQ2 == 1 ) ) {
971 CandQ[CandidatsN] = 0; 940 CandQ[CandidatsN] = 0;
972 } 941 }
973 ++CandidatsN; << 942 CandidatsN++;
974 } 943 }
975 } 944 }
976 } 945 }
977 } 946 }
978 } 947 }
979 } 948 }
980 } 949 }
981 950
982 if ( CandidatsN != 0 ) { 951 if ( CandidatsN != 0 ) {
983 G4int SampledCase = (G4int)G4RandFlat::sho << 952 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) );
984 LeftAQ = common.AQ[ CandAQ[SampledCase] ]; 953 LeftAQ = common.AQ[ CandAQ[SampledCase] ];
985 LeftQ = common.Q[ CandQ[SampledCase] ]; 954 LeftQ = common.Q[ CandQ[SampledCase] ];
986 955
987 // Set the string properties 956 // Set the string properties
>> 957 projectile->SplitUp();
988 projectile->SetFirstParton( LeftQ ); 958 projectile->SetFirstParton( LeftQ );
989 projectile->SetSecondParton( LeftAQ ); 959 projectile->SetSecondParton( LeftAQ );
990 projectile->SetStatus( 0 ); 960 projectile->SetStatus( 0 );
991 G4int aAQ = std::abs( LeftAQ ), aQ = std:: 961 G4int aAQ = std::abs( LeftAQ ), aQ = std::abs( LeftQ );
992 G4int NewCode = 0; 962 G4int NewCode = 0;
993 G4double aKsi = G4UniformRand(); 963 G4double aKsi = G4UniformRand();
994 // The string can have the quantum number 964 // 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 965 // of the PDG code is, respectively, 1 and 3). For simplicity only scalar is considered here.
996 if ( aAQ == aQ ) { 966 if ( aAQ == aQ ) {
997 if ( aAQ != 3 ) { 967 if ( aAQ != 3 ) {
998 NewCode = 111; // Pi0-meson 968 NewCode = 111; // Pi0-meson
999 if ( aKsi < 0.5 ) { 969 if ( aKsi < 0.5 ) {
1000 NewCode = 221; // Eta -meson 970 NewCode = 221; // Eta -meson
1001 if ( aKsi < 0.25 ) { 971 if ( aKsi < 0.25 ) {
1002 NewCode = 331; // Eta'-meson 972 NewCode = 331; // Eta'-meson
1003 } 973 }
1004 } 974 }
1005 } else { 975 } else {
1006 NewCode = 221; // Eta -meson 976 NewCode = 221; // Eta -meson
1007 if ( aKsi < 0.5 ) { 977 if ( aKsi < 0.5 ) {
1008 NewCode = 331; // Eta'-meson 978 NewCode = 331; // Eta'-meson
1009 } 979 }
1010 } 980 }
1011 } else { 981 } else {
1012 if ( aAQ > aQ ) { 982 if ( aAQ > aQ ) {
1013 NewCode = aAQ*100 + aQ*10 + 1; NewCod 983 NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ;
1014 } else { 984 } else {
1015 NewCode = aQ*100 + aAQ*10 + 1; NewCod 985 NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ;
1016 } 986 }
1017 } 987 }
1018 988
1019 G4ParticleDefinition* TestParticle = G4Pa 989 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode );
1020 if ( ! TestParticle ) return false; 990 if ( ! TestParticle ) return false;
1021 projectile->SetDefinition( TestParticle ) 991 projectile->SetDefinition( TestParticle );
1022 theParameters->SetProjMinDiffMass( 0.5 ); 992 theParameters->SetProjMinDiffMass( 0.5 );
1023 theParameters->SetProjMinNonDiffMass( 0.5 993 theParameters->SetProjMinNonDiffMass( 0.5 );
1024 994
1025 target->SetStatus( 4 ); // The target nu 995 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4
1026 common.Pprojectile.setPx( 0.0 ); 996 common.Pprojectile.setPx( 0.0 );
1027 common.Pprojectile.setPy( 0.0 ); 997 common.Pprojectile.setPy( 0.0 );
1028 common.Pprojectile.setPz( 0.0 ); 998 common.Pprojectile.setPz( 0.0 );
1029 common.Pprojectile.setE( common.SqrtS ); 999 common.Pprojectile.setE( common.SqrtS );
1030 <<
1031 common.Pprojectile.transform( common.toLa 1000 common.Pprojectile.transform( common.toLab );
1032 1001
1033 G4LorentzVector Pquark = G4LorentzVector << 1002 G4LorentzVector Pquark = G4LorentzVector( 0.0, 0.0, common.SqrtS/2.0, common.SqrtS/2.0 );
1034 G4LorentzVector Paquark = G4LorentzVector << 1003 G4LorentzVector Paquark = G4LorentzVector( 0.0, 0.0, -common.SqrtS/2.0, common.SqrtS/2.0 );
1035 <<
1036 if ( common.RotateStrings ) { 1004 if ( common.RotateStrings ) {
1037 Pquark *= common.RandomRotation; Paquar 1005 Pquark *= common.RandomRotation; Paquark *= common.RandomRotation;
1038 } 1006 }
1039 Pquark.transform(common.toLab); projecti 1007 Pquark.transform(common.toLab); projectile->GetNextParton()->Set4Momentum(Pquark);
1040 Paquark.transform(common.toLab); projecti 1008 Paquark.transform(common.toLab); projectile->GetNextAntiParton()->Set4Momentum(Paquark);
1041 1009
1042 projectile->Splitting(); <<
1043 <<
1044 // Calculation of the creation time 1010 // Calculation of the creation time
1045 // Creation time and position of target n 1011 // Creation time and position of target nucleon were determined in ReggeonCascade() of G4FTFModel
1046 projectile->SetTimeOfCreation( target->Ge 1012 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1047 projectile->SetPosition( target->GetPosit 1013 projectile->SetPosition( target->GetPosition() );
1048 projectile->Set4Momentum( common.Pproject 1014 projectile->Set4Momentum( common.Pprojectile );
1049 <<
1050 projectile->IncrementCollisionCount( 1 ); 1015 projectile->IncrementCollisionCount( 1 );
1051 target->IncrementCollisionCount( 1 ); 1016 target->IncrementCollisionCount( 1 );
1052 1017
1053 return true; 1018 return true;
1054 } // End of if ( CandidatsN != 0 ) 1019 } // End of if ( CandidatsN != 0 )
1055 1020
1056 return true; 1021 return true;
1057 } 1022 }
1058 1023
1059 1024
1060 //=========================================== 1025 //============================================================================
1061 1026
1062 G4double G4FTFAnnihilation::ChooseX( G4double 1027 G4double G4FTFAnnihilation::ChooseX( G4double /* Alpha */, G4double /* Beta */ ) const {
1063 // If for sampling Xs other values of Alfa 1028 // If for sampling Xs other values of Alfa and Beta instead of 0.5 will be
1064 // chosen the method will be implemented 1029 // chosen the method will be implemented
1065 //G4double tmp = Alpha*Beta; 1030 //G4double tmp = Alpha*Beta;
1066 //tmp *= 1.0; 1031 //tmp *= 1.0;
1067 return 0.5; 1032 return 0.5;
1068 } 1033 }
1069 1034
1070 1035
1071 //=========================================== 1036 //============================================================================
1072 1037
1073 G4ThreeVector G4FTFAnnihilation::GaussianPt( 1038 G4ThreeVector G4FTFAnnihilation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1074 // @@ this method is used in FTFModel as w 1039 // @@ this method is used in FTFModel as well. Should go somewhere common!
1075 G4double Pt2 = 0.0; 1040 G4double Pt2 = 0.0;
1076 if ( AveragePt2 <= 0.0 ) { 1041 if ( AveragePt2 <= 0.0 ) {
1077 Pt2 = 0.0; 1042 Pt2 = 0.0;
1078 } else { 1043 } else {
1079 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor 1044 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1080 ( G4E 1045 ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
1081 } 1046 }
1082 G4double Pt = std::sqrt( Pt2 ); 1047 G4double Pt = std::sqrt( Pt2 );
1083 G4double phi = G4UniformRand() * twopi; 1048 G4double phi = G4UniformRand() * twopi;
1084 return G4ThreeVector ( Pt*std::cos( phi ), 1049 return G4ThreeVector ( Pt*std::cos( phi ), Pt*std::sin( phi ), 0.0 );
1085 } 1050 }
1086 1051
1087 1052
1088 //=========================================== 1053 //============================================================================
1089 1054
1090 void G4FTFAnnihilation::UnpackBaryon( G4int I 1055 void G4FTFAnnihilation::UnpackBaryon( G4int IdPDG, G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1091 G4int AbsId = std::abs( IdPDG ); 1056 G4int AbsId = std::abs( IdPDG );
1092 Q1 = AbsId / 1000; 1057 Q1 = AbsId / 1000;
1093 Q2 = ( AbsId % 1000 ) / 100; 1058 Q2 = ( AbsId % 1000 ) / 100;
1094 Q3 = ( AbsId % 100 ) / 10; 1059 Q3 = ( AbsId % 100 ) / 10;
1095 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = 1060 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = -Q3; } // Anti-baryon
1096 return; 1061 return;
1097 } 1062 }
1098 1063
1099 1064
1100 //=========================================== 1065 //============================================================================
1101 1066
1102 G4FTFAnnihilation::G4FTFAnnihilation( const G 1067 G4FTFAnnihilation::G4FTFAnnihilation( const G4FTFAnnihilation& ) {
1103 throw G4HadronicException( __FILE__, __LINE 1068 throw G4HadronicException( __FILE__, __LINE__,
1104 "G4FTFAnnihilati 1069 "G4FTFAnnihilation copy constructor not meant to be called" );
1105 } 1070 }
1106 1071
1107 1072
1108 //=========================================== 1073 //============================================================================
1109 1074
1110 const G4FTFAnnihilation & G4FTFAnnihilation:: 1075 const G4FTFAnnihilation & G4FTFAnnihilation::operator=( const G4FTFAnnihilation& ) {
1111 throw G4HadronicException( __FILE__, __LINE 1076 throw G4HadronicException( __FILE__, __LINE__,
1112 "G4FTFAnnihilati 1077 "G4FTFAnnihilation = operator not meant to be called" );
1113 } 1078 }
1114 1079
1115 1080
1116 //=========================================== 1081 //============================================================================
1117 1082
1118 G4bool G4FTFAnnihilation::operator==( const G 1083 G4bool G4FTFAnnihilation::operator==( const G4FTFAnnihilation& ) const {
1119 throw G4HadronicException( __FILE__, __LINE 1084 throw G4HadronicException( __FILE__, __LINE__,
1120 "G4FTFAnnihilati 1085 "G4FTFAnnihilation == operator not meant to be called" );
1121 } 1086 }
1122 1087
1123 1088
1124 //=========================================== 1089 //============================================================================
1125 1090
1126 G4bool G4FTFAnnihilation::operator!=( const G 1091 G4bool G4FTFAnnihilation::operator!=( const G4FTFAnnihilation& ) const {
1127 throw G4HadronicException( __FILE__, __LINE 1092 throw G4HadronicException( __FILE__, __LINE__,
1128 "G4DiffractiveEx 1093 "G4DiffractiveExcitation != operator not meant to be called" );
1129 } 1094 }
>> 1095
1130 1096