Geant4  10.00.p02
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 74627 2013-10-17 07:04:38Z 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, Uzhi
46 // ---------------------------------------------------------------------
47 
48 #include "globals.hh"
49 #include "Randomize.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 
54 #include "G4FTFParameters.hh"
55 #include "G4ElasticHNScattering.hh"
56 
57 #include "G4LorentzRotation.hh"
58 #include "G4RotationMatrix.hh"
59 #include "G4ThreeVector.hh"
60 #include "G4ParticleDefinition.hh"
61 #include "G4VSplitableHadron.hh"
62 #include "G4ExcitedString.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4Neutron.hh"
65 #include "G4ParticleDefinition.hh"
66 
67 //#include "G4ios.hh"
68 //#include "UZHI_diffraction.hh"
69 
70 
71 //============================================================================
72 
73 //#define debugFTFexictation
74 
75 
76 //============================================================================
77 
79 
80 
81 //============================================================================
82 
84 
85 
86 //============================================================================
87 
89  G4VSplitableHadron* target,
90  G4FTFParameters* theParameters,
91  G4ElasticHNScattering* theElastic ) const {
92 
93  #ifdef debugFTFexictation
94  G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
95  #endif
96 
97  // Projectile parameters
98  G4LorentzVector Pprojectile = projectile->Get4Momentum();
99  if ( Pprojectile.z() < 0.0 ) return false;
100 
101  G4double ProjectileRapidity = Pprojectile.rapidity();
102  G4int ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
103  G4int absProjectilePDGcode = std::abs( ProjectilePDGcode );
104 
105  G4bool PutOnMassShell( false );
106  //G4double M0projectile = projectile->GetDefinition()->GetPDGMass(); // With de-excitation
107  G4double M0projectile = Pprojectile.mag(); // Without de-excitation
108  if ( M0projectile < projectile->GetDefinition()->GetPDGMass() ) {
109  PutOnMassShell = true;
110  M0projectile = projectile->GetDefinition()->GetPDGMass();
111  }
112  G4double M0projectile2 = M0projectile * M0projectile;
113  G4double ProjectileDiffStateMinMass( 0.0 ), ProjectileNonDiffStateMinMass( 0.0 );
114  if ( M0projectile > projectile->GetDefinition()->GetPDGMass() ) { // Uzhi 31.05.2013
115  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV;
116  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV;
117  } else {
118  ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass();
119  ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass();
120  }
121 
122  // Target parameters
123  G4int TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
124  G4int absTargetPDGcode = std::abs( TargetPDGcode );
125  G4LorentzVector Ptarget = target->Get4Momentum();
126  G4double M0target = Ptarget.mag();
127  //G4double TargetRapidity = Ptarget.rapidity();
128 
129  #ifdef debugFTFexictation
130  G4cout << "Proj Targ PDGcodes " << ProjectilePDGcode << " " << TargetPDGcode << G4endl
131  << "M0projectile Y " << M0projectile << " " << ProjectileRapidity << G4endl;
132  //G4cout << "M0target Y " << M0target << " " << TargetRapidity << G4endl;
133  G4cout << "Pproj " << Pprojectile << G4endl << "Ptarget " << Ptarget << G4endl;
134  #endif
135 
136  if ( M0target < target->GetDefinition()->GetPDGMass() ) {
137  PutOnMassShell = true;
138  M0target=target->GetDefinition()->GetPDGMass();
139  }
140  G4double M0target2 = M0target * M0target;
141  G4double TargetDiffStateMinMass( 0.0 ), TargetNonDiffStateMinMass( 0.0 ); // Uzhi 31.05.2013
142  if ( M0target > target->GetDefinition()->GetPDGMass() ) {
143  TargetDiffStateMinMass = M0target + 220.0*MeV;
144  TargetNonDiffStateMinMass = M0target + 220.0*MeV;
145  } else {
146  TargetDiffStateMinMass = theParameters->GetTarMinDiffMass();
147  TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass();
148  }
149 
150  G4double AveragePt2 = theParameters->GetAveragePt2();
151  G4double ProbLogDistr = theParameters->GetProbLogDistr(); // 21.05.2012
152  G4double SumMasses = M0projectile + M0target + 220.0*MeV; // 200->220 7 June 2011
153 
154  // Kinematical properties of the interactions
155  G4LorentzVector Psum = Pprojectile + Ptarget; // 4-momentum in CMS
156  G4double S = Psum.mag2();
157  //Uzhi_SqrtS = std::sqrt( S );
158 
159  // Transform momenta to cms and then rotate parallel to z axis;
160  G4LorentzRotation toCms( -1 * Psum.boostVector() );
161 
162  G4LorentzVector Ptmp = toCms * Pprojectile;
163  if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
164 
165  toCms.rotateZ( -1*Ptmp.phi() );
166  toCms.rotateY( -1*Ptmp.theta() );
167  G4LorentzRotation toLab(toCms.inverse());
168  Pprojectile.transform( toCms );
169  Ptarget.transform( toCms );
170 
171  G4double PZcms2, PZcms;
172  G4double SqrtS = std::sqrt( S );
173 
174  #ifdef debugFTFexictation
175  G4cout << "SqrtS " << SqrtS << G4endl << "M0pr M0tr SumM+220 " << M0projectile << " "
176  << M0target << " " << SumMasses << G4endl;
177  #endif
178 
179  if ( SqrtS < M0projectile + M0target ) return false;
180  if ( SqrtS < SumMasses ) return false;
181  // The model cannot work at low energy
182 
183  PZcms2 = ( S*S + M0projectile2*M0projectile2 + M0target2*M0target2
184  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
185 
186  #ifdef debugFTFexictation
187  G4cout << "PZcms2 after PutOnMassShell " << PZcms2 << G4endl;
188  #endif
189 
190  if ( PZcms2 < 0 ) return false;
191  // It can be in an interaction with off-shell nuclear nucleon
192 
193  PZcms = std::sqrt( PZcms2 );
194  if ( PutOnMassShell ) {
195  if ( Pprojectile.z() > 0.0 ) {
196  Pprojectile.setPz( PZcms );
197  Ptarget.setPz( -PZcms );
198  } else {
199  Pprojectile.setPz( -PZcms );
200  Ptarget.setPz( PZcms );
201  };
202  Pprojectile.setE( std::sqrt( M0projectile2 +
203  Pprojectile.x()*Pprojectile.x() +
204  Pprojectile.y()*Pprojectile.y() +
205  PZcms2 ) );
206  Ptarget.setE( std::sqrt( M0target2 +
207  Ptarget.x()*Ptarget.x() +
208  Ptarget.y()*Ptarget.y() +
209  PZcms2 ) );
210  }
211 
212  G4double maxPtSquare; // = PZcms2;
213 
214  //Uzhi_QEnex = 0;
215  //Uzhi_QEexc = 0;
216  //Uzhi_targetdiffraction = 0;
217  //Uzhi_projectilediffraction = 0;
218  //Uzhi_nondiffraction = 0;
219  //G4int UzhiPrD( 0 ), UzhiTrD( 0 ), UzhiND( 0 );
220 
221  #ifdef debugFTFexictation
222  G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << M0projectile
223  << " " << ProjectileDiffStateMinMass << " " << ProjectileNonDiffStateMinMass << G4endl
224  << "Targ M0 Mdif Mndif " << M0target << " " << TargetDiffStateMinMass << " "
225  << TargetNonDiffStateMinMass << G4endl << "SqrtS " << SqrtS << G4endl
226  << "Proj CMS " << Pprojectile << G4endl << "Targ CMS " << Ptarget << G4endl;
227  #endif
228 
229  // Charge exchange can be possible
230  // Getting the values needed for exchange
231  // Check for possible quark exchange
232  G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity );
233  G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity );
234  G4double ProbProjectileDiffraction = theParameters->GetProcProb( 2, ProjectileRapidity );
235  G4double ProbTargetDiffraction = theParameters->GetProcProb( 3, ProjectileRapidity );
236  G4double ProbExc( 0.0 );
237  if ( QeExc + QeNoExc != 0.0 ) ProbExc = QeExc/(QeExc + QeNoExc);
238  G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
240 
241  #ifdef debugFTFexictation
242  G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " " << ProbProjectileDiffraction
243  << " " << ProbTargetDiffraction << G4endl
244  << "ProjectileRapidity " << ProjectileRapidity << G4endl;
245  //G4int Uzhi; G4cin >> Uzhi;
246  #endif
247 
248  if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
249  ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
250  ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
251  }
252 
253  if ( G4UniformRand() < QeExc + QeNoExc ) {
254 
255  #ifdef debugFTFexictation
256  G4cout << "Q exchange --------------------------" << G4endl;
257  #endif
258 
259  G4int NewProjCode( 0 ), NewTargCode( 0 );
260  G4int ProjQ1( 0 ), ProjQ2( 0 ), ProjQ3( 0 );
261  // Projectile unpacking
262  if ( absProjectilePDGcode < 1000 ) { // projectile is meson
263  UnpackMeson( ProjectilePDGcode, ProjQ1, ProjQ2 );
264  } else { // projectile is baryon
265  UnpackBaryon( ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
266  }
267  // Target unpacking
268  G4int TargQ1( 0 ), TargQ2( 0 ), TargQ3( 0 );
269  UnpackBaryon( TargetPDGcode, TargQ1, TargQ2, TargQ3 );
270 
271  #ifdef debugFTFexictation
272  G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
273  << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
274  #endif
275 
276  // Sampling of exchanged quarks
277  G4int ProjExchangeQ( 0 );
278  G4int TargExchangeQ( 0 );
279 
280  if ( absProjectilePDGcode < 1000 ) { // projectile is meson
281 
282  if ( ProjQ1 > 0 ) { // ProjQ1 is quark
283  ProjExchangeQ = ProjQ1;
284  G4int Nsampled = G4RandFlat::shootInt( G4long( 3 ) ) + 1;
285  if ( Nsampled == 1 ) {
286  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ1 = TargExchangeQ;
287  } else if ( Nsampled == 2 ) {
288  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ1 = TargExchangeQ;
289  } else {
290  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ1 = TargExchangeQ;
291  }
292  } else { // ProjQ2 is quark
293  ProjExchangeQ = ProjQ2;
294  G4int Nsampled = G4RandFlat::shootInt( G4long( 3 ) ) + 1;
295  if ( Nsampled == 1 ) {
296  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ2 = TargExchangeQ;
297  } else if ( Nsampled == 2 ) {
298  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ2 = TargExchangeQ;
299  } else {
300  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ2 = TargExchangeQ;
301  }
302  }
303 
304  #ifdef debugFTFexictation
305  G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
306  #endif
307 
308  G4int aProjQ1 = std::abs( ProjQ1 );
309  G4int aProjQ2 = std::abs( ProjQ2 );
310  if ( aProjQ1 == aProjQ2 ) {
311  NewProjCode = 111; // Pi0-meson
312  } else {
313  if ( aProjQ1 > aProjQ2 ) {
314  NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
315  } else {
316  NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
317  }
318  }
319 
320  #ifdef debugFTFexictation
321  G4cout << "NewProjCode " << NewProjCode << G4endl;
322  #endif
323 
324  G4bool ProjExcited = false;
325  if ( G4UniformRand() < 0.5 ) {
326  NewProjCode += 2; // Excited Pi0-meson
327  ProjExcited = true;
328  }
329  if ( aProjQ1 != aProjQ2 ) NewProjCode *= ( ProjectilePDGcode / absProjectilePDGcode );
330 
331  #ifdef debugFTFexictation
332  G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
333  #endif
334 
335  G4ParticleDefinition* TestParticle =
337 
338  if ( TestParticle ) {
339  G4double MtestPart =
341 
342  #ifdef debugFTFexictation
343  G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()
344  << G4endl << "MtestPart M0projectile projectile->GetDefinition()->GetPDGMass() "
345  << MtestPart << " " << M0projectile << " "
346  << projectile->GetDefinition()->GetPDGMass() << G4endl;
347  #endif
348 
349  if ( MtestPart > M0projectile ) {
350  M0projectile = MtestPart;
351  } else {
352  if ( std::abs( M0projectile - projectile->GetDefinition()->GetPDGMass() )
353  < 140.0*MeV ) {
354  M0projectile = MtestPart;
355  }
356  }
357 
358  #ifdef debugFTFexictation
359  G4cout << "M0projectile After check " << M0projectile << G4endl;
360  #endif
361 
362  M0projectile2 = M0projectile * M0projectile;
363  ProjectileDiffStateMinMass = M0projectile + 210.0*MeV; //210 MeV=m_pi+70 MeV
364  ProjectileNonDiffStateMinMass = M0projectile + 210.0*MeV; //210 MeV=m_pi+70 MeV
365  } else {
366  return false;
367  }
368 
369  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
370 
371  #ifdef debugFTFexictation
372  G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
373  << "NewTargCode " << NewTargCode << G4endl;
374  #endif
375 
376  // Lambda or Sigma0 ???
377  // if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) {
378  // if ( G4UniformRand() < 0.5 ) NewTargCode=
379  // if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 && SqrtS > M0projectile + DeltaMass ) {
380  // //Create Delta isobar
381  // NewTargCode +=2;
382 
383  if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
384  NewTargCode += 2; ProjExcited = true; //Create Delta isobar
385  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
386  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
387  NewTargCode += 2; ProjExcited = true; // Save Delta isobar
388  } else {
389  } // De-excite initial Delta isobar
390  } else if ( ! ProjExcited &&
391  G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
392  SqrtS > M0projectile + DeltaMass ) { // Create Delta isobar
393  NewTargCode +=2; // Save initial nucleon
394  } else {
395  }
396 
397  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
398 
399  #ifdef debugFTFexictation
400  G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
401  #endif
402 
403  if ( TestParticle ) {
404  G4double MtestPart =
406  if ( MtestPart > M0target ) {
407  M0target = MtestPart; M0target2 = sqr( M0target );
408  } else {
409  if ( std::abs( M0target - target->GetDefinition()->GetPDGMass() ) < 140.0*MeV ) {
410  M0target = MtestPart;
411  }
412  }
413 
414  #ifdef debugFTFexictation
415  G4cout << "New targ M0 M0^2 " << M0target << " " << M0target2 << G4endl;
416  #endif
417 
418  TargetDiffStateMinMass = M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
419  TargetNonDiffStateMinMass = M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
420 
421  if ( ProjExcited ) {
422  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0; // Uzhi May 2013
423  } else {
424  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0; // Uzhi May 2013
425  }
426  } else {
427  return false;
428  }
429 
430  } else { // of the if ( absProjectilePDGcode < 1000 ) ; i.e. the projectile is baryon
431 
432  G4double Same = theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0.
433  G4bool ProjDeltaHasCreated( false );
434  G4bool TargDeltaHasCreated( false );
435 
436  G4double Ksi = G4UniformRand();
437  if ( G4UniformRand() < 0.5 ) { // Sampling exchange quark from proj. or targ.
438  // Sampling exchanged quark from the projectile
439 
440  if ( Ksi < 0.333333 ) {
441  ProjExchangeQ = ProjQ1;
442  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
443  ProjExchangeQ = ProjQ2;
444  } else {
445  ProjExchangeQ = ProjQ3;
446  }
447 
448  if ( ProjExchangeQ != TargQ1 || G4UniformRand() < Same ) {
449  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
450  } else {
451  if ( ProjExchangeQ != TargQ2 || G4UniformRand() < Same ) {
452  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
453  } else {
454  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
455  }
456  }
457 
458  #ifdef debugFTFexictation
459  G4cout << "Exchange Qs Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
460  #endif
461 
462  if ( Ksi < 0.333333 ) {
463  ProjQ1 = ProjExchangeQ;
464  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
465  ProjQ2 = ProjExchangeQ;
466  } else {
467  ProjQ3 = ProjExchangeQ;
468  }
469 
470  } else { // Sampling exchanged quark from the target
471 
472  if ( Ksi < 0.333333 ) {
473  TargExchangeQ = TargQ1;
474  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
475  TargExchangeQ = TargQ2;
476  } else {
477  TargExchangeQ = TargQ3;
478  }
479  if ( TargExchangeQ != ProjQ1 || G4UniformRand() < Same ) {
480  ProjExchangeQ = ProjQ1; ProjQ1 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
481  } else {
482  if ( TargExchangeQ != ProjQ2 || G4UniformRand() < Same ) {
483  ProjExchangeQ = ProjQ2; ProjQ2 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
484  } else {
485  ProjExchangeQ = ProjQ3; ProjQ3 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
486  }
487  }
488 
489  if ( Ksi < 0.333333 ) {
490  TargQ1 = TargExchangeQ;
491  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
492  TargQ2 = TargExchangeQ;
493  } else {
494  TargQ3 = TargExchangeQ;
495  }
496 
497  } // End of sampling baryon
498 
499  NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
500 
501  if ( ProjQ1 == ProjQ2 && ProjQ1 == ProjQ3 ) {
502  NewProjCode += 2; ProjDeltaHasCreated = true;
503  } else if ( projectile->GetDefinition()->GetPDGiIsospin() == 3 ) { // Projectile was Delta
504  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
505  NewProjCode += 2; ProjDeltaHasCreated = true;
506  } else {
507  NewProjCode += 0; ProjDeltaHasCreated = false;
508  }
509  } else { // Projectile was Nucleon
510  if ( G4UniformRand() < DeltaProbAtQuarkExchange && SqrtS > DeltaMass + M0target ) {
511  NewProjCode += 2; ProjDeltaHasCreated = true;
512  } else {
513  NewProjCode += 0; ProjDeltaHasCreated = false;
514  }
515  }
516 
517  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
518 
519  if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
520  NewTargCode += 2; TargDeltaHasCreated = true;
521  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Target was Delta
522  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
523  NewTargCode += 2; TargDeltaHasCreated = true;
524  } else {
525  NewTargCode += 0; TargDeltaHasCreated = false;
526  }
527  } else { // Target was Nucleon
528  if ( G4UniformRand() < DeltaProbAtQuarkExchange && SqrtS > M0projectile + DeltaMass ) {
529  NewTargCode += 2; TargDeltaHasCreated = true;
530  } else {
531  NewTargCode += 0; TargDeltaHasCreated = false;
532  }
533  }
534 
535  #ifdef debugFTFexictation
536  G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
537  G4int Uzhi; G4cin >> Uzhi;
538  #endif
539 
540  if ( absProjectilePDGcode == NewProjCode && absTargetPDGcode == NewTargCode ) {
541  } // Nothing was changed! It is not right!?
542 
543  // Forming baryons
544 
545  if ( ProjDeltaHasCreated ) { // 8 Apr. 2013
546  if ( G4UniformRand() > 0.5 ) {
547  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
548  } else {
549  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
550  }
551  }
552 
553  if ( TargDeltaHasCreated ) {
554  if ( G4UniformRand() > 0.5 ) {
555  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
556  } else {
557  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
558  }
559  }
560 
561  if ( ProjDeltaHasCreated ) {
562  G4double MtestPart = // 31.05.2012
564  if ( MtestPart >= M0projectile ) { // 31.05.2012
565  M0projectile = MtestPart; // 31.05.2012
566  M0projectile2 = M0projectile * M0projectile; // 31.05.2012
567  } // 31.05.2012
568  ProjectileDiffStateMinMass = M0projectile + 210.0*MeV; //210 MeV=m_pi+70 MeV
569  ProjectileNonDiffStateMinMass = M0projectile + 210.0*MeV; //210 MeV=m_pi+70 MeV
570  }
571 
572  if ( TargDeltaHasCreated ) {
573  G4double MtestPart = // 31.05.2012
575  if ( MtestPart >= M0target ) { // 31.05.2012
576  M0target = MtestPart; // 31.05.2012
577  M0target2 = M0target * M0target; // 31.05.2012
578  } // 31.05.2012
579  TargetDiffStateMinMass = M0target + 210.0*MeV; //210 MeV=m_pi+70 MeV;
580  TargetNonDiffStateMinMass = M0target + 210.0*MeV; //210 MeV=m_pi+70 MeV;
581  }
582 
583  } // End of if ( absProjectilePDGcode < 1000 )
584 
585  // If we assume that final state hadrons after the charge exchange will be
586  // in the ground states, we have to put
587  if ( SqrtS < M0projectile + M0target ) return false;
588 
589  PZcms2 = ( S*S + M0projectile2*M0projectile2 + M0target2*M0target2
590  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
591 
592  #ifdef debugFTFexictation
593  G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
594  << "At the end// NewTargCode " << NewTargCode << G4endl
595  << "M0pr M0tr SqS " << M0projectile << " " << M0target << " " << SqrtS << G4endl
596  << "M0pr2 M0tr2 SqS " << M0projectile2 << " " << M0target2 << " " << SqrtS << G4endl
597  << "PZcms2 after the change " << PZcms2 << G4endl << G4endl;
598  #endif
599 
600  if ( PZcms2 < 0 ) return false; // It can be if energy is not sufficient for Delta
601 
602  projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
603  target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
604 
605  PZcms = std::sqrt( PZcms2 );
606  Pprojectile.setPz( PZcms );
607  Pprojectile.setE( std::sqrt( M0projectile2 + PZcms2 ) );
608  Ptarget.setPz( -PZcms );
609  Ptarget.setE( std::sqrt( M0target2 + PZcms2 ) );
610 
611  #ifdef debugFTFexictation
612  G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << Pprojectile << G4endl << Ptarget
613  << G4endl << Pprojectile + Ptarget << G4endl;
614  #endif
615 
616  if ( absProjectilePDGcode < 1000 ) { // For projectile meson
617  if ( G4UniformRand() > ProbExc ) { // Make elastic scattering
618 
619  #ifdef debugFTFexictation
620  G4cout << "Make elastic scattering of new hadrons" << G4endl;
621  #endif
622 
623  Pprojectile.transform( toLab );
624  Ptarget.transform( toLab );
625  // Uzhi 9.11 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
626  // Uzhi 9.11 projectile->SetPosition( target->GetPosition() );
627  projectile->Set4Momentum( Pprojectile );
628  target->Set4Momentum( Ptarget );
629 
630  G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
631 
632  #ifdef debugFTFexictation
633  G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
634  << G4endl << Pprojectile << G4endl << Ptarget << G4endl
635  << Pprojectile + Ptarget << " " << (Pprojectile + Ptarget).mag() << G4endl;
636  #endif
637 
638  //Uzhi_QEnex++;
639  return Result;
640  }
641  //Uzhi_QEexc++;
642  } else { // For projectile baryon
643  if ( G4UniformRand() > ProbExc ) { // Make elastic scattering
644 
645  #ifdef debugFTFexictation
646  G4cout << "Make elastic scattering of new hadrons" << G4endl;
647  #endif
648 
649  Pprojectile.transform( toLab );
650  Ptarget.transform( toLab );
651  // Uzhi 9.11 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
652  // Uzhi 9.11 projectile->SetPosition( target->GetPosition() );
653  projectile->Set4Momentum( Pprojectile );
654  target->Set4Momentum( Ptarget );
655 
656  G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
657 
658  //Uzhi_QEnex++;
659  return Result;
660  }
661  //Uzhi_QEexc++;
662  }
663 
664  #ifdef debugFTFexictation
665  G4cout << "Make excitation of new hadrons" << G4endl;
666  #endif
667 
668  } // End of if ( G4UniformRand() < QeExc + QeNoExc ) , i.e. of the charge exchange part
669 
670  G4double ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction;
671 
672  #ifdef debugFTFexictation
673  G4cout << "Excitation --------------------" << G4endl
674  << "Proj M0 MdMin MndMin " << M0projectile << " " << ProjectileDiffStateMinMass << " "
675  << ProjectileNonDiffStateMinMass << G4endl
676  << "Targ M0 MdMin MndMin " << M0target << " " << TargetDiffStateMinMass << " "
677  << TargetNonDiffStateMinMass << G4endl << "SqrtS " << SqrtS << G4endl
678  << "Prob: ProjDiff TargDiff + Sum " << ProbProjectileDiffraction << " "
679  << ProbTargetDiffraction << " " << ProbOfDiffraction << G4endl;
680  #endif
681 
682  if ( ProbOfDiffraction != 0.0 ) {
683  ProbProjectileDiffraction /= ProbOfDiffraction;
684  } else {
685  ProbProjectileDiffraction = 0.0;
686  }
687 
688  #ifdef debugFTFexictation
689  G4cout << "Prob: ProjDiff TargDiff + Sum " << ProbProjectileDiffraction << " "
690  << ProbTargetDiffraction << " " << ProbOfDiffraction << G4endl;
691  #endif
692 
693  G4double ProjectileDiffStateMinMass2 = sqr( ProjectileDiffStateMinMass );
694  G4double ProjectileNonDiffStateMinMass2 = sqr( ProjectileNonDiffStateMinMass );
695  G4double TargetDiffStateMinMass2 = sqr( TargetDiffStateMinMass );
696  G4double TargetNonDiffStateMinMass2 = sqr( TargetNonDiffStateMinMass );
697 
698  G4double Pt2;
699  G4double ProjMassT2, ProjMassT;
700  G4double TargMassT2, TargMassT;
701  G4double PMinusMin, PMinusMax;
702  //G4double PPlusMin , PPlusMax;
703  G4double TPlusMin, TPlusMax;
704  G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
705  G4LorentzVector Qmomentum;
706  G4double Qminus, Qplus;
707  G4int whilecount = 0;
708  //ProbOfDiffraction = 1.0; // Uzhi Difr For testing purposes
709  //ProbProjectileDiffraction = 1.0;
710 
711  // Choose a process
712  if ( G4UniformRand() < ProbOfDiffraction ) {
713 
714  if ( G4UniformRand() < ProbProjectileDiffraction ) { // projectile diffraction
715 
716  #ifdef debugFTFexictation
717  G4cout << "projectile diffraction" << G4endl;
718  #endif
719 
720  //UzhiPrD++;
721 
722  do { // while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileDiffStateMinMass2 )
723 
724  //Uzhi_projectilediffraction = 1;
725  //Uzhi_targetdiffraction = 0;
726  //Uzhi_Mx2 = 1.0;
727 
728  // Generate pt and mass of projectile
729 
730  whilecount++;
731  if ( whilecount > 1000 ) {
732  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
733  return false; // Ignore this interaction
734  };
735 
736  // Check that the interaction is possible
737  ProjMassT2 = ProjectileDiffStateMinMass2;
738  ProjMassT = ProjectileDiffStateMinMass;
739  TargMassT2 = M0target2;
740  TargMassT = M0target;
741  if ( SqrtS < ProjMassT + TargMassT ) return false;
742 
743  PZcms2 =( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
744  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
745 
746  if ( PZcms2 < 0 ) return false;
747 
748  maxPtSquare = PZcms2;
749 
750  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
751 
752  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
753  ProjMassT2 = ProjectileDiffStateMinMass2 + Pt2;
754  ProjMassT = std::sqrt( ProjMassT2 );
755  TargMassT2 = M0target2 + Pt2;
756  TargMassT = std::sqrt( TargMassT2 );
757  if ( SqrtS < ProjMassT + TargMassT ) continue;
758 
759  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
760  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
761 
762  if ( PZcms2 < 0 ) continue;
763 
764  PZcms = std::sqrt( PZcms2 );
765  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
766  PMinusMax = SqrtS - TargMassT;
767 
768  PMinusNew = ChooseP( PMinusMin, PMinusMax );
769 
770  // An attempt to use different distributions
771  //PMinusNew = 1.0/sqrt( 1.0/PMinusMin - G4UniformRand()*( 1.0/PMinusMin - 1.0/PMinusMax ) );
772  //PMinusNew = 1.0/sqr( 1.0/std::sqrt( PMinusMin ) -
773  // G4UniformRand()*( 1.0/std::sqrt( PMinusMin ) -
774  // 1.0/std::sqrt( PMinusMax ) ) );
775 
776  TMinusNew = SqrtS - PMinusNew;
777  Qminus = Ptarget.minus() - TMinusNew;
778  TPlusNew = TargMassT2 / TMinusNew;
779  Qplus = Ptarget.plus() - TPlusNew;
780  Qmomentum.setPz( (Qplus - Qminus)/2 );
781  Qmomentum.setE( (Qplus + Qminus)/2 );
782 
783  } while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileDiffStateMinMass2 );
784  // Repeat the sampling because there was not any excitation
785 
786  projectile->SetStatus( 1*projectile->GetStatus() ); // VU 10.04.2012
787 
788  } else { // Target diffraction
789 
790  #ifdef debugFTFexictation
791  G4cout << "Target diffraction" << G4endl;
792  #endif
793 
794  //UzhiTrD++;
795 
796  do { // while ( ( Ptarget - Qmomentum ).mag2() < TargetDiffStateMinMass2 )
797 
798  //Uzhi_projectilediffraction = 0;
799  //Uzhi_targetdiffraction = 1;
800  //Uzhi_Mx2 = 1.0;
801 
802  // Generate pt and target mass
803 
804  whilecount++;
805  if ( whilecount > 1000 ) {
806  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
807  return false; // Ignore this interaction
808  };
809 
810  // Check that the interaction is possible
811  ProjMassT2 = M0projectile2;
812  ProjMassT = M0projectile;
813  TargMassT2 = TargetDiffStateMinMass2;
814  TargMassT = TargetDiffStateMinMass;
815  if ( SqrtS < ProjMassT + TargMassT ) return false;
816 
817  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
818  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
819 
820  if ( PZcms2 < 0 ) return false;
821 
822  maxPtSquare = PZcms2;
823 
824  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
825 
826  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
827  ProjMassT2 = M0projectile2 + Pt2;
828  ProjMassT = std::sqrt( ProjMassT2 );
829  TargMassT2 = TargetDiffStateMinMass2 + Pt2;
830  TargMassT = std::sqrt( TargMassT2 );
831  if ( SqrtS < ProjMassT + TargMassT ) continue;
832 
833  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
834  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
835 
836  if ( PZcms2 < 0 ) continue;
837 
838  PZcms = std::sqrt( PZcms2 );
839  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
840  TPlusMax = SqrtS - ProjMassT;
841 
842  TPlusNew = ChooseP( TPlusMin, TPlusMax );
843 
844  //TPlusNew = 1.0/sqr( 1.0/std::sqrt( TPlusMin ) -
845  // G4UniformRand()*( 1.0/std::sqrt( TPlusMin ) -
846  // 1.0/std::sqrt( TPlusMax ) ) ); // 1/M^3/2
847  //TPlusNew = TPlusMax; // For extreme case of maximum mass
848 
849  PPlusNew = SqrtS - TPlusNew;
850  Qplus = PPlusNew - Pprojectile.plus();
851  PMinusNew = ProjMassT2 / PPlusNew;
852  Qminus = PMinusNew - Pprojectile.minus();
853  Qmomentum.setPz( (Qplus - Qminus)/2 );
854  Qmomentum.setE( (Qplus + Qminus)/2 );
855 
856  } while ( ( Ptarget - Qmomentum ).mag2() < TargetDiffStateMinMass2 );
857  // Repeat the sampling because there was not any excitation
858 
859  target->SetStatus( 1*target->GetStatus() ); // VU 10.04.2012
860 
861  } // End of if ( G4UniformRand() < ProbProjectileDiffraction )
862 
863  } else { // Non-diffraction process
864 
865  #ifdef debugFTFexictation
866  G4cout << "Non-diffraction process" << G4endl;
867  #endif
868 
869  //UzhiND++;
870 
871  do { // while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileNonDiffStateMinMass2 || ...
872 
873  //Uzhi_projectilediffraction = 0;
874  //Uzhi_targetdiffraction = 0;
875  //Uzhi_Mx2 = 1.0;
876 
877  // Generate pt and masses
878 
879  whilecount++;
880  if ( whilecount > 1000 ) {
881  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
882  return false; // Ignore this interaction
883  };
884 
885  // Check that the interaction is possible
886  ProjMassT2 = ProjectileNonDiffStateMinMass2;
887  ProjMassT = ProjectileNonDiffStateMinMass;
888  TargMassT2 = TargetNonDiffStateMinMass2;
889  TargMassT = TargetNonDiffStateMinMass;
890  if ( SqrtS < ProjMassT + TargMassT ) return false;
891 
892  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
893  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
894 
895  if ( PZcms2 < 0 ) return false;
896 
897  maxPtSquare = PZcms2;
898 
899  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
900 
901  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
902  ProjMassT2 = ProjectileNonDiffStateMinMass2 + Pt2;
903  ProjMassT = std::sqrt( ProjMassT2 );
904  TargMassT2 = TargetNonDiffStateMinMass2 + Pt2;
905  TargMassT = std::sqrt( TargMassT2 );
906  if ( SqrtS < ProjMassT + TargMassT ) continue;
907 
908  PZcms2 =( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
909  -2.0*S*ProjMassT2 - 2.0*S*TargMassT2 -2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
910 
911  if ( PZcms2 < 0 ) continue;
912 
913  PZcms = std::sqrt( PZcms2 );
914  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
915  PMinusMax = SqrtS - TargMassT;
916  if ( G4UniformRand() < ProbLogDistr ) { // Uzhi 25.04.2012
917  PMinusNew = ChooseP( PMinusMin, PMinusMax );
918  } else {
919  PMinusNew = ( PMinusMax - PMinusMin )*G4UniformRand() + PMinusMin;
920  }
921  Qminus = PMinusNew - Pprojectile.minus();
922  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
923  TPlusMax = SqrtS - PMinusNew;
924  //TPlusMax = SqrtS - ProjMassT; // To study an extreme case
925 
926  if ( G4UniformRand() < 0.5 ) { //ProbLogDistr) // Uzhi 29.05.2012 0.5)
927  TPlusNew = ChooseP( TPlusMin, TPlusMax );
928  } else {
929  TPlusNew = ( TPlusMax - TPlusMin )*G4UniformRand() + TPlusMin;
930  }
931  Qplus = -( TPlusNew - Ptarget.plus() );
932  Qmomentum.setPz( (Qplus - Qminus)/2 );
933  Qmomentum.setE( (Qplus + Qminus)/2 );
934 
935  #ifdef debugFTFexictation
936  G4cout << ( Pprojectile + Qmomentum ).mag2() << " " << ProjectileNonDiffStateMinMass2
937  << G4endl << ( Ptarget - Qmomentum ).mag2() << " "
938  << TargetNonDiffStateMinMass2 << G4endl;
939  G4int Uzhi; G4cin >> Uzhi;
940  #endif
941 
942  } while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileNonDiffStateMinMass2 || //No double Diffraction
943  ( Ptarget - Qmomentum ).mag2() < TargetNonDiffStateMinMass2 );
944 
945  projectile->SetStatus( 0*projectile->GetStatus() ); // VU 10.04.2012
946  target->SetStatus( 0*target->GetStatus() ); // VU 10.04.2012
947 
948  } // End of if ( G4UniformRand() < ProbOfDiffraction )
949 
950  Pprojectile += Qmomentum;
951  Ptarget -= Qmomentum;
952 
953  // Transform back and update SplitableHadron Participant.
954  Pprojectile.transform( toLab );
955  Ptarget.transform( toLab );
956 
957  // Calculation of the creation time
958  //Uzhi 9.11 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
959  //Uzhi 9.11 projectile->SetPosition( target->GetPosition() );
960  // Creation time and position of target nucleon were determined in
961  // ReggeonCascade() of G4FTFModel
962  //
963  //if ( Uzhi_projectilediffraction != 0 ) {
964  // Uzhi_Mx2 = Pprojectile.mag2(); Uzhi_modT = ( target->Get4Momentum() - Ptarget ).mag2();
965  //}
966  //if ( Uzhi_targetdiffraction != 0 ) {
967  // Uzhi_Mx2 = Ptarget.mag2(); Uzhi_modT = ( projectile->Get4Momentum() - Pprojectile ).mag2();
968  //}
969  //if ( Uzhi_QE != 0 ) {
970  // Uzhi_projectilediffraction = 0;
971  // Uzhi_targetdiffraction = 0;
972  // Uzhi_Mx2 = 1.0;
973  //}
974 
975  #ifdef debugFTFexictation
976  G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl;
977  #endif
978 
979  projectile->Set4Momentum( Pprojectile );
980  target->Set4Momentum( Ptarget );
981  projectile->IncrementCollisionCount( 1 );
982  target->IncrementCollisionCount( 1 );
983 
984  //Uzhi_projectilediffraction = UzhiPrD;
985  //Uzhi_targetdiffraction = UzhiTrD;
986  //Uzhi_nondiffraction = UzhiND;
987  //G4cout << Uzhi_projectilediffraction << " " << Uzhi_targetdiffraction << " "
988  // << Uzhi_nondiffraction << G4endl;
989 
990  return true;
991 }
992 
993 
994 //============================================================================
995 
997  G4bool isProjectile,
998  G4ExcitedString*& FirstString,
999  G4ExcitedString*& SecondString,
1000  G4FTFParameters* theParameters ) const {
1001 
1002  //G4cout << "Create Strings SplitUp " << hadron << G4endl
1003  // << "Defin " << hadron->GetDefinition() << G4endl
1004  // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
1005 
1006  hadron->SplitUp();
1007 
1008  G4Parton* start = hadron->GetNextParton();
1009  if ( start == NULL ) {
1010  G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1011  FirstString = 0; SecondString = 0;
1012  return;
1013  }
1014 
1015  G4Parton* end = hadron->GetNextParton();
1016  if ( end == NULL ) {
1017  G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1018  FirstString = 0; SecondString = 0;
1019  return;
1020  }
1021 
1022  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1023  // << G4endl
1024  // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1025 
1026  G4LorentzVector Phadron = hadron->Get4Momentum();
1027  //G4cout << "String mom " << Phadron << G4endl;
1028  G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1029  G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1030  G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1031  G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1032  G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1033 
1034  G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1035  G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1036  //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1037 
1038  G4double Wmin( 0.0 );
1039  if ( isProjectile ) {
1040  Wmin = theParameters->GetProjMinDiffMass();
1041  } else {
1042  Wmin = theParameters->GetTarMinDiffMass();
1043  }
1044 
1045  G4double W = hadron->Get4Momentum().mag();
1046  //G4cout << "Wmin W " << Wmin << " " << W << G4endl;
1047  //G4int Uzhi; G4cin >> Uzhi;
1048  G4double W2 = W*W;
1049  G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); // x2( 0.0 )
1050  G4bool Kink = false;
1051 
1052  if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1053  end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1054  ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1055  end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1056  // Kinky strings are allowed only for qq-q strings;
1057  // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1058  // according to the analysis of Pbar P interactions
1059 
1060  if ( W > Wmin ) { // Kink is possible
1061  if ( hadron->GetStatus() == 0 ) { // VU 10.04.2012
1062  G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1063  Pt = std::sqrt( Pt2kink * ( std::pow( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );
1064  } else {
1065  Pt = 0.0;
1066  }
1067 
1068  if ( Pt > 500.0*MeV ) {
1069  G4double Ymax = std::log( W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1070  G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
1071  x1 = 1.0 - Pt/W * std::exp( Y );
1072  x3 = 1.0 - Pt/W * std::exp(-Y );
1073  //x2 = 2.0 - x1 - x3;
1074 
1075  G4double Mass_startQ = 650.0*MeV;
1076  if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV;
1077  if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV;
1078  if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV;
1079  G4double Mass_endQ = 650.0*MeV;
1080  if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV;
1081  if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV;
1082  if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV;
1083 
1084  G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1085  G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1086  G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
1087  if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1088  Kink = false;
1089  } else {
1090  G4double P_1 = std::sqrt( P2_1 );
1091  G4double P_2 = std::sqrt( P2_2 );
1092  G4double P_3 = std::sqrt( P2_3 );
1093  G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1094  G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1095  //Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated 11.12.09
1096 
1097  if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1098  Kink = false;
1099  } else {
1100  Kink = true;
1101  Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated 11.12.09
1102  Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1103  Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1104  Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1105  Pstart.setE( x3*W/2.0 );
1106  Pkink.setE( Pkink.vect().mag() );
1107  Pend.setE( x1*W/2.0 );
1108 
1109  G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1110  if ( Pkink.getZ() > 0.0 ) {
1111  if ( XkQ > 0.5 ) {
1112  PkinkQ1 = XkQ*Pkink;
1113  } else {
1114  PkinkQ1 = (1.0 - XkQ)*Pkink;
1115  }
1116  } else {
1117  if ( XkQ > 0.5 ) {
1118  PkinkQ1 = (1.0 - XkQ)*Pkink;
1119  } else {
1120  PkinkQ1 = XkQ*Pkink;
1121  }
1122  }
1123 
1124  PkinkQ2 = Pkink - PkinkQ1;
1125  // Minimizing Pt1^2+Pt3^2
1126  G4double Cos2Psi = ( sqr(x1) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1127  std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
1128  G4double Psi = std::acos( Cos2Psi );
1129 
1130  G4LorentzRotation Rotate;
1131  if ( isProjectile ) {
1132  Rotate.rotateY( Psi );
1133  } else {
1134  Rotate.rotateY( pi - Psi );
1135  }
1136  Rotate.rotateZ( twopi * G4UniformRand() );
1137  Pstart *= Rotate;
1138  Pkink *= Rotate;
1139  PkinkQ1 *= Rotate;
1140  PkinkQ2 *= Rotate;
1141  Pend *= Rotate;
1142  }
1143  } // End of if ( P2_1 <= 0.0 || P2_3 <= 0.0 )
1144  } // End of if ( Pt > 500.0*MeV )
1145  } // End of if ( W > Wmin ) : check for a kink
1146  } // end of qq-q string selection
1147 
1148  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetParticleSubType() << " "
1149  // << end->GetDefinition()->GetParticleSubType() << G4endl;
1150  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetPDGEncoding() << " "
1151  // << end->GetDefinition()->GetPDGEncoding() << G4endl;
1152  //G4int Uzhi; G4cin >> Uzhi;
1153 
1154  if ( Kink ) { // Kink is possible
1155 
1156  //G4cout << "Kink is sampled!" << G4endl;
1157  std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1158  theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1159 
1160  G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1161  for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1162  //G4cout << "Iq " << Iq << G4endl;
1163  if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1164  }
1165  //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1166  G4Parton* Gquark = new G4Parton( QuarkInGluon );
1167  G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1168  //G4cout << "Lorentz " << G4endl;
1169 
1170  G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1171  G4LorentzRotation toLab( toCMS.inverse() );
1172  //G4cout << "Pstart " << Pstart << G4endl;
1173  //G4cout << "Pend " << Pend << G4endl;
1174  Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1175  PkinkQ1.transform( toLab );
1176  PkinkQ2.transform( toLab );
1177  Pend.transform( toLab ); end->Set4Momentum( Pend );
1178  //G4cout << "Pstart " << Pstart << G4endl;
1179  //G4cout << "Pend " << Pend << G4endl;
1180  //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1181  //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1182 
1183  //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1184  G4int absPDGcode = 1500; // 23 Dec
1185  if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1186  end->GetDefinition()->GetParticleSubType() == "quark" ) {
1187  absPDGcode = 110;
1188  }
1189  //G4cout << "absPDGcode " << absPDGcode << G4endl;
1190 
1191  if ( absPDGcode < 1000 ) { // meson
1192  if ( isProjectile ) { // Projectile
1193  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1194  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1195  SecondString = new G4ExcitedString( Gquark, start , +1 );
1196  Ganti_quark->Set4Momentum( PkinkQ1 );
1197  Gquark->Set4Momentum( PkinkQ2 );
1198  } else { // Anti_Quark on the end
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  if ( isProjectile ) { // Projectile
1219  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1220  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1221  FirstString = new G4ExcitedString( end , Gquark, +1 );
1222  SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1223  Gquark->Set4Momentum( PkinkQ1 );
1224  Ganti_quark->Set4Momentum( PkinkQ2 );
1225  } else { // Anti_DiQuark on the end or quark
1226  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1227  SecondString = new G4ExcitedString( Gquark, start , +1 );
1228  Ganti_quark->Set4Momentum( PkinkQ1 );
1229  Gquark->Set4Momentum( PkinkQ2 );
1230  }
1231  } else { // Target
1232  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1233  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1234  FirstString = new G4ExcitedString( Gquark, end , -1 );
1235  SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1236  Gquark->Set4Momentum( PkinkQ1 );
1237  Ganti_quark->Set4Momentum( PkinkQ2 );
1238  } else { // Anti_DiQuark on the end or Q
1239  FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1240  SecondString = new G4ExcitedString( start , Gquark, -1 );
1241  Gquark->Set4Momentum( PkinkQ2 );
1242  Ganti_quark->Set4Momentum( PkinkQ1 );
1243  }
1244  }
1245  }
1246 
1247  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1248  FirstString->SetPosition( hadron->GetPosition() );
1249  SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1250  SecondString->SetPosition( hadron->GetPosition() );
1251 
1252  } else { // End of kink is possible: Kink is impossible
1253 
1254  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1255  // << G4endl;
1256  if ( isProjectile ) {
1257  FirstString = new G4ExcitedString( end, start, +1 );
1258  } else {
1259  FirstString = new G4ExcitedString( start, end, -1 );
1260  }
1261  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1262  FirstString->SetPosition( hadron->GetPosition() );
1263  SecondString = 0;
1264 
1265  // momenta of string ends
1266  G4double Momentum = hadron->Get4Momentum().vect().mag();
1267  G4double Plus = hadron->Get4Momentum().e() + Momentum;
1268  G4double Minus = hadron->Get4Momentum().e() - Momentum;
1269  G4ThreeVector tmp;
1270  if ( Momentum > 0.0 ) {
1271  tmp.set( hadron->Get4Momentum().px(),
1272  hadron->Get4Momentum().py(),
1273  hadron->Get4Momentum().pz() );
1274  tmp /= Momentum;
1275  } else {
1276  tmp.set( 0.0, 0.0, 1.0 );
1277  }
1278  G4LorentzVector Pstart1( tmp, 0.0 );
1279  G4LorentzVector Pend1( tmp, 0.0 );
1280  if ( isProjectile ) {
1281  Pstart1 *= (-1.0)*Minus/2.0;
1282  Pend1 *= (+1.0)*Plus /2.0;
1283  } else {
1284  Pstart1 *= (+1.0)*Plus/ 2.0;
1285  Pend1 *= (-1.0)*Minus/2.0;
1286  }
1287  Momentum = -Pstart1.mag();
1288  Pstart1.setT( Momentum ); // It is assumed that quark has m=0.
1289  Momentum = -Pend1.mag();
1290  Pend1.setT( Momentum ); // It is assumed that di-quark has m=0.
1291  start->Set4Momentum( Pstart1 );
1292  end->Set4Momentum( Pend1 );
1293  SecondString = 0;
1294 
1295  } // End of kink is impossible
1296 
1297  //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1298  // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1299  // << FirstString << " " << SecondString << G4endl;
1300 
1301  #ifdef G4_FTFDEBUG
1302  G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1303  << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1304  << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1305  << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1306  << end->Get4Momentum().mag() << G4endl << " sum of ends "
1307  << Pstart + Pend << G4endl << " Original "
1308  << hadron->Get4Momentum() << G4endl;
1309  #endif
1310 
1311  return;
1312 }
1313 
1314 
1315 //============================================================================
1316 
1318  // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1319  // To be improved...
1320  G4double range = Pmax - Pmin;
1321  if ( Pmin <= 0.0 || range <= 0.0 ) {
1322  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1323  throw G4HadronicException( __FILE__, __LINE__,
1324  "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1325  }
1326  G4double P = Pmin * std::pow( Pmax/Pmin, G4UniformRand() );
1327  //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1328  return P;
1329 }
1330 
1331 
1332 //============================================================================
1333 
1335  // @@ this method is used in FTFModel as well. Should go somewhere common!
1336  G4double Pt2( 0.0 );
1337  if ( AveragePt2 <= 0.0 ) {
1338  Pt2 = 0.0;
1339  } else {
1340  Pt2 = -AveragePt2 * std::log( 1.0 + G4UniformRand() *
1341  ( std::exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1342  }
1343  G4double Pt = std::sqrt( Pt2 );
1344  G4double phi = G4UniformRand() * twopi;
1345  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1346 }
1347 
1348 
1349 //============================================================================
1350 
1352  G4double z, yf;
1353  do {
1354  z = zmin + G4UniformRand() * (zmax - zmin);
1355  yf = z*z + sqr(1.0 - z);
1356  } while ( G4UniformRand() > yf );
1357  return z;
1358 }
1359 
1360 
1361 //============================================================================
1362 
1363 void G4DiffractiveExcitation::UnpackMeson( const G4int IdPDG, G4int& Q1, G4int& Q2 ) const {
1364  G4int absIdPDG = std::abs( IdPDG );
1365  Q1 = absIdPDG / 100;
1366  Q2 = (absIdPDG % 100) / 10;
1367  G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1368  if ( IdPDG < 0 ) anti *= -1;
1369  Q1 *= anti;
1370  Q2 *= -1 * anti;
1371  return;
1372 }
1373 
1374 
1375 //============================================================================
1376 
1378  G4int& Q1, G4int& Q2, G4int& Q3 ) const {
1379  Q1 = IdPDG / 1000;
1380  Q2 = (IdPDG % 1000) / 100;
1381  Q3 = (IdPDG % 100) / 10;
1382  return;
1383 }
1384 
1385 
1386 //============================================================================
1387 
1389  G4int TmpQ( 0 );
1390  if ( Q3 > Q2 ) {
1391  TmpQ = Q2;
1392  Q2 = Q3;
1393  Q3 = TmpQ;
1394  } else if ( Q3 > Q1 ) {
1395  TmpQ = Q1;
1396  Q1 = Q3;
1397  Q3 = TmpQ;
1398  }
1399  if ( Q2 > Q1 ) {
1400  TmpQ = Q1;
1401  Q1 = Q2;
1402  Q2 = TmpQ;
1403  }
1404  G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1405  return NewCode;
1406 }
1407 
1408 
1409 //============================================================================
1410 
1412  throw G4HadronicException( __FILE__, __LINE__,
1413  "G4DiffractiveExcitation copy contructor not meant to be called" );
1414 }
1415 
1416 
1417 //============================================================================
1418 
1420  throw G4HadronicException( __FILE__, __LINE__,
1421  "G4DiffractiveExcitation = operator not meant to be called" );
1422  return *this;
1423 }
1424 
1425 
1426 //============================================================================
1427 
1429  throw G4HadronicException( __FILE__, __LINE__,
1430  "G4DiffractiveExcitation == operator not meant to be called" );
1431 }
1432 
1433 
1434 //============================================================================
1435 
1437  throw G4HadronicException( __FILE__, __LINE__,
1438  "G4DiffractiveExcitation != operator not meant to be called" );
1439 }
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
static const double MeV
Definition: G4SIunits.hh:193
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4int GetPDGcode() const
Definition: G4Parton.hh:124
G4double GetProjMinDiffMass()
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
G4int NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const
int operator==(const G4DiffractiveExcitation &right) const
G4double z
Definition: TRTMaterials.hh:39
G4double GetProbLogDistr()
const G4double pi
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:140
G4double GetTarMinDiffMass()
long G4long
Definition: G4Types.hh:80
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:145
G4ParticleDefinition * GetDefinition()
Definition: G4Parton.hh:158
G4double GetPt2Kink()
const G4String & GetParticleSubType() const
int G4int
Definition: G4Types.hh:78
void SetStatus(const G4int aStatus)
void SetDefinition(G4ParticleDefinition *aDefinition)
const G4String & GetParticleName() const
virtual void SplitUp()=0
G4ParticleDefinition * GetDefinition() const
#define G4cin
Definition: G4ios.hh:60
void SetPosition(const G4ThreeVector &aPosition)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
void SetTimeOfCreation(G4double aTime)
bool G4bool
Definition: G4Types.hh:79
virtual G4Parton * GetNextParton()=0
int operator!=(const G4DiffractiveExcitation &right) const
void UnpackMeson(G4int IdPDG, G4int &Q1, G4int &Q2) const
G4double GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
void IncrementCollisionCount(G4int aCount)
const G4String & GetParticleType() const
G4double GetProcProb(const G4int ProcN, const G4double y)
G4double GetAveragePt2()
const G4LorentzVector & Get4Momentum() const
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
G4double ChooseP(G4double Pmin, G4double Pmax) const
G4double GetPDGMass() const
G4double GetTarMinNonDiffMass()
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4ThreeVector & GetPosition() const
G4double GetDeltaProbAtQuarkExchange()
const G4DiffractiveExcitation & operator=(const G4DiffractiveExcitation &right)
G4double GetProbOfSameQuarkExchange()
#define G4endl
Definition: G4ios.hh:61
G4double GetProjMinNonDiffMass()
void Set4Momentum(const G4LorentzVector &a4Momentum)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
void UnpackBaryon(G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
CLHEP::HepLorentzVector G4LorentzVector