94 #ifdef debugFTFexcitation
98 CommonVariables common;
102 if ( common.Pprojectile.z() < 0.0 )
return false;
104 common.absProjectilePDGcode = std::abs( common.ProjectilePDGcode );
106 G4double ProjectileRapidity = common.Pprojectile.rapidity();
111 common.absTargetPDGcode = std::abs( common.TargetPDGcode );
113 G4double TargetRapidity = common.Ptarget.rapidity();
117 common.S = Psum.
mag2();
118 common.SqrtS = std::sqrt( common.S );
121 G4bool toBePutOnMassShell =
true;
122 common.MminProjectile = common.BrW.GetMinimumMass( projectile->
GetDefinition() );
131 common.M0projectile2 = common.M0projectile * common.M0projectile;
134 if ( common.M0projectile > common.ProjectileDiffStateMinMass ) {
135 common.ProjectileDiffStateMinMass = common.MminProjectile + 220.0*MeV;
136 common.ProjectileNonDiffStateMinMass = common.MminProjectile + 220.0*MeV;
137 if ( common.absProjectilePDGcode > 3000 ) {
138 common.ProjectileDiffStateMinMass += 140.0*MeV;
139 common.ProjectileNonDiffStateMinMass += 140.0*MeV;
142 common.MminTarget = common.BrW.GetMinimumMass( target->
GetDefinition() );
151 common.M0target2 = common.M0target * common.M0target;
154 if ( common.M0target > common.TargetDiffStateMinMass ) {
155 common.TargetDiffStateMinMass = common.MminTarget + 220.0*MeV;
156 common.TargetNonDiffStateMinMass = common.MminTarget + 220.0*MeV;
157 if ( common.absTargetPDGcode > 3000 ) {
158 common.TargetDiffStateMinMass += 140.0*MeV;
159 common.TargetNonDiffStateMinMass += 140.0*MeV;
162 #ifdef debugFTFexcitation
163 G4cout <<
"Proj Targ PDGcodes " << common.ProjectilePDGcode <<
" " << common.TargetPDGcode <<
G4endl
164 <<
"Mprojectile Y " << common.Pprojectile.mag() <<
" " << ProjectileRapidity <<
G4endl
165 <<
"M0projectile Y " << common.M0projectile <<
" " << ProjectileRapidity <<
G4endl;
166 G4cout <<
"Mtarget Y " << common.Ptarget.mag() <<
" " << TargetRapidity <<
G4endl
167 <<
"M0target Y " << common.M0target <<
" " << TargetRapidity <<
G4endl;
168 G4cout <<
"Pproj " << common.Pprojectile <<
G4endl <<
"Ptarget " << common.Ptarget <<
G4endl;
174 if ( Ptmp.
pz() <= 0.0 )
return false;
177 common.toLab = common.toCms.inverse();
178 common.Pprojectile.
transform( common.toCms );
179 common.Ptarget.
transform( common.toCms );
181 G4double SumMasses = common.M0projectile + common.M0target;
182 #ifdef debugFTFexcitation
183 G4cout <<
"SqrtS " << common.SqrtS <<
G4endl <<
"M0pr M0tr SumM " << common.M0projectile
184 <<
" " << common.M0target <<
" " << SumMasses <<
G4endl;
186 if ( common.SqrtS < SumMasses )
return false;
188 common.PZcms2 = (
sqr( common.S ) +
sqr( common.M0projectile2 ) +
sqr( common.M0target2 )
189 - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
190 + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
191 #ifdef debugFTFexcitation
192 G4cout <<
"PZcms2 after toBePutOnMassShell " << common.PZcms2 <<
G4endl;
194 if ( common.PZcms2 < 0.0 )
return false;
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 );
202 common.Pprojectile.setPz( -common.PZcms );
203 common.Ptarget.setPz( common.PZcms );
205 common.Pprojectile.setE( std::sqrt( common.M0projectile2
206 + common.Pprojectile.x() * common.Pprojectile.x()
207 + common.Pprojectile.y() * common.Pprojectile.y()
209 common.Ptarget.setE( std::sqrt( common.M0target2
210 + common.Ptarget.x() * common.Ptarget.x()
211 + common.Ptarget.y() * common.Ptarget.y()
214 #ifdef debugFTFexcitation
215 G4cout <<
"Start --------------------" <<
G4endl <<
"Proj M0 Mdif Mndif " << common.M0projectile
216 <<
" " << common.ProjectileDiffStateMinMass <<
" " << common.ProjectileNonDiffStateMinMass
218 <<
"Targ M0 Mdif Mndif " << common.M0target <<
" " << common.TargetDiffStateMinMass
219 <<
" " << common.TargetNonDiffStateMinMass <<
G4endl <<
"SqrtS " << common.SqrtS <<
G4endl
220 <<
"Proj CMS " << common.Pprojectile <<
G4endl <<
"Targ CMS " << common.Ptarget <<
G4endl;
224 ProjectileRapidity = common.Pprojectile.rapidity();
225 TargetRapidity = common.Ptarget.rapidity();
228 theParameters->
GetProcProb( 4, ProjectileRapidity - TargetRapidity );
229 common.ProbProjectileDiffraction =
230 theParameters->
GetProcProb( 2, ProjectileRapidity - TargetRapidity );
231 common.ProbTargetDiffraction =
232 theParameters->
GetProcProb( 3, ProjectileRapidity - TargetRapidity );
233 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
234 #ifdef debugFTFexcitation
235 G4cout <<
"Proc Probs " << QeNoExc <<
" " << QeExc <<
" "
236 << common.ProbProjectileDiffraction <<
" " << common.ProbTargetDiffraction <<
G4endl
237 <<
"ProjectileRapidity " << ProjectileRapidity <<
G4endl;
240 if ( QeNoExc + QeExc + common.ProbProjectileDiffraction + common.ProbTargetDiffraction > 1.0 ) {
241 QeNoExc = 1.0 - QeExc - common.ProbProjectileDiffraction - common.ProbTargetDiffraction;
243 if ( QeExc + QeNoExc != 0.0 ) {
244 common.ProbExc = QeExc / ( QeExc + QeNoExc );
246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
247 common.ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
248 common.ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
250 #ifdef debugFTFexcitation
251 G4cout <<
"Proc Probs " << QeNoExc <<
" " << QeExc <<
" "
252 << common.ProbProjectileDiffraction <<
" " << common.ProbTargetDiffraction <<
G4endl
253 <<
"ProjectileRapidity " << ProjectileRapidity <<
G4endl;
257 G4int returnCode = 1;
259 returnCode = ExciteParticipants_doChargeExchange( projectile, target, theParameters,
260 theElastic, common );
263 G4bool returnResult =
false;
264 if ( returnCode == 0 ) {
266 }
else if ( returnCode == 1 ) {
268 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
269 #ifdef debugFTFexcitation
271 <<
"Proj M0 MdMin MndMin " << common.M0projectile <<
" "
272 << common.ProjectileDiffStateMinMass <<
" " << common.ProjectileNonDiffStateMinMass
274 <<
"Targ M0 MdMin MndMin " << common.M0target <<
" " << common.TargetDiffStateMinMass
275 <<
" " << common.TargetNonDiffStateMinMass <<
G4endl <<
"SqrtS " << common.SqrtS
277 <<
"Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction <<
" "
278 << common.ProbTargetDiffraction <<
" " << common.ProbOfDiffraction <<
G4endl;
280 if ( common.ProbOfDiffraction != 0.0 ) {
281 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
283 common.ProbProjectileDiffraction = 0.0;
285 #ifdef debugFTFexcitation
286 G4cout <<
"Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction <<
" "
287 << common.ProbTargetDiffraction <<
" " << common.ProbOfDiffraction <<
G4endl;
289 common.ProjectileDiffStateMinMass2 =
sqr( common.ProjectileDiffStateMinMass );
290 common.ProjectileNonDiffStateMinMass2 =
sqr( common.ProjectileNonDiffStateMinMass );
291 common.TargetDiffStateMinMass2 =
sqr( common.TargetDiffStateMinMass );
292 common.TargetNonDiffStateMinMass2 =
sqr( common.TargetNonDiffStateMinMass );
296 returnResult = ExciteParticipants_doDiffraction( projectile, target, theParameters, common );
300 returnResult = ExciteParticipants_doNonDiffraction( projectile, target, theParameters, common );
303 if ( returnResult ) {
304 common.Pprojectile += common.Qmomentum;
305 common.Ptarget -= common.Qmomentum;
307 common.Pprojectile.transform( common.toLab );
308 common.Ptarget.transform( common.toLab );
309 #ifdef debugFTFexcitation
310 G4cout <<
"Mproj " << common.Pprojectile.mag() <<
G4endl <<
"Mtarg " << common.Ptarget.mag()
325G4int G4DiffractiveExcitation::
330 G4DiffractiveExcitation::CommonVariables& common )
const {
337 G4int returnCode = 99;
341 G4double MtestPr = 0.0, MtestTr = 0.0;
343 #ifdef debugFTFexcitation
344 G4cout <<
"Q exchange --------------------------" <<
G4endl;
357 G4int NewProjCode = 0, NewTargCode = 0, ProjQ1 = 0, ProjQ2 = 0, ProjQ3 = 0;
359 if ( common.absProjectilePDGcode < 1000 ) {
360 UnpackMeson( common.ProjectilePDGcode, ProjQ1, ProjQ2 );
362 UnpackBaryon( common.ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;
366 UnpackBaryon( common.TargetPDGcode, TargQ1, TargQ2, TargQ3 );
367 #ifdef debugFTFexcitation
368 G4cout <<
"Proj Quarks " << ProjQ1 <<
" " << ProjQ2 <<
" " << ProjQ3 <<
G4endl
369 <<
"Targ Quarks " << TargQ1 <<
" " << TargQ2 <<
" " << TargQ3 <<
G4endl;
373 G4int ProjExchangeQ = 0, TargExchangeQ = 0;
374 if ( common.absProjectilePDGcode < 1000 ) {
376 G4bool isProjQ1Quark =
false;
377 ProjExchangeQ = ProjQ2;
379 isProjQ1Quark =
true;
380 ProjExchangeQ = ProjQ1;
383 G4int NpossibleStates = 0;
384 if ( ProjExchangeQ != TargQ1 ) NpossibleStates++;
385 if ( ProjExchangeQ != TargQ2 ) NpossibleStates++;
386 if ( ProjExchangeQ != TargQ3 ) NpossibleStates++;
387 G4int Nsampled = (
G4int)G4RandFlat::shootInt(
G4long( NpossibleStates ) )+1;
389 if ( ProjExchangeQ != TargQ1 ) {
390 if ( ++NpossibleStates == Nsampled ) {
391 TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ;
392 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
395 if ( ProjExchangeQ != TargQ2 ) {
396 if ( ++NpossibleStates == Nsampled ) {
397 TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ;
398 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
401 if ( ProjExchangeQ != TargQ3 ) {
402 if ( ++NpossibleStates == Nsampled ) {
403 TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ;
404 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
407 #ifdef debugFTFexcitation
408 G4cout <<
"Exchanged Qs in Pr Tr " << ProjExchangeQ <<
" " << TargExchangeQ <<
G4endl;
411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ2 = std::abs( ProjQ2 );
412 G4bool ProjExcited =
false;
413 const G4int maxNumberOfAttempts = 50;
415 while ( attempts++ < maxNumberOfAttempts ) {
420 if ( aProjQ1 == aProjQ2 ) {
430 }
else if ( aProjQ1 == 3 ) {
435 }
else if ( aProjQ1 == 4 ) {
437 }
else if ( aProjQ1 == 5 ) {
446 }
else if ( aProjQ1 == 3 ) {
448 }
else if ( aProjQ1 == 4 ) {
450 }
else if ( aProjQ1 == 5 ) {
455 if ( aProjQ1 > aProjQ2 ) {
456 NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
458 NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
461 #ifdef debugFTFexcitation
468 if ( aProjQ1 <= 3 && aProjQ2 <= 3 &&
G4UniformRand() < 0.5 ) {
475 G4int value = ProjQ1, absValue = aProjQ1, Qquarks = 0;
476 for (
G4int iQ = 0; iQ < 2; ++iQ ) {
478 value = ProjQ2; absValue = aProjQ2;
480 if ( absValue == 2 || absValue == 4 ) {
481 Qquarks += 2*value/absValue;
483 Qquarks -= value/absValue;
502 if ( Qquarks < 0 || ( Qquarks == 0 && aProjQ1 != aProjQ2 && aProjQ1%2 == 0 ) ) {
505 #ifdef debugFTFexcitation
506 G4cout <<
"NewProjCode +2 or 0 " << NewProjCode <<
G4endl;
507 G4cout<<
"+++++++++++++++++++++++++++++++++++++++"<<
G4endl;
509 G4cout<<
"+++++++++++++++++++++++++++++++++++++++"<<
G4endl;
514 if ( ! TestParticle )
continue;
515 common.MminProjectile = common.BrW.GetMinimumMass( TestParticle );
516 if ( common.SqrtS - common.M0target < common.MminProjectile )
continue;
517 MtestPr = common.BrW.SampleMass( TestParticle, TestParticle->
GetPDGMass()
519 #ifdef debugFTFexcitation
522 <<
"MtestPart MtestPart0 "<<MtestPr<<
" "<<TestParticle->
GetPDGMass()<<
G4endl
523 <<
"M0projectile projectile PDGMass " << common.M0projectile <<
" "
528 NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
529 #ifdef debugFTFexcitation
530 G4cout <<
"New TrQ " << TargQ1 <<
" " << TargQ2 <<
" " << TargQ3 <<
G4endl
531 <<
"NewTargCode " << NewTargCode <<
G4endl;
536 if ( TargQ1 <= 3 && TargQ2 <= 3 && TargQ3 <= 3 ) {
537 if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) {
543 }
else if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
544 NewTargCode += 2; ProjExcited =
true;
547 NewTargCode += 2; ProjExcited =
true;
549 }
else if ( ! ProjExcited &&
551 common.SqrtS > common.M0projectile +
561 if ( NewTargCode == 3124 ||
562 NewTargCode == 3224 ||
563 NewTargCode == 3214 ||
564 NewTargCode == 3114 ||
565 NewTargCode == 3324 ||
566 NewTargCode == 3314 ) {
574 #ifdef debug_heavyHadrons
575 G4int initialNewTargCode = NewTargCode;
577 if ( NewTargCode == 4322 ) NewTargCode = 4232;
578 else if ( NewTargCode == 4312 ) NewTargCode = 4132;
579 else if ( NewTargCode == 5312 ) NewTargCode = 5132;
580 else if ( NewTargCode == 5322 ) NewTargCode = 5232;
581 #ifdef debug_heavyHadrons
582 if ( NewTargCode != initialNewTargCode ) {
583 G4cout <<
"G4DiffractiveExcitation::ExciteParticipants_doChargeExchange : forcing (inexisting in G4)" <<
G4endl
584 <<
"\t target heavy baryon with pdgCode=" << initialNewTargCode
585 <<
" into pdgCode=" << NewTargCode <<
G4endl;
590 if ( ! TestParticle )
continue;
591 #ifdef debugFTFexcitation
594 common.MminTarget = common.BrW.GetMinimumMass( TestParticle );
595 if ( common.SqrtS - MtestPr < common.MminTarget )
continue;
596 MtestTr = common.BrW.SampleMass( TestParticle, TestParticle->
GetPDGMass()
598 if ( common.SqrtS > MtestPr + MtestTr )
break;
601 if ( attempts >= maxNumberOfAttempts )
return returnCode;
603 if ( MtestPr >= common.Pprojectile.mag() || projectile->
GetStatus() != 0 ) {
604 common.M0projectile = MtestPr;
606 #ifdef debugFTFexcitation
607 G4cout <<
"M0projectile After check " << common.M0projectile <<
G4endl;
609 common.M0projectile2 = common.M0projectile * common.M0projectile;
610 common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV;
611 common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV;
612 if ( MtestTr >= common.Ptarget.mag() || target->
GetStatus() != 0 ) {
613 common.M0target = MtestTr;
615 common.M0target2 = common.M0target * common.M0target;
616 #ifdef debugFTFexcitation
617 G4cout <<
"New targ M0 M0^2 " << common.M0target <<
" " << common.M0target2 <<
G4endl;
619 common.TargetDiffStateMinMass = common.M0target + 220.0*MeV;
620 common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV;
627 if ( common.ProjectilePDGcode < 0 )
return 1;
631 G4bool isProjectileExchangedQ =
false;
632 G4int firstQ = TargQ1, secondQ = TargQ2, thirdQ = TargQ3;
633 G4int otherFirstQ = ProjQ1, otherSecondQ = ProjQ2, otherThirdQ = ProjQ3;
635 isProjectileExchangedQ =
true;
636 firstQ = ProjQ1; secondQ = ProjQ2; thirdQ = ProjQ3;
637 otherFirstQ = TargQ1; otherSecondQ = TargQ2; otherThirdQ = TargQ3;
641 G4int exchangedQ = 0;
643 if ( Ksi < 0.333333 ) {
645 }
else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
646 exchangedQ = secondQ;
650 #ifdef debugFTFexcitation
651 G4cout <<
"Exchange Qs isProjectile Q " << isProjectileExchangedQ <<
" " << exchangedQ <<
" ";
657 const G4int MaxCount = 100;
658 G4int count = 0, otherExchangedQ = 0;
660 if ( exchangedQ != otherFirstQ ||
G4UniformRand() < probSame ) {
661 otherExchangedQ = otherFirstQ; otherFirstQ = exchangedQ; exchangedQ = otherExchangedQ;
663 if ( exchangedQ != otherSecondQ ||
G4UniformRand() < probSame ) {
664 otherExchangedQ = otherSecondQ; otherSecondQ = exchangedQ; exchangedQ = otherExchangedQ;
666 if ( exchangedQ != otherThirdQ ||
G4UniformRand() < probSame ) {
667 otherExchangedQ = otherThirdQ; otherThirdQ = exchangedQ; exchangedQ = otherExchangedQ;
671 }
while ( otherExchangedQ == 0 && ++count < MaxCount );
672 if ( count >= MaxCount )
return returnCode;
675 if ( Ksi < 0.333333 ) {
677 }
else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
678 secondQ = exchangedQ;
682 if ( isProjectileExchangedQ ) {
683 ProjQ1 = firstQ; ProjQ2 = secondQ; ProjQ3 = thirdQ;
684 TargQ1 = otherFirstQ; TargQ2 = otherSecondQ; TargQ3 = otherThirdQ;
686 TargQ1 = firstQ; TargQ2 = secondQ; TargQ3 = thirdQ;
687 ProjQ1 = otherFirstQ; ProjQ2 = otherSecondQ; ProjQ3 = otherThirdQ;
689 #ifdef debugFTFexcitation
690 G4cout <<
"Exchange Qs Pr Tr " << ( isProjectileExchangedQ ? exchangedQ : otherExchangedQ )
691 <<
" " << ( isProjectileExchangedQ ? otherExchangedQ : exchangedQ ) <<
G4endl;
694 NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
695 NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
703 for (
G4int iHadron = 0; iHadron < 2; iHadron++ ) {
705 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, codeQ3 = ProjQ3, newHadCode = NewProjCode;
706 G4double massConstraint = common.M0target;
708 if ( iHadron == 1 ) {
709 codeQ1 = TargQ1, codeQ2 = TargQ2, codeQ3 = TargQ3, newHadCode = NewTargCode;
710 massConstraint = common.M0projectile;
713 if ( codeQ1 > 3 || codeQ2 > 3 || codeQ3 > 3 )
continue;
714 if ( codeQ1 == codeQ2 && codeQ1 == codeQ3 ) {
716 }
else if ( isHadronADelta ) {
731 if ( iHadron == 0 ) {
732 NewProjCode = newHadCode;
734 NewTargCode = newHadCode;
737 #ifdef debugFTFexcitation
738 G4cout <<
"NewProjCode NewTargCode " << NewProjCode <<
" " << NewTargCode <<
G4endl;
745 if ( NewProjCode == 3124 ||
746 NewProjCode == 3224 ||
747 NewProjCode == 3214 ||
748 NewProjCode == 3114 ||
749 NewProjCode == 3324 ||
750 NewProjCode == 3314 ) {
755 if ( NewTargCode == 3124 ||
756 NewTargCode == 3224 ||
757 NewTargCode == 3214 ||
758 NewTargCode == 3114 ||
759 NewTargCode == 3324 ||
760 NewTargCode == 3314 ) {
768 #ifdef debug_heavyHadrons
769 G4int initialNewProjCode = NewProjCode, initialNewTargCode = NewTargCode;
771 if ( NewProjCode == 4322 ) NewProjCode = 4232;
772 else if ( NewProjCode == 4312 ) NewProjCode = 4132;
773 else if ( NewProjCode == 5312 ) NewProjCode = 5132;
774 else if ( NewProjCode == 5322 ) NewProjCode = 5232;
775 if ( NewTargCode == 4322 ) NewTargCode = 4232;
776 else if ( NewTargCode == 4312 ) NewTargCode = 4132;
777 else if ( NewTargCode == 5312 ) NewTargCode = 5132;
778 else if ( NewTargCode == 5322 ) NewTargCode = 5232;
779 #ifdef debug_heavyHadrons
780 if ( NewProjCode != initialNewProjCode || NewTargCode != initialNewTargCode ) {
781 G4cout <<
"G4DiffractiveExcitation::ExciteParticipants_doChargeExchange : forcing (inexisting in G4)" <<
G4endl
782 <<
"\t heavy baryon into an existing one:" <<
G4endl;
783 if ( NewProjCode != initialNewProjCode ) {
784 G4cout <<
"\t \t projectile baryon with pdgCode=" << initialNewProjCode
785 <<
" into pdgCode=" << NewProjCode <<
G4endl;
787 if ( NewTargCode != initialNewTargCode ) {
788 G4cout <<
"\t \t target baryon with pdgCode=" << initialNewTargCode
789 <<
" into pdgCode=" << NewTargCode <<
G4endl;
801 G4int firstHadronCode = NewTargCode, secondHadronCode = NewProjCode;
802 G4double massConstraint = common.M0projectile;
803 G4bool isFirstTarget =
true;
806 firstHadron = projectile; secondHadron = target;
807 firstHadronCode = NewProjCode; secondHadronCode = NewTargCode;
808 massConstraint = common.M0target;
809 isFirstTarget =
false;
811 G4double Mtest1st = 0.0, Mtest2nd = 0.0, Mmin1st = 0.0, Mmin2nd = 0.0;
812 for (
int iSamplingCase = 0; iSamplingCase < 2; iSamplingCase++ ) {
814 G4int aHadronCode = firstHadronCode;
815 if ( iSamplingCase == 1 ) {
816 aHadron = secondHadron; aHadronCode = secondHadronCode; massConstraint = Mtest1st;
818 G4double MtestHadron = 0.0, MminHadron = 0.0;
821 if ( ! TestParticle )
return returnCode;
822 MminHadron = common.BrW.GetMinimumMass( TestParticle );
823 if ( common.SqrtS - massConstraint < MminHadron )
return returnCode;
825 MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->
GetPDGMass() );
827 const G4int maxNumberOfAttempts = 50;
829 while ( attempts < maxNumberOfAttempts ) {
831 MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->
GetPDGMass()
833 if ( common.SqrtS < MtestHadron + massConstraint ) {
839 if ( attempts >= maxNumberOfAttempts )
return returnCode;
842 if ( iSamplingCase == 0 ) {
843 Mtest1st = MtestHadron; Mmin1st = MminHadron;
845 Mtest2nd = MtestHadron; Mmin2nd = MminHadron;
848 if ( isFirstTarget ) {
849 MtestTr = Mtest1st; MtestPr = Mtest2nd;
850 common.MminTarget = Mmin1st; common.MminProjectile = Mmin2nd;
852 MtestTr = Mtest2nd; MtestPr = Mtest1st;
853 common.MminTarget = Mmin2nd; common.MminProjectile = Mmin1st;
856 if ( MtestPr != 0.0 ) {
857 common.M0projectile = MtestPr;
858 common.M0projectile2 = common.M0projectile * common.M0projectile;
859 common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV;
860 common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV;
862 if ( MtestTr != 0.0 ) {
863 common.M0target = MtestTr;
864 common.M0target2 = common.M0target * common.M0target;
865 common.TargetDiffStateMinMass = common.M0target + 220.0*MeV;
866 common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV;
873 if ( common.SqrtS < common.M0projectile + common.M0target )
return returnCode;
875 common.PZcms2 = (
sqr( common.S ) +
sqr( common.M0projectile2 ) +
sqr( common.M0target2 )
876 - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
877 + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
878 #ifdef debugFTFexcitation
879 G4cout <<
"At the end// NewProjCode " << NewProjCode <<
G4endl
880 <<
"At the end// NewTargCode " << NewTargCode <<
G4endl
881 <<
"M0pr M0tr SqS " << common.M0projectile <<
" " << common.M0target <<
" "
883 <<
"M0pr2 M0tr2 SqS " << common.M0projectile2 <<
" " << common.M0target2 <<
" "
885 <<
"PZcms2 after the change " << common.PZcms2 <<
G4endl <<
G4endl;
887 if ( common.PZcms2 < 0.0 )
return returnCode;
891 common.PZcms = std::sqrt( common.PZcms2 );
892 common.Pprojectile.setPz( common.PZcms );
893 common.Pprojectile.setE( std::sqrt( common.M0projectile2 + common.PZcms2 ) );
894 common.Ptarget.setPz( -common.PZcms );
895 common.Ptarget.setE( std::sqrt( common.M0target2 + common.PZcms2 ) );
896 #ifdef debugFTFexcitation
898 << common.Ptarget <<
G4endl << common.Pprojectile + common.Ptarget <<
G4endl;
905 if ( common.SqrtS < common.M0projectile + common.TargetDiffStateMinMass ||
906 common.SqrtS < common.ProjectileDiffStateMinMass + common.M0target ||
907 common.ProbOfDiffraction == 0.0 ) common.ProbExc = 0.0;
910 #ifdef debugFTFexcitation
911 G4cout <<
"Make elastic scattering of new hadrons" <<
G4endl;
913 common.Pprojectile.transform( common.toLab );
914 common.Ptarget.transform( common.toLab );
918 #ifdef debugFTFexcitation
919 G4cout <<
"Result of el. scatt " << Result <<
G4endl <<
"Proj Targ and Proj+Targ in Lab"
924 if ( Result ) returnCode = 0;
928 #ifdef debugFTFexcitation
933 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
934 if ( common.ProbOfDiffraction != 0.0 ) {
935 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
936 common.ProbTargetDiffraction /= common.ProbOfDiffraction;
939 return returnCode = 1;
944G4bool G4DiffractiveExcitation::
947 G4DiffractiveExcitation::CommonVariables& common )
const {
951 G4bool isProjectileDiffraction =
false;
953 isProjectileDiffraction =
true;
954 #ifdef debugFTFexcitation
957 common.ProjMassT2 = common.ProjectileDiffStateMinMass2;
958 common.ProjMassT = common.ProjectileDiffStateMinMass;
959 common.TargMassT2 = common.M0target2;
960 common.TargMassT = common.M0target;
962 #ifdef debugFTFexcitation
965 common.ProjMassT2 = common.M0projectile2;
966 common.ProjMassT = common.M0projectile;
967 common.TargMassT2 = common.TargetDiffStateMinMass2;
968 common.TargMassT = common.TargetDiffStateMinMass;
972 if ( common.SqrtS < common.ProjMassT + common.TargMassT )
return false;
974 common.PZcms2 = (
sqr( common.S ) +
sqr( common.ProjMassT2 ) +
sqr( common.TargMassT2 )
975 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
976 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
977 if ( common.PZcms2 < 0.0 )
return false;
978 common.maxPtSquare = common.PZcms2;
981 G4bool loopCondition =
true;
982 G4int whilecount = 0;
986 if ( whilecount > 1000 ) {
991 common.Qmomentum =
G4LorentzVector( GaussianPt( DiffrAveragePt2, common.maxPtSquare ), 0 );
993 if ( isProjectileDiffraction ) {
994 common.ProjMassT2 = common.ProjectileDiffStateMinMass2 + common.Pt2;
995 common.TargMassT2 = common.M0target2 + common.Pt2;
997 common.ProjMassT2 = common.M0projectile2 + common.Pt2;
998 common.TargMassT2 = common.TargetDiffStateMinMass2 + common.Pt2;
1000 common.ProjMassT = std::sqrt( common.ProjMassT2 );
1001 common.TargMassT = std::sqrt( common.TargMassT2 );
1002 if ( common.SqrtS < common.ProjMassT + common.TargMassT )
continue;
1004 common.PZcms2 = (
sqr( common.S ) +
sqr( common.ProjMassT2 ) +
sqr( common.TargMassT2 )
1005 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1006 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1007 if ( common.PZcms2 < 0.0 )
continue;
1009 common.PZcms = std::sqrt( common.PZcms2 );
1010 if ( isProjectileDiffraction ) {
1011 common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
1012 common.PMinusMax = common.SqrtS - common.TargMassT;
1013 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1014 common.TMinusNew = common.SqrtS - common.PMinusNew;
1015 common.Qminus = common.Ptarget.minus() - common.TMinusNew;
1016 common.TPlusNew = common.TargMassT2 / common.TMinusNew;
1017 common.Qplus = common.Ptarget.plus() - common.TPlusNew;
1018 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1019 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1020 loopCondition = ( ( common.Pprojectile + common.Qmomentum ).mag2() <
1021 common.ProjectileDiffStateMinMass2 );
1023 common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
1024 common.TPlusMax = common.SqrtS - common.ProjMassT;
1025 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1026 common.PPlusNew = common.SqrtS - common.TPlusNew;
1027 common.Qplus = common.PPlusNew - common.Pprojectile.plus();
1028 common.PMinusNew = common.ProjMassT2 / common.PPlusNew;
1029 common.Qminus = common.PMinusNew - common.Pprojectile.minus();
1030 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1031 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1032 loopCondition = ( ( common.Ptarget - common.Qmomentum ).mag2() <
1033 common.TargetDiffStateMinMass2 );
1036 }
while ( loopCondition );
1039 if ( isProjectileDiffraction ) {
1052G4bool G4DiffractiveExcitation::
1056 G4DiffractiveExcitation::CommonVariables& common )
const {
1060 #ifdef debugFTFexcitation
1065 common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2;
1066 common.ProjMassT = common.ProjectileNonDiffStateMinMass;
1067 common.TargMassT2 = common.TargetNonDiffStateMinMass2;
1068 common.TargMassT = common.TargetNonDiffStateMinMass;
1069 if ( common.SqrtS < common.ProjMassT + common.TargMassT )
return false;
1071 common.PZcms2 = (
sqr( common.S ) +
sqr( common.ProjMassT2 ) +
sqr( common.TargMassT2 )
1072 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1073 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1074 common.maxPtSquare = common.PZcms2;
1076 G4int whilecount = 0;
1080 if ( whilecount > 1000 ) {
1086 common.maxPtSquare ), 0 );
1088 common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2 + common.Pt2;
1089 common.ProjMassT = std::sqrt( common.ProjMassT2 );
1090 common.TargMassT2 = common.TargetNonDiffStateMinMass2 + common.Pt2;
1091 common.TargMassT = std::sqrt( common.TargMassT2 );
1092 if ( common.SqrtS < common.ProjMassT + common.TargMassT )
continue;
1094 common.PZcms2 =(
sqr( common.S ) +
sqr( common.ProjMassT2 ) +
sqr( common.TargMassT2 )
1095 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1096 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1097 if ( common.PZcms2 < 0.0 )
continue;
1099 common.PZcms = std::sqrt( common.PZcms2 );
1100 common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
1101 common.PMinusMax = common.SqrtS - common.TargMassT;
1102 common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
1103 common.TPlusMax = common.SqrtS - common.ProjMassT;
1107 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1109 common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*
G4UniformRand() + common.PMinusMin;
1112 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1114 common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*
G4UniformRand() + common.TPlusMin;
1118 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1120 common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*
G4UniformRand() + common.TPlusMin;
1123 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1125 common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*
G4UniformRand() + common.PMinusMin;
1129 common.Qminus = common.PMinusNew - common.Pprojectile.minus();
1130 common.Qplus = -( common.TPlusNew - common.Ptarget.plus() );
1131 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1132 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1133 #ifdef debugFTFexcitation
1135 << ( common.Pprojectile + common.Qmomentum ).mag() <<
" "
1136 << common.ProjectileNonDiffStateMinMass <<
G4endl
1137 << ( common.Ptarget - common.Qmomentum ).mag() <<
" "
1138 << common.TargetNonDiffStateMinMass <<
G4endl;
1141 }
while ( ( common.Pprojectile + common.Qmomentum ).mag2() <
1142 common.ProjectileNonDiffStateMinMass2 ||
1143 ( common.Ptarget - common.Qmomentum ).mag2() <
1144 common.TargetNonDiffStateMinMass2 );
1165 if( ! HadronIsString ) hadron->
SplitUp();
1168 if ( start ==
nullptr ) {
1169 G4cout <<
" G4FTFModel::String() Error: No start parton found" <<
G4endl;
1170 FirstString = 0; SecondString = 0;
1175 if ( end ==
nullptr ) {
1176 G4cout <<
" G4FTFModel::String() Error: No end parton found" <<
G4endl;
1177 FirstString = 0; SecondString = 0;
1185 if ( HadronIsString ) {
1186 if ( isProjectile ) {
1210 if ( isProjectile ) {
1218 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 );
1241 if ( Pt > 500.0*MeV ) {
1242 G4double Ymax =
G4Log(
W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1249 if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV;
1250 if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV;
1251 if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV;
1252 if ( PDGcode_startQ == 5 ) Mass_startQ = 4500.0*MeV;
1254 if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV;
1255 if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV;
1256 if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV;
1257 if ( PDGcode_endQ == 5 ) Mass_endQ = 4500.0*MeV;
1259 G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1260 G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1262 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1268 G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1269 G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1271 if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1275 Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 );
1279 Pstart.
setE( x3*
W/2.0 );
1281 Pend.
setE( x1*
W/2.0 );
1283 G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1284 if ( Pkink.
getZ() > 0.0 ) {
1286 PkinkQ1 = XkQ*Pkink;
1288 PkinkQ1 = (1.0 - XkQ)*Pkink;
1292 PkinkQ1 = (1.0 - XkQ)*Pkink;
1294 PkinkQ1 = XkQ*Pkink;
1298 PkinkQ2 = Pkink - PkinkQ1;
1301 std::sqrt(
sqr(
sqr(x1) -
sqr(x3) ) +
sqr( 2.0*x1*x3*CosT13 ) );
1302 G4double Psi = std::acos( Cos2Psi );
1305 if ( isProjectile ) {
1325 std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1329 for (
unsigned int Iq = 0; Iq < 3; Iq++ ) {
1331 if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1356 G4int absPDGcode = 1500;
1363 if ( absPDGcode < 1000 ) {
1364 if ( isProjectile ) {
1390 if ( isProjectile ) {
1426 if ( isProjectile ) {
1442 G4double Exp = std::sqrt(
sqr(Pz) + (
sqr(Mt2) - 4.0*
sqr(Eh)*Pt2/4.0 )/Mt2 )/2.0;
1449 Pstart = Pstart1; Pend = Pend1;
1461 <<
" generated string momenta: Diquark " << end->
Get4Momentum() <<
"mass : "
1463 << Pstart + Pend <<
G4endl <<
" Original "
1477 if ( Pmin <= 0.0 || range <= 0.0 ) {
1478 G4cout <<
" Pmin, range : " << Pmin <<
" , " << range <<
G4endl;
1480 "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1493 if ( AveragePt2 <= 0.0 ) {
1497 (
G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1501 return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1509 const G4int maxNumberOfLoops = 10000;
1510 G4int loopCounter = 0;
1513 yf = z*z +
sqr(1.0 - z);
1515 ++loopCounter < maxNumberOfLoops );
1516 if ( loopCounter >= maxNumberOfLoops ) {
1517 z = 0.5*(zmin + zmax);
1525void G4DiffractiveExcitation::UnpackMeson(
const G4int IdPDG,
G4int& Q1,
G4int& Q2 )
const {
1526 G4int absIdPDG = std::abs( IdPDG );
1527 if ( ! ( absIdPDG == 111 || absIdPDG == 221 || absIdPDG == 331 ||
1528 absIdPDG == 441 || absIdPDG == 443 || absIdPDG == 553 ) ) {
1530 Q1 = absIdPDG / 100;
1531 Q2 = (absIdPDG % 100) / 10;
1532 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1533 if ( IdPDG < 0 ) anti *= -1;
1537 if ( absIdPDG == 441 || absIdPDG == 443 ) {
1539 }
else if ( absIdPDG == 553 ) {
1555void G4DiffractiveExcitation::UnpackBaryon(
G4int IdPDG,
1558 Q2 = (IdPDG % 1000) / 100;
1559 Q3 = (IdPDG % 100) / 10;
1573 }
else if ( Q3 > Q1 ) {
1584 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1593 "G4DiffractiveExcitation copy constructor not meant to be called" );
1601 "G4DiffractiveExcitation = operator not meant to be called" );
1610 "G4DiffractiveExcitation == operator not meant to be called" );
1618 "G4DiffractiveExcitation != operator not meant to be called" );
G4double Y(G4double density)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
CLHEP::HepLorentzRotation G4LorentzRotation
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
Hep3Vector boostVector() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)
G4DiffractiveExcitation()
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
virtual ~G4DiffractiveExcitation()
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
void SetPosition(const G4ThreeVector &aPosition)
void SetTimeOfCreation(G4double aTime)
G4double GetProbLogDistrPrD()
G4double GetProbLogDistr()
G4double GetProjMinNonDiffMass()
G4double GetAvaragePt2ofElasticScattering()
G4double GetTarMinDiffMass()
G4double GetTarMinNonDiffMass()
G4double GetProcProb(const G4int ProcN, const G4double y)
G4double GetDeltaProbAtQuarkExchange()
G4double GetProjMinDiffMass()
G4double GetProbOfSameQuarkExchange()
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
G4int GetPDGiIsospin() const
const G4String & GetParticleType() const
G4double GetPDGMass() const
G4int GetPDGEncoding() const
G4double GetPDGWidth() const
const G4String & GetParticleName() const
const G4String & GetParticleSubType() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * GetDefinition()
void Set4Momentum(const G4LorentzVector &aMomentum)
const G4LorentzVector & Get4Momentum() const
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
void SetStatus(const G4int aStatus)
const G4ParticleDefinition * GetDefinition() const
G4double GetTimeOfCreation()
void Set4Momentum(const G4LorentzVector &a4Momentum)
virtual G4Parton * GetNextParton()=0
void SetDefinition(const G4ParticleDefinition *aDefinition)
const G4LorentzVector & Get4Momentum() const
const G4ThreeVector & GetPosition() const
void IncrementCollisionCount(G4int aCount)
G4int GetSoftCollisionCount()