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