Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Classes | 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 1138 of file G4DiffractiveExcitation.cc.

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