Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FTFAnnihilation.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 // ---------------- G4FTFAnnihilation --------------
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 "G4DiffractiveSplitableHadron.hh"
52 #include "G4DiffractiveExcitation.hh"
53 #include "G4FTFParameters.hh"
54 #include "G4ElasticHNScattering.hh"
55 #include "G4FTFAnnihilation.hh"
56 
57 #include "G4LorentzRotation.hh"
58 #include "G4RotationMatrix.hh"
59 #include "G4ThreeVector.hh"
60 #include "G4ParticleDefinition.hh"
61 #include "G4VSplitableHadron.hh"
62 #include "G4ExcitedString.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4Neutron.hh"
65 #include "G4ParticleDefinition.hh"
66 
67 //#include "G4ios.hh"
68 //#include "UZHI_diffraction.hh"
69 
70 G4FTFAnnihilation::G4FTFAnnihilation()
71 {
72 }
73 
74 // ---------------------------------------------------------------------
75 G4bool G4FTFAnnihilation::
76  Annihilate(G4VSplitableHadron *projectile,
77  G4VSplitableHadron *target,
78  G4VSplitableHadron *&AdditionalString,
79  G4FTFParameters *theParameters) const
80 {
81 // -------------------- Projectile parameters -----------------------
82  G4LorentzVector Pprojectile=projectile->Get4Momentum();
83 //G4cout<<"---------------------------- Annihilation----------------"<<G4endl;
84 //G4cout<<"Pprojectile "<<Pprojectile<<G4endl;
85 //G4cout<<"Pprojectile.mag2 "<<Pprojectile.mag2()<<G4endl;
86  G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
87  if(ProjectilePDGcode > 0)
88  {
89  target->SetStatus(2);
90  return false;
91  }
92 
93 // G4double M0projectile = Pprojectile.mag();
94 // G4double M0projectile2= projectile->GetDefinition()->GetPDGMass()*
95 // projectile->GetDefinition()->GetPDGMass();
96  G4double M0projectile2=Pprojectile.mag2();
97 // -------------------- Target parameters -------------------------
98  G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
99 
100  G4LorentzVector Ptarget=target->Get4Momentum();
101 //G4cout<<"Ptarget "<<Ptarget<<G4endl;
102 //G4cout<<"Ptarget.mag2 "<<Ptarget.mag2()<<G4endl;
103 
104 // G4double M0target = Ptarget.mag();
105 // G4double M0target2= target->GetDefinition()->GetPDGMass()*
106 // target->GetDefinition()->GetPDGMass();
107  G4double M0target2=Ptarget.mag2();
108 
109 //G4cout<<"Annihilate "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
110 //G4cout<<"Pprojec "<<Pprojectile<<" "<<Pprojectile.mag2()<<G4endl;
111 //G4cout<<"Ptarget "<<Ptarget <<" "<<Ptarget.mag2() <<G4endl;
112 //G4cout<<"M0 proj target "<<M0projectile<<" "<<M0target<<G4endl;
113 
114  G4double AveragePt2=theParameters->GetAveragePt2();
115 
116 // Kinematical properties of the interactions --------------
117  G4LorentzVector Psum; // 4-momentum in CMS
118  Psum=Pprojectile+Ptarget;
119  G4double S=Psum.mag2();
120 //G4cout<<"Psum S"<<Psum<<" "<<S<<G4endl;
121 // Transform momenta to cms and then rotate parallel to z axis;
122  G4LorentzRotation toCms(-1*Psum.boostVector());
123 //G4cout<<"G4LorentzRotation toCms(-1*Psum.boostVector());"<<G4endl;
124  G4LorentzVector Ptmp=toCms*Pprojectile;
125 
126 /* // For anti-baryons it is not needed !
127  if ( Ptmp.pz() <= 0. )
128  {
129  target->SetStatus(2);
130  // "String" moving backwards in CMS, abort collision !!
131  return false;
132  }
133 */
134 
135  toCms.rotateZ(-1*Ptmp.phi());
136  toCms.rotateY(-1*Ptmp.theta());
137 
138  G4LorentzRotation toLab(toCms.inverse());
139 
140  G4double SqrtS=std::sqrt(S);
141 
142  G4double maxPtSquare;
143 
144 //G4cout<<"M0projectile+M0target Sqrt(S) (GeV) "<<M0projectile2/GeV<<" "<<M0target2/GeV<<" "<<(M0projectile2+M0target2)/GeV<<" "<<SqrtS/GeV<<G4endl;
145 
146  G4double X_a(0.), X_b(0.), X_c(0.), X_d(0.);
147  G4double MesonProdThreshold=projectile->GetDefinition()->GetPDGMass()+
148  target->GetDefinition()->GetPDGMass()+
149  (2.*140.+16.)*MeV; // 2 Mpi +DeltaE
150 
151  G4double Prel2= S*S + M0projectile2*M0projectile2 + M0target2*M0target2 -
152  2.*S*M0projectile2 - 2.*S*M0target2 - 2.*M0projectile2*M0target2;
153  Prel2/=S;
154 
155  if(Prel2 < 0. ) // *MeV*MeV 1600.
156  { // Annihilation at rest! Values are copied from Paratemets.
157  X_a= 625.1; // mb // 3-shirt diagram
158  X_b= 9.780; // mb // anti-quark-quark annihilation
159  X_c= 49.989; // mb
160  X_d= 6.614; // mb
161  }
162  else
163  { // Annihilation in flight!
164  G4double FlowF=1./std::sqrt(Prel2)*GeV;
165 
166 //G4cout<<"Annig FlowF "<<FlowF<<" sqrt "<<SqrtS/GeV<<G4endl;
167 
168 // Process cross sections ---------------------------------------------------
169  X_a=25.*FlowF; // mb 3-shirt diagram
170 
171  // mb anti-quark-quark annihilation
172  if(SqrtS < MesonProdThreshold)
173  {
174  X_b=3.13+140.*std::pow((MesonProdThreshold - SqrtS)/GeV,2.5);
175  }
176  else
177  {
178  X_b=6.8*GeV/SqrtS;
179  }
180  if(projectile->GetDefinition()->GetPDGMass()+
181  target->GetDefinition()->GetPDGMass() > SqrtS) {X_b=0.;}
182 // This can be in an interaction of low energy anti-baryon with off-shell nuclear nucleon
183 
184 // ????????????????????????????????????????
185  X_c=2.*FlowF*sqr(projectile->GetDefinition()->GetPDGMass()+
186  target->GetDefinition()->GetPDGMass())/S;
187  // mb re-arrangement of 2 quarks and 2 anti-quarks
188 // ????????????????????????????????????????
189  X_d=23.3*GeV*GeV/S; // mb anti-quark-quark string creation
190  } // end of if(Prel2 < 1600. ) // *MeV*MeV
191 
192 //G4cout<<"Annih X a b c d "<<X_a<<" "<<X_b<<" "<<X_c<<" "<<X_d<<G4endl;
193 
194  if((ProjectilePDGcode == -2212)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
195  {X_b*=5.; X_c*=5.; X_d*=6.;} // Pbar P
196  else if((ProjectilePDGcode == -2212)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
197  {X_b*=4.; X_c*=4.; X_d*=4.;} // Pbar N
198  else if((ProjectilePDGcode == -2112)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
199  {X_b*=4.; X_c*=4.; X_d*=4.;} // NeutrBar P
200  else if((ProjectilePDGcode == -2112)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
201  {X_b*=5.; X_c*=5.; X_d*=6.;} // NeutrBar N
202  else if((ProjectilePDGcode == -3122)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
203  {X_b*=3.; X_c*=3.; X_d*=2.;} // LambdaBar P
204  else if((ProjectilePDGcode == -3122)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
205  {X_b*=3.; X_c*=3.; X_d*=2.;} // LambdaBar N
206  else if((ProjectilePDGcode == -3112)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
207  {X_b*=2.; X_c*=2.; X_d*=0.;} // Sigma-Bar P
208  else if((ProjectilePDGcode == -3112)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
209  {X_b*=4.; X_c*=4.; X_d*=2.;} // Sigma-Bar N
210  else if((ProjectilePDGcode == -3212)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
211  {X_b*=3.; X_c*=3.; X_d*=2.;} // Sigma0Bar P
212  else if((ProjectilePDGcode == -3212)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
213  {X_b*=3.; X_c*=3.; X_d*=2.;} // Sigma0Bar N
214  else if((ProjectilePDGcode == -3222)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
215  {X_b*=4.; X_c*=4.; X_d*=2.;} // Sigma+Bar P
216  else if((ProjectilePDGcode == -3222)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
217  {X_b*=2.; X_c*=2.; X_d*=0.;} // Sigma+Bar N
218  else if((ProjectilePDGcode == -3312)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
219  {X_b*=1.; X_c*=1.; X_d*=0.;} // Xi-Bar P
220  else if((ProjectilePDGcode == -3312)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
221  {X_b*=2.; X_c*=2.; X_d*=0.;} // Xi-Bar N
222  else if((ProjectilePDGcode == -3322)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
223  {X_b*=2.; X_c*=2.; X_d*=0.;} // Xi0Bar P
224  else if((ProjectilePDGcode == -3322)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
225  {X_b*=1.; X_c*=1.; X_d*=0.;} // Xi0Bar N
226  else if((ProjectilePDGcode == -3334)&&((TargetPDGcode == 2212)||(TargetPDGcode == 2214)))
227  {X_b*=0.; X_c*=0.; X_d*=0.;} // Omega-Bar P
228  else if((ProjectilePDGcode == -3334)&&((TargetPDGcode == 2112)||(TargetPDGcode == 2114)))
229  {X_b*=0.; X_c*=0.; X_d*=0.;} // Omega-Bar N
230  else {G4cout<<"Unknown anti-baryon for FTF annihilation: PDGcodes - "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;}
231 
232 //G4cout<<"Annih X a b c d "<<X_a<<" "<<X_b<<" "<<X_c<<" "<<X_d<<G4endl;
233 //=========================================
234 //X_a=0.;
235 //X_b=0.;
236 //X_c=0.;
237 //X_d=0.;
238 //G4cout<<"Annih X a b c d "<<X_a<<" "<<X_b<<" "<<X_c<<" "<<X_d<<G4endl;
239 //=========================================
240  G4double Xannihilation=X_a+X_b+X_c+X_d;
241 // ------------------------------------------------------
242 // ------ Projectile unpacking --------------------------
243  G4int AQ[3];
244  UnpackBaryon(ProjectilePDGcode, AQ[0], AQ[1], AQ[2]);
245 
246 // ------ Target unpacking ------------------------------
247  G4int Q[3];
248  UnpackBaryon(TargetPDGcode, Q[0], Q[1], Q[2]);
249 
250 // ------------------------------------------------------
251  G4double Ksi=G4UniformRand();
252 
253  if(Ksi < X_a/Xannihilation)
254  {
255 //============================================================
256 // Simulation of 3 anti-quark-quark strings creation
257 // Sampling of anti-quark order in projectile
258 //G4cout<<"Process a"<<G4endl;
259  G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(6));
260 
261  G4int Tmp1(0), Tmp2(0);
262  if(SampledCase == 0) { }
263  if(SampledCase == 1) {Tmp1=AQ[1]; AQ[1]=AQ[2]; AQ[2]=Tmp1;}
264  if(SampledCase == 2) {Tmp1=AQ[0]; AQ[0]=AQ[1]; AQ[1]=Tmp1;}
265  if(SampledCase == 3) {Tmp1=AQ[0]; Tmp2=AQ[1]; AQ[0]=AQ[2]; AQ[1]=Tmp1; AQ[2]=Tmp2;}
266  if(SampledCase == 4) {Tmp1=AQ[0]; Tmp2=AQ[1]; AQ[0]=Tmp2; AQ[1]=AQ[2]; AQ[2]=Tmp1;}
267  if(SampledCase == 5) {Tmp1=AQ[0]; Tmp2=AQ[1]; AQ[0]=AQ[2]; AQ[1]=Tmp2; AQ[2]=Tmp1;}
268 
269 // --------------- Set the string properties ---------------
270 //G4cout<<"String 1 "<<AQ[0]<<" "<<Q[0]<<G4endl;
271  projectile->SplitUp();
272 
273  projectile->SetFirstParton(AQ[0]);
274  projectile->SetSecondParton(Q[0]);
275  projectile->SetStatus(1);
276 
277 //G4cout<<"String 2 "<<Q[1]<<" "<<AQ[1]<<G4endl;
278  target->SplitUp();
279 
280  target->SetFirstParton(Q[1]);
281  target->SetSecondParton(AQ[1]);
282  target->SetStatus(1);
283 
284 //G4cout<<"String 3 "<<AQ[2]<<" "<<Q[2]<<G4endl;
285  AdditionalString=new G4DiffractiveSplitableHadron();
286  AdditionalString->SplitUp();
287  AdditionalString->SetFirstParton(AQ[2]);
288  AdditionalString->SetSecondParton(Q[2]);
289  AdditionalString->SetStatus(1);
290 //G4cout<<G4endl<<"*AdditionalString in Annih"<<AdditionalString<<G4endl;
291 
292 // Sampling kinematical properties
293 // 1 string AQ[0]-Q[0]// 2 string AQ[1]-Q[1]// 3 string AQ[2]-Q[2]
294 
295  G4ThreeVector Quark_Mom[6];
296  G4double ModMom2[6]; //ModMom[6],
297 
298 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
299  AveragePt2=200.*200.; maxPtSquare=S;
300 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
301 
302  G4double SumMt(0.);
303 
304  G4double MassQ2=0.; //100.*100.*MeV*MeV;
305 
306  G4int NumberOfTries(0);
307  G4double ScaleFactor(1.);
308  do // while(SumMt >SqrtS)
309  {
310  NumberOfTries++;
311 
312  if(NumberOfTries == 100*(NumberOfTries/100))
313  { // At large number of tries it would be better to reduce the values of <Pt^2>
314  ScaleFactor/=2.;
315  AveragePt2 *=ScaleFactor;
316  }
317 
318  G4ThreeVector PtSum(0.,0.,0.);
319  for(G4int i=0; i<6; i++)
320  {
321  Quark_Mom[i]=GaussianPt(AveragePt2, maxPtSquare);
322  PtSum+=Quark_Mom[i];
323  }
324 
325  PtSum/=6.;
326 
327  SumMt=0.;
328  for(G4int i=0; i<6; i++)
329  {
330  Quark_Mom[i]-=PtSum;
331 // ModMom[i] =Quark_Mom[i].mag();
332  ModMom2[i]=Quark_Mom[i].mag2();
333  SumMt+=std::sqrt(ModMom2[i]+MassQ2);
334  }
335  } while(SumMt > SqrtS);
336 
337  G4double WminusTarget(0.), WplusProjectile(0.);
338 /* //--------------------- Closed is variant with sampling of Xs at minimum
339  G4double SumMod_anti=ModMom[0]+ModMom[1]+ModMom[2];
340  Quark_Mom[0].setZ(ModMom[0]/SumMod_anti);
341  Quark_Mom[1].setZ(ModMom[1]/SumMod_anti);
342  Quark_Mom[2].setZ(ModMom[2]/SumMod_anti);
343 
344  G4double SumMod_bary=ModMom[3]+ModMom[4]+ModMom[5];
345  Quark_Mom[3].setZ(ModMom[3]/SumMod_bary);
346  Quark_Mom[4].setZ(ModMom[4]/SumMod_bary);
347  Quark_Mom[5].setZ(ModMom[5]/SumMod_bary);
348 
349  G4double Alfa=SumMod_anti*SumMod_anti;
350  G4double Beta=SumMod_bary*SumMod_bary;
351 //------------------------------------
352  G4double DecayMomentum2=S*S + Alfa*Alfa + Beta*Beta
353  - 2.*S*Alfa - 2.*S*Beta - 2.*Alfa*Beta;
354 
355  WminusTarget=(S-Alfa+Beta+std::sqrt(DecayMomentum2))/2./SqrtS;
356  WplusProjectile=SqrtS-Beta/WminusTarget;
357 */ //--------------------- Closed is variant with sampling of Xs at minimum
358 //
359 // // ------------------------------------------------
360 // Sampling X's of anti-baryon -------
361  G4double Alfa_R=0.5;
362 
363  NumberOfTries=0;
364  ScaleFactor=1.;
365 
366  G4bool Succes(true);
367  do // while(!Succes)
368  {
369  Succes=true;
370  NumberOfTries++;
371 
372  if(NumberOfTries == 100*(NumberOfTries/100))
373  { // At large number of tries it would be better to reduce the values of Pt's
374  ScaleFactor/=2.;
375  }
376 
377  if(Alfa_R == 1.)
378  {
379  G4double Xaq1=1.-std::sqrt(G4UniformRand());
380  G4double Xaq2=(1.-Xaq1)*G4UniformRand();
381  G4double Xaq3=1.-Xaq1-Xaq2;
382 
383  Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2); Quark_Mom[2].setZ(Xaq3);
384  }
385  else
386  {
387  G4double Xaq1=sqr(G4UniformRand());
388  G4double Xaq2=(1.-Xaq1)*sqr(std::sin(pi/2.*G4UniformRand()));
389  G4double Xaq3=1.-Xaq1-Xaq2;
390 
391  Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2); Quark_Mom[2].setZ(Xaq3);
392  } // end of if(Alfa_R == 0.)
393 
394 // Sampling X's of baryon ------------
395  if(Alfa_R == 1.)
396  {
397  G4double Xq1=1.-std::sqrt(G4UniformRand());
398  G4double Xq2=(1.-Xq1)*G4UniformRand();
399  G4double Xq3=1.-Xq1-Xq2;
400 
401  Quark_Mom[3].setZ(Xq1); Quark_Mom[4].setZ(Xq2); Quark_Mom[5].setZ(Xq3);
402  }
403  else
404  {
405  G4double Xq1=sqr(G4UniformRand());
406  G4double Xq2=(1.-Xq1)*sqr(std::sin(pi/2.*G4UniformRand()));
407  G4double Xq3=1.-Xq1-Xq2;
408 
409  Quark_Mom[3].setZ(Xq1); Quark_Mom[4].setZ(Xq2); Quark_Mom[5].setZ(Xq3);
410  } // end of if(Alfa_R == 0.)
411 //
412  G4double Alfa(0.), Beta(0.);
413 
414  for(G4int i=0; i<3; i++) // For Anti-baryon
415  {
416  if(Quark_Mom[i].getZ() != 0.)
417  {Alfa+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
418  else {Succes=false;}
419  }
420 
421  for(G4int i=3; i<6; i++) // For baryon
422  {
423  if(Quark_Mom[i].getZ() != 0.)
424  {Beta+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
425  else {Succes=false;}
426  }
427 
428  if(!Succes) continue;
429 
430  if(std::sqrt(Alfa)+std::sqrt(Beta) > SqrtS) {Succes=false; continue;}
431 
432  G4double DecayMomentum2=S*S + Alfa*Alfa + Beta*Beta
433  - 2.*S*Alfa - 2.*S*Beta - 2.*Alfa*Beta;
434 
435  WminusTarget=(S-Alfa+Beta+std::sqrt(DecayMomentum2))/2./SqrtS;
436  WplusProjectile=SqrtS-Beta/WminusTarget;
437 
438  } while(!Succes);
439 // //--------------------------------------------------
440 
441  G4double SqrtScaleF=std::sqrt(ScaleFactor);
442 
443  for(G4int i=0; i<3; i++)
444  {
445  G4double Pz=WplusProjectile*Quark_Mom[i].getZ()/2.-
446  (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WplusProjectile*Quark_Mom[i].getZ());
447  Quark_Mom[i].setZ(Pz);
448 
449  if(ScaleFactor != 1.)
450  {
451  Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
452  Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
453  }
454  }
455 
456  for(G4int i=3; i<6; i++)
457  {
458  G4double Pz=-WminusTarget*Quark_Mom[i].getZ()/2.+
459  (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WminusTarget*Quark_Mom[i].getZ());
460  Quark_Mom[i].setZ(Pz);
461 
462  if(ScaleFactor != 1.)
463  {
464  Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
465  Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
466  }
467  }
468 //G4cout<<"Sum AQ "<<Quark_Mom[0]+Quark_Mom[1]+Quark_Mom[2]<<G4endl;
469 //G4cout<<"Sum Q "<<Quark_Mom[3]+Quark_Mom[4]+Quark_Mom[5]<<G4endl;
470 //-------------------------------------
471 
472  G4ThreeVector tmp=Quark_Mom[0]+Quark_Mom[3];
473  G4LorentzVector Pstring1(tmp,std::sqrt(Quark_Mom[0].mag2()+MassQ2)+
474  std::sqrt(Quark_Mom[3].mag2()+MassQ2));
475  G4double Ystring1=Pstring1.rapidity();
476 /*
477 G4cout<<"Mom 1 string "<<G4endl;
478 G4cout<<Quark_Mom[0]<<G4endl;
479 G4cout<<Quark_Mom[3]<<G4endl;
480 G4cout<<tmp<<" "<<tmp.mag()<<G4endl;
481 */
482 //G4cout<<"1 str "<<Pstring1<<" "<<Pstring1.mag()<<" "<<Ystring1<<G4endl;
483 
484  tmp=Quark_Mom[1]+Quark_Mom[4];
485  G4LorentzVector Pstring2(tmp,std::sqrt(Quark_Mom[1].mag2()+MassQ2)+
486  std::sqrt(Quark_Mom[4].mag2()+MassQ2));
487  G4double Ystring2=Pstring2.rapidity();
488 /*
489 G4cout<<"Mom 2 string "<<G4endl;
490 G4cout<<Quark_Mom[1]<<G4endl;
491 G4cout<<Quark_Mom[4]<<G4endl;
492 G4cout<<tmp<<" "<<tmp.mag()<<G4endl;
493 */
494 //G4cout<<"2 str "<<Pstring2<<" "<<Pstring2.mag()<<" "<<Ystring2<<G4endl;
495 
496  tmp=Quark_Mom[2]+Quark_Mom[5];
497  G4LorentzVector Pstring3(tmp,std::sqrt(Quark_Mom[2].mag2()+MassQ2)+
498  std::sqrt(Quark_Mom[5].mag2()+MassQ2));
499  G4double Ystring3=Pstring3.rapidity();
500 /*
501 G4cout<<"Mom 3 string "<<G4endl;
502 G4cout<<Quark_Mom[2]<<G4endl;
503 G4cout<<Quark_Mom[5]<<G4endl;
504 G4cout<<tmp<<" "<<tmp.mag()<<G4endl;
505 */
506 //G4cout<<"3 str "<<Pstring3<<" "<<Pstring3.mag()<<" "<<Ystring3<<G4endl;
507 //G4cout<<"SumE "<<Pstring1.e()+Pstring2.e()+Pstring3.e()<<G4endl;
508 //G4cout<<Pstring1.mag()<<" "<<Pstring2.mag()<<" "<<Pstring3.mag()<<G4endl;
509 //G4int Uzhi; G4cin>>Uzhi;
510 //--------------------------------
511  G4LorentzVector LeftString(0.,0.,0.,0.);
512 //-----
513  if((Ystring1 > Ystring2)&&(Ystring2 > Ystring3))
514  {
515  Pprojectile=Pstring1;
516  LeftString =Pstring2;
517  Ptarget =Pstring3;
518  }
519 
520  if((Ystring1 > Ystring3)&&(Ystring3 > Ystring2))
521  {
522  Pprojectile=Pstring1;
523  LeftString =Pstring3;
524  Ptarget =Pstring2;
525  }
526 //-----
527  if((Ystring2 > Ystring1)&&(Ystring1 > Ystring3))
528  {
529  Pprojectile=Pstring2;
530  LeftString =Pstring1;
531  Ptarget =Pstring3;
532  }
533 
534  if((Ystring2 > Ystring3)&&(Ystring3 > Ystring1))
535  {
536  Pprojectile=Pstring2;
537  LeftString =Pstring3;
538  Ptarget =Pstring1;
539  }
540 //-----
541  if((Ystring3 > Ystring1)&&(Ystring1 > Ystring2))
542  {
543  Pprojectile=Pstring3;
544  LeftString =Pstring1;
545  Ptarget =Pstring2;
546  }
547 
548  if((Ystring3 > Ystring2)&&(Ystring2 > Ystring1))
549  {
550  Pprojectile=Pstring3;
551  LeftString =Pstring2;
552  Ptarget =Pstring1;
553  }
554 
555 //-------------------------------------------------------
556 //G4cout<<"SumP "<<Pprojectile+LeftString+Ptarget<<" "<<SqrtS<<G4endl;
557 
558  Pprojectile.transform(toLab);
559  LeftString.transform(toLab);
560  Ptarget.transform(toLab);
561 //G4cout<<"SumP "<<Pprojectile+LeftString+Ptarget<<" "<<SqrtS<<G4endl;
562 
563 // Calculation of the creation time ---------------------
564  projectile->SetTimeOfCreation(target->GetTimeOfCreation());
565  projectile->SetPosition(target->GetPosition());
566 
567  AdditionalString->SetTimeOfCreation(target->GetTimeOfCreation());
568  AdditionalString->SetPosition(target->GetPosition());
569 // Creation time and position of target nucleon were determined at
570 // ReggeonCascade() of G4FTFModel
571 // ------------------------------------------------------
572 
573 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
574 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
575  projectile->Set4Momentum(Pprojectile);
576  AdditionalString->Set4Momentum(LeftString);
577  target->Set4Momentum(Ptarget);
578 
579  projectile->IncrementCollisionCount(1);
580  AdditionalString->IncrementCollisionCount(1);
581  target->IncrementCollisionCount(1);
582 
583  return true;
584  }
585 
586 //============================================================
587 // Simulation of anti-diquark-diquark string creation
588 //
589  if(Ksi < (X_a+X_b)/Xannihilation)
590  {
591 //G4cout<<"Process b"<<G4endl;
592  G4int CandidatsN(0), CandAQ[9][2], CandQ[9][2];
593  G4int LeftAQ1(0), LeftAQ2(0), LeftQ1(0), LeftQ2(0);
594 //------------------------------------------------------------
595  for(G4int iAQ=0; iAQ<3; iAQ++)
596  {
597  for(G4int iQ=0; iQ<3; iQ++)
598  {
599  if(-AQ[iAQ] == Q[iQ])
600  {
601  if(iAQ == 0) {CandAQ[CandidatsN][0]=1; CandAQ[CandidatsN][1]=2;}
602  if(iAQ == 1) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=2;}
603  if(iAQ == 2) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=1;}
604  if(iQ == 0) {CandQ[CandidatsN][0] =1; CandQ[CandidatsN][1]=2;}
605  if(iQ == 1) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=2;}
606  if(iQ == 2) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=1;}
607  CandidatsN++;
608  } //end of if(-AQ[i] == Q[j])
609  } //end of cycle on targ. quarks
610  } //end of cycle on proj. anti-quarks
611 //------------------------------------------------------------
612 //G4cout<<"CandidatsN "<<CandidatsN<<G4endl;
613 
614  if(CandidatsN != 0)
615  {
616  G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(CandidatsN));
617 
618  LeftAQ1=AQ[CandAQ[SampledCase][0]];
619  LeftAQ2=AQ[CandAQ[SampledCase][1]];
620 
621  LeftQ1=Q[CandQ[SampledCase][0]];
622  LeftQ2=Q[CandQ[SampledCase][1]];
623 
624 // -------- Build anti-diquark and diquark
625  G4int Anti_DQ(0), DQ(0);
626 
627  if(std::abs(LeftAQ1) > std::abs(LeftAQ2))
628  {
629  Anti_DQ=1000*LeftAQ1+100*LeftAQ2-3; // 1
630  } else
631  {
632  Anti_DQ=1000*LeftAQ2+100*LeftAQ1-3; // 1
633  }
634 // if(G4UniformRand() > 0.5) Anti_DQ-=2;
635 
636  if(std::abs(LeftQ1) > std::abs(LeftQ2))
637  {
638  DQ=1000*LeftQ1+100*LeftQ2+3; // 1
639  } else
640  {
641  DQ=1000*LeftQ2+100*LeftQ1+3; // 1
642  }
643 // if(G4UniformRand() > 0.5) DQ+=2;
644 
645 // --------------- Set the string properties ---------------
646 //G4cout<<"Left ADiQ DiQ "<<Anti_DQ<<" "<<DQ<<G4endl;
647 
648  projectile->SplitUp();
649 
650 // projectile->SetFirstParton(Anti_DQ);
651 // projectile->SetSecondParton(DQ);
652  projectile->SetFirstParton(DQ);
653  projectile->SetSecondParton(Anti_DQ);
654 
655  projectile->SetStatus(1);
656  target->SetStatus(3); // The target nucleon has annihilated
657 
658  Pprojectile.setPx(0.); // VU Mar1
659  Pprojectile.setPy(0.); // VU Mar1
660  Pprojectile.setPz(0.);
661  Pprojectile.setE(SqrtS);
662  Pprojectile.transform(toLab);
663 
664 // Calculation of the creation time ---------------------
665  projectile->SetTimeOfCreation(target->GetTimeOfCreation());
666  projectile->SetPosition(target->GetPosition());
667 // Creation time and position of target nucleon were determined at
668 // ReggeonCascade() of G4FTFModel
669 // ------------------------------------------------------
670 
671 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
672 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
673  projectile->Set4Momentum(Pprojectile);
674 
675  projectile->IncrementCollisionCount(1);
676 
677  return true;
678  } // end of if(CandidatsN != 0)
679  } // if(Ksi < (X_a+X_b)/Xannihilation)
680 
681 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
682 
683  if(Ksi < (X_a+X_b+X_c)/Xannihilation)
684  {
685 //============================================================
686 // Simulation of 2 anti-quark-quark strings creation
687 //G4cout<<"Process c"<<G4endl;
688  G4int CandidatsN(0), CandAQ[9][2], CandQ[9][2];
689  G4int LeftAQ1(0), LeftAQ2(0), LeftQ1(0), LeftQ2(0);
690 //------------------------------------------------------------
691  for(G4int iAQ=0; iAQ<3; iAQ++)
692  {
693  for(G4int iQ=0; iQ<3; iQ++)
694  {
695  if(-AQ[iAQ] == Q[iQ])
696  {
697  if(iAQ == 0) {CandAQ[CandidatsN][0]=1; CandAQ[CandidatsN][1]=2;}
698  if(iAQ == 1) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=2;}
699  if(iAQ == 2) {CandAQ[CandidatsN][0]=0; CandAQ[CandidatsN][1]=1;}
700  if(iQ == 0) {CandQ[CandidatsN][0] =1; CandQ[CandidatsN][1]=2;}
701  if(iQ == 1) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=2;}
702  if(iQ == 2) {CandQ[CandidatsN][0] =0; CandQ[CandidatsN][1]=1;}
703  CandidatsN++;
704  } //end of if(-AQ[i] == Q[j])
705  } //end of cycle on targ. quarks
706  } //end of cycle on proj. anti-quarks
707 //------------------------------------------------------------
708 //G4cout<<"CandidatsN "<<CandidatsN<<G4endl;
709 
710  if(CandidatsN != 0)
711  {
712  G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(CandidatsN));
713 
714  LeftAQ1=AQ[CandAQ[SampledCase][0]];
715  LeftAQ2=AQ[CandAQ[SampledCase][1]];
716 
717  if(G4UniformRand() < 0.5)
718  {
719  LeftQ1=Q[CandQ[SampledCase][0]];
720  LeftQ2=Q[CandQ[SampledCase][1]];
721  } else
722  {
723  LeftQ2=Q[CandQ[SampledCase][0]];
724  LeftQ1=Q[CandQ[SampledCase][1]];
725  }
726 
727 // --------------- Set the string properties ---------------
728 //G4cout<<"String 1 "<<LeftAQ1<<" "<<LeftQ1<<G4endl;
729  projectile->SplitUp();
730 
731  projectile->SetFirstParton(LeftAQ1);
732  projectile->SetSecondParton(LeftQ1);
733  projectile->SetStatus(1);
734 
735 //G4cout<<"String 2 "<<LeftAQ2<<" "<<LeftQ2<<G4endl;
736  target->SplitUp();
737 
738  target->SetFirstParton(LeftQ2);
739  target->SetSecondParton(LeftAQ2);
740  target->SetStatus(1);
741 
742 // Sampling kinematical properties
743 // 1 string LeftAQ1-LeftQ1// 2 string LeftAQ2-LeftQ2
744 
745  G4ThreeVector Quark_Mom[4];
746  G4double ModMom2[4]; //ModMom[4],
747 
748 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
749  AveragePt2=200.*200.; maxPtSquare=S;
750 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
751 
752  G4double SumMt(0.);
753 
754  G4double MassQ2=0.; //100.*100.*MeV*MeV;
755 
756  G4int NumberOfTries(0);
757  G4double ScaleFactor(1.);
758  do // while(SumMt >SqrtS)
759  {
760  NumberOfTries++;
761 
762  if(NumberOfTries == 100*(NumberOfTries/100))
763  { // At large number of tries it would be better to reduce the values of <Pt^2>
764  ScaleFactor/=2.;
765  AveragePt2 *=ScaleFactor;
766  }
767 
768  G4ThreeVector PtSum(0.,0.,0.);
769  for(G4int i=0; i<4; i++)
770  {
771  Quark_Mom[i]=GaussianPt(AveragePt2, maxPtSquare);
772  PtSum+=Quark_Mom[i];
773  }
774 
775  PtSum/=4.;
776 
777  SumMt=0.;
778  for(G4int i=0; i<4; i++)
779  {
780  Quark_Mom[i]-=PtSum;
781 // ModMom[i] =Quark_Mom[i].mag();
782  ModMom2[i]=Quark_Mom[i].mag2();
783  SumMt+=std::sqrt(ModMom2[i]+MassQ2);
784  }
785  } while(SumMt > SqrtS);
786 
787  G4double WminusTarget(0.), WplusProjectile(0.);
788 
789 // Sampling X's of anti-baryon -------
790  G4double Alfa_R=0.5;
791 
792  NumberOfTries=0;
793  ScaleFactor=1.;
794 
795  G4bool Succes(true);
796  do // while(!Succes)
797  {
798  Succes=true;
799  NumberOfTries++;
800 
801  if(NumberOfTries == 100*(NumberOfTries/100))
802  { // At large number of tries it would be better to reduce the values of Pt's
803  ScaleFactor/=2.;
804  }
805 
806  if(Alfa_R == 1.)
807  {
808  G4double Xaq1=std::sqrt(G4UniformRand());
809  G4double Xaq2=1.-Xaq1;
810 
811  Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2);
812  }
813  else
814  {
815  G4double Xaq1=sqr(std::sin(pi/2.*G4UniformRand()));
816  G4double Xaq2=1.-Xaq1;
817 
818  Quark_Mom[0].setZ(Xaq1); Quark_Mom[1].setZ(Xaq2);
819  } // end of if(Alfa_R == 0.)
820 
821 // Sampling X's of baryon ------------
822  if(Alfa_R == 1.)
823  {
824  G4double Xq1=1.-std::sqrt(G4UniformRand());
825  G4double Xq2=1.-Xq1;
826 
827  Quark_Mom[2].setZ(Xq1); Quark_Mom[3].setZ(Xq2);
828  }
829  else
830  {
831  G4double Xq1=sqr(std::sin(pi/2.*G4UniformRand()));
832  G4double Xq2=1.-Xq1;
833 
834  Quark_Mom[2].setZ(Xq1); Quark_Mom[3].setZ(Xq2);
835  } // end of if(Alfa_R == 0.)
836 //
837  G4double Alfa(0.), Beta(0.);
838 
839  for(G4int i=0; i<2; i++) // For Anti-baryon
840  {
841  if(Quark_Mom[i].getZ() != 0.)
842  {Alfa+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
843  else {Succes=false;}
844  }
845 
846  for(G4int i=2; i<4; i++) // For baryon
847  {
848  if(Quark_Mom[i].getZ() != 0.)
849  {Beta+=(ScaleFactor*ModMom2[i]+MassQ2)/Quark_Mom[i].getZ();}
850  else {Succes=false;}
851  }
852 
853  if(!Succes) continue;
854 
855  if(std::sqrt(Alfa)+std::sqrt(Beta) > SqrtS) {Succes=false; continue;}
856 
857  G4double DecayMomentum2=S*S + Alfa*Alfa + Beta*Beta
858  - 2.*S*Alfa - 2.*S*Beta - 2.*Alfa*Beta;
859 
860  WminusTarget=(S-Alfa+Beta+std::sqrt(DecayMomentum2))/2./SqrtS;
861  WplusProjectile=SqrtS-Beta/WminusTarget;
862 
863  } while(!Succes);
864 // //--------------------------------------------------
865 
866  G4double SqrtScaleF=std::sqrt(ScaleFactor);
867 
868  for(G4int i=0; i<2; i++)
869  {
870  G4double Pz=WplusProjectile*Quark_Mom[i].getZ()/2.-
871  (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WplusProjectile*Quark_Mom[i].getZ());
872  Quark_Mom[i].setZ(Pz);
873 
874  if(ScaleFactor != 1.)
875  {
876  Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
877  Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
878  }
879 //G4cout<<"Anti Q "<<i<<" "<<Quark_Mom[i]<<G4endl;
880  }
881 
882  for(G4int i=2; i<4; i++)
883  {
884  G4double Pz=-WminusTarget*Quark_Mom[i].getZ()/2.+
885  (ScaleFactor*ModMom2[i]+MassQ2)/(2.*WminusTarget*Quark_Mom[i].getZ());
886  Quark_Mom[i].setZ(Pz);
887 
888  if(ScaleFactor != 1.)
889  {
890  Quark_Mom[i].setX(SqrtScaleF*Quark_Mom[i].getX());
891  Quark_Mom[i].setY(SqrtScaleF*Quark_Mom[i].getY());
892  }
893 //G4cout<<"Bary Q "<<i<<" "<<Quark_Mom[i]<<G4endl;
894  }
895 //G4cout<<"Sum AQ "<<Quark_Mom[0]+Quark_Mom[1]<<G4endl;
896 //G4cout<<"Sum Q "<<Quark_Mom[2]+Quark_Mom[3]<<G4endl;
897 //-------------------------------------
898 
899  G4ThreeVector tmp=Quark_Mom[0]+Quark_Mom[2];
900  G4LorentzVector Pstring1(tmp,std::sqrt(Quark_Mom[0].mag2()+MassQ2)+
901  std::sqrt(Quark_Mom[2].mag2()+MassQ2));
902  G4double Ystring1=Pstring1.rapidity();
903 /*
904 G4cout<<"Mom 1 string "<<G4endl;
905 G4cout<<Quark_Mom[0]<<G4endl;
906 G4cout<<Quark_Mom[2]<<G4endl;
907 G4cout<<tmp<<" "<<tmp.mag()<<G4endl;
908 //G4cout<<"1 str "<<Pstring1<<" "<<Pstring1.mag()<<" "<<Ystring1<<G4endl;
909 */
910 
911  tmp=Quark_Mom[1]+Quark_Mom[3];
912  G4LorentzVector Pstring2(tmp,std::sqrt(Quark_Mom[1].mag2()+MassQ2)+
913  std::sqrt(Quark_Mom[3].mag2()+MassQ2));
914  G4double Ystring2=Pstring2.rapidity();
915 /*
916 G4cout<<"Mom 2 string "<<G4endl;
917 G4cout<<Quark_Mom[1]<<G4endl;
918 G4cout<<Quark_Mom[3]<<G4endl;
919 G4cout<<tmp<<" "<<tmp.mag()<<G4endl;
920 G4cout<<"2 str "<<Pstring2<<" "<<Pstring2.mag()<<" "<<Ystring2<<G4endl;
921 */
922 //--------------------------------
923  if(Ystring1 > Ystring2)
924  {
925  Pprojectile=Pstring1;
926  Ptarget =Pstring2;
927  } else
928  {
929  Pprojectile=Pstring2;
930  Ptarget =Pstring1;
931  }
932 
933 //-------------------------------------------------------
934 //G4cout<<"SumP CMS "<<Pprojectile+Ptarget<<" "<<SqrtS<<G4endl;
935 
936  Pprojectile.transform(toLab);
937  Ptarget.transform(toLab);
938 //G4cout<<"SumP Lab "<<Pprojectile+Ptarget<<" "<<SqrtS<<G4endl;
939 
940 // Calculation of the creation time ---------------------
941  projectile->SetTimeOfCreation(target->GetTimeOfCreation());
942  projectile->SetPosition(target->GetPosition());
943 
944 // Creation time and position of target nucleon were determined at
945 // ReggeonCascade() of G4FTFModel
946 // ------------------------------------------------------
947 
948 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
949 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
950  projectile->Set4Momentum(Pprojectile);
951 
952  target->Set4Momentum(Ptarget);
953 
954  projectile->IncrementCollisionCount(1);
955  target->IncrementCollisionCount(1);
956 
957  return true;
958  } // End of if(CandidatsN != 0)
959  }
960 
961 //============================================================
962 // Simulation of anti-quark-quark string creation
963 //
964  if(Ksi < (X_a+X_b+X_c+X_d)/Xannihilation)
965  {
966 //G4cout<<"Process d"<<G4endl;
967  G4int CandidatsN(0), CandAQ[9], CandQ[9];
968  G4int LeftAQ(0), LeftQ(0);
969 //------------------------------------------------------------
970  for(G4int iAQ1=0; iAQ1<3; iAQ1++)
971  {
972  for(G4int iAQ2=0; iAQ2<3; iAQ2++)
973  {
974  if(iAQ1 != iAQ2)
975  {
976  for(G4int iQ1=0; iQ1<3; iQ1++)
977  {
978  for(G4int iQ2=0; iQ2<3; iQ2++)
979  {
980  if(iQ1 != iQ2)
981  {
982  if((-AQ[iAQ1] == Q[iQ1]) && (-AQ[iAQ2] == Q[iQ2]))
983  {
984  if((iAQ1 == 0) && (iAQ2 == 1)){CandAQ[CandidatsN]=2;}
985  if((iAQ1 == 1) && (iAQ2 == 0)){CandAQ[CandidatsN]=2;}
986 
987  if((iAQ1 == 0) && (iAQ2 == 2)){CandAQ[CandidatsN]=1;}
988  if((iAQ1 == 2) && (iAQ2 == 0)){CandAQ[CandidatsN]=1;}
989 
990  if((iAQ1 == 1) && (iAQ2 == 2)){CandAQ[CandidatsN]=0;}
991  if((iAQ1 == 2) && (iAQ2 == 1)){CandAQ[CandidatsN]=0;}
992 //----------------------------------------------------------------
993  if((iQ1 == 0) && (iQ2 == 1)){CandQ[CandidatsN]=2;}
994  if((iQ1 == 1) && (iQ2 == 0)){CandQ[CandidatsN]=2;}
995 
996  if((iQ1 == 0) && (iQ2 == 2)){CandQ[CandidatsN]=1;}
997  if((iQ1 == 2) && (iQ2 == 0)){CandQ[CandidatsN]=1;}
998 
999  if((iQ1 == 1) && (iQ2 == 2)){CandQ[CandidatsN]=0;}
1000  if((iQ1 == 2) && (iQ2 == 1)){CandQ[CandidatsN]=0;}
1001  CandidatsN++;
1002  }//--------------------------
1003  } //end of if(jQ1 != jQ2)
1004  } //end of for(G4int jQ2=0; j<3; j++)
1005  } //end of for(G4int jQ=0; j<3; j++)
1006  } //end of if(iAQ1 != iAQ2)
1007  } //end of for(G4int iAQ2=0; i<3; i++)
1008  } //end of for(G4int iAQ1=0; i<3; i++)
1009 //------------------------------------------------------------
1010 
1011  if(CandidatsN != 0)
1012  {
1013  G4int SampledCase=CLHEP::RandFlat::shootInt(G4long(CandidatsN));
1014 
1015  LeftAQ=AQ[CandAQ[SampledCase]];
1016 
1017  LeftQ =Q[CandQ[SampledCase]];
1018 
1019 // --------------- Set the string properties ---------------
1020 //G4cout<<"Left Aq Q "<<LeftAQ<<" "<<LeftQ<<G4endl;
1021 
1022  projectile->SplitUp();
1023 
1024 // projectile->SetFirstParton(LeftAQ);
1025 // projectile->SetSecondParton(LeftQ);
1026  projectile->SetFirstParton(LeftQ);
1027  projectile->SetSecondParton(LeftAQ);
1028 
1029  projectile->SetStatus(1);
1030  target->SetStatus(3); // The target nucleon has annihilated
1031 
1032  Pprojectile.setPx(0.); // VU Mar1
1033  Pprojectile.setPy(0.); // Vu Mar1
1034  Pprojectile.setPz(0.);
1035  Pprojectile.setE(SqrtS);
1036  Pprojectile.transform(toLab);
1037 
1038 // Calculation of the creation time ---------------------
1039  projectile->SetTimeOfCreation(target->GetTimeOfCreation());
1040  projectile->SetPosition(target->GetPosition());
1041 // Creation time and position of target nucleon were determined at
1042 // ReggeonCascade() of G4FTFModel
1043 // ------------------------------------------------------
1044 
1045 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
1046 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
1047  projectile->Set4Momentum(Pprojectile);
1048 
1049  projectile->IncrementCollisionCount(1);
1050  return true;
1051  } // end of if(CandidatsN != 0)
1052  } // if(Ksi < (X_a+X_b+X_c+X_d/Xannihilation)
1053 
1054 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1055 //G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
1056 return true;
1057 }
1058 
1059 // ---------------------------------------------------------------------
1060 G4double G4FTFAnnihilation::ChooseX(G4double Alpha, G4double Beta) const
1061 {
1062 // If for sampling Xs other values of Alfa and Beta instead of 0.5 will be choose
1063 // the method will be implemented
1064  G4double tmp=Alpha*Beta;
1065  tmp*=1.;
1066  return 0.5;
1067 }
1068 
1069 // ---------------------------------------------------------------------
1070 G4ThreeVector G4FTFAnnihilation::GaussianPt(G4double AveragePt2,
1071  G4double maxPtSquare) const
1072 { // @@ this method is used in FTFModel as well. Should go somewhere common!
1073 
1074  G4double Pt2(0.);
1075  if(AveragePt2 <= 0.) {Pt2=0.;}
1076  else
1077  {
1078  Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
1079  (std::exp(-maxPtSquare/AveragePt2)-1.));
1080  }
1081  G4double Pt=std::sqrt(Pt2);
1082  G4double phi=G4UniformRand() * twopi;
1083  return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
1084 }
1085 
1086 
1087 // ---------------------------------------------------------------------
1088 void G4FTFAnnihilation::UnpackBaryon(G4int IdPDG,
1089  G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09
1090  {
1091  G4int AbsId=std::abs(IdPDG);
1092 
1093  Q1 = AbsId / 1000;
1094  Q2 = (AbsId % 1000) / 100;
1095  Q3 = (AbsId % 100) / 10;
1096 
1097  if(IdPDG < 0 ) {Q1=-Q1; Q2=-Q2; Q3=-Q3;} // Anti-baryon
1098 
1099  return;
1100  }
1101 
1102 // ---------------------------------------------------------------------
1103 G4FTFAnnihilation::G4FTFAnnihilation(const G4FTFAnnihilation &)
1104 {
1105  throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation copy contructor not meant to be called");
1106 }
1107 
1108 
1109 G4FTFAnnihilation::~G4FTFAnnihilation()
1110 {
1111 }
1112 
1113 
1114 const G4FTFAnnihilation & G4FTFAnnihilation::operator=(const G4FTFAnnihilation &)
1115 {
1116  throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation = operator not meant to be called");
1117  //return *this; //A.R. 25-Jul-2012 : fix Coverity
1118 }
1119 
1120 
1121 int G4FTFAnnihilation::operator==(const G4FTFAnnihilation &) const
1122 {
1123  throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation == operator not meant to be called");
1124  //return false; //A.R. 25-Jul-2012 : fix Coverity
1125 }
1126 
1127 int G4FTFAnnihilation::operator!=(const G4FTFAnnihilation &) const
1128 {
1129  throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator not meant to be called");
1130  //return true; //A.R. 25-Jul-2012 : fix Coverity
1131 }
Generated on Sat May 25 2013 14:34:05 for Geant4 by   doxygen 1.8.4