Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4DiffractiveExcitation Class Reference

#include <G4DiffractiveExcitation.hh>

Public Member Functions

 G4DiffractiveExcitation ()
 
virtual ~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
 

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 1154 of file G4DiffractiveExcitation.cc.

1158 {
1159
1160 //G4cout << "Create Strings SplitUp " << hadron << G4endl
1161 // << "Defin " << hadron->GetDefinition() << G4endl
1162 // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
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: No start parton found" << G4endl;
1170 FirstString = 0; SecondString = 0;
1171 return;
1172 }
1173
1174 G4Parton* end = hadron->GetNextParton();
1175 if ( end == nullptr ) {
1176 G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1177 FirstString = 0; SecondString = 0;
1178 return;
1179 }
1180
1181 //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1182 // << G4endl
1183 // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1184
1185 if ( HadronIsString ) {
1186 if ( isProjectile ) {
1187 FirstString = new G4ExcitedString( end, start, +1 );
1188 } else {
1189 FirstString = new G4ExcitedString( end, start, -1 );
1190 }
1191 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1192 FirstString->SetPosition( hadron->GetPosition() );
1193 SecondString = 0;
1194 return;
1195 }
1196
1197 G4LorentzVector Phadron = hadron->Get4Momentum();
1198 //G4cout << "String mom " << Phadron << G4endl;
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->GetDefinition()->GetPDGEncoding() );
1206 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1207 //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
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 ); // x2( 0.0 )
1219 G4bool Kink = false;
1220
1221 if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1222 end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1223 ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1224 end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1225 // Kinky strings are allowed only for qq-q strings;
1226 // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1227 // according to the analysis of Pbar P interactions
1228
1229 if ( W > Wmin ) { // Kink is possible
1230 if ( hadron->GetStatus() == 0 ) {
1231 G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1232 if ( Pt2kink ) {
1233 Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );
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::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1243 G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
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_startQ = 325.0*MeV; // For constituent up or down quark
1250 if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV; // For constituent strange quark
1251 if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV; // For constituent charm quark
1252 if ( PDGcode_startQ == 5 ) Mass_startQ = 4500.0*MeV; // For constituent bottom quark
1253 G4double Mass_endQ = 650.0*MeV;
1254 if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV; // For constituent up or down quark
1255 if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV; // For constituent strange quark
1256 if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV; // For constituent charm quark
1257 if ( PDGcode_endQ == 5 ) Mass_endQ = 4500.0*MeV; // For constituent bottom quark
1258
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;
1261 G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
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 - P2_2 ) / (2.0*P_1*P_2);
1269 G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1270
1271 if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1272 Kink = false;
1273 } else {
1274 Kink = true;
1275 Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated
1276 Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1277 Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1278 Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1279 Pstart.setE( x3*W/2.0 );
1280 Pkink.setE( Pkink.vect().mag() );
1281 Pend.setE( x1*W/2.0 );
1282
1283 G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
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) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1301 std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
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 * G4UniformRand() );
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_3 <= 0.0 )
1318 } // End of if ( Pt > 500.0*MeV )
1319 } // End of if ( W > Wmin ) : check for a kink
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 > QuarkProbabilitiesAtGluonSplitUp =
1326 theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1327
1328 G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1329 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1330 //G4cout << "Iq " << Iq << G4endl;
1331 if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1332 }
1333 //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1334 G4Parton* Gquark = new G4Parton( QuarkInGluon );
1335 G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1336 //G4cout << "Lorentz " << G4endl;
1337
1338 G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
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<<G4endl;
1345
1346 Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1347 PkinkQ1.transform( toLab );
1348 PkinkQ2.transform( toLab );
1349 Pend.transform( toLab ); end->Set4Momentum( Pend );
1350 //G4cout << "Pstart " << Pstart << G4endl;
1351 //G4cout << "Pend " << Pend << G4endl;
1352 //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1353 //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1354
1355 //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1356 G4int absPDGcode = 1500;
1357 if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1358 end->GetDefinition()->GetParticleSubType() == "quark" ) {
1359 absPDGcode = 110;
1360 }
1361 //G4cout << "absPDGcode " << absPDGcode << G4endl;
1362
1363 if ( absPDGcode < 1000 ) { // meson
1364 if ( isProjectile ) { // Projectile
1365 if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1366 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1367 SecondString = new G4ExcitedString( Gquark, start , +1 );
1368 Ganti_quark->Set4Momentum( PkinkQ1 );
1369 Gquark->Set4Momentum( PkinkQ2 );
1370 } else { // Anti_Quark on the end
1371 FirstString = new G4ExcitedString( end , Gquark, +1 );
1372 SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1373 Gquark->Set4Momentum( PkinkQ1 );
1374 Ganti_quark->Set4Momentum( PkinkQ2 );
1375 }
1376 } else { // Target
1377 if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1378 FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1379 SecondString = new G4ExcitedString( start , Gquark, -1 );
1380 Ganti_quark->Set4Momentum( PkinkQ2 );
1381 Gquark->Set4Momentum( PkinkQ1 );
1382 } else { // Anti_Quark on the end
1383 FirstString = new G4ExcitedString( Gquark, end , -1 );
1384 SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1385 Gquark->Set4Momentum( PkinkQ2 );
1386 Ganti_quark->Set4Momentum( PkinkQ1 );
1387 }
1388 }
1389 } else { // Baryon/AntiBaryon
1390 if ( isProjectile ) { // Projectile
1391 if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1392 end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1393 FirstString = new G4ExcitedString( end , Gquark, +1 );
1394 SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1395 Gquark->Set4Momentum( PkinkQ1 );
1396 Ganti_quark->Set4Momentum( PkinkQ2 );
1397 } else { // Anti_DiQuark on the end or quark
1398 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1399 SecondString = new G4ExcitedString( Gquark, start , +1 );
1400 Ganti_quark->Set4Momentum( PkinkQ1 );
1401 Gquark->Set4Momentum( PkinkQ2 );
1402 }
1403 } else { // Target
1404 if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1405 end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1406 Gquark->Set4Momentum( PkinkQ1 );
1407 Ganti_quark->Set4Momentum( PkinkQ2 );
1408 FirstString = new G4ExcitedString( end, Gquark, -1 );
1409 SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1410 } else { // Anti_DiQuark on the end or Q
1411 FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1412 SecondString = new G4ExcitedString( start , Gquark, -1 );
1413 Gquark->Set4Momentum( PkinkQ2 );
1414 Ganti_quark->Set4Momentum( PkinkQ1 );
1415 }
1416 }
1417 }
1418
1419 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1420 FirstString->SetPosition( hadron->GetPosition() );
1421 SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1422 SecondString->SetPosition( hadron->GetPosition() );
1423
1424 } else { // Kink is impossible
1425
1426 if ( isProjectile ) {
1427 FirstString = new G4ExcitedString( end, start, +1 );
1428 } else {
1429 FirstString = new G4ExcitedString( end, start, -1 );
1430 }
1431 FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1432 FirstString->SetPosition( hadron->GetPosition() );
1433 SecondString = 0;
1434 // momenta of string ends
1435 G4LorentzVector HadronMom = hadron->Get4Momentum();
1436 G4LorentzVector Pstart1 = G4LorentzVector( HadronMom.px()/2.0, HadronMom.py()/2.0, 0.0, 0.0 ); // Quark momentum
1437 G4LorentzVector Pend1 = G4LorentzVector( HadronMom.px()/2.0, HadronMom.py()/2.0, 0.0, 0.0 ); // Di-quark momentum
1438 G4double Pz = HadronMom.pz();
1439 G4double Eh = HadronMom.e();
1440 G4double Pt2 = sqr( HadronMom.px() ) + sqr( HadronMom.py() );
1441 G4double Mt2 = HadronMom.mt2();
1442 G4double Exp = std::sqrt( sqr(Pz) + ( sqr(Mt2) - 4.0*sqr(Eh)*Pt2/4.0 )/Mt2 )/2.0;
1443 G4double Pzq = Pz/2.0 - Exp; Pstart1.setZ( Pzq );
1444 G4double Eq = std::sqrt( sqr(Pzq) + Pt2/4.0 ); Pstart1.setE( Eq );
1445 G4double Pzqq = Pz/2.0 + Exp; Pend1.setZ(Pzqq);
1446 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt2/4.0 ); Pend1.setE(Eqq);
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 creation" << FirstString->GetRightParton()->GetPDGcode()
1454 // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1455 // << FirstString << " " << SecondString << G4endl;
1456
1457 #ifdef G4_FTFDEBUG
1458 G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1459 << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1460 << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1461 << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1462 << end->Get4Momentum().mag() << G4endl << " sum of ends "
1463 << Pstart + Pend << G4endl << " Original "
1464 << hadron->Get4Momentum() << " "<<hadron->Get4Momentum().mag() << G4endl;
1465 #endif
1466
1467 return;
1468}
Y
G4double Y(G4double density)
Definition G4DNAElectronHoleRecombination.cc:57
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:227
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 ) const
virtual

Definition at line 89 of file G4DiffractiveExcitation.cc.

92 {
93
94 #ifdef debugFTFexcitation
95 G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
96 #endif
97
98 CommonVariables common;
99
100 // Projectile parameters
101 common.Pprojectile = projectile->Get4Momentum();
102 if ( common.Pprojectile.z() < 0.0 ) return false;
103 common.ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
104 common.absProjectilePDGcode = std::abs( common.ProjectilePDGcode );
105 common.M0projectile = projectile->GetDefinition()->GetPDGMass(); //Uzhi Aug.2019 common.Pprojectile.mag();
106 G4double ProjectileRapidity = common.Pprojectile.rapidity();
107
108 // Target parameters
109 common.Ptarget = target->Get4Momentum();
110 common.TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
111 common.absTargetPDGcode = std::abs( common.TargetPDGcode );
112 common.M0target = target->GetDefinition()->GetPDGMass(); //Uzhi Aug.2019 common.Ptarget.mag();
113 G4double TargetRapidity = common.Ptarget.rapidity();
114
115 // Kinematical properties of the interactions
116 G4LorentzVector Psum = common.Pprojectile + common.Ptarget; // 4-momentum in Lab.
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 Aug.2019 false;
122 common.MminProjectile = common.BrW.GetMinimumMass( projectile->GetDefinition() );
123 /* Uzhi Aug.2019
124 if ( common.M0projectile < common.MminProjectile ) {
125 toBePutOnMassShell = true;
126 common.M0projectile = common.BrW.SampleMass( projectile->GetDefinition(),
127 projectile->GetDefinition()->GetPDGMass()
128 + 5.0*projectile->GetDefinition()->GetPDGWidth() );
129 }
130 */
131 common.M0projectile2 = common.M0projectile * common.M0projectile;
132 common.ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass();
133 common.ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass();
134 if ( common.M0projectile > common.ProjectileDiffStateMinMass ) {
135 common.ProjectileDiffStateMinMass = common.MminProjectile + 220.0*MeV; //Uzhi Aug.2019 common.M0projectile + 220.0*MeV;
136 common.ProjectileNonDiffStateMinMass = common.MminProjectile + 220.0*MeV; //Uzhi Aug.2019 common.M0projectile + 220.0*MeV;
137 if ( common.absProjectilePDGcode > 3000 ) { // Strange baryon
138 common.ProjectileDiffStateMinMass += 140.0*MeV;
139 common.ProjectileNonDiffStateMinMass += 140.0*MeV;
140 }
141 }
142 common.MminTarget = common.BrW.GetMinimumMass( target->GetDefinition() );
143 /* Uzhi Aug.2019
144 if ( common.M0target < common.MminTarget ) {
145 toBePutOnMassShell = true;
146 common.M0target = common.BrW.SampleMass( target->GetDefinition(),
147 target->GetDefinition()->GetPDGMass()
148 + 5.0*target->GetDefinition()->GetPDGWidth() );
149 }
150 */
151 common.M0target2 = common.M0target * common.M0target;
152 common.TargetDiffStateMinMass = theParameters->GetTarMinDiffMass();
153 common.TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass();
154 if ( common.M0target > common.TargetDiffStateMinMass ) {
155 common.TargetDiffStateMinMass = common.MminTarget + 220.0*MeV; //Uzhi Aug.2019 common.M0target + 220.0*MeV;
156 common.TargetNonDiffStateMinMass = common.MminTarget + 220.0*MeV; //Uzhi Aug.2019 common.M0target + 220.0*MeV;
157 if ( common.absTargetPDGcode > 3000 ) { // Strange baryon
158 common.TargetDiffStateMinMass += 140.0*MeV;
159 common.TargetNonDiffStateMinMass += 140.0*MeV;
160 }
161 };
162 #ifdef debugFTFexcitation
163 G4cout << "Proj Targ PDGcodes " << common.ProjectilePDGcode << " " << common.TargetPDGcode << G4endl
164 << "Mprojectile Y " << common.Pprojectile.mag() << " " << ProjectileRapidity << G4endl // Uzhi Aug.2019
165 << "M0projectile Y " << common.M0projectile << " " << ProjectileRapidity << G4endl;
166 G4cout << "Mtarget Y " << common.Ptarget.mag() << " " << TargetRapidity << G4endl // Uzhi Aug.2019
167 << "M0target Y " << common.M0target << " " << TargetRapidity << G4endl;
168 G4cout << "Pproj " << common.Pprojectile << G4endl << "Ptarget " << common.Ptarget << G4endl;
169 #endif
170
171 // Transform momenta to cms and then rotate parallel to z axis;
172 common.toCms = G4LorentzRotation( -1 * Psum.boostVector() );
173 G4LorentzVector Ptmp = common.toCms * common.Pprojectile;
174 if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
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 + common.M0target; // + 220.0*MeV;
182 #ifdef debugFTFexcitation
183 G4cout << "SqrtS " << common.SqrtS << G4endl << "M0pr M0tr SumM " << common.M0projectile
184 << " " << common.M0target << " " << SumMasses << G4endl;
185 #endif
186 if ( common.SqrtS < SumMasses ) return false; // The model cannot work at low energy
187
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;
193 #endif
194 if ( common.PZcms2 < 0.0 ) return false; // It can be in an interaction with off-shell nuclear nucleon
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.M0projectile2
206 + common.Pprojectile.x() * common.Pprojectile.x()
207 + common.Pprojectile.y() * common.Pprojectile.y()
208 + common.PZcms2 ) );
209 common.Ptarget.setE( std::sqrt( common.M0target2
210 + common.Ptarget.x() * common.Ptarget.x()
211 + common.Ptarget.y() * common.Ptarget.y()
212 + common.PZcms2 ) );
213 }
214 #ifdef debugFTFexcitation
215 G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << common.M0projectile
216 << " " << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
217 << G4endl
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;
221 #endif
222
223 // Check for possible quark exchange
224 ProjectileRapidity = common.Pprojectile.rapidity();
225 TargetRapidity = common.Ptarget.rapidity();
226 G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity - TargetRapidity );
227 G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity - TargetRapidity ) *
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;
238 #endif
239
240 if ( QeNoExc + QeExc + common.ProbProjectileDiffraction + common.ProbTargetDiffraction > 1.0 ) {
241 QeNoExc = 1.0 - QeExc - common.ProbProjectileDiffraction - common.ProbTargetDiffraction;
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 - QeExc - QeNoExc );
248 common.ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
249 }
250 #ifdef debugFTFexcitation
251 G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " "
252 << common.ProbProjectileDiffraction << " " << common.ProbTargetDiffraction << G4endl
253 << "ProjectileRapidity " << ProjectileRapidity << G4endl;
254 #endif
255
256 // Try out charge exchange
257 G4int returnCode = 1;
258 if ( G4UniformRand() < QeExc + QeNoExc ) {
259 returnCode = ExciteParticipants_doChargeExchange( projectile, target, theParameters,
260 theElastic, common );
261 }
262
263 G4bool returnResult = false;
264 if ( returnCode == 0 ) {
265 returnResult = true; // Successfully ended: no need of extra work
266 } else if ( returnCode == 1 ) {
267
268 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
269 #ifdef debugFTFexcitation
270 G4cout << "Excitation --------------------" << G4endl
271 << "Proj M0 MdMin MndMin " << common.M0projectile << " "
272 << common.ProjectileDiffStateMinMass << " " << common.ProjectileNonDiffStateMinMass
273 << G4endl
274 << "Targ M0 MdMin MndMin " << common.M0target << " " << common.TargetDiffStateMinMass
275 << " " << common.TargetNonDiffStateMinMass << G4endl << "SqrtS " << common.SqrtS
276 << G4endl
277 << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
278 << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
279 #endif
280 if ( common.ProbOfDiffraction != 0.0 ) {
281 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
282 } else {
283 common.ProbProjectileDiffraction = 0.0;
284 }
285 #ifdef debugFTFexcitation
286 G4cout << "Prob: ProjDiff TargDiff + Sum " << common.ProbProjectileDiffraction << " "
287 << common.ProbTargetDiffraction << " " << common.ProbOfDiffraction << G4endl;
288 #endif
289 common.ProjectileDiffStateMinMass2 = sqr( common.ProjectileDiffStateMinMass );
290 common.ProjectileNonDiffStateMinMass2 = sqr( common.ProjectileNonDiffStateMinMass );
291 common.TargetDiffStateMinMass2 = sqr( common.TargetDiffStateMinMass );
292 common.TargetNonDiffStateMinMass2 = sqr( common.TargetNonDiffStateMinMass );
293 // Choose between diffraction and non-diffraction process
294 if ( G4UniformRand() < common.ProbOfDiffraction ) {
295
296 returnResult = ExciteParticipants_doDiffraction( projectile, target, theParameters, common );
297
298 } else {
299
300 returnResult = ExciteParticipants_doNonDiffraction( projectile, target, theParameters, common );
301
302 }
303 if ( returnResult ) {
304 common.Pprojectile += common.Qmomentum;
305 common.Ptarget -= common.Qmomentum;
306 // Transform back and update SplitableHadron Participant.
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()
311 << G4endl;
312 #endif
313 projectile->Set4Momentum( common.Pprojectile );
314 target->Set4Momentum( common.Ptarget );
315 projectile->IncrementCollisionCount( 1 );
316 target->IncrementCollisionCount( 1 );
317 }
318 }
319
320 return returnResult;
321}
G4LorentzRotation
CLHEP::HepLorentzRotation G4LorentzRotation
Definition G4LorentzRotation.hh:33
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: