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