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