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