ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DiffractiveExcitation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DiffractiveExcitation.cc
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 
76 
77 //============================================================================
78 
79 G4DiffractiveExcitation::G4DiffractiveExcitation() {}
80 
81 
82 //============================================================================
83 
84 G4DiffractiveExcitation::~G4DiffractiveExcitation() {}
85 
86 
87 //============================================================================
88 
89 G4bool G4DiffractiveExcitation::ExciteParticipants( G4VSplitableHadron* projectile,
90  G4VSplitableHadron* target,
91  G4FTFParameters* theParameters,
92  G4ElasticHNScattering* theElastic ) const {
93 
94  #ifdef 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 
325 G4int G4DiffractiveExcitation::
326 ExciteParticipants_doChargeExchange( G4VSplitableHadron* projectile,
327  G4VSplitableHadron* target,
328  G4FTFParameters* theParameters,
329  G4ElasticHNScattering* theElastic,
330  G4DiffractiveExcitation::CommonVariables& common ) const {
331  // First of the three utility methods used only by ExciteParticipants:
332  // it does the sampling for the charge-exchange case.
333  // This method returns an integer code - instead of a boolean, with the following meaning:
334  // "0" : successfully ended and nothing else needs to be done;
335  // "1" : successfully completed, but the work needs to be continued;
336  // "99" : unsuccessfully ended, nothing else can be done.
337  G4int returnCode = 99;
338 
339  G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
340  G4ParticleDefinition* TestParticle = 0;
341  G4double MtestPr = 0.0, MtestTr = 0.0;
342 
343  #ifdef debugFTFexictation
344  G4cout << "Q exchange --------------------------" << G4endl;
345  #endif
346  G4int NewProjCode = 0, NewTargCode = 0, ProjQ1 = 0, ProjQ2 = 0, ProjQ3 = 0;
347  // Projectile unpacking
348  if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
349  UnpackMeson( common.ProjectilePDGcode, ProjQ1, ProjQ2 );
350  } else { // projectile is baryon
351  UnpackBaryon( common.ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
352  }
353  // Target unpacking
354  G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;
355  UnpackBaryon( common.TargetPDGcode, TargQ1, TargQ2, TargQ3 );
356  #ifdef debugFTFexictation
357  G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
358  << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
359  #endif
360 
361  // Sampling of exchanged quarks
362  G4int ProjExchangeQ = 0, TargExchangeQ = 0;
363  if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
364 
365  G4bool isProjQ1Quark = false;
366  ProjExchangeQ = ProjQ2;
367  if ( ProjQ1 > 0 ) { // ProjQ1 is a quark
368  isProjQ1Quark = true;
369  ProjExchangeQ = ProjQ1;
370  }
371  // Exchange of non-identical quarks is allowed
372  G4int NpossibleStates = 0;
373  if ( ProjExchangeQ != TargQ1 ) NpossibleStates++;
374  if ( ProjExchangeQ != TargQ2 ) NpossibleStates++;
375  if ( ProjExchangeQ != TargQ3 ) NpossibleStates++;
376  G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
377  NpossibleStates = 0;
378  if ( ProjExchangeQ != TargQ1 ) {
379  if ( ++NpossibleStates == Nsampled ) {
380  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ;
381  isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
382  }
383  }
384  if ( ProjExchangeQ != TargQ2 ) {
385  if ( ++NpossibleStates == Nsampled ) {
386  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ;
387  isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
388  }
389  }
390  if ( ProjExchangeQ != TargQ3 ) {
391  if ( ++NpossibleStates == Nsampled ) {
392  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ;
393  isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
394  }
395  }
396  #ifdef debugFTFexictation
397  G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
398  #endif
399 
400  G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ2 = std::abs( ProjQ2 );
401  G4bool ProjExcited = false;
402  const G4int maxNumberOfAttempts = 50;
403  G4int attempts = 0;
404  while ( attempts++ < maxNumberOfAttempts ) { /* Loop checking, 10.08.2015, A.Ribon */
405 
406  // Determination of a new projectile ID which satisfies energy-momentum conservation
407  G4double Ksi = G4UniformRand();
408  if ( aProjQ1 == aProjQ2 ) {
409  if ( aProjQ1 != 3 ) {
410  NewProjCode = 111; // Pi0-meson
411  if ( Ksi < 0.5 ) {
412  NewProjCode = 221; // Eta -meson
413  if ( Ksi < 0.25 ) {
414  NewProjCode = 331; // Eta'-meson
415  }
416  }
417  } else {
418  NewProjCode = 221; // Eta -meson
419  if ( Ksi < 0.5 ) {
420  NewProjCode = 331; // Eta'-meson
421  }
422  }
423  } else {
424  if ( aProjQ1 > aProjQ2 ) {
425  NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
426  } else {
427  NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
428  }
429  }
430  #ifdef debugFTFexictation
431  G4cout << "NewProjCode " << NewProjCode << G4endl;
432  #endif
433  ProjExcited = false;
434  if ( G4UniformRand() < 0.5 ) {
435  NewProjCode += 2; // Excited meson (J=1 -> 2*J+1=2+1=3, last digit of the PDG code)
436  ProjExcited = true;
437  }
438  // So far we have used the absolute values of the PDG codes of the two constituent quarks:
439  // now look at their signed values to set properly the signed of the meson's PDG code.
440  G4int value = ProjQ1, absValue = aProjQ1, Qquarks = 0;
441  for ( G4int iQuark = 0; iQuark < 2; iQuark++ ) {
442  if ( iQuark == 1 ) {
443  value = ProjQ2; absValue = aProjQ2;
444  }
445  if ( absValue == 2 ) {
446  Qquarks += value; // u or ubar : u-quark is positively charged +2 (e/3 unit)
447  } else {
448  Qquarks -= value/absValue; // d or dbar or s or sbar : d- or s-quark is negatively charged -1 (e/3 unit)
449  }
450  }
451  if ( Qquarks < 0 ) NewProjCode *= -1;
452  #ifdef debugFTFexictation
453  G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
454  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
455  G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquarks<<G4endl;
456  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
457  #endif
458 
459  // Proj
460  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
461  if ( ! TestParticle ) continue;
462  common.MminProjectile = common.BrW.GetMinimumMass( TestParticle );
463  if ( common.SqrtS - common.M0target < common.MminProjectile ) continue;
464  MtestPr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
465  + 5.0*TestParticle->GetPDGWidth() );
466  #ifdef debugFTFexictation
467  G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()
468  << G4endl
469  << "MtestPart MtestPart0 "<<MtestPr<<" "<<TestParticle->GetPDGMass()<<G4endl
470  << "M0projectile projectile PDGMass " << common.M0projectile << " "
471  << projectile->GetDefinition()->GetPDGMass() << G4endl;
472  #endif
473 
474  // Targ
475  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
476  #ifdef debugFTFexictation
477  G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
478  << "NewTargCode " << NewTargCode << G4endl;
479  #endif
480  if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) { // Lambda or Sigma0 ?
481  if ( G4UniformRand() < 0.5 ) {
482  NewTargCode += 2;
483  } else if ( G4UniformRand() < 0.75 ) {
484  NewTargCode = 3122;
485  }
486  } else if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
487  NewTargCode += 2; ProjExcited = true; // Create Delta isobar
488  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
489  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
490  NewTargCode += 2; ProjExcited = true;
491  }
492  } else if ( ! ProjExcited &&
493  G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
494  common.SqrtS > common.M0projectile + // Delta mass
495  G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass() ) {
496  NewTargCode += 2; // Create Delta isobar
497  }
498  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
499  if ( ! TestParticle ) continue;
500  #ifdef debugFTFexictation
501  G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
502  #endif
503  common.MminTarget = common.BrW.GetMinimumMass( TestParticle );
504  if ( common.SqrtS - MtestPr < common.MminTarget ) continue;
505  MtestTr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
506  + 5.0*TestParticle->GetPDGWidth() );
507  if ( common.SqrtS > MtestPr + MtestTr ) break;
508 
509  } // End of while loop
510  if ( attempts >= maxNumberOfAttempts ) return returnCode; // unsuccessfully ended, nothing else can be done
511 
512  if ( MtestPr >= common.Pprojectile.mag() || projectile->GetStatus() != 0 ) {
513  common.M0projectile = MtestPr;
514  }
515  #ifdef debugFTFexictation
516  G4cout << "M0projectile After check " << common.M0projectile << G4endl;
517  #endif
518  common.M0projectile2 = common.M0projectile * common.M0projectile;
519  common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
520  common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
521  if ( MtestTr >= common.Ptarget.mag() || target->GetStatus() != 0 ) {
522  common.M0target = MtestTr;
523  }
524  common.M0target2 = common.M0target * common.M0target;
525  #ifdef debugFTFexictation
526  G4cout << "New targ M0 M0^2 " << common.M0target << " " << common.M0target2 << G4endl;
527  #endif
528  common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
529  common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
530 
531  } else { // of the if ( common.absProjectilePDGcode < 1000 ) ; the projectile is baryon now
532 
533  // Choose randomly, with equal probability, whether to consider the quarks of the
534  // projectile or target hadron for selecting the flavour of the exchanged quark.
535  G4bool isProjectileExchangedQ = false;
536  G4int firstQ = TargQ1, secondQ = TargQ2, thirdQ = TargQ3;
537  G4int otherFirstQ = ProjQ1, otherSecondQ = ProjQ2, otherThirdQ = ProjQ3;
538  if ( G4UniformRand() < 0.5 ) {
539  isProjectileExchangedQ = true;
540  firstQ = ProjQ1; secondQ = ProjQ2; thirdQ = ProjQ3;
541  otherFirstQ = TargQ1; otherSecondQ = TargQ2; otherThirdQ = TargQ3;
542  }
543  // Choose randomly, with equal probability, which of the three quarks of the
544  // selected (projectile or target) hadron is the exhanged quark.
545  G4int exchangedQ = 0;
546  G4double Ksi = G4UniformRand();
547  if ( Ksi < 0.333333 ) {
548  exchangedQ = firstQ;
549  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
550  exchangedQ = secondQ;
551  } else {
552  exchangedQ = thirdQ;
553  }
554  #ifdef debugFTFexictation
555  G4cout << "Exchange Qs isProjectile Q " << isProjectileExchangedQ << " " << exchangedQ << " ";
556  #endif
557  // The exchanged quarks (one of the projectile hadron and one of the target hadron)
558  // are always accepted if they have different flavour, else (i.e. when they have the
559  // same flavour) they are accepted only with a specified probability.
560  G4double probSame = theParameters->GetProbOfSameQuarkExchange();
561  const G4int MaxCount = 100;
562  G4int count = 0, otherExchangedQ = 0;
563  do {
564  if ( exchangedQ != otherFirstQ || G4UniformRand() < probSame ) {
565  otherExchangedQ = otherFirstQ; otherFirstQ = exchangedQ; exchangedQ = otherExchangedQ;
566  } else {
567  if ( exchangedQ != otherSecondQ || G4UniformRand() < probSame ) {
568  otherExchangedQ = otherSecondQ; otherSecondQ = exchangedQ; exchangedQ = otherExchangedQ;
569  } else {
570  //ALB if ( exchangedQ != otherThirdQ || G4UniformRand() < probSame ) {
571  otherExchangedQ = otherThirdQ; otherThirdQ = exchangedQ; exchangedQ = otherExchangedQ;
572  //ALB }
573  }
574  }
575  } while ( otherExchangedQ == 0 && ++count < MaxCount );
576  if ( count >= MaxCount ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
577  // Swap (between projectile and target hadron) the two quarks that have been sampled
578  // as "exchanged" quarks.
579  if ( Ksi < 0.333333 ) {
580  firstQ = exchangedQ;
581  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
582  secondQ = exchangedQ;
583  } else {
584  thirdQ = exchangedQ;
585  }
586  if ( isProjectileExchangedQ ) {
587  ProjQ1 = firstQ; ProjQ2 = secondQ; ProjQ3 = thirdQ;
588  TargQ1 = otherFirstQ; TargQ2 = otherSecondQ; TargQ3 = otherThirdQ;
589  } else {
590  TargQ1 = firstQ; TargQ2 = secondQ; TargQ3 = thirdQ;
591  ProjQ1 = otherFirstQ; ProjQ2 = otherSecondQ; ProjQ3 = otherThirdQ;
592  }
593  #ifdef debugFTFexictation
594  G4cout << "Exchange Qs Pr Tr " << ( isProjectileExchangedQ ? exchangedQ : otherExchangedQ )
595  << " " << ( isProjectileExchangedQ ? otherExchangedQ : exchangedQ ) << G4endl;
596  #endif
597 
598  NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
599  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
600  // Decide whether the new projectile hadron is a Delta particle;
601  // then decide whether the new target hadron is a Delta particle.
602  // Notice that a Delta particle has the last PDG digit "4" (because its spin is 3/2),
603  // whereas a nucleon has "2" (because its spin is 1/2).
604  for ( G4int iHadron = 0; iHadron < 2; iHadron++ ) {
605  // First projectile hadron, then target hadron
606  G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, codeQ3 = ProjQ3, newHadCode = NewProjCode;
607  G4double massConstraint = common.M0target;
608  G4bool isHadronADelta = ( projectile->GetDefinition()->GetPDGiIsospin() == 3 );
609  if ( iHadron == 1 ) { // Target hadron
610  codeQ1 = TargQ1, codeQ2 = TargQ2, codeQ3 = TargQ3, newHadCode = NewTargCode;
611  massConstraint = common.M0projectile;
612  isHadronADelta = ( target->GetDefinition()->GetPDGiIsospin() == 3 );
613  }
614  if ( codeQ1 == codeQ2 && codeQ1 == codeQ3 ) { // The three quarks are the same
615  newHadCode += 2; // Delta++ (uuu) or Delta- (ddd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
616  } else if ( isHadronADelta ) { // Hadron (projectile or target) was Delta
617  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
618  newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
619  } else {
620  newHadCode += 0; // No delta (so the last PDG digit remains 2)
621  }
622  } else { // Hadron (projectile or target) was Nucleon
623  if ( G4UniformRand() < DeltaProbAtQuarkExchange &&
624  common.SqrtS > G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass()
625  + massConstraint ) {
626  newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
627  } else {
628  newHadCode += 0; // No delta (so the last PDG digit remains 2)
629  }
630  }
631  if ( iHadron == 0 ) { // Projectile hadron
632  NewProjCode = newHadCode;
633  } else { // Target hadron
634  NewTargCode = newHadCode;
635  }
636  }
637  #ifdef debugFTFexictation
638  G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
639  #endif
640 
641  if ( common.absProjectilePDGcode == NewProjCode && common.absTargetPDGcode == NewTargCode ) {
642  } // Nothing was changed! It is not right!?
643 
644  // Sampling of the masses of the projectile and target nucleons.
645  // Because of energy conservation, the ordering of the sampling matters:
646  // randomly, half of the time we sample first the target nucleon mass and
647  // then the projectile nucleon mass, and the other half of the time we
648  // sample first the projectile nucleon mass and then the target nucleon mass.
649  G4VSplitableHadron* firstHadron = target;
650  G4VSplitableHadron* secondHadron = projectile;
651  G4int firstHadronCode = NewTargCode, secondHadronCode = NewProjCode;
652  G4double massConstraint = common.M0projectile;
653  G4bool isFirstTarget = true;
654  if ( G4UniformRand() < 0.5 ) {
655  // Sample first the projectile nucleon mass, then the target nucleon mass.
656  firstHadron = projectile; secondHadron = target;
657  firstHadronCode = NewProjCode; secondHadronCode = NewTargCode;
658  massConstraint = common.M0target;
659  isFirstTarget = false;
660  }
661  G4double Mtest1st = 0.0, Mtest2nd = 0.0, Mmin1st = 0.0, Mmin2nd = 0.0;
662  for ( int iSamplingCase = 0; iSamplingCase < 2; iSamplingCase++ ) {
663  G4VSplitableHadron* aHadron = firstHadron;
664  G4int aHadronCode = firstHadronCode;
665  if ( iSamplingCase == 1 ) { // Second nucleon mass sampling
666  aHadron = secondHadron; aHadronCode = secondHadronCode; massConstraint = Mtest1st;
667  }
668  G4double MtestHadron = 0.0, MminHadron = 0.0;
669  if ( aHadron->GetStatus() == 1 || aHadron->GetStatus() == 2 ) {
670  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( aHadronCode );
671  if ( ! TestParticle ) return returnCode; // Not possible to find such a hadron: unsuccessfully ended, nothing else can be done
672  MminHadron = common.BrW.GetMinimumMass( TestParticle );
673  if ( common.SqrtS - massConstraint < MminHadron ) return returnCode; // Kinematically impossible: unsuccessfully ended, nothing else can be done
674  if ( TestParticle->GetPDGWidth() == 0.0 ) {
675  MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass() );
676  } else {
677  const G4int maxNumberOfAttempts = 50;
678  G4int attempts = 0;
679  while ( attempts < maxNumberOfAttempts ) {
680  attempts++;
681  MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
682  + 5.0*TestParticle->GetPDGWidth() );
683  if ( common.SqrtS < MtestHadron + massConstraint ) {
684  continue; // Kinematically unacceptable: try again
685  } else {
686  break; // Kinematically acceptable: the mass sampling is successful
687  }
688  }
689  if ( attempts >= maxNumberOfAttempts ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
690  }
691  }
692  if ( iSamplingCase == 0 ) {
693  Mtest1st = MtestHadron; Mmin1st = MminHadron;
694  } else {
695  Mtest2nd = MtestHadron; Mmin2nd = MminHadron;
696  }
697  } // End for loop on the two sampling cases (1st and 2nd)
698  if ( isFirstTarget ) {
699  MtestTr = Mtest1st; MtestPr = Mtest2nd;
700  common.MminTarget = Mmin1st; common.MminProjectile = Mmin2nd;
701  } else {
702  MtestTr = Mtest2nd; MtestPr = Mtest1st;
703  common.MminTarget = Mmin2nd; common.MminProjectile = Mmin1st;
704  }
705 
706  if ( MtestPr != 0.0 ) {
707  common.M0projectile = MtestPr;
708  common.M0projectile2 = common.M0projectile * common.M0projectile;
709  common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
710  common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
711  }
712  if ( MtestTr != 0.0 ) {
713  common.M0target = MtestTr;
714  common.M0target2 = common.M0target * common.M0target;
715  common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
716  common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
717  }
718 
719  } // End of if ( common.absProjectilePDGcode < 1000 )
720 
721  // If we assume that final state hadrons after the charge exchange will be
722  // in the ground states, we have to put
723  if ( common.SqrtS < common.M0projectile + common.M0target ) return returnCode; // unsuccessfully ended, nothing else can be done
724 
725  common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
726  - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
727  + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
728  #ifdef debugFTFexictation
729  G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
730  << "At the end// NewTargCode " << NewTargCode << G4endl
731  << "M0pr M0tr SqS " << common.M0projectile << " " << common.M0target << " "
732  << common.SqrtS << G4endl
733  << "M0pr2 M0tr2 SqS " << common.M0projectile2 << " " << common.M0target2 << " "
734  << common.SqrtS << G4endl
735  << "PZcms2 after the change " << common.PZcms2 << G4endl << G4endl;
736  #endif
737  if ( common.PZcms2 < 0.0 ) return returnCode; // It can be if energy is not sufficient for Delta
738  // unsuccessfully ended, nothing else can be done
739  projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
740  target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
741  common.PZcms = std::sqrt( common.PZcms2 );
742  common.Pprojectile.setPz( common.PZcms );
743  common.Pprojectile.setE( std::sqrt( common.M0projectile2 + common.PZcms2 ) );
744  common.Ptarget.setPz( -common.PZcms );
745  common.Ptarget.setE( std::sqrt( common.M0target2 + common.PZcms2 ) );
746  #ifdef debugFTFexictation
747  G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << common.Pprojectile << G4endl
748  << common.Ptarget << G4endl << common.Pprojectile + common.Ptarget << G4endl;
749  #endif
750 
751  if ( projectile->GetStatus() != 0 ) projectile->SetStatus( 2 );
752  if ( target->GetStatus() != 0 ) target->SetStatus( 2 );
753 
754  // Check for possible excitation of the participants
755  if ( common.SqrtS < common.M0projectile + common.TargetDiffStateMinMass ||
756  common.SqrtS < common.ProjectileDiffStateMinMass + common.M0target ||
757  common.ProbOfDiffraction == 0.0 ) common.ProbExc = 0.0;
758 
759  if ( G4UniformRand() > common.ProbExc ) { // Make elastic scattering
760  #ifdef debugFTFexictation
761  G4cout << "Make elastic scattering of new hadrons" << G4endl;
762  #endif
763  common.Pprojectile.transform( common.toLab );
764  common.Ptarget.transform( common.toLab );
765  projectile->Set4Momentum( common.Pprojectile );
766  target->Set4Momentum( common.Ptarget );
767  G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
768  #ifdef debugFTFexictation
769  G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
770  << G4endl << projectile->Get4Momentum() << G4endl << target->Get4Momentum() << G4endl
771  << projectile->Get4Momentum() + target->Get4Momentum() << " "
772  << (projectile->Get4Momentum() + target->Get4Momentum()).mag() << G4endl;
773  #endif
774  if ( Result ) returnCode = 0; // successfully ended and nothing else needs to be done
775  return returnCode;
776  }
777 
778  #ifdef debugFTFexictation
779  G4cout << "Make excitation of new hadrons" << G4endl;
780  #endif
781 
782  // Redefinition of ProbOfDiffraction because the probabilities are changed due to quark exchange
783  common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
784  if ( common.ProbOfDiffraction != 0.0 ) {
785  common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
786  common.ProbTargetDiffraction /= common.ProbOfDiffraction;
787  }
788 
789  return returnCode = 1; // successfully completed, but the work needs to be continued
790 }
791 
792 //-----------------------------------------------------------------------------
793 
794 G4bool G4DiffractiveExcitation::
795 ExciteParticipants_doDiffraction( G4VSplitableHadron* projectile, G4VSplitableHadron* target,
796  G4FTFParameters* theParameters,
797  G4DiffractiveExcitation::CommonVariables& common ) const {
798  // Second of the three utility methods used only by ExciteParticipants:
799  // it does the sampling for the diffraction case, either projectile or target diffraction.
800 
801  G4bool isProjectileDiffraction = false;
802  if ( G4UniformRand() < common.ProbProjectileDiffraction ) { // projectile diffraction
803  isProjectileDiffraction = true;
804  #ifdef debugFTFexictation
805  G4cout << "projectile diffraction" << G4endl;
806  #endif
807  common.ProjMassT2 = common.ProjectileDiffStateMinMass2;
808  common.ProjMassT = common.ProjectileDiffStateMinMass;
809  common.TargMassT2 = common.M0target2;
810  common.TargMassT = common.M0target;
811  } else { // target diffraction
812  #ifdef debugFTFexictation
813  G4cout << "Target diffraction" << G4endl;
814  #endif
815  common.ProjMassT2 = common.M0projectile2;
816  common.ProjMassT = common.M0projectile;
817  common.TargMassT2 = common.TargetDiffStateMinMass2;
818  common.TargMassT = common.TargetDiffStateMinMass;
819  }
820 
821  G4double DiffrAveragePt2 = theParameters->GetAvaragePt2ofElasticScattering()*1.2;
822  G4bool loopCondition = true;
823  G4int whilecount = 0;
824  do { // Generate pt and mass of projectile
825 
826  whilecount++;
827  if ( whilecount > 1000 ) {
828  common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
829  return false; // Ignore this interaction
830  };
831 
832  // Check that the interaction is possible
833  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
834 
835  common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
836  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
837  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
838  if ( common.PZcms2 < 0.0 ) return false;
839 
840  common.maxPtSquare = common.PZcms2;
841  common.Qmomentum = G4LorentzVector( GaussianPt( DiffrAveragePt2, common.maxPtSquare ), 0 );
842  common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
843  if ( isProjectileDiffraction ) { // projectile diffraction
844  common.ProjMassT2 = common.ProjectileDiffStateMinMass2 + common.Pt2;
845  common.TargMassT2 = common.M0target2 + common.Pt2;
846  } else { // target diffraction
847  common.ProjMassT2 = common.M0projectile2 + common.Pt2;
848  common.TargMassT2 = common.TargetDiffStateMinMass2 + common.Pt2;
849  }
850  common.ProjMassT = std::sqrt( common.ProjMassT2 );
851  common.TargMassT = std::sqrt( common.TargMassT2 );
852  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
853 
854  common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
855  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
856  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
857  if ( common.PZcms2 < 0.0 ) continue;
858 
859  common.PZcms = std::sqrt( common.PZcms2 );
860  if ( isProjectileDiffraction ) { // projectile diffraction
861  common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
862  common.PMinusMax = common.SqrtS - common.TargMassT;
863  common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
864  common.TMinusNew = common.SqrtS - common.PMinusNew;
865  common.Qminus = common.Ptarget.minus() - common.TMinusNew;
866  common.TPlusNew = common.TargMassT2 / common.TMinusNew;
867  common.Qplus = common.Ptarget.plus() - common.TPlusNew;
868  common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
869  common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
870  loopCondition = ( ( common.Pprojectile + common.Qmomentum ).mag2() <
871  common.ProjectileDiffStateMinMass2 );
872  } else { // target diffraction
873  common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
874  common.TPlusMax = common.SqrtS - common.ProjMassT;
875  common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
876  common.PPlusNew = common.SqrtS - common.TPlusNew;
877  common.Qplus = common.PPlusNew - common.Pprojectile.plus();
878  common.PMinusNew = common.ProjMassT2 / common.PPlusNew;
879  common.Qminus = common.PMinusNew - common.Pprojectile.minus();
880  common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
881  common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
882  loopCondition = ( ( common.Ptarget - common.Qmomentum ).mag2() <
883  common.TargetDiffStateMinMass2 );
884  }
885 
886  } while ( loopCondition ); /* Loop checking, 10.08.2015, A.Ribon */
887  // Repeat the sampling because there was not any excitation
888 
889  if ( isProjectileDiffraction ) { // projectile diffraction
890  projectile->SetStatus( 0 );
891  if ( projectile->GetStatus() == 2 ) projectile->SetStatus( 1 );
892  if ( target->GetStatus() == 1 && target->GetSoftCollisionCount() == 0 ) target->SetStatus( 2 );
893  } else { // target diffraction
894  target->SetStatus( 0 );
895  }
896 
897  return true;
898 }
899 
900 //-----------------------------------------------------------------------------
901 
902 G4bool G4DiffractiveExcitation::
903 ExciteParticipants_doNonDiffraction( G4VSplitableHadron* projectile,
904  G4VSplitableHadron* target,
905  G4FTFParameters* theParameters,
906  G4DiffractiveExcitation::CommonVariables& common ) const {
907  // Third of the three utility methods used only by ExciteParticipants:
908  // it does the sampling for the non-diffraction case.
909 
910  #ifdef debugFTFexictation
911  G4cout << "Non-diffraction process" << G4endl;
912  #endif
913  G4int whilecount = 0;
914  do { // Generate pt and masses
915 
916  whilecount++;
917  if ( whilecount > 1000 ) {
918  common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
919  return false; // Ignore this interaction
920  };
921 
922  // Check that the interaction is possible
923  common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2;
924  common.ProjMassT = common.ProjectileNonDiffStateMinMass;
925  common.TargMassT2 = common.TargetNonDiffStateMinMass2;
926  common.TargMassT = common.TargetNonDiffStateMinMass;
927  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
928 
929  common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
930  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
931  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
932  if ( common.PZcms2 < 0.0 ) return false;
933 
934  common.maxPtSquare = common.PZcms2;
935  common.Qmomentum = G4LorentzVector( GaussianPt( theParameters->GetAveragePt2(),
936  common.maxPtSquare ), 0 );
937  common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
938  common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2 + common.Pt2;
939  common.ProjMassT = std::sqrt( common.ProjMassT2 );
940  common.TargMassT2 = common.TargetNonDiffStateMinMass2 + common.Pt2;
941  common.TargMassT = std::sqrt( common.TargMassT2 );
942  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
943 
944  common.PZcms2 =( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
945  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
946  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
947  if ( common.PZcms2 < 0.0 ) continue;
948 
949  common.PZcms = std::sqrt( common.PZcms2 );
950  common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
951  common.PMinusMax = common.SqrtS - common.TargMassT;
952  common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
953  common.TPlusMax = common.SqrtS - common.ProjMassT;
954  if ( G4UniformRand() < theParameters->GetProbLogDistrPrD() ) {
955  common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
956  } else {
957  common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*G4UniformRand() + common.PMinusMin;
958  }
959  if ( G4UniformRand() < theParameters->GetProbLogDistr() ) {
960  common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
961  } else {
962  common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*G4UniformRand() + common.TPlusMin;
963  }
964  common.Qminus = common.PMinusNew - common.Pprojectile.minus();
965  common.Qplus = -( common.TPlusNew - common.Ptarget.plus() );
966  common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
967  common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
968  #ifdef debugFTFexictation
969  G4cout <<"Sampled: Mpr, MdifPr, Mtr, MdifTr "<<G4endl
970  << ( common.Pprojectile + common.Qmomentum ).mag() << " "
971  << common.ProjectileNonDiffStateMinMass << G4endl
972  << ( common.Ptarget - common.Qmomentum ).mag() << " "
973  << common.TargetNonDiffStateMinMass << G4endl;
974  #endif
975 
976  } while ( ( common.Pprojectile + common.Qmomentum ).mag2() <
977  common.ProjectileNonDiffStateMinMass2 || // No double Diffraction
978  ( common.Ptarget - common.Qmomentum ).mag2() <
979  common.TargetNonDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
980 
981  projectile->SetStatus( 0 );
982  target->SetStatus( 0 );
983  return true;
984 }
985 
986 
987 //============================================================================
988 
989 void G4DiffractiveExcitation::CreateStrings( G4VSplitableHadron* hadron,
990  G4bool isProjectile,
991  G4ExcitedString*& FirstString,
992  G4ExcitedString*& SecondString,
993  G4FTFParameters* theParameters ) const {
994 
995  //G4cout << "Create Strings SplitUp " << hadron << G4endl
996  // << "Defin " << hadron->GetDefinition() << G4endl
997  // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
998 
999  hadron->SplitUp();
1000 
1001  G4Parton* start = hadron->GetNextParton();
1002  if ( start == NULL ) {
1003  G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1004  FirstString = 0; SecondString = 0;
1005  return;
1006  }
1007 
1008  G4Parton* end = hadron->GetNextParton();
1009  if ( end == NULL ) {
1010  G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1011  FirstString = 0; SecondString = 0;
1012  return;
1013  }
1014 
1015  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1016  // << G4endl
1017  // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1018 
1019  G4LorentzVector Phadron = hadron->Get4Momentum();
1020  //G4cout << "String mom " << Phadron << G4endl;
1021  G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1022  G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1023  G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1024  G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1025  G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1026 
1027  G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1028  G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1029  //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1030 
1031  G4double Wmin( 0.0 );
1032  if ( isProjectile ) {
1033  Wmin = theParameters->GetProjMinDiffMass();
1034  } else {
1035  Wmin = theParameters->GetTarMinDiffMass();
1036  }
1037 
1038  G4double W = hadron->Get4Momentum().mag();
1039  G4double W2 = W*W;
1040  G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); // x2( 0.0 )
1041  G4bool Kink = false;
1042 
1043  if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1044  end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1045  ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1046  end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1047  // Kinky strings are allowed only for qq-q strings;
1048  // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1049  // according to the analysis of Pbar P interactions
1050 
1051  if ( W > Wmin ) { // Kink is possible
1052  if ( hadron->GetStatus() == 0 ) {
1053  G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1054  if ( Pt2kink ) {
1055  Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );
1056  } else {
1057  Pt = 0.0;
1058  }
1059  } else {
1060  Pt = 0.0;
1061  }
1062 
1063  if ( Pt > 500.0*MeV ) {
1064  G4double Ymax = G4Log( W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1065  G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
1066  x1 = 1.0 - Pt/W * G4Exp( Y );
1067  x3 = 1.0 - Pt/W * G4Exp(-Y );
1068  //x2 = 2.0 - x1 - x3;
1069 
1070  G4double Mass_startQ = 650.0*MeV;
1071  if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV;
1072  if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV;
1073  if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV;
1074  G4double Mass_endQ = 650.0*MeV;
1075  if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV;
1076  if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV;
1077  if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV;
1078 
1079  G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1080  G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1081  G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
1082  if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1083  Kink = false;
1084  } else {
1085  G4double P_1 = std::sqrt( P2_1 );
1086  G4double P_2 = std::sqrt( P2_2 );
1087  G4double P_3 = std::sqrt( P2_3 );
1088  G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1089  G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1090 
1091  if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1092  Kink = false;
1093  } else {
1094  Kink = true;
1095  Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated
1096  Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1097  Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1098  Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1099  Pstart.setE( x3*W/2.0 );
1100  Pkink.setE( Pkink.vect().mag() );
1101  Pend.setE( x1*W/2.0 );
1102 
1103  G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1104  if ( Pkink.getZ() > 0.0 ) {
1105  if ( XkQ > 0.5 ) {
1106  PkinkQ1 = XkQ*Pkink;
1107  } else {
1108  PkinkQ1 = (1.0 - XkQ)*Pkink;
1109  }
1110  } else {
1111  if ( XkQ > 0.5 ) {
1112  PkinkQ1 = (1.0 - XkQ)*Pkink;
1113  } else {
1114  PkinkQ1 = XkQ*Pkink;
1115  }
1116  }
1117 
1118  PkinkQ2 = Pkink - PkinkQ1;
1119  // Minimizing Pt1^2+Pt3^2
1120  G4double Cos2Psi = ( sqr(x1) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1121  std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
1122  G4double Psi = std::acos( Cos2Psi );
1123 
1124  G4LorentzRotation Rotate;
1125  if ( isProjectile ) {
1126  Rotate.rotateY( Psi );
1127  } else {
1128  Rotate.rotateY( pi + Psi );
1129  }
1130  Rotate.rotateZ( twopi * G4UniformRand() );
1131  Pstart *= Rotate;
1132  Pkink *= Rotate;
1133  PkinkQ1 *= Rotate;
1134  PkinkQ2 *= Rotate;
1135  Pend *= Rotate;
1136  }
1137  } // End of if ( P2_1 <= 0.0 || P2_3 <= 0.0 )
1138  } // End of if ( Pt > 500.0*MeV )
1139  } // End of if ( W > Wmin ) : check for a kink
1140  } // end of qq-q string selection
1141 
1142  if ( Kink ) { // Kink is possible
1143 
1144  //G4cout << "Kink is sampled!" << G4endl;
1145  std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1146  theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1147 
1148  G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1149  for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1150  //G4cout << "Iq " << Iq << G4endl;
1151  if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1152  }
1153  //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1154  G4Parton* Gquark = new G4Parton( QuarkInGluon );
1155  G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1156  //G4cout << "Lorentz " << G4endl;
1157 
1158  G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1159  G4LorentzRotation toLab( toCMS.inverse() );
1160  //G4cout << "Pstart " << Pstart << G4endl;
1161  //G4cout << "Pend " << Pend << G4endl;
1162  //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1163  //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1164  //G4cout << "Pstart " << Pstart << G4endl<<G4endl;
1165 
1166  Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1167  PkinkQ1.transform( toLab );
1168  PkinkQ2.transform( toLab );
1169  Pend.transform( toLab ); end->Set4Momentum( Pend );
1170  //G4cout << "Pstart " << Pstart << G4endl;
1171  //G4cout << "Pend " << Pend << G4endl;
1172  //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1173  //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1174 
1175  //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1176  G4int absPDGcode = 1500; // 23 Dec
1177  if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1178  end->GetDefinition()->GetParticleSubType() == "quark" ) {
1179  absPDGcode = 110;
1180  }
1181  //G4cout << "absPDGcode " << absPDGcode << G4endl;
1182 
1183  if ( absPDGcode < 1000 ) { // meson
1184  if ( isProjectile ) { // Projectile
1185  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1186 
1187  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1188  SecondString = new G4ExcitedString( Gquark, start , +1 );
1189  Ganti_quark->Set4Momentum( PkinkQ1 );
1190  Gquark->Set4Momentum( PkinkQ2 );
1191  } else { // Anti_Quark on the end
1192 
1193  FirstString = new G4ExcitedString( end , Gquark, +1 );
1194  SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1195  Gquark->Set4Momentum( PkinkQ1 );
1196  Ganti_quark->Set4Momentum( PkinkQ2 );
1197  }
1198  } else { // Target
1199  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1200  FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1201  SecondString = new G4ExcitedString( start , Gquark, -1 );
1202  Ganti_quark->Set4Momentum( PkinkQ2 );
1203  Gquark->Set4Momentum( PkinkQ1 );
1204  } else { // Anti_Quark on the end
1205  FirstString = new G4ExcitedString( Gquark, end , -1 );
1206  SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1207  Gquark->Set4Momentum( PkinkQ2 );
1208  Ganti_quark->Set4Momentum( PkinkQ1 );
1209  }
1210  }
1211  } else { // Baryon/AntiBaryon
1212 
1213  if ( isProjectile ) { // Projectile
1214  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1215  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1216  FirstString = new G4ExcitedString( end , Gquark, +1 );
1217  SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1218  Gquark->Set4Momentum( PkinkQ1 );
1219  Ganti_quark->Set4Momentum( PkinkQ2 );
1220 
1221  } else { // Anti_DiQuark on the end or quark
1222  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1223  SecondString = new G4ExcitedString( Gquark, start , +1 );
1224  Ganti_quark->Set4Momentum( PkinkQ1 );
1225  Gquark->Set4Momentum( PkinkQ2 );
1226  }
1227  } else { // Target
1228  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1229  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1230 
1231  Gquark->Set4Momentum( PkinkQ1 );
1232  Ganti_quark->Set4Momentum( PkinkQ2 );
1233 
1234  FirstString = new G4ExcitedString( end, Gquark, -1 );
1235  SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1236 
1237  } else { // Anti_DiQuark on the end or Q
1238  FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1239  SecondString = new G4ExcitedString( start , Gquark, -1 );
1240  Gquark->Set4Momentum( PkinkQ2 );
1241  Ganti_quark->Set4Momentum( PkinkQ1 );
1242  }
1243  }
1244  }
1245 
1246  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1247  FirstString->SetPosition( hadron->GetPosition() );
1248  SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1249  SecondString->SetPosition( hadron->GetPosition() );
1250 
1251  } else { // End of kink is possible: Kink is impossible
1252 
1253  FirstString = new G4ExcitedString( end, start, +1 );
1254 
1255  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1256  FirstString->SetPosition( hadron->GetPosition() );
1257  SecondString = 0;
1258 
1259  if ( ! (end->Get4Momentum().e() != 0.) ) {
1260  // momenta of string ends
1261  G4double Momentum = hadron->Get4Momentum().vect().mag();
1262  G4double Plus = hadron->Get4Momentum().e() + Momentum;
1263  G4double Minus = hadron->Get4Momentum().e() - Momentum;
1264  G4ThreeVector tmp;
1265  if ( Momentum > 0.0 ) {
1266  tmp.set( hadron->Get4Momentum().px(),
1267  hadron->Get4Momentum().py(),
1268  hadron->Get4Momentum().pz() );
1269  tmp /= Momentum;
1270  } else {
1271  tmp.set( 0.0, 0.0, 1.0 );
1272  }
1273  G4LorentzVector Pstart1( tmp, 0.0 );
1274  G4LorentzVector Pend1( tmp, 0.0 );
1275  if ( isProjectile ) {
1276  Pstart1 *= (-1.0)*Minus/2.0;
1277  Pend1 *= (+1.0)*Plus /2.0;
1278  } else {
1279  Pstart1 *= (+1.0)*Plus/ 2.0;
1280  Pend1 *= (-1.0)*Minus/2.0;
1281  }
1282 
1283  Momentum = Pstart1.vect().mag();
1284 
1285  Pstart1.setT( Momentum ); // It is assumed that quark has m=0.
1286 
1287  Momentum = Pend1.vect().mag();
1288 
1289  Pend1.setT( Momentum ); // It is assumed that di-quark has m=0.
1290  start->Set4Momentum( Pstart1 );
1291  end->Set4Momentum( Pend1 );
1292  SecondString = 0;
1293  }
1294  } // End of kink is impossible
1295 
1296  //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1297  // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1298  // << FirstString << " " << SecondString << G4endl;
1299 
1300  #ifdef G4_FTFDEBUG
1301  G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1302  << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1303  << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1304  << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1305  << end->Get4Momentum().mag() << G4endl << " sum of ends "
1306  << Pstart + Pend << G4endl << " Original "
1307  << hadron->Get4Momentum() << " "<<hadron->Get4Momentum().mag() << G4endl;
1308  #endif
1309 
1310  return;
1311 }
1312 
1313 
1314 //============================================================================
1315 
1316 G4double G4DiffractiveExcitation::ChooseP( G4double Pmin, G4double Pmax ) const {
1317  // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1318  // To be improved...
1319  G4double range = Pmax - Pmin;
1320  if ( Pmin <= 0.0 || range <= 0.0 ) {
1321  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1322  throw G4HadronicException( __FILE__, __LINE__,
1323  "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1324  }
1325  G4double P = Pmin * G4Pow::GetInstance()->powA( Pmax/Pmin, G4UniformRand() );
1326  //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1327  return P;
1328 }
1329 
1330 
1331 //============================================================================
1332 
1333 G4ThreeVector G4DiffractiveExcitation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1334  // @@ this method is used in FTFModel as well. Should go somewhere common!
1335  G4double Pt2( 0.0 );
1336  if ( AveragePt2 <= 0.0 ) {
1337  Pt2 = 0.0;
1338  } else {
1339  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1340  ( G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1341  }
1342  G4double Pt = std::sqrt( Pt2 );
1343  G4double phi = G4UniformRand() * twopi;
1344  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1345 }
1346 
1347 
1348 //============================================================================
1349 
1350 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink( G4double zmin, G4double zmax ) const {
1351  G4double z, yf;
1352  const G4int maxNumberOfLoops = 10000;
1353  G4int loopCounter = 0;
1354  do {
1355  z = zmin + G4UniformRand() * (zmax - zmin);
1356  yf = z*z + sqr(1.0 - z);
1357  } while ( ( G4UniformRand() > yf ) &&
1358  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1359  if ( loopCounter >= maxNumberOfLoops ) {
1360  z = 0.5*(zmin + zmax); // Just something acceptable, without any physics consideration.
1361  }
1362  return z;
1363 }
1364 
1365 
1366 //============================================================================
1367 
1368 void G4DiffractiveExcitation::UnpackMeson( const G4int IdPDG, G4int& Q1, G4int& Q2 ) const {
1369  G4int absIdPDG = std::abs( IdPDG );
1370  if ( ! ( absIdPDG == 111 || absIdPDG == 221 || absIdPDG == 331 ) ) {
1371  // Ordinary mesons
1372  Q1 = absIdPDG / 100;
1373  Q2 = (absIdPDG % 100) / 10;
1374  G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1375  if ( IdPDG < 0 ) anti *= -1;
1376  Q1 *= anti;
1377  Q2 *= -1 * anti;
1378  } else {
1379  // Pi0, Eta, Eta'
1380  if ( G4UniformRand() < 0.5 ) { Q1 = 1; Q2 = -1; }
1381  else { Q1 = 2; Q2 = -2; }
1382  }
1383  return;
1384 }
1385 
1386 
1387 //============================================================================
1388 
1389 void G4DiffractiveExcitation::UnpackBaryon( G4int IdPDG,
1390  G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1391  Q1 = IdPDG / 1000;
1392  Q2 = (IdPDG % 1000) / 100;
1393  Q3 = (IdPDG % 100) / 10;
1394  return;
1395 }
1396 
1397 
1398 //============================================================================
1399 
1400 G4int G4DiffractiveExcitation::NewNucleonId( G4int Q1, G4int Q2, G4int Q3 ) const {
1401  G4int TmpQ( 0 );
1402  if ( Q3 > Q2 ) {
1403  TmpQ = Q2;
1404  Q2 = Q3;
1405  Q3 = TmpQ;
1406  } else if ( Q3 > Q1 ) {
1407  TmpQ = Q1;
1408  Q1 = Q3;
1409  Q3 = TmpQ;
1410  }
1411  if ( Q2 > Q1 ) {
1412  TmpQ = Q1;
1413  Q1 = Q2;
1414  Q2 = TmpQ;
1415  }
1416  G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1417  return NewCode;
1418 }
1419 
1420 
1421 //============================================================================
1422 
1423 G4DiffractiveExcitation::G4DiffractiveExcitation( const G4DiffractiveExcitation& ) {
1424  throw G4HadronicException( __FILE__, __LINE__,
1425  "G4DiffractiveExcitation copy constructor not meant to be called" );
1426 }
1427 
1428 
1429 //============================================================================
1430 
1431 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=( const G4DiffractiveExcitation& ) {
1432  throw G4HadronicException( __FILE__, __LINE__,
1433  "G4DiffractiveExcitation = operator not meant to be called" );
1434  return *this;
1435 }
1436 
1437 
1438 //============================================================================
1439 
1440 G4bool G4DiffractiveExcitation::operator==( const G4DiffractiveExcitation& ) const {
1441  throw G4HadronicException( __FILE__, __LINE__,
1442  "G4DiffractiveExcitation == operator not meant to be called" );
1443 }
1444 
1445 
1446 //============================================================================
1447 
1448 G4bool G4DiffractiveExcitation::operator!= ( const G4DiffractiveExcitation& ) const {
1449  throw G4HadronicException( __FILE__, __LINE__,
1450  "G4DiffractiveExcitation != operator not meant to be called" );
1451 }
1452