Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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
70G4FTFAnnihilation::G4FTFAnnihilation()
71{
72}
73
74// ---------------------------------------------------------------------
75G4bool 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/*
477G4cout<<"Mom 1 string "<<G4endl;
478G4cout<<Quark_Mom[0]<<G4endl;
479G4cout<<Quark_Mom[3]<<G4endl;
480G4cout<<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/*
489G4cout<<"Mom 2 string "<<G4endl;
490G4cout<<Quark_Mom[1]<<G4endl;
491G4cout<<Quark_Mom[4]<<G4endl;
492G4cout<<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/*
501G4cout<<"Mom 3 string "<<G4endl;
502G4cout<<Quark_Mom[2]<<G4endl;
503G4cout<<Quark_Mom[5]<<G4endl;
504G4cout<<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/*
904G4cout<<"Mom 1 string "<<G4endl;
905G4cout<<Quark_Mom[0]<<G4endl;
906G4cout<<Quark_Mom[2]<<G4endl;
907G4cout<<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/*
916G4cout<<"Mom 2 string "<<G4endl;
917G4cout<<Quark_Mom[1]<<G4endl;
918G4cout<<Quark_Mom[3]<<G4endl;
919G4cout<<tmp<<" "<<tmp.mag()<<G4endl;
920G4cout<<"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;
1056return true;
1057}
1058
1059// ---------------------------------------------------------------------
1060G4double 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// ---------------------------------------------------------------------
1070G4ThreeVector 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// ---------------------------------------------------------------------
1088void 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// ---------------------------------------------------------------------
1103G4FTFAnnihilation::G4FTFAnnihilation(const G4FTFAnnihilation &)
1104{
1105 throw G4HadronicException(__FILE__, __LINE__, "G4FTFAnnihilation copy contructor not meant to be called");
1106}
1107
1108
1109G4FTFAnnihilation::~G4FTFAnnihilation()
1110{
1111}
1112
1113
1114const 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
1121int 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
1127int 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}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4double
double G4double
Definition: G4Types.hh:64
G4long
long G4long
Definition: G4Types.hh:68
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::getZ
double getZ() const
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:162
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:174
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
CLHEP::RandFlat::shootInt
static long shootInt(long n)
G4FTFAnnihilation::Annihilate
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
Definition: G4FTFAnnihilation.cc:76
G4FTFParameters::GetAveragePt2
G4double GetAveragePt2()
Definition: G4FTFParameters.hh:377
G4VSplitableHadron::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition: G4VSplitableHadron.hh:152
G4VSplitableHadron::SetStatus
void SetStatus(const G4int aStatus)
Definition: G4VSplitableHadron.hh:172
G4VSplitableHadron::SplitUp
virtual void SplitUp()=0
Definition: G4VSplitableHadron.cc:114
G4VSplitableHadron::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4VSplitableHadron.hh:142
G4VSplitableHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4VSplitableHadron.hh:127
G4VSplitableHadron::SetSecondParton
virtual void SetSecondParton(G4int PDGcode)=0
G4VSplitableHadron::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4VSplitableHadron.hh:132
G4VSplitableHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4VSplitableHadron.hh:167
G4VSplitableHadron::IncrementCollisionCount
void IncrementCollisionCount(G4int aCount)
Definition: G4VSplitableHadron.hh:147
G4VSplitableHadron::SetFirstParton
virtual void SetFirstParton(G4int PDGcode)=0
G4VSplitableHadron::SetPosition
void SetPosition(const G4ThreeVector &aPosition)
Definition: G4VSplitableHadron.hh:162
sqr
T sqr(const T &x)
Definition: templates.hh:145