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