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