Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DiffractiveExcitation.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4DiffractiveExcitation.cc 100828 2016-11-02 15:25:59Z 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 
55 #include "G4DiffractiveExcitation.hh"
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 
84 G4DiffractiveExcitation::G4DiffractiveExcitation() {}
85 
86 
87 //============================================================================
88 
89 G4DiffractiveExcitation::~G4DiffractiveExcitation() {}
90 
91 
92 //============================================================================
93 
94 G4bool G4DiffractiveExcitation::ExciteParticipants( G4VSplitableHadron* projectile,
95  G4VSplitableHadron* target,
96  G4FTFParameters* theParameters,
97  G4ElasticHNScattering* theElastic ) const {
98 
99  #ifdef debugFTFexictation
100  G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
101  #endif
102 
103  // Projectile parameters
104  G4LorentzVector Pprojectile = projectile->Get4Momentum();
105  if ( Pprojectile.z() < 0.0 ) return false;
106 
107  G4int ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
108  G4int absProjectilePDGcode = std::abs( ProjectilePDGcode );
109  G4double M0projectile = Pprojectile.mag();
110  G4double ProjectileRapidity = Pprojectile.rapidity();
111 
112  // Target parameters
113  G4LorentzVector Ptarget = target->Get4Momentum();
114  G4int TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
115  G4int absTargetPDGcode = std::abs( TargetPDGcode );
116  G4double M0target = Ptarget.mag();
117  //G4double TargetRapidity = Ptarget.rapidity();
118 
119  // Kinematical properties of the interactions
120  G4LorentzVector Psum = Pprojectile + Ptarget; // 4-momentum in CMS
121  G4double S = Psum.mag2();
122  G4double SqrtS = std::sqrt( S );
123  //Uzhi_SqrtS = std::sqrt( S );
124 
125  // Check off-shellness of the participants
126  G4SampleResonance BrW;
127 
128  G4bool PutOnMassShell( false );
129 
130  G4double MminProjectile( 0.0 );
131  MminProjectile = BrW.GetMinimumMass(projectile->GetDefinition());
132  //G4double M0projectile = MminProjectile; // With de-excitation
133  //G4double M0projectile = Pprojectile.mag(); // Without de-excitation, see above
134 
135  if ( M0projectile < MminProjectile ) {
136  PutOnMassShell = true;
137  M0projectile = BrW.SampleMass( projectile->GetDefinition(),
138  projectile->GetDefinition()->GetPDGMass() + 5.0*projectile->GetDefinition()->GetPDGWidth() );
139  }
140 
141  G4double M0projectile2 = M0projectile * M0projectile;
142  G4double ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass();
143  G4double ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass();
144  if ( M0projectile > ProjectileDiffStateMinMass ) {
145  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV;
146  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV;
147  if ( absProjectilePDGcode > 3000 ) { // Strange baryon
148  ProjectileDiffStateMinMass += 140.0*MeV;
149  ProjectileNonDiffStateMinMass += 140.0*MeV;
150  }
151  }
152 
153  G4double MminTarget( 0.0 );
154  MminTarget = BrW.GetMinimumMass( target->GetDefinition() );
155  if ( M0target < MminTarget ) {
156  PutOnMassShell = true;
157  M0target = BrW.SampleMass( target->GetDefinition(),
158  target->GetDefinition()->GetPDGMass() + 5.0*target->GetDefinition()->GetPDGWidth() );
159  }
160 
161  G4double M0target2 = M0target * M0target;
162  G4double TargetDiffStateMinMass = theParameters->GetTarMinDiffMass();
163  G4double TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass();
164  if ( M0target > TargetDiffStateMinMass ) {
165  TargetDiffStateMinMass = M0target + 220.0*MeV;
166  TargetNonDiffStateMinMass = M0target + 220.0*MeV;
167  if ( absTargetPDGcode > 3000 ) { // Strange baryon
168  TargetDiffStateMinMass += 140.0*MeV;
169  TargetNonDiffStateMinMass += 140.0*MeV;
170  }
171  };
172 
173  #ifdef debugFTFexictation
174  G4cout << "Proj Targ PDGcodes " << ProjectilePDGcode << " " << TargetPDGcode << G4endl
175  << "M0projectile Y " << M0projectile << " " << ProjectileRapidity << G4endl;
176  //G4cout << "M0target Y " << M0target << " " << TargetRapidity << G4endl;
177  G4cout << "Pproj " << Pprojectile << G4endl << "Ptarget " << Ptarget << G4endl;
178  #endif
179 
180  G4double AveragePt2 = theParameters->GetAveragePt2();
181  //G4double ProbLogDistrPrD = theParameters->GetProbLogDistrPrD();
182  G4double ProbLogDistr = theParameters->GetProbLogDistr();
183  G4double SumMasses = M0projectile + M0target; // + 220.0*MeV;
184 
185  // Transform momenta to cms and then rotate parallel to z axis;
186  G4LorentzRotation toCms( -1 * Psum.boostVector() );
187 
188  G4LorentzVector Ptmp = toCms * Pprojectile;
189  if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
190 
191  toCms.rotateZ( -1*Ptmp.phi() );
192  toCms.rotateY( -1*Ptmp.theta() );
193  G4LorentzRotation toLab(toCms.inverse());
194  Pprojectile.transform( toCms );
195  Ptarget.transform( toCms );
196 
197  G4double PZcms2, PZcms;
198 
199  #ifdef debugFTFexictation
200  G4cout << "SqrtS " << SqrtS << G4endl << "M0pr M0tr SumM+220 " << M0projectile << " "
201  << M0target << " " << SumMasses << G4endl;
202  #endif
203 
204  if ( SqrtS < M0projectile + M0target ) return false;
205  if ( SqrtS < SumMasses ) return false;
206  // The model cannot work at low energy
207 
208  PZcms2 = ( S*S + M0projectile2*M0projectile2 + M0target2*M0target2
209  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
210 
211  #ifdef debugFTFexictation
212  G4cout << "PZcms2 after PutOnMassShell " << PZcms2 << G4endl;
213  #endif
214 
215  if ( PZcms2 < 0 ) return false;
216  // It can be in an interaction with off-shell nuclear nucleon
217 
218  PZcms = std::sqrt( PZcms2 );
219  if ( PutOnMassShell ) {
220  if ( Pprojectile.z() > 0.0 ) {
221  Pprojectile.setPz( PZcms );
222  Ptarget.setPz( -PZcms );
223  } else {
224  Pprojectile.setPz( -PZcms );
225  Ptarget.setPz( PZcms );
226  }
227  Pprojectile.setE( std::sqrt( M0projectile2 +
228  Pprojectile.x()*Pprojectile.x() +
229  Pprojectile.y()*Pprojectile.y() +
230  PZcms2 ) );
231  Ptarget.setE( std::sqrt( M0target2 +
232  Ptarget.x()*Ptarget.x() +
233  Ptarget.y()*Ptarget.y() +
234  PZcms2 ) );
235  }
236 
237  G4double maxPtSquare( 0.0 ); // = PZcms2;
238 
239  //Uzhi_QEnex = 0;
240  //Uzhi_QEexc = 0;
241  //Uzhi_targetdiffraction = 0;
242  //Uzhi_projectilediffraction = 0;
243  //Uzhi_nondiffraction = 0;
244  //G4int UzhiPrD( 0 ), UzhiTrD( 0 ), UzhiND( 0 );
245 
246  #ifdef debugFTFexictation
247  G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << M0projectile
248  << " " << ProjectileDiffStateMinMass << " " << ProjectileNonDiffStateMinMass << G4endl
249  << "Targ M0 Mdif Mndif " << M0target << " " << TargetDiffStateMinMass << " "
250  << TargetNonDiffStateMinMass << G4endl << "SqrtS " << SqrtS << G4endl
251  << "Proj CMS " << Pprojectile << G4endl << "Targ CMS " << Ptarget << G4endl;
252  #endif
253 
254  // Charge exchange can be possible
255  // Getting the values needed for exchange
256  // Check for possible quark exchange
257  G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity );
258  G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity )*theParameters->GetProcProb( 4, ProjectileRapidity );
259  G4double ProbProjectileDiffraction = theParameters->GetProcProb( 2, ProjectileRapidity );
260  G4double ProbTargetDiffraction = theParameters->GetProcProb( 3, ProjectileRapidity );
261 
262  //QeNoExc = 1.0;
263  //QeExc = 0.0;
264  //ProbProjectileDiffraction = 0.0;
265  //ProbTargetDiffraction = 0.0;
266  //AveragePt2 = 0.0;
267 
268  if ( QeNoExc + QeExc + ProbProjectileDiffraction + ProbTargetDiffraction > 1.0 ) {
269  QeNoExc = 1.0 - QeExc - ProbProjectileDiffraction - ProbTargetDiffraction;
270  }
271  //QeNoExc = 0.0;
272  G4double ProbExc( 0.0 );
273  if ( QeExc + QeNoExc != 0.0 ) ProbExc = QeExc/(QeExc + QeNoExc);
274  G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
275  G4double DeltaMass = G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass();
276 
277  //ProbProjectileDiffraction = 0.5;
278  //ProbTargetDiffraction = 0.5;
279  G4double ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction;
280 
281  #ifdef debugFTFexictation
282  G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " " << ProbProjectileDiffraction
283  << " " << ProbTargetDiffraction << G4endl
284  << "ProjectileRapidity " << ProjectileRapidity << G4endl;
285  //G4int Uzhi; G4cin >> Uzhi;
286  #endif
287 
288  G4ParticleDefinition* TestParticle( 0 );
289  G4double MtestPr( 0.0 ), MtestTr( 0.0 );
290 
291  if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
292  ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
293  ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
294  }
295 
296  if ( G4UniformRand() < QeExc + QeNoExc ) {
297 
298  #ifdef debugFTFexictation
299  G4cout << "Q exchange --------------------------" << G4endl;
300  #endif
301 
302  G4int NewProjCode( 0 ), NewTargCode( 0 );
303  G4int ProjQ1( 0 ), ProjQ2( 0 ), ProjQ3( 0 );
304 
305  // Projectile unpacking
306  if ( absProjectilePDGcode < 1000 ) { // projectile is meson
307  UnpackMeson( ProjectilePDGcode, ProjQ1, ProjQ2 );
308  } else { // projectile is baryon
309  UnpackBaryon( ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
310  }
311 
312  // Target unpacking
313  G4int TargQ1( 0 ), TargQ2( 0 ), TargQ3( 0 );
314  UnpackBaryon( TargetPDGcode, TargQ1, TargQ2, TargQ3 );
315 
316  #ifdef debugFTFexictation
317  G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
318  << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
319  #endif
320 
321  // Sampling of exchanged quarks
322  G4int ProjExchangeQ( 0 );
323  G4int TargExchangeQ( 0 );
324 
325  if ( absProjectilePDGcode < 1000 ) { // projectile is meson
326 
327  if ( ProjQ1 > 0 ) { // ProjQ1 is quark
328  ProjExchangeQ = ProjQ1;
329  // Exchange of non-identical quarks is allowed
330  G4int NpossibleStates = 3;
331  NpossibleStates = 0;
332  if ( ProjQ1 != TargQ1 ) NpossibleStates++;
333  if ( ProjQ1 != TargQ2 ) NpossibleStates++;
334  if ( ProjQ1 != TargQ3 ) NpossibleStates++;
335 
336  G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
337  //G4cout<<"NpossibleStates Nsampled "<<NpossibleStates<<" "<<Nsampled<<G4endl;
338 
339  NpossibleStates=0;
340  if ( ProjQ1 != TargQ1 ) {
341  NpossibleStates++;
342  if ( NpossibleStates == Nsampled ) {
343  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ1 = TargExchangeQ;
344  }
345  }
346  if ( ProjQ1 != TargQ2 ) {
347  NpossibleStates++;
348  if ( NpossibleStates == Nsampled ) {
349  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ1 = TargExchangeQ;
350  }
351  }
352  if ( ProjQ1 != TargQ3 ) {
353  NpossibleStates++;
354  if ( NpossibleStates == Nsampled ) {
355  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ1 = TargExchangeQ;
356  }
357  }
358 
359  //if ( Nsampled == 1 ) { TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ1 = TargExchangeQ; }
360  //else if ( Nsampled == 2 ) { TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ1 = TargExchangeQ; }
361  //else { TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ1 = TargExchangeQ; }
362 
363  } else { // ProjQ2 is quark
364  ProjExchangeQ = ProjQ2;
365  // Exchange of non-identical quarks is allowed
366  G4int NpossibleStates = 3;
367  NpossibleStates = 0;
368  if ( ProjQ2 != TargQ1 ) NpossibleStates++;
369  if ( ProjQ2 != TargQ2 ) NpossibleStates++;
370  if ( ProjQ2 != TargQ3 ) NpossibleStates++;
371 
372  G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
373 
374  NpossibleStates = 0;
375  if ( ProjQ2 != TargQ1 ) {
376  NpossibleStates++;
377  if ( NpossibleStates == Nsampled ) {
378  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ2 = TargExchangeQ;
379  }
380  }
381  if ( ProjQ2 != TargQ2 ) {
382  NpossibleStates++;
383  if ( NpossibleStates == Nsampled ) {
384  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ2 = TargExchangeQ;
385  }
386  }
387  if ( ProjQ2 != TargQ3 ) {
388  NpossibleStates++;
389  if ( NpossibleStates == Nsampled ) {
390  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ2 = TargExchangeQ;
391  }
392  }
393 
394  //if ( Nsampled == 1 ) { TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ2 = TargExchangeQ; }
395  //else if ( Nsampled == 2 ) { TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ2 = TargExchangeQ; }
396  //else { TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ2 = TargExchangeQ; }
397 
398  } // End of if ( ProjQ1 > 0 )
399 
400  #ifdef debugFTFexictation
401  G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
402  #endif
403 
404  G4int aProjQ1 = std::abs( ProjQ1 );
405  G4int aProjQ2 = std::abs( ProjQ2 );
406 
407  G4bool ProjExcited = false;
408 
409  G4int attempts = 0;
410  while ( attempts < 50 ) { /* Loop checking, 10.08.2015, A.Ribon */
411  // Determination of a new projectile ID which garanty energy-momentum conservation
412  attempts++;
413  G4double Ksi = G4UniformRand();
414  if ( aProjQ1 == aProjQ2 ) {
415  if ( aProjQ1 != 3 ) {
416  NewProjCode = 111; // Pi0-meson
417  if ( Ksi < 0.5 ) {
418  NewProjCode = 221; // Eta -meson
419  if ( Ksi < 0.25 ) {
420  NewProjCode = 331; // Eta'-meson
421  }
422  }
423  } else {
424  NewProjCode = 221; // Eta -meson
425  if ( Ksi < 0.5 ) {
426  NewProjCode = 331; // Eta'-meson
427  }
428  }
429  } else {
430  if ( aProjQ1 > aProjQ2 ) {
431  NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
432  } else {
433  NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
434  }
435  }
436 
437  #ifdef debugFTFexictation
438  G4cout << "NewProjCode " << NewProjCode << G4endl;
439  #endif
440 
441  ProjExcited = false;
442  if ( G4UniformRand() < 0.5 ) {
443  NewProjCode += 2; // Excited meson
444  ProjExcited = true;
445  }
446  //if ( aProjQ1 != aProjQ2 ) NewProjCode *= ( ProjectilePDGcode / absProjectilePDGcode );
447  G4int Qquarks=0;
448  if ( aProjQ1 == 1 ) { Qquarks -= ProjQ1; }
449  else if ( aProjQ1 == 2 ) { Qquarks += ProjQ1; }
450  else { Qquarks -= ProjQ1/aProjQ1; }
451 
452  if ( aProjQ2 == 1 ) { Qquarks -= ProjQ2; }
453  else if ( aProjQ2 == 2 ) { Qquarks += ProjQ2; }
454  else { Qquarks -= ProjQ2/aProjQ2; }
455 
456  if ( Qquarks < 0 ) NewProjCode *= -1;
457 
458  #ifdef debugFTFexictation
459  G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
460  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
461  G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquarks<<G4endl;
462  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
463  #endif
464 
465  // Proj
466  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
467  if ( ! TestParticle ) continue;
468 
469  //MminProjectile = TestParticle->GetPDGMass();
470  MminProjectile = BrW.GetMinimumMass( TestParticle );
471 
472  if ( SqrtS - M0target < MminProjectile ) continue;
473 
474  MtestPr = BrW.SampleMass( TestParticle,
475  TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth() );
476  //G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode )->GetPDGMass();
477 
478  #ifdef debugFTFexictation
479  G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()<< G4endl;
480  G4cout << "MtestPart MtestPart0 "<<MtestPr<<" "<<TestParticle->GetPDGMass()<<G4endl;
481  G4cout << "M0projectile projectile PDGMass " << M0projectile << " "
482  << projectile->GetDefinition()->GetPDGMass() << G4endl;
483  #endif
484 
485  // Targ
486  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
487 
488  #ifdef debugFTFexictation
489  G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
490  << "NewTargCode " << NewTargCode << G4endl;
491  #endif
492 
493  if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) { // Lambda or Sigma0 ?
494  if ( G4UniformRand() < 0.5 ) {
495  NewTargCode += 2;
496  } else {
497  if ( G4UniformRand() < 0.75 ) NewTargCode = 3122;
498  }
499  } else if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
500  NewTargCode += 2; ProjExcited = true; //Create Delta isobar
501  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
502  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
503  NewTargCode += 2; ProjExcited = true; // Save Delta isobar
504  } else { // De-excite initial Delta isobar
505  }
506  } else if ( ! ProjExcited &&
507  G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
508  SqrtS > M0projectile + DeltaMass ) { // Create Delta isobar
509  NewTargCode += 2; // Save initial nucleon
510  } else {
511  }
512 
513  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
514 
515  if ( ! TestParticle ) continue;
516 
517  #ifdef debugFTFexictation
518  G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
519  #endif
520 
521  MminTarget = BrW.GetMinimumMass( TestParticle );
522 
523  if ( SqrtS - MtestPr < MminTarget ) continue;
524 
525  MtestTr = BrW.SampleMass( TestParticle,
526  TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth() );
527 
528  if ( SqrtS > MtestPr + MtestTr ) break;
529 
530  } // End of while ( attempts < 50 )
531 
532  if ( attempts >= 50 ) return false;
533 
534  /*
535  if ( MtestPr > Pprojectile.mag() ) {
536  M0projectile = MtestPr;
537  } else {
538  if ( std::abs( MtestPr - M0projectile ) // projectile->GetDefinition()->GetPDGMass() )
539  < 140.0*MeV ) {
540  M0projectile = MtestPr;
541  }
542  }
543  */
544 
545  if ( MtestPr >= Pprojectile.mag() ) {
546  M0projectile = MtestPr;
547  } else if ( projectile->GetStatus() != 0 ) {
548  M0projectile = MtestPr;
549  }
550 
551  #ifdef debugFTFexictation
552  G4cout << "M0projectile After check " << M0projectile << G4endl;
553  #endif
554 
555  M0projectile2 = M0projectile * M0projectile;
556  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
557  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
558 
559  if ( MtestTr >= Ptarget.mag() ) {
560  M0target = MtestTr;
561  } else if ( target->GetStatus() != 0 ) {
562  M0target = MtestTr;
563  }
564 
565  /*
566  M0target = MtestTr;
567  if ( MtestTr > Ptarget.mag() ) {
568  M0target = MtestTr;
569  } else {
570  if ( std::abs( MtestTr - M0target ) // target->GetDefinition()->GetPDGMass() )
571  < 140.0*MeV ) {
572  M0target = MtestTr;
573  }
574  }
575  */
576 
577  M0target2 = M0target * M0target;
578 
579  #ifdef debugFTFexictation
580  G4cout << "New targ M0 M0^2 " << M0target << " " << M0target2 << G4endl;
581  #endif
582 
583  TargetDiffStateMinMass = M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
584  TargetNonDiffStateMinMass = M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
585 
586  } else { // of the if ( absProjectilePDGcode < 1000 ) ;
587  // The projectile is baryon now
588 
589  G4double Same = theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0.
590  //G4bool ProjDeltaHasCreated( false );
591  //G4bool TargDeltaHasCreated( false );
592 
593  G4double Ksi = G4UniformRand();
594  if ( G4UniformRand() < 0.5 ) { // Sampling exchange quark from proj. or targ.
595 
596  // Sampling exchanged quark from the projectile
597  if ( Ksi < 0.333333 ) {
598  ProjExchangeQ = ProjQ1;
599  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
600  ProjExchangeQ = ProjQ2;
601  } else {
602  ProjExchangeQ = ProjQ3;
603  }
604 
605  if ( ProjExchangeQ != TargQ1 || G4UniformRand() < Same ) {
606  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
607  } else {
608  if ( ProjExchangeQ != TargQ2 || G4UniformRand() < Same ) {
609  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
610  } else {
611  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
612  }
613  }
614 
615  #ifdef debugFTFexictation
616  G4cout << "Exchange Qs Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
617  #endif
618 
619  if ( Ksi < 0.333333 ) {
620  ProjQ1 = ProjExchangeQ;
621  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
622  ProjQ2 = ProjExchangeQ;
623  } else {
624  ProjQ3 = ProjExchangeQ;
625  }
626 
627  } else { // Sampling exchanged quark from the target
628 
629  if ( Ksi < 0.333333 ) {
630  TargExchangeQ = TargQ1;
631  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
632  TargExchangeQ = TargQ2;
633  } else {
634  TargExchangeQ = TargQ3;
635  }
636  if ( TargExchangeQ != ProjQ1 || G4UniformRand() < Same ) {
637  ProjExchangeQ = ProjQ1; ProjQ1 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
638  } else {
639  if ( TargExchangeQ != ProjQ2 || G4UniformRand() < Same ) {
640  ProjExchangeQ = ProjQ2; ProjQ2 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
641  } else {
642  ProjExchangeQ = ProjQ3; ProjQ3 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
643  }
644  }
645 
646  if ( Ksi < 0.333333 ) {
647  TargQ1 = TargExchangeQ;
648  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
649  TargQ2 = TargExchangeQ;
650  } else {
651  TargQ3 = TargExchangeQ;
652  }
653 
654  } // End of quark sampling for the baryons
655 
656  NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
657  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
658 
659  G4int attempts = 0;
660  while ( attempts < 50 ) { /* Loop checking, 10.08.2015, A.Ribon */
661  // Determination of a new projectile ID which garanty energy-momentum conservation
662  attempts++;
663 
664  if ( ProjQ1 == ProjQ2 && ProjQ1 == ProjQ3 ) {
665  NewProjCode += 2; // ProjDeltaHasCreated = true;
666  } else if ( projectile->GetDefinition()->GetPDGiIsospin() == 3 ) { // Projectile was Delta
667  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
668  NewProjCode += 2; //ProjDeltaHasCreated = true;
669  } else {
670  NewProjCode += 0; //ProjDeltaHasCreated = false;
671  }
672  } else { // Projectile was Nucleon
673  if ( G4UniformRand() < DeltaProbAtQuarkExchange && SqrtS > DeltaMass + M0target ) {
674  NewProjCode += 2; //ProjDeltaHasCreated = true;
675  } else {
676  NewProjCode += 0; //ProjDeltaHasCreated = false;
677  }
678  }
679 
680  if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
681  NewTargCode += 2; //TargDeltaHasCreated = true;
682  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Target was Delta
683  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
684  NewTargCode += 2; //TargDeltaHasCreated = true;
685  } else {
686  NewTargCode += 0; //TargDeltaHasCreated = false;
687  }
688  } else { // Target was Nucleon
689  if ( G4UniformRand() < DeltaProbAtQuarkExchange && SqrtS > M0projectile + DeltaMass ) {
690  NewTargCode += 2; //TargDeltaHasCreated = true;
691  } else {
692  NewTargCode += 0; //TargDeltaHasCreated = false;
693  }
694  }
695 
696  #ifdef debugFTFexictation
697  G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
698  //G4int Uzhi; G4cin >> Uzhi;
699  #endif
700 
701  if ( absProjectilePDGcode == NewProjCode && absTargetPDGcode == NewTargCode ) {
702  } // Nothing was changed! It is not right!?
703 
704  // Forming baryons
705  /*
706  if ( G4UniformRand() > 0.5 ) {
707  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
708  } else {
709  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
710  }
711 
712  if ( ProjDeltaHasCreated ) {
713  if ( G4UniformRand() > 0.5 ) {
714  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
715  } else {
716  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
717  }
718  }
719 
720  if ( TargDeltaHasCreated ) {
721  if ( G4UniformRand() > 0.5 ) {
722  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
723  } else {
724  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
725  }
726  }
727  */
728 
729  // Proj
730  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
731  if ( ! TestParticle ) continue;
732 
733  MminProjectile = BrW.GetMinimumMass( TestParticle );
734 
735  if ( SqrtS - M0target < MminProjectile ) continue;
736 
737  MtestPr = BrW.SampleMass( TestParticle,
738  TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth() );
739 
740  // Targ
741  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
742  if ( ! TestParticle ) continue;
743 
744  MminTarget = BrW.GetMinimumMass( TestParticle );
745 
746  if ( SqrtS - MtestPr < MminTarget ) continue;
747 
748  MtestTr = BrW.SampleMass( TestParticle,
749  TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth() );
750 
751  if ( SqrtS > MtestPr + MtestTr ) break;
752  } // End of while ( attempts < 50 )
753 
754  if ( attempts >= 50 ) return false;
755 
756  if ( MtestPr >= Pprojectile.mag() ) {
757  M0projectile = MtestPr;
758  } else if ( projectile->GetStatus() != 0 ) {
759  M0projectile = MtestPr;
760  }
761  M0projectile2 = M0projectile * M0projectile;
762  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
763  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
764 
765  if ( MtestTr >= Ptarget.mag() ) {
766  M0target = MtestTr;
767  } else if ( target->GetStatus() != 0 ) {
768  M0target = MtestTr;
769  }
770  M0target2 = M0target * M0target;
771  TargetDiffStateMinMass = M0target + 220.0*MeV; //220 MeV=m_pi+80 MeV;
772  TargetNonDiffStateMinMass = M0target + 220.0*MeV; //220 MeV=m_pi+80 MeV;
773 
774  } // End of if ( absProjectilePDGcode < 1000 )
775 
776  // If we assume that final state hadrons after the charge exchange will be
777  // in the ground states, we have to put
778  if ( SqrtS < M0projectile + M0target ) return false;
779 
780  PZcms2 = ( S*S + M0projectile2*M0projectile2 + M0target2*M0target2
781  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
782 
783  #ifdef debugFTFexictation
784  G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
785  << "At the end// NewTargCode " << NewTargCode << G4endl
786  << "M0pr M0tr SqS " << M0projectile << " " << M0target << " " << SqrtS << G4endl
787  << "M0pr2 M0tr2 SqS " << M0projectile2 << " " << M0target2 << " " << SqrtS << G4endl
788  << "PZcms2 after the change " << PZcms2 << G4endl << G4endl;
789  #endif
790 
791  if ( PZcms2 < 0 ) return false; // It can be if energy is not sufficient for Delta
792 
793  projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
794  target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
795 
796  PZcms = std::sqrt( PZcms2 );
797  Pprojectile.setPz( PZcms );
798  Pprojectile.setE( std::sqrt( M0projectile2 + PZcms2 ) );
799  Ptarget.setPz( -PZcms );
800  Ptarget.setE( std::sqrt( M0target2 + PZcms2 ) );
801 
802  if ( projectile->GetStatus() != 0 ) projectile->SetStatus( 2 );
803  if ( target->GetStatus() != 0 ) target->SetStatus( 2 );
804 
805  #ifdef debugFTFexictation
806  G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << Pprojectile << G4endl << Ptarget
807  << G4endl << Pprojectile + Ptarget << G4endl;
808  #endif
809 
810  // Check for possible excitation of the participants
811  if ( SqrtS < M0projectile + TargetDiffStateMinMass ||
812  SqrtS < ProjectileDiffStateMinMass + M0target ||
813  ProbOfDiffraction == 0.0 ) ProbExc = 0.0;
814 
815  if ( G4UniformRand() > ProbExc ) { // Make elastic scattering
816 
817  #ifdef debugFTFexictation
818  G4cout << "Make elastic scattering of new hadrons" << G4endl;
819  #endif
820 
821  Pprojectile.transform( toLab );
822  Ptarget.transform( toLab );
823 
824  projectile->Set4Momentum( Pprojectile );
825  target->Set4Momentum( Ptarget );
826 
827  G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
828 
829  #ifdef debugFTFexictation
830  G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
831  << G4endl << projectile->Get4Momentum() << G4endl << target->Get4Momentum() << G4endl
832  << projectile->Get4Momentum() + target->Get4Momentum() << " " << (projectile->Get4Momentum() + target->Get4Momentum()).mag() << G4endl;
833  #endif
834 
835  //Uzhi_QEnex++;
836  return Result;
837  }
838  //Uzhi_QEexc++;
839 
840  #ifdef debugFTFexictation
841  G4cout << "Make excitation of new hadrons" << G4endl;
842  #endif
843 
844  // Redefinition of ProbOfDiffraction because the probabilities are changed due to quark exchange
845 
846  ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction;
847  if ( ProbOfDiffraction != 0.0 ) {
848  ProbProjectileDiffraction /= ProbOfDiffraction;
849  ProbTargetDiffraction /= ProbOfDiffraction;
850  }
851 
852  //Uzhi_QEnex++;
853  } // End of if ( G4UniformRand() < QeExc + QeNoExc ) , i.e. of the charge exchange part
854 
855  ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction;
856 
857  #ifdef debugFTFexictation
858  G4cout << "Excitation --------------------" << G4endl
859  << "Proj M0 MdMin MndMin " << M0projectile << " " << ProjectileDiffStateMinMass << " "
860  << ProjectileNonDiffStateMinMass << G4endl
861  << "Targ M0 MdMin MndMin " << M0target << " " << TargetDiffStateMinMass << " "
862  << TargetNonDiffStateMinMass << G4endl << "SqrtS " << SqrtS << G4endl
863  << "Prob: ProjDiff TargDiff + Sum " << ProbProjectileDiffraction << " "
864  << ProbTargetDiffraction << " " << ProbOfDiffraction << G4endl;
865  #endif
866 
867  if ( ProbOfDiffraction != 0.0 ) {
868  ProbProjectileDiffraction /= ProbOfDiffraction;
869  } else {
870  ProbProjectileDiffraction = 0.0;
871  }
872 
873  #ifdef debugFTFexictation
874  G4cout << "Prob: ProjDiff TargDiff + Sum " << ProbProjectileDiffraction << " "
875  << ProbTargetDiffraction << " " << ProbOfDiffraction << G4endl;
876  #endif
877 
878  G4double ProjectileDiffStateMinMass2 = sqr( ProjectileDiffStateMinMass );
879  G4double ProjectileNonDiffStateMinMass2 = sqr( ProjectileNonDiffStateMinMass );
880  G4double TargetDiffStateMinMass2 = sqr( TargetDiffStateMinMass );
881  G4double TargetNonDiffStateMinMass2 = sqr( TargetNonDiffStateMinMass );
882 
883  G4double Pt2;
884  G4double ProjMassT2, ProjMassT;
885  G4double TargMassT2, TargMassT;
886  G4double PMinusMin, PMinusMax;
887  //G4double PPlusMin , PPlusMax;
888  G4double TPlusMin, TPlusMax;
889  G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
890  G4LorentzVector Qmomentum;
891  G4double Qminus, Qplus;
892  G4int whilecount = 0;
893 
894  // Choose a process
895  if ( G4UniformRand() < ProbOfDiffraction ) {
896 
897  if ( G4UniformRand() < ProbProjectileDiffraction ) { // projectile diffraction
898 
899  #ifdef debugFTFexictation
900  G4cout << "projectile diffraction" << G4endl;
901  #endif
902 
903  //UzhiPrD++;
904 
905  do {
906 
907  //Uzhi_projectilediffraction = 1;
908  //Uzhi_targetdiffraction = 0;
909  //Uzhi_Mx2 = 1.0;
910 
911  // Generate pt and mass of projectile
912 
913  whilecount++;
914  if ( whilecount > 1000 ) {
915  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
916  return false; // Ignore this interaction
917  };
918 
919  // Check that the interaction is possible
920  ProjMassT2 = ProjectileDiffStateMinMass2;
921  ProjMassT = ProjectileDiffStateMinMass;
922  TargMassT2 = M0target2;
923  TargMassT = M0target;
924  if ( SqrtS < ProjMassT + TargMassT ) return false;
925 
926  PZcms2 =( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
927  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
928 
929  if ( PZcms2 < 0 ) return false;
930 
931  maxPtSquare = PZcms2;
932 
933  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
934 
935  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
936  ProjMassT2 = ProjectileDiffStateMinMass2 + Pt2;
937  ProjMassT = std::sqrt( ProjMassT2 );
938  TargMassT2 = M0target2 + Pt2;
939  TargMassT = std::sqrt( TargMassT2 );
940  if ( SqrtS < ProjMassT + TargMassT ) continue;
941 
942  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
943  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
944 
945  if ( PZcms2 < 0 ) continue;
946 
947  PZcms = std::sqrt( PZcms2 );
948  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
949  PMinusMax = SqrtS - TargMassT;
950 
951  //PMinusNew = PMinusMax; //ChooseP( PMinusMin, PMinusMax );
952  PMinusNew = ChooseP( PMinusMin, PMinusMax );
953 
954  TMinusNew = SqrtS - PMinusNew;
955  Qminus = Ptarget.minus() - TMinusNew;
956  TPlusNew = TargMassT2 / TMinusNew;
957  Qplus = Ptarget.plus() - TPlusNew;
958  Qmomentum.setPz( (Qplus - Qminus)/2 );
959  Qmomentum.setE( (Qplus + Qminus)/2 );
960 
961  } while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
962  // Repeat the sampling because there was not any excitation
963 
964  //projectile->SetStatus( 1*projectile->GetStatus() );
965 
966  if ( projectile->GetStatus() == 2 ) projectile->SetStatus( 1 );
967  if ( target->GetStatus() == 1 && target->GetSoftCollisionCount() == 0 ) target->SetStatus( 2 );
968 
969  } else { // Target diffraction
970 
971  #ifdef debugFTFexictation
972  G4cout << "Target diffraction" << G4endl;
973  #endif
974 
975  //UzhiTrD++;
976 
977  do {
978 
979  //Uzhi_projectilediffraction = 0;
980  //Uzhi_targetdiffraction = 1;
981  //Uzhi_Mx2 = 1.0;
982 
983  // Generate pt and target mass
984 
985  whilecount++;
986  if ( whilecount > 1000 ) {
987  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
988  return false; // Ignore this interaction
989  };
990 
991  // Check that the interaction is possible
992  ProjMassT2 = M0projectile2;
993  ProjMassT = M0projectile;
994 
995  TargMassT2 = TargetDiffStateMinMass2;
996  TargMassT = TargetDiffStateMinMass;
997 
998  if ( SqrtS < ProjMassT + TargMassT ) return false;
999 
1000  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1001  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1002 
1003  if ( PZcms2 < 0 ) return false;
1004 
1005  maxPtSquare = PZcms2;
1006 
1007  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
1008 
1009  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
1010  ProjMassT2 = M0projectile2 + Pt2;
1011  ProjMassT = std::sqrt( ProjMassT2 );
1012  TargMassT2 = TargetDiffStateMinMass2 + Pt2;
1013  TargMassT = std::sqrt( TargMassT2 );
1014  if ( SqrtS < ProjMassT + TargMassT ) continue;
1015 
1016  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1017  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1018 
1019  if ( PZcms2 < 0 ) continue;
1020 
1021  PZcms = std::sqrt( PZcms2 );
1022  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
1023  //TPlusMax = std::sqrt( TargMassT2 + PZcms2 ) + PZcms;
1024  TPlusMax = SqrtS - ProjMassT;
1025 
1026  TPlusNew = ChooseP( TPlusMin, TPlusMax );
1027  //TPlusNew = TPlusMin;
1028 
1029  PPlusNew = SqrtS - TPlusNew;
1030  Qplus = PPlusNew - Pprojectile.plus();
1031  PMinusNew = ProjMassT2 / PPlusNew;
1032  Qminus = PMinusNew - Pprojectile.minus();
1033  Qmomentum.setPz( (Qplus - Qminus)/2 );
1034  Qmomentum.setE( (Qplus + Qminus)/2 );
1035 
1036  } while ( ( Ptarget - Qmomentum ).mag2() < TargetDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
1037  // Repeat the sampling because there was not any excitation
1038 
1039  //target->SetStatus( 1*target->GetStatus() );
1040 
1041  if ( projectile->GetStatus() == 1 && projectile->GetSoftCollisionCount() == 0 ) projectile->SetStatus( 2 );
1042  if ( target->GetStatus() == 2 ) target->SetStatus( 1 );
1043 
1044  } // End of if ( G4UniformRand() < ProbProjectileDiffraction )
1045 
1046  } else { // Non-diffraction process
1047 
1048  #ifdef debugFTFexictation
1049  G4cout << "Non-diffraction process" << G4endl;
1050  #endif
1051 
1052  //UzhiND++;
1053  //Uzhi_QEnex++;
1054  do {
1055 
1056  //Uzhi_projectilediffraction = 0;
1057  //Uzhi_targetdiffraction = 0;
1058  //Uzhi_Mx2 = 1.0;
1059 
1060  // Generate pt and masses
1061 
1062  whilecount++;
1063  if ( whilecount > 1000 ) {
1064  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1065  return false; // Ignore this interaction
1066  };
1067 
1068  // Check that the interaction is possible
1069  ProjMassT2 = ProjectileNonDiffStateMinMass2;
1070  ProjMassT = ProjectileNonDiffStateMinMass;
1071  TargMassT2 = TargetNonDiffStateMinMass2;
1072  TargMassT = TargetNonDiffStateMinMass;
1073  if ( SqrtS < ProjMassT + TargMassT ) return false;
1074 
1075  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1076  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1077 
1078  if ( PZcms2 < 0 ) return false;
1079 
1080  maxPtSquare = PZcms2;
1081 
1082  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
1083 
1084  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
1085  ProjMassT2 = ProjectileNonDiffStateMinMass2 + Pt2;
1086  ProjMassT = std::sqrt( ProjMassT2 );
1087  TargMassT2 = TargetNonDiffStateMinMass2 + Pt2;
1088  TargMassT = std::sqrt( TargMassT2 );
1089  if ( SqrtS < ProjMassT + TargMassT ) continue;
1090 
1091  PZcms2 =( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1092  -2.0*S*ProjMassT2 - 2.0*S*TargMassT2 -2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1093 
1094  if ( PZcms2 < 0 ) continue;
1095 
1096  PZcms = std::sqrt( PZcms2 );
1097  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
1098  //PMinusMax = std::sqrt( ProjMassT2 + PZcms2 ) + PZcms;
1099  PMinusMax = SqrtS - TargMassT;
1100 
1101  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
1102  //TPlusMax = std::sqrt( TargMassT2 + PZcms2 ) + PZcms;
1103  TPlusMax = SqrtS - ProjMassT;
1104 
1105  /*
1106  if ( G4UniformRand() < ProbLogDistrPrD ) {
1107  PMinusNew = ChooseP( PMinusMin, PMinusMax );
1108  } else {
1109  PMinusNew = ( PMinusMax - PMinusMin )*G4UniformRand() + PMinusMin;
1110  }
1111  //Qminus = PMinusNew - Pprojectile.minus();
1112  */
1113 
1114  //TPlusMax = SqrtS - PMinusNew;
1115 
1116  if ( G4UniformRand() < ProbLogDistr ) {
1117  PMinusNew = ChooseP( PMinusMin, PMinusMax );
1118  TPlusNew = ChooseP( TPlusMin, TPlusMax );
1119  } else {
1120  PMinusNew = ( PMinusMax - PMinusMin )*G4UniformRand() + PMinusMin;
1121  TPlusNew = ( TPlusMax - TPlusMin )*G4UniformRand() + TPlusMin;
1122  }
1123 
1124  Qminus = PMinusNew - Pprojectile.minus();
1125 
1126  Qplus = -( TPlusNew - Ptarget.plus() );
1127  Qmomentum.setPz( (Qplus - Qminus)/2 );
1128  Qmomentum.setE( (Qplus + Qminus)/2 );
1129 
1130  #ifdef debugFTFexictation
1131  G4cout << ( Pprojectile + Qmomentum ).mag2() << " " << ProjectileNonDiffStateMinMass2
1132  << G4endl << ( Ptarget - Qmomentum ).mag2() << " "
1133  << TargetNonDiffStateMinMass2 << G4endl;
1134  G4cout<<"To continue - enter any integer"<<G4endl;
1135  // G4int Uzhi; G4cin >> Uzhi;
1136  #endif
1137 
1138  } while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileNonDiffStateMinMass2 || //No double Diffraction
1139  ( Ptarget - Qmomentum ).mag2() < TargetNonDiffStateMinMass2 || // ); //
1140  ( Pprojectile + Qmomentum ).pz() < 0.); /* Loop checking, 10.08.2015, A.Ribon */
1141 
1142  projectile->SetStatus( 0*projectile->GetStatus() );
1143  target->SetStatus( 0*target->GetStatus() );
1144 
1145  } // End of if ( G4UniformRand() < ProbOfDiffraction )
1146 
1147  Pprojectile += Qmomentum;
1148  Ptarget -= Qmomentum;
1149 
1150  // Transform back and update SplitableHadron Participant.
1151  Pprojectile.transform( toLab );
1152  Ptarget.transform( toLab );
1153 
1154  // Calculation of the creation time
1155  //projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1156  //projectile->SetPosition( target->GetPosition() );
1157  // Creation time and position of target nucleon were determined in
1158  // ReggeonCascade() of G4FTFModel
1159  //
1160  //if ( Uzhi_projectilediffraction != 0 ) {
1161  // Uzhi_Mx2 = Pprojectile.mag2(); Uzhi_modT = ( target->Get4Momentum() - Ptarget ).mag2();
1162  //}
1163  //if ( Uzhi_targetdiffraction != 0 ) {
1164  // Uzhi_Mx2 = Ptarget.mag2(); Uzhi_modT = ( projectile->Get4Momentum() - Pprojectile ).mag2();
1165  //}
1166  //if ( Uzhi_QE != 0 ) {
1167  // Uzhi_projectilediffraction = 0;
1168  // Uzhi_targetdiffraction = 0;
1169  // Uzhi_Mx2 = 1.0;
1170  //}
1171 
1172  #ifdef debugFTFexictation
1173  G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl;
1174  #endif
1175 
1176  projectile->Set4Momentum( Pprojectile );
1177  target->Set4Momentum( Ptarget );
1178  projectile->IncrementCollisionCount( 1 );
1179  target->IncrementCollisionCount( 1 );
1180 
1181  //Uzhi_projectilediffraction = UzhiPrD;
1182  //Uzhi_targetdiffraction = UzhiTrD;
1183  //Uzhi_nondiffraction = UzhiND;
1184  //G4cout << Uzhi_projectilediffraction << " " << Uzhi_targetdiffraction << " "
1185  // << Uzhi_nondiffraction << G4endl;
1186 
1187  return true;
1188 }
1189 
1190 
1191 //============================================================================
1192 
1193 void G4DiffractiveExcitation::CreateStrings( G4VSplitableHadron* hadron,
1194  G4bool isProjectile,
1195  G4ExcitedString*& FirstString,
1196  G4ExcitedString*& SecondString,
1197  G4FTFParameters* theParameters ) const {
1198 
1199  //G4cout << "Create Strings SplitUp " << hadron << G4endl
1200  // << "Defin " << hadron->GetDefinition() << G4endl
1201  // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
1202 
1203  hadron->SplitUp();
1204 
1205  G4Parton* start = hadron->GetNextParton();
1206  if ( start == NULL ) {
1207  G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1208  FirstString = 0; SecondString = 0;
1209  return;
1210  }
1211 
1212  G4Parton* end = hadron->GetNextParton();
1213  if ( end == NULL ) {
1214  G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1215  FirstString = 0; SecondString = 0;
1216  return;
1217  }
1218 
1219  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1220  // << G4endl
1221  // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1222 
1223  G4LorentzVector Phadron = hadron->Get4Momentum();
1224  //G4cout << "String mom " << Phadron << G4endl;
1225  G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1226  G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1227  G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1228  G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1229  G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1230 
1231  G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1232  G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1233  //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1234 
1235  G4double Wmin( 0.0 );
1236  if ( isProjectile ) {
1237  Wmin = theParameters->GetProjMinDiffMass();
1238  } else {
1239  Wmin = theParameters->GetTarMinDiffMass();
1240  }
1241 
1242  G4double W = hadron->Get4Momentum().mag();
1243  //G4cout << "Wmin W " << Wmin << " " << W << G4endl;
1244  //G4int Uzhi; G4cin >> Uzhi;
1245  G4double W2 = W*W;
1246  G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); // x2( 0.0 )
1247  G4bool Kink = false;
1248 
1249  if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1250  end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1251  ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1252  end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1253  // Kinky strings are allowed only for qq-q strings;
1254  // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1255  // according to the analysis of Pbar P interactions
1256 
1257  if ( W > Wmin ) { // Kink is possible
1258  if ( hadron->GetStatus() == 0 ) {
1259  G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1260  if ( Pt2kink ) {
1261  Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );
1262  } else {
1263  Pt = 0.0;
1264  }
1265  } else {
1266  Pt = 0.0;
1267  }
1268 
1269  if ( Pt > 500.0*MeV ) {
1270  G4double Ymax = G4Log( W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1271  G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
1272  x1 = 1.0 - Pt/W * G4Exp( Y );
1273  x3 = 1.0 - Pt/W * G4Exp(-Y );
1274  //x2 = 2.0 - x1 - x3;
1275 
1276  G4double Mass_startQ = 650.0*MeV;
1277  if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV;
1278  if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV;
1279  if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV;
1280  G4double Mass_endQ = 650.0*MeV;
1281  if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV;
1282  if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV;
1283  if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV;
1284 
1285  G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1286  G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1287  G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
1288  if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1289  Kink = false;
1290  } else {
1291  G4double P_1 = std::sqrt( P2_1 );
1292  G4double P_2 = std::sqrt( P2_2 );
1293  G4double P_3 = std::sqrt( P2_3 );
1294  G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1295  G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1296  //Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated
1297 
1298  if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1299  Kink = false;
1300  } else {
1301  Kink = true;
1302  Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated
1303  Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1304  Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1305  Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1306  Pstart.setE( x3*W/2.0 );
1307  Pkink.setE( Pkink.vect().mag() );
1308  Pend.setE( x1*W/2.0 );
1309 
1310  G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1311  if ( Pkink.getZ() > 0.0 ) {
1312  if ( XkQ > 0.5 ) {
1313  PkinkQ1 = XkQ*Pkink;
1314  } else {
1315  PkinkQ1 = (1.0 - XkQ)*Pkink;
1316  }
1317  } else {
1318  if ( XkQ > 0.5 ) {
1319  PkinkQ1 = (1.0 - XkQ)*Pkink;
1320  } else {
1321  PkinkQ1 = XkQ*Pkink;
1322  }
1323  }
1324 
1325  PkinkQ2 = Pkink - PkinkQ1;
1326  // Minimizing Pt1^2+Pt3^2
1327  G4double Cos2Psi = ( sqr(x1) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1328  std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
1329  G4double Psi = std::acos( Cos2Psi );
1330 
1331  G4LorentzRotation Rotate;
1332  if ( isProjectile ) {
1333  Rotate.rotateY( Psi );
1334  } else {
1335  Rotate.rotateY( pi + Psi );
1336  }
1337  Rotate.rotateZ( twopi * G4UniformRand() );
1338  Pstart *= Rotate;
1339  Pkink *= Rotate;
1340  PkinkQ1 *= Rotate;
1341  PkinkQ2 *= Rotate;
1342  Pend *= Rotate;
1343  }
1344  } // End of if ( P2_1 <= 0.0 || P2_3 <= 0.0 )
1345  } // End of if ( Pt > 500.0*MeV )
1346  } // End of if ( W > Wmin ) : check for a kink
1347  } // end of qq-q string selection
1348 
1349  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetParticleSubType() << " "
1350  // << end->GetDefinition()->GetParticleSubType() << G4endl;
1351  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetPDGEncoding() << " "
1352  // << end->GetDefinition()->GetPDGEncoding() << G4endl;
1353  //G4int Uzhi; G4cin >> Uzhi;
1354 
1355  if ( Kink ) { // Kink is possible
1356 
1357  //G4cout << "Kink is sampled!" << G4endl;
1358  std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1359  theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1360 
1361  G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1362  for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1363  //G4cout << "Iq " << Iq << G4endl;
1364  if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1365  }
1366  //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1367  G4Parton* Gquark = new G4Parton( QuarkInGluon );
1368  G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1369  //G4cout << "Lorentz " << G4endl;
1370 
1371  G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1372  G4LorentzRotation toLab( toCMS.inverse() );
1373  //G4cout << "Pstart " << Pstart << G4endl;
1374  //G4cout << "Pend " << Pend << G4endl;
1375  //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1376  //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1377  //G4cout << "Pstart " << Pstart << G4endl<<G4endl;
1378 
1379  Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1380  PkinkQ1.transform( toLab );
1381  PkinkQ2.transform( toLab );
1382  Pend.transform( toLab ); end->Set4Momentum( Pend );
1383  //G4cout << "Pstart " << Pstart << G4endl;
1384  //G4cout << "Pend " << Pend << G4endl;
1385  //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1386  //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1387 
1388  //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1389  G4int absPDGcode = 1500; // 23 Dec
1390  if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1391  end->GetDefinition()->GetParticleSubType() == "quark" ) {
1392  absPDGcode = 110;
1393  }
1394  //G4cout << "absPDGcode " << absPDGcode << G4endl;
1395 
1396  if ( absPDGcode < 1000 ) { // meson
1397  if ( isProjectile ) { // Projectile
1398  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1399  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1400  SecondString = new G4ExcitedString( Gquark, start , +1 );
1401  Ganti_quark->Set4Momentum( PkinkQ1 );
1402  Gquark->Set4Momentum( PkinkQ2 );
1403  } else { // Anti_Quark on the end
1404  FirstString = new G4ExcitedString( end , Gquark, +1 );
1405  SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1406  Gquark->Set4Momentum( PkinkQ1 );
1407  Ganti_quark->Set4Momentum( PkinkQ2 );
1408  }
1409  } else { // Target
1410  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1411  FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1412  SecondString = new G4ExcitedString( start , Gquark, -1 );
1413  Ganti_quark->Set4Momentum( PkinkQ2 );
1414  Gquark->Set4Momentum( PkinkQ1 );
1415  } else { // Anti_Quark on the end
1416  FirstString = new G4ExcitedString( Gquark, end , -1 );
1417  SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1418  Gquark->Set4Momentum( PkinkQ2 );
1419  Ganti_quark->Set4Momentum( PkinkQ1 );
1420  }
1421  }
1422  } else { // Baryon/AntiBaryon
1423  //G4cout<<"isProjectile "<<isProjectile<<G4endl;
1424  //G4cout<<"end "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1425  //G4cout<<PkinkQ1<<G4endl;
1426  //G4cout<<PkinkQ2<<G4endl;
1427  //G4cout<<"start "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1428  //G4int Uzhi; G4cin>>Uzhi;
1429  if ( isProjectile ) { // Projectile
1430  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1431  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1432  FirstString = new G4ExcitedString( end , Gquark, +1 ); // Open Uzhi
1433  SecondString = new G4ExcitedString( Ganti_quark, start , +1 ); // Open Uzhi
1434  Gquark->Set4Momentum( PkinkQ1 );
1435  Ganti_quark->Set4Momentum( PkinkQ2 );
1436  //FirstString = new G4ExcitedString( Gquark, end, +1 );
1437  //SecondString = new G4ExcitedString( start, Ganti_quark, +1 );
1438  } else { // Anti_DiQuark on the end or quark
1439  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1440  SecondString = new G4ExcitedString( Gquark, start , +1 );
1441  Ganti_quark->Set4Momentum( PkinkQ1 );
1442  Gquark->Set4Momentum( PkinkQ2 );
1443  }
1444  } else { // Target
1445  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1446  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1447  //FirstString = new G4ExcitedString( Gquark, end , -1 );
1448  //SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1449  Gquark->Set4Momentum( PkinkQ1 );
1450  Ganti_quark->Set4Momentum( PkinkQ2 );
1451 
1452  FirstString = new G4ExcitedString( end, Gquark, -1 );
1453  SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1454 
1455  } else { // Anti_DiQuark on the end or Q
1456  FirstString = new G4ExcitedString( Ganti_quark, end , -1 ); // Uzhi ?????
1457  SecondString = new G4ExcitedString( start , Gquark, -1 );
1458  Gquark->Set4Momentum( PkinkQ2 );
1459  Ganti_quark->Set4Momentum( PkinkQ1 );
1460  }
1461  }
1462  }
1463 
1464  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1465  FirstString->SetPosition( hadron->GetPosition() );
1466  SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1467  SecondString->SetPosition( hadron->GetPosition() );
1468 
1469  } else { // End of kink is possible: Kink is impossible
1470 
1471  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1472  // << G4endl;
1473  /*
1474  if ( isProjectile ) {
1475  FirstString = new G4ExcitedString( end, start, +1 );
1476  } else {
1477  FirstString = new G4ExcitedString( start, end, -1 );
1478  }
1479  */
1480 
1481  FirstString = new G4ExcitedString( end, start, +1 );
1482 
1483  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1484  FirstString->SetPosition( hadron->GetPosition() );
1485  SecondString = 0;
1486 
1487  // momenta of string ends
1488  G4double Momentum = hadron->Get4Momentum().vect().mag();
1489  G4double Plus = hadron->Get4Momentum().e() + Momentum;
1490  G4double Minus = hadron->Get4Momentum().e() - Momentum;
1491  G4ThreeVector tmp;
1492  if ( Momentum > 0.0 ) {
1493  tmp.set( hadron->Get4Momentum().px(),
1494  hadron->Get4Momentum().py(),
1495  hadron->Get4Momentum().pz() );
1496  tmp /= Momentum;
1497  } else {
1498  tmp.set( 0.0, 0.0, 1.0 );
1499  }
1500  G4LorentzVector Pstart1( tmp, 0.0 );
1501  G4LorentzVector Pend1( tmp, 0.0 );
1502  if ( isProjectile ) {
1503  Pstart1 *= (-1.0)*Minus/2.0;
1504  Pend1 *= (+1.0)*Plus /2.0;
1505  } else {
1506  Pstart1 *= (+1.0)*Plus/ 2.0;
1507  Pend1 *= (-1.0)*Minus/2.0;
1508  }
1509  Momentum = -Pstart1.mag();
1510  Pstart1.setT( Momentum ); // It is assumed that quark has m=0.
1511  Momentum = -Pend1.mag();
1512  Pend1.setT( Momentum ); // It is assumed that di-quark has m=0.
1513  start->Set4Momentum( Pstart1 );
1514  end->Set4Momentum( Pend1 );
1515  SecondString = 0;
1516 
1517  } // End of kink is impossible
1518 
1519  //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1520  // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1521  // << FirstString << " " << SecondString << G4endl;
1522 
1523  #ifdef G4_FTFDEBUG
1524  G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1525  << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1526  << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1527  << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1528  << end->Get4Momentum().mag() << G4endl << " sum of ends "
1529  << Pstart + Pend << G4endl << " Original "
1530  << hadron->Get4Momentum() << G4endl;
1531  #endif
1532 
1533  return;
1534 }
1535 
1536 
1537 //============================================================================
1538 
1539 G4double G4DiffractiveExcitation::ChooseP( G4double Pmin, G4double Pmax ) const {
1540  // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1541  // To be improved...
1542  G4double range = Pmax - Pmin;
1543  if ( Pmin <= 0.0 || range <= 0.0 ) {
1544  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1545  throw G4HadronicException( __FILE__, __LINE__,
1546  "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1547  }
1548  G4double P = Pmin * G4Pow::GetInstance()->powA( Pmax/Pmin, G4UniformRand() );
1549  //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1550  return P;
1551 }
1552 
1553 
1554 //============================================================================
1555 
1556 G4ThreeVector G4DiffractiveExcitation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1557  // @@ this method is used in FTFModel as well. Should go somewhere common!
1558  G4double Pt2( 0.0 );
1559  if ( AveragePt2 <= 0.0 ) {
1560  Pt2 = 0.0;
1561  } else {
1562  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1563  ( G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1564  }
1565  G4double Pt = std::sqrt( Pt2 );
1566  G4double phi = G4UniformRand() * twopi;
1567  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1568 }
1569 
1570 
1571 //============================================================================
1572 
1573 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink( G4double zmin, G4double zmax ) const {
1574  G4double z, yf;
1575  const G4int maxNumberOfLoops = 10000;
1576  G4int loopCounter = 0;
1577  do {
1578  z = zmin + G4UniformRand() * (zmax - zmin);
1579  yf = z*z + sqr(1.0 - z);
1580  } while ( ( G4UniformRand() > yf ) &&
1581  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1582  if ( loopCounter >= maxNumberOfLoops ) {
1583  z = 0.5*(zmin + zmax); // Just something acceptable, without any physics consideration.
1584  }
1585  return z;
1586 }
1587 
1588 
1589 //============================================================================
1590 
1591 void G4DiffractiveExcitation::UnpackMeson( const G4int IdPDG, G4int& Q1, G4int& Q2 ) const {
1592  G4int absIdPDG = std::abs( IdPDG );
1593  if ( ! ( absIdPDG == 111 || absIdPDG == 221 || absIdPDG == 331 ) ) {
1594  // Ordinary mesons
1595  Q1 = absIdPDG / 100;
1596  Q2 = (absIdPDG % 100) / 10;
1597  G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1598  if ( IdPDG < 0 ) anti *= -1;
1599  Q1 *= anti;
1600  Q2 *= -1 * anti;
1601  } else {
1602  // Pi0, Eta, Eta'
1603  if ( G4UniformRand() < 0.5 ) { Q1 = 1; Q2 = -1; }
1604  else { Q1 = 2; Q2 = -2; }
1605  }
1606  return;
1607 }
1608 
1609 
1610 //============================================================================
1611 
1612 void G4DiffractiveExcitation::UnpackBaryon( G4int IdPDG,
1613  G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1614  Q1 = IdPDG / 1000;
1615  Q2 = (IdPDG % 1000) / 100;
1616  Q3 = (IdPDG % 100) / 10;
1617  return;
1618 }
1619 
1620 
1621 //============================================================================
1622 
1623 G4int G4DiffractiveExcitation::NewNucleonId( G4int Q1, G4int Q2, G4int Q3 ) const {
1624  G4int TmpQ( 0 );
1625  if ( Q3 > Q2 ) {
1626  TmpQ = Q2;
1627  Q2 = Q3;
1628  Q3 = TmpQ;
1629  } else if ( Q3 > Q1 ) {
1630  TmpQ = Q1;
1631  Q1 = Q3;
1632  Q3 = TmpQ;
1633  }
1634  if ( Q2 > Q1 ) {
1635  TmpQ = Q1;
1636  Q1 = Q2;
1637  Q2 = TmpQ;
1638  }
1639  G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1640  return NewCode;
1641 }
1642 
1643 
1644 //============================================================================
1645 
1646 G4DiffractiveExcitation::G4DiffractiveExcitation( const G4DiffractiveExcitation& ) {
1647  throw G4HadronicException( __FILE__, __LINE__,
1648  "G4DiffractiveExcitation copy contructor not meant to be called" );
1649 }
1650 
1651 
1652 //============================================================================
1653 
1654 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=( const G4DiffractiveExcitation& ) {
1655  throw G4HadronicException( __FILE__, __LINE__,
1656  "G4DiffractiveExcitation = operator not meant to be called" );
1657  return *this;
1658 }
1659 
1660 
1661 //============================================================================
1662 
1663 int G4DiffractiveExcitation::operator==( const G4DiffractiveExcitation& ) const {
1664  throw G4HadronicException( __FILE__, __LINE__,
1665  "G4DiffractiveExcitation == operator not meant to be called" );
1666 }
1667 
1668 
1669 //============================================================================
1670 
1671 int G4DiffractiveExcitation::operator!= ( const G4DiffractiveExcitation& ) const {
1672  throw G4HadronicException( __FILE__, __LINE__,
1673  "G4DiffractiveExcitation != operator not meant to be called" );
1674 }
1675 
G4Exp.hh
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4SampleResonance::SampleMass
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
Definition: G4SampleResonance.cc:110
G4ElasticHNScattering::ElasticScattering
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
Definition: G4ElasticHNScattering.cc:64
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4ExcitedString.hh
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
G4FTFParameters.hh
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:566
G4Parton::GetPDGcode
G4int GetPDGcode() const
Definition: G4Parton.hh:127
Y
double Y(double density)
Definition: G4DNAElectronHoleRecombination.cc:55
target
const XML_Char * target
Definition: expat.h:268
G4FTFParameters::GetProjMinDiffMass
G4double GetProjMinDiffMass()
Definition: G4FTFParameters.hh:429
G4ElasticHNScattering.hh
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:77
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4DiffractiveExcitation::ExciteParticipants
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
Definition: G4DiffractiveExcitation.cc:94
G4FTFParameters::GetProbLogDistr
G4double GetProbLogDistr()
Definition: G4FTFParameters.hh:453
G4Parton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
G4FTFParameters::GetTarMinDiffMass
G4double GetTarMinDiffMass()
Definition: G4FTFParameters.hh:437
G4long
long G4long
Definition: G4Types.hh:80
G4Parton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
G4Parton::GetDefinition
G4ParticleDefinition * GetDefinition()
Definition: G4Parton.hh:161
G4FTFParameters::GetPt2Kink
G4double GetPt2Kink()
Definition: G4FTFParameters.hh:459
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:142
G4SampleResonance.hh
G4VSplitableHadron::SetDefinition
void SetDefinition(const G4ParticleDefinition *aDefinition)
Definition: G4VSplitableHadron.hh:133
CLHEP::HepLorentzVector::rotateZ
HepLorentzVector & rotateZ(double)
Definition: LorentzVector.cc:80
G4int
int G4int
Definition: G4Types.hh:78
G4VSplitableHadron::SetStatus
void SetStatus(const G4int aStatus)
Definition: G4VSplitableHadron.hh:168
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:162
P
static double P[]
Definition: nf_gammaFunctions.cc:89
G4VSplitableHadron::SplitUp
virtual void SplitUp()=0
Definition: G4VSplitableHadron.cc:113
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4VSplitableHadron::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4VSplitableHadron.hh:138
G4ParticleDefinition::GetPDGWidth
G4double GetPDGWidth() const
Definition: G4ParticleDefinition.hh:123
CLHEP::HepLorentzVector::rapidity
double rapidity() const
Definition: LorentzVectorK.cc:128
G4ExcitedString::SetPosition
void SetPosition(const G4ThreeVector &aPosition)
Definition: G4ExcitedString.hh:155
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4VSplitableHadron.hh
G4ExcitedString::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition: G4ExcitedString.hh:143
pyG4EmCalculator::range
const G4ParticleDefinition const G4Material *G4double range
Definition: pyG4EmCalculator.cc:185
G4ParticleTable.hh
G4Neutron.hh
CLHEP::HepLorentzVector::mag
double mag() const
CLHEP::HepLorentzVector::py
double py() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VSplitableHadron::GetNextParton
virtual G4Parton * GetNextParton()=0
CLHEP::HepLorentzVector::px
double px() const
CLHEP::HepLorentzRotation::transform
HepLorentzRotation & transform(const HepBoost &b)
G4VSplitableHadron::IncrementCollisionCount
void IncrementCollisionCount(G4int aCount)
Definition: G4VSplitableHadron.hh:143
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:141
G4ParticleDefinition.hh
Randomize.hh
G4DiffractiveExcitation.hh
G4FTFParameters::GetProcProb
G4double GetProcProb(const G4int ProcN, const G4double y)
Definition: G4FTFParameters.cc:758
G4Log.hh
G4FTFParameters::GetAveragePt2
G4double GetAveragePt2()
Definition: G4FTFParameters.hh:445
G4PhysicalConstants.hh
G4VSplitableHadron::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4VSplitableHadron.hh:128
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:174
G4SampleResonance::GetMinimumMass
G4double GetMinimumMass(const G4ParticleDefinition *p) const
Definition: G4SampleResonance.cc:46
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4DiffractiveExcitation::CreateStrings
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
Definition: G4DiffractiveExcitation.cc:1193
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4FTFParameters::GetTarMinNonDiffMass
G4double GetTarMinNonDiffMass()
Definition: G4FTFParameters.hh:441
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4RotationMatrix.hh
CLHEP::HepLorentzVector::mag2
double mag2() const
CLHEP::HepLorentzVector::rotateY
HepLorentzVector & rotateY(double)
Definition: LorentzVector.cc:76
test.z
tuple z
Definition: test.py:28
CLHEP::HepLorentzVector::pz
double pz() const
G4VSplitableHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4VSplitableHadron.hh:163
G4FTFParameters::GetDeltaProbAtQuarkExchange
G4double GetDeltaProbAtQuarkExchange()
Definition: G4FTFParameters.hh:421
G4FTFParameters::GetProbOfSameQuarkExchange
G4double GetProbOfSameQuarkExchange()
Definition: G4FTFParameters.hh:425
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4FTFParameters::GetProjMinNonDiffMass
G4double GetProjMinNonDiffMass()
Definition: G4FTFParameters.hh:433
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4VSplitableHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4VSplitableHadron.hh:123
CLHEP::HepLorentzVector::getZ
double getZ() const
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4ThreeVector.hh
G4Pow.hh
CLHEP::Hep3Vector::mag
double mag() const
G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
Definition: G4FTFParameters.hh:463
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4LorentzRotation.hh