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 // ---------------- G4DiffractiveExcitati 33 // by Gunter Folger, October 1998. 34 // diffractive Excitation used by strin 35 // Take a projectile and a target 36 // excite the projectile and targe 37 // Essential changed by V. Uzhinsky in Novemb 38 // in order to put it in a correspondence wit 39 // model. Variant of FRITIOF with nucleon de- 40 // Other changes by V.Uzhinsky in May 2007 we 41 // meson-nucleon interactions. Additional cha 42 // were introduced in December 2006. They tre 43 // processes more exactly. 44 // Correct treatment of the diffraction disso 45 // Mass distributions for resonances and uu-d 46 // and dd-diquarks suppression in neutrons we 47 // ------------------------------------------- 48 49 #include "globals.hh" 50 #include "Randomize.hh" 51 #include "G4PhysicalConstants.hh" 52 #include "G4SystemOfUnits.hh" 53 54 #include "G4DiffractiveExcitation.hh" 55 #include "G4FTFParameters.hh" 56 #include "G4ElasticHNScattering.hh" 57 58 #include "G4RotationMatrix.hh" 59 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleTable.hh" 61 #include "G4SampleResonance.hh" 62 #include "G4VSplitableHadron.hh" 63 #include "G4ExcitedString.hh" 64 #include "G4Neutron.hh" 65 66 #include "G4Exp.hh" 67 #include "G4Log.hh" 68 #include "G4Pow.hh" 69 70 //#include "G4ios.hh" 71 72 //============================================ 73 74 //#define debugFTFexcitation 75 //#define debug_heavyHadrons 76 77 //============================================ 78 79 G4DiffractiveExcitation::G4DiffractiveExcitati 80 81 82 //============================================ 83 84 G4DiffractiveExcitation::~G4DiffractiveExcitat 85 86 87 //============================================ 88 89 G4bool G4DiffractiveExcitation::ExciteParticip 90 91 92 93 94 #ifdef debugFTFexcitation 95 G4cout << G4endl << "FTF ExciteParticipants 96 #endif 97 98 CommonVariables common; 99 100 // Projectile parameters 101 common.Pprojectile = projectile->Get4Momentu 102 if ( common.Pprojectile.z() < 0.0 ) return f 103 common.ProjectilePDGcode = projectile->GetDe 104 common.absProjectilePDGcode = std::abs( comm 105 common.M0projectile = projectile->GetDefinit 106 G4double ProjectileRapidity = common.Pprojec 107 108 // Target parameters 109 common.Ptarget = target->Get4Momentum(); 110 common.TargetPDGcode = target->GetDefinition 111 common.absTargetPDGcode = std::abs( common.T 112 common.M0target = target->GetDefinition()->G 113 G4double TargetRapidity = common.Ptarget.rap 114 115 // Kinematical properties of the interaction 116 G4LorentzVector Psum = common.Pprojectile + 117 common.S = Psum.mag2(); 118 common.SqrtS = std::sqrt( common.S ); 119 120 // Check off-shellness of the participants 121 G4bool toBePutOnMassShell = true; //Uzhi Au 122 common.MminProjectile = common.BrW.GetMinimu 123 /* Uzhi Aug.2019 124 if ( common.M0projectile < common.MminProjec 125 toBePutOnMassShell = true; 126 common.M0projectile = common.BrW.SampleMas 127 128 129 } 130 */ 131 common.M0projectile2 = common.M0projectile * 132 common.ProjectileDiffStateMinMass = thePa 133 common.ProjectileNonDiffStateMinMass = thePa 134 if ( common.M0projectile > common.Projectile 135 common.ProjectileDiffStateMinMass = com 136 common.ProjectileNonDiffStateMinMass = com 137 if ( common.absProjectilePDGcode > 3000 ) 138 common.ProjectileDiffStateMinMass += 139 common.ProjectileNonDiffStateMinMass += 140 } 141 } 142 common.MminTarget = common.BrW.GetMinimumMas 143 /* Uzhi Aug.2019 144 if ( common.M0target < common.MminTarget ) { 145 toBePutOnMassShell = true; 146 common.M0target = common.BrW.SampleMass( t 147 t 148 + 149 } 150 */ 151 common.M0target2 = common.M0target * common. 152 common.TargetDiffStateMinMass = theParame 153 common.TargetNonDiffStateMinMass = theParame 154 if ( common.M0target > common.TargetDiffStat 155 common.TargetDiffStateMinMass = common. 156 common.TargetNonDiffStateMinMass = common. 157 if ( common.absTargetPDGcode > 3000 ) { / 158 common.TargetDiffStateMinMass += 140. 159 common.TargetNonDiffStateMinMass += 140. 160 } 161 }; 162 #ifdef debugFTFexcitation 163 G4cout << "Proj Targ PDGcodes " << common.Pr 164 << "Mprojectile Y " << common.Pproje 165 << "M0projectile Y " << common.M0proj 166 G4cout << "Mtarget Y " << common.Ptarge 167 << "M0target Y " << common.M0targ 168 G4cout << "Pproj " << common.Pprojectile << 169 #endif 170 171 // Transform momenta to cms and then rotate 172 common.toCms = G4LorentzRotation( -1 * Psum. 173 G4LorentzVector Ptmp = common.toCms * common 174 if ( Ptmp.pz() <= 0.0 ) return false; // "S 175 common.toCms.rotateZ( -1*Ptmp.phi() ); 176 common.toCms.rotateY( -1*Ptmp.theta() ); 177 common.toLab = common.toCms.inverse(); 178 common.Pprojectile.transform( common.toCms ) 179 common.Ptarget.transform( common.toCms ); 180 181 G4double SumMasses = common.M0projectile + c 182 #ifdef debugFTFexcitation 183 G4cout << "SqrtS " << common.SqrtS << G4 184 << " " << common.M0target << " " << S 185 #endif 186 if ( common.SqrtS < SumMasses ) return false 187 188 common.PZcms2 = ( sqr( common.S ) + sqr( com 189 - 2.0 * ( common.S * ( com 190 + common.M0proje 191 #ifdef debugFTFexcitation 192 G4cout << "PZcms2 after toBePutOnMassShell " 193 #endif 194 if ( common.PZcms2 < 0.0 ) return false; // 195 196 common.PZcms = std::sqrt( common.PZcms2 ); 197 if ( toBePutOnMassShell ) { 198 if ( common.Pprojectile.z() > 0.0 ) { 199 common.Pprojectile.setPz( common.PZcms 200 common.Ptarget.setPz( -common.PZcms 201 } else { 202 common.Pprojectile.setPz( -common.PZcms 203 common.Ptarget.setPz( common.PZcms 204 } 205 common.Pprojectile.setE( std::sqrt( common 206 + comm 207 + comm 208 + comm 209 common.Ptarget.setE( std::sqrt( common.M0t 210 + common.P 211 + common.P 212 + common.P 213 } 214 #ifdef debugFTFexcitation 215 G4cout << "Start --------------------" << G4 216 << " " << common.ProjectileDiffStateM 217 << G4endl 218 << "Targ M0 Mdif Mndif " << common.M0 219 << " " << common.TargetNonDiffStateMi 220 << "Proj CMS " << common.Pprojectile 221 #endif 222 223 // Check for possible quark exchange 224 ProjectileRapidity = common.Pprojectile.rapi 225 TargetRapidity = common.Ptarget.rapidity(); 226 G4double QeNoExc = theParameters->GetProcPro 227 G4double QeExc = theParameters->GetProcPro 228 theParameters->GetProcPro 229 common.ProbProjectileDiffraction = 230 theParameters->GetProcProb( 2, ProjectileR 231 common.ProbTargetDiffraction = 232 theParameters->GetProcProb( 3, ProjectileR 233 common.ProbOfDiffraction = common.ProbProjec 234 #ifdef debugFTFexcitation 235 G4cout << "Proc Probs " << QeNoExc << " " << 236 << common.ProbProjectileDiffraction < 237 << "ProjectileRapidity " << Projectil 238 #endif 239 240 if ( QeNoExc + QeExc + common.ProbProjectile 241 QeNoExc = 1.0 - QeExc - common.ProbProject 242 } 243 if ( QeExc + QeNoExc != 0.0 ) { 244 common.ProbExc = QeExc / ( QeExc + QeNoExc 245 } 246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) { 247 common.ProbProjectileDiffraction /= ( 1.0 248 common.ProbTargetDiffraction /= ( 1.0 249 } 250 #ifdef debugFTFexcitation 251 G4cout << "Proc Probs " << QeNoExc << " " << 252 << common.ProbProjectileDiffraction < 253 << "ProjectileRapidity " << Projectil 254 #endif 255 256 // Try out charge exchange 257 G4int returnCode = 1; 258 if ( G4UniformRand() < QeExc + QeNoExc ) { 259 returnCode = ExciteParticipants_doChargeEx 260 261 } 262 263 G4bool returnResult = false; 264 if ( returnCode == 0 ) { 265 returnResult = true; // Successfully ende 266 } else if ( returnCode == 1 ) { 267 268 common.ProbOfDiffraction = common.ProbProj 269 #ifdef debugFTFexcitation 270 G4cout << "Excitation -------------------- 271 << "Proj M0 MdMin MndMin " << commo 272 << common.ProjectileDiffStateMinMas 273 << G4endl 274 << "Targ M0 MdMin MndMin " << commo 275 << " " << common.TargetNonDiffState 276 << G4endl 277 << "Prob: ProjDiff TargDiff + Sum " 278 << common.ProbTargetDiffraction << 279 #endif 280 if ( common.ProbOfDiffraction != 0.0 ) { 281 common.ProbProjectileDiffraction /= comm 282 } else { 283 common.ProbProjectileDiffraction = 0.0; 284 } 285 #ifdef debugFTFexcitation 286 G4cout << "Prob: ProjDiff TargDiff + Sum " 287 << common.ProbTargetDiffraction << 288 #endif 289 common.ProjectileDiffStateMinMass2 = sq 290 common.ProjectileNonDiffStateMinMass2 = sq 291 common.TargetDiffStateMinMass2 = sq 292 common.TargetNonDiffStateMinMass2 = sq 293 // Choose between diffraction and non-diff 294 if ( G4UniformRand() < common.ProbOfDiffra 295 296 returnResult = ExciteParticipants_doDiff 297 298 } else { 299 300 returnResult = ExciteParticipants_doNonD 301 302 } 303 if ( returnResult ) { 304 common.Pprojectile += common.Qmomentum; 305 common.Ptarget -= common.Qmomentum; 306 // Transform back and update SplitableHa 307 common.Pprojectile.transform( common.toL 308 common.Ptarget.transform( common.toLab ) 309 #ifdef debugFTFexcitation 310 G4cout << "Mproj " << common.Pprojectile 311 << G4endl; 312 #endif 313 projectile->Set4Momentum( common.Pprojec 314 target->Set4Momentum( common.Ptarget ); 315 projectile->IncrementCollisionCount( 1 ) 316 target->IncrementCollisionCount( 1 ); 317 } 318 } 319 320 return returnResult; 321 } 322 323 //-------------------------------------------- 324 325 G4int G4DiffractiveExcitation:: 326 ExciteParticipants_doChargeExchange( G4VSplita 327 G4VSplita 328 G4FTFPara 329 G4Elastic 330 G4Diffrac 331 // First of the three utility methods used o 332 // it does the sampling for the charge-excha 333 // This method returns an integer code - ins 334 // "0" : successfully ended and nothing el 335 // "1" : successfully completed, but the w 336 // "99" : unsuccessfully ended, nothing els 337 G4int returnCode = 99; 338 339 G4double DeltaProbAtQuarkExchange = theParam 340 G4ParticleDefinition* TestParticle = 0; 341 G4double MtestPr = 0.0, MtestTr = 0.0; 342 343 #ifdef debugFTFexcitation 344 G4cout << "Q exchange ---------------------- 345 #endif 346 347 // The target hadron is always a nucleon (i. 348 // never an antinucleon), therefore only a q 349 // exchanged between the projectile hadron a 350 // we could get a quark-quark-antiquark syst 351 // This implies that any projectile meson or 352 // a constituent quark in all cases - can ha 353 // hadron. Instead, any projectile anti-bary 354 // with a target hadron (because it has only 355 // projectile baryons, instead can have char 356 357 G4int NewProjCode = 0, NewTargCode = 0, Proj 358 // Projectile unpacking 359 if ( common.absProjectilePDGcode < 1000 ) { 360 UnpackMeson( common.ProjectilePDGcode, Pr 361 } else { 362 UnpackBaryon( common.ProjectilePDGcode, Pr 363 } 364 // Target unpacking 365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0; 366 UnpackBaryon( common.TargetPDGcode, TargQ1, 367 #ifdef debugFTFexcitation 368 G4cout << "Proj Quarks " << ProjQ1 << " " << 369 << "Targ Quarks " << TargQ1 << " " << 370 #endif 371 372 // Sampling of exchanged quarks 373 G4int ProjExchangeQ = 0, TargExchangeQ = 0; 374 if ( common.absProjectilePDGcode < 1000 ) { 375 376 G4bool isProjQ1Quark = false; 377 ProjExchangeQ = ProjQ2; 378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark 379 isProjQ1Quark = true; 380 ProjExchangeQ = ProjQ1; 381 } 382 // Exchange of non-identical quarks is all 383 G4int NpossibleStates = 0; 384 if ( ProjExchangeQ != TargQ1 ) NpossibleSt 385 if ( ProjExchangeQ != TargQ2 ) NpossibleSt 386 if ( ProjExchangeQ != TargQ3 ) NpossibleSt 387 G4int Nsampled = (G4int)G4RandFlat::shootI 388 NpossibleStates = 0; 389 if ( ProjExchangeQ != TargQ1 ) { 390 if ( ++NpossibleStates == Nsampled ) { 391 TargExchangeQ = TargQ1; TargQ1 = ProjE 392 isProjQ1Quark ? ProjQ1 = TargExchangeQ 393 } 394 } 395 if ( ProjExchangeQ != TargQ2 ) { 396 if ( ++NpossibleStates == Nsampled ) { 397 TargExchangeQ = TargQ2; TargQ2 = ProjE 398 isProjQ1Quark ? ProjQ1 = TargExchangeQ 399 } 400 } 401 if ( ProjExchangeQ != TargQ3 ) { 402 if ( ++NpossibleStates == Nsampled ) { 403 TargExchangeQ = TargQ3; TargQ3 = ProjE 404 isProjQ1Quark ? ProjQ1 = TargExchangeQ 405 } 406 } 407 #ifdef debugFTFexcitation 408 G4cout << "Exchanged Qs in Pr Tr " << Proj 409 #endif 410 411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ 412 G4bool ProjExcited = false; 413 const G4int maxNumberOfAttempts = 50; 414 G4int attempts = 0; 415 while ( attempts++ < maxNumberOfAttempts ) 416 417 // Determination of a new projectile ID 418 G4double ProbSpin0 = 0.5; 419 G4double Ksi = G4UniformRand(); 420 if ( aProjQ1 == aProjQ2 ) { 421 if ( G4UniformRand() < ProbSpin0 ) { 422 if ( aProjQ1 < 3 ) { 423 NewProjCode = 111; 424 if ( Ksi < 0.5 ) { 425 NewProjCode = 221; 426 if ( Ksi < 0.25 ) { 427 NewProjCode = 331; 428 } 429 } 430 } else if ( aProjQ1 == 3 ) { 431 NewProjCode = 221; 432 if ( Ksi < 0.5 ) { 433 NewProjCode = 331; 434 } 435 } else if ( aProjQ1 == 4 ) { 436 NewProjCode = 441; // eta 437 } else if ( aProjQ1 == 5 ) { 438 NewProjCode = 551; // eta 439 } 440 } else { 441 if ( aProjQ1 < 3 ) { 442 NewProjCode = 113; 443 if ( Ksi < 0.5 ) { 444 NewProjCode = 223; 445 } 446 } else if ( aProjQ1 == 3 ) { 447 NewProjCode = 333; 448 } else if ( aProjQ1 == 4 ) { 449 NewProjCode = 443; // J/p 450 } else if ( aProjQ1 == 5 ) { 451 NewProjCode = 553; // Ups 452 } 453 } 454 } else { 455 if ( aProjQ1 > aProjQ2 ) { 456 NewProjCode = aProjQ1*100 + aProjQ2* 457 } else { 458 NewProjCode = aProjQ2*100 + aProjQ1* 459 } 460 } 461 #ifdef debugFTFexcitation 462 G4cout << "NewProjCode " << NewProjCode 463 #endif 464 // Decide (with 50% probability) whether 465 // but not in the case of charmed and bo 466 // there are no excited charmed and bott 467 ProjExcited = false; 468 if ( aProjQ1 <= 3 && aProjQ2 <= 3 && 469 NewProjCode += 2; // Excited meson (J 470 ProjExcited = true; 471 } 472 473 // So far we have used the absolute valu 474 // now look at their signed values to se 475 G4int value = ProjQ1, absValue = aProjQ1 476 for ( G4int iQ = 0; iQ < 2; ++iQ ) { // 477 if ( iQ == 1 ) { 478 value = ProjQ2; absValue = aProjQ2; 479 } 480 if ( absValue == 2 || absValue == 4 481 Qquarks += 2*value/absValue; // u, 482 } else { 483 Qquarks -= value/absValue; // d, 484 } 485 } 486 // If Qquarks is positive, the sign of N 487 // If Qquarks is negative, then the sign 488 // If Qquarks is zero, then we need to d 489 // 1. If aProjQ1 and aProjQ2 are the sam 490 // (because the antiparticle is the s 491 // 2. If aProjQ1 and aProjQ2 are not the 492 // we have only two possibilities: 493 // a. aProjQ1 and aProjQ2 are two dif 494 // (s,d) or (b,d), or (b,s). In th 495 // is fine (because the heaviest o 496 // to be anti-quark belonging to t 497 // implies a meson with positive P 498 // b. aProjQ1 and aProjQ2 are two dif 499 // The heaviest of the two (c) has 500 // in the projectile, therefore th 501 // reverse: 421 -> -421 anti_D0 (c 502 if ( Qquarks < 0 || ( Qquarks == 0 && 503 NewProjCode *= -1; 504 } 505 #ifdef debugFTFexcitation 506 G4cout << "NewProjCode +2 or 0 " << NewP 507 G4cout<<"+++++++++++++++++++++++++++++++ 508 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquark 509 G4cout<<"+++++++++++++++++++++++++++++++ 510 #endif 511 512 // Proj 513 TestParticle = G4ParticleTable::GetParti 514 if ( ! TestParticle ) continue; 515 common.MminProjectile = common.BrW.GetMi 516 if ( common.SqrtS - common.M0target < co 517 MtestPr = common.BrW.SampleMass( TestPar 518 519 #ifdef debugFTFexcitation 520 G4cout << "TestParticle Name " << NewPro 521 << G4endl 522 << "MtestPart MtestPart0 "<<Mtest 523 << "M0projectile projectile PDGMa 524 << projectile->GetDefinition()->G 525 #endif 526 527 // Targ 528 NewTargCode = NewNucleonId( TargQ1, Targ 529 #ifdef debugFTFexcitation 530 G4cout << "New TrQ " << TargQ1 << " " << 531 << "NewTargCode " << NewTargCode 532 #endif 533 534 // If the target has not heavy (charm or 535 // see whether a Delta isobar can be cre 536 if ( TargQ1 <= 3 && TargQ2 <= 3 && T 537 if ( TargQ1 != TargQ2 && TargQ1 != T 538 if ( G4UniformRand() < 0.5 ) { 539 NewTargCode += 2; 540 } else if ( G4UniformRand() < 0.75 ) 541 NewTargCode = 3122; // Lambda 542 } 543 } else if ( TargQ1 == TargQ2 && Targ 544 NewTargCode += 2; ProjExcited = true 545 } else if ( target->GetDefinition()->G 546 if ( G4UniformRand() > DeltaProbAtQu 547 NewTargCode += 2; ProjExcited = tr 548 } 549 } else if ( ! ProjExcited && 550 G4UniformRand() < DeltaPro 551 common.SqrtS > common.M0pr 552 G4ParticleTable::GetP 553 NewTargCode += 2; // Create Delta i 554 } 555 } 556 557 // Protection against: 558 // - Lambda* (i.e. excited Lambda state) 559 // - Sigma* (i.e. excited Sigma hyperon 560 // - Xi* (i.e. excited Xi hyperon st 561 if ( NewTargCode == 3124 || // Lambda 562 NewTargCode == 3224 || // Sigma*+ NOT 563 NewTargCode == 3214 || // Sigma*0 NOT 564 NewTargCode == 3114 || // Sigma*- NOT 565 NewTargCode == 3324 || // Xi*0 NOT 566 NewTargCode == 3314 ) { // Xi*- NOT 567 //G4cout << "G4DiffractiveExcitation::Excite 568 // << NewTargCode << G4endl; 569 NewTargCode -= 2; // Corresponding ground-s 570 } 571 572 // Special treatment for charmed and bot 573 // so we need to transform them by hand 574 #ifdef debug_heavyHadrons 575 G4int initialNewTargCode = NewTargCode; 576 #endif 577 if ( NewTargCode == 4322 ) NewTargC 578 else if ( NewTargCode == 4312 ) NewTargC 579 else if ( NewTargCode == 5312 ) NewTargC 580 else if ( NewTargCode == 5322 ) NewTargC 581 #ifdef debug_heavyHadrons 582 if ( NewTargCode != initialNewTargCode ) 583 G4cout << "G4DiffractiveExcitation::ExcitePa 584 << "\t target heavy baryon with pdgCo 585 << " into pdgCode=" << NewTargCode << 586 } 587 #endif 588 589 TestParticle = G4ParticleTable::GetParti 590 if ( ! TestParticle ) continue; 591 #ifdef debugFTFexcitation 592 G4cout << "New targ " << NewTargCode << 593 #endif 594 common.MminTarget = common.BrW.GetMinimu 595 if ( common.SqrtS - MtestPr < common.Mmi 596 MtestTr = common.BrW.SampleMass( TestPar 597 598 if ( common.SqrtS > MtestPr + MtestTr ) 599 600 } // End of while loop 601 if ( attempts >= maxNumberOfAttempts ) ret 602 603 if ( MtestPr >= common.Pprojectile.mag() 604 common.M0projectile = MtestPr; 605 } 606 #ifdef debugFTFexcitation 607 G4cout << "M0projectile After check " << c 608 #endif 609 common.M0projectile2 = common.M0projectile 610 common.ProjectileDiffStateMinMass = com 611 common.ProjectileNonDiffStateMinMass = com 612 if ( MtestTr >= common.Ptarget.mag() || 613 common.M0target = MtestTr; 614 } 615 common.M0target2 = common.M0target * commo 616 #ifdef debugFTFexcitation 617 G4cout << "New targ M0 M0^2 " << common.M0 618 #endif 619 common.TargetDiffStateMinMass = common. 620 common.TargetNonDiffStateMinMass = common. 621 622 } else { // of the if ( common.absProjectil 623 624 // If it is a projectile anti-baryon, no q 625 // therefore returns immediately 1 (which 626 // needs to be continued"). 627 if ( common.ProjectilePDGcode < 0 ) return 628 629 // Choose randomly, with equal probability 630 // projectile or target hadron for selecti 631 G4bool isProjectileExchangedQ = false; 632 G4int firstQ = TargQ1, secondQ = 633 G4int otherFirstQ = ProjQ1, otherSecondQ = 634 if ( G4UniformRand() < 0.5 ) { 635 isProjectileExchangedQ = true; 636 firstQ = ProjQ1; secondQ = Pro 637 otherFirstQ = TargQ1; otherSecondQ = Tar 638 } 639 // Choose randomly, with equal probability 640 // selected (projectile or target) hadron 641 G4int exchangedQ = 0; 642 G4double Ksi = G4UniformRand(); 643 if ( Ksi < 0.333333 ) { 644 exchangedQ = firstQ; 645 } else if ( 0.333333 <= Ksi && Ksi < 0.6 646 exchangedQ = secondQ; 647 } else { 648 exchangedQ = thirdQ; 649 } 650 #ifdef debugFTFexcitation 651 G4cout << "Exchange Qs isProjectile Q " << 652 #endif 653 // The exchanged quarks (one of the projec 654 // are always accepted if they have differ 655 // same flavour) they are accepted only wi 656 G4double probSame = theParameters->GetProb 657 const G4int MaxCount = 100; 658 G4int count = 0, otherExchangedQ = 0; 659 do { 660 if ( exchangedQ != otherFirstQ || G4Un 661 otherExchangedQ = otherFirstQ; otherFi 662 } else { 663 if ( exchangedQ != otherSecondQ || G 664 otherExchangedQ = otherSecondQ; othe 665 } else { 666 if ( exchangedQ != otherThirdQ || 667 otherExchangedQ = otherThirdQ; oth 668 } 669 } 670 } 671 } while ( otherExchangedQ == 0 && ++coun 672 if ( count >= MaxCount ) return returnCode 673 // Swap (between projectile and target had 674 // as "exchanged" quarks. 675 if ( Ksi < 0.333333 ) { 676 firstQ = exchangedQ; 677 } else if ( 0.333333 <= Ksi && Ksi < 0.6 678 secondQ = exchangedQ; 679 } else { 680 thirdQ = exchangedQ; 681 } 682 if ( isProjectileExchangedQ ) { 683 ProjQ1 = firstQ; ProjQ2 = secondQ; 684 TargQ1 = otherFirstQ; TargQ2 = otherSeco 685 } else { 686 TargQ1 = firstQ; TargQ2 = secondQ; 687 ProjQ1 = otherFirstQ; ProjQ2 = otherSeco 688 } 689 #ifdef debugFTFexcitation 690 G4cout << "Exchange Qs Pr Tr " << ( isPro 691 << " " << ( isProjectileExchangedQ 692 #endif 693 694 NewProjCode = NewNucleonId( ProjQ1, ProjQ2 695 NewTargCode = NewNucleonId( TargQ1, TargQ2 696 // Decide whether the new projectile hadro 697 // then decide whether the new target hadr 698 // Notice that a Delta particle has the la 699 // whereas a nucleon has "2" (because its 700 // If the new projectile hadron or the new 701 // constituent quark, then skip this part 702 // excited charm and bottom hadrons). 703 for ( G4int iHadron = 0; iHadron < 2; iHad 704 // First projectile hadron, then target 705 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, 706 G4double massConstraint = common.M0targe 707 G4bool isHadronADelta = ( projectile->Ge 708 if ( iHadron == 1 ) { // Target hadron 709 codeQ1 = TargQ1, codeQ2 = TargQ2, code 710 massConstraint = common.M0projectile; 711 isHadronADelta = ( target->GetDefiniti 712 } 713 if ( codeQ1 > 3 || codeQ2 > 3 || cod 714 if ( codeQ1 == codeQ2 && codeQ1 == cod 715 newHadCode += 2; // Delta++ (uuu) or 716 } else if ( isHadronADelta ) { // Hadro 717 if ( G4UniformRand() > DeltaProbAtQuar 718 newHadCode += 2; // Delta+ (uud) or 719 } else { 720 newHadCode += 0; // No delta (so th 721 } 722 } else { // Hadron (projectile or targe 723 if ( G4UniformRand() < DeltaProbAtQuar 724 common.SqrtS > G4ParticleTable::G 725 + massConstraint ) 726 newHadCode += 2; // Delta+ (uud) or 727 } else { 728 newHadCode += 0; // No delta (so th 729 } 730 } 731 if ( iHadron == 0 ) { // Projectile had 732 NewProjCode = newHadCode; 733 } else { // Target hadron 734 NewTargCode = newHadCode; 735 } 736 } 737 #ifdef debugFTFexcitation 738 G4cout << "NewProjCode NewTargCode " << Ne 739 #endif 740 741 // Protection against: 742 // - Lambda* (i.e. excited Lambda state) 743 // - Sigma* (i.e. excited Sigma hyperon s 744 // - Xi* (i.e. excited Xi hyperon stat 745 if ( NewProjCode == 3124 || // Lambda* 746 NewProjCode == 3224 || // Sigma*+ NOT ex 747 NewProjCode == 3214 || // Sigma*0 NOT ex 748 NewProjCode == 3114 || // Sigma*- NOT ex 749 NewProjCode == 3324 || // Xi*0 NOT ex 750 NewProjCode == 3314 ) { // Xi*- NOT ex 751 //G4cout << "G4DiffractiveExcitation::Ex 752 // << NewProjCode << G4endl; 753 NewProjCode -= 2; // Corresponding grou 754 } 755 if ( NewTargCode == 3124 || // Lambda* 756 NewTargCode == 3224 || // Sigma*+ NOT ex 757 NewTargCode == 3214 || // Sigma*0 NOT ex 758 NewTargCode == 3114 || // Sigma*- NOT ex 759 NewTargCode == 3324 || // Xi*0 NOT ex 760 NewTargCode == 3314 ) { // Xi*- NOT ex 761 //G4cout << "G4DiffractiveExcitation::Ex 762 // << NewTargCode << G4endl; 763 NewTargCode -= 2; // Corresponding grou 764 } 765 766 // Special treatment for charmed and botto 767 // so we need to transform them by hand to 768 #ifdef debug_heavyHadrons 769 G4int initialNewProjCode = NewProjCode, in 770 #endif 771 if ( NewProjCode == 4322 ) NewProjCod 772 else if ( NewProjCode == 4312 ) NewProjCod 773 else if ( NewProjCode == 5312 ) NewProjCod 774 else if ( NewProjCode == 5322 ) NewProjCod 775 if ( NewTargCode == 4322 ) NewTargCod 776 else if ( NewTargCode == 4312 ) NewTargCod 777 else if ( NewTargCode == 5312 ) NewTargCod 778 else if ( NewTargCode == 5322 ) NewTargCod 779 #ifdef debug_heavyHadrons 780 if ( NewProjCode != initialNewProjCode || 781 G4cout << "G4DiffractiveExcitation::Exci 782 << "\t heavy baryon into an exist 783 if ( NewProjCode != initialNewProjCode ) 784 G4cout << "\t \t projectile baryon with pdgC 785 << " into pdgCode=" << NewProjCode << 786 } 787 if ( NewTargCode != initialNewTargCode ) 788 G4cout << "\t \t target baryon with pd 789 << " into pdgCode=" << NewTargCode << 790 } 791 } 792 #endif 793 794 // Sampling of the masses of the projectil 795 // Because of energy conservation, the ord 796 // randomly, half of the time we sample fi 797 // then the projectile nucleon mass, and t 798 // sample first the projectile nucleon mas 799 G4VSplitableHadron* firstHadron = target; 800 G4VSplitableHadron* secondHadron = project 801 G4int firstHadronCode = NewTargCode, secon 802 G4double massConstraint = common.M0project 803 G4bool isFirstTarget = true; 804 if ( G4UniformRand() < 0.5 ) { 805 // Sample first the projectile nucleon m 806 firstHadron = projectile; secondHa 807 firstHadronCode = NewProjCode; secondHa 808 massConstraint = common.M0target; 809 isFirstTarget = false; 810 } 811 G4double Mtest1st = 0.0, Mtest2nd = 0.0, M 812 for ( int iSamplingCase = 0; iSamplingCase 813 G4VSplitableHadron* aHadron = firstHadro 814 G4int aHadronCode = firstHadronCode; 815 if ( iSamplingCase == 1 ) { // Second n 816 aHadron = secondHadron; aHadronCode = 817 } 818 G4double MtestHadron = 0.0, MminHadron = 819 if ( aHadron->GetStatus() == 1 || aHad 820 TestParticle = G4ParticleTable::GetPar 821 if ( ! TestParticle ) return returnCod 822 MminHadron = common.BrW.GetMinimumMass 823 if ( common.SqrtS - massConstraint < M 824 if ( TestParticle->GetPDGWidth() == 0. 825 MtestHadron = common.BrW.SampleMass( 826 } else { 827 const G4int maxNumberOfAttempts = 50 828 G4int attempts = 0; 829 while ( attempts < maxNumberOfAttemp 830 attempts++; 831 MtestHadron = common.BrW.SampleMas 832 833 if ( common.SqrtS < MtestHadron + 834 continue; // Kinematically unac 835 } else { 836 break; // Kinematically acce 837 } 838 } 839 if ( attempts >= maxNumberOfAttempts 840 } 841 } 842 if ( iSamplingCase == 0 ) { 843 Mtest1st = MtestHadron; Mmin1st = Mmi 844 } else { 845 Mtest2nd = MtestHadron; Mmin2nd = Mmi 846 } 847 } // End for loop on the two sampling cas 848 if ( isFirstTarget ) { 849 MtestTr = Mtest1st; MtestPr = Mtest2n 850 common.MminTarget = Mmin1st; common.Mmi 851 } else { 852 MtestTr = Mtest2nd; MtestPr = Mtest1s 853 common.MminTarget = Mmin2nd; common.Mmi 854 } 855 856 if ( MtestPr != 0.0 ) { 857 common.M0projectile = MtestPr; 858 common.M0projectile2 = common.M0projecti 859 common.ProjectileDiffStateMinMass = c 860 common.ProjectileNonDiffStateMinMass = c 861 } 862 if ( MtestTr != 0.0 ) { 863 common.M0target = MtestTr; 864 common.M0target2 = common.M0target * com 865 common.TargetDiffStateMinMass = commo 866 common.TargetNonDiffStateMinMass = commo 867 } 868 869 } // End of if ( common.absProjectilePDGcod 870 871 // If we assume that final state hadrons aft 872 // in the ground states, we have to put 873 if ( common.SqrtS < common.M0projectile + co 874 875 common.PZcms2 = ( sqr( common.S ) + sqr( com 876 - 2.0 * ( common.S * ( com 877 + common.M0proje 878 #ifdef debugFTFexcitation 879 G4cout << "At the end// NewProjCode " << New 880 << "At the end// NewTargCode " << New 881 << "M0pr M0tr SqS " << common.M0pro 882 << common.SqrtS << G4endl 883 << "M0pr2 M0tr2 SqS " << common.M0pro 884 << common.SqrtS << G4endl 885 << "PZcms2 after the change " << comm 886 #endif 887 if ( common.PZcms2 < 0.0 ) return returnCode 888 889 projectile->SetDefinition( G4ParticleTable:: 890 target->SetDefinition( G4ParticleTable::GetP 891 common.PZcms = std::sqrt( common.PZcms2 ); 892 common.Pprojectile.setPz( common.PZcms ); 893 common.Pprojectile.setE( std::sqrt( common.M 894 common.Ptarget.setPz( -common.PZcms ); 895 common.Ptarget.setE( std::sqrt( common.M0tar 896 #ifdef debugFTFexcitation 897 G4cout << "Proj Targ and Proj+Targ in CMS" < 898 << common.Ptarget << G4endl << common 899 #endif 900 901 if ( projectile->GetStatus() != 0 ) projecti 902 if ( target->GetStatus() != 0 ) target-> 903 904 // Check for possible excitation of the part 905 if ( common.SqrtS < common.M0projectile + co 906 common.SqrtS < common.ProjectileDiffSta 907 common.ProbOfDiffraction == 0.0 ) commo 908 909 if ( G4UniformRand() > common.ProbExc ) { / 910 #ifdef debugFTFexcitation 911 G4cout << "Make elastic scattering of new 912 #endif 913 common.Pprojectile.transform( common.toLab 914 common.Ptarget.transform( common.toLab ); 915 projectile->Set4Momentum( common.Pprojecti 916 target->Set4Momentum( common.Ptarget ); 917 G4bool Result = theElastic->ElasticScatter 918 #ifdef debugFTFexcitation 919 G4cout << "Result of el. scatt " << Result 920 << G4endl << projectile->Get4Moment 921 << projectile->Get4Momentum() + tar 922 << (projectile->Get4Momentum() + ta 923 #endif 924 if ( Result ) returnCode = 0; // successf 925 return returnCode; 926 } 927 928 #ifdef debugFTFexcitation 929 G4cout << "Make excitation of new hadrons" < 930 #endif 931 932 // Redefinition of ProbOfDiffraction because 933 common.ProbOfDiffraction = common.ProbProjec 934 if ( common.ProbOfDiffraction != 0.0 ) { 935 common.ProbProjectileDiffraction /= common 936 common.ProbTargetDiffraction /= common 937 } 938 939 return returnCode = 1; // successfully comp 940 } 941 942 //-------------------------------------------- 943 944 G4bool G4DiffractiveExcitation:: 945 ExciteParticipants_doDiffraction( G4VSplitable 946 G4FTFParamet 947 G4Diffractiv 948 // Second of the three utility methods used 949 // it does the sampling for the diffraction 950 951 G4bool isProjectileDiffraction = false; 952 if ( G4UniformRand() < common.ProbProjectile 953 isProjectileDiffraction = true; 954 #ifdef debugFTFexcitation 955 G4cout << "projectile diffraction" << G4en 956 #endif 957 common.ProjMassT2 = common.ProjectileDiffS 958 common.ProjMassT = common.ProjectileDiffS 959 common.TargMassT2 = common.M0target2; 960 common.TargMassT = common.M0target; 961 } else { 962 #ifdef debugFTFexcitation 963 G4cout << "Target diffraction" << G4endl; 964 #endif 965 common.ProjMassT2 = common.M0projectile2; 966 common.ProjMassT = common.M0projectile; 967 common.TargMassT2 = common.TargetDiffState 968 common.TargMassT = common.TargetDiffState 969 } 970 971 // Check whether the interaction is possible 972 if ( common.SqrtS < common.ProjMassT + commo 973 974 common.PZcms2 = ( sqr( common.S ) + sqr( com 975 - 2.0 * ( common.S * ( com 976 + common.ProjMassT2 977 if ( common.PZcms2 < 0.0 ) return false; 978 common.maxPtSquare = common.PZcms2; 979 980 G4double DiffrAveragePt2 = theParameters->Ge 981 G4bool loopCondition = true; 982 G4int whilecount = 0; 983 do { // Generate pt and mass of projectile 984 985 whilecount++; 986 if ( whilecount > 1000 ) { 987 common.Qmomentum = G4LorentzVector( 0.0, 988 return false; // Ignore this interacti 989 }; 990 991 common.Qmomentum = G4LorentzVector( Gaussi 992 common.Pt2 = G4ThreeVector( common.Qmoment 993 if ( isProjectileDiffraction ) { // proje 994 common.ProjMassT2 = common.ProjectileDif 995 common.TargMassT2 = common.M0target2 + c 996 } else { // targe 997 common.ProjMassT2 = common.M0projectile2 998 common.TargMassT2 = common.TargetDiffSta 999 } 1000 common.ProjMassT = std::sqrt( common.Proj 1001 common.TargMassT = std::sqrt( common.Targ 1002 if ( common.SqrtS < common.ProjMassT + co 1003 1004 common.PZcms2 = ( sqr( common.S ) + sqr( 1005 - 2.0 * ( common.S * ( 1006 + common.Proj 1007 if ( common.PZcms2 < 0.0 ) continue; 1008 1009 common.PZcms = std::sqrt( common.PZcms2 ) 1010 if ( isProjectileDiffraction ) { // proj 1011 common.PMinusMin = std::sqrt( common.Pr 1012 common.PMinusMax = common.SqrtS - commo 1013 common.PMinusNew = ChooseP( common.PMin 1014 common.TMinusNew = common.SqrtS - commo 1015 common.Qminus = common.Ptarget.minus() 1016 common.TPlusNew = common.TargMassT2 / c 1017 common.Qplus = common.Ptarget.plus() - 1018 common.Qmomentum.setPz( (common.Qplus - 1019 common.Qmomentum.setE( (common.Qplus + 1020 loopCondition = ( ( common.Pprojectile 1021 common.ProjectileDiff 1022 } else { // targ 1023 common.TPlusMin = std::sqrt( common.Tar 1024 common.TPlusMax = common.SqrtS - common 1025 common.TPlusNew = ChooseP( common.TPlus 1026 common.PPlusNew = common.SqrtS - common 1027 common.Qplus = common.PPlusNew - common 1028 common.PMinusNew = common.ProjMassT2 / 1029 common.Qminus = common.PMinusNew - comm 1030 common.Qmomentum.setPz( (common.Qplus - 1031 common.Qmomentum.setE( (common.Qplus + 1032 loopCondition = ( ( common.Ptarget - c 1033 common.TargetDiffSta 1034 } 1035 1036 } while ( loopCondition ); /* Loop checkin 1037 // Repeat the sampling because ther 1038 1039 if ( isProjectileDiffraction ) { // projec 1040 projectile->SetStatus( 0 ); 1041 if ( projectile->GetStatus() == 2 ) proje 1042 if ( target->GetStatus() == 1 && target 1043 } else { // target 1044 target->SetStatus( 0 ); 1045 } 1046 1047 return true; 1048 } 1049 1050 //------------------------------------------- 1051 1052 G4bool G4DiffractiveExcitation:: 1053 ExciteParticipants_doNonDiffraction( G4VSplit 1054 G4VSplit 1055 G4FTFPar 1056 G4Diffra 1057 // Third of the three utility methods used 1058 // it does the sampling for the non-diffrac 1059 1060 #ifdef debugFTFexcitation 1061 G4cout << "Non-diffraction process" << G4en 1062 #endif 1063 1064 // Check whether the interaction is possibl 1065 common.ProjMassT2 = common.ProjectileNonDif 1066 common.ProjMassT = common.ProjectileNonDif 1067 common.TargMassT2 = common.TargetNonDiffSta 1068 common.TargMassT = common.TargetNonDiffSta 1069 if ( common.SqrtS < common.ProjMassT + comm 1070 1071 common.PZcms2 = ( sqr( common.S ) + sqr( co 1072 - 2.0 * ( common.S * ( comm 1073 + common.ProjMassT2 * 1074 common.maxPtSquare = common.PZcms2; 1075 1076 G4int whilecount = 0; 1077 do { // Generate pt and masses 1078 1079 whilecount++; 1080 if ( whilecount > 1000 ) { 1081 common.Qmomentum = G4LorentzVector( 0.0 1082 return false; // Ignore this interacti 1083 }; 1084 1085 common.Qmomentum = G4LorentzVector( Gauss 1086 1087 common.Pt2 = G4ThreeVector( common.Qmomen 1088 common.ProjMassT2 = common.ProjectileNonD 1089 common.ProjMassT = std::sqrt( common.Pro 1090 common.TargMassT2 = common.TargetNonDiffS 1091 common.TargMassT = std::sqrt( common.Tar 1092 if ( common.SqrtS < common.ProjMassT + co 1093 1094 common.PZcms2 =( sqr( common.S ) + sqr( c 1095 - 2.0 * ( common.S * ( c 1096 + common.ProjM 1097 if ( common.PZcms2 < 0.0 ) continue; 1098 1099 common.PZcms = std::sqrt( common.PZcms2 ) 1100 common.PMinusMin = std::sqrt( common.Proj 1101 common.PMinusMax = common.SqrtS - common. 1102 common.TPlusMin = std::sqrt( common.TargM 1103 common.TPlusMax = common.SqrtS - common.P 1104 1105 if ( G4UniformRand() <= 0.5 ) { // Rando 1106 if ( G4UniformRand() < theParameters->G 1107 common.PMinusNew = ChooseP( common.PM 1108 } else { 1109 common.PMinusNew = ( common.PMinusMax 1110 } 1111 if ( G4UniformRand() < theParameters->G 1112 common.TPlusNew = ChooseP( common.TPl 1113 } else { 1114 common.TPlusNew = ( common.TPlusMax - 1115 } 1116 } else { 1117 if ( G4UniformRand() < theParameters->G 1118 common.TPlusNew = ChooseP( common.TPl 1119 } else { 1120 common.TPlusNew = ( common.TPlusMax - 1121 } 1122 if ( G4UniformRand() < theParameters->G 1123 common.PMinusNew = ChooseP( common.PM 1124 } else { 1125 common.PMinusNew = ( common.PMinusMax 1126 } 1127 } 1128 1129 common.Qminus = common.PMinusNew - common 1130 common.Qplus = -( common.TPlusNew - commo 1131 common.Qmomentum.setPz( (common.Qplus - c 1132 common.Qmomentum.setE( (common.Qplus + c 1133 #ifdef debugFTFexcitation 1134 G4cout <<"Sampled: Mpr, MdifPr, Mtr, Mdif 1135 << ( common.Pprojectile + common.Q 1136 << common.ProjectileNonDiffStateMi 1137 << ( common.Ptarget - common.Qmome 1138 << common.TargetNonDiffStateMinMas 1139 #endif 1140 1141 } while ( ( common.Pprojectile + common.Qmo 1142 common.ProjectileNonDiffStateMinM 1143 ( common.Ptarget - common.Qmoment 1144 common.TargetNonDiffStateMinMass2 1145 1146 projectile->SetStatus( 0 ); 1147 target->SetStatus( 0 ); 1148 return true; 1149 } 1150 1151 1152 //=========================================== 1153 1154 void G4DiffractiveExcitation::CreateStrings( 1155 1156 1157 1158 1159 1160 //G4cout << "Create Strings SplitUp " << ha 1161 // << "Defin " << hadron->GetDefiniti 1162 // << "Defin " << hadron->GetDefiniti 1163 1164 G4bool HadronIsString = hadron->IsSplit(); 1165 if( ! HadronIsString ) hadron->SplitUp(); 1166 1167 G4Parton* start = hadron->GetNextParton(); 1168 if ( start == nullptr ) { 1169 G4cout << " G4FTFModel::String() Error: N 1170 FirstString = 0; SecondString = 0; 1171 return; 1172 } 1173 1174 G4Parton* end = hadron->GetNextParton(); 1175 if ( end == nullptr ) { 1176 G4cout << " G4FTFModel::String() Error: N 1177 FirstString = 0; SecondString = 0; 1178 return; 1179 } 1180 1181 //G4cout << start << " " << start->GetPDGco 1182 // << G4endl 1183 // << "Create string " << start->GetP 1184 1185 if ( HadronIsString ) { 1186 if ( isProjectile ) { 1187 FirstString = new G4ExcitedString( end, 1188 } else { 1189 FirstString = new G4ExcitedString( end, 1190 } 1191 FirstString->SetTimeOfCreation( hadron->G 1192 FirstString->SetPosition( hadron->GetPosi 1193 SecondString = 0; 1194 return; 1195 } 1196 1197 G4LorentzVector Phadron = hadron->Get4Momen 1198 //G4cout << "String mom " << Phadron << G4e 1199 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 1200 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 ); 1201 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 ) 1202 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 1203 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 1204 1205 G4int PDGcode_startQ = std::abs( start->Get 1206 G4int PDGcode_endQ = std::abs( end->Get 1207 //G4cout << "PDGcode_startQ " << PDGcode_st 1208 1209 G4double Wmin( 0.0 ); 1210 if ( isProjectile ) { 1211 Wmin = theParameters->GetProjMinDiffMass( 1212 } else { 1213 Wmin = theParameters->GetTarMinDiffMass() 1214 } 1215 1216 G4double W = hadron->Get4Momentum().mag(); 1217 G4double W2 = W*W; 1218 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); 1219 G4bool Kink = false; 1220 1221 if ( ! ( ( start->GetDefinition()->GetParti 1222 end->GetDefinition()->GetParti 1223 ( start->GetDefinition()->GetParti 1224 end->GetDefinition()->GetParti 1225 // Kinky strings are allowed only for qq- 1226 // Kinky strings are impossible for other 1227 // according to the analysis of Pbar P in 1228 1229 if ( W > Wmin ) { // Kink is possible 1230 if ( hadron->GetStatus() == 0 ) { 1231 G4double Pt2kink = theParameters->Get 1232 if ( Pt2kink ) { 1233 Pt = std::sqrt( Pt2kink * ( G4Pow:: 1234 } else { 1235 Pt = 0.0; 1236 } 1237 } else { 1238 Pt = 0.0; 1239 } 1240 1241 if ( Pt > 500.0*MeV ) { 1242 G4double Ymax = G4Log( W/2.0/Pt + std 1243 G4double Y = Ymax*( 1.0 - 2.0*G4Unifo 1244 x1 = 1.0 - Pt/W * G4Exp( Y ); 1245 x3 = 1.0 - Pt/W * G4Exp(-Y ); 1246 //x2 = 2.0 - x1 - x3; 1247 1248 G4double Mass_startQ = 650.0*MeV; 1249 if ( PDGcode_startQ < 3 ) Mass_start 1250 if ( PDGcode_startQ == 3 ) Mass_start 1251 if ( PDGcode_startQ == 4 ) Mass_start 1252 if ( PDGcode_startQ == 5 ) Mass_start 1253 G4double Mass_endQ = 650.0*MeV; 1254 if ( PDGcode_endQ < 3 ) Mass_endQ = 1255 if ( PDGcode_endQ == 3 ) Mass_endQ = 1256 if ( PDGcode_endQ == 4 ) Mass_endQ = 1257 if ( PDGcode_endQ == 5 ) Mass_endQ = 1258 1259 G4double P2_1 = W2*x1*x1/4.0 - Mass_e 1260 G4double P2_3 = W2*x3*x3/4.0 - Mass_s 1261 G4double P2_2 = sqr( (2.0 - x1 - x3)* 1262 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) { 1263 Kink = false; 1264 } else { 1265 G4double P_1 = std::sqrt( P2_1 ); 1266 G4double P_2 = std::sqrt( P2_2 ); 1267 G4double P_3 = std::sqrt( P2_3 ); 1268 G4double CosT12 = ( P2_3 - P2_1 - P 1269 G4double CosT13 = ( P2_2 - P2_1 - P 1270 1271 if ( std::abs( CosT12 ) > 1.0 || 1272 Kink = false; 1273 } else { 1274 Kink = true; 1275 Pt = P_2 * std::sqrt( 1.0 - CosT1 1276 Pstart.setPx( -Pt ); Pstart.setPy 1277 Pend.setPx( 0.0 ); Pend.setPy( 1278 Pkink.setPx( Pt ); Pkink.setPy( 1279 Pstart.setE( x3*W/2.0 ); 1280 Pkink.setE( Pkink.vect().mag() ); 1281 Pend.setE( x1*W/2.0 ); 1282 1283 G4double XkQ = GetQuarkFractionOf 1284 if ( Pkink.getZ() > 0.0 ) { 1285 if ( XkQ > 0.5 ) { 1286 PkinkQ1 = XkQ*Pkink; 1287 } else { 1288 PkinkQ1 = (1.0 - XkQ)*Pkink; 1289 } 1290 } else { 1291 if ( XkQ > 0.5 ) { 1292 PkinkQ1 = (1.0 - XkQ)*Pkink; 1293 } else { 1294 PkinkQ1 = XkQ*Pkink; 1295 } 1296 } 1297 1298 PkinkQ2 = Pkink - PkinkQ1; 1299 // Minimizing Pt1^2+Pt3^2 1300 G4double Cos2Psi = ( sqr(x1) - sq 1301 std::sqrt( sqr 1302 G4double Psi = std::acos( Cos2Psi 1303 1304 G4LorentzRotation Rotate; 1305 if ( isProjectile ) { 1306 Rotate.rotateY( Psi ); 1307 } else { 1308 Rotate.rotateY( pi + Psi ); 1309 } 1310 Rotate.rotateZ( twopi * G4Uniform 1311 Pstart *= Rotate; 1312 Pkink *= Rotate; 1313 PkinkQ1 *= Rotate; 1314 PkinkQ2 *= Rotate; 1315 Pend *= Rotate; 1316 } 1317 } // End of if ( P2_1 <= 0.0 || P2 1318 } // End of if ( Pt > 500.0*MeV ) 1319 } // End of if ( W > Wmin ) : check for a 1320 } // end of qq-q string selection 1321 1322 if ( Kink ) { // Kink is possible 1323 1324 //G4cout << "Kink is sampled!" << G4endl; 1325 std::vector< G4double > QuarkProbabilitie 1326 theParameters->GetQuarkProbabilitiesA 1327 1328 G4int QuarkInGluon( 1 ); G4double Ksi = G 1329 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) 1330 //G4cout << "Iq " << Iq << G4endl; 1331 if ( Ksi > QuarkProbabilitiesAtGluonSpl 1332 } 1333 //G4cout << "Last Iq " << QuarkInGluon << 1334 G4Parton* Gquark = new G4Parton( QuarkInG 1335 G4Parton* Ganti_quark = new G4Parton( -Qu 1336 //G4cout << "Lorentz " << G4endl; 1337 1338 G4LorentzRotation toCMS( -1 * Phadron.boo 1339 G4LorentzRotation toLab( toCMS.inverse() 1340 //G4cout << "Pstart " << Pstart << G4endl 1341 //G4cout << "Pend " << Pend << G4endl; 1342 //G4cout << "Kink1 " <<PkinkQ1<<G4endl; 1343 //G4cout << "Kink2 " <<PkinkQ2<<G4endl; 1344 //G4cout << "Pstart " << Pstart << G4endl 1345 1346 Pstart.transform( toLab ); start->Set4Mo 1347 PkinkQ1.transform( toLab ); 1348 PkinkQ2.transform( toLab ); 1349 Pend.transform( toLab ); end->Set4Mome 1350 //G4cout << "Pstart " << Pstart << G4endl 1351 //G4cout << "Pend " << Pend << G4endl; 1352 //G4cout << "Defin " << hadron->GetDefini 1353 //G4cout << "Defin " << hadron->GetDefini 1354 1355 //G4int absPDGcode = std::abs( hadron->Ge 1356 G4int absPDGcode = 1500; 1357 if ( start->GetDefinition()->GetParticleS 1358 end->GetDefinition()->GetParticleSub 1359 absPDGcode = 110; 1360 } 1361 //G4cout << "absPDGcode " << absPDGcode < 1362 1363 if ( absPDGcode < 1000 ) { // meson 1364 if ( isProjectile ) { // Projectile 1365 if ( end->GetDefinition()->GetPDGEnco 1366 FirstString = new G4ExcitedString( 1367 SecondString = new G4ExcitedString( 1368 Ganti_quark->Set4Momentum( PkinkQ1 1369 Gquark->Set4Momentum( PkinkQ2 ); 1370 } else { // Anti_Quark on the end 1371 FirstString = new G4ExcitedString( 1372 SecondString = new G4ExcitedString( 1373 Gquark->Set4Momentum( PkinkQ1 ); 1374 Ganti_quark->Set4Momentum( PkinkQ2 1375 } 1376 } else { // Target 1377 if ( end->GetDefinition()->GetPDGEnco 1378 FirstString = new G4ExcitedString( 1379 SecondString = new G4ExcitedString( 1380 Ganti_quark->Set4Momentum( PkinkQ2 1381 Gquark->Set4Momentum( PkinkQ1 ); 1382 } else { // Anti_Quark on the end 1383 FirstString = new G4ExcitedString( 1384 SecondString = new G4ExcitedString( 1385 Gquark->Set4Momentum( PkinkQ2 ); 1386 Ganti_quark->Set4Momentum( PkinkQ1 1387 } 1388 } 1389 } else { // Baryon/AntiBaryon 1390 if ( isProjectile ) { // Projectile 1391 if ( end->GetDefinition()->GetParticl 1392 end->GetDefinition()->GetPDGEnco 1393 FirstString = new G4ExcitedString( 1394 SecondString = new G4ExcitedString( 1395 Gquark->Set4Momentum( PkinkQ1 ); 1396 Ganti_quark->Set4Momentum( PkinkQ2 1397 } else { / 1398 FirstString = new G4ExcitedString( 1399 SecondString = new G4ExcitedString( 1400 Ganti_quark->Set4Momentum( PkinkQ1 1401 Gquark->Set4Momentum( PkinkQ2 ); 1402 } 1403 } else { // Target 1404 if ( end->GetDefinition()->GetParticl 1405 end->GetDefinition()->GetPDGEnco 1406 Gquark->Set4Momentum( PkinkQ1 ); 1407 Ganti_quark->Set4Momentum( PkinkQ2 1408 FirstString = new G4ExcitedString( 1409 SecondString = new G4ExcitedString( 1410 } else { // Anti_DiQuark on the end 1411 FirstString = new G4ExcitedString( 1412 SecondString = new G4ExcitedString( 1413 Gquark->Set4Momentum( PkinkQ2 ); 1414 Ganti_quark->Set4Momentum( PkinkQ1 1415 } 1416 } 1417 } 1418 1419 FirstString->SetTimeOfCreation( hadron->G 1420 FirstString->SetPosition( hadron->GetPosi 1421 SecondString->SetTimeOfCreation( hadron-> 1422 SecondString->SetPosition( hadron->GetPos 1423 1424 } else { // Kink is impossible 1425 1426 if ( isProjectile ) { 1427 FirstString = new G4ExcitedString( end, 1428 } else { 1429 FirstString = new G4ExcitedString( end, 1430 } 1431 FirstString->SetTimeOfCreation( hadron->G 1432 FirstString->SetPosition( hadron->GetPosi 1433 SecondString = 0; 1434 // momenta of string ends 1435 G4LorentzVector HadronMom = hadron->Get4M 1436 G4LorentzVector Pstart1 = G4LorentzVector 1437 G4LorentzVector Pend1 = G4LorentzVector 1438 G4double Pz = HadronMom.pz(); 1439 G4double Eh = HadronMom.e(); 1440 G4double Pt2 = sqr( HadronMom.px() ) + sq 1441 G4double Mt2 = HadronMom.mt2(); 1442 G4double Exp = std::sqrt( sqr(Pz) + ( sqr 1443 G4double Pzq = Pz/2.0 - Exp; 1444 G4double Eq = std::sqrt( sqr(Pzq) + Pt2 1445 G4double Pzqq = Pz/2.0 + Exp; 1446 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt 1447 start->Set4Momentum( Pstart1 ); 1448 end->Set4Momentum( Pend1 ); 1449 Pstart = Pstart1; Pend = Pend1; 1450 1451 } // End of "if (Kink)" 1452 1453 //G4cout << "Quarks in the string at creati 1454 // << " " << FirstString->GetLeftPart 1455 // << FirstString << " " << SecondStr 1456 1457 #ifdef G4_FTFDEBUG 1458 G4cout << " generated string flavors 1459 << end->GetPDGcode() << G4endl << " 1460 << start->Get4Momentum() << "mass : 1461 << " generated string momenta: Diqua 1462 << end->Get4Momentum().mag() << G4en 1463 << Pstart + Pend << G4endl << " Orig 1464 << hadron->Get4Momentum() << " "<<ha 1465 #endif 1466 1467 return; 1468 } 1469 1470 1471 //=========================================== 1472 1473 G4double G4DiffractiveExcitation::ChooseP( G4 1474 // Choose an x between Xmin and Xmax with P 1475 // To be improved... 1476 G4double range = Pmax - Pmin; 1477 if ( Pmin <= 0.0 || range <= 0.0 ) { 1478 G4cout << " Pmin, range : " << Pmin << " 1479 throw G4HadronicException( __FILE__, __LI 1480 "G4Diffractive 1481 } 1482 G4double P = Pmin * G4Pow::GetInstance()->p 1483 //G4double P = (Pmax - Pmin) * G4UniformRan 1484 return P; 1485 } 1486 1487 1488 //=========================================== 1489 1490 G4ThreeVector G4DiffractiveExcitation::Gaussi 1491 // @@ this method is used in FTFModel as w 1492 G4double Pt2( 0.0 ); 1493 if ( AveragePt2 <= 0.0 ) { 1494 Pt2 = 0.0; 1495 } else { 1496 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor 1497 ( G4Ex 1498 } 1499 G4double Pt = std::sqrt( Pt2 ); 1500 G4double phi = G4UniformRand() * twopi; 1501 return G4ThreeVector( Pt * std::cos( phi ), 1502 } 1503 1504 1505 //=========================================== 1506 1507 G4double G4DiffractiveExcitation::GetQuarkFra 1508 G4double z, yf; 1509 const G4int maxNumberOfLoops = 10000; 1510 G4int loopCounter = 0; 1511 do { 1512 z = zmin + G4UniformRand() * (zmax - zmin 1513 yf = z*z + sqr(1.0 - z); 1514 } while ( ( G4UniformRand() > yf ) && 1515 ++loopCounter < maxNumberOfLoops 1516 if ( loopCounter >= maxNumberOfLoops ) { 1517 z = 0.5*(zmin + zmax); // Just something 1518 } 1519 return z; 1520 } 1521 1522 1523 //=========================================== 1524 1525 void G4DiffractiveExcitation::UnpackMeson( co 1526 G4int absIdPDG = std::abs( IdPDG ); 1527 if ( ! ( absIdPDG == 111 || absIdPDG == 221 1528 absIdPDG == 441 || absIdPDG == 443 1529 // All other projectile mesons (including 1530 Q1 = absIdPDG / 100; 1531 Q2 = (absIdPDG % 100) / 10; 1532 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) 1533 if ( IdPDG < 0 ) anti *= -1; 1534 Q1 *= anti; 1535 Q2 *= -1 * anti; 1536 } else { 1537 if ( absIdPDG == 441 || absIdPDG == 443 ) 1538 Q1 = 4; Q2 = -4; 1539 } else if ( absIdPDG == 553 ) { 1540 Q1 = 5; Q2 = -5; 1541 } else { 1542 if ( G4UniformRand() < 0.5 ) { 1543 Q1 = 1; Q2 = -1; 1544 } else { 1545 Q1 = 2; Q2 = -2; 1546 } 1547 } 1548 } 1549 return; 1550 } 1551 1552 1553 //=========================================== 1554 1555 void G4DiffractiveExcitation::UnpackBaryon( G 1556 G 1557 Q1 = IdPDG / 1000; 1558 Q2 = (IdPDG % 1000) / 100; 1559 Q3 = (IdPDG % 100) / 10; 1560 return; 1561 } 1562 1563 1564 //=========================================== 1565 1566 G4int G4DiffractiveExcitation::NewNucleonId( 1567 // Order the three integers in such a way t 1568 G4int TmpQ( 0 ); 1569 if ( Q3 > Q2 ) { 1570 TmpQ = Q2; 1571 Q2 = Q3; 1572 Q3 = TmpQ; 1573 } else if ( Q3 > Q1 ) { 1574 TmpQ = Q1; 1575 Q1 = Q3; 1576 Q3 = TmpQ; 1577 } 1578 if ( Q2 > Q1 ) { 1579 TmpQ = Q1; 1580 Q1 = Q2; 1581 Q2 = TmpQ; 1582 } 1583 // By now Q1 >= Q2 >= Q3 1584 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 1585 return NewCode; 1586 } 1587 1588 1589 //=========================================== 1590 1591 G4DiffractiveExcitation::G4DiffractiveExcitat 1592 throw G4HadronicException( __FILE__, __LINE 1593 "G4DiffractiveEx 1594 } 1595 1596 1597 //=========================================== 1598 1599 const G4DiffractiveExcitation & G4Diffractive 1600 throw G4HadronicException( __FILE__, __LINE 1601 "G4DiffractiveEx 1602 return *this; 1603 } 1604 1605 1606 //=========================================== 1607 1608 G4bool G4DiffractiveExcitation::operator==( c 1609 throw G4HadronicException( __FILE__, __LINE 1610 "G4DiffractiveEx 1611 } 1612 1613 1614 //=========================================== 1615 1616 G4bool G4DiffractiveExcitation::operator!= ( 1617 throw G4HadronicException( __FILE__, __LINE 1618 "G4DiffractiveEx 1619 } 1620 1621