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