Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DiffractiveExcitation Class Reference

#include <G4DiffractiveExcitation.hh>

Public Member Functions

 G4DiffractiveExcitation ()
virtual ~G4DiffractiveExcitation ()
virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic, G4bool &isDiffractive) const
virtual void CreateStrings (G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const

Detailed Description

Definition at line 52 of file G4DiffractiveExcitation.hh.

Constructor & Destructor Documentation

◆ G4DiffractiveExcitation()

G4DiffractiveExcitation::G4DiffractiveExcitation ( )

Definition at line 79 of file G4DiffractiveExcitation.cc.

79{}

◆ ~G4DiffractiveExcitation()

G4DiffractiveExcitation::~G4DiffractiveExcitation ( )
virtual

Definition at line 84 of file G4DiffractiveExcitation.cc.

84{}

Member Function Documentation

◆ CreateStrings()

void G4DiffractiveExcitation::CreateStrings ( G4VSplitableHadron * aHadron,
G4bool isProjectile,
G4ExcitedString *& FirstString,
G4ExcitedString *& SecondString,
G4FTFParameters * theParameters ) const
virtual

Definition at line 1157 of file G4DiffractiveExcitation.cc.

1161 {
1162
1163 //G4cout << "Create Strings SplitUp " << hadron << G4endl
1164 // << "Defin " << hadron->GetDefinition() << G4endl
1165 // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
1166
1167 G4bool HadronIsString = hadron->IsSplit();
1168 if( ! HadronIsString ) hadron->SplitUp();
1169
1170 G4Parton* start = hadron->GetNextParton();
1171 if ( start == nullptr ) {
1172 G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1173 FirstString = 0; SecondString = 0;
1174 return;
1175 }
1176
1177 G4Parton* end = hadron->GetNextParton();
1178 if ( end == nullptr ) {
1179 G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1180 FirstString = 0; SecondString = 0;
1181 return;
1182 }
1183
1184 //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1185 // << G4endl
1186 // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1187
1188 if ( HadronIsString ) {
1189 if ( isProjectile ) {
1190 FirstString = new G4ExcitedString( end, start, +1 );
1191 } else {
1192 FirstString = new G4ExcitedString( end, start, -1 );
1193 }
1194 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1195 FirstString->SetPosition( hadron->GetPosition() );
1196 SecondString = 0;
1197 return;
1198 }
1199
1200 G4LorentzVector Phadron = hadron->Get4Momentum();
1201 //G4cout << "String mom " << Phadron << G4endl;
1202 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1203 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1204 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1205 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1206 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1207
1208 G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1209 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1210 //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1211
1212 G4double Wmin( 0.0 );
1213 if ( isProjectile ) {
1214 Wmin = theParameters->GetProjMinDiffMass();
1215 } else {
1216 Wmin = theParameters->GetTarMinDiffMass();
1217 }
1218
1219 G4double W = hadron->Get4Momentum().mag();
1220 G4double W2 = W*W;
1221 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); // x2( 0.0 )
1222 G4bool Kink = false;
1223
1224 if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1225 end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1226 ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1227 end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1228 // Kinky strings are allowed only for qq-q strings;
1229 // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1230 // according to the analysis of Pbar P interactions
1231
1232 if ( W > Wmin ) { // Kink is possible
1233 if ( hadron->GetStatus() == 0 ) {
1234 G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1235 if ( Pt2kink ) {
1236 Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );
1237 } else {
1238 Pt = 0.0;
1239 }
1240 } else {
1241 Pt = 0.0;
1242 }
1243
1244 if ( Pt > 500.0*MeV ) {
1245 G4double Ymax = G4Log( W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1246 G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
1247 x1 = 1.0 - Pt/W * G4Exp( Y );
1248 x3 = 1.0 - Pt/W * G4Exp(-Y );
1249 //x2 = 2.0 - x1 - x3;
1250
1251 G4double Mass_startQ = 650.0*MeV;
1252 if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV; // For constituent up or down quark
1253 if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV; // For constituent strange quark
1254 if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV; // For constituent charm quark
1255 if ( PDGcode_startQ == 5 ) Mass_startQ = 4500.0*MeV; // For constituent bottom quark
1256 G4double Mass_endQ = 650.0*MeV;
1257 if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV; // For constituent up or down quark
1258 if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV; // For constituent strange quark
1259 if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV; // For constituent charm quark
1260 if ( PDGcode_endQ == 5 ) Mass_endQ = 4500.0*MeV; // For constituent bottom quark
1261
1262 G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1263 G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1264 G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
1265 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1266 Kink = false;
1267 } else {
1268 G4double P_1 = std::sqrt( P2_1 );
1269 G4double P_2 = std::sqrt( P2_2 );
1270 G4double P_3 = std::sqrt( P2_3 );
1271 G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1272 G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1273
1274 if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1275 Kink = false;
1276 } else {
1277 Kink = true;
1278 Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated
1279 Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1280 Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1281 Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1282 Pstart.setE( x3*W/2.0 );
1283 Pkink.setE( Pkink.vect().mag() );
1284 Pend.setE( x1*W/2.0 );
1285
1286 G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1287 if ( Pkink.getZ() > 0.0 ) {
1288 if ( XkQ > 0.5 ) {
1289 PkinkQ1 = XkQ*Pkink;
1290 } else {
1291 PkinkQ1 = (1.0 - XkQ)*Pkink;
1292 }
1293 } else {
1294 if ( XkQ > 0.5 ) {
1295 PkinkQ1 = (1.0 - XkQ)*Pkink;
1296 } else {
1297 PkinkQ1 = XkQ*Pkink;
1298 }
1299 }
1300
1301 PkinkQ2 = Pkink - PkinkQ1;
1302 // Minimizing Pt1^2+Pt3^2
1303 G4double Cos2Psi = ( sqr(x1) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1304 std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
1305 G4double Psi = std::acos( Cos2Psi );
1306
1307 G4LorentzRotation Rotate;
1308 if ( isProjectile ) {
1309 Rotate.rotateY( Psi );
1310 } else {
1311 Rotate.rotateY( pi + Psi );
1312 }
1313 Rotate.rotateZ( twopi * G4UniformRand() );
1314 Pstart *= Rotate;
1315 Pkink *= Rotate;
1316 PkinkQ1 *= Rotate;
1317 PkinkQ2 *= Rotate;
1318 Pend *= Rotate;
1319 }
1320 } // End of if ( P2_1 <= 0.0 || P2_3 <= 0.0 )
1321 } // End of if ( Pt > 500.0*MeV )
1322 } // End of if ( W > Wmin ) : check for a kink
1323 } // end of qq-q string selection
1324
1325 if ( Kink ) { // Kink is possible
1326
1327 //G4cout << "Kink is sampled!" << G4endl;
1328 std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1329 theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1330
1331 G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1332 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1333 //G4cout << "Iq " << Iq << G4endl;
1334 if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1335 }
1336 //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1337 G4Parton* Gquark = new G4Parton( QuarkInGluon );
1338 G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1339 //G4cout << "Lorentz " << G4endl;
1340
1341 G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1342 G4LorentzRotation toLab( toCMS.inverse() );
1343 //G4cout << "Pstart " << Pstart << G4endl;
1344 //G4cout << "Pend " << Pend << G4endl;
1345 //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1346 //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1347 //G4cout << "Pstart " << Pstart << G4endl<<G4endl;
1348
1349 Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1350 PkinkQ1.transform( toLab );
1351 PkinkQ2.transform( toLab );
1352 Pend.transform( toLab ); end->Set4Momentum( Pend );
1353 //G4cout << "Pstart " << Pstart << G4endl;
1354 //G4cout << "Pend " << Pend << G4endl;
1355 //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1356 //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1357
1358 //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1359 G4int absPDGcode = 1500;
1360 if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1361 end->GetDefinition()->GetParticleSubType() == "quark" ) {
1362 absPDGcode = 110;
1363 }
1364 //G4cout << "absPDGcode " << absPDGcode << G4endl;
1365
1366 if ( absPDGcode < 1000 ) { // meson
1367 if ( isProjectile ) { // Projectile
1368 if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1369 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1370 SecondString = new G4ExcitedString( Gquark, start , +1 );
1371 Ganti_quark->Set4Momentum( PkinkQ1 );
1372 Gquark->Set4Momentum( PkinkQ2 );
1373 } else { // Anti_Quark on the end
1374 FirstString = new G4ExcitedString( end , Gquark, +1 );
1375 SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1376 Gquark->Set4Momentum( PkinkQ1 );
1377 Ganti_quark->Set4Momentum( PkinkQ2 );
1378 }
1379 } else { // Target
1380 if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1381 FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1382 SecondString = new G4ExcitedString( start , Gquark, -1 );
1383 Ganti_quark->Set4Momentum( PkinkQ2 );
1384 Gquark->Set4Momentum( PkinkQ1 );
1385 } else { // Anti_Quark on the end
1386 FirstString = new G4ExcitedString( Gquark, end , -1 );
1387 SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1388 Gquark->Set4Momentum( PkinkQ2 );
1389 Ganti_quark->Set4Momentum( PkinkQ1 );
1390 }
1391 }
1392 } else { // Baryon/AntiBaryon
1393 if ( isProjectile ) { // Projectile
1394 if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1395 end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1396 FirstString = new G4ExcitedString( end , Gquark, +1 );
1397 SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1398 Gquark->Set4Momentum( PkinkQ1 );
1399 Ganti_quark->Set4Momentum( PkinkQ2 );
1400 } else { // Anti_DiQuark on the end or quark
1401 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1402 SecondString = new G4ExcitedString( Gquark, start , +1 );
1403 Ganti_quark->Set4Momentum( PkinkQ1 );
1404 Gquark->Set4Momentum( PkinkQ2 );
1405 }
1406 } else { // Target
1407 if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1408 end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1409 Gquark->Set4Momentum( PkinkQ1 );
1410 Ganti_quark->Set4Momentum( PkinkQ2 );
1411 FirstString = new G4ExcitedString( end, Gquark, -1 );
1412 SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1413 } else { // Anti_DiQuark on the end or Q
1414 FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1415 SecondString = new G4ExcitedString( start , Gquark, -1 );
1416 Gquark->Set4Momentum( PkinkQ2 );
1417 Ganti_quark->Set4Momentum( PkinkQ1 );
1418 }
1419 }
1420 }
1421
1422 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1423 FirstString->SetPosition( hadron->GetPosition() );
1424 SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1425 SecondString->SetPosition( hadron->GetPosition() );
1426
1427 } else { // Kink is impossible
1428
1429 if ( isProjectile ) {
1430 FirstString = new G4ExcitedString( end, start, +1 );
1431 } else {
1432 FirstString = new G4ExcitedString( end, start, -1 );
1433 }
1434 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1435 FirstString->SetPosition( hadron->GetPosition() );
1436 SecondString = 0;
1437 // momenta of string ends
1438 G4LorentzVector HadronMom = hadron->Get4Momentum();
1439 G4LorentzVector Pstart1 = G4LorentzVector( HadronMom.px()/2.0, HadronMom.py()/2.0, 0.0, 0.0 ); // Quark momentum
1440 G4LorentzVector Pend1 = G4LorentzVector( HadronMom.px()/2.0, HadronMom.py()/2.0, 0.0, 0.0 ); // Di-quark momentum
1441 G4double Pz = HadronMom.pz();
1442 G4double Eh = HadronMom.e();
1443 G4double Pt2 = sqr( HadronMom.px() ) + sqr( HadronMom.py() );
1444 G4double Mt2 = HadronMom.mt2();
1445 G4double Exp = std::sqrt( sqr(Pz) + ( sqr(Mt2) - 4.0*sqr(Eh)*Pt2/4.0 )/Mt2 )/2.0;
1446 G4double Pzq = Pz/2.0 - Exp; Pstart1.setZ( Pzq );
1447 G4double Eq = std::sqrt( sqr(Pzq) + Pt2/4.0 ); Pstart1.setE( Eq );
1448 G4double Pzqq = Pz/2.0 + Exp; Pend1.setZ(Pzqq);
1449 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt2/4.0 ); Pend1.setE(Eqq);
1450 start->Set4Momentum( Pstart1 );
1451 end->Set4Momentum( Pend1 );
1452 Pstart = Pstart1; Pend = Pend1;
1453
1454 } // End of "if (Kink)"
1455
1456 //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1457 // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1458 // << FirstString << " " << SecondString << G4endl;
1459
1460 #ifdef G4_FTFDEBUG
1461 G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1462 << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1463 << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1464 << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1465 << end->Get4Momentum().mag() << G4endl << " sum of ends "
1466 << Pstart + Pend << G4endl << " Original "
1467 << hadron->Get4Momentum() << " "<<hadron->Get4Momentum().mag() << G4endl;
1468 #endif
1469
1470 return;
1471}
Y
G4double Y(G4double density)
Definition G4DNAElectronHoleRecombination.cc:57
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:169
G4LorentzRotation
CLHEP::HepLorentzRotation G4LorentzRotation
Definition G4LorentzRotation.hh:33
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition G4LorentzVector.hh:33
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition LorentzRotation.cc:158
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition LorentzRotation.cc:170
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition LorentzVector.cc:173
G4ExcitedString::SetPosition
void SetPosition(const G4ThreeVector &aPosition)
Definition G4ExcitedString.hh:155
G4ExcitedString::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition G4ExcitedString.hh:143
G4FTFParameters::GetTarMinDiffMass
G4double GetTarMinDiffMass()
Definition G4FTFParameters.hh:456
G4FTFParameters::GetProjMinDiffMass
G4double GetProjMinDiffMass()
Definition G4FTFParameters.hh:448
G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
Definition G4FTFParameters.hh:482
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition G4ParticleDefinition.hh:118
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition G4ParticleDefinition.hh:119
G4Parton::GetDefinition
G4ParticleDefinition * GetDefinition()
Definition G4Parton.hh:161
G4Parton::GetPDGcode
G4int GetPDGcode() const
Definition G4Parton.hh:127
G4Parton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition G4Parton.hh:148
G4Parton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition G4Parton.hh:143
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition G4Pow.cc:41
W
#define W
Definition crc32.c:85
sqr
T sqr(const T &x)
Definition templates.hh:128

◆ ExciteParticipants()

G4bool G4DiffractiveExcitation::ExciteParticipants ( G4VSplitableHadron * aPartner,
G4VSplitableHadron * bPartner,
G4FTFParameters * theParameters,
G4ElasticHNScattering * theElastic,
G4bool & isDiffractive ) const
virtual

Definition at line 89 of file G4DiffractiveExcitation.cc.

93 {
94
95 #ifdef debugFTFexcitation
96 G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
97 #endif
98
99 CommonVariables common;
100
101 // Projectile parameters
102 common.Pprojectile = projectile->Get4Momentum();
103 if ( common.Pprojectile.z() < 0.0 ) return false;
104 common.ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
105 common.absProjectilePDGcode = std::abs( common.ProjectilePDGcode );
106 common.M0projectile = projectile->GetDefinition()->GetPDGMass(); //Uzhi Aug.2019 common.Pprojectile.mag();
107 G4double ProjectileRapidity = common.Pprojectile.rapidity();
108
109 // Target parameters
110 common.Ptarget = target->Get4Momentum();
111 common.TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
112 common.absTargetPDGcode = std::abs( common.TargetPDGcode );
113 common.M0target = target->GetDefinition()->GetPDGMass(); //Uzhi Aug.2019 common.Ptarget.mag();
114 G4double TargetRapidity = common.Ptarget.rapidity();
115
116 // Kinematical properties of the interactions
117 G4LorentzVector Psum = common.Pprojectile + common.Ptarget; // 4-momentum in Lab.
118 common.S = Psum.mag2();
119 common.SqrtS = std::sqrt( common.S );
120
121 // Check off-shellness of the participants
122 G4bool toBePutOnMassShell = true; //Uzhi Aug.2019 false;
123 common.MminProjectile = common.BrW.GetMinimumMass( projectile->GetDefinition() );
124 /* Uzhi Aug.2019
125 if ( common.M0projectile < common.MminProjectile ) {
126 toBePutOnMassShell = true;
127 common.M0projectile = common.BrW.SampleMass( projectile->GetDefinition(),
128 projectile->GetDefinition()->GetPDGMass()
129 + 5.0*projectile->GetDefinition()->GetPDGWidth() );
130 }
131 */
132 common.M0projectile2 = common.M0projectile * common.M0projectile;
133 common.ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass();
134 common.ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass();
135 if ( common.M0projectile > common.ProjectileDiffStateMinMass ) {
136 common.ProjectileDiffStateMinMass = common.MminProjectile + 220.0*MeV; //Uzhi Aug.2019 common.M0projectile + 220.0*MeV;
137 common.ProjectileNonDiffStateMinMass = common.MminProjectile + 220.0*MeV; //Uzhi Aug.2019 common.M0projectile + 220.0*MeV;
138 if ( common.absProjectilePDGcode > 3000 ) { // Strange baryon
139 common.ProjectileDiffStateMinMass += 140.0*MeV;
140 common.ProjectileNonDiffStateMinMass += 140.0*MeV;
141 }
142 }
143 common.MminTarget = common.BrW.GetMinimumMass( target->GetDefinition() );
144 /* Uzhi Aug.2019
145 if ( common.M0target < common.MminTarget ) {
146 toBePutOnMassShell = true;
147 common.M0target = common.BrW.SampleMass( target->GetDefinition(),
148 target->GetDefinition()->GetPDGMass()
149 + 5.0*target->GetDefinition()->GetPDGWidth() );
150 }
151 */
152 common.M0target2 = common.M0target * common.M0target;
153 common.TargetDiffStateMinMass = theParameters->GetTarMinDiffMass();
154 common.TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass();
155 if ( common.M0target > common.TargetDiffStateMinMass ) {
156 common.TargetDiffStateMinMass = common.MminTarget + 220.0*MeV; //Uzhi Aug.2019 common.M0target + 220.0*MeV;
157 common.TargetNonDiffStateMinMass = common.MminTarget + 220.0*MeV; //Uzhi Aug.2019 common.M0target + 220.0*MeV;
158 if ( common.absTargetPDGcode > 3000 ) { // Strange baryon
159 common.TargetDiffStateMinMass += 140.0*MeV;
160 common.TargetNonDiffStateMinMass += 140.0*MeV;
161 }
162 };
163 #ifdef debugFTFexcitation
164 G4cout << "Proj Targ PDGcodes " << common.ProjectilePDGcode << " " << common.TargetPDGcode << G4endl
165 << "Mprojectile Y " << common.Pprojectile.mag() << " " << ProjectileRapidity << G4endl // Uzhi Aug.2019
166 << "M0projectile Y " << common.M0projectile << " " << ProjectileRapidity << G4endl;
167 G4cout << "Mtarget Y " << common.Ptarget.mag() << " " << TargetRapidity << G4endl // Uzhi Aug.2019
168 << "M0target Y " << common.M0target << " " << TargetRapidity << G4endl;
169 G4cout << "Pproj " << common.Pprojectile << G4endl << "Ptarget " << common.Ptarget << G4endl;
170 #endif
171
172 // Transform momenta to cms and then rotate parallel to z axis;
173 common.toCms = G4LorentzRotation( -1 * Psum.boostVector() );
174 G4LorentzVector Ptmp = common.toCms * common.Pprojectile;
175 if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
176 common.toCms.rotateZ( -1*Ptmp.phi() );
177 common.toCms.rotateY( -1*Ptmp.theta() );
178 common.toLab = common.toCms.inverse();
179 common.Pprojectile.transform( common.toCms );
180 common.Ptarget.transform( common.toCms );
181
182 G4double SumMasses = common.M0projectile + common.M0target; // + 220.0*MeV;
183 #ifdef debugFTFexcitation
184 G4cout << "SqrtS " << common.SqrtS << G4endl << "M0pr M0tr SumM " << common.M0projectile
185 << " " << common.M0target << " " << SumMasses << G4endl;
186 #endif
187 if ( common.SqrtS < SumMasses ) return false; // The model cannot work at low energy
188
189 common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
190 - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
191 + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
192 #ifdef debugFTFexcitation
193 G4cout << "PZcms2 after toBePutOnMassShell " << common.PZcms2 << G4endl;
194 #endif
195 if ( common.PZcms2 < 0.0 ) return false; // It can be in an interaction with off-shell nuclear nucleon
196
197 common.PZcms = std::sqrt( common.PZcms2 );
198 if ( toBePutOnMassShell ) {
199 if ( common.Pprojectile.z() > 0.0 ) {
200 common.Pprojectile.setPz( common.PZcms );
201 common.Ptarget.setPz( -common.PZcms );
202 } else {
203 common.Pprojectile.setPz( -common.PZcms );
204 common.Ptarget.setPz( common.PZcms );
205 }
206 common.Pprojectile.setE( std::sqrt( common.M0projectile2
207 + common.Pprojectile.x() * common.Pprojectile.x()
208 + common.Pprojectile.y() * common.Pprojectile.y()
209 + common.PZcms2 ) );
210 common.Ptarget.setE( std::sqrt( common.M0target2
211 + common.Ptarget.x() * common.Ptarget.x()
212 + common.Ptarget.y() * common.Ptarget.y()
213 + common.PZcms2 ) );
214 }
215 #ifdef debugFTFexcitation
216 G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << common.M0projectile
217 << " " << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
218 << G4endl
219 << "Targ M0 Mdif Mndif " << common.M0target << " " << common.TargetDiffStateMinMass
220 << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS << G4endl
221 << "Proj CMS " << common.Pprojectile << G4endl << "Targ CMS " << common.Ptarget << G4endl;
222 #endif
223
224 // Check for possible quark exchange
225 ProjectileRapidity = common.Pprojectile.rapidity();
226 TargetRapidity = common.Ptarget.rapidity();
227 G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity - TargetRapidity );
228 G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity - TargetRapidity ) *
229 theParameters->GetProcProb( 4, ProjectileRapidity - TargetRapidity );
230 common.ProbProjectileDiffraction =
231 theParameters->GetProcProb( 2, ProjectileRapidity - TargetRapidity );
232 common.ProbTargetDiffraction =
233 theParameters->GetProcProb( 3, ProjectileRapidity - TargetRapidity );
234 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
235 #ifdef debugFTFexcitation
236 G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
237 << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
238 << "ProjectileRapidity " << ProjectileRapidity << G4endl;
239 #endif
240
241 if ( QeNoExc + QeExc + common.ProbProjectileDiffraction + common.ProbTargetDiffraction > 1.0 ) {
242 QeNoExc = 1.0 - QeExc - common.ProbProjectileDiffraction - common.ProbTargetDiffraction;
243 }
244 if ( QeExc + QeNoExc != 0.0 ) {
245 common.ProbExc = QeExc / ( QeExc + QeNoExc );
246 }
247 if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
248 common.ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
249 common.ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
250 }
251 #ifdef debugFTFexcitation
252 G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
253 << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
254 << "ProjectileRapidity " << ProjectileRapidity << G4endl;
255 #endif
256
257 // Try out charge exchange
258 G4int returnCode = 1;
259 if ( G4UniformRand() < QeExc + QeNoExc ) {
260 returnCode = ExciteParticipants_doChargeExchange( projectile, target, theParameters,
261 theElastic, common );
262 }
263
264 G4bool returnResult = false;
265 if ( returnCode == 0 ) {
266 returnResult = true; // Successfully ended: no need of extra work
267 } else if ( returnCode == 1 ) {
268
269 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
270 #ifdef debugFTFexcitation
271 G4cout << "Excitation --------------------" << G4endl
272 << "Proj M0 MdMin MndMin " << common.M0projectile << " "
273 << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
274 << G4endl
275 << "Targ M0 MdMin MndMin " << common.M0target << " " << common.TargetDiffStateMinMass
276 << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS
277 << G4endl
278 << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
279 << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
280 #endif
281 if ( common.ProbOfDiffraction != 0.0 ) {
282 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
283 } else {
284 common.ProbProjectileDiffraction = 0.0;
285 }
286 #ifdef debugFTFexcitation
287 G4cout << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
288 << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
289 #endif
290 common.ProjectileDiffStateMinMass2 = sqr( common.ProjectileDiffStateMinMass );
291 common.ProjectileNonDiffStateMinMass2 = sqr( common.ProjectileNonDiffStateMinMass );
292 common.TargetDiffStateMinMass2 = sqr( common.TargetDiffStateMinMass );
293 common.TargetNonDiffStateMinMass2 = sqr( common.TargetNonDiffStateMinMass );
294 // Choose between diffraction and non-diffraction process
295 if ( G4UniformRand() < common.ProbOfDiffraction ) {
296
297 returnResult = ExciteParticipants_doDiffraction( projectile, target, theParameters, common );
298
299 if ( returnResult ) isDiffractive = true;
300
301 } else {
302
303 returnResult = ExciteParticipants_doNonDiffraction( projectile, target, theParameters, common );
304
305 }
306 if ( returnResult ) {
307 common.Pprojectile += common.Qmomentum;
308 common.Ptarget -= common.Qmomentum;
309 // Transform back and update SplitableHadron Participant.
310 common.Pprojectile.transform( common.toLab );
311 common.Ptarget.transform( common.toLab );
312 #ifdef debugFTFexcitation
313 G4cout << "Mproj " << common.Pprojectile.mag() << G4endl << "Mtarg " << common.Ptarget.mag()
314 << G4endl;
315 #endif
316 projectile->Set4Momentum( common.Pprojectile );
317 target->Set4Momentum( common.Ptarget );
318 projectile->IncrementCollisionCount( 1 );
319 target->IncrementCollisionCount( 1 );
320 }
321 }
322
323 return returnResult;
324}
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::rotateZ
HepLorentzVector & rotateZ(double)
Definition LorentzVector.cc:76
CLHEP::HepLorentzVector::rotateY
HepLorentzVector & rotateY(double)
Definition LorentzVector.cc:72
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4FTFParameters::GetProjMinNonDiffMass
G4double GetProjMinNonDiffMass()
Definition G4FTFParameters.hh:452
G4FTFParameters::GetTarMinNonDiffMass
G4double GetTarMinNonDiffMass()
Definition G4FTFParameters.hh:460
G4FTFParameters::GetProcProb
G4double GetProcProb(const G4int ProcN, const G4double y)
Definition G4FTFParameters.cc:904
The documentation for this class was generated from the following files: