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