Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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