| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28
29 // -------------------------------------------
30 // GEANT 4 class implemetation file
31 //
32 // ---------------- G4FTFAnnihilation ---
33 // by V. Uzhinsky, Spring 2011.
34 // Take a projectile and a targ
35 // make annihilation or re-orangement o
36 // Ideas of Quark-Gluon-String model my A.
37 // are implemented.
38 // -------------------------------------------
39
40 #include "globals.hh"
41 #include "Randomize.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44
45 #include "G4DiffractiveSplitableHadron.hh"
46 #include "G4DiffractiveExcitation.hh"
47 #include "G4FTFParameters.hh"
48 #include "G4ElasticHNScattering.hh"
49 #include "G4FTFAnnihilation.hh"
50
51 #include "G4LorentzRotation.hh"
52 #include "G4RotationMatrix.hh"
53 #include "G4ThreeVector.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4VSplitableHadron.hh"
56 #include "G4ExcitedString.hh"
57 #include "G4ParticleTable.hh"
58 #include "G4Neutron.hh"
59 #include "G4ParticleDefinition.hh"
60
61 #include "G4Exp.hh"
62 #include "G4Log.hh"
63 #include "G4Pow.hh"
64
65 //#include "UZHI_diffraction.hh"
66
67 #include "G4ParticleTable.hh"
68
69 //============================================
70
71 //#define debugFTFannih
72
73
74 //============================================
75
76 G4FTFAnnihilation::G4FTFAnnihilation() {}
77
78
79 //============================================
80
81 G4FTFAnnihilation::~G4FTFAnnihilation() {}
82
83
84 //============================================
85
86 G4bool G4FTFAnnihilation::Annihilate( G4VSplit
87 G4VSplit
88 G4VSplit
89 G4FTFPar
90
91 #ifdef debugFTFannih
92 G4cout << "---------------------------- Anni
93 #endif
94
95 CommonVariables common;
96
97 // Projectile parameters
98 common.Pprojectile = projectile->Get4Momentu
99 G4int ProjectilePDGcode = projectile->GetDef
100 if ( ProjectilePDGcode > 0 ) {
101 target->SetStatus( 3 );
102 return false;
103 }
104 G4double M0projectile2 = common.Pprojectile.
105
106 // Target parameters
107 G4int TargetPDGcode = target->GetDefinition(
108 common.Ptarget = target->Get4Momentum();
109 G4double M0target2 = common.Ptarget.mag2();
110
111 #ifdef debugFTFannih
112 G4cout << "PDG codes " << ProjectilePDGcode
113 << "Pprojec " << common.Pprojectile <
114 << "Ptarget " << common.Ptarget <<
115 << "M0 proj target " << std::sqrt( M0
116 << " " << std::sqrt( M0target2 ) << G
117 #endif
118
119 // Kinematical properties of the interaction
120 G4LorentzVector Psum = common.Pprojectile +
121 common.S = Psum.mag2();
122 common.SqrtS = std::sqrt( common.S );
123 #ifdef debugFTFannih
124 G4cout << "Psum SqrtS S " << Psum << " " <<
125 #endif
126
127 // Transform momenta to cms and then rotate
128 G4LorentzRotation toCms( -1*Psum.boostVector
129 G4LorentzVector Ptmp( toCms*common.Pprojecti
130 toCms.rotateZ( -1*Ptmp.phi() );
131 toCms.rotateY( -1*Ptmp.theta() );
132 common.toLab = toCms.inverse();
133
134 if ( G4UniformRand() <= G4Pow::GetInstance()
135 common.RotateStrings = true;
136 common.RandomRotation.rotateZ( 2.0*pi*G4Un
137 common.RandomRotation.rotateY( std::acos(
138 common.RandomRotation.rotateZ( 2.0*pi*G4Un
139 }
140
141 G4double MesonProdThreshold = projectile->Ge
142 target->GetDef
143 ( 2.0*140.0 +
144 G4double Prel2 = sqr(common.S) + sqr(M0proje
145 - 2.0*( common.S*(M0project
146 Prel2 /= common.S;
147 G4double X_a = 0.0, X_b = 0.0, X_c = 0.0, X_
148 if ( Prel2 <= 0.0 ) {
149 // Annihilation at rest! Values are copied
150 X_a = 625.1; // mb // 3-shirt diagram
151 X_b = 0.0; // mb // anti-quark-quark
152 X_c = 49.989; // mb // 2 Q-Qbar string
153 X_d = 6.614; // mb // One Q-Qbar strin
154 #ifdef debugFTFannih
155 G4cout << "Annih at Rest X a b c d " << X_
156 << G4endl;
157 #endif
158 } else { // Annihilation in flight!
159 G4double FlowF = 1.0 / std::sqrt( Prel2 )*
160 // Process cross sections
161 X_a = 25.0*FlowF; // mb 3-shirt diagram
162 if ( common.SqrtS < MesonProdThreshold ) {
163 X_b = 3.13 + 140.0*G4Pow::GetInstance()-
164 } else {
165 X_b = 6.8*GeV / common.SqrtS; // mb ant
166 }
167 if ( projectile->GetDefinition()->GetPDGMa
168 > common.SqrtS ) {
169 X_b = 0.0;
170 }
171 // This can be in an interaction of low en
172 X_c = 2.0 * FlowF * sqr( projectile->GetDe
173 target->GetDefini
174
175 X_d = 23.3*GeV*GeV / common.S; // mb anti-
176 #ifdef debugFTFannih
177 G4cout << "Annih in Flight X a b c d " <<
178 << G4endl << "SqrtS MesonProdThresh
179 << G4endl;
180 #endif
181 }
182
183 G4bool isUnknown = false;
184 if ( TargetPDGcode == 2212 || TargetPDGcod
185 if ( ProjectilePDGcode == -2212 ||
186 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; //
187 } else if ( ProjectilePDGcode == -2112 ||
188 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; //
189 } else if ( ProjectilePDGcode == -3122 ) {
190 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; //
191 } else if ( ProjectilePDGcode == -3112 ) {
192 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; //
193 } else if ( ProjectilePDGcode == -3212 ) {
194 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; //
195 } else if ( ProjectilePDGcode == -3222 ) {
196 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; //
197 } else if ( ProjectilePDGcode == -3312 ) {
198 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; //
199 } else if ( ProjectilePDGcode == -3322 ) {
200 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; //
201 } else if ( ProjectilePDGcode == -3334 ) {
202 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; //
203 } else {
204 isUnknown = true;
205 }
206 } else if ( TargetPDGcode == 2112 || Targe
207 if ( ProjectilePDGcode == -2212 ||
208 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; //
209 } else if ( ProjectilePDGcode == -2112 ||
210 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; //
211 } else if ( ProjectilePDGcode == -3122 ) {
212 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; //
213 } else if ( ProjectilePDGcode == -3112 ) {
214 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; //
215 } else if ( ProjectilePDGcode == -3212 ) {
216 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; //
217 } else if ( ProjectilePDGcode == -3222 ) {
218 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; //
219 } else if ( ProjectilePDGcode == -3312 ) {
220 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; //
221 } else if ( ProjectilePDGcode == -3322 ) {
222 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; //
223 } else if ( ProjectilePDGcode == -3334 ) {
224 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; //
225 } else {
226 isUnknown = true;
227 }
228 } else {
229 isUnknown = true;
230 }
231 if ( isUnknown ) {
232 G4cout << "Unknown anti-baryon for FTF ann
233 << ProjectilePDGcode << " " << Targ
234 }
235
236 #ifdef debugFTFannih
237 G4cout << "Annih Actual X a b c d " << X_a <
238 #endif
239
240 G4double Xannihilation = X_a + X_b + X_c + X
241
242 // Projectile unpacking
243 UnpackBaryon( ProjectilePDGcode, common.AQ[0
244
245 // Target unpacking
246 UnpackBaryon( TargetPDGcode, common.Q[0], co
247
248 G4double Ksi = G4UniformRand();
249
250 if ( Ksi < X_a / Xannihilation ) {
251 return Create3QuarkAntiQuarkStrings( proje
252 }
253
254 G4int resultCode = 99;
255 if ( Ksi < (X_a + X_b) / Xannihilation ) {
256 resultCode = Create1DiquarkAntiDiquarkStri
257 if ( resultCode == 0 ) {
258 return true;
259 } else if ( resultCode == 99 ) {
260 return false;
261 }
262 }
263
264 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilati
265 resultCode = Create2QuarkAntiQuarkStrings(
266 if ( resultCode == 0 ) {
267 return true;
268 } else if ( resultCode == 99 ) {
269 return false;
270 }
271 }
272
273 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xanni
274 return Create1QuarkAntiQuarkString( projec
275 }
276
277 return true;
278 }
279
280
281 //--------------------------------------------
282
283 G4bool G4FTFAnnihilation::
284 Create3QuarkAntiQuarkStrings( G4VSplitableHadr
285 G4VSplitableHadr
286 G4VSplitableHadr
287 G4FTFParameters*
288 G4FTFAnnihilatio
289 // Simulation of 3 anti-quark - quark string
290
291 #ifdef debugFTFannih
292 G4cout << "Process a, 3 shirt diagram" << G4
293 #endif
294
295 // Sampling kinematical properties of quark.
296
297 const G4int maxNumberOfLoops = 1000;
298 G4double MassQ2 = 0.0; // Simp
299 // In p
300 G4double Quark_Xs[6];
301 G4ThreeVector Quark_Mom[6];
302
303 G4double Alfa_R = 0.5;
304 G4double AveragePt2 = 200.0*200.0, maxPtSqua
305 G4double ScaleFactor = 1.0;
306 G4double Alfa = 0.0, Beta = 0.0;
307
308 G4int NumberOfTries = 0, loopCounter = 0;
309
310 do {
311 // Sampling X's of anti-baryon and baryon
312 G4double x1 = 0.0, x2 = 0.0, x3 = 0.0;
313 G4double Product = 1.0;
314 for ( G4int iCase = 0; iCase < 2; ++iCase
315
316 G4double r1 = G4UniformRand(), r2 = G4Un
317 if ( Alfa_R == 1.0 ) {
318 x1 = 1.0 - std::sqrt( r1 );
319 x2 = (1.0 - x1) * r2;
320 } else {
321 x1 = sqr( r1 );
322 x2 = (1.0 - x1) * sqr( std::sin( pi/2.
323 }
324 x3 = 1.0 - x1 - x2;
325
326 G4int index = iCase*3; // 0 for anti-ba
327 Quark_Xs[index] = x1; Quark_Xs[index+1]
328 Product *= (x1*x2*x3);
329 }
330
331 if ( Product == 0.0 ) continue;
332
333 ++NumberOfTries;
334 if ( NumberOfTries == 100*(NumberOfTries/1
335 // After a large number of tries, it is
336 ScaleFactor /= 2.0;
337 AveragePt2 *= ScaleFactor;
338 }
339
340 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
341 for ( G4int i = 0; i < 6; ++i ) {
342 Quark_Mom [i] = GaussianPt( AveragePt2,
343 PtSum += Quark_Mom[i];
344 }
345
346 PtSum /= 6.0;
347 Alfa = 0.0; Beta = 0.0;
348
349 for ( G4int i = 0; i < 6; ++i ) { // Loop
350 Quark_Mom[i] -= PtSum;
351
352 G4double val = ( Quark_Mom[i].mag2() + M
353 if ( i < 3 ) { // anti-baryon
354 Alfa += val;
355 } else { // baryon (iCase == 1)
356 Beta += val;
357 }
358 }
359
360 } while ( ( std::sqrt( Alfa ) + std::sqrt( B
361 ++loopCounter < maxNumberOfLoops )
362
363 if ( loopCounter >= maxNumberOfLoops ) {
364 return false;
365 }
366
367 G4double DecayMomentum2 = sqr(common.S) + sq
368 - 2.0*( common.S*(
369
370 G4double WminusTarget = 0.0, WplusProjectile
371 WminusTarget = ( common.S - Alfa + Beta + st
372 WplusProjectile = common.SqrtS - Beta/Wminus
373
374 for ( G4int iCase = 0; iCase < 2; ++iCase )
375 G4int index = iCase*3; //
376 G4double w = WplusProjectile; //
377 if ( iCase == 1 ) w = - WminusTarget; //
378 for ( G4int i = 0; i < 3; ++i ) {
379 G4double Pz = w * Quark_Xs[index+i] / 2.
380 ( Quark_Mom[index+i].mag2(
381 ( 2.0 * w * Quark_Xs[index
382 Quark_Mom[index+i].setZ( Pz );
383 }
384 }
385
386 // Sampling of anti-quark order in projectil
387 G4int SampledCase = (G4int)G4RandFlat::shoot
388 G4int Tmp1 = 0, Tmp2 = 0;
389 switch ( SampledCase ) {
390 case 1 : Tmp1 = common.AQ[1]; common.AQ[1]
391 case 2 : Tmp1 = common.AQ[0]; common.AQ[0]
392 case 3 : Tmp1 = common.AQ[0]; Tmp2
393 common.AQ[1] = Tmp1; comm
394 case 4 : Tmp1 = common.AQ[0]; Tmp2
395 common.AQ[1] = common.AQ[2]; comm
396 case 5 : Tmp1 = common.AQ[0]; Tmp2
397 common.AQ[1] = Tmp2; comm
398 }
399
400 // Set the string properties
401 // An anti quark - quark pair can have the q
402 // or a vector meson: the last digit of the
403 // For simplicity only scalar is considered
404 G4int NewCode = 0, antiQuark = 0, quark = 0;
405 G4ParticleDefinition* TestParticle = nullptr
406 for ( G4int iString = 0; iString < 3; ++iStr
407 if ( iString == 0 ) {
408 antiQuark = common.AQ[0]; quark = commo
409 projectile->SetFirstParton( antiQuark );
410 projectile->SetSecondParton( quark );
411 projectile->SetStatus( 0 );
412 } else if ( iString == 1 ) {
413 quark = common.Q[1]; antiQuark = common
414 target->SetFirstParton( quark );
415 target->SetSecondParton( antiQuark );
416 target->SetStatus( 0 );
417 } else { // iString == 2
418 antiQuark = common.AQ[2]; quark = commo
419 }
420 G4int absAntiQuark = std::abs( antiQuark )
421 G4double aKsi = G4UniformRand();
422 if ( absAntiQuark == absQuark ) {
423 if ( absAntiQuark != 3 ) { // Not yet c
424 NewCode = 111; // Pi0-meson
425 if ( aKsi < 0.5 ) {
426 NewCode = 221; // Eta -meso
427 if ( aKsi < 0.25 ) {
428 NewCode = 331; // Eta'-meso
429 }
430 }
431 } else {
432 NewCode = 221; // Eta -meso
433 if ( aKsi < 0.5 ) {
434 NewCode = 331; // Eta'-meso
435 }
436 }
437 } else { // Vector mesons - rho, omega, p
438 if ( absAntiQuark > absQuark ) {
439 NewCode = absAntiQuark*100 + absQuark*
440 } else {
441 NewCode = absQuark*100 + absAntiQuark*
442 }
443 }
444 if ( iString == 2 ) AdditionalString = new
445 TestParticle = G4ParticleTable::GetParticl
446 if ( ! TestParticle ) return false;
447 if ( iString == 0 ) {
448 projectile->SetDefinition( TestParticle
449 theParameters->SetProjMinDiffMass( 0.5 )
450 theParameters->SetProjMinNonDiffMass( 0.
451 } else if ( iString == 1 ) {
452 target->SetDefinition( TestParticle );
453 theParameters->SetTarMinDiffMass( 0.5 );
454 theParameters->SetTarMinNonDiffMass( 0.5
455 } else { // iString == 2
456 AdditionalString->SetDefinition( TestPar
457 AdditionalString->SetFirstParton( common
458 AdditionalString->SetSecondParton( commo
459 AdditionalString->SetStatus( 0 );
460 }
461 } // End of the for loop over the 3 string
462
463 // 1st string AQ[0]-Q[0], 2nd string AQ[1]-Q
464
465 G4LorentzVector Pstring1, Pstring2, Pstring3
466 G4int QuarkOrder[3] = { 0 };
467 G4double YstringMax = 0.0, YstringMin = 0.0;
468 for ( G4int i = 0; i < 3; ++i ) {
469 G4ThreeVector tmp = Quark_Mom[i] + Quark_M
470 G4LorentzVector Pstring( tmp, std::sqrt( Q
471 std::sqrt( Q
472 // Add protection for rapidity = 0.5*ln(
473 G4double Ystring = 0.0;
474 if ( Pstring.e() > 1.0e-30 ) {
475 if ( Pstring.e() + Pstring.pz() < 1.0e-3
476 Ystring = -1.0e30; // A very large n
477 if ( Pstring.e() - Pstring.pz() < 1.0e
478 Ystring = 1.0e30; // A very large p
479 } else { // Normal case
480 Ystring = Pstring.rapidity();
481 }
482 }
483 }
484 // Keep ordering in rapidity: "1" highest,
485 if ( i == 0 ) {
486 Pstring1 = Pstring; YstringMax = Yst
487 QuarkOrder[0] = 0;
488 } else if ( i == 1 ) {
489 if ( Ystring > YstringMax ) {
490 Pstring2 = Pstring1; YstringMin = Yst
491 Pstring1 = Pstring; YstringMax = Yst
492 QuarkOrder[0] = 1; QuarkOrder[1] = 0;
493 } else {
494 Pstring2 = Pstring; YstringMin = Yst
495 QuarkOrder[1] = 1;
496 }
497 } else { // i == 2
498 if ( Ystring > YstringMax ) {
499 Pstring3 = Pstring2;
500 Pstring2 = Pstring1;
501 Pstring1 = Pstring;
502 QuarkOrder[1] = QuarkOrder[0];
503 QuarkOrder[2] = QuarkOrder[1];
504 QuarkOrder[0] = 2;
505 } else if ( Ystring > YstringMin ) {
506 Pstring3 = Pstring2;
507 Pstring2 = Pstring;
508 } else {
509 Pstring3 = Pstring;
510 QuarkOrder[2] = 2;
511 }
512 }
513 }
514
515 G4LorentzVector Quark_4Mom[6];
516 for ( G4int i = 0; i < 6; ++i ) {
517 Quark_4Mom[i] = G4LorentzVector( Quark_Mom
518 if ( common.RotateStrings ) Quark_4Mom[i]
519 Quark_4Mom[i].transform( common.toLab );
520 }
521
522 projectile->Splitting();
523 projectile->GetNextAntiParton()->Set4Momentu
524 projectile->GetNextParton()->Set4Momentum( Q
525
526 target->Splitting();
527 target->GetNextParton()->Set4Momentum( Quark
528 target->GetNextAntiParton()->Set4Momentum( Q
529
530 AdditionalString->Splitting();
531 AdditionalString->GetNextAntiParton()->Set4M
532 AdditionalString->GetNextParton()->Set4Momen
533
534 common.Pprojectile = Pstring1; //
535 common.Ptarget = Pstring3; //
536 G4LorentzVector LeftString( Pstring2 ); //
537
538 if ( common.RotateStrings ) {
539 common.Pprojectile *= common.RandomRotatio
540 common.Ptarget *= common.RandomRotatio
541 LeftString *= common.RandomRotatio
542 }
543
544 common.Pprojectile.transform( common.toLab )
545 common.Ptarget.transform( common.toLab );
546 LeftString.transform( common.toLab );
547
548 // Calculation of the creation time
549 // Creation time and position of target nucl
550 projectile->SetTimeOfCreation( target->GetTi
551 projectile->SetPosition( target->GetPosition
552 AdditionalString->SetTimeOfCreation( target-
553 AdditionalString->SetPosition( target->GetPo
554
555 projectile->Set4Momentum( common.Pprojectile
556 AdditionalString->Set4Momentum( LeftString )
557 target->Set4Momentum( common.Ptarget );
558
559 projectile->IncrementCollisionCount( 1 );
560 AdditionalString->IncrementCollisionCount( 1
561 target->IncrementCollisionCount( 1 );
562
563 return true;
564 }
565
566
567 //--------------------------------------------
568
569 G4int G4FTFAnnihilation::
570 Create1DiquarkAntiDiquarkString( G4VSplitableH
571 G4VSplitableH
572 G4FTFAnnihila
573 // Simulation of anti-diquark-diquark string
574 // This method returns an integer code - ins
575 // "0" : successfully ended and nothing el
576 // "1" : successfully completed, but the w
577 // "99" : unsuccessfully ended, nothing els
578
579 #ifdef debugFTFannih
580 G4cout << "Process b, quark - anti-quark ann
581 #endif
582
583 G4int CandidatsN = 0, CandAQ[9][2] = {}, Can
584 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { //
585 for ( G4int iQ = 0; iQ < 3; ++iQ ) { //
586 if ( -common.AQ[iAQ] == common.Q[iQ] ) {
587 // Here "0", "1", "2" means, respectiv
588 // of the (anti-baryon) projectile or
589 if ( iAQ == 0 ) { CandAQ[CandidatsN][0
590 if ( iAQ == 1 ) { CandAQ[CandidatsN][0
591 if ( iAQ == 2 ) { CandAQ[CandidatsN][0
592 if ( iQ == 0 ) { CandQ[CandidatsN][0]
593 if ( iQ == 1 ) { CandQ[CandidatsN][0]
594 if ( iQ == 2 ) { CandQ[CandidatsN][0]
595 ++CandidatsN;
596 }
597 }
598 }
599
600 // Remaining two (anti-)quarks that form the
601 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0,
602 if ( CandidatsN != 0 ) {
603 G4int SampledCase = (G4int)G4RandFlat::sho
604 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0
605 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1
606 LeftQ1 = common.Q[ CandQ[SampledCase][0]
607 LeftQ2 = common.Q[ CandQ[SampledCase][1]
608
609 // Build anti-diquark and diquark : the la
610 // of anti-quark - anti-quark and quark -
611 // or quarks are different. For simplicity
612 G4int Anti_DQ = 0, DQ = 0;
613 if ( std::abs( LeftAQ1 ) > std::abs( LeftA
614 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3
615 } else {
616 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3
617 }
618 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2
619 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3;
620 } else {
621 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3;
622 }
623
624 // Set the string properties
625 projectile->SetFirstParton( DQ );
626 projectile->SetSecondParton( Anti_DQ );
627
628 // It is assumed that quark and di-quark m
629 G4LorentzVector Pquark = G4LorentzVector(
630 G4LorentzVector Paquark = G4LorentzVector(
631
632 if ( common.RotateStrings ) {
633 Pquark *= common.RandomRotation;
634 Paquark *= common.RandomRotation;
635 }
636
637 Pquark.transform( common.toLab );
638 Paquark.transform( common.toLab );
639
640 projectile->GetNextParton()->Set4Momentum(
641 projectile->GetNextAntiParton()->Set4Momen
642
643 projectile->Splitting();
644
645 projectile->SetStatus( 0 );
646 target->SetStatus( 4 ); // The target nuc
647 common.Pprojectile.setPx( 0.0 );
648 common.Pprojectile.setPy( 0.0 );
649 common.Pprojectile.setPz( 0.0 );
650 common.Pprojectile.setE( common.SqrtS );
651 common.Pprojectile.transform( common.toLab
652
653 // Calculation of the creation time
654 // Creation time and position of target nu
655 projectile->SetTimeOfCreation( target->Get
656 projectile->SetPosition( target->GetPositi
657 projectile->Set4Momentum( common.Pprojecti
658
659 projectile->IncrementCollisionCount( 1 );
660 target->IncrementCollisionCount( 1 );
661
662 return 0; // Completed successfully: noth
663 } // End of if ( CandidatsN != 0 )
664
665 // If we allow the string to interact with o
666 // set up MinDiffrMass in Parameters, and as
667
668 return 1; // Successfully ended, but the wo
669 }
670
671
672 //--------------------------------------------
673
674 G4int G4FTFAnnihilation::
675 Create2QuarkAntiQuarkStrings( G4VSplitableHadr
676 G4VSplitableHadr
677 G4FTFParameters*
678 G4FTFAnnihilatio
679 // Simulation of 2 anti-quark-quark strings
680 // This method returns an integer code - ins
681 // "0" : successfully ended and nothing el
682 // "1" : successfully completed, but the w
683 // "99" : unsuccessfully ended, nothing els
684
685 #ifdef debugFTFannih
686 G4cout << "Process c, quark - anti-quark and
687 << G4endl;
688 #endif
689
690 // Sampling kinematical properties: 1st stri
691 G4ThreeVector Quark_Mom[4];
692 G4double Quark_Xs[4];
693 G4double AveragePt2 = 200.0*200.0, maxPtSqua
694 G4int NumberOfTries = 0, loopCounter = 0;
695 const G4int maxNumberOfLoops = 1000;
696 G4double Alfa = 0.0, Beta = 0.0;
697 G4double WminusTarget = 0.0, WplusProjectile
698 do {
699 // Sampling X's of the 2 quarks and 2 anti
700
701 G4double Product = 1.0;
702 for ( G4int iCase = 0; iCase < 2; ++iCase
703 G4double x = 0.0, r = G4UniformRand();
704 if ( Alfa_R == 1.0 ) {
705 if ( iCase == 0 ) { // first string
706 x = std::sqrt( r );
707 } else { // second string
708 x = 1.0 - std::sqrt( r );
709 }
710 } else {
711 x = sqr( std::sin( pi/2.0*r ) );
712 }
713 G4int index = iCase*2; // 0 for the fir
714 Quark_Xs[index] = x ; Quark_Xs[index+1]
715 Product *= x*(1.0-x);
716 }
717
718 if ( Product == 0.0 ) continue;
719
720 ++NumberOfTries;
721 if ( NumberOfTries == 100*(NumberOfTries/1
722 // After a large number of tries, it is
723 ScaleFactor /= 2.0;
724 AveragePt2 *= ScaleFactor;
725 }
726
727 G4ThreeVector PtSum( 0.0, 0.0, 0.0 );
728 for( G4int i = 0; i < 4; ++i ) {
729 Quark_Mom[i] = GaussianPt( AveragePt2, m
730 PtSum += Quark_Mom[i];
731 }
732
733 PtSum /= 4.0;
734 for ( G4int i = 0; i < 4; ++i ) {
735 Quark_Mom[i] -= PtSum;
736 }
737
738 Alfa = 0.0; Beta = 0.0;
739 for ( G4int iCase = 0; iCase < 2; ++iCase
740 G4int index = iCase * 2;
741 for ( G4int i = 0; i < 2; ++i ) {
742 G4double val = ( Quark_Mom[index+i].
743 if ( iCase == 0 ) { // first string
744 Alfa += val;
745 } else { // second strin
746 Beta += val;
747 }
748 }
749 }
750
751 } while ( ( std::sqrt( Alfa ) + std::sqrt( B
752 ++loopCounter < maxNumberOfLoops )
753
754 if ( loopCounter >= maxNumberOfLoops ) {
755 return 99; // unsuccessfully ended, nothi
756 }
757
758 G4double DecayMomentum2 = sqr(common.S) + sq
759 - 2.0*( common.S*(
760 WminusTarget = ( common.S - Alfa + Beta + st
761 WplusProjectile = common.SqrtS - Beta/Wminus
762
763 for ( G4int iCase = 0; iCase < 2; ++iCase )
764 G4int index = iCase*2; // 0 for the first
765 for ( G4int i = 0; i < 2; ++i ) {
766 G4double w = WplusProjectile; /
767 if ( iCase == 1 ) w = - WminusTarget; /
768 G4double Pz = w * Quark_Xs[index+i] / 2.
769 - ( Quark_Mom[index+i].mag
770 ( 2.0 * w * Quark_Xs[ind
771 Quark_Mom[index+i].setZ( Pz );
772 }
773 }
774
775 G4int CandidatsN = 0, CandAQ[9][2] = {}, Can
776 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0,
777 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { //
778 for ( G4int iQ = 0; iQ < 3; ++iQ ) { //
779 if ( -common.AQ[iAQ] == common.Q[iQ] ) {
780 // Here "0", "1", "2" means, respectiv
781 // of the (anti-baryon) projectile or
782 if ( iAQ == 0 ) { CandAQ[CandidatsN][0
783 if ( iAQ == 1 ) { CandAQ[CandidatsN][0
784 if ( iAQ == 2 ) { CandAQ[CandidatsN][0
785 if ( iQ == 0 ) { CandQ[CandidatsN][0]
786 if ( iQ == 1 ) { CandQ[CandidatsN][0]
787 if ( iQ == 2 ) { CandQ[CandidatsN][0]
788 ++CandidatsN;
789 }
790 }
791 }
792
793 if ( CandidatsN != 0 ) {
794 G4int SampledCase = (G4int)G4RandFlat::sho
795 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0
796 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1
797 if ( G4UniformRand() < 0.5 ) {
798 LeftQ1 = common.Q[ CandQ[SampledCase][0]
799 LeftQ2 = common.Q[ CandQ[SampledCase][1]
800 } else {
801 LeftQ2 = common.Q[ CandQ[SampledCase][0]
802 LeftQ1 = common.Q[ CandQ[SampledCase][1]
803 }
804
805 // Set the string properties
806 // An anti quark - quark pair can have the
807 // or a vector meson: the last digit of th
808 // For simplicity only scalar is considere
809 G4int NewCode = 0, antiQuark = 0, quark =
810 G4ParticleDefinition* TestParticle = nullp
811 for ( G4int iString = 0; iString < 2; ++iS
812 if ( iString == 0 ) {
813 antiQuark = LeftAQ1; quark = LeftQ1;
814 projectile->SetFirstParton( antiQuark
815 projectile->SetSecondParton( quark );
816 projectile->SetStatus( 0 );
817 } else { // iString == 1
818 quark = LeftQ2; antiQuark = LeftAQ2;
819 target->SetFirstParton( quark );
820 target->SetSecondParton( antiQuark );
821 target->SetStatus( 0 );
822 }
823 G4int absAntiQuark = std::abs( antiQuark
824 G4double aKsi = G4UniformRand();
825 if ( absAntiQuark == absQuark ) {
826 if ( absAntiQuark != 3 ) {
827 NewCode = 111; // Pi0-meson
828 if ( aKsi < 0.5 ) {
829 NewCode = 221; // Eta -meso
830 if ( aKsi < 0.25 ) {
831 NewCode = 331; // Eta'-meso
832 }
833 }
834 } else {
835 NewCode = 221; // Eta -meso
836 if ( aKsi < 0.5 ) {
837 NewCode = 331; // Eta'-meso
838 }
839 }
840 } else {
841 if ( absAntiQuark > absQuark ) {
842 NewCode = absAntiQuark*100 + absQuar
843 } else {
844 NewCode = absQuark*100 + absAntiQuar
845 }
846 }
847 TestParticle = G4ParticleTable::GetParti
848 if ( ! TestParticle ) return 99; // uns
849 if ( iString == 0 ) {
850 projectile->SetDefinition( TestParticl
851 theParameters->SetProjMinDiffMass( 0.5
852 theParameters->SetProjMinNonDiffMass(
853 } else { // iString == 1
854 target->SetDefinition( TestParticle );
855 theParameters->SetTarMinDiffMass( 0.5
856 theParameters->SetTarMinNonDiffMass( 0
857 }
858 } // End of loop over the 2 string cases
859
860 G4int QuarkOrder[2];
861 G4LorentzVector Pstring1, Pstring2;
862 G4double Ystring1 = 0.0, Ystring2 = 0.0;
863
864 for ( G4int iCase = 0; iCase < 2; ++iCase
865 G4ThreeVector tmp = Quark_Mom[iCase] + Q
866 G4LorentzVector Pstring( tmp, std::sqrt(
867 std::sqrt(
868 // Add protection for rapidity = 0.5*ln
869 G4double Ystring = 0.0;
870 if ( Pstring.e() > 1.0e-30 ) {
871 if ( Pstring.e() + Pstring.pz() < 1.0e
872 Ystring = -1.0e30; // A very large
873 if ( Pstring.e() - Pstring.pz() < 1.
874 Ystring = 1.0e30; // A very large
875 } else { // Normal case
876 Ystring = Pstring.rapidity();
877 }
878 }
879 }
880 if ( iCase == 0 ) { // For the first st
881 Pstring1 = Pstring; Ystring1 = Ystring
882 } else { // For the second s
883 Pstring2 = Pstring; Ystring2 = Ystring
884 }
885 }
886 if ( Ystring1 > Ystring2 ) {
887 common.Pprojectile = Pstring1; common.P
888 QuarkOrder[0] = 0; QuarkOrder[1] = 1;
889 } else {
890 common.Pprojectile = Pstring2; common.P
891 QuarkOrder[0] = 1; QuarkOrder[1] = 0;
892 }
893
894 if ( common.RotateStrings ) {
895 common.Pprojectile *= common.RandomRotat
896 common.Ptarget *= common.RandomRotat
897 }
898
899 common.Pprojectile.transform( common.toLab
900 common.Ptarget.transform( common.toLab );
901
902 G4LorentzVector Quark_4Mom[4];
903 for ( G4int i = 0; i < 4; ++i ) {
904 Quark_4Mom[i] = G4LorentzVector( Quark_M
905 if ( common.RotateStrings ) Quark_4Mom[i
906 Quark_4Mom[i].transform( common.toLab );
907 }
908
909 projectile->Splitting();
910 projectile->GetNextAntiParton()->Set4Momen
911 projectile->GetNextParton()->Set4Momentum(
912
913 target->Splitting();
914 target->GetNextParton()->Set4Momentum( Qua
915 target->GetNextAntiParton()->Set4Momentum(
916
917 // Calculation of the creation time
918 // Creation time and position of target nu
919 projectile->SetTimeOfCreation( target->Get
920 projectile->SetPosition( target->GetPositi
921 projectile->Set4Momentum( common.Pprojecti
922 target->Set4Momentum( common.Ptarget );
923
924 projectile->IncrementCollisionCount( 1 );
925 target->IncrementCollisionCount( 1 );
926
927 return 0; // Completed successfully: noth
928 } // End of if ( CandidatsN != 0 )
929
930 return 1; // Successfully ended, but the wo
931 }
932
933
934 //--------------------------------------------
935
936 G4bool G4FTFAnnihilation::
937 Create1QuarkAntiQuarkString( G4VSplitableHadro
938 G4VSplitableHadro
939 G4FTFParameters*
940 G4FTFAnnihilation
941 // Simulation of anti-quark - quark string c
942
943 #ifdef debugFTFannih
944 G4cout << "Process d, only 1 quark - anti-qu
945 #endif
946
947 // Determine the set of candidates anti-quar
948 // Here "0", "1", "2" means, respectively, "
949 // of the (anti-baryon) projectile or (nucle
950 G4int CandidatsN = 0, CandAQ[36], CandQ[36];
951 G4int LeftAQ = 0, LeftQ = 0;
952 for ( G4int iAQ1 = 0; iAQ1 < 3; ++iAQ1 ) {
953 for ( G4int iAQ2 = 0; iAQ2 < 3; ++iAQ2 ) {
954 if ( iAQ1 != iAQ2 ) {
955 for ( G4int iQ1 = 0; iQ1 < 3; ++iQ1 )
956 for ( G4int iQ2 = 0; iQ2 < 3; ++iQ2
957 if ( iQ1 != iQ2 ) {
958 if ( -common.AQ[iAQ1] == common.
959 if ( ( iAQ1 == 0 && i
960 CandAQ[CandidatsN] = 2;
961 } else if ( ( iAQ1 == 0 && i
962 CandAQ[CandidatsN] = 1;
963 } else if ( ( iAQ1 == 1 && i
964 CandAQ[CandidatsN] = 0;
965 }
966 if ( ( iQ1 == 0 &&
967 CandQ[CandidatsN] = 2;
968 } else if ( ( iQ1 == 0 &&
969 CandQ[CandidatsN] = 1;
970 } else if ( ( iQ1 == 1 &&
971 CandQ[CandidatsN] = 0;
972 }
973 ++CandidatsN;
974 }
975 }
976 }
977 }
978 }
979 }
980 }
981
982 if ( CandidatsN != 0 ) {
983 G4int SampledCase = (G4int)G4RandFlat::sho
984 LeftAQ = common.AQ[ CandAQ[SampledCase] ];
985 LeftQ = common.Q[ CandQ[SampledCase] ];
986
987 // Set the string properties
988 projectile->SetFirstParton( LeftQ );
989 projectile->SetSecondParton( LeftAQ );
990 projectile->SetStatus( 0 );
991 G4int aAQ = std::abs( LeftAQ ), aQ = std::
992 G4int NewCode = 0;
993 G4double aKsi = G4UniformRand();
994 // The string can have the quantum number
995 // of the PDG code is, respectively, 1 and
996 if ( aAQ == aQ ) {
997 if ( aAQ != 3 ) {
998 NewCode = 111; // Pi0-meson
999 if ( aKsi < 0.5 ) {
1000 NewCode = 221; // Eta -meson
1001 if ( aKsi < 0.25 ) {
1002 NewCode = 331; // Eta'-meson
1003 }
1004 }
1005 } else {
1006 NewCode = 221; // Eta -meson
1007 if ( aKsi < 0.5 ) {
1008 NewCode = 331; // Eta'-meson
1009 }
1010 }
1011 } else {
1012 if ( aAQ > aQ ) {
1013 NewCode = aAQ*100 + aQ*10 + 1; NewCod
1014 } else {
1015 NewCode = aQ*100 + aAQ*10 + 1; NewCod
1016 }
1017 }
1018
1019 G4ParticleDefinition* TestParticle = G4Pa
1020 if ( ! TestParticle ) return false;
1021 projectile->SetDefinition( TestParticle )
1022 theParameters->SetProjMinDiffMass( 0.5 );
1023 theParameters->SetProjMinNonDiffMass( 0.5
1024
1025 target->SetStatus( 4 ); // The target nu
1026 common.Pprojectile.setPx( 0.0 );
1027 common.Pprojectile.setPy( 0.0 );
1028 common.Pprojectile.setPz( 0.0 );
1029 common.Pprojectile.setE( common.SqrtS );
1030
1031 common.Pprojectile.transform( common.toLa
1032
1033 G4LorentzVector Pquark = G4LorentzVector
1034 G4LorentzVector Paquark = G4LorentzVector
1035
1036 if ( common.RotateStrings ) {
1037 Pquark *= common.RandomRotation; Paquar
1038 }
1039 Pquark.transform(common.toLab); projecti
1040 Paquark.transform(common.toLab); projecti
1041
1042 projectile->Splitting();
1043
1044 // Calculation of the creation time
1045 // Creation time and position of target n
1046 projectile->SetTimeOfCreation( target->Ge
1047 projectile->SetPosition( target->GetPosit
1048 projectile->Set4Momentum( common.Pproject
1049
1050 projectile->IncrementCollisionCount( 1 );
1051 target->IncrementCollisionCount( 1 );
1052
1053 return true;
1054 } // End of if ( CandidatsN != 0 )
1055
1056 return true;
1057 }
1058
1059
1060 //===========================================
1061
1062 G4double G4FTFAnnihilation::ChooseX( G4double
1063 // If for sampling Xs other values of Alfa
1064 // chosen the method will be implemented
1065 //G4double tmp = Alpha*Beta;
1066 //tmp *= 1.0;
1067 return 0.5;
1068 }
1069
1070
1071 //===========================================
1072
1073 G4ThreeVector G4FTFAnnihilation::GaussianPt(
1074 // @@ this method is used in FTFModel as w
1075 G4double Pt2 = 0.0;
1076 if ( AveragePt2 <= 0.0 ) {
1077 Pt2 = 0.0;
1078 } else {
1079 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor
1080 ( G4E
1081 }
1082 G4double Pt = std::sqrt( Pt2 );
1083 G4double phi = G4UniformRand() * twopi;
1084 return G4ThreeVector ( Pt*std::cos( phi ),
1085 }
1086
1087
1088 //===========================================
1089
1090 void G4FTFAnnihilation::UnpackBaryon( G4int I
1091 G4int AbsId = std::abs( IdPDG );
1092 Q1 = AbsId / 1000;
1093 Q2 = ( AbsId % 1000 ) / 100;
1094 Q3 = ( AbsId % 100 ) / 10;
1095 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 =
1096 return;
1097 }
1098
1099
1100 //===========================================
1101
1102 G4FTFAnnihilation::G4FTFAnnihilation( const G
1103 throw G4HadronicException( __FILE__, __LINE
1104 "G4FTFAnnihilati
1105 }
1106
1107
1108 //===========================================
1109
1110 const G4FTFAnnihilation & G4FTFAnnihilation::
1111 throw G4HadronicException( __FILE__, __LINE
1112 "G4FTFAnnihilati
1113 }
1114
1115
1116 //===========================================
1117
1118 G4bool G4FTFAnnihilation::operator==( const G
1119 throw G4HadronicException( __FILE__, __LINE
1120 "G4FTFAnnihilati
1121 }
1122
1123
1124 //===========================================
1125
1126 G4bool G4FTFAnnihilation::operator!=( const G
1127 throw G4HadronicException( __FILE__, __LINE
1128 "G4DiffractiveEx
1129 }
1130