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