Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DiffractiveExcitation.cc
Go to the documentation of this file.
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"
52#include "G4SystemOfUnits.hh"
53
55#include "G4FTFParameters.hh"
57
58#include "G4RotationMatrix.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
80
81
82//============================================================================
83
85
86
87//============================================================================
88
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
325G4int G4DiffractiveExcitation::
326ExciteParticipants_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
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 ) );
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
944G4bool G4DiffractiveExcitation::
945ExciteParticipants_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
1052G4bool G4DiffractiveExcitation::
1053ExciteParticipants_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
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
1473G4double 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
1490G4ThreeVector 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
1507G4double 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
1525void 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
1555void 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
1566G4int 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
1592 throw G4HadronicException( __FILE__, __LINE__,
1593 "G4DiffractiveExcitation copy constructor not meant to be called" );
1594}
1595
1596
1597//============================================================================
1598
1599const 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
1608G4bool G4DiffractiveExcitation::operator==( const G4DiffractiveExcitation& ) const {
1609 throw G4HadronicException( __FILE__, __LINE__,
1610 "G4DiffractiveExcitation == operator not meant to be called" );
1611}
1612
1613
1614//============================================================================
1615
1616G4bool G4DiffractiveExcitation::operator!= ( const G4DiffractiveExcitation& ) const {
1617 throw G4HadronicException( __FILE__, __LINE__,
1618 "G4DiffractiveExcitation != operator not meant to be called" );
1619}
1620
G4double Y(G4double density)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double G4Log(G4double x)
Definition: G4Log.hh:227
CLHEP::HepLorentzRotation G4LorentzRotation
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:83
long G4long
Definition: G4Types.hh:87
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
double mag2() const
double mag() const
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
double theta() const
Hep3Vector boostVector() const
Hep3Vector vect() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
void SetPosition(const G4ThreeVector &aPosition)
void SetTimeOfCreation(G4double aTime)
G4double GetProbLogDistrPrD()
G4double GetAveragePt2()
G4double GetProbLogDistr()
G4double GetProjMinNonDiffMass()
G4double GetAvaragePt2ofElasticScattering()
G4double GetTarMinDiffMass()
G4double GetTarMinNonDiffMass()
G4double GetProcProb(const G4int ProcN, const G4double y)
G4double GetPt2Kink()
G4double GetDeltaProbAtQuarkExchange()
G4double GetProjMinDiffMass()
G4double GetProbOfSameQuarkExchange()
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
const G4String & GetParticleType() const
G4double GetPDGWidth() const
const G4String & GetParticleName() const
const G4String & GetParticleSubType() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * GetDefinition()
Definition: G4Parton.hh:161
G4int GetPDGcode() const
Definition: G4Parton.hh:127
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
void SetStatus(const G4int aStatus)
virtual void SplitUp()=0
const G4ParticleDefinition * GetDefinition() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
virtual G4Parton * GetNextParton()=0
void SetDefinition(const G4ParticleDefinition *aDefinition)
const G4LorentzVector & Get4Momentum() const
const G4ThreeVector & GetPosition() const
void IncrementCollisionCount(G4int aCount)
#define W
Definition: crc32.c:84
T sqr(const T &x)
Definition: templates.hh:128