Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4DiffractiveExcitation.cc 110870 2018-06-22 12:14:16Z gcosmo $
28 //
29 
30 // ------------------------------------------------------------
31 // GEANT 4 class implemetation file
32 //
33 // ---------------- G4DiffractiveExcitation --------------
34 // by Gunter Folger, October 1998.
35 // diffractive Excitation used by strings models
36 // Take a projectile and a target
37 // excite the projectile and target
38 // Essential changed by V. Uzhinsky in November - December 2006
39 // in order to put it in a correspondence with original FRITIOF
40 // model. Variant of FRITIOF with nucleon de-excitation is implemented.
41 // Other changes by V.Uzhinsky in May 2007 were introduced to fit
42 // meson-nucleon interactions. Additional changes by V. Uzhinsky
43 // were introduced in December 2006. They treat diffraction dissociation
44 // processes more exactly.
45 // Correct treatment of the diffraction dissociation - 2012, V. Uzhinsky
46 // Mass distributions for resonances and uu-diquark suppression in protons,
47 // and dd-diquarks suppression in neutrons were introduced by V. Uzhinsky, 2014
48 // ---------------------------------------------------------------------
49 
50 #include "globals.hh"
51 #include "Randomize.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 
56 #include "G4FTFParameters.hh"
57 #include "G4ElasticHNScattering.hh"
58 
59 #include "G4LorentzRotation.hh"
60 #include "G4RotationMatrix.hh"
61 #include "G4ThreeVector.hh"
62 #include "G4ParticleDefinition.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4SampleResonance.hh"
65 #include "G4VSplitableHadron.hh"
66 #include "G4ExcitedString.hh"
67 #include "G4Neutron.hh"
68 
69 #include "G4Exp.hh"
70 #include "G4Log.hh"
71 #include "G4Pow.hh"
72 
73 //#include "G4ios.hh"
74 //#include "UZHI_diffraction.hh"
75 
76 
77 //============================================================================
78 
79 //#define debugFTFexictation
80 
81 
82 //============================================================================
83 
85 
86 
87 //============================================================================
88 
90 
91 
92 //============================================================================
93 
96  G4FTFParameters* theParameters,
97  G4ElasticHNScattering* theElastic ) const {
98 
99  #ifdef debugFTFexictation
100  G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
101  #endif
102 
104 
105  // Projectile parameters
106  common.Pprojectile = projectile->Get4Momentum();
107  if ( common.Pprojectile.z() < 0.0 ) return false;
108  common.ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
109  common.absProjectilePDGcode = std::abs( common.ProjectilePDGcode );
110  common.M0projectile = common.Pprojectile.mag();
111  G4double ProjectileRapidity = common.Pprojectile.rapidity();
112 
113  // Target parameters
114  common.Ptarget = target->Get4Momentum();
115  common.TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
116  common.absTargetPDGcode = std::abs( common.TargetPDGcode );
117  common.M0target = common.Ptarget.mag();
118  G4double TargetRapidity = common.Ptarget.rapidity();
119 
120  // Kinematical properties of the interactions
121  G4LorentzVector Psum = common.Pprojectile + common.Ptarget; // 4-momentum in Lab.
122  common.S = Psum.mag2();
123  common.SqrtS = std::sqrt( common.S );
124 
125  // Check off-shellness of the participants
126  G4bool toBePutOnMassShell = false;
127  common.MminProjectile = common.BrW.GetMinimumMass( projectile->GetDefinition() );
128  if ( common.M0projectile < common.MminProjectile ) {
129  toBePutOnMassShell = true;
130  common.M0projectile = common.BrW.SampleMass( projectile->GetDefinition(),
131  projectile->GetDefinition()->GetPDGMass()
132  + 5.0*projectile->GetDefinition()->GetPDGWidth() );
133  }
134  common.M0projectile2 = common.M0projectile * common.M0projectile;
135  common.ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass();
136  common.ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass();
137  if ( common.M0projectile > common.ProjectileDiffStateMinMass ) {
138  common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV;
139  common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV;
140  if ( common.absProjectilePDGcode > 3000 ) { // Strange baryon
141  common.ProjectileDiffStateMinMass += 140.0*MeV;
142  common.ProjectileNonDiffStateMinMass += 140.0*MeV;
143  }
144  }
145  common.MminTarget = common.BrW.GetMinimumMass( target->GetDefinition() );
146  if ( common.M0target < common.MminTarget ) {
147  toBePutOnMassShell = true;
148  common.M0target = common.BrW.SampleMass( target->GetDefinition(),
149  target->GetDefinition()->GetPDGMass()
150  + 5.0*target->GetDefinition()->GetPDGWidth() );
151  }
152  common.M0target2 = common.M0target * common.M0target;
153  common.TargetDiffStateMinMass = theParameters->GetTarMinDiffMass();
154  common.TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass();
155  if ( common.M0target > common.TargetDiffStateMinMass ) {
156  common.TargetDiffStateMinMass = common.M0target + 220.0*MeV;
157  common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV;
158  if ( common.absTargetPDGcode > 3000 ) { // Strange baryon
159  common.TargetDiffStateMinMass += 140.0*MeV;
160  common.TargetNonDiffStateMinMass += 140.0*MeV;
161  }
162  };
163  #ifdef debugFTFexictation
164  G4cout << "Proj Targ PDGcodes " << common.ProjectilePDGcode << " " << common.TargetPDGcode << G4endl
165  << "M0projectile Y " << common.M0projectile << " " << ProjectileRapidity << G4endl;
166  //G4cout << "M0target Y " << common.M0target << " " << TargetRapidity << G4endl;
167  G4cout << "Pproj " << common.Pprojectile << G4endl << "Ptarget " << common.Ptarget << G4endl;
168  #endif
169 
170  // Transform momenta to cms and then rotate parallel to z axis;
171  common.toCms = G4LorentzRotation( -1 * Psum.boostVector() );
172  G4LorentzVector Ptmp = common.toCms * common.Pprojectile;
173  if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
174  common.toCms.rotateZ( -1*Ptmp.phi() );
175  common.toCms.rotateY( -1*Ptmp.theta() );
176  common.toLab = common.toCms.inverse();
177  common.Pprojectile.transform( common.toCms );
178  common.Ptarget.transform( common.toCms );
179 
180  G4double SumMasses = common.M0projectile + common.M0target; // + 220.0*MeV;
181  #ifdef debugFTFexictation
182  G4cout << "SqrtS " << common.SqrtS << G4endl << "M0pr M0tr SumM " << common.M0projectile
183  << " " << common.M0target << " " << SumMasses << G4endl;
184  #endif
185  if ( common.SqrtS < SumMasses ) return false; // The model cannot work at low energy
186 
187  common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
188  - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
189  + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
190  #ifdef debugFTFexictation
191  G4cout << "PZcms2 after toBePutOnMassShell " << common.PZcms2 << G4endl;
192  #endif
193  if ( common.PZcms2 < 0.0 ) return false; // It can be in an interaction with off-shell nuclear nucleon
194 
195  common.PZcms = std::sqrt( common.PZcms2 );
196  if ( toBePutOnMassShell ) {
197  if ( common.Pprojectile.z() > 0.0 ) {
198  common.Pprojectile.setPz( common.PZcms );
199  common.Ptarget.setPz( -common.PZcms );
200  } else {
201  common.Pprojectile.setPz( -common.PZcms );
202  common.Ptarget.setPz( common.PZcms );
203  }
204  common.Pprojectile.setE( std::sqrt( common.M0projectile2
205  + common.Pprojectile.x() * common.Pprojectile.x()
206  + common.Pprojectile.y() * common.Pprojectile.y()
207  + common.PZcms2 ) );
208  common.Ptarget.setE( std::sqrt( common.M0target2
209  + common.Ptarget.x() * common.Ptarget.x()
210  + common.Ptarget.y() * common.Ptarget.y()
211  + common.PZcms2 ) );
212  }
213  #ifdef debugFTFexictation
214  G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << common.M0projectile
215  << " " << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
216  << G4endl
217  << "Targ M0 Mdif Mndif " << common.M0target << " " << common.TargetDiffStateMinMass
218  << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS << G4endl
219  << "Proj CMS " << common.Pprojectile << G4endl << "Targ CMS " << common.Ptarget << G4endl;
220  #endif
221 
222  // Check for possible quark exchange
223  ProjectileRapidity = common.Pprojectile.rapidity();
224  TargetRapidity = common.Ptarget.rapidity();
225  G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity - TargetRapidity );
226  G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity - TargetRapidity ) *
227  theParameters->GetProcProb( 4, ProjectileRapidity - TargetRapidity );
228  common.ProbProjectileDiffraction =
229  theParameters->GetProcProb( 2, ProjectileRapidity - TargetRapidity );
230  common.ProbTargetDiffraction =
231  theParameters->GetProcProb( 3, ProjectileRapidity - TargetRapidity );
233  #ifdef debugFTFexictation
234  G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
235  << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
236  << "ProjectileRapidity " << ProjectileRapidity << G4endl;
237  #endif
238  if ( QeNoExc + QeExc + common.ProbProjectileDiffraction + common.ProbTargetDiffraction > 1.0 ) {
239  QeNoExc = 1.0 - QeExc - common.ProbProjectileDiffraction - common.ProbTargetDiffraction;
240  }
241  if ( QeExc + QeNoExc != 0.0 ) {
242  common.ProbExc = QeExc / ( QeExc + QeNoExc );
243  }
244  if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
245  common.ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
246  common.ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
247  }
248  #ifdef debugFTFexictation
249  G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
250  << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
251  << "ProjectileRapidity " << ProjectileRapidity << G4endl;
252  #endif
253 
254  // Try out charge exchange
255  G4int returnCode = 1;
256  if ( G4UniformRand() < QeExc + QeNoExc ) {
257  returnCode = ExciteParticipants_doChargeExchange( projectile, target, theParameters,
258  theElastic, common );
259  }
260 
261  G4bool returnResult = false;
262  if ( returnCode == 0 ) {
263  returnResult = true; // Successfully ended: no need of extra work
264  } else if ( returnCode == 1 ) {
265 
267  #ifdef debugFTFexictation
268  G4cout << "Excitation --------------------" << G4endl
269  << "Proj M0 MdMin MndMin " << common.M0projectile << " "
271  << G4endl
272  << "Targ M0 MdMin MndMin " << common.M0target << " " << common.TargetDiffStateMinMass
273  << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS
274  << G4endl
275  << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
276  << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
277  #endif
278  if ( common.ProbOfDiffraction != 0.0 ) {
280  } else {
281  common.ProbProjectileDiffraction = 0.0;
282  }
283  #ifdef debugFTFexictation
284  G4cout << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
285  << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
286  #endif
291  // Choose between diffraction and non-diffraction process
292  if ( G4UniformRand() < common.ProbOfDiffraction ) {
293 
294  returnResult = ExciteParticipants_doDiffraction( projectile, target, theParameters, common );
295 
296  } else {
297 
298  returnResult = ExciteParticipants_doNonDiffraction( projectile, target, theParameters, common );
299 
300  }
301  if ( returnResult ) {
302  common.Pprojectile += common.Qmomentum;
303  common.Ptarget -= common.Qmomentum;
304  // Transform back and update SplitableHadron Participant.
305  common.Pprojectile.transform( common.toLab );
306  common.Ptarget.transform( common.toLab );
307  #ifdef debugFTFexictation
308  G4cout << "Mproj " << common.Pprojectile.mag() << G4endl << "Mtarg " << common.Ptarget.mag()
309  << G4endl;
310  #endif
311  projectile->Set4Momentum( common.Pprojectile );
312  target->Set4Momentum( common.Ptarget );
313  projectile->IncrementCollisionCount( 1 );
314  target->IncrementCollisionCount( 1 );
315  }
316  }
317 
318  return returnResult;
319 }
320 
321 //-----------------------------------------------------------------------------
322 
326  G4FTFParameters* theParameters,
327  G4ElasticHNScattering* theElastic,
329  // First of the three utility methods used only by ExciteParticipants:
330  // it does the sampling for the charge-exchange case.
331  // This method returns an integer code - instead of a boolean, with the following meaning:
332  // "0" : successfully ended and nothing else needs to be done;
333  // "1" : successfully completed, but the work needs to be continued;
334  // "99" : unsuccessfully ended, nothing else can be done.
335  G4int returnCode = 99;
336 
337  G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
338  G4ParticleDefinition* TestParticle = 0;
339  G4double MtestPr = 0.0, MtestTr = 0.0;
340 
341  #ifdef debugFTFexictation
342  G4cout << "Q exchange --------------------------" << G4endl;
343  #endif
344  G4int NewProjCode = 0, NewTargCode = 0, ProjQ1 = 0, ProjQ2 = 0, ProjQ3 = 0;
345  // Projectile unpacking
346  if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
347  UnpackMeson( common.ProjectilePDGcode, ProjQ1, ProjQ2 );
348  } else { // projectile is baryon
349  UnpackBaryon( common.ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
350  }
351  // Target unpacking
352  G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;
353  UnpackBaryon( common.TargetPDGcode, TargQ1, TargQ2, TargQ3 );
354  #ifdef debugFTFexictation
355  G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
356  << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
357  #endif
358 
359  // Sampling of exchanged quarks
360  G4int ProjExchangeQ = 0, TargExchangeQ = 0;
361  if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
362 
363  G4bool isProjQ1Quark = false;
364  ProjExchangeQ = ProjQ2;
365  if ( ProjQ1 > 0 ) { // ProjQ1 is a quark
366  isProjQ1Quark = true;
367  ProjExchangeQ = ProjQ1;
368  }
369  // Exchange of non-identical quarks is allowed
370  G4int NpossibleStates = 0;
371  if ( ProjExchangeQ != TargQ1 ) NpossibleStates++;
372  if ( ProjExchangeQ != TargQ2 ) NpossibleStates++;
373  if ( ProjExchangeQ != TargQ3 ) NpossibleStates++;
374  G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
375  NpossibleStates = 0;
376  if ( ProjExchangeQ != TargQ1 ) {
377  if ( ++NpossibleStates == Nsampled ) {
378  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ;
379  isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
380  }
381  }
382  if ( ProjExchangeQ != TargQ2 ) {
383  if ( ++NpossibleStates == Nsampled ) {
384  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ;
385  isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
386  }
387  }
388  if ( ProjExchangeQ != TargQ3 ) {
389  if ( ++NpossibleStates == Nsampled ) {
390  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ;
391  isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
392  }
393  }
394  #ifdef debugFTFexictation
395  G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
396  #endif
397 
398  G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ2 = std::abs( ProjQ2 );
399  G4bool ProjExcited = false;
400  const G4int maxNumberOfAttempts = 50;
401  G4int attempts = 0;
402  while ( attempts++ < maxNumberOfAttempts ) { /* Loop checking, 10.08.2015, A.Ribon */
403 
404  // Determination of a new projectile ID which satisfies energy-momentum conservation
405  G4double Ksi = G4UniformRand();
406  if ( aProjQ1 == aProjQ2 ) {
407  if ( aProjQ1 != 3 ) {
408  NewProjCode = 111; // Pi0-meson
409  if ( Ksi < 0.5 ) {
410  NewProjCode = 221; // Eta -meson
411  if ( Ksi < 0.25 ) {
412  NewProjCode = 331; // Eta'-meson
413  }
414  }
415  } else {
416  NewProjCode = 221; // Eta -meson
417  if ( Ksi < 0.5 ) {
418  NewProjCode = 331; // Eta'-meson
419  }
420  }
421  } else {
422  if ( aProjQ1 > aProjQ2 ) {
423  NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
424  } else {
425  NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
426  }
427  }
428  #ifdef debugFTFexictation
429  G4cout << "NewProjCode " << NewProjCode << G4endl;
430  #endif
431  ProjExcited = false;
432  if ( G4UniformRand() < 0.5 ) {
433  NewProjCode += 2; // Excited meson (J=1 -> 2*J+1=2+1=3, last digit of the PDG code)
434  ProjExcited = true;
435  }
436  // So far we have used the absolute values of the PDG codes of the two constituent quarks:
437  // now look at their signed values to set properly the signed of the meson's PDG code.
438  G4int value = ProjQ1, absValue = aProjQ1, Qquarks = 0;
439  for ( G4int iQuark = 0; iQuark < 2; iQuark++ ) {
440  if ( iQuark == 1 ) {
441  value = ProjQ2; absValue = aProjQ2;
442  }
443  if ( absValue == 2 ) {
444  Qquarks += value; // u or ubar : u-quark is positively charged +2 (e/3 unit)
445  } else {
446  Qquarks -= value/absValue; // d or dbar or s or sbar : d- or s-quark is negatively charged -1 (e/3 unit)
447  }
448  }
449  if ( Qquarks < 0 ) NewProjCode *= -1;
450  #ifdef debugFTFexictation
451  G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
452  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
453  G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquarks<<G4endl;
454  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
455  #endif
456 
457  // Proj
458  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
459  if ( ! TestParticle ) continue;
460  common.MminProjectile = common.BrW.GetMinimumMass( TestParticle );
461  if ( common.SqrtS - common.M0target < common.MminProjectile ) continue;
462  MtestPr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
463  + 5.0*TestParticle->GetPDGWidth() );
464  #ifdef debugFTFexictation
465  G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()
466  << G4endl
467  << "MtestPart MtestPart0 "<<MtestPr<<" "<<TestParticle->GetPDGMass()<<G4endl
468  << "M0projectile projectile PDGMass " << common.M0projectile << " "
469  << projectile->GetDefinition()->GetPDGMass() << G4endl;
470  #endif
471 
472  // Targ
473  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
474  #ifdef debugFTFexictation
475  G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
476  << "NewTargCode " << NewTargCode << G4endl;
477  #endif
478  if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) { // Lambda or Sigma0 ?
479  if ( G4UniformRand() < 0.5 ) {
480  NewTargCode += 2;
481  } else if ( G4UniformRand() < 0.75 ) {
482  NewTargCode = 3122;
483  }
484  } else if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
485  NewTargCode += 2; ProjExcited = true; // Create Delta isobar
486  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
487  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
488  NewTargCode += 2; ProjExcited = true;
489  }
490  } else if ( ! ProjExcited &&
491  G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
492  common.SqrtS > common.M0projectile + // Delta mass
494  NewTargCode += 2; // Create Delta isobar
495  }
496  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
497  if ( ! TestParticle ) continue;
498  #ifdef debugFTFexictation
499  G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
500  #endif
501  common.MminTarget = common.BrW.GetMinimumMass( TestParticle );
502  if ( common.SqrtS - MtestPr < common.MminTarget ) continue;
503  MtestTr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
504  + 5.0*TestParticle->GetPDGWidth() );
505  if ( common.SqrtS > MtestPr + MtestTr ) break;
506 
507  } // End of while loop
508  if ( attempts >= maxNumberOfAttempts ) return returnCode; // unsuccessfully ended, nothing else can be done
509 
510  if ( MtestPr >= common.Pprojectile.mag() || projectile->GetStatus() != 0 ) {
511  common.M0projectile = MtestPr;
512  }
513  #ifdef debugFTFexictation
514  G4cout << "M0projectile After check " << common.M0projectile << G4endl;
515  #endif
516  common.M0projectile2 = common.M0projectile * common.M0projectile;
517  common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
518  common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
519  if ( MtestTr >= common.Ptarget.mag() || target->GetStatus() != 0 ) {
520  common.M0target = MtestTr;
521  }
522  common.M0target2 = common.M0target * common.M0target;
523  #ifdef debugFTFexictation
524  G4cout << "New targ M0 M0^2 " << common.M0target << " " << common.M0target2 << G4endl;
525  #endif
526  common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
527  common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
528 
529  } else { // of the if ( common.absProjectilePDGcode < 1000 ) ; the projectile is baryon now
530 
531  // Choose randomly, with equal probability, whether to consider the quarks of the
532  // projectile or target hadron for selecting the flavour of the exchanged quark.
533  G4bool isProjectileExchangedQ = false;
534  G4int firstQ = TargQ1, secondQ = TargQ2, thirdQ = TargQ3;
535  G4int otherFirstQ = ProjQ1, otherSecondQ = ProjQ2, otherThirdQ = ProjQ3;
536  if ( G4UniformRand() < 0.5 ) {
537  isProjectileExchangedQ = true;
538  firstQ = ProjQ1; secondQ = ProjQ2; thirdQ = ProjQ3;
539  otherFirstQ = TargQ1; otherSecondQ = TargQ2; otherThirdQ = TargQ3;
540  }
541  // Choose randomly, with equal probability, which of the three quarks of the
542  // selected (projectile or target) hadron is the exhanged quark.
543  G4int exchangedQ = 0;
544  G4double Ksi = G4UniformRand();
545  if ( Ksi < 0.333333 ) {
546  exchangedQ = firstQ;
547  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
548  exchangedQ = secondQ;
549  } else {
550  exchangedQ = thirdQ;
551  }
552  #ifdef debugFTFexictation
553  G4cout << "Exchange Qs isProjectile Q " << isProjectileExchangedQ << " " << exchangedQ << " ";
554  #endif
555  // The exchanged quarks (one of the projectile hadron and one of the target hadron)
556  // are always accepted if they have different flavour, else (i.e. when they have the
557  // same flavour) they are accepted only with a specified probability.
558  G4double probSame = theParameters->GetProbOfSameQuarkExchange();
559  const G4int MaxCount = 100;
560  G4int count = 0, otherExchangedQ = 0;
561  do {
562  if ( exchangedQ != otherFirstQ || G4UniformRand() < probSame ) {
563  otherExchangedQ = otherFirstQ; otherFirstQ = exchangedQ; exchangedQ = otherExchangedQ;
564  } else {
565  if ( exchangedQ != otherSecondQ || G4UniformRand() < probSame ) {
566  otherExchangedQ = otherSecondQ; otherSecondQ = exchangedQ; exchangedQ = otherExchangedQ;
567  } else {
568  //ALB if ( exchangedQ != otherThirdQ || G4UniformRand() < probSame ) {
569  otherExchangedQ = otherThirdQ; otherThirdQ = exchangedQ; exchangedQ = otherExchangedQ;
570  //ALB }
571  }
572  }
573  } while ( otherExchangedQ == 0 && ++count < MaxCount );
574  if ( count >= MaxCount ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
575  // Swap (between projectile and target hadron) the two quarks that have been sampled
576  // as "exchanged" quarks.
577  if ( Ksi < 0.333333 ) {
578  firstQ = exchangedQ;
579  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
580  secondQ = exchangedQ;
581  } else {
582  thirdQ = exchangedQ;
583  }
584  if ( isProjectileExchangedQ ) {
585  ProjQ1 = firstQ; ProjQ2 = secondQ; ProjQ3 = thirdQ;
586  TargQ1 = otherFirstQ; TargQ2 = otherSecondQ; TargQ3 = otherThirdQ;
587  } else {
588  TargQ1 = firstQ; TargQ2 = secondQ; TargQ3 = thirdQ;
589  ProjQ1 = otherFirstQ; ProjQ2 = otherSecondQ; ProjQ3 = otherThirdQ;
590  }
591  #ifdef debugFTFexictation
592  G4cout << "Exchange Qs Pr Tr " << ( isProjectileExchangedQ ? exchangedQ : otherExchangedQ )
593  << " " << ( isProjectileExchangedQ ? otherExchangedQ : exchangedQ ) << G4endl;
594  #endif
595 
596  NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
597  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
598  // Decide whether the new projectile hadron is a Delta particle;
599  // then decide whether the new target hadron is a Delta particle.
600  // Notice that a Delta particle has the last PDG digit "4" (because its spin is 3/2),
601  // whereas a nucleon has "2" (because its spin is 1/2).
602  for ( G4int iHadron = 0; iHadron < 2; iHadron++ ) {
603  // First projectile hadron, then target hadron
604  G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, codeQ3 = ProjQ3, newHadCode = NewProjCode;
605  G4double massConstraint = common.M0target;
606  G4bool isHadronADelta = ( projectile->GetDefinition()->GetPDGiIsospin() == 3 );
607  if ( iHadron == 1 ) { // Target hadron
608  codeQ1 = TargQ1, codeQ2 = TargQ2, codeQ3 = TargQ3, newHadCode = NewTargCode;
609  massConstraint = common.M0projectile;
610  isHadronADelta = ( target->GetDefinition()->GetPDGiIsospin() == 3 );
611  }
612  if ( codeQ1 == codeQ2 && codeQ1 == codeQ3 ) { // The three quarks are the same
613  newHadCode += 2; // Delta++ (uuu) or Delta- (ddd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
614  } else if ( isHadronADelta ) { // Hadron (projectile or target) was Delta
615  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
616  newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
617  } else {
618  newHadCode += 0; // No delta (so the last PDG digit remains 2)
619  }
620  } else { // Hadron (projectile or target) was Nucleon
621  if ( G4UniformRand() < DeltaProbAtQuarkExchange &&
623  + massConstraint ) {
624  newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
625  } else {
626  newHadCode += 0; // No delta (so the last PDG digit remains 2)
627  }
628  }
629  if ( iHadron == 0 ) { // Projectile hadron
630  NewProjCode = newHadCode;
631  } else { // Target hadron
632  NewTargCode = newHadCode;
633  }
634  }
635  #ifdef debugFTFexictation
636  G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
637  #endif
638 
639  if ( common.absProjectilePDGcode == NewProjCode && common.absTargetPDGcode == NewTargCode ) {
640  } // Nothing was changed! It is not right!?
641 
642  // Sampling of the masses of the projectile and target nucleons.
643  // Because of energy conservation, the ordering of the sampling matters:
644  // randomly, half of the time we sample first the target nucleon mass and
645  // then the projectile nucleon mass, and the other half of the time we
646  // sample first the projectile nucleon mass and then the target nucleon mass.
647  G4VSplitableHadron* firstHadron = target;
648  G4VSplitableHadron* secondHadron = projectile;
649  G4int firstHadronCode = NewTargCode, secondHadronCode = NewProjCode;
650  G4double massConstraint = common.M0projectile;
651  G4bool isFirstTarget = true;
652  if ( G4UniformRand() < 0.5 ) {
653  // Sample first the projectile nucleon mass, then the target nucleon mass.
654  firstHadron = projectile; secondHadron = target;
655  firstHadronCode = NewProjCode; secondHadronCode = NewTargCode;
656  massConstraint = common.M0target;
657  isFirstTarget = false;
658  }
659  G4double Mtest1st = 0.0, Mtest2nd = 0.0, Mmin1st = 0.0, Mmin2nd = 0.0;
660  for ( int iSamplingCase = 0; iSamplingCase < 2; iSamplingCase++ ) {
661  G4VSplitableHadron* aHadron = firstHadron;
662  G4int aHadronCode = firstHadronCode;
663  if ( iSamplingCase == 1 ) { // Second nucleon mass sampling
664  aHadron = secondHadron; aHadronCode = secondHadronCode; massConstraint = Mtest1st;
665  }
666  G4double MtestHadron = 0.0, MminHadron = 0.0;
667  if ( aHadron->GetStatus() == 1 || aHadron->GetStatus() == 2 ) {
668  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( aHadronCode );
669  if ( ! TestParticle ) return returnCode; // Not possible to find such a hadron: unsuccessfully ended, nothing else can be done
670  MminHadron = common.BrW.GetMinimumMass( TestParticle );
671  if ( common.SqrtS - massConstraint < MminHadron ) return returnCode; // Kinematically impossible: unsuccessfully ended, nothing else can be done
672  if ( TestParticle->GetPDGWidth() == 0.0 ) {
673  MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass() );
674  } else {
675  const G4int maxNumberOfAttempts = 50;
676  G4int attempts = 0;
677  while ( attempts < maxNumberOfAttempts ) {
678  attempts++;
679  MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
680  + 5.0*TestParticle->GetPDGWidth() );
681  if ( common.SqrtS < MtestHadron + massConstraint ) {
682  continue; // Kinematically unacceptable: try again
683  } else {
684  break; // Kinematically acceptable: the mass sampling is successful
685  }
686  }
687  if ( attempts >= maxNumberOfAttempts ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
688  }
689  }
690  if ( iSamplingCase == 0 ) {
691  Mtest1st = MtestHadron; Mmin1st = MminHadron;
692  } else {
693  Mtest2nd = MtestHadron; Mmin2nd = MminHadron;
694  }
695  } // End for loop on the two sampling cases (1st and 2nd)
696  if ( isFirstTarget ) {
697  MtestTr = Mtest1st; MtestPr = Mtest2nd;
698  common.MminTarget = Mmin1st; common.MminProjectile = Mmin2nd;
699  } else {
700  MtestTr = Mtest2nd; MtestPr = Mtest1st;
701  common.MminTarget = Mmin2nd; common.MminProjectile = Mmin1st;
702  }
703 
704  if ( MtestPr != 0.0 ) {
705  common.M0projectile = MtestPr;
706  common.M0projectile2 = common.M0projectile * common.M0projectile;
707  common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
708  common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
709  }
710  if ( MtestTr != 0.0 ) {
711  common.M0target = MtestTr;
712  common.M0target2 = common.M0target * common.M0target;
713  common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
714  common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
715  }
716 
717  } // End of if ( common.absProjectilePDGcode < 1000 )
718 
719  // If we assume that final state hadrons after the charge exchange will be
720  // in the ground states, we have to put
721  if ( common.SqrtS < common.M0projectile + common.M0target ) return returnCode; // unsuccessfully ended, nothing else can be done
722 
723  common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
724  - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
725  + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
726  #ifdef debugFTFexictation
727  G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
728  << "At the end// NewTargCode " << NewTargCode << G4endl
729  << "M0pr M0tr SqS " << common.M0projectile << " " << common.M0target << " "
730  << common.SqrtS << G4endl
731  << "M0pr2 M0tr2 SqS " << common.M0projectile2 << " " << common.M0target2 << " "
732  << common.SqrtS << G4endl
733  << "PZcms2 after the change " << common.PZcms2 << G4endl << G4endl;
734  #endif
735  if ( common.PZcms2 < 0.0 ) return returnCode; // It can be if energy is not sufficient for Delta
736  // unsuccessfully ended, nothing else can be done
737  projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
738  target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
739  common.PZcms = std::sqrt( common.PZcms2 );
740  common.Pprojectile.setPz( common.PZcms );
741  common.Pprojectile.setE( std::sqrt( common.M0projectile2 + common.PZcms2 ) );
742  common.Ptarget.setPz( -common.PZcms );
743  common.Ptarget.setE( std::sqrt( common.M0target2 + common.PZcms2 ) );
744  #ifdef debugFTFexictation
745  G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << common.Pprojectile << G4endl
746  << common.Ptarget << G4endl << common.Pprojectile + common.Ptarget << G4endl;
747  #endif
748 
749  if ( projectile->GetStatus() != 0 ) projectile->SetStatus( 2 );
750  if ( target->GetStatus() != 0 ) target->SetStatus( 2 );
751 
752  // Check for possible excitation of the participants
753  if ( common.SqrtS < common.M0projectile + common.TargetDiffStateMinMass ||
754  common.SqrtS < common.ProjectileDiffStateMinMass + common.M0target ||
755  common.ProbOfDiffraction == 0.0 ) common.ProbExc = 0.0;
756 
757  if ( G4UniformRand() > common.ProbExc ) { // Make elastic scattering
758  #ifdef debugFTFexictation
759  G4cout << "Make elastic scattering of new hadrons" << G4endl;
760  #endif
761  common.Pprojectile.transform( common.toLab );
762  common.Ptarget.transform( common.toLab );
763  projectile->Set4Momentum( common.Pprojectile );
764  target->Set4Momentum( common.Ptarget );
765  G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
766  #ifdef debugFTFexictation
767  G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
768  << G4endl << projectile->Get4Momentum() << G4endl << target->Get4Momentum() << G4endl
769  << projectile->Get4Momentum() + target->Get4Momentum() << " "
770  << (projectile->Get4Momentum() + target->Get4Momentum()).mag() << G4endl;
771  #endif
772  if ( Result ) returnCode = 0; // successfully ended and nothing else needs to be done
773  return returnCode;
774  }
775 
776  #ifdef debugFTFexictation
777  G4cout << "Make excitation of new hadrons" << G4endl;
778  #endif
779 
780  // Redefinition of ProbOfDiffraction because the probabilities are changed due to quark exchange
782  if ( common.ProbOfDiffraction != 0.0 ) {
784  common.ProbTargetDiffraction /= common.ProbOfDiffraction;
785  }
786 
787  return returnCode = 1; // successfully completed, but the work needs to be continued
788 }
789 
790 //-----------------------------------------------------------------------------
791 
794  G4FTFParameters* theParameters,
796  // Second of the three utility methods used only by ExciteParticipants:
797  // it does the sampling for the diffraction case, either projectile or target diffraction.
798 
799  G4bool isProjectileDiffraction = false;
800  if ( G4UniformRand() < common.ProbProjectileDiffraction ) { // projectile diffraction
801  isProjectileDiffraction = true;
802  #ifdef debugFTFexictation
803  G4cout << "projectile diffraction" << G4endl;
804  #endif
805  common.ProjMassT2 = common.ProjectileDiffStateMinMass2;
806  common.ProjMassT = common.ProjectileDiffStateMinMass;
807  common.TargMassT2 = common.M0target2;
808  common.TargMassT = common.M0target;
809  } else { // target diffraction
810  #ifdef debugFTFexictation
811  G4cout << "Target diffraction" << G4endl;
812  #endif
813  common.ProjMassT2 = common.M0projectile2;
814  common.ProjMassT = common.M0projectile;
815  common.TargMassT2 = common.TargetDiffStateMinMass2;
816  common.TargMassT = common.TargetDiffStateMinMass;
817  }
818 
819  G4double DiffrAveragePt2 = theParameters->GetAvaragePt2ofElasticScattering()*1.2;
820  G4bool loopCondition = true;
821  G4int whilecount = 0;
822  do { // Generate pt and mass of projectile
823 
824  whilecount++;
825  if ( whilecount > 1000 ) {
826  common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
827  return false; // Ignore this interaction
828  };
829 
830  // Check that the interaction is possible
831  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
832 
833  common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
834  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
835  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
836  if ( common.PZcms2 < 0.0 ) return false;
837 
838  common.maxPtSquare = common.PZcms2;
839  common.Qmomentum = G4LorentzVector( GaussianPt( DiffrAveragePt2, common.maxPtSquare ), 0 );
840  common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
841  if ( isProjectileDiffraction ) { // projectile diffraction
842  common.ProjMassT2 = common.ProjectileDiffStateMinMass2 + common.Pt2;
843  common.TargMassT2 = common.M0target2 + common.Pt2;
844  } else { // target diffraction
845  common.ProjMassT2 = common.M0projectile2 + common.Pt2;
846  common.TargMassT2 = common.TargetDiffStateMinMass2 + common.Pt2;
847  }
848  common.ProjMassT = std::sqrt( common.ProjMassT2 );
849  common.TargMassT = std::sqrt( common.TargMassT2 );
850  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
851 
852  common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
853  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
854  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
855  if ( common.PZcms2 < 0.0 ) continue;
856 
857  common.PZcms = std::sqrt( common.PZcms2 );
858  if ( isProjectileDiffraction ) { // projectile diffraction
859  common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
860  common.PMinusMax = common.SqrtS - common.TargMassT;
861  common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
862  common.TMinusNew = common.SqrtS - common.PMinusNew;
863  common.Qminus = common.Ptarget.minus() - common.TMinusNew;
864  common.TPlusNew = common.TargMassT2 / common.TMinusNew;
865  common.Qplus = common.Ptarget.plus() - common.TPlusNew;
866  common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
867  common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
868  loopCondition = ( ( common.Pprojectile + common.Qmomentum ).mag2() <
869  common.ProjectileDiffStateMinMass2 );
870  } else { // target diffraction
871  common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
872  common.TPlusMax = common.SqrtS - common.ProjMassT;
873  common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
874  common.PPlusNew = common.SqrtS - common.TPlusNew;
875  common.Qplus = common.PPlusNew - common.Pprojectile.plus();
876  common.PMinusNew = common.ProjMassT2 / common.PPlusNew;
877  common.Qminus = common.PMinusNew - common.Pprojectile.minus();
878  common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
879  common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
880  loopCondition = ( ( common.Ptarget - common.Qmomentum ).mag2() <
881  common.TargetDiffStateMinMass2 );
882  }
883 
884  } while ( loopCondition ); /* Loop checking, 10.08.2015, A.Ribon */
885  // Repeat the sampling because there was not any excitation
886 
887  if ( isProjectileDiffraction ) { // projectile diffraction
888  projectile->SetStatus( 0 );
889  if ( projectile->GetStatus() == 2 ) projectile->SetStatus( 1 );
890  if ( target->GetStatus() == 1 && target->GetSoftCollisionCount() == 0 ) target->SetStatus( 2 );
891  } else { // target diffraction
892  target->SetStatus( 0 );
893  }
894 
895  return true;
896 }
897 
898 //-----------------------------------------------------------------------------
899 
903  G4FTFParameters* theParameters,
905  // Third of the three utility methods used only by ExciteParticipants:
906  // it does the sampling for the non-diffraction case.
907 
908  #ifdef debugFTFexictation
909  G4cout << "Non-diffraction process" << G4endl;
910  #endif
911  G4int whilecount = 0;
912  do { // Generate pt and masses
913 
914  whilecount++;
915  if ( whilecount > 1000 ) {
916  common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
917  return false; // Ignore this interaction
918  };
919 
920  // Check that the interaction is possible
923  common.TargMassT2 = common.TargetNonDiffStateMinMass2;
924  common.TargMassT = common.TargetNonDiffStateMinMass;
925  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
926 
927  common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
928  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
929  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
930  if ( common.PZcms2 < 0.0 ) return false;
931 
932  common.maxPtSquare = common.PZcms2;
933  common.Qmomentum = G4LorentzVector( GaussianPt( theParameters->GetAveragePt2(),
934  common.maxPtSquare ), 0 );
935  common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
936  common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2 + common.Pt2;
937  common.ProjMassT = std::sqrt( common.ProjMassT2 );
938  common.TargMassT2 = common.TargetNonDiffStateMinMass2 + common.Pt2;
939  common.TargMassT = std::sqrt( common.TargMassT2 );
940  if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
941 
942  common.PZcms2 =( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
943  - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
944  + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
945  if ( common.PZcms2 < 0.0 ) continue;
946 
947  common.PZcms = std::sqrt( common.PZcms2 );
948  common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
949  common.PMinusMax = common.SqrtS - common.TargMassT;
950  common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
951  common.TPlusMax = common.SqrtS - common.ProjMassT;
952  if ( G4UniformRand() < theParameters->GetProbLogDistrPrD() ) {
953  common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
954  } else {
955  common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*G4UniformRand() + common.PMinusMin;
956  }
957  if ( G4UniformRand() < theParameters->GetProbLogDistr() ) {
958  common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
959  } else {
960  common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*G4UniformRand() + common.TPlusMin;
961  }
962  common.Qminus = common.PMinusNew - common.Pprojectile.minus();
963  common.Qplus = -( common.TPlusNew - common.Ptarget.plus() );
964  common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
965  common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
966  #ifdef debugFTFexictation
967  G4cout <<"Sampled: Mpr, MdifPr, Mtr, MdifTr "<<G4endl
968  << ( common.Pprojectile + common.Qmomentum ).mag() << " "
969  << common.ProjectileNonDiffStateMinMass << G4endl
970  << ( common.Ptarget - common.Qmomentum ).mag() << " "
971  << common.TargetNonDiffStateMinMass << G4endl;
972  #endif
973 
974  } while ( ( common.Pprojectile + common.Qmomentum ).mag2() <
975  common.ProjectileNonDiffStateMinMass2 || // No double Diffraction
976  ( common.Ptarget - common.Qmomentum ).mag2() <
977  common.TargetNonDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
978 
979  projectile->SetStatus( 0 );
980  target->SetStatus( 0 );
981  return true;
982 }
983 
984 
985 //============================================================================
986 
988  G4bool isProjectile,
989  G4ExcitedString*& FirstString,
990  G4ExcitedString*& SecondString,
991  G4FTFParameters* theParameters ) const {
992 
993  //G4cout << "Create Strings SplitUp " << hadron << G4endl
994  // << "Defin " << hadron->GetDefinition() << G4endl
995  // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
996 
997  hadron->SplitUp();
998 
999  G4Parton* start = hadron->GetNextParton();
1000  if ( start == NULL ) {
1001  G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1002  FirstString = 0; SecondString = 0;
1003  return;
1004  }
1005 
1006  G4Parton* end = hadron->GetNextParton();
1007  if ( end == NULL ) {
1008  G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1009  FirstString = 0; SecondString = 0;
1010  return;
1011  }
1012 
1013  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1014  // << G4endl
1015  // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1016 
1017  G4LorentzVector Phadron = hadron->Get4Momentum();
1018  //G4cout << "String mom " << Phadron << G4endl;
1019  G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1020  G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1021  G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1022  G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1023  G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1024 
1025  G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1026  G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1027  //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1028 
1029  G4double Wmin( 0.0 );
1030  if ( isProjectile ) {
1031  Wmin = theParameters->GetProjMinDiffMass();
1032  } else {
1033  Wmin = theParameters->GetTarMinDiffMass();
1034  }
1035 
1036  G4double W = hadron->Get4Momentum().mag();
1037  //G4cout << "Wmin W " << Wmin << " " << W << G4endl;
1038  //G4int Uzhi; G4cin >> Uzhi;
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  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetParticleSubType() << " "
1143  // << end->GetDefinition()->GetParticleSubType() << G4endl;
1144  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetPDGEncoding() << " "
1145  // << end->GetDefinition()->GetPDGEncoding() << G4endl;
1146  //G4int Uzhi; G4cin >> Uzhi;
1147 
1148  if ( Kink ) { // Kink is possible
1149 
1150  //G4cout << "Kink is sampled!" << G4endl;
1151  std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1152  theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1153 
1154  G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1155  for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1156  //G4cout << "Iq " << Iq << G4endl;
1157  if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1158  }
1159  //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1160  G4Parton* Gquark = new G4Parton( QuarkInGluon );
1161  G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1162  //G4cout << "Lorentz " << G4endl;
1163 
1164  G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1165  G4LorentzRotation toLab( toCMS.inverse() );
1166  //G4cout << "Pstart " << Pstart << G4endl;
1167  //G4cout << "Pend " << Pend << G4endl;
1168  //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1169  //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1170  //G4cout << "Pstart " << Pstart << G4endl<<G4endl;
1171 
1172  Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1173  PkinkQ1.transform( toLab );
1174  PkinkQ2.transform( toLab );
1175  Pend.transform( toLab ); end->Set4Momentum( Pend );
1176  //G4cout << "Pstart " << Pstart << G4endl;
1177  //G4cout << "Pend " << Pend << G4endl;
1178  //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1179  //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1180 
1181  //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1182  G4int absPDGcode = 1500; // 23 Dec
1183  if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1184  end->GetDefinition()->GetParticleSubType() == "quark" ) {
1185  absPDGcode = 110;
1186  }
1187  //G4cout << "absPDGcode " << absPDGcode << G4endl;
1188 
1189  if ( absPDGcode < 1000 ) { // meson
1190  if ( isProjectile ) { // Projectile
1191  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1192 
1193  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1194  SecondString = new G4ExcitedString( Gquark, start , +1 );
1195  Ganti_quark->Set4Momentum( PkinkQ1 );
1196  Gquark->Set4Momentum( PkinkQ2 );
1197  } else { // Anti_Quark on the end
1198 
1199  FirstString = new G4ExcitedString( end , Gquark, +1 );
1200  SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1201  Gquark->Set4Momentum( PkinkQ1 );
1202  Ganti_quark->Set4Momentum( PkinkQ2 );
1203  }
1204  } else { // Target
1205  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1206  FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1207  SecondString = new G4ExcitedString( start , Gquark, -1 );
1208  Ganti_quark->Set4Momentum( PkinkQ2 );
1209  Gquark->Set4Momentum( PkinkQ1 );
1210  } else { // Anti_Quark on the end
1211  FirstString = new G4ExcitedString( Gquark, end , -1 );
1212  SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1213  Gquark->Set4Momentum( PkinkQ2 );
1214  Ganti_quark->Set4Momentum( PkinkQ1 );
1215  }
1216  }
1217  } else { // Baryon/AntiBaryon
1218 
1219  if ( isProjectile ) { // Projectile
1220  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1221  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1222  FirstString = new G4ExcitedString( end , Gquark, +1 ); // Open Uzhi
1223  SecondString = new G4ExcitedString( Ganti_quark, start , +1 ); // Open Uzhi
1224  Gquark->Set4Momentum( PkinkQ1 );
1225  Ganti_quark->Set4Momentum( PkinkQ2 );
1226 
1227  } else { // Anti_DiQuark on the end or quark
1228  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1229  SecondString = new G4ExcitedString( Gquark, start , +1 );
1230  Ganti_quark->Set4Momentum( PkinkQ1 );
1231  Gquark->Set4Momentum( PkinkQ2 );
1232  }
1233  } else { // Target
1234  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1235  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1236 
1237  Gquark->Set4Momentum( PkinkQ1 );
1238  Ganti_quark->Set4Momentum( PkinkQ2 );
1239 
1240  FirstString = new G4ExcitedString( end, Gquark, -1 );
1241  SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1242 
1243  } else { // Anti_DiQuark on the end or Q
1244  FirstString = new G4ExcitedString( Ganti_quark, end , -1 ); // Uzhi ?????
1245  SecondString = new G4ExcitedString( start , Gquark, -1 );
1246  Gquark->Set4Momentum( PkinkQ2 );
1247  Ganti_quark->Set4Momentum( PkinkQ1 );
1248  }
1249  }
1250  }
1251 
1252  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1253  FirstString->SetPosition( hadron->GetPosition() );
1254  SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1255  SecondString->SetPosition( hadron->GetPosition() );
1256 
1257  } else { // End of kink is possible: Kink is impossible
1258 
1259  FirstString = new G4ExcitedString( end, start, +1 );
1260 
1261  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1262  FirstString->SetPosition( hadron->GetPosition() );
1263  SecondString = 0;
1264 
1265  if ( ! (end->Get4Momentum().e() != 0.) ) { //AR-Oct2017
1266  // momenta of string ends
1267  G4double Momentum = hadron->Get4Momentum().vect().mag();
1268  G4double Plus = hadron->Get4Momentum().e() + Momentum;
1269  G4double Minus = hadron->Get4Momentum().e() - Momentum;
1271  if ( Momentum > 0.0 ) {
1272  tmp.set( hadron->Get4Momentum().px(),
1273  hadron->Get4Momentum().py(),
1274  hadron->Get4Momentum().pz() );
1275  tmp /= Momentum;
1276  } else {
1277  tmp.set( 0.0, 0.0, 1.0 );
1278  }
1279  G4LorentzVector Pstart1( tmp, 0.0 );
1280  G4LorentzVector Pend1( tmp, 0.0 );
1281  if ( isProjectile ) {
1282  Pstart1 *= (-1.0)*Minus/2.0;
1283  Pend1 *= (+1.0)*Plus /2.0;
1284  } else {
1285  Pstart1 *= (+1.0)*Plus/ 2.0;
1286  Pend1 *= (-1.0)*Minus/2.0;
1287  }
1288 
1289  //AR-Oct2017 Momentum = -Pstart1.mag();
1290  Momentum = Pstart1.vect().mag();
1291 
1292  Pstart1.setT( Momentum ); // It is assumed that quark has m=0.
1293 
1294  //AR-Oct2017 Momentum = -Pend1.mag();
1295  Momentum = Pend1.vect().mag();
1296 
1297  Pend1.setT( Momentum ); // It is assumed that di-quark has m=0.
1298  start->Set4Momentum( Pstart1 );
1299  end->Set4Momentum( Pend1 );
1300  SecondString = 0;
1301  } //AR-Oct2017
1302  } // End of kink is impossible
1303 
1304  //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1305  // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1306  // << FirstString << " " << SecondString << G4endl;
1307 
1308  #ifdef G4_FTFDEBUG
1309  G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1310  << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1311  << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1312  << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1313  << end->Get4Momentum().mag() << G4endl << " sum of ends "
1314  << Pstart + Pend << G4endl << " Original "
1315  << hadron->Get4Momentum() << " "<<hadron->Get4Momentum().mag() << G4endl;
1316  #endif
1317 
1318  return;
1319 }
1320 
1321 
1322 //============================================================================
1323 
1325  // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1326  // To be improved...
1327  G4double range = Pmax - Pmin;
1328  if ( Pmin <= 0.0 || range <= 0.0 ) {
1329  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1330  throw G4HadronicException( __FILE__, __LINE__,
1331  "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1332  }
1333  G4double P = Pmin * G4Pow::GetInstance()->powA( Pmax/Pmin, G4UniformRand() );
1334  //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1335  return P;
1336 }
1337 
1338 
1339 //============================================================================
1340 
1342  // @@ this method is used in FTFModel as well. Should go somewhere common!
1343  G4double Pt2( 0.0 );
1344  if ( AveragePt2 <= 0.0 ) {
1345  Pt2 = 0.0;
1346  } else {
1347  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1348  ( G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1349  }
1350  G4double Pt = std::sqrt( Pt2 );
1351  G4double phi = G4UniformRand() * twopi;
1352  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1353 }
1354 
1355 
1356 //============================================================================
1357 
1359  G4double z, yf;
1360  const G4int maxNumberOfLoops = 10000;
1361  G4int loopCounter = 0;
1362  do {
1363  z = zmin + G4UniformRand() * (zmax - zmin);
1364  yf = z*z + sqr(1.0 - z);
1365  } while ( ( G4UniformRand() > yf ) &&
1366  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1367  if ( loopCounter >= maxNumberOfLoops ) {
1368  z = 0.5*(zmin + zmax); // Just something acceptable, without any physics consideration.
1369  }
1370  return z;
1371 }
1372 
1373 
1374 //============================================================================
1375 
1376 void G4DiffractiveExcitation::UnpackMeson( const G4int IdPDG, G4int& Q1, G4int& Q2 ) const {
1377  G4int absIdPDG = std::abs( IdPDG );
1378  if ( ! ( absIdPDG == 111 || absIdPDG == 221 || absIdPDG == 331 ) ) {
1379  // Ordinary mesons
1380  Q1 = absIdPDG / 100;
1381  Q2 = (absIdPDG % 100) / 10;
1382  G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1383  if ( IdPDG < 0 ) anti *= -1;
1384  Q1 *= anti;
1385  Q2 *= -1 * anti;
1386  } else {
1387  // Pi0, Eta, Eta'
1388  if ( G4UniformRand() < 0.5 ) { Q1 = 1; Q2 = -1; }
1389  else { Q1 = 2; Q2 = -2; }
1390  }
1391  return;
1392 }
1393 
1394 
1395 //============================================================================
1396 
1398  G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1399  Q1 = IdPDG / 1000;
1400  Q2 = (IdPDG % 1000) / 100;
1401  Q3 = (IdPDG % 100) / 10;
1402  return;
1403 }
1404 
1405 
1406 //============================================================================
1407 
1409  G4int TmpQ( 0 );
1410  if ( Q3 > Q2 ) {
1411  TmpQ = Q2;
1412  Q2 = Q3;
1413  Q3 = TmpQ;
1414  } else if ( Q3 > Q1 ) {
1415  TmpQ = Q1;
1416  Q1 = Q3;
1417  Q3 = TmpQ;
1418  }
1419  if ( Q2 > Q1 ) {
1420  TmpQ = Q1;
1421  Q1 = Q2;
1422  Q2 = TmpQ;
1423  }
1424  G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1425  return NewCode;
1426 }
1427 
1428 
1429 //============================================================================
1430 
1432  throw G4HadronicException( __FILE__, __LINE__,
1433  "G4DiffractiveExcitation copy contructor not meant to be called" );
1434 }
1435 
1436 
1437 //============================================================================
1438 
1440  throw G4HadronicException( __FILE__, __LINE__,
1441  "G4DiffractiveExcitation = operator not meant to be called" );
1442  return *this;
1443 }
1444 
1445 
1446 //============================================================================
1447 
1449  throw G4HadronicException( __FILE__, __LINE__,
1450  "G4DiffractiveExcitation == operator not meant to be called" );
1451 }
1452 
1453 
1454 //============================================================================
1455 
1457  throw G4HadronicException( __FILE__, __LINE__,
1458  "G4DiffractiveExcitation != operator not meant to be called" );
1459 }
1460 
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
void SetStatus(const G4int aStatus)
G4double GetProcProb(const G4int ProcN, const G4double y)
void set(double x, double y, double z)
double mag2() const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetPDGWidth() const
const G4ParticleDefinition * GetDefinition() const
G4double GetProbOfSameQuarkExchange()
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
const G4DiffractiveExcitation & operator=(const G4DiffractiveExcitation &right)
HepLorentzVector & transform(const HepRotation &)
G4double GetAvaragePt2ofElasticScattering()
CLHEP::Hep3Vector G4ThreeVector
G4double GetProjMinDiffMass()
int operator!=(const G4DiffractiveExcitation &right) const
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
G4double ChooseP(G4double Pmin, G4double Pmax) const
double plus() const
const XML_Char * target
Definition: expat.h:268
void UnpackBaryon(G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * GetDefinition()
Definition: G4Parton.hh:161
static constexpr double MeV
Definition: G4SIunits.hh:214
double py() const
Float_t x1[n_points_granero]
Definition: compare.C:5
const G4String & GetParticleSubType() const
#define G4endl
Definition: G4ios.hh:61
double px() const
G4double GetTarMinNonDiffMass()
Double_t z
Float_t Y
const G4String & GetParticleName() const
G4double GetProbLogDistrPrD()
const G4String & GetParticleType() const
G4double GetProbLogDistr()
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
static int FASTCALL common(PROLOG_STATE *state, int tok)
Definition: xmlrole.cc:1309
Float_t tmp
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
virtual G4Parton * GetNextParton()=0
G4double GetDeltaProbAtQuarkExchange()
G4int NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const
G4double GetPDGMass() const
double rapidity() const
int operator==(const G4DiffractiveExcitation &right) const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
void SetTimeOfCreation(G4double aTime)
int Ymax
G4bool ExciteParticipants_doDiffraction(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, CommonVariables &common) const
CLHEP::HepLorentzRotation G4LorentzRotation
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4double GetAveragePt2()
if(nlines<=0)
G4double GetProjMinNonDiffMass()
G4double GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
HepLorentzRotation & rotateZ(double delta)
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ThreeVector & GetPosition() const
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
long G4long
Definition: G4Types.hh:80
const G4ParticleDefinition const G4Material *G4double range
double minus() const
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
HepLorentzRotation inverse() const
virtual void SplitUp()=0
static double P[]
double mag() const
void SetPosition(const G4ThreeVector &aPosition)
G4int GetPDGcode() const
Definition: G4Parton.hh:127
double mag() const
int G4int
Definition: G4Types.hh:78
void SetDefinition(const G4ParticleDefinition *aDefinition)
double pz() const
G4bool ExciteParticipants_doNonDiffraction(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, CommonVariables &common) const
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
Hep3Vector vect() const
G4int ExciteParticipants_doChargeExchange(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic, CommonVariables &common) const
T sqr(const T &x)
Definition: templates.hh:145
static constexpr double pi
Definition: G4SIunits.hh:75
CLHEP::HepLorentzVector G4LorentzVector
HepLorentzRotation & rotateY(double delta)
G4double GetPt2Kink()
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
void UnpackMeson(G4int IdPDG, G4int &Q1, G4int &Q2) const
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
double getZ() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
void IncrementCollisionCount(G4int aCount)
G4double GetTarMinDiffMass()
G4double GetMinimumMass(const G4ParticleDefinition *p) const