| 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 // ---------------- G4DiffractiveExcitati 32 // ---------------- G4DiffractiveExcitation --------------
33 // by Gunter Folger, October 1998. 33 // by Gunter Folger, October 1998.
34 // diffractive Excitation used by strin 34 // diffractive Excitation used by strings models
35 // Take a projectile and a target 35 // Take a projectile and a target
36 // excite the projectile and targe 36 // excite the projectile and target
37 // Essential changed by V. Uzhinsky in Novemb 37 // Essential changed by V. Uzhinsky in November - December 2006
38 // in order to put it in a correspondence wit 38 // in order to put it in a correspondence with original FRITIOF
39 // model. Variant of FRITIOF with nucleon de- 39 // model. Variant of FRITIOF with nucleon de-excitation is implemented.
40 // Other changes by V.Uzhinsky in May 2007 we 40 // Other changes by V.Uzhinsky in May 2007 were introduced to fit
41 // meson-nucleon interactions. Additional cha 41 // meson-nucleon interactions. Additional changes by V. Uzhinsky
42 // were introduced in December 2006. They tre 42 // were introduced in December 2006. They treat diffraction dissociation
43 // processes more exactly. 43 // processes more exactly.
44 // Correct treatment of the diffraction disso 44 // Correct treatment of the diffraction dissociation - 2012, V. Uzhinsky
45 // Mass distributions for resonances and uu-d 45 // Mass distributions for resonances and uu-diquark suppression in protons,
46 // and dd-diquarks suppression in neutrons we 46 // and dd-diquarks suppression in neutrons were introduced by V. Uzhinsky, 2014
47 // ------------------------------------------- 47 // ---------------------------------------------------------------------
48 48
49 #include "globals.hh" 49 #include "globals.hh"
50 #include "Randomize.hh" 50 #include "Randomize.hh"
51 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh" 52 #include "G4SystemOfUnits.hh"
53 53
54 #include "G4DiffractiveExcitation.hh" 54 #include "G4DiffractiveExcitation.hh"
55 #include "G4FTFParameters.hh" 55 #include "G4FTFParameters.hh"
56 #include "G4ElasticHNScattering.hh" 56 #include "G4ElasticHNScattering.hh"
57 57
58 #include "G4RotationMatrix.hh" 58 #include "G4RotationMatrix.hh"
59 #include "G4ParticleDefinition.hh" 59 #include "G4ParticleDefinition.hh"
60 #include "G4ParticleTable.hh" 60 #include "G4ParticleTable.hh"
61 #include "G4SampleResonance.hh" 61 #include "G4SampleResonance.hh"
62 #include "G4VSplitableHadron.hh" 62 #include "G4VSplitableHadron.hh"
63 #include "G4ExcitedString.hh" 63 #include "G4ExcitedString.hh"
64 #include "G4Neutron.hh" 64 #include "G4Neutron.hh"
65 65
66 #include "G4Exp.hh" 66 #include "G4Exp.hh"
67 #include "G4Log.hh" 67 #include "G4Log.hh"
68 #include "G4Pow.hh" 68 #include "G4Pow.hh"
69 69
70 //#include "G4ios.hh" 70 //#include "G4ios.hh"
71 71
72 //============================================ 72 //============================================================================
73 73
74 //#define debugFTFexcitation << 74 //#define debugFTFexictation
75 //#define debug_heavyHadrons 75 //#define debug_heavyHadrons
76 76
77 //============================================ 77 //============================================================================
78 78
79 G4DiffractiveExcitation::G4DiffractiveExcitati 79 G4DiffractiveExcitation::G4DiffractiveExcitation() {}
80 80
81 81
82 //============================================ 82 //============================================================================
83 83
84 G4DiffractiveExcitation::~G4DiffractiveExcitat 84 G4DiffractiveExcitation::~G4DiffractiveExcitation() {}
85 85
86 86
87 //============================================ 87 //============================================================================
88 88
89 G4bool G4DiffractiveExcitation::ExciteParticip 89 G4bool G4DiffractiveExcitation::ExciteParticipants( G4VSplitableHadron* projectile,
90 90 G4VSplitableHadron* target,
91 91 G4FTFParameters* theParameters,
92 92 G4ElasticHNScattering* theElastic ) const {
93 93
94 #ifdef debugFTFexcitation << 94 #ifdef debugFTFexictation
95 G4cout << G4endl << "FTF ExciteParticipants 95 G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
96 #endif 96 #endif
97 97
98 CommonVariables common; 98 CommonVariables common;
99 99
100 // Projectile parameters 100 // Projectile parameters
101 common.Pprojectile = projectile->Get4Momentu 101 common.Pprojectile = projectile->Get4Momentum();
102 if ( common.Pprojectile.z() < 0.0 ) return f 102 if ( common.Pprojectile.z() < 0.0 ) return false;
103 common.ProjectilePDGcode = projectile->GetDe 103 common.ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
104 common.absProjectilePDGcode = std::abs( comm 104 common.absProjectilePDGcode = std::abs( common.ProjectilePDGcode );
105 common.M0projectile = projectile->GetDefinit 105 common.M0projectile = projectile->GetDefinition()->GetPDGMass(); //Uzhi Aug.2019 common.Pprojectile.mag();
106 G4double ProjectileRapidity = common.Pprojec 106 G4double ProjectileRapidity = common.Pprojectile.rapidity();
107 107
108 // Target parameters 108 // Target parameters
109 common.Ptarget = target->Get4Momentum(); 109 common.Ptarget = target->Get4Momentum();
110 common.TargetPDGcode = target->GetDefinition 110 common.TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
111 common.absTargetPDGcode = std::abs( common.T 111 common.absTargetPDGcode = std::abs( common.TargetPDGcode );
112 common.M0target = target->GetDefinition()->G 112 common.M0target = target->GetDefinition()->GetPDGMass(); //Uzhi Aug.2019 common.Ptarget.mag();
113 G4double TargetRapidity = common.Ptarget.rap 113 G4double TargetRapidity = common.Ptarget.rapidity();
114 114
115 // Kinematical properties of the interaction 115 // Kinematical properties of the interactions
116 G4LorentzVector Psum = common.Pprojectile + 116 G4LorentzVector Psum = common.Pprojectile + common.Ptarget; // 4-momentum in Lab.
117 common.S = Psum.mag2(); 117 common.S = Psum.mag2();
118 common.SqrtS = std::sqrt( common.S ); 118 common.SqrtS = std::sqrt( common.S );
119 119
120 // Check off-shellness of the participants 120 // Check off-shellness of the participants
121 G4bool toBePutOnMassShell = true; //Uzhi Au 121 G4bool toBePutOnMassShell = true; //Uzhi Aug.2019 false;
122 common.MminProjectile = common.BrW.GetMinimu 122 common.MminProjectile = common.BrW.GetMinimumMass( projectile->GetDefinition() );
123 /* Uzhi Aug.2019 123 /* Uzhi Aug.2019
124 if ( common.M0projectile < common.MminProjec 124 if ( common.M0projectile < common.MminProjectile ) {
125 toBePutOnMassShell = true; 125 toBePutOnMassShell = true;
126 common.M0projectile = common.BrW.SampleMas 126 common.M0projectile = common.BrW.SampleMass( projectile->GetDefinition(),
127 127 projectile->GetDefinition()->GetPDGMass()
128 128 + 5.0*projectile->GetDefinition()->GetPDGWidth() );
129 } 129 }
130 */ 130 */
131 common.M0projectile2 = common.M0projectile * 131 common.M0projectile2 = common.M0projectile * common.M0projectile;
132 common.ProjectileDiffStateMinMass = thePa 132 common.ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass();
133 common.ProjectileNonDiffStateMinMass = thePa 133 common.ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass();
134 if ( common.M0projectile > common.Projectile 134 if ( common.M0projectile > common.ProjectileDiffStateMinMass ) {
135 common.ProjectileDiffStateMinMass = com 135 common.ProjectileDiffStateMinMass = common.MminProjectile + 220.0*MeV; //Uzhi Aug.2019 common.M0projectile + 220.0*MeV;
136 common.ProjectileNonDiffStateMinMass = com 136 common.ProjectileNonDiffStateMinMass = common.MminProjectile + 220.0*MeV; //Uzhi Aug.2019 common.M0projectile + 220.0*MeV;
137 if ( common.absProjectilePDGcode > 3000 ) 137 if ( common.absProjectilePDGcode > 3000 ) { // Strange baryon
138 common.ProjectileDiffStateMinMass += 138 common.ProjectileDiffStateMinMass += 140.0*MeV;
139 common.ProjectileNonDiffStateMinMass += 139 common.ProjectileNonDiffStateMinMass += 140.0*MeV;
140 } 140 }
141 } 141 }
142 common.MminTarget = common.BrW.GetMinimumMas 142 common.MminTarget = common.BrW.GetMinimumMass( target->GetDefinition() );
143 /* Uzhi Aug.2019 143 /* Uzhi Aug.2019
144 if ( common.M0target < common.MminTarget ) { 144 if ( common.M0target < common.MminTarget ) {
145 toBePutOnMassShell = true; 145 toBePutOnMassShell = true;
146 common.M0target = common.BrW.SampleMass( t 146 common.M0target = common.BrW.SampleMass( target->GetDefinition(),
147 t 147 target->GetDefinition()->GetPDGMass()
148 + 148 + 5.0*target->GetDefinition()->GetPDGWidth() );
149 } 149 }
150 */ 150 */
151 common.M0target2 = common.M0target * common. 151 common.M0target2 = common.M0target * common.M0target;
152 common.TargetDiffStateMinMass = theParame 152 common.TargetDiffStateMinMass = theParameters->GetTarMinDiffMass();
153 common.TargetNonDiffStateMinMass = theParame 153 common.TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass();
154 if ( common.M0target > common.TargetDiffStat 154 if ( common.M0target > common.TargetDiffStateMinMass ) {
155 common.TargetDiffStateMinMass = common. 155 common.TargetDiffStateMinMass = common.MminTarget + 220.0*MeV; //Uzhi Aug.2019 common.M0target + 220.0*MeV;
156 common.TargetNonDiffStateMinMass = common. 156 common.TargetNonDiffStateMinMass = common.MminTarget + 220.0*MeV; //Uzhi Aug.2019 common.M0target + 220.0*MeV;
157 if ( common.absTargetPDGcode > 3000 ) { / 157 if ( common.absTargetPDGcode > 3000 ) { // Strange baryon
158 common.TargetDiffStateMinMass += 140. 158 common.TargetDiffStateMinMass += 140.0*MeV;
159 common.TargetNonDiffStateMinMass += 140. 159 common.TargetNonDiffStateMinMass += 140.0*MeV;
160 } 160 }
161 }; 161 };
162 #ifdef debugFTFexcitation << 162 #ifdef debugFTFexictation
163 G4cout << "Proj Targ PDGcodes " << common.Pr 163 G4cout << "Proj Targ PDGcodes " << common.ProjectilePDGcode << " " << common.TargetPDGcode << G4endl
164 << "Mprojectile Y " << common.Pproje 164 << "Mprojectile Y " << common.Pprojectile.mag() << " " << ProjectileRapidity << G4endl // Uzhi Aug.2019
165 << "M0projectile Y " << common.M0proj 165 << "M0projectile Y " << common.M0projectile << " " << ProjectileRapidity << G4endl;
166 G4cout << "Mtarget Y " << common.Ptarge 166 G4cout << "Mtarget Y " << common.Ptarget.mag() << " " << TargetRapidity << G4endl // Uzhi Aug.2019
167 << "M0target Y " << common.M0targ 167 << "M0target Y " << common.M0target << " " << TargetRapidity << G4endl;
168 G4cout << "Pproj " << common.Pprojectile << 168 G4cout << "Pproj " << common.Pprojectile << G4endl << "Ptarget " << common.Ptarget << G4endl;
169 #endif 169 #endif
170 170
171 // Transform momenta to cms and then rotate 171 // Transform momenta to cms and then rotate parallel to z axis;
172 common.toCms = G4LorentzRotation( -1 * Psum. 172 common.toCms = G4LorentzRotation( -1 * Psum.boostVector() );
173 G4LorentzVector Ptmp = common.toCms * common 173 G4LorentzVector Ptmp = common.toCms * common.Pprojectile;
174 if ( Ptmp.pz() <= 0.0 ) return false; // "S 174 if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
175 common.toCms.rotateZ( -1*Ptmp.phi() ); 175 common.toCms.rotateZ( -1*Ptmp.phi() );
176 common.toCms.rotateY( -1*Ptmp.theta() ); 176 common.toCms.rotateY( -1*Ptmp.theta() );
177 common.toLab = common.toCms.inverse(); 177 common.toLab = common.toCms.inverse();
178 common.Pprojectile.transform( common.toCms ) 178 common.Pprojectile.transform( common.toCms );
179 common.Ptarget.transform( common.toCms ); 179 common.Ptarget.transform( common.toCms );
180 180
181 G4double SumMasses = common.M0projectile + c 181 G4double SumMasses = common.M0projectile + common.M0target; // + 220.0*MeV;
182 #ifdef debugFTFexcitation << 182 #ifdef debugFTFexictation
183 G4cout << "SqrtS " << common.SqrtS << G4 183 G4cout << "SqrtS " << common.SqrtS << G4endl << "M0pr M0tr SumM " << common.M0projectile
184 << " " << common.M0target << " " << S 184 << " " << common.M0target << " " << SumMasses << G4endl;
185 #endif 185 #endif
186 if ( common.SqrtS < SumMasses ) return false 186 if ( common.SqrtS < SumMasses ) return false; // The model cannot work at low energy
187 187
188 common.PZcms2 = ( sqr( common.S ) + sqr( com 188 common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
189 - 2.0 * ( common.S * ( com 189 - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
190 + common.M0proje 190 + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
191 #ifdef debugFTFexcitation << 191 #ifdef debugFTFexictation
192 G4cout << "PZcms2 after toBePutOnMassShell " 192 G4cout << "PZcms2 after toBePutOnMassShell " << common.PZcms2 << G4endl;
193 #endif 193 #endif
194 if ( common.PZcms2 < 0.0 ) return false; // 194 if ( common.PZcms2 < 0.0 ) return false; // It can be in an interaction with off-shell nuclear nucleon
195 195
196 common.PZcms = std::sqrt( common.PZcms2 ); 196 common.PZcms = std::sqrt( common.PZcms2 );
197 if ( toBePutOnMassShell ) { 197 if ( toBePutOnMassShell ) {
198 if ( common.Pprojectile.z() > 0.0 ) { 198 if ( common.Pprojectile.z() > 0.0 ) {
199 common.Pprojectile.setPz( common.PZcms 199 common.Pprojectile.setPz( common.PZcms );
200 common.Ptarget.setPz( -common.PZcms 200 common.Ptarget.setPz( -common.PZcms );
201 } else { 201 } else {
202 common.Pprojectile.setPz( -common.PZcms 202 common.Pprojectile.setPz( -common.PZcms );
203 common.Ptarget.setPz( common.PZcms 203 common.Ptarget.setPz( common.PZcms );
204 } 204 }
205 common.Pprojectile.setE( std::sqrt( common 205 common.Pprojectile.setE( std::sqrt( common.M0projectile2
206 + comm 206 + common.Pprojectile.x() * common.Pprojectile.x()
207 + comm 207 + common.Pprojectile.y() * common.Pprojectile.y()
208 + comm 208 + common.PZcms2 ) );
209 common.Ptarget.setE( std::sqrt( common.M0t 209 common.Ptarget.setE( std::sqrt( common.M0target2
210 + common.P 210 + common.Ptarget.x() * common.Ptarget.x()
211 + common.P 211 + common.Ptarget.y() * common.Ptarget.y()
212 + common.P 212 + common.PZcms2 ) );
213 } 213 }
214 #ifdef debugFTFexcitation << 214 #ifdef debugFTFexictation
215 G4cout << "Start --------------------" << G4 215 G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << common.M0projectile
216 << " " << common.ProjectileDiffStateM 216 << " " << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
217 << G4endl 217 << G4endl
218 << "Targ M0 Mdif Mndif " << common.M0 218 << "Targ M0 Mdif Mndif " << common.M0target << " " << common.TargetDiffStateMinMass
219 << " " << common.TargetNonDiffStateMi 219 << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS << G4endl
220 << "Proj CMS " << common.Pprojectile 220 << "Proj CMS " << common.Pprojectile << G4endl << "Targ CMS " << common.Ptarget << G4endl;
221 #endif 221 #endif
222 222
223 // Check for possible quark exchange 223 // Check for possible quark exchange
224 ProjectileRapidity = common.Pprojectile.rapi 224 ProjectileRapidity = common.Pprojectile.rapidity();
225 TargetRapidity = common.Ptarget.rapidity(); 225 TargetRapidity = common.Ptarget.rapidity();
226 G4double QeNoExc = theParameters->GetProcPro 226 G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity - TargetRapidity );
227 G4double QeExc = theParameters->GetProcPro 227 G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity - TargetRapidity ) *
228 theParameters->GetProcPro 228 theParameters->GetProcProb( 4, ProjectileRapidity - TargetRapidity );
229 common.ProbProjectileDiffraction = 229 common.ProbProjectileDiffraction =
230 theParameters->GetProcProb( 2, ProjectileR 230 theParameters->GetProcProb( 2, ProjectileRapidity - TargetRapidity );
231 common.ProbTargetDiffraction = 231 common.ProbTargetDiffraction =
232 theParameters->GetProcProb( 3, ProjectileR 232 theParameters->GetProcProb( 3, ProjectileRapidity - TargetRapidity );
233 common.ProbOfDiffraction = common.ProbProjec 233 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
234 #ifdef debugFTFexcitation << 234 #ifdef debugFTFexictation
235 G4cout << "Proc Probs " << QeNoExc << " " << 235 G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
236 << common.ProbProjectileDiffraction < 236 << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
237 << "ProjectileRapidity " << Projectil 237 << "ProjectileRapidity " << ProjectileRapidity << G4endl;
238 #endif 238 #endif
239 239
240 if ( QeNoExc + QeExc + common.ProbProjectile 240 if ( QeNoExc + QeExc + common.ProbProjectileDiffraction + common.ProbTargetDiffraction > 1.0 ) {
241 QeNoExc = 1.0 - QeExc - common.ProbProject 241 QeNoExc = 1.0 - QeExc - common.ProbProjectileDiffraction - common.ProbTargetDiffraction;
242 } 242 }
243 if ( QeExc + QeNoExc != 0.0 ) { 243 if ( QeExc + QeNoExc != 0.0 ) {
244 common.ProbExc = QeExc / ( QeExc + QeNoExc 244 common.ProbExc = QeExc / ( QeExc + QeNoExc );
245 } 245 }
246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) { 246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
247 common.ProbProjectileDiffraction /= ( 1.0 247 common.ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
248 common.ProbTargetDiffraction /= ( 1.0 248 common.ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
249 } 249 }
250 #ifdef debugFTFexcitation << 250 #ifdef debugFTFexictation
251 G4cout << "Proc Probs " << QeNoExc << " " << 251 G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
252 << common.ProbProjectileDiffraction < 252 << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
253 << "ProjectileRapidity " << Projectil 253 << "ProjectileRapidity " << ProjectileRapidity << G4endl;
254 #endif 254 #endif
255 255
256 // Try out charge exchange 256 // Try out charge exchange
257 G4int returnCode = 1; 257 G4int returnCode = 1;
258 if ( G4UniformRand() < QeExc + QeNoExc ) { 258 if ( G4UniformRand() < QeExc + QeNoExc ) {
259 returnCode = ExciteParticipants_doChargeEx 259 returnCode = ExciteParticipants_doChargeExchange( projectile, target, theParameters,
260 260 theElastic, common );
261 } 261 }
262 262
263 G4bool returnResult = false; 263 G4bool returnResult = false;
264 if ( returnCode == 0 ) { 264 if ( returnCode == 0 ) {
265 returnResult = true; // Successfully ende 265 returnResult = true; // Successfully ended: no need of extra work
266 } else if ( returnCode == 1 ) { 266 } else if ( returnCode == 1 ) {
267 267
268 common.ProbOfDiffraction = common.ProbProj 268 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
269 #ifdef debugFTFexcitation << 269 #ifdef debugFTFexictation
270 G4cout << "Excitation -------------------- 270 G4cout << "Excitation --------------------" << G4endl
271 << "Proj M0 MdMin MndMin " << commo 271 << "Proj M0 MdMin MndMin " << common.M0projectile << " "
272 << common.ProjectileDiffStateMinMas 272 << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
273 << G4endl 273 << G4endl
274 << "Targ M0 MdMin MndMin " << commo 274 << "Targ M0 MdMin MndMin " << common.M0target << " " << common.TargetDiffStateMinMass
275 << " " << common.TargetNonDiffState 275 << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS
276 << G4endl 276 << G4endl
277 << "Prob: ProjDiff TargDiff + Sum " 277 << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
278 << common.ProbTargetDiffraction << 278 << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
279 #endif 279 #endif
280 if ( common.ProbOfDiffraction != 0.0 ) { 280 if ( common.ProbOfDiffraction != 0.0 ) {
281 common.ProbProjectileDiffraction /= comm 281 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
282 } else { 282 } else {
283 common.ProbProjectileDiffraction = 0.0; 283 common.ProbProjectileDiffraction = 0.0;
284 } 284 }
285 #ifdef debugFTFexcitation << 285 #ifdef debugFTFexictation
286 G4cout << "Prob: ProjDiff TargDiff + Sum " 286 G4cout << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
287 << common.ProbTargetDiffraction << 287 << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
288 #endif 288 #endif
289 common.ProjectileDiffStateMinMass2 = sq 289 common.ProjectileDiffStateMinMass2 = sqr( common.ProjectileDiffStateMinMass );
290 common.ProjectileNonDiffStateMinMass2 = sq 290 common.ProjectileNonDiffStateMinMass2 = sqr( common.ProjectileNonDiffStateMinMass );
291 common.TargetDiffStateMinMass2 = sq 291 common.TargetDiffStateMinMass2 = sqr( common.TargetDiffStateMinMass );
292 common.TargetNonDiffStateMinMass2 = sq 292 common.TargetNonDiffStateMinMass2 = sqr( common.TargetNonDiffStateMinMass );
293 // Choose between diffraction and non-diff 293 // Choose between diffraction and non-diffraction process
294 if ( G4UniformRand() < common.ProbOfDiffra 294 if ( G4UniformRand() < common.ProbOfDiffraction ) {
295 295
296 returnResult = ExciteParticipants_doDiff 296 returnResult = ExciteParticipants_doDiffraction( projectile, target, theParameters, common );
297 297
298 } else { 298 } else {
299 299
300 returnResult = ExciteParticipants_doNonD 300 returnResult = ExciteParticipants_doNonDiffraction( projectile, target, theParameters, common );
301 301
302 } 302 }
303 if ( returnResult ) { 303 if ( returnResult ) {
304 common.Pprojectile += common.Qmomentum; 304 common.Pprojectile += common.Qmomentum;
305 common.Ptarget -= common.Qmomentum; 305 common.Ptarget -= common.Qmomentum;
306 // Transform back and update SplitableHa 306 // Transform back and update SplitableHadron Participant.
307 common.Pprojectile.transform( common.toL 307 common.Pprojectile.transform( common.toLab );
308 common.Ptarget.transform( common.toLab ) 308 common.Ptarget.transform( common.toLab );
309 #ifdef debugFTFexcitation << 309 #ifdef debugFTFexictation
310 G4cout << "Mproj " << common.Pprojectile 310 G4cout << "Mproj " << common.Pprojectile.mag() << G4endl << "Mtarg " << common.Ptarget.mag()
311 << G4endl; 311 << G4endl;
312 #endif 312 #endif
313 projectile->Set4Momentum( common.Pprojec 313 projectile->Set4Momentum( common.Pprojectile );
314 target->Set4Momentum( common.Ptarget ); 314 target->Set4Momentum( common.Ptarget );
315 projectile->IncrementCollisionCount( 1 ) 315 projectile->IncrementCollisionCount( 1 );
316 target->IncrementCollisionCount( 1 ); 316 target->IncrementCollisionCount( 1 );
317 } 317 }
318 } 318 }
319 319
320 return returnResult; 320 return returnResult;
321 } 321 }
322 322
323 //-------------------------------------------- 323 //-----------------------------------------------------------------------------
324 324
325 G4int G4DiffractiveExcitation:: 325 G4int G4DiffractiveExcitation::
326 ExciteParticipants_doChargeExchange( G4VSplita 326 ExciteParticipants_doChargeExchange( G4VSplitableHadron* projectile,
327 G4VSplita 327 G4VSplitableHadron* target,
328 G4FTFPara 328 G4FTFParameters* theParameters,
329 G4Elastic 329 G4ElasticHNScattering* theElastic,
330 G4Diffrac 330 G4DiffractiveExcitation::CommonVariables& common ) const {
331 // First of the three utility methods used o 331 // First of the three utility methods used only by ExciteParticipants:
332 // it does the sampling for the charge-excha 332 // it does the sampling for the charge-exchange case.
333 // This method returns an integer code - ins 333 // This method returns an integer code - instead of a boolean, with the following meaning:
334 // "0" : successfully ended and nothing el 334 // "0" : successfully ended and nothing else needs to be done;
335 // "1" : successfully completed, but the w 335 // "1" : successfully completed, but the work needs to be continued;
336 // "99" : unsuccessfully ended, nothing els 336 // "99" : unsuccessfully ended, nothing else can be done.
337 G4int returnCode = 99; 337 G4int returnCode = 99;
338 338
339 G4double DeltaProbAtQuarkExchange = theParam 339 G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
340 G4ParticleDefinition* TestParticle = 0; 340 G4ParticleDefinition* TestParticle = 0;
341 G4double MtestPr = 0.0, MtestTr = 0.0; 341 G4double MtestPr = 0.0, MtestTr = 0.0;
342 342
343 #ifdef debugFTFexcitation << 343 #ifdef debugFTFexictation
344 G4cout << "Q exchange ---------------------- 344 G4cout << "Q exchange --------------------------" << G4endl;
345 #endif 345 #endif
346 346
347 // The target hadron is always a nucleon (i. 347 // The target hadron is always a nucleon (i.e. either a proton or a neutron,
348 // never an antinucleon), therefore only a q 348 // never an antinucleon), therefore only a quark (not an antiquark) can be
349 // exchanged between the projectile hadron a 349 // exchanged between the projectile hadron and the target hadron (otherwise
350 // we could get a quark-quark-antiquark syst 350 // we could get a quark-quark-antiquark system which cannot be a bound state).
351 // This implies that any projectile meson or 351 // This implies that any projectile meson or anti-meson - given that it has
352 // a constituent quark in all cases - can ha 352 // a constituent quark in all cases - can have charge exchange with a target
353 // hadron. Instead, any projectile anti-bary 353 // hadron. Instead, any projectile anti-baryon can never have charge exchange
354 // with a target hadron (because it has only 354 // with a target hadron (because it has only constituent anti-quarks);
355 // projectile baryons, instead can have char 355 // projectile baryons, instead can have charge exchange with a target hadron.
356 356
357 G4int NewProjCode = 0, NewTargCode = 0, Proj 357 G4int NewProjCode = 0, NewTargCode = 0, ProjQ1 = 0, ProjQ2 = 0, ProjQ3 = 0;
358 // Projectile unpacking 358 // Projectile unpacking
359 if ( common.absProjectilePDGcode < 1000 ) { 359 if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
360 UnpackMeson( common.ProjectilePDGcode, Pr 360 UnpackMeson( common.ProjectilePDGcode, ProjQ1, ProjQ2 );
361 } else { 361 } else { // projectile is baryon
362 UnpackBaryon( common.ProjectilePDGcode, Pr 362 UnpackBaryon( common.ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
363 } 363 }
364 // Target unpacking 364 // Target unpacking
365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0; 365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;
366 UnpackBaryon( common.TargetPDGcode, TargQ1, 366 UnpackBaryon( common.TargetPDGcode, TargQ1, TargQ2, TargQ3 );
367 #ifdef debugFTFexcitation << 367 #ifdef debugFTFexictation
368 G4cout << "Proj Quarks " << ProjQ1 << " " << 368 G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
369 << "Targ Quarks " << TargQ1 << " " << 369 << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
370 #endif 370 #endif
371 371
372 // Sampling of exchanged quarks 372 // Sampling of exchanged quarks
373 G4int ProjExchangeQ = 0, TargExchangeQ = 0; 373 G4int ProjExchangeQ = 0, TargExchangeQ = 0;
374 if ( common.absProjectilePDGcode < 1000 ) { 374 if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
375 375
376 G4bool isProjQ1Quark = false; 376 G4bool isProjQ1Quark = false;
377 ProjExchangeQ = ProjQ2; 377 ProjExchangeQ = ProjQ2;
378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark 378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark
379 isProjQ1Quark = true; 379 isProjQ1Quark = true;
380 ProjExchangeQ = ProjQ1; 380 ProjExchangeQ = ProjQ1;
381 } 381 }
382 // Exchange of non-identical quarks is all 382 // Exchange of non-identical quarks is allowed
383 G4int NpossibleStates = 0; 383 G4int NpossibleStates = 0;
384 if ( ProjExchangeQ != TargQ1 ) NpossibleSt 384 if ( ProjExchangeQ != TargQ1 ) NpossibleStates++;
385 if ( ProjExchangeQ != TargQ2 ) NpossibleSt 385 if ( ProjExchangeQ != TargQ2 ) NpossibleStates++;
386 if ( ProjExchangeQ != TargQ3 ) NpossibleSt 386 if ( ProjExchangeQ != TargQ3 ) NpossibleStates++;
387 G4int Nsampled = (G4int)G4RandFlat::shootI << 387 G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
388 NpossibleStates = 0; 388 NpossibleStates = 0;
389 if ( ProjExchangeQ != TargQ1 ) { 389 if ( ProjExchangeQ != TargQ1 ) {
390 if ( ++NpossibleStates == Nsampled ) { 390 if ( ++NpossibleStates == Nsampled ) {
391 TargExchangeQ = TargQ1; TargQ1 = ProjE 391 TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ;
392 isProjQ1Quark ? ProjQ1 = TargExchangeQ 392 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
393 } 393 }
394 } 394 }
395 if ( ProjExchangeQ != TargQ2 ) { 395 if ( ProjExchangeQ != TargQ2 ) {
396 if ( ++NpossibleStates == Nsampled ) { 396 if ( ++NpossibleStates == Nsampled ) {
397 TargExchangeQ = TargQ2; TargQ2 = ProjE 397 TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ;
398 isProjQ1Quark ? ProjQ1 = TargExchangeQ 398 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
399 } 399 }
400 } 400 }
401 if ( ProjExchangeQ != TargQ3 ) { 401 if ( ProjExchangeQ != TargQ3 ) {
402 if ( ++NpossibleStates == Nsampled ) { 402 if ( ++NpossibleStates == Nsampled ) {
403 TargExchangeQ = TargQ3; TargQ3 = ProjE 403 TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ;
404 isProjQ1Quark ? ProjQ1 = TargExchangeQ 404 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
405 } 405 }
406 } 406 }
407 #ifdef debugFTFexcitation << 407 #ifdef debugFTFexictation
408 G4cout << "Exchanged Qs in Pr Tr " << Proj 408 G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
409 #endif 409 #endif
410 410
411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ 411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ2 = std::abs( ProjQ2 );
412 G4bool ProjExcited = false; 412 G4bool ProjExcited = false;
413 const G4int maxNumberOfAttempts = 50; 413 const G4int maxNumberOfAttempts = 50;
414 G4int attempts = 0; 414 G4int attempts = 0;
415 while ( attempts++ < maxNumberOfAttempts ) 415 while ( attempts++ < maxNumberOfAttempts ) { /* Loop checking, 10.08.2015, A.Ribon */
416 416
417 // Determination of a new projectile ID 417 // Determination of a new projectile ID which satisfies energy-momentum conservation
418 G4double ProbSpin0 = 0.5; <<
419 G4double Ksi = G4UniformRand(); 418 G4double Ksi = G4UniformRand();
420 if ( aProjQ1 == aProjQ2 ) { 419 if ( aProjQ1 == aProjQ2 ) {
421 if ( G4UniformRand() < ProbSpin0 ) { << 420 if ( aProjQ1 < 3 ) {
422 if ( aProjQ1 < 3 ) { << 421 NewProjCode = 111; // Pi0-meson
423 NewProjCode = 111; << 422 if ( Ksi < 0.5 ) {
424 if ( Ksi < 0.5 ) { << 423 NewProjCode = 221; // Eta-meson
425 NewProjCode = 221; << 424 if ( Ksi < 0.25 ) {
426 if ( Ksi < 0.25 ) { << 425 NewProjCode = 331; // Eta'-meson
427 NewProjCode = 331; << 426 }
428 } << 427 }
429 } << 428 } else if ( aProjQ1 == 3 ) {
430 } else if ( aProjQ1 == 3 ) { << 429 NewProjCode = 221; // Eta-meson
431 NewProjCode = 221; << 430 if ( Ksi < 0.5 ) {
432 if ( Ksi < 0.5 ) { << 431 NewProjCode = 331; // Eta'-meson
433 NewProjCode = 331; << 432 }
434 } << 433 } else if ( aProjQ1 == 4 ) {
435 } else if ( aProjQ1 == 4 ) { << 434 NewProjCode = 441; // Eta_c
436 NewProjCode = 441; // eta << 435 } else if ( aProjQ1 == 5 ) {
437 } else if ( aProjQ1 == 5 ) { << 436 NewProjCode = 553; // Upsilon
438 NewProjCode = 551; // eta << 437 }
439 } <<
440 } else { <<
441 if ( aProjQ1 < 3 ) { <<
442 NewProjCode = 113; <<
443 if ( Ksi < 0.5 ) { <<
444 NewProjCode = 223; <<
445 } <<
446 } else if ( aProjQ1 == 3 ) { <<
447 NewProjCode = 333; <<
448 } else if ( aProjQ1 == 4 ) { <<
449 NewProjCode = 443; // J/p <<
450 } else if ( aProjQ1 == 5 ) { <<
451 NewProjCode = 553; // Ups <<
452 } <<
453 } <<
454 } else { 438 } else {
455 if ( aProjQ1 > aProjQ2 ) { 439 if ( aProjQ1 > aProjQ2 ) {
456 NewProjCode = aProjQ1*100 + aProjQ2* 440 NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
457 } else { 441 } else {
458 NewProjCode = aProjQ2*100 + aProjQ1* 442 NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
459 } 443 }
460 } 444 }
461 #ifdef debugFTFexcitation << 445 #ifdef debugFTFexictation
462 G4cout << "NewProjCode " << NewProjCode 446 G4cout << "NewProjCode " << NewProjCode << G4endl;
463 #endif 447 #endif
464 // Decide (with 50% probability) whether 448 // Decide (with 50% probability) whether the projectile hadrons is excited,
465 // but not in the case of charmed and bo 449 // but not in the case of charmed and bottom hadrons (because in Geant4
466 // there are no excited charmed and bott 450 // there are no excited charmed and bottom states).
467 ProjExcited = false; 451 ProjExcited = false;
468 if ( aProjQ1 <= 3 && aProjQ2 <= 3 && 452 if ( aProjQ1 <= 3 && aProjQ2 <= 3 && G4UniformRand() < 0.5 ) {
469 NewProjCode += 2; // Excited meson (J 453 NewProjCode += 2; // Excited meson (J=1 -> 2*J+1=2+1=3, last digit of the PDG code)
470 ProjExcited = true; 454 ProjExcited = true;
471 } 455 }
472 456
473 // So far we have used the absolute valu 457 // So far we have used the absolute values of the PDG codes of the two constituent quarks:
474 // now look at their signed values to se 458 // now look at their signed values to set properly the signed of the meson's PDG code.
475 G4int value = ProjQ1, absValue = aProjQ1 459 G4int value = ProjQ1, absValue = aProjQ1, Qquarks = 0;
476 for ( G4int iQ = 0; iQ < 2; ++iQ ) { // 460 for ( G4int iQ = 0; iQ < 2; ++iQ ) { // 0 : first quark; 1 : second quark
477 if ( iQ == 1 ) { 461 if ( iQ == 1 ) {
478 value = ProjQ2; absValue = aProjQ2; 462 value = ProjQ2; absValue = aProjQ2;
479 } 463 }
480 if ( absValue == 2 || absValue == 4 464 if ( absValue == 2 || absValue == 4 ) { // u or c
481 Qquarks += 2*value/absValue; // u, 465 Qquarks += 2*value/absValue; // u, c : positively charged +2 (e/3 unit)
482 } else { 466 } else {
483 Qquarks -= value/absValue; // d, 467 Qquarks -= value/absValue; // d, s, b : negatively charged -1 (e/3 unit)
484 } 468 }
485 } 469 }
486 // If Qquarks is positive, the sign of N 470 // If Qquarks is positive, the sign of NewProjCode is fine.
487 // If Qquarks is negative, then the sign 471 // If Qquarks is negative, then the sign of NewProjCode needs to be reversed.
488 // If Qquarks is zero, then we need to d 472 // If Qquarks is zero, then we need to distinguish between 3 cases:
489 // 1. If aProjQ1 and aProjQ2 are the sam 473 // 1. If aProjQ1 and aProjQ2 are the same, then the sign of NewProjCode is fine
490 // (because the antiparticle is the s 474 // (because the antiparticle is the same as the particle, e.g. eta, eta').
491 // 2. If aProjQ1 and aProjQ2 are not the 475 // 2. If aProjQ1 and aProjQ2 are not the same, given that Qquarks is zero,
492 // we have only two possibilities: 476 // we have only two possibilities:
493 // a. aProjQ1 and aProjQ2 are two dif 477 // a. aProjQ1 and aProjQ2 are two different down-type quarks, i.e.
494 // (s,d) or (b,d), or (b,s). In th 478 // (s,d) or (b,d), or (b,s). In this case, the sign of NewProjCode
495 // is fine (because the heaviest o 479 // is fine (because the heaviest of the two down-type quarks has
496 // to be anti-quark belonging to t 480 // to be anti-quark belonging to the initial projectile, which
497 // implies a meson with positive P 481 // implies a meson with positive PDG code, e.g. B0 (bbar,d), Bs (bbar,s).
498 // b. aProjQ1 and aProjQ2 are two dif 482 // b. aProjQ1 and aProjQ2 are two different up-type quarks, i.e. (u,c).
499 // The heaviest of the two (c) has 483 // The heaviest of the two (c) has to be an anti-quark (cbar) left
500 // in the projectile, therefore th 484 // in the projectile, therefore the sign of NewProjCode needs to be
501 // reverse: 421 -> -421 anti_D0 (c 485 // reverse: 421 -> -421 anti_D0 (cbar,u)
502 if ( Qquarks < 0 || ( Qquarks == 0 && 486 if ( Qquarks < 0 || ( Qquarks == 0 && aProjQ1 != aProjQ2 && aProjQ1%2 == 0 ) ) {
503 NewProjCode *= -1; 487 NewProjCode *= -1;
504 } 488 }
505 #ifdef debugFTFexcitation << 489 #ifdef debugFTFexictation
506 G4cout << "NewProjCode +2 or 0 " << NewP 490 G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
507 G4cout<<"+++++++++++++++++++++++++++++++ 491 G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
508 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquark 492 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquarks<<G4endl;
509 G4cout<<"+++++++++++++++++++++++++++++++ 493 G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
510 #endif 494 #endif
511 495
512 // Proj 496 // Proj
513 TestParticle = G4ParticleTable::GetParti 497 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
514 if ( ! TestParticle ) continue; 498 if ( ! TestParticle ) continue;
515 common.MminProjectile = common.BrW.GetMi 499 common.MminProjectile = common.BrW.GetMinimumMass( TestParticle );
516 if ( common.SqrtS - common.M0target < co 500 if ( common.SqrtS - common.M0target < common.MminProjectile ) continue;
517 MtestPr = common.BrW.SampleMass( TestPar 501 MtestPr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
518 502 + 5.0*TestParticle->GetPDGWidth() );
519 #ifdef debugFTFexcitation << 503 #ifdef debugFTFexictation
520 G4cout << "TestParticle Name " << NewPro 504 G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()
521 << G4endl 505 << G4endl
522 << "MtestPart MtestPart0 "<<Mtest 506 << "MtestPart MtestPart0 "<<MtestPr<<" "<<TestParticle->GetPDGMass()<<G4endl
523 << "M0projectile projectile PDGMa 507 << "M0projectile projectile PDGMass " << common.M0projectile << " "
524 << projectile->GetDefinition()->G 508 << projectile->GetDefinition()->GetPDGMass() << G4endl;
525 #endif 509 #endif
526 510
527 // Targ 511 // Targ
528 NewTargCode = NewNucleonId( TargQ1, Targ 512 NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
529 #ifdef debugFTFexcitation << 513 #ifdef debugFTFexictation
530 G4cout << "New TrQ " << TargQ1 << " " << 514 G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
531 << "NewTargCode " << NewTargCode 515 << "NewTargCode " << NewTargCode << G4endl;
532 #endif 516 #endif
533 517
534 // If the target has not heavy (charm or 518 // If the target has not heavy (charm or bottom) constituent quark,
535 // see whether a Delta isobar can be cre 519 // see whether a Delta isobar can be created.
536 if ( TargQ1 <= 3 && TargQ2 <= 3 && T 520 if ( TargQ1 <= 3 && TargQ2 <= 3 && TargQ3 <= 3 ) {
537 if ( TargQ1 != TargQ2 && TargQ1 != T 521 if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) { // Lambda or Sigma0 ?
538 if ( G4UniformRand() < 0.5 ) { 522 if ( G4UniformRand() < 0.5 ) {
539 NewTargCode += 2; 523 NewTargCode += 2;
540 } else if ( G4UniformRand() < 0.75 ) 524 } else if ( G4UniformRand() < 0.75 ) {
541 NewTargCode = 3122; // Lambda 525 NewTargCode = 3122; // Lambda
542 } 526 }
543 } else if ( TargQ1 == TargQ2 && Targ 527 } else if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
544 NewTargCode += 2; ProjExcited = true 528 NewTargCode += 2; ProjExcited = true; // Create Delta isobar
545 } else if ( target->GetDefinition()->G 529 } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
546 if ( G4UniformRand() > DeltaProbAtQu 530 if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
547 NewTargCode += 2; ProjExcited = tr 531 NewTargCode += 2; ProjExcited = true;
548 } 532 }
549 } else if ( ! ProjExcited && 533 } else if ( ! ProjExcited &&
550 G4UniformRand() < DeltaPro 534 G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
551 common.SqrtS > common.M0pr 535 common.SqrtS > common.M0projectile + // Delta mass
552 G4ParticleTable::GetP 536 G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass() ) {
553 NewTargCode += 2; // Create Delta i 537 NewTargCode += 2; // Create Delta isobar
554 } 538 }
555 } 539 }
556 540
557 // Protection against: 541 // Protection against:
558 // - Lambda* (i.e. excited Lambda state) 542 // - Lambda* (i.e. excited Lambda state) NOT existing in PDG , -> Lambda
559 // - Sigma* (i.e. excited Sigma hyperon 543 // - Sigma* (i.e. excited Sigma hyperon states) NOT existing in Geant4 -> Sigma
560 // - Xi* (i.e. excited Xi hyperon st 544 // - Xi* (i.e. excited Xi hyperon states) NOT existing in Geant4 -> Xi
561 if ( NewTargCode == 3124 || // Lambda 545 if ( NewTargCode == 3124 || // Lambda* NOT existing in PDG !
562 NewTargCode == 3224 || // Sigma*+ NOT 546 NewTargCode == 3224 || // Sigma*+ NOT existing in Geant4
563 NewTargCode == 3214 || // Sigma*0 NOT 547 NewTargCode == 3214 || // Sigma*0 NOT existing in Geant4
564 NewTargCode == 3114 || // Sigma*- NOT 548 NewTargCode == 3114 || // Sigma*- NOT existing in Geant4
565 NewTargCode == 3324 || // Xi*0 NOT 549 NewTargCode == 3324 || // Xi*0 NOT existing in Geant4
566 NewTargCode == 3314 ) { // Xi*- NOT 550 NewTargCode == 3314 ) { // Xi*- NOT existing in Geant4
567 //G4cout << "G4DiffractiveExcitation::Excite 551 //G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange INEXISTING TARGET PDGcode="
568 // << NewTargCode << G4endl; 552 // << NewTargCode << G4endl;
569 NewTargCode -= 2; // Corresponding ground-s 553 NewTargCode -= 2; // Corresponding ground-state hyperon
570 } 554 }
571 555
572 // Special treatment for charmed and bot 556 // Special treatment for charmed and bottom baryons : in Geant4 there are no Xi_c' and Xi_b'
573 // so we need to transform them by hand 557 // so we need to transform them by hand to Xi_c and Xi_b, respectively.
574 #ifdef debug_heavyHadrons 558 #ifdef debug_heavyHadrons
575 G4int initialNewTargCode = NewTargCode; 559 G4int initialNewTargCode = NewTargCode;
576 #endif 560 #endif
577 if ( NewTargCode == 4322 ) NewTargC 561 if ( NewTargCode == 4322 ) NewTargCode = 4232; // Xi_c'+ -> Xi_c+
578 else if ( NewTargCode == 4312 ) NewTargC 562 else if ( NewTargCode == 4312 ) NewTargCode = 4132; // Xi_c'0 -> Xi_c0
579 else if ( NewTargCode == 5312 ) NewTargC 563 else if ( NewTargCode == 5312 ) NewTargCode = 5132; // Xi_b'- -> Xi_b-
580 else if ( NewTargCode == 5322 ) NewTargC 564 else if ( NewTargCode == 5322 ) NewTargCode = 5232; // Xi_b'0 -> Xi_b0
581 #ifdef debug_heavyHadrons 565 #ifdef debug_heavyHadrons
582 if ( NewTargCode != initialNewTargCode ) 566 if ( NewTargCode != initialNewTargCode ) {
583 G4cout << "G4DiffractiveExcitation::ExcitePa 567 G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange : forcing (inexisting in G4)" << G4endl
584 << "\t target heavy baryon with pdgCo 568 << "\t target heavy baryon with pdgCode=" << initialNewTargCode
585 << " into pdgCode=" << NewTargCode << 569 << " into pdgCode=" << NewTargCode << G4endl;
586 } 570 }
587 #endif 571 #endif
588 572
589 TestParticle = G4ParticleTable::GetParti 573 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
590 if ( ! TestParticle ) continue; 574 if ( ! TestParticle ) continue;
591 #ifdef debugFTFexcitation << 575 #ifdef debugFTFexictation
592 G4cout << "New targ " << NewTargCode << 576 G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
593 #endif 577 #endif
594 common.MminTarget = common.BrW.GetMinimu 578 common.MminTarget = common.BrW.GetMinimumMass( TestParticle );
595 if ( common.SqrtS - MtestPr < common.Mmi 579 if ( common.SqrtS - MtestPr < common.MminTarget ) continue;
596 MtestTr = common.BrW.SampleMass( TestPar 580 MtestTr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
597 581 + 5.0*TestParticle->GetPDGWidth() );
598 if ( common.SqrtS > MtestPr + MtestTr ) 582 if ( common.SqrtS > MtestPr + MtestTr ) break;
599 583
600 } // End of while loop 584 } // End of while loop
601 if ( attempts >= maxNumberOfAttempts ) ret 585 if ( attempts >= maxNumberOfAttempts ) return returnCode; // unsuccessfully ended, nothing else can be done
602 586
603 if ( MtestPr >= common.Pprojectile.mag() 587 if ( MtestPr >= common.Pprojectile.mag() || projectile->GetStatus() != 0 ) {
604 common.M0projectile = MtestPr; 588 common.M0projectile = MtestPr;
605 } 589 }
606 #ifdef debugFTFexcitation << 590 #ifdef debugFTFexictation
607 G4cout << "M0projectile After check " << c 591 G4cout << "M0projectile After check " << common.M0projectile << G4endl;
608 #endif 592 #endif
609 common.M0projectile2 = common.M0projectile 593 common.M0projectile2 = common.M0projectile * common.M0projectile;
610 common.ProjectileDiffStateMinMass = com 594 common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
611 common.ProjectileNonDiffStateMinMass = com 595 common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
612 if ( MtestTr >= common.Ptarget.mag() || 596 if ( MtestTr >= common.Ptarget.mag() || target->GetStatus() != 0 ) {
613 common.M0target = MtestTr; 597 common.M0target = MtestTr;
614 } 598 }
615 common.M0target2 = common.M0target * commo 599 common.M0target2 = common.M0target * common.M0target;
616 #ifdef debugFTFexcitation << 600 #ifdef debugFTFexictation
617 G4cout << "New targ M0 M0^2 " << common.M0 601 G4cout << "New targ M0 M0^2 " << common.M0target << " " << common.M0target2 << G4endl;
618 #endif 602 #endif
619 common.TargetDiffStateMinMass = common. 603 common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
620 common.TargetNonDiffStateMinMass = common. 604 common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
621 605
622 } else { // of the if ( common.absProjectil 606 } else { // of the if ( common.absProjectilePDGcode < 1000 ) ; the projectile is baryon now
623 607
624 // If it is a projectile anti-baryon, no q 608 // If it is a projectile anti-baryon, no quark exchange is possible with a target hadron,
625 // therefore returns immediately 1 (which 609 // therefore returns immediately 1 (which means "successfully completed, but the work
626 // needs to be continued"). 610 // needs to be continued").
627 if ( common.ProjectilePDGcode < 0 ) return 611 if ( common.ProjectilePDGcode < 0 ) return 1;
628 612
629 // Choose randomly, with equal probability 613 // Choose randomly, with equal probability, whether to consider the quarks of the
630 // projectile or target hadron for selecti 614 // projectile or target hadron for selecting the flavour of the exchanged quark.
631 G4bool isProjectileExchangedQ = false; 615 G4bool isProjectileExchangedQ = false;
632 G4int firstQ = TargQ1, secondQ = 616 G4int firstQ = TargQ1, secondQ = TargQ2, thirdQ = TargQ3;
633 G4int otherFirstQ = ProjQ1, otherSecondQ = 617 G4int otherFirstQ = ProjQ1, otherSecondQ = ProjQ2, otherThirdQ = ProjQ3;
634 if ( G4UniformRand() < 0.5 ) { 618 if ( G4UniformRand() < 0.5 ) {
635 isProjectileExchangedQ = true; 619 isProjectileExchangedQ = true;
636 firstQ = ProjQ1; secondQ = Pro 620 firstQ = ProjQ1; secondQ = ProjQ2; thirdQ = ProjQ3;
637 otherFirstQ = TargQ1; otherSecondQ = Tar 621 otherFirstQ = TargQ1; otherSecondQ = TargQ2; otherThirdQ = TargQ3;
638 } 622 }
639 // Choose randomly, with equal probability 623 // Choose randomly, with equal probability, which of the three quarks of the
640 // selected (projectile or target) hadron 624 // selected (projectile or target) hadron is the exhanged quark.
641 G4int exchangedQ = 0; 625 G4int exchangedQ = 0;
642 G4double Ksi = G4UniformRand(); 626 G4double Ksi = G4UniformRand();
643 if ( Ksi < 0.333333 ) { 627 if ( Ksi < 0.333333 ) {
644 exchangedQ = firstQ; 628 exchangedQ = firstQ;
645 } else if ( 0.333333 <= Ksi && Ksi < 0.6 629 } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
646 exchangedQ = secondQ; 630 exchangedQ = secondQ;
647 } else { 631 } else {
648 exchangedQ = thirdQ; 632 exchangedQ = thirdQ;
649 } 633 }
650 #ifdef debugFTFexcitation << 634 #ifdef debugFTFexictation
651 G4cout << "Exchange Qs isProjectile Q " << 635 G4cout << "Exchange Qs isProjectile Q " << isProjectileExchangedQ << " " << exchangedQ << " ";
652 #endif 636 #endif
653 // The exchanged quarks (one of the projec 637 // The exchanged quarks (one of the projectile hadron and one of the target hadron)
654 // are always accepted if they have differ 638 // are always accepted if they have different flavour, else (i.e. when they have the
655 // same flavour) they are accepted only wi 639 // same flavour) they are accepted only with a specified probability.
656 G4double probSame = theParameters->GetProb 640 G4double probSame = theParameters->GetProbOfSameQuarkExchange();
657 const G4int MaxCount = 100; 641 const G4int MaxCount = 100;
658 G4int count = 0, otherExchangedQ = 0; 642 G4int count = 0, otherExchangedQ = 0;
659 do { 643 do {
660 if ( exchangedQ != otherFirstQ || G4Un 644 if ( exchangedQ != otherFirstQ || G4UniformRand() < probSame ) {
661 otherExchangedQ = otherFirstQ; otherFi 645 otherExchangedQ = otherFirstQ; otherFirstQ = exchangedQ; exchangedQ = otherExchangedQ;
662 } else { 646 } else {
663 if ( exchangedQ != otherSecondQ || G 647 if ( exchangedQ != otherSecondQ || G4UniformRand() < probSame ) {
664 otherExchangedQ = otherSecondQ; othe 648 otherExchangedQ = otherSecondQ; otherSecondQ = exchangedQ; exchangedQ = otherExchangedQ;
665 } else { 649 } else {
666 if ( exchangedQ != otherThirdQ || 650 if ( exchangedQ != otherThirdQ || G4UniformRand() < probSame ) {
667 otherExchangedQ = otherThirdQ; oth 651 otherExchangedQ = otherThirdQ; otherThirdQ = exchangedQ; exchangedQ = otherExchangedQ;
668 } 652 }
669 } 653 }
670 } 654 }
671 } while ( otherExchangedQ == 0 && ++coun 655 } while ( otherExchangedQ == 0 && ++count < MaxCount );
672 if ( count >= MaxCount ) return returnCode 656 if ( count >= MaxCount ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
673 // Swap (between projectile and target had 657 // Swap (between projectile and target hadron) the two quarks that have been sampled
674 // as "exchanged" quarks. 658 // as "exchanged" quarks.
675 if ( Ksi < 0.333333 ) { 659 if ( Ksi < 0.333333 ) {
676 firstQ = exchangedQ; 660 firstQ = exchangedQ;
677 } else if ( 0.333333 <= Ksi && Ksi < 0.6 661 } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
678 secondQ = exchangedQ; 662 secondQ = exchangedQ;
679 } else { 663 } else {
680 thirdQ = exchangedQ; 664 thirdQ = exchangedQ;
681 } 665 }
682 if ( isProjectileExchangedQ ) { 666 if ( isProjectileExchangedQ ) {
683 ProjQ1 = firstQ; ProjQ2 = secondQ; 667 ProjQ1 = firstQ; ProjQ2 = secondQ; ProjQ3 = thirdQ;
684 TargQ1 = otherFirstQ; TargQ2 = otherSeco 668 TargQ1 = otherFirstQ; TargQ2 = otherSecondQ; TargQ3 = otherThirdQ;
685 } else { 669 } else {
686 TargQ1 = firstQ; TargQ2 = secondQ; 670 TargQ1 = firstQ; TargQ2 = secondQ; TargQ3 = thirdQ;
687 ProjQ1 = otherFirstQ; ProjQ2 = otherSeco 671 ProjQ1 = otherFirstQ; ProjQ2 = otherSecondQ; ProjQ3 = otherThirdQ;
688 } 672 }
689 #ifdef debugFTFexcitation << 673 #ifdef debugFTFexictation
690 G4cout << "Exchange Qs Pr Tr " << ( isPro 674 G4cout << "Exchange Qs Pr Tr " << ( isProjectileExchangedQ ? exchangedQ : otherExchangedQ )
691 << " " << ( isProjectileExchangedQ 675 << " " << ( isProjectileExchangedQ ? otherExchangedQ : exchangedQ ) << G4endl;
692 #endif 676 #endif
693 677
694 NewProjCode = NewNucleonId( ProjQ1, ProjQ2 678 NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
695 NewTargCode = NewNucleonId( TargQ1, TargQ2 679 NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
696 // Decide whether the new projectile hadro 680 // Decide whether the new projectile hadron is a Delta particle;
697 // then decide whether the new target hadr 681 // then decide whether the new target hadron is a Delta particle.
698 // Notice that a Delta particle has the la 682 // Notice that a Delta particle has the last PDG digit "4" (because its spin is 3/2),
699 // whereas a nucleon has "2" (because its 683 // whereas a nucleon has "2" (because its spin is 1/2).
700 // If the new projectile hadron or the new 684 // If the new projectile hadron or the new target hadron has a heavy (c or b)
701 // constituent quark, then skip this part 685 // constituent quark, then skip this part (because Geant4 does not have
702 // excited charm and bottom hadrons). 686 // excited charm and bottom hadrons).
703 for ( G4int iHadron = 0; iHadron < 2; iHad 687 for ( G4int iHadron = 0; iHadron < 2; iHadron++ ) {
704 // First projectile hadron, then target 688 // First projectile hadron, then target hadron
705 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, 689 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, codeQ3 = ProjQ3, newHadCode = NewProjCode;
706 G4double massConstraint = common.M0targe 690 G4double massConstraint = common.M0target;
707 G4bool isHadronADelta = ( projectile->Ge 691 G4bool isHadronADelta = ( projectile->GetDefinition()->GetPDGiIsospin() == 3 );
708 if ( iHadron == 1 ) { // Target hadron 692 if ( iHadron == 1 ) { // Target hadron
709 codeQ1 = TargQ1, codeQ2 = TargQ2, code 693 codeQ1 = TargQ1, codeQ2 = TargQ2, codeQ3 = TargQ3, newHadCode = NewTargCode;
710 massConstraint = common.M0projectile; 694 massConstraint = common.M0projectile;
711 isHadronADelta = ( target->GetDefiniti 695 isHadronADelta = ( target->GetDefinition()->GetPDGiIsospin() == 3 );
712 } 696 }
713 if ( codeQ1 > 3 || codeQ2 > 3 || cod 697 if ( codeQ1 > 3 || codeQ2 > 3 || codeQ3 > 3 ) continue; // No excited charm or bottom states in Geant4
714 if ( codeQ1 == codeQ2 && codeQ1 == cod 698 if ( codeQ1 == codeQ2 && codeQ1 == codeQ3 ) { // The three quarks are the same
715 newHadCode += 2; // Delta++ (uuu) or 699 newHadCode += 2; // Delta++ (uuu) or Delta- (ddd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
716 } else if ( isHadronADelta ) { // Hadro 700 } else if ( isHadronADelta ) { // Hadron (projectile or target) was Delta
717 if ( G4UniformRand() > DeltaProbAtQuar 701 if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
718 newHadCode += 2; // Delta+ (uud) or 702 newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
719 } else { 703 } else {
720 newHadCode += 0; // No delta (so th 704 newHadCode += 0; // No delta (so the last PDG digit remains 2)
721 } 705 }
722 } else { // Hadron (projectile or targe 706 } else { // Hadron (projectile or target) was Nucleon
723 if ( G4UniformRand() < DeltaProbAtQuar 707 if ( G4UniformRand() < DeltaProbAtQuarkExchange &&
724 common.SqrtS > G4ParticleTable::G 708 common.SqrtS > G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass()
725 + massConstraint ) 709 + massConstraint ) {
726 newHadCode += 2; // Delta+ (uud) or 710 newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
727 } else { 711 } else {
728 newHadCode += 0; // No delta (so th 712 newHadCode += 0; // No delta (so the last PDG digit remains 2)
729 } 713 }
730 } 714 }
731 if ( iHadron == 0 ) { // Projectile had 715 if ( iHadron == 0 ) { // Projectile hadron
732 NewProjCode = newHadCode; 716 NewProjCode = newHadCode;
733 } else { // Target hadron 717 } else { // Target hadron
734 NewTargCode = newHadCode; 718 NewTargCode = newHadCode;
735 } 719 }
736 } 720 }
737 #ifdef debugFTFexcitation << 721 #ifdef debugFTFexictation
738 G4cout << "NewProjCode NewTargCode " << Ne 722 G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
739 #endif 723 #endif
740 724
741 // Protection against: 725 // Protection against:
742 // - Lambda* (i.e. excited Lambda state) 726 // - Lambda* (i.e. excited Lambda state) NOT existing in PDG , -> Lambda
743 // - Sigma* (i.e. excited Sigma hyperon s 727 // - Sigma* (i.e. excited Sigma hyperon states) NOT existing in Geant4 -> Sigma
744 // - Xi* (i.e. excited Xi hyperon stat 728 // - Xi* (i.e. excited Xi hyperon states) NOT existing in Geant4 -> Xi
745 if ( NewProjCode == 3124 || // Lambda* 729 if ( NewProjCode == 3124 || // Lambda* NOT existing in PDG !
746 NewProjCode == 3224 || // Sigma*+ NOT ex 730 NewProjCode == 3224 || // Sigma*+ NOT existing in Geant4
747 NewProjCode == 3214 || // Sigma*0 NOT ex 731 NewProjCode == 3214 || // Sigma*0 NOT existing in Geant4
748 NewProjCode == 3114 || // Sigma*- NOT ex 732 NewProjCode == 3114 || // Sigma*- NOT existing in Geant4
749 NewProjCode == 3324 || // Xi*0 NOT ex 733 NewProjCode == 3324 || // Xi*0 NOT existing in Geant4
750 NewProjCode == 3314 ) { // Xi*- NOT ex 734 NewProjCode == 3314 ) { // Xi*- NOT existing in Geant4
751 //G4cout << "G4DiffractiveExcitation::Ex 735 //G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange INEXISTING PROJECTILE PDGcode="
752 // << NewProjCode << G4endl; 736 // << NewProjCode << G4endl;
753 NewProjCode -= 2; // Corresponding grou 737 NewProjCode -= 2; // Corresponding ground-state hyperon
754 } 738 }
755 if ( NewTargCode == 3124 || // Lambda* 739 if ( NewTargCode == 3124 || // Lambda* NOT existing in PDG !
756 NewTargCode == 3224 || // Sigma*+ NOT ex 740 NewTargCode == 3224 || // Sigma*+ NOT existing in Geant4
757 NewTargCode == 3214 || // Sigma*0 NOT ex 741 NewTargCode == 3214 || // Sigma*0 NOT existing in Geant4
758 NewTargCode == 3114 || // Sigma*- NOT ex 742 NewTargCode == 3114 || // Sigma*- NOT existing in Geant4
759 NewTargCode == 3324 || // Xi*0 NOT ex 743 NewTargCode == 3324 || // Xi*0 NOT existing in Geant4
760 NewTargCode == 3314 ) { // Xi*- NOT ex 744 NewTargCode == 3314 ) { // Xi*- NOT existing in Geant4
761 //G4cout << "G4DiffractiveExcitation::Ex 745 //G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange INEXISTING TARGET PDGcode="
762 // << NewTargCode << G4endl; 746 // << NewTargCode << G4endl;
763 NewTargCode -= 2; // Corresponding grou 747 NewTargCode -= 2; // Corresponding ground-state hyperon
764 } 748 }
765 749
766 // Special treatment for charmed and botto 750 // Special treatment for charmed and bottom baryons : in Geant4 there are no Xi_c' and Xi_b'
767 // so we need to transform them by hand to 751 // so we need to transform them by hand to the, respectively, Xi_c and Xi_b.
768 #ifdef debug_heavyHadrons 752 #ifdef debug_heavyHadrons
769 G4int initialNewProjCode = NewProjCode, in 753 G4int initialNewProjCode = NewProjCode, initialNewTargCode = NewTargCode;
770 #endif 754 #endif
771 if ( NewProjCode == 4322 ) NewProjCod 755 if ( NewProjCode == 4322 ) NewProjCode = 4232; // Xi_c'+ -> Xi_c+
772 else if ( NewProjCode == 4312 ) NewProjCod 756 else if ( NewProjCode == 4312 ) NewProjCode = 4132; // Xi_c'0 -> Xi_c0
773 else if ( NewProjCode == 5312 ) NewProjCod 757 else if ( NewProjCode == 5312 ) NewProjCode = 5132; // Xi_b'- -> Xi_b-
774 else if ( NewProjCode == 5322 ) NewProjCod 758 else if ( NewProjCode == 5322 ) NewProjCode = 5232; // Xi_b'0 -> Xi_b0
775 if ( NewTargCode == 4322 ) NewTargCod 759 if ( NewTargCode == 4322 ) NewTargCode = 4232; // Xi_c'+ -> Xi_c+
776 else if ( NewTargCode == 4312 ) NewTargCod 760 else if ( NewTargCode == 4312 ) NewTargCode = 4132; // Xi_c'0 -> Xi_c0
777 else if ( NewTargCode == 5312 ) NewTargCod 761 else if ( NewTargCode == 5312 ) NewTargCode = 5132; // Xi_b'- -> Xi_b-
778 else if ( NewTargCode == 5322 ) NewTargCod 762 else if ( NewTargCode == 5322 ) NewTargCode = 5232; // Xi_b'0 -> Xi_b0
779 #ifdef debug_heavyHadrons 763 #ifdef debug_heavyHadrons
780 if ( NewProjCode != initialNewProjCode || 764 if ( NewProjCode != initialNewProjCode || NewTargCode != initialNewTargCode ) {
781 G4cout << "G4DiffractiveExcitation::Exci 765 G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange : forcing (inexisting in G4)" << G4endl
782 << "\t heavy baryon into an exist 766 << "\t heavy baryon into an existing one:" << G4endl;
783 if ( NewProjCode != initialNewProjCode ) 767 if ( NewProjCode != initialNewProjCode ) {
784 G4cout << "\t \t projectile baryon with pdgC 768 G4cout << "\t \t projectile baryon with pdgCode=" << initialNewProjCode
785 << " into pdgCode=" << NewProjCode << 769 << " into pdgCode=" << NewProjCode << G4endl;
786 } 770 }
787 if ( NewTargCode != initialNewTargCode ) 771 if ( NewTargCode != initialNewTargCode ) {
788 G4cout << "\t \t target baryon with pd 772 G4cout << "\t \t target baryon with pdgCode=" << initialNewTargCode
789 << " into pdgCode=" << NewTargCode << 773 << " into pdgCode=" << NewTargCode << G4endl;
790 } 774 }
791 } 775 }
792 #endif 776 #endif
793 777
794 // Sampling of the masses of the projectil 778 // Sampling of the masses of the projectile and target nucleons.
795 // Because of energy conservation, the ord 779 // Because of energy conservation, the ordering of the sampling matters:
796 // randomly, half of the time we sample fi 780 // randomly, half of the time we sample first the target nucleon mass and
797 // then the projectile nucleon mass, and t 781 // then the projectile nucleon mass, and the other half of the time we
798 // sample first the projectile nucleon mas 782 // sample first the projectile nucleon mass and then the target nucleon mass.
799 G4VSplitableHadron* firstHadron = target; 783 G4VSplitableHadron* firstHadron = target;
800 G4VSplitableHadron* secondHadron = project 784 G4VSplitableHadron* secondHadron = projectile;
801 G4int firstHadronCode = NewTargCode, secon 785 G4int firstHadronCode = NewTargCode, secondHadronCode = NewProjCode;
802 G4double massConstraint = common.M0project 786 G4double massConstraint = common.M0projectile;
803 G4bool isFirstTarget = true; 787 G4bool isFirstTarget = true;
804 if ( G4UniformRand() < 0.5 ) { 788 if ( G4UniformRand() < 0.5 ) {
805 // Sample first the projectile nucleon m 789 // Sample first the projectile nucleon mass, then the target nucleon mass.
806 firstHadron = projectile; secondHa 790 firstHadron = projectile; secondHadron = target;
807 firstHadronCode = NewProjCode; secondHa 791 firstHadronCode = NewProjCode; secondHadronCode = NewTargCode;
808 massConstraint = common.M0target; 792 massConstraint = common.M0target;
809 isFirstTarget = false; 793 isFirstTarget = false;
810 } 794 }
811 G4double Mtest1st = 0.0, Mtest2nd = 0.0, M 795 G4double Mtest1st = 0.0, Mtest2nd = 0.0, Mmin1st = 0.0, Mmin2nd = 0.0;
812 for ( int iSamplingCase = 0; iSamplingCase 796 for ( int iSamplingCase = 0; iSamplingCase < 2; iSamplingCase++ ) {
813 G4VSplitableHadron* aHadron = firstHadro 797 G4VSplitableHadron* aHadron = firstHadron;
814 G4int aHadronCode = firstHadronCode; 798 G4int aHadronCode = firstHadronCode;
815 if ( iSamplingCase == 1 ) { // Second n 799 if ( iSamplingCase == 1 ) { // Second nucleon mass sampling
816 aHadron = secondHadron; aHadronCode = 800 aHadron = secondHadron; aHadronCode = secondHadronCode; massConstraint = Mtest1st;
817 } 801 }
818 G4double MtestHadron = 0.0, MminHadron = 802 G4double MtestHadron = 0.0, MminHadron = 0.0;
819 if ( aHadron->GetStatus() == 1 || aHad 803 if ( aHadron->GetStatus() == 1 || aHadron->GetStatus() == 2 ) {
820 TestParticle = G4ParticleTable::GetPar 804 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( aHadronCode );
821 if ( ! TestParticle ) return returnCod 805 if ( ! TestParticle ) return returnCode; // Not possible to find such a hadron: unsuccessfully ended, nothing else can be done
822 MminHadron = common.BrW.GetMinimumMass 806 MminHadron = common.BrW.GetMinimumMass( TestParticle );
823 if ( common.SqrtS - massConstraint < M 807 if ( common.SqrtS - massConstraint < MminHadron ) return returnCode; // Kinematically impossible: unsuccessfully ended, nothing else can be done
824 if ( TestParticle->GetPDGWidth() == 0. 808 if ( TestParticle->GetPDGWidth() == 0.0 ) {
825 MtestHadron = common.BrW.SampleMass( 809 MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass() );
826 } else { 810 } else {
827 const G4int maxNumberOfAttempts = 50 811 const G4int maxNumberOfAttempts = 50;
828 G4int attempts = 0; 812 G4int attempts = 0;
829 while ( attempts < maxNumberOfAttemp 813 while ( attempts < maxNumberOfAttempts ) {
830 attempts++; 814 attempts++;
831 MtestHadron = common.BrW.SampleMas 815 MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
832 816 + 5.0*TestParticle->GetPDGWidth() );
833 if ( common.SqrtS < MtestHadron + 817 if ( common.SqrtS < MtestHadron + massConstraint ) {
834 continue; // Kinematically unac 818 continue; // Kinematically unacceptable: try again
835 } else { 819 } else {
836 break; // Kinematically acce 820 break; // Kinematically acceptable: the mass sampling is successful
837 } 821 }
838 } 822 }
839 if ( attempts >= maxNumberOfAttempts 823 if ( attempts >= maxNumberOfAttempts ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
840 } 824 }
841 } 825 }
842 if ( iSamplingCase == 0 ) { 826 if ( iSamplingCase == 0 ) {
843 Mtest1st = MtestHadron; Mmin1st = Mmi 827 Mtest1st = MtestHadron; Mmin1st = MminHadron;
844 } else { 828 } else {
845 Mtest2nd = MtestHadron; Mmin2nd = Mmi 829 Mtest2nd = MtestHadron; Mmin2nd = MminHadron;
846 } 830 }
847 } // End for loop on the two sampling cas 831 } // End for loop on the two sampling cases (1st and 2nd)
848 if ( isFirstTarget ) { 832 if ( isFirstTarget ) {
849 MtestTr = Mtest1st; MtestPr = Mtest2n 833 MtestTr = Mtest1st; MtestPr = Mtest2nd;
850 common.MminTarget = Mmin1st; common.Mmi 834 common.MminTarget = Mmin1st; common.MminProjectile = Mmin2nd;
851 } else { 835 } else {
852 MtestTr = Mtest2nd; MtestPr = Mtest1s 836 MtestTr = Mtest2nd; MtestPr = Mtest1st;
853 common.MminTarget = Mmin2nd; common.Mmi 837 common.MminTarget = Mmin2nd; common.MminProjectile = Mmin1st;
854 } 838 }
855 839
856 if ( MtestPr != 0.0 ) { 840 if ( MtestPr != 0.0 ) {
857 common.M0projectile = MtestPr; 841 common.M0projectile = MtestPr;
858 common.M0projectile2 = common.M0projecti 842 common.M0projectile2 = common.M0projectile * common.M0projectile;
859 common.ProjectileDiffStateMinMass = c 843 common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
860 common.ProjectileNonDiffStateMinMass = c 844 common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
861 } 845 }
862 if ( MtestTr != 0.0 ) { 846 if ( MtestTr != 0.0 ) {
863 common.M0target = MtestTr; 847 common.M0target = MtestTr;
864 common.M0target2 = common.M0target * com 848 common.M0target2 = common.M0target * common.M0target;
865 common.TargetDiffStateMinMass = commo 849 common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
866 common.TargetNonDiffStateMinMass = commo 850 common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
867 } 851 }
868 852
869 } // End of if ( common.absProjectilePDGcod 853 } // End of if ( common.absProjectilePDGcode < 1000 )
870 854
871 // If we assume that final state hadrons aft 855 // If we assume that final state hadrons after the charge exchange will be
872 // in the ground states, we have to put 856 // in the ground states, we have to put
873 if ( common.SqrtS < common.M0projectile + co 857 if ( common.SqrtS < common.M0projectile + common.M0target ) return returnCode; // unsuccessfully ended, nothing else can be done
874 858
875 common.PZcms2 = ( sqr( common.S ) + sqr( com 859 common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
876 - 2.0 * ( common.S * ( com 860 - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
877 + common.M0proje 861 + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
878 #ifdef debugFTFexcitation << 862 #ifdef debugFTFexictation
879 G4cout << "At the end// NewProjCode " << New 863 G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
880 << "At the end// NewTargCode " << New 864 << "At the end// NewTargCode " << NewTargCode << G4endl
881 << "M0pr M0tr SqS " << common.M0pro 865 << "M0pr M0tr SqS " << common.M0projectile << " " << common.M0target << " "
882 << common.SqrtS << G4endl 866 << common.SqrtS << G4endl
883 << "M0pr2 M0tr2 SqS " << common.M0pro 867 << "M0pr2 M0tr2 SqS " << common.M0projectile2 << " " << common.M0target2 << " "
884 << common.SqrtS << G4endl 868 << common.SqrtS << G4endl
885 << "PZcms2 after the change " << comm 869 << "PZcms2 after the change " << common.PZcms2 << G4endl << G4endl;
886 #endif 870 #endif
887 if ( common.PZcms2 < 0.0 ) return returnCode 871 if ( common.PZcms2 < 0.0 ) return returnCode; // It can be if energy is not sufficient for Delta;
888 872 // unsuccessfully ended, nothing else can be done
889 projectile->SetDefinition( G4ParticleTable:: 873 projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
890 target->SetDefinition( G4ParticleTable::GetP 874 target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
891 common.PZcms = std::sqrt( common.PZcms2 ); 875 common.PZcms = std::sqrt( common.PZcms2 );
892 common.Pprojectile.setPz( common.PZcms ); 876 common.Pprojectile.setPz( common.PZcms );
893 common.Pprojectile.setE( std::sqrt( common.M 877 common.Pprojectile.setE( std::sqrt( common.M0projectile2 + common.PZcms2 ) );
894 common.Ptarget.setPz( -common.PZcms ); 878 common.Ptarget.setPz( -common.PZcms );
895 common.Ptarget.setE( std::sqrt( common.M0tar 879 common.Ptarget.setE( std::sqrt( common.M0target2 + common.PZcms2 ) );
896 #ifdef debugFTFexcitation << 880 #ifdef debugFTFexictation
897 G4cout << "Proj Targ and Proj+Targ in CMS" < 881 G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << common.Pprojectile << G4endl
898 << common.Ptarget << G4endl << common 882 << common.Ptarget << G4endl << common.Pprojectile + common.Ptarget << G4endl;
899 #endif 883 #endif
900 884
901 if ( projectile->GetStatus() != 0 ) projecti 885 if ( projectile->GetStatus() != 0 ) projectile->SetStatus( 2 );
902 if ( target->GetStatus() != 0 ) target-> 886 if ( target->GetStatus() != 0 ) target->SetStatus( 2 );
903 887
904 // Check for possible excitation of the part 888 // Check for possible excitation of the participants
905 if ( common.SqrtS < common.M0projectile + co 889 if ( common.SqrtS < common.M0projectile + common.TargetDiffStateMinMass ||
906 common.SqrtS < common.ProjectileDiffSta 890 common.SqrtS < common.ProjectileDiffStateMinMass + common.M0target ||
907 common.ProbOfDiffraction == 0.0 ) commo 891 common.ProbOfDiffraction == 0.0 ) common.ProbExc = 0.0;
908 892
909 if ( G4UniformRand() > common.ProbExc ) { / 893 if ( G4UniformRand() > common.ProbExc ) { // Make elastic scattering
910 #ifdef debugFTFexcitation << 894 #ifdef debugFTFexictation
911 G4cout << "Make elastic scattering of new 895 G4cout << "Make elastic scattering of new hadrons" << G4endl;
912 #endif 896 #endif
913 common.Pprojectile.transform( common.toLab 897 common.Pprojectile.transform( common.toLab );
914 common.Ptarget.transform( common.toLab ); 898 common.Ptarget.transform( common.toLab );
915 projectile->Set4Momentum( common.Pprojecti 899 projectile->Set4Momentum( common.Pprojectile );
916 target->Set4Momentum( common.Ptarget ); 900 target->Set4Momentum( common.Ptarget );
917 G4bool Result = theElastic->ElasticScatter 901 G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
918 #ifdef debugFTFexcitation << 902 #ifdef debugFTFexictation
919 G4cout << "Result of el. scatt " << Result 903 G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
920 << G4endl << projectile->Get4Moment 904 << G4endl << projectile->Get4Momentum() << G4endl << target->Get4Momentum() << G4endl
921 << projectile->Get4Momentum() + tar 905 << projectile->Get4Momentum() + target->Get4Momentum() << " "
922 << (projectile->Get4Momentum() + ta 906 << (projectile->Get4Momentum() + target->Get4Momentum()).mag() << G4endl;
923 #endif 907 #endif
924 if ( Result ) returnCode = 0; // successf 908 if ( Result ) returnCode = 0; // successfully ended and nothing else needs to be done
925 return returnCode; 909 return returnCode;
926 } 910 }
927 911
928 #ifdef debugFTFexcitation << 912 #ifdef debugFTFexictation
929 G4cout << "Make excitation of new hadrons" < 913 G4cout << "Make excitation of new hadrons" << G4endl;
930 #endif 914 #endif
931 915
932 // Redefinition of ProbOfDiffraction because 916 // Redefinition of ProbOfDiffraction because the probabilities are changed due to quark exchange
933 common.ProbOfDiffraction = common.ProbProjec 917 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
934 if ( common.ProbOfDiffraction != 0.0 ) { 918 if ( common.ProbOfDiffraction != 0.0 ) {
935 common.ProbProjectileDiffraction /= common 919 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
936 common.ProbTargetDiffraction /= common 920 common.ProbTargetDiffraction /= common.ProbOfDiffraction;
937 } 921 }
938 922
939 return returnCode = 1; // successfully comp 923 return returnCode = 1; // successfully completed, but the work needs to be continued
940 } 924 }
941 925
942 //-------------------------------------------- 926 //-----------------------------------------------------------------------------
943 927
944 G4bool G4DiffractiveExcitation:: 928 G4bool G4DiffractiveExcitation::
945 ExciteParticipants_doDiffraction( G4VSplitable 929 ExciteParticipants_doDiffraction( G4VSplitableHadron* projectile, G4VSplitableHadron* target,
946 G4FTFParamet 930 G4FTFParameters* theParameters,
947 G4Diffractiv 931 G4DiffractiveExcitation::CommonVariables& common ) const {
948 // Second of the three utility methods used 932 // Second of the three utility methods used only by ExciteParticipants:
949 // it does the sampling for the diffraction 933 // it does the sampling for the diffraction case, either projectile or target diffraction.
950 934
951 G4bool isProjectileDiffraction = false; 935 G4bool isProjectileDiffraction = false;
952 if ( G4UniformRand() < common.ProbProjectile 936 if ( G4UniformRand() < common.ProbProjectileDiffraction ) { // projectile diffraction
953 isProjectileDiffraction = true; 937 isProjectileDiffraction = true;
954 #ifdef debugFTFexcitation << 938 #ifdef debugFTFexictation
955 G4cout << "projectile diffraction" << G4en 939 G4cout << "projectile diffraction" << G4endl;
956 #endif 940 #endif
957 common.ProjMassT2 = common.ProjectileDiffS 941 common.ProjMassT2 = common.ProjectileDiffStateMinMass2;
958 common.ProjMassT = common.ProjectileDiffS 942 common.ProjMassT = common.ProjectileDiffStateMinMass;
959 common.TargMassT2 = common.M0target2; 943 common.TargMassT2 = common.M0target2;
960 common.TargMassT = common.M0target; 944 common.TargMassT = common.M0target;
961 } else { 945 } else { // target diffraction
962 #ifdef debugFTFexcitation << 946 #ifdef debugFTFexictation
963 G4cout << "Target diffraction" << G4endl; 947 G4cout << "Target diffraction" << G4endl;
964 #endif 948 #endif
965 common.ProjMassT2 = common.M0projectile2; 949 common.ProjMassT2 = common.M0projectile2;
966 common.ProjMassT = common.M0projectile; 950 common.ProjMassT = common.M0projectile;
967 common.TargMassT2 = common.TargetDiffState 951 common.TargMassT2 = common.TargetDiffStateMinMass2;
968 common.TargMassT = common.TargetDiffState 952 common.TargMassT = common.TargetDiffStateMinMass;
969 } 953 }
970 954
971 // Check whether the interaction is possible 955 // Check whether the interaction is possible
972 if ( common.SqrtS < common.ProjMassT + commo 956 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
973 957
974 common.PZcms2 = ( sqr( common.S ) + sqr( com 958 common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
975 - 2.0 * ( common.S * ( com 959 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
976 + common.ProjMassT2 960 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
977 if ( common.PZcms2 < 0.0 ) return false; 961 if ( common.PZcms2 < 0.0 ) return false;
978 common.maxPtSquare = common.PZcms2; 962 common.maxPtSquare = common.PZcms2;
979 963
980 G4double DiffrAveragePt2 = theParameters->Ge 964 G4double DiffrAveragePt2 = theParameters->GetAvaragePt2ofElasticScattering()*1.2;
981 G4bool loopCondition = true; 965 G4bool loopCondition = true;
982 G4int whilecount = 0; 966 G4int whilecount = 0;
983 do { // Generate pt and mass of projectile 967 do { // Generate pt and mass of projectile
984 968
985 whilecount++; 969 whilecount++;
986 if ( whilecount > 1000 ) { 970 if ( whilecount > 1000 ) {
987 common.Qmomentum = G4LorentzVector( 0.0, 971 common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
988 return false; // Ignore this interacti 972 return false; // Ignore this interaction
989 }; 973 };
990 974
991 common.Qmomentum = G4LorentzVector( Gaussi 975 common.Qmomentum = G4LorentzVector( GaussianPt( DiffrAveragePt2, common.maxPtSquare ), 0 );
992 common.Pt2 = G4ThreeVector( common.Qmoment 976 common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
993 if ( isProjectileDiffraction ) { // proje 977 if ( isProjectileDiffraction ) { // projectile diffraction
994 common.ProjMassT2 = common.ProjectileDif 978 common.ProjMassT2 = common.ProjectileDiffStateMinMass2 + common.Pt2;
995 common.TargMassT2 = common.M0target2 + c 979 common.TargMassT2 = common.M0target2 + common.Pt2;
996 } else { // targe 980 } else { // target diffraction
997 common.ProjMassT2 = common.M0projectile2 981 common.ProjMassT2 = common.M0projectile2 + common.Pt2;
998 common.TargMassT2 = common.TargetDiffSta 982 common.TargMassT2 = common.TargetDiffStateMinMass2 + common.Pt2;
999 } 983 }
1000 common.ProjMassT = std::sqrt( common.Proj 984 common.ProjMassT = std::sqrt( common.ProjMassT2 );
1001 common.TargMassT = std::sqrt( common.Targ 985 common.TargMassT = std::sqrt( common.TargMassT2 );
1002 if ( common.SqrtS < common.ProjMassT + co 986 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
1003 987
1004 common.PZcms2 = ( sqr( common.S ) + sqr( 988 common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
1005 - 2.0 * ( common.S * ( 989 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1006 + common.Proj 990 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1007 if ( common.PZcms2 < 0.0 ) continue; 991 if ( common.PZcms2 < 0.0 ) continue;
1008 992
1009 common.PZcms = std::sqrt( common.PZcms2 ) 993 common.PZcms = std::sqrt( common.PZcms2 );
1010 if ( isProjectileDiffraction ) { // proj 994 if ( isProjectileDiffraction ) { // projectile diffraction
1011 common.PMinusMin = std::sqrt( common.Pr 995 common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
1012 common.PMinusMax = common.SqrtS - commo 996 common.PMinusMax = common.SqrtS - common.TargMassT;
1013 common.PMinusNew = ChooseP( common.PMin 997 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1014 common.TMinusNew = common.SqrtS - commo 998 common.TMinusNew = common.SqrtS - common.PMinusNew;
1015 common.Qminus = common.Ptarget.minus() 999 common.Qminus = common.Ptarget.minus() - common.TMinusNew;
1016 common.TPlusNew = common.TargMassT2 / c 1000 common.TPlusNew = common.TargMassT2 / common.TMinusNew;
1017 common.Qplus = common.Ptarget.plus() - 1001 common.Qplus = common.Ptarget.plus() - common.TPlusNew;
1018 common.Qmomentum.setPz( (common.Qplus - 1002 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1019 common.Qmomentum.setE( (common.Qplus + 1003 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1020 loopCondition = ( ( common.Pprojectile 1004 loopCondition = ( ( common.Pprojectile + common.Qmomentum ).mag2() <
1021 common.ProjectileDiff 1005 common.ProjectileDiffStateMinMass2 );
1022 } else { // targ 1006 } else { // target diffraction
1023 common.TPlusMin = std::sqrt( common.Tar 1007 common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
1024 common.TPlusMax = common.SqrtS - common 1008 common.TPlusMax = common.SqrtS - common.ProjMassT;
1025 common.TPlusNew = ChooseP( common.TPlus 1009 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1026 common.PPlusNew = common.SqrtS - common 1010 common.PPlusNew = common.SqrtS - common.TPlusNew;
1027 common.Qplus = common.PPlusNew - common 1011 common.Qplus = common.PPlusNew - common.Pprojectile.plus();
1028 common.PMinusNew = common.ProjMassT2 / 1012 common.PMinusNew = common.ProjMassT2 / common.PPlusNew;
1029 common.Qminus = common.PMinusNew - comm 1013 common.Qminus = common.PMinusNew - common.Pprojectile.minus();
1030 common.Qmomentum.setPz( (common.Qplus - 1014 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1031 common.Qmomentum.setE( (common.Qplus + 1015 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1032 loopCondition = ( ( common.Ptarget - c 1016 loopCondition = ( ( common.Ptarget - common.Qmomentum ).mag2() <
1033 common.TargetDiffSta 1017 common.TargetDiffStateMinMass2 );
1034 } 1018 }
1035 1019
1036 } while ( loopCondition ); /* Loop checkin 1020 } while ( loopCondition ); /* Loop checking, 10.08.2015, A.Ribon */
1037 // Repeat the sampling because ther 1021 // Repeat the sampling because there was not any excitation
1038 1022
1039 if ( isProjectileDiffraction ) { // projec 1023 if ( isProjectileDiffraction ) { // projectile diffraction
1040 projectile->SetStatus( 0 ); 1024 projectile->SetStatus( 0 );
1041 if ( projectile->GetStatus() == 2 ) proje 1025 if ( projectile->GetStatus() == 2 ) projectile->SetStatus( 1 );
1042 if ( target->GetStatus() == 1 && target 1026 if ( target->GetStatus() == 1 && target->GetSoftCollisionCount() == 0 ) target->SetStatus( 2 );
1043 } else { // target 1027 } else { // target diffraction
1044 target->SetStatus( 0 ); 1028 target->SetStatus( 0 );
1045 } 1029 }
1046 1030
1047 return true; 1031 return true;
1048 } 1032 }
1049 1033
1050 //------------------------------------------- 1034 //-----------------------------------------------------------------------------
1051 1035
1052 G4bool G4DiffractiveExcitation:: 1036 G4bool G4DiffractiveExcitation::
1053 ExciteParticipants_doNonDiffraction( G4VSplit 1037 ExciteParticipants_doNonDiffraction( G4VSplitableHadron* projectile,
1054 G4VSplit 1038 G4VSplitableHadron* target,
1055 G4FTFPar 1039 G4FTFParameters* theParameters,
1056 G4Diffra 1040 G4DiffractiveExcitation::CommonVariables& common ) const {
1057 // Third of the three utility methods used 1041 // Third of the three utility methods used only by ExciteParticipants:
1058 // it does the sampling for the non-diffrac 1042 // it does the sampling for the non-diffraction case.
1059 1043
1060 #ifdef debugFTFexcitation << 1044 #ifdef debugFTFexictation
1061 G4cout << "Non-diffraction process" << G4en 1045 G4cout << "Non-diffraction process" << G4endl;
1062 #endif 1046 #endif
1063 1047
1064 // Check whether the interaction is possibl 1048 // Check whether the interaction is possible
1065 common.ProjMassT2 = common.ProjectileNonDif 1049 common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2;
1066 common.ProjMassT = common.ProjectileNonDif 1050 common.ProjMassT = common.ProjectileNonDiffStateMinMass;
1067 common.TargMassT2 = common.TargetNonDiffSta 1051 common.TargMassT2 = common.TargetNonDiffStateMinMass2;
1068 common.TargMassT = common.TargetNonDiffSta 1052 common.TargMassT = common.TargetNonDiffStateMinMass;
1069 if ( common.SqrtS < common.ProjMassT + comm 1053 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
1070 1054
1071 common.PZcms2 = ( sqr( common.S ) + sqr( co 1055 common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
1072 - 2.0 * ( common.S * ( comm 1056 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1073 + common.ProjMassT2 * 1057 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1074 common.maxPtSquare = common.PZcms2; 1058 common.maxPtSquare = common.PZcms2;
1075 1059
1076 G4int whilecount = 0; 1060 G4int whilecount = 0;
1077 do { // Generate pt and masses 1061 do { // Generate pt and masses
1078 1062
1079 whilecount++; 1063 whilecount++;
1080 if ( whilecount > 1000 ) { 1064 if ( whilecount > 1000 ) {
1081 common.Qmomentum = G4LorentzVector( 0.0 1065 common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1082 return false; // Ignore this interacti 1066 return false; // Ignore this interaction
1083 }; 1067 };
1084 1068
1085 common.Qmomentum = G4LorentzVector( Gauss 1069 common.Qmomentum = G4LorentzVector( GaussianPt( theParameters->GetAveragePt2(),
1086 1070 common.maxPtSquare ), 0 );
1087 common.Pt2 = G4ThreeVector( common.Qmomen 1071 common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
1088 common.ProjMassT2 = common.ProjectileNonD 1072 common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2 + common.Pt2;
1089 common.ProjMassT = std::sqrt( common.Pro 1073 common.ProjMassT = std::sqrt( common.ProjMassT2 );
1090 common.TargMassT2 = common.TargetNonDiffS 1074 common.TargMassT2 = common.TargetNonDiffStateMinMass2 + common.Pt2;
1091 common.TargMassT = std::sqrt( common.Tar 1075 common.TargMassT = std::sqrt( common.TargMassT2 );
1092 if ( common.SqrtS < common.ProjMassT + co 1076 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
1093 1077
1094 common.PZcms2 =( sqr( common.S ) + sqr( c 1078 common.PZcms2 =( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
1095 - 2.0 * ( common.S * ( c 1079 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1096 + common.ProjM 1080 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1097 if ( common.PZcms2 < 0.0 ) continue; 1081 if ( common.PZcms2 < 0.0 ) continue;
1098 1082
1099 common.PZcms = std::sqrt( common.PZcms2 ) 1083 common.PZcms = std::sqrt( common.PZcms2 );
1100 common.PMinusMin = std::sqrt( common.Proj 1084 common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
1101 common.PMinusMax = common.SqrtS - common. 1085 common.PMinusMax = common.SqrtS - common.TargMassT;
1102 common.TPlusMin = std::sqrt( common.TargM 1086 common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
1103 common.TPlusMax = common.SqrtS - common.P 1087 common.TPlusMax = common.SqrtS - common.ProjMassT;
1104 1088
1105 if ( G4UniformRand() <= 0.5 ) { // Rando 1089 if ( G4UniformRand() <= 0.5 ) { // Random choice projectile or target sampling
1106 if ( G4UniformRand() < theParameters->G 1090 if ( G4UniformRand() < theParameters->GetProbLogDistrPrD() ) {
1107 common.PMinusNew = ChooseP( common.PM 1091 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1108 } else { 1092 } else {
1109 common.PMinusNew = ( common.PMinusMax 1093 common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*G4UniformRand() + common.PMinusMin;
1110 } 1094 }
1111 if ( G4UniformRand() < theParameters->G 1095 if ( G4UniformRand() < theParameters->GetProbLogDistr() ) {
1112 common.TPlusNew = ChooseP( common.TPl 1096 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1113 } else { 1097 } else {
1114 common.TPlusNew = ( common.TPlusMax - 1098 common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*G4UniformRand() + common.TPlusMin;
1115 } 1099 }
1116 } else { 1100 } else {
1117 if ( G4UniformRand() < theParameters->G 1101 if ( G4UniformRand() < theParameters->GetProbLogDistr() ) {
1118 common.TPlusNew = ChooseP( common.TPl 1102 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1119 } else { 1103 } else {
1120 common.TPlusNew = ( common.TPlusMax - 1104 common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*G4UniformRand() + common.TPlusMin;
1121 } 1105 }
1122 if ( G4UniformRand() < theParameters->G 1106 if ( G4UniformRand() < theParameters->GetProbLogDistrPrD() ) {
1123 common.PMinusNew = ChooseP( common.PM 1107 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1124 } else { 1108 } else {
1125 common.PMinusNew = ( common.PMinusMax 1109 common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*G4UniformRand() + common.PMinusMin;
1126 } 1110 }
1127 } 1111 }
1128 1112
1129 common.Qminus = common.PMinusNew - common 1113 common.Qminus = common.PMinusNew - common.Pprojectile.minus();
1130 common.Qplus = -( common.TPlusNew - commo 1114 common.Qplus = -( common.TPlusNew - common.Ptarget.plus() );
1131 common.Qmomentum.setPz( (common.Qplus - c 1115 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1132 common.Qmomentum.setE( (common.Qplus + c 1116 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1133 #ifdef debugFTFexcitation << 1117 #ifdef debugFTFexictation
1134 G4cout <<"Sampled: Mpr, MdifPr, Mtr, Mdif 1118 G4cout <<"Sampled: Mpr, MdifPr, Mtr, MdifTr "<<G4endl
1135 << ( common.Pprojectile + common.Q 1119 << ( common.Pprojectile + common.Qmomentum ).mag() << " "
1136 << common.ProjectileNonDiffStateMi 1120 << common.ProjectileNonDiffStateMinMass << G4endl
1137 << ( common.Ptarget - common.Qmome 1121 << ( common.Ptarget - common.Qmomentum ).mag() << " "
1138 << common.TargetNonDiffStateMinMas 1122 << common.TargetNonDiffStateMinMass << G4endl;
1139 #endif 1123 #endif
1140 1124
1141 } while ( ( common.Pprojectile + common.Qmo 1125 } while ( ( common.Pprojectile + common.Qmomentum ).mag2() <
1142 common.ProjectileNonDiffStateMinM 1126 common.ProjectileNonDiffStateMinMass2 || // No double Diffraction
1143 ( common.Ptarget - common.Qmoment 1127 ( common.Ptarget - common.Qmomentum ).mag2() <
1144 common.TargetNonDiffStateMinMass2 1128 common.TargetNonDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
1145 1129
1146 projectile->SetStatus( 0 ); 1130 projectile->SetStatus( 0 );
1147 target->SetStatus( 0 ); 1131 target->SetStatus( 0 );
1148 return true; 1132 return true;
1149 } 1133 }
1150 1134
1151 1135
1152 //=========================================== 1136 //============================================================================
1153 1137
1154 void G4DiffractiveExcitation::CreateStrings( 1138 void G4DiffractiveExcitation::CreateStrings( G4VSplitableHadron* hadron,
1155 1139 G4bool isProjectile,
1156 1140 G4ExcitedString*& FirstString,
1157 1141 G4ExcitedString*& SecondString,
1158 1142 G4FTFParameters* theParameters ) const {
1159 1143
1160 //G4cout << "Create Strings SplitUp " << ha 1144 //G4cout << "Create Strings SplitUp " << hadron << G4endl
1161 // << "Defin " << hadron->GetDefiniti 1145 // << "Defin " << hadron->GetDefinition() << G4endl
1162 // << "Defin " << hadron->GetDefiniti 1146 // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
1163 1147
1164 G4bool HadronIsString = hadron->IsSplit(); << 1148 hadron->SplitUp();
1165 if( ! HadronIsString ) hadron->SplitUp(); <<
1166 1149
1167 G4Parton* start = hadron->GetNextParton(); 1150 G4Parton* start = hadron->GetNextParton();
1168 if ( start == nullptr ) { << 1151 if ( start == NULL ) {
1169 G4cout << " G4FTFModel::String() Error: N 1152 G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1170 FirstString = 0; SecondString = 0; 1153 FirstString = 0; SecondString = 0;
1171 return; 1154 return;
1172 } 1155 }
1173 1156
1174 G4Parton* end = hadron->GetNextParton(); 1157 G4Parton* end = hadron->GetNextParton();
1175 if ( end == nullptr ) { << 1158 if ( end == NULL ) {
1176 G4cout << " G4FTFModel::String() Error: N 1159 G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1177 FirstString = 0; SecondString = 0; 1160 FirstString = 0; SecondString = 0;
1178 return; 1161 return;
1179 } 1162 }
1180 1163
1181 //G4cout << start << " " << start->GetPDGco 1164 //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1182 // << G4endl 1165 // << G4endl
1183 // << "Create string " << start->GetP 1166 // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1184 <<
1185 if ( HadronIsString ) { <<
1186 if ( isProjectile ) { <<
1187 FirstString = new G4ExcitedString( end, <<
1188 } else { <<
1189 FirstString = new G4ExcitedString( end, <<
1190 } <<
1191 FirstString->SetTimeOfCreation( hadron->G <<
1192 FirstString->SetPosition( hadron->GetPosi <<
1193 SecondString = 0; <<
1194 return; <<
1195 } <<
1196 1167
1197 G4LorentzVector Phadron = hadron->Get4Momen 1168 G4LorentzVector Phadron = hadron->Get4Momentum();
1198 //G4cout << "String mom " << Phadron << G4e 1169 //G4cout << "String mom " << Phadron << G4endl;
1199 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 1170 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1200 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 ); 1171 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1201 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 ) 1172 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1202 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 1173 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1203 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 1174 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1204 1175
1205 G4int PDGcode_startQ = std::abs( start->Get 1176 G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1206 G4int PDGcode_endQ = std::abs( end->Get 1177 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1207 //G4cout << "PDGcode_startQ " << PDGcode_st 1178 //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1208 1179
1209 G4double Wmin( 0.0 ); 1180 G4double Wmin( 0.0 );
1210 if ( isProjectile ) { 1181 if ( isProjectile ) {
1211 Wmin = theParameters->GetProjMinDiffMass( 1182 Wmin = theParameters->GetProjMinDiffMass();
1212 } else { 1183 } else {
1213 Wmin = theParameters->GetTarMinDiffMass() 1184 Wmin = theParameters->GetTarMinDiffMass();
1214 } 1185 }
1215 1186
1216 G4double W = hadron->Get4Momentum().mag(); 1187 G4double W = hadron->Get4Momentum().mag();
1217 G4double W2 = W*W; 1188 G4double W2 = W*W;
1218 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); 1189 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); // x2( 0.0 )
1219 G4bool Kink = false; 1190 G4bool Kink = false;
1220 1191
1221 if ( ! ( ( start->GetDefinition()->GetParti 1192 if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1222 end->GetDefinition()->GetParti 1193 end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1223 ( start->GetDefinition()->GetParti 1194 ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1224 end->GetDefinition()->GetParti 1195 end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1225 // Kinky strings are allowed only for qq- 1196 // Kinky strings are allowed only for qq-q strings;
1226 // Kinky strings are impossible for other 1197 // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1227 // according to the analysis of Pbar P in 1198 // according to the analysis of Pbar P interactions
1228 1199
1229 if ( W > Wmin ) { // Kink is possible 1200 if ( W > Wmin ) { // Kink is possible
1230 if ( hadron->GetStatus() == 0 ) { 1201 if ( hadron->GetStatus() == 0 ) {
1231 G4double Pt2kink = theParameters->Get 1202 G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1232 if ( Pt2kink ) { 1203 if ( Pt2kink ) {
1233 Pt = std::sqrt( Pt2kink * ( G4Pow:: 1204 Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );
1234 } else { 1205 } else {
1235 Pt = 0.0; 1206 Pt = 0.0;
1236 } 1207 }
1237 } else { 1208 } else {
1238 Pt = 0.0; 1209 Pt = 0.0;
1239 } 1210 }
1240 1211
1241 if ( Pt > 500.0*MeV ) { 1212 if ( Pt > 500.0*MeV ) {
1242 G4double Ymax = G4Log( W/2.0/Pt + std 1213 G4double Ymax = G4Log( W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1243 G4double Y = Ymax*( 1.0 - 2.0*G4Unifo 1214 G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
1244 x1 = 1.0 - Pt/W * G4Exp( Y ); 1215 x1 = 1.0 - Pt/W * G4Exp( Y );
1245 x3 = 1.0 - Pt/W * G4Exp(-Y ); 1216 x3 = 1.0 - Pt/W * G4Exp(-Y );
1246 //x2 = 2.0 - x1 - x3; 1217 //x2 = 2.0 - x1 - x3;
1247 1218
1248 G4double Mass_startQ = 650.0*MeV; 1219 G4double Mass_startQ = 650.0*MeV;
1249 if ( PDGcode_startQ < 3 ) Mass_start 1220 if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV; // For constituent up or down quark
1250 if ( PDGcode_startQ == 3 ) Mass_start 1221 if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV; // For constituent strange quark
1251 if ( PDGcode_startQ == 4 ) Mass_start 1222 if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV; // For constituent charm quark
1252 if ( PDGcode_startQ == 5 ) Mass_start 1223 if ( PDGcode_startQ == 5 ) Mass_startQ = 4500.0*MeV; // For constituent bottom quark
1253 G4double Mass_endQ = 650.0*MeV; 1224 G4double Mass_endQ = 650.0*MeV;
1254 if ( PDGcode_endQ < 3 ) Mass_endQ = 1225 if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV; // For constituent up or down quark
1255 if ( PDGcode_endQ == 3 ) Mass_endQ = 1226 if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV; // For constituent strange quark
1256 if ( PDGcode_endQ == 4 ) Mass_endQ = 1227 if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV; // For constituent charm quark
1257 if ( PDGcode_endQ == 5 ) Mass_endQ = 1228 if ( PDGcode_endQ == 5 ) Mass_endQ = 4500.0*MeV; // For constituent bottom quark
1258 1229
1259 G4double P2_1 = W2*x1*x1/4.0 - Mass_e 1230 G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1260 G4double P2_3 = W2*x3*x3/4.0 - Mass_s 1231 G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1261 G4double P2_2 = sqr( (2.0 - x1 - x3)* 1232 G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
1262 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) { 1233 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1263 Kink = false; 1234 Kink = false;
1264 } else { 1235 } else {
1265 G4double P_1 = std::sqrt( P2_1 ); 1236 G4double P_1 = std::sqrt( P2_1 );
1266 G4double P_2 = std::sqrt( P2_2 ); 1237 G4double P_2 = std::sqrt( P2_2 );
1267 G4double P_3 = std::sqrt( P2_3 ); 1238 G4double P_3 = std::sqrt( P2_3 );
1268 G4double CosT12 = ( P2_3 - P2_1 - P 1239 G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1269 G4double CosT13 = ( P2_2 - P2_1 - P 1240 G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1270 1241
1271 if ( std::abs( CosT12 ) > 1.0 || 1242 if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1272 Kink = false; 1243 Kink = false;
1273 } else { 1244 } else {
1274 Kink = true; 1245 Kink = true;
1275 Pt = P_2 * std::sqrt( 1.0 - CosT1 1246 Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated
1276 Pstart.setPx( -Pt ); Pstart.setPy 1247 Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1277 Pend.setPx( 0.0 ); Pend.setPy( 1248 Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1278 Pkink.setPx( Pt ); Pkink.setPy( 1249 Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1279 Pstart.setE( x3*W/2.0 ); 1250 Pstart.setE( x3*W/2.0 );
1280 Pkink.setE( Pkink.vect().mag() ); 1251 Pkink.setE( Pkink.vect().mag() );
1281 Pend.setE( x1*W/2.0 ); 1252 Pend.setE( x1*W/2.0 );
1282 1253
1283 G4double XkQ = GetQuarkFractionOf 1254 G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1284 if ( Pkink.getZ() > 0.0 ) { 1255 if ( Pkink.getZ() > 0.0 ) {
1285 if ( XkQ > 0.5 ) { 1256 if ( XkQ > 0.5 ) {
1286 PkinkQ1 = XkQ*Pkink; 1257 PkinkQ1 = XkQ*Pkink;
1287 } else { 1258 } else {
1288 PkinkQ1 = (1.0 - XkQ)*Pkink; 1259 PkinkQ1 = (1.0 - XkQ)*Pkink;
1289 } 1260 }
1290 } else { 1261 } else {
1291 if ( XkQ > 0.5 ) { 1262 if ( XkQ > 0.5 ) {
1292 PkinkQ1 = (1.0 - XkQ)*Pkink; 1263 PkinkQ1 = (1.0 - XkQ)*Pkink;
1293 } else { 1264 } else {
1294 PkinkQ1 = XkQ*Pkink; 1265 PkinkQ1 = XkQ*Pkink;
1295 } 1266 }
1296 } 1267 }
1297 1268
1298 PkinkQ2 = Pkink - PkinkQ1; 1269 PkinkQ2 = Pkink - PkinkQ1;
1299 // Minimizing Pt1^2+Pt3^2 1270 // Minimizing Pt1^2+Pt3^2
1300 G4double Cos2Psi = ( sqr(x1) - sq 1271 G4double Cos2Psi = ( sqr(x1) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1301 std::sqrt( sqr 1272 std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
1302 G4double Psi = std::acos( Cos2Psi 1273 G4double Psi = std::acos( Cos2Psi );
1303 1274
1304 G4LorentzRotation Rotate; 1275 G4LorentzRotation Rotate;
1305 if ( isProjectile ) { 1276 if ( isProjectile ) {
1306 Rotate.rotateY( Psi ); 1277 Rotate.rotateY( Psi );
1307 } else { 1278 } else {
1308 Rotate.rotateY( pi + Psi ); 1279 Rotate.rotateY( pi + Psi );
1309 } 1280 }
1310 Rotate.rotateZ( twopi * G4Uniform 1281 Rotate.rotateZ( twopi * G4UniformRand() );
1311 Pstart *= Rotate; 1282 Pstart *= Rotate;
1312 Pkink *= Rotate; 1283 Pkink *= Rotate;
1313 PkinkQ1 *= Rotate; 1284 PkinkQ1 *= Rotate;
1314 PkinkQ2 *= Rotate; 1285 PkinkQ2 *= Rotate;
1315 Pend *= Rotate; 1286 Pend *= Rotate;
1316 } 1287 }
1317 } // End of if ( P2_1 <= 0.0 || P2 1288 } // End of if ( P2_1 <= 0.0 || P2_3 <= 0.0 )
1318 } // End of if ( Pt > 500.0*MeV ) 1289 } // End of if ( Pt > 500.0*MeV )
1319 } // End of if ( W > Wmin ) : check for a 1290 } // End of if ( W > Wmin ) : check for a kink
1320 } // end of qq-q string selection 1291 } // end of qq-q string selection
1321 1292
1322 if ( Kink ) { // Kink is possible 1293 if ( Kink ) { // Kink is possible
1323 1294
1324 //G4cout << "Kink is sampled!" << G4endl; 1295 //G4cout << "Kink is sampled!" << G4endl;
1325 std::vector< G4double > QuarkProbabilitie 1296 std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1326 theParameters->GetQuarkProbabilitiesA 1297 theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1327 1298
1328 G4int QuarkInGluon( 1 ); G4double Ksi = G 1299 G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1329 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) 1300 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1330 //G4cout << "Iq " << Iq << G4endl; 1301 //G4cout << "Iq " << Iq << G4endl;
1331 if ( Ksi > QuarkProbabilitiesAtGluonSpl 1302 if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1332 } 1303 }
1333 //G4cout << "Last Iq " << QuarkInGluon << 1304 //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1334 G4Parton* Gquark = new G4Parton( QuarkInG 1305 G4Parton* Gquark = new G4Parton( QuarkInGluon );
1335 G4Parton* Ganti_quark = new G4Parton( -Qu 1306 G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1336 //G4cout << "Lorentz " << G4endl; 1307 //G4cout << "Lorentz " << G4endl;
1337 1308
1338 G4LorentzRotation toCMS( -1 * Phadron.boo 1309 G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1339 G4LorentzRotation toLab( toCMS.inverse() 1310 G4LorentzRotation toLab( toCMS.inverse() );
1340 //G4cout << "Pstart " << Pstart << G4endl 1311 //G4cout << "Pstart " << Pstart << G4endl;
1341 //G4cout << "Pend " << Pend << G4endl; 1312 //G4cout << "Pend " << Pend << G4endl;
1342 //G4cout << "Kink1 " <<PkinkQ1<<G4endl; 1313 //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1343 //G4cout << "Kink2 " <<PkinkQ2<<G4endl; 1314 //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1344 //G4cout << "Pstart " << Pstart << G4endl 1315 //G4cout << "Pstart " << Pstart << G4endl<<G4endl;
1345 1316
1346 Pstart.transform( toLab ); start->Set4Mo 1317 Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1347 PkinkQ1.transform( toLab ); 1318 PkinkQ1.transform( toLab );
1348 PkinkQ2.transform( toLab ); 1319 PkinkQ2.transform( toLab );
1349 Pend.transform( toLab ); end->Set4Mome 1320 Pend.transform( toLab ); end->Set4Momentum( Pend );
1350 //G4cout << "Pstart " << Pstart << G4endl 1321 //G4cout << "Pstart " << Pstart << G4endl;
1351 //G4cout << "Pend " << Pend << G4endl; 1322 //G4cout << "Pend " << Pend << G4endl;
1352 //G4cout << "Defin " << hadron->GetDefini 1323 //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1353 //G4cout << "Defin " << hadron->GetDefini 1324 //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1354 1325
1355 //G4int absPDGcode = std::abs( hadron->Ge 1326 //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1356 G4int absPDGcode = 1500; 1327 G4int absPDGcode = 1500;
1357 if ( start->GetDefinition()->GetParticleS 1328 if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1358 end->GetDefinition()->GetParticleSub 1329 end->GetDefinition()->GetParticleSubType() == "quark" ) {
1359 absPDGcode = 110; 1330 absPDGcode = 110;
1360 } 1331 }
1361 //G4cout << "absPDGcode " << absPDGcode < 1332 //G4cout << "absPDGcode " << absPDGcode << G4endl;
1362 1333
1363 if ( absPDGcode < 1000 ) { // meson 1334 if ( absPDGcode < 1000 ) { // meson
1364 if ( isProjectile ) { // Projectile 1335 if ( isProjectile ) { // Projectile
1365 if ( end->GetDefinition()->GetPDGEnco 1336 if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1366 FirstString = new G4ExcitedString( 1337 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1367 SecondString = new G4ExcitedString( 1338 SecondString = new G4ExcitedString( Gquark, start , +1 );
1368 Ganti_quark->Set4Momentum( PkinkQ1 1339 Ganti_quark->Set4Momentum( PkinkQ1 );
1369 Gquark->Set4Momentum( PkinkQ2 ); 1340 Gquark->Set4Momentum( PkinkQ2 );
1370 } else { // Anti_Quark on the end 1341 } else { // Anti_Quark on the end
1371 FirstString = new G4ExcitedString( 1342 FirstString = new G4ExcitedString( end , Gquark, +1 );
1372 SecondString = new G4ExcitedString( 1343 SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1373 Gquark->Set4Momentum( PkinkQ1 ); 1344 Gquark->Set4Momentum( PkinkQ1 );
1374 Ganti_quark->Set4Momentum( PkinkQ2 1345 Ganti_quark->Set4Momentum( PkinkQ2 );
1375 } 1346 }
1376 } else { // Target 1347 } else { // Target
1377 if ( end->GetDefinition()->GetPDGEnco 1348 if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1378 FirstString = new G4ExcitedString( 1349 FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1379 SecondString = new G4ExcitedString( 1350 SecondString = new G4ExcitedString( start , Gquark, -1 );
1380 Ganti_quark->Set4Momentum( PkinkQ2 1351 Ganti_quark->Set4Momentum( PkinkQ2 );
1381 Gquark->Set4Momentum( PkinkQ1 ); 1352 Gquark->Set4Momentum( PkinkQ1 );
1382 } else { // Anti_Quark on the end 1353 } else { // Anti_Quark on the end
1383 FirstString = new G4ExcitedString( 1354 FirstString = new G4ExcitedString( Gquark, end , -1 );
1384 SecondString = new G4ExcitedString( 1355 SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1385 Gquark->Set4Momentum( PkinkQ2 ); 1356 Gquark->Set4Momentum( PkinkQ2 );
1386 Ganti_quark->Set4Momentum( PkinkQ1 1357 Ganti_quark->Set4Momentum( PkinkQ1 );
1387 } 1358 }
1388 } 1359 }
1389 } else { // Baryon/AntiBaryon 1360 } else { // Baryon/AntiBaryon
1390 if ( isProjectile ) { // Projectile 1361 if ( isProjectile ) { // Projectile
1391 if ( end->GetDefinition()->GetParticl 1362 if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1392 end->GetDefinition()->GetPDGEnco 1363 end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1393 FirstString = new G4ExcitedString( 1364 FirstString = new G4ExcitedString( end , Gquark, +1 );
1394 SecondString = new G4ExcitedString( 1365 SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1395 Gquark->Set4Momentum( PkinkQ1 ); 1366 Gquark->Set4Momentum( PkinkQ1 );
1396 Ganti_quark->Set4Momentum( PkinkQ2 1367 Ganti_quark->Set4Momentum( PkinkQ2 );
1397 } else { / 1368 } else { // Anti_DiQuark on the end or quark
1398 FirstString = new G4ExcitedString( 1369 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1399 SecondString = new G4ExcitedString( 1370 SecondString = new G4ExcitedString( Gquark, start , +1 );
1400 Ganti_quark->Set4Momentum( PkinkQ1 1371 Ganti_quark->Set4Momentum( PkinkQ1 );
1401 Gquark->Set4Momentum( PkinkQ2 ); 1372 Gquark->Set4Momentum( PkinkQ2 );
1402 } 1373 }
1403 } else { // Target 1374 } else { // Target
1404 if ( end->GetDefinition()->GetParticl 1375 if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1405 end->GetDefinition()->GetPDGEnco 1376 end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1406 Gquark->Set4Momentum( PkinkQ1 ); 1377 Gquark->Set4Momentum( PkinkQ1 );
1407 Ganti_quark->Set4Momentum( PkinkQ2 1378 Ganti_quark->Set4Momentum( PkinkQ2 );
1408 FirstString = new G4ExcitedString( 1379 FirstString = new G4ExcitedString( end, Gquark, -1 );
1409 SecondString = new G4ExcitedString( 1380 SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1410 } else { // Anti_DiQuark on the end 1381 } else { // Anti_DiQuark on the end or Q
1411 FirstString = new G4ExcitedString( 1382 FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1412 SecondString = new G4ExcitedString( 1383 SecondString = new G4ExcitedString( start , Gquark, -1 );
1413 Gquark->Set4Momentum( PkinkQ2 ); 1384 Gquark->Set4Momentum( PkinkQ2 );
1414 Ganti_quark->Set4Momentum( PkinkQ1 1385 Ganti_quark->Set4Momentum( PkinkQ1 );
1415 } 1386 }
1416 } 1387 }
1417 } 1388 }
1418 1389
1419 FirstString->SetTimeOfCreation( hadron->G 1390 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1420 FirstString->SetPosition( hadron->GetPosi 1391 FirstString->SetPosition( hadron->GetPosition() );
1421 SecondString->SetTimeOfCreation( hadron-> 1392 SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1422 SecondString->SetPosition( hadron->GetPos 1393 SecondString->SetPosition( hadron->GetPosition() );
1423 1394
1424 } else { // Kink is impossible 1395 } else { // Kink is impossible
1425 1396
1426 if ( isProjectile ) { 1397 if ( isProjectile ) {
1427 FirstString = new G4ExcitedString( end, 1398 FirstString = new G4ExcitedString( end, start, +1 );
1428 } else { 1399 } else {
1429 FirstString = new G4ExcitedString( end, 1400 FirstString = new G4ExcitedString( end, start, -1 );
1430 } 1401 }
1431 FirstString->SetTimeOfCreation( hadron->G 1402 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1432 FirstString->SetPosition( hadron->GetPosi 1403 FirstString->SetPosition( hadron->GetPosition() );
1433 SecondString = 0; 1404 SecondString = 0;
1434 // momenta of string ends 1405 // momenta of string ends
1435 G4LorentzVector HadronMom = hadron->Get4M 1406 G4LorentzVector HadronMom = hadron->Get4Momentum();
1436 G4LorentzVector Pstart1 = G4LorentzVector 1407 G4LorentzVector Pstart1 = G4LorentzVector( HadronMom.px()/2.0, HadronMom.py()/2.0, 0.0, 0.0 ); // Quark momentum
1437 G4LorentzVector Pend1 = G4LorentzVector 1408 G4LorentzVector Pend1 = G4LorentzVector( HadronMom.px()/2.0, HadronMom.py()/2.0, 0.0, 0.0 ); // Di-quark momentum
1438 G4double Pz = HadronMom.pz(); 1409 G4double Pz = HadronMom.pz();
1439 G4double Eh = HadronMom.e(); 1410 G4double Eh = HadronMom.e();
1440 G4double Pt2 = sqr( HadronMom.px() ) + sq 1411 G4double Pt2 = sqr( HadronMom.px() ) + sqr( HadronMom.py() );
1441 G4double Mt2 = HadronMom.mt2(); 1412 G4double Mt2 = HadronMom.mt2();
1442 G4double Exp = std::sqrt( sqr(Pz) + ( sqr 1413 G4double Exp = std::sqrt( sqr(Pz) + ( sqr(Mt2) - 4.0*sqr(Eh)*Pt2/4.0 )/Mt2 )/2.0;
1443 G4double Pzq = Pz/2.0 - Exp; 1414 G4double Pzq = Pz/2.0 - Exp; Pstart1.setZ( Pzq );
1444 G4double Eq = std::sqrt( sqr(Pzq) + Pt2 1415 G4double Eq = std::sqrt( sqr(Pzq) + Pt2/4.0 ); Pstart1.setE( Eq );
1445 G4double Pzqq = Pz/2.0 + Exp; 1416 G4double Pzqq = Pz/2.0 + Exp; Pend1.setZ(Pzqq);
1446 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt 1417 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt2/4.0 ); Pend1.setE(Eqq);
1447 start->Set4Momentum( Pstart1 ); 1418 start->Set4Momentum( Pstart1 );
1448 end->Set4Momentum( Pend1 ); 1419 end->Set4Momentum( Pend1 );
1449 Pstart = Pstart1; Pend = Pend1; 1420 Pstart = Pstart1; Pend = Pend1;
1450 1421
1451 } // End of "if (Kink)" 1422 } // End of "if (Kink)"
1452 1423
1453 //G4cout << "Quarks in the string at creati 1424 //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1454 // << " " << FirstString->GetLeftPart 1425 // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1455 // << FirstString << " " << SecondStr 1426 // << FirstString << " " << SecondString << G4endl;
1456 1427
1457 #ifdef G4_FTFDEBUG 1428 #ifdef G4_FTFDEBUG
1458 G4cout << " generated string flavors 1429 G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1459 << end->GetPDGcode() << G4endl << " 1430 << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1460 << start->Get4Momentum() << "mass : 1431 << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1461 << " generated string momenta: Diqua 1432 << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1462 << end->Get4Momentum().mag() << G4en 1433 << end->Get4Momentum().mag() << G4endl << " sum of ends "
1463 << Pstart + Pend << G4endl << " Orig 1434 << Pstart + Pend << G4endl << " Original "
1464 << hadron->Get4Momentum() << " "<<ha 1435 << hadron->Get4Momentum() << " "<<hadron->Get4Momentum().mag() << G4endl;
1465 #endif 1436 #endif
1466 1437
1467 return; 1438 return;
1468 } 1439 }
1469 1440
1470 1441
1471 //=========================================== 1442 //============================================================================
1472 1443
1473 G4double G4DiffractiveExcitation::ChooseP( G4 1444 G4double G4DiffractiveExcitation::ChooseP( G4double Pmin, G4double Pmax ) const {
1474 // Choose an x between Xmin and Xmax with P 1445 // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1475 // To be improved... 1446 // To be improved...
1476 G4double range = Pmax - Pmin; 1447 G4double range = Pmax - Pmin;
1477 if ( Pmin <= 0.0 || range <= 0.0 ) { 1448 if ( Pmin <= 0.0 || range <= 0.0 ) {
1478 G4cout << " Pmin, range : " << Pmin << " 1449 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1479 throw G4HadronicException( __FILE__, __LI 1450 throw G4HadronicException( __FILE__, __LINE__,
1480 "G4Diffractive 1451 "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1481 } 1452 }
1482 G4double P = Pmin * G4Pow::GetInstance()->p 1453 G4double P = Pmin * G4Pow::GetInstance()->powA( Pmax/Pmin, G4UniformRand() );
1483 //G4double P = (Pmax - Pmin) * G4UniformRan 1454 //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1484 return P; 1455 return P;
1485 } 1456 }
1486 1457
1487 1458
1488 //=========================================== 1459 //============================================================================
1489 1460
1490 G4ThreeVector G4DiffractiveExcitation::Gaussi 1461 G4ThreeVector G4DiffractiveExcitation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1491 // @@ this method is used in FTFModel as w 1462 // @@ this method is used in FTFModel as well. Should go somewhere common!
1492 G4double Pt2( 0.0 ); 1463 G4double Pt2( 0.0 );
1493 if ( AveragePt2 <= 0.0 ) { 1464 if ( AveragePt2 <= 0.0 ) {
1494 Pt2 = 0.0; 1465 Pt2 = 0.0;
1495 } else { 1466 } else {
1496 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor 1467 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1497 ( G4Ex 1468 ( G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1498 } 1469 }
1499 G4double Pt = std::sqrt( Pt2 ); 1470 G4double Pt = std::sqrt( Pt2 );
1500 G4double phi = G4UniformRand() * twopi; 1471 G4double phi = G4UniformRand() * twopi;
1501 return G4ThreeVector( Pt * std::cos( phi ), 1472 return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1502 } 1473 }
1503 1474
1504 1475
1505 //=========================================== 1476 //============================================================================
1506 1477
1507 G4double G4DiffractiveExcitation::GetQuarkFra 1478 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink( G4double zmin, G4double zmax ) const {
1508 G4double z, yf; 1479 G4double z, yf;
1509 const G4int maxNumberOfLoops = 10000; 1480 const G4int maxNumberOfLoops = 10000;
1510 G4int loopCounter = 0; 1481 G4int loopCounter = 0;
1511 do { 1482 do {
1512 z = zmin + G4UniformRand() * (zmax - zmin 1483 z = zmin + G4UniformRand() * (zmax - zmin);
1513 yf = z*z + sqr(1.0 - z); 1484 yf = z*z + sqr(1.0 - z);
1514 } while ( ( G4UniformRand() > yf ) && 1485 } while ( ( G4UniformRand() > yf ) &&
1515 ++loopCounter < maxNumberOfLoops 1486 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1516 if ( loopCounter >= maxNumberOfLoops ) { 1487 if ( loopCounter >= maxNumberOfLoops ) {
1517 z = 0.5*(zmin + zmax); // Just something 1488 z = 0.5*(zmin + zmax); // Just something acceptable, without any physics consideration.
1518 } 1489 }
1519 return z; 1490 return z;
1520 } 1491 }
1521 1492
1522 1493
1523 //=========================================== 1494 //============================================================================
1524 1495
1525 void G4DiffractiveExcitation::UnpackMeson( co 1496 void G4DiffractiveExcitation::UnpackMeson( const G4int IdPDG, G4int& Q1, G4int& Q2 ) const {
1526 G4int absIdPDG = std::abs( IdPDG ); 1497 G4int absIdPDG = std::abs( IdPDG );
1527 if ( ! ( absIdPDG == 111 || absIdPDG == 221 1498 if ( ! ( absIdPDG == 111 || absIdPDG == 221 || absIdPDG == 331 || // Pi0 , Eta , Eta'
1528 absIdPDG == 441 || absIdPDG == 443 1499 absIdPDG == 441 || absIdPDG == 443 || absIdPDG == 553 ) ) { // Etac , J/psi , Upsilon
1529 // All other projectile mesons (including 1500 // All other projectile mesons (including charmed and bottom ones)
1530 Q1 = absIdPDG / 100; 1501 Q1 = absIdPDG / 100;
1531 Q2 = (absIdPDG % 100) / 10; 1502 Q2 = (absIdPDG % 100) / 10;
1532 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) 1503 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1533 if ( IdPDG < 0 ) anti *= -1; 1504 if ( IdPDG < 0 ) anti *= -1;
1534 Q1 *= anti; 1505 Q1 *= anti;
1535 Q2 *= -1 * anti; 1506 Q2 *= -1 * anti;
1536 } else { 1507 } else {
1537 if ( absIdPDG == 441 || absIdPDG == 443 ) 1508 if ( absIdPDG == 441 || absIdPDG == 443 ) { // Etac , J/psi
1538 Q1 = 4; Q2 = -4; 1509 Q1 = 4; Q2 = -4;
1539 } else if ( absIdPDG == 553 ) { 1510 } else if ( absIdPDG == 553 ) { // Upsilon
1540 Q1 = 5; Q2 = -5; 1511 Q1 = 5; Q2 = -5;
1541 } else { 1512 } else { // Pi0 , Eta , Eta'
1542 if ( G4UniformRand() < 0.5 ) { 1513 if ( G4UniformRand() < 0.5 ) {
1543 Q1 = 1; Q2 = -1; 1514 Q1 = 1; Q2 = -1;
1544 } else { 1515 } else {
1545 Q1 = 2; Q2 = -2; 1516 Q1 = 2; Q2 = -2;
1546 } 1517 }
1547 } 1518 }
1548 } 1519 }
1549 return; 1520 return;
1550 } 1521 }
1551 1522
1552 1523
1553 //=========================================== 1524 //============================================================================
1554 1525
1555 void G4DiffractiveExcitation::UnpackBaryon( G 1526 void G4DiffractiveExcitation::UnpackBaryon( G4int IdPDG,
1556 G 1527 G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1557 Q1 = IdPDG / 1000; 1528 Q1 = IdPDG / 1000;
1558 Q2 = (IdPDG % 1000) / 100; 1529 Q2 = (IdPDG % 1000) / 100;
1559 Q3 = (IdPDG % 100) / 10; 1530 Q3 = (IdPDG % 100) / 10;
1560 return; 1531 return;
1561 } 1532 }
1562 1533
1563 1534
1564 //=========================================== 1535 //============================================================================
1565 1536
1566 G4int G4DiffractiveExcitation::NewNucleonId( 1537 G4int G4DiffractiveExcitation::NewNucleonId( G4int Q1, G4int Q2, G4int Q3 ) const {
1567 // Order the three integers in such a way t 1538 // Order the three integers in such a way that Q1 >= Q2 >= Q3
1568 G4int TmpQ( 0 ); 1539 G4int TmpQ( 0 );
1569 if ( Q3 > Q2 ) { 1540 if ( Q3 > Q2 ) {
1570 TmpQ = Q2; 1541 TmpQ = Q2;
1571 Q2 = Q3; 1542 Q2 = Q3;
1572 Q3 = TmpQ; 1543 Q3 = TmpQ;
1573 } else if ( Q3 > Q1 ) { 1544 } else if ( Q3 > Q1 ) {
1574 TmpQ = Q1; 1545 TmpQ = Q1;
1575 Q1 = Q3; 1546 Q1 = Q3;
1576 Q3 = TmpQ; 1547 Q3 = TmpQ;
1577 } 1548 }
1578 if ( Q2 > Q1 ) { 1549 if ( Q2 > Q1 ) {
1579 TmpQ = Q1; 1550 TmpQ = Q1;
1580 Q1 = Q2; 1551 Q1 = Q2;
1581 Q2 = TmpQ; 1552 Q2 = TmpQ;
1582 } 1553 }
1583 // By now Q1 >= Q2 >= Q3 1554 // By now Q1 >= Q2 >= Q3
1584 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 1555 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1585 return NewCode; 1556 return NewCode;
1586 } 1557 }
1587 1558
1588 1559
1589 //=========================================== 1560 //============================================================================
1590 1561
1591 G4DiffractiveExcitation::G4DiffractiveExcitat 1562 G4DiffractiveExcitation::G4DiffractiveExcitation( const G4DiffractiveExcitation& ) {
1592 throw G4HadronicException( __FILE__, __LINE 1563 throw G4HadronicException( __FILE__, __LINE__,
1593 "G4DiffractiveEx 1564 "G4DiffractiveExcitation copy constructor not meant to be called" );
1594 } 1565 }
1595 1566
1596 1567
1597 //=========================================== 1568 //============================================================================
1598 1569
1599 const G4DiffractiveExcitation & G4Diffractive 1570 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=( const G4DiffractiveExcitation& ) {
1600 throw G4HadronicException( __FILE__, __LINE 1571 throw G4HadronicException( __FILE__, __LINE__,
1601 "G4DiffractiveEx 1572 "G4DiffractiveExcitation = operator not meant to be called" );
1602 return *this; 1573 return *this;
1603 } 1574 }
1604 1575
1605 1576
1606 //=========================================== 1577 //============================================================================
1607 1578
1608 G4bool G4DiffractiveExcitation::operator==( c 1579 G4bool G4DiffractiveExcitation::operator==( const G4DiffractiveExcitation& ) const {
1609 throw G4HadronicException( __FILE__, __LINE 1580 throw G4HadronicException( __FILE__, __LINE__,
1610 "G4DiffractiveEx 1581 "G4DiffractiveExcitation == operator not meant to be called" );
1611 } 1582 }
1612 1583
1613 1584
1614 //=========================================== 1585 //============================================================================
1615 1586
1616 G4bool G4DiffractiveExcitation::operator!= ( 1587 G4bool G4DiffractiveExcitation::operator!= ( const G4DiffractiveExcitation& ) const {
1617 throw G4HadronicException( __FILE__, __LINE 1588 throw G4HadronicException( __FILE__, __LINE__,
1618 "G4DiffractiveEx 1589 "G4DiffractiveExcitation != operator not meant to be called" );
1619 } 1590 }
1620 1591
1621 1592