Geant4-11
Data Structures | Public Member Functions | Private Member Functions
G4DiffractiveExcitation Class Reference

#include <G4DiffractiveExcitation.hh>

Data Structures

struct  CommonVariables
 

Public Member Functions

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

Private Member Functions

G4double ChooseP (G4double Pmin, G4double Pmax) const
 
G4int ExciteParticipants_doChargeExchange (G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic, CommonVariables &common) const
 
G4bool ExciteParticipants_doDiffraction (G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, CommonVariables &common) const
 
G4bool ExciteParticipants_doNonDiffraction (G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, CommonVariables &common) const
 
 G4DiffractiveExcitation (const G4DiffractiveExcitation &right)
 
G4ThreeVector GaussianPt (G4double AveragePt2, G4double maxPtSquare) const
 
G4double GetQuarkFractionOfKink (G4double zmin, G4double zmax) const
 
G4double LambdaF (G4double sqrM, G4double sqrM1, G4double sqrM2) const
 
G4int NewNucleonId (G4int Q1, G4int Q2, G4int Q3) const
 
G4bool operator!= (const G4DiffractiveExcitation &right) const
 
const G4DiffractiveExcitation & operator= (const G4DiffractiveExcitation &right)
 
G4bool operator== (const G4DiffractiveExcitation &right) const
 
void UnpackBaryon (G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const
 
void UnpackMeson (G4int IdPDG, G4int &Q1, G4int &Q2) const
 

Detailed Description

Definition at line 52 of file G4DiffractiveExcitation.hh.

Constructor & Destructor Documentation

◆ G4DiffractiveExcitation() [1/2]

G4DiffractiveExcitation::G4DiffractiveExcitation ( )

Definition at line 79 of file G4DiffractiveExcitation.cc.

79{}

◆ ~G4DiffractiveExcitation()

G4DiffractiveExcitation::~G4DiffractiveExcitation ( )
virtual

Definition at line 84 of file G4DiffractiveExcitation.cc.

84{}

◆ G4DiffractiveExcitation() [2/2]

G4DiffractiveExcitation::G4DiffractiveExcitation ( const G4DiffractiveExcitation &  right)
private

Definition at line 1562 of file G4DiffractiveExcitation.cc.

1562 {
1563 throw G4HadronicException( __FILE__, __LINE__,
1564 "G4DiffractiveExcitation copy constructor not meant to be called" );
1565}

Member Function Documentation

◆ ChooseP()

G4double G4DiffractiveExcitation::ChooseP ( G4double  Pmin,
G4double  Pmax 
) const
private

Definition at line 1444 of file G4DiffractiveExcitation.cc.

1444 {
1445 // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1446 // To be improved...
1447 G4double range = Pmax - Pmin;
1448 if ( Pmin <= 0.0 || range <= 0.0 ) {
1449 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1450 throw G4HadronicException( __FILE__, __LINE__,
1451 "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1452 }
1453 G4double P = Pmin * G4Pow::GetInstance()->powA( Pmax/Pmin, G4UniformRand() );
1454 //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1455 return P;
1456}
G4double
double G4double
Definition: G4Types.hh:83
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
anonymous_namespace{G4ChipsKaonMinusInelasticXS.cc}::Pmin
const G4double Pmin
Definition: G4ChipsKaonMinusInelasticXS.cc:56
anonymous_namespace{G4ChipsKaonMinusInelasticXS.cc}::Pmax
const G4double Pmax
Definition: G4ChipsKaonMinusInelasticXS.cc:57
P
static double P[]
Definition: nf_gammaFunctions.cc:89

References G4cout, G4endl, G4UniformRand, G4Pow::GetInstance(), P, anonymous_namespace{G4ChipsKaonMinusInelasticXS.cc}::Pmax, anonymous_namespace{G4ChipsKaonMinusInelasticXS.cc}::Pmin, and G4Pow::powA().

Referenced by ExciteParticipants_doDiffraction(), and ExciteParticipants_doNonDiffraction().

◆ 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
G4double Y(G4double density)
Definition: G4DNAElectronHoleRecombination.cc:54
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
twopi
static constexpr double twopi
Definition: G4SIunits.hh:56
MeV
static constexpr double MeV
Definition: G4SIunits.hh:200
pi
static constexpr double pi
Definition: G4SIunits.hh:55
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
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
G4DiffractiveExcitation::GetQuarkFractionOfKink
G4double GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
Definition: G4DiffractiveExcitation.cc:1478
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:710
G4FTFParameters::GetPt2Kink
G4double GetPt2Kink()
Definition: G4FTFParameters.hh:732
G4FTFParameters::GetProjMinDiffMass
G4double GetProjMinDiffMass()
Definition: G4FTFParameters.hh:702
G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
Definition: G4FTFParameters.hh:736
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:130
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:131
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
sqr
T sqr(const T &x)
Definition: templates.hh:128

References CLHEP::HepLorentzVector::boostVector(), CLHEP::HepLorentzVector::e(), G4cout, G4endl, G4Exp(), G4Log(), G4UniformRand, G4VSplitableHadron::Get4Momentum(), G4Parton::Get4Momentum(), G4Parton::GetDefinition(), G4Pow::GetInstance(), G4VSplitableHadron::GetNextParton(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetParticleType(), G4Parton::GetPDGcode(), G4ParticleDefinition::GetPDGEncoding(), G4VSplitableHadron::GetPosition(), G4FTFParameters::GetProjMinDiffMass(), G4FTFParameters::GetPt2Kink(), GetQuarkFractionOfKink(), G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp(), G4VSplitableHadron::GetStatus(), G4FTFParameters::GetTarMinDiffMass(), G4VSplitableHadron::GetTimeOfCreation(), CLHEP::HepLorentzVector::getZ(), CLHEP::HepLorentzRotation::inverse(), CLHEP::HepLorentzVector::mag(), CLHEP::Hep3Vector::mag(), MeV, CLHEP::HepLorentzVector::mt2(), pi, CLHEP::HepLorentzVector::px(), CLHEP::HepLorentzVector::py(), CLHEP::HepLorentzVector::pz(), CLHEP::HepLorentzRotation::rotateY(), CLHEP::HepLorentzRotation::rotateZ(), G4Parton::Set4Momentum(), CLHEP::HepLorentzVector::setE(), G4ExcitedString::SetPosition(), CLHEP::HepLorentzVector::setPx(), CLHEP::HepLorentzVector::setPy(), CLHEP::HepLorentzVector::setPz(), G4ExcitedString::SetTimeOfCreation(), CLHEP::HepLorentzVector::setZ(), G4VSplitableHadron::SplitUp(), sqr(), CLHEP::HepLorentzVector::transform(), twopi, CLHEP::HepLorentzVector::vect(), and Y().

Referenced by G4FTFModel::BuildStrings().

◆ 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
G4DiffractiveExcitation::ExciteParticipants_doDiffraction
G4bool ExciteParticipants_doDiffraction(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, CommonVariables &common) const
Definition: G4DiffractiveExcitation.cc:929
G4DiffractiveExcitation::ExciteParticipants_doChargeExchange
G4int ExciteParticipants_doChargeExchange(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic, CommonVariables &common) const
Definition: G4DiffractiveExcitation.cc:326
G4DiffractiveExcitation::ExciteParticipants_doNonDiffraction
G4bool ExciteParticipants_doNonDiffraction(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4FTFParameters *theParameters, CommonVariables &common) const
Definition: G4DiffractiveExcitation.cc:1037
G4FTFParameters::GetProjMinNonDiffMass
G4double GetProjMinNonDiffMass()
Definition: G4FTFParameters.hh:706
G4FTFParameters::GetTarMinNonDiffMass
G4double GetTarMinNonDiffMass()
Definition: G4FTFParameters.hh:714
G4FTFParameters::GetProcProb
G4double GetProcProb(const G4int ProcN, const G4double y)
Definition: G4FTFParameters.cc:848
common
static int FASTCALL common(PROLOG_STATE *state, int tok)
Definition: xmlrole.cc:1309

References CLHEP::HepLorentzVector::boostVector(), common(), ExciteParticipants_doChargeExchange(), ExciteParticipants_doDiffraction(), ExciteParticipants_doNonDiffraction(), G4cout, G4endl, G4UniformRand, G4VSplitableHadron::Get4Momentum(), G4VSplitableHadron::GetDefinition(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4FTFParameters::GetProcProb(), G4FTFParameters::GetProjMinDiffMass(), G4FTFParameters::GetProjMinNonDiffMass(), G4FTFParameters::GetTarMinDiffMass(), G4FTFParameters::GetTarMinNonDiffMass(), G4VSplitableHadron::IncrementCollisionCount(), CLHEP::HepLorentzVector::mag2(), MeV, CLHEP::HepLorentzVector::phi(), CLHEP::HepLorentzVector::pz(), G4VSplitableHadron::Set4Momentum(), sqr(), and CLHEP::HepLorentzVector::theta().

Referenced by G4FTFModel::ExciteParticipants().

◆ ExciteParticipants_doChargeExchange()

G4int G4DiffractiveExcitation::ExciteParticipants_doChargeExchange ( G4VSplitableHadron *  projectile,
G4VSplitableHadron *  target,
G4FTFParameters *  theParameters,
G4ElasticHNScattering *  theElastic,
G4DiffractiveExcitation::CommonVariables &  common 
) const
private

Definition at line 325 of file G4DiffractiveExcitation.cc.

330 {
331 // First of the three utility methods used only by ExciteParticipants:
332 // it does the sampling for the charge-exchange case.
333 // This method returns an integer code - instead of a boolean, with the following meaning:
334 // "0" : successfully ended and nothing else needs to be done;
335 // "1" : successfully completed, but the work needs to be continued;
336 // "99" : unsuccessfully ended, nothing else can be done.
337 G4int returnCode = 99;
338
339 G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
340 G4ParticleDefinition* TestParticle = 0;
341 G4double MtestPr = 0.0, MtestTr = 0.0;
342
343 #ifdef debugFTFexictation
344 G4cout << "Q exchange --------------------------" << G4endl;
345 #endif
346
347 // The target hadron is always a nucleon (i.e. either a proton or a neutron,
348 // never an antinucleon), therefore only a quark (not an antiquark) can be
349 // exchanged between the projectile hadron and the target hadron (otherwise
350 // we could get a quark-quark-antiquark system which cannot be a bound state).
351 // This implies that any projectile meson or anti-meson - given that it has
352 // a constituent quark in all cases - can have charge exchange with a target
353 // hadron. Instead, any projectile anti-baryon can never have charge exchange
354 // with a target hadron (because it has only constituent anti-quarks);
355 // projectile baryons, instead can have charge exchange with a target hadron.
356
357 G4int NewProjCode = 0, NewTargCode = 0, ProjQ1 = 0, ProjQ2 = 0, ProjQ3 = 0;
358 // Projectile unpacking
359 if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
360 UnpackMeson( common.ProjectilePDGcode, ProjQ1, ProjQ2 );
361 } else { // projectile is baryon
362 UnpackBaryon( common.ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
363 }
364 // Target unpacking
365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;
366 UnpackBaryon( common.TargetPDGcode, TargQ1, TargQ2, TargQ3 );
367 #ifdef debugFTFexictation
368 G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
369 << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
370 #endif
371
372 // Sampling of exchanged quarks
373 G4int ProjExchangeQ = 0, TargExchangeQ = 0;
374 if ( common.absProjectilePDGcode < 1000 ) { // projectile is meson
375
376 G4bool isProjQ1Quark = false;
377 ProjExchangeQ = ProjQ2;
378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark
379 isProjQ1Quark = true;
380 ProjExchangeQ = ProjQ1;
381 }
382 // Exchange of non-identical quarks is allowed
383 G4int NpossibleStates = 0;
384 if ( ProjExchangeQ != TargQ1 ) NpossibleStates++;
385 if ( ProjExchangeQ != TargQ2 ) NpossibleStates++;
386 if ( ProjExchangeQ != TargQ3 ) NpossibleStates++;
387 G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
388 NpossibleStates = 0;
389 if ( ProjExchangeQ != TargQ1 ) {
390 if ( ++NpossibleStates == Nsampled ) {
391 TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ;
392 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
393 }
394 }
395 if ( ProjExchangeQ != TargQ2 ) {
396 if ( ++NpossibleStates == Nsampled ) {
397 TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ;
398 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
399 }
400 }
401 if ( ProjExchangeQ != TargQ3 ) {
402 if ( ++NpossibleStates == Nsampled ) {
403 TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ;
404 isProjQ1Quark ? ProjQ1 = TargExchangeQ : ProjQ2 = TargExchangeQ;
405 }
406 }
407 #ifdef debugFTFexictation
408 G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
409 #endif
410
411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ2 = std::abs( ProjQ2 );
412 G4bool ProjExcited = false;
413 const G4int maxNumberOfAttempts = 50;
414 G4int attempts = 0;
415 while ( attempts++ < maxNumberOfAttempts ) { /* Loop checking, 10.08.2015, A.Ribon */
416
417 // Determination of a new projectile ID which satisfies energy-momentum conservation
418 G4double Ksi = G4UniformRand();
419 if ( aProjQ1 == aProjQ2 ) {
420 if ( aProjQ1 < 3 ) {
421 NewProjCode = 111; // Pi0-meson
422 if ( Ksi < 0.5 ) {
423 NewProjCode = 221; // Eta-meson
424 if ( Ksi < 0.25 ) {
425 NewProjCode = 331; // Eta'-meson
426 }
427 }
428 } else if ( aProjQ1 == 3 ) {
429 NewProjCode = 221; // Eta-meson
430 if ( Ksi < 0.5 ) {
431 NewProjCode = 331; // Eta'-meson
432 }
433 } else if ( aProjQ1 == 4 ) {
434 NewProjCode = 441; // Eta_c
435 } else if ( aProjQ1 == 5 ) {
436 NewProjCode = 553; // Upsilon
437 }
438 } else {
439 if ( aProjQ1 > aProjQ2 ) {
440 NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
441 } else {
442 NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
443 }
444 }
445 #ifdef debugFTFexictation
446 G4cout << "NewProjCode " << NewProjCode << G4endl;
447 #endif
448 // Decide (with 50% probability) whether the projectile hadrons is excited,
449 // but not in the case of charmed and bottom hadrons (because in Geant4
450 // there are no excited charmed and bottom states).
451 ProjExcited = false;
452 if ( aProjQ1 <= 3 && aProjQ2 <= 3 && G4UniformRand() < 0.5 ) {
453 NewProjCode += 2; // Excited meson (J=1 -> 2*J+1=2+1=3, last digit of the PDG code)
454 ProjExcited = true;
455 }
456
457 // So far we have used the absolute values of the PDG codes of the two constituent quarks:
458 // now look at their signed values to set properly the signed of the meson's PDG code.
459 G4int value = ProjQ1, absValue = aProjQ1, Qquarks = 0;
460 for ( G4int iQ = 0; iQ < 2; ++iQ ) { // 0 : first quark; 1 : second quark
461 if ( iQ == 1 ) {
462 value = ProjQ2; absValue = aProjQ2;
463 }
464 if ( absValue == 2 || absValue == 4 ) { // u or c
465 Qquarks += 2*value/absValue; // u, c : positively charged +2 (e/3 unit)
466 } else {
467 Qquarks -= value/absValue; // d, s, b : negatively charged -1 (e/3 unit)
468 }
469 }
470 // If Qquarks is positive, the sign of NewProjCode is fine.
471 // If Qquarks is negative, then the sign of NewProjCode needs to be reversed.
472 // If Qquarks is zero, then we need to distinguish between 3 cases:
473 // 1. If aProjQ1 and aProjQ2 are the same, then the sign of NewProjCode is fine
474 // (because the antiparticle is the same as the particle, e.g. eta, eta').
475 // 2. If aProjQ1 and aProjQ2 are not the same, given that Qquarks is zero,
476 // we have only two possibilities:
477 // a. aProjQ1 and aProjQ2 are two different down-type quarks, i.e.
478 // (s,d) or (b,d), or (b,s). In this case, the sign of NewProjCode
479 // is fine (because the heaviest of the two down-type quarks has
480 // to be anti-quark belonging to the initial projectile, which
481 // implies a meson with positive PDG code, e.g. B0 (bbar,d), Bs (bbar,s).
482 // b. aProjQ1 and aProjQ2 are two different up-type quarks, i.e. (u,c).
483 // The heaviest of the two (c) has to be an anti-quark (cbar) left
484 // in the projectile, therefore the sign of NewProjCode needs to be
485 // reverse: 421 -> -421 anti_D0 (cbar,u)
486 if ( Qquarks < 0 || ( Qquarks == 0 && aProjQ1 != aProjQ2 && aProjQ1%2 == 0 ) ) {
487 NewProjCode *= -1;
488 }
489 #ifdef debugFTFexictation
490 G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
491 G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
492 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquarks<<G4endl;
493 G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
494 #endif
495
496 // Proj
497 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
498 if ( ! TestParticle ) continue;
499 common.MminProjectile = common.BrW.GetMinimumMass( TestParticle );
500 if ( common.SqrtS - common.M0target < common.MminProjectile ) continue;
501 MtestPr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
502 + 5.0*TestParticle->GetPDGWidth() );
503 #ifdef debugFTFexictation
504 G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()
505 << G4endl
506 << "MtestPart MtestPart0 "<<MtestPr<<" "<<TestParticle->GetPDGMass()<<G4endl
507 << "M0projectile projectile PDGMass " << common.M0projectile << " "
508 << projectile->GetDefinition()->GetPDGMass() << G4endl;
509 #endif
510
511 // Targ
512 NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
513 #ifdef debugFTFexictation
514 G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
515 << "NewTargCode " << NewTargCode << G4endl;
516 #endif
517
518 // If the target has not heavy (charm or bottom) constituent quark,
519 // see whether a Delta isobar can be created.
520 if ( TargQ1 <= 3 && TargQ2 <= 3 && TargQ3 <= 3 ) {
521 if ( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 ) { // Lambda or Sigma0 ?
522 if ( G4UniformRand() < 0.5 ) {
523 NewTargCode += 2;
524 } else if ( G4UniformRand() < 0.75 ) {
525 NewTargCode = 3122; // Lambda
526 }
527 } else if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
528 NewTargCode += 2; ProjExcited = true; // Create Delta isobar
529 } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
530 if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
531 NewTargCode += 2; ProjExcited = true;
532 }
533 } else if ( ! ProjExcited &&
534 G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
535 common.SqrtS > common.M0projectile + // Delta mass
536 G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass() ) {
537 NewTargCode += 2; // Create Delta isobar
538 }
539 }
540
541 // Protection against:
542 // - Lambda* (i.e. excited Lambda state) NOT existing in PDG , -> Lambda
543 // - Sigma* (i.e. excited Sigma hyperon states) NOT existing in Geant4 -> Sigma
544 // - Xi* (i.e. excited Xi hyperon states) NOT existing in Geant4 -> Xi
545 if ( NewTargCode == 3124 || // Lambda* NOT existing in PDG !
546 NewTargCode == 3224 || // Sigma*+ NOT existing in Geant4
547 NewTargCode == 3214 || // Sigma*0 NOT existing in Geant4
548 NewTargCode == 3114 || // Sigma*- NOT existing in Geant4
549 NewTargCode == 3324 || // Xi*0 NOT existing in Geant4
550 NewTargCode == 3314 ) { // Xi*- NOT existing in Geant4
551 //G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange INEXISTING TARGET PDGcode="
552 // << NewTargCode << G4endl;
553 NewTargCode -= 2; // Corresponding ground-state hyperon
554 }
555
556 // Special treatment for charmed and bottom baryons : in Geant4 there are no Xi_c' and Xi_b'
557 // so we need to transform them by hand to Xi_c and Xi_b, respectively.
558 #ifdef debug_heavyHadrons
559 G4int initialNewTargCode = NewTargCode;
560 #endif
561 if ( NewTargCode == 4322 ) NewTargCode = 4232; // Xi_c'+ -> Xi_c+
562 else if ( NewTargCode == 4312 ) NewTargCode = 4132; // Xi_c'0 -> Xi_c0
563 else if ( NewTargCode == 5312 ) NewTargCode = 5132; // Xi_b'- -> Xi_b-
564 else if ( NewTargCode == 5322 ) NewTargCode = 5232; // Xi_b'0 -> Xi_b0
565 #ifdef debug_heavyHadrons
566 if ( NewTargCode != initialNewTargCode ) {
567 G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange : forcing (inexisting in G4)" << G4endl
568 << "\t target heavy baryon with pdgCode=" << initialNewTargCode
569 << " into pdgCode=" << NewTargCode << G4endl;
570 }
571 #endif
572
573 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
574 if ( ! TestParticle ) continue;
575 #ifdef debugFTFexictation
576 G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
577 #endif
578 common.MminTarget = common.BrW.GetMinimumMass( TestParticle );
579 if ( common.SqrtS - MtestPr < common.MminTarget ) continue;
580 MtestTr = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
581 + 5.0*TestParticle->GetPDGWidth() );
582 if ( common.SqrtS > MtestPr + MtestTr ) break;
583
584 } // End of while loop
585 if ( attempts >= maxNumberOfAttempts ) return returnCode; // unsuccessfully ended, nothing else can be done
586
587 if ( MtestPr >= common.Pprojectile.mag() || projectile->GetStatus() != 0 ) {
588 common.M0projectile = MtestPr;
589 }
590 #ifdef debugFTFexictation
591 G4cout << "M0projectile After check " << common.M0projectile << G4endl;
592 #endif
593 common.M0projectile2 = common.M0projectile * common.M0projectile;
594 common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
595 common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
596 if ( MtestTr >= common.Ptarget.mag() || target->GetStatus() != 0 ) {
597 common.M0target = MtestTr;
598 }
599 common.M0target2 = common.M0target * common.M0target;
600 #ifdef debugFTFexictation
601 G4cout << "New targ M0 M0^2 " << common.M0target << " " << common.M0target2 << G4endl;
602 #endif
603 common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
604 common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
605
606 } else { // of the if ( common.absProjectilePDGcode < 1000 ) ; the projectile is baryon now
607
608 // If it is a projectile anti-baryon, no quark exchange is possible with a target hadron,
609 // therefore returns immediately 1 (which means "successfully completed, but the work
610 // needs to be continued").
611 if ( common.ProjectilePDGcode < 0 ) return 1;
612
613 // Choose randomly, with equal probability, whether to consider the quarks of the
614 // projectile or target hadron for selecting the flavour of the exchanged quark.
615 G4bool isProjectileExchangedQ = false;
616 G4int firstQ = TargQ1, secondQ = TargQ2, thirdQ = TargQ3;
617 G4int otherFirstQ = ProjQ1, otherSecondQ = ProjQ2, otherThirdQ = ProjQ3;
618 if ( G4UniformRand() < 0.5 ) {
619 isProjectileExchangedQ = true;
620 firstQ = ProjQ1; secondQ = ProjQ2; thirdQ = ProjQ3;
621 otherFirstQ = TargQ1; otherSecondQ = TargQ2; otherThirdQ = TargQ3;
622 }
623 // Choose randomly, with equal probability, which of the three quarks of the
624 // selected (projectile or target) hadron is the exhanged quark.
625 G4int exchangedQ = 0;
626 G4double Ksi = G4UniformRand();
627 if ( Ksi < 0.333333 ) {
628 exchangedQ = firstQ;
629 } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
630 exchangedQ = secondQ;
631 } else {
632 exchangedQ = thirdQ;
633 }
634 #ifdef debugFTFexictation
635 G4cout << "Exchange Qs isProjectile Q " << isProjectileExchangedQ << " " << exchangedQ << " ";
636 #endif
637 // The exchanged quarks (one of the projectile hadron and one of the target hadron)
638 // are always accepted if they have different flavour, else (i.e. when they have the
639 // same flavour) they are accepted only with a specified probability.
640 G4double probSame = theParameters->GetProbOfSameQuarkExchange();
641 const G4int MaxCount = 100;
642 G4int count = 0, otherExchangedQ = 0;
643 do {
644 if ( exchangedQ != otherFirstQ || G4UniformRand() < probSame ) {
645 otherExchangedQ = otherFirstQ; otherFirstQ = exchangedQ; exchangedQ = otherExchangedQ;
646 } else {
647 if ( exchangedQ != otherSecondQ || G4UniformRand() < probSame ) {
648 otherExchangedQ = otherSecondQ; otherSecondQ = exchangedQ; exchangedQ = otherExchangedQ;
649 } else {
650 if ( exchangedQ != otherThirdQ || G4UniformRand() < probSame ) {
651 otherExchangedQ = otherThirdQ; otherThirdQ = exchangedQ; exchangedQ = otherExchangedQ;
652 }
653 }
654 }
655 } while ( otherExchangedQ == 0 && ++count < MaxCount );
656 if ( count >= MaxCount ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
657 // Swap (between projectile and target hadron) the two quarks that have been sampled
658 // as "exchanged" quarks.
659 if ( Ksi < 0.333333 ) {
660 firstQ = exchangedQ;
661 } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
662 secondQ = exchangedQ;
663 } else {
664 thirdQ = exchangedQ;
665 }
666 if ( isProjectileExchangedQ ) {
667 ProjQ1 = firstQ; ProjQ2 = secondQ; ProjQ3 = thirdQ;
668 TargQ1 = otherFirstQ; TargQ2 = otherSecondQ; TargQ3 = otherThirdQ;
669 } else {
670 TargQ1 = firstQ; TargQ2 = secondQ; TargQ3 = thirdQ;
671 ProjQ1 = otherFirstQ; ProjQ2 = otherSecondQ; ProjQ3 = otherThirdQ;
672 }
673 #ifdef debugFTFexictation
674 G4cout << "Exchange Qs Pr Tr " << ( isProjectileExchangedQ ? exchangedQ : otherExchangedQ )
675 << " " << ( isProjectileExchangedQ ? otherExchangedQ : exchangedQ ) << G4endl;
676 #endif
677
678 NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
679 NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
680 // Decide whether the new projectile hadron is a Delta particle;
681 // then decide whether the new target hadron is a Delta particle.
682 // Notice that a Delta particle has the last PDG digit "4" (because its spin is 3/2),
683 // whereas a nucleon has "2" (because its spin is 1/2).
684 // If the new projectile hadron or the new target hadron has a heavy (c or b)
685 // constituent quark, then skip this part (because Geant4 does not have
686 // excited charm and bottom hadrons).
687 for ( G4int iHadron = 0; iHadron < 2; iHadron++ ) {
688 // First projectile hadron, then target hadron
689 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, codeQ3 = ProjQ3, newHadCode = NewProjCode;
690 G4double massConstraint = common.M0target;
691 G4bool isHadronADelta = ( projectile->GetDefinition()->GetPDGiIsospin() == 3 );
692 if ( iHadron == 1 ) { // Target hadron
693 codeQ1 = TargQ1, codeQ2 = TargQ2, codeQ3 = TargQ3, newHadCode = NewTargCode;
694 massConstraint = common.M0projectile;
695 isHadronADelta = ( target->GetDefinition()->GetPDGiIsospin() == 3 );
696 }
697 if ( codeQ1 > 3 || codeQ2 > 3 || codeQ3 > 3 ) continue; // No excited charm or bottom states in Geant4
698 if ( codeQ1 == codeQ2 && codeQ1 == codeQ3 ) { // The three quarks are the same
699 newHadCode += 2; // Delta++ (uuu) or Delta- (ddd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
700 } else if ( isHadronADelta ) { // Hadron (projectile or target) was Delta
701 if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
702 newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
703 } else {
704 newHadCode += 0; // No delta (so the last PDG digit remains 2)
705 }
706 } else { // Hadron (projectile or target) was Nucleon
707 if ( G4UniformRand() < DeltaProbAtQuarkExchange &&
708 common.SqrtS > G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass()
709 + massConstraint ) {
710 newHadCode += 2; // Delta+ (uud) or Delta0 (udd) : spin 3/2, last PDG digit = 4 (so +2 wrt p or n)
711 } else {
712 newHadCode += 0; // No delta (so the last PDG digit remains 2)
713 }
714 }
715 if ( iHadron == 0 ) { // Projectile hadron
716 NewProjCode = newHadCode;
717 } else { // Target hadron
718 NewTargCode = newHadCode;
719 }
720 }
721 #ifdef debugFTFexictation
722 G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
723 #endif
724
725 // Protection against:
726 // - Lambda* (i.e. excited Lambda state) NOT existing in PDG , -> Lambda
727 // - Sigma* (i.e. excited Sigma hyperon states) NOT existing in Geant4 -> Sigma
728 // - Xi* (i.e. excited Xi hyperon states) NOT existing in Geant4 -> Xi
729 if ( NewProjCode == 3124 || // Lambda* NOT existing in PDG !
730 NewProjCode == 3224 || // Sigma*+ NOT existing in Geant4
731 NewProjCode == 3214 || // Sigma*0 NOT existing in Geant4
732 NewProjCode == 3114 || // Sigma*- NOT existing in Geant4
733 NewProjCode == 3324 || // Xi*0 NOT existing in Geant4
734 NewProjCode == 3314 ) { // Xi*- NOT existing in Geant4
735 //G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange INEXISTING PROJECTILE PDGcode="
736 // << NewProjCode << G4endl;
737 NewProjCode -= 2; // Corresponding ground-state hyperon
738 }
739 if ( NewTargCode == 3124 || // Lambda* NOT existing in PDG !
740 NewTargCode == 3224 || // Sigma*+ NOT existing in Geant4
741 NewTargCode == 3214 || // Sigma*0 NOT existing in Geant4
742 NewTargCode == 3114 || // Sigma*- NOT existing in Geant4
743 NewTargCode == 3324 || // Xi*0 NOT existing in Geant4
744 NewTargCode == 3314 ) { // Xi*- NOT existing in Geant4
745 //G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange INEXISTING TARGET PDGcode="
746 // << NewTargCode << G4endl;
747 NewTargCode -= 2; // Corresponding ground-state hyperon
748 }
749
750 // Special treatment for charmed and bottom baryons : in Geant4 there are no Xi_c' and Xi_b'
751 // so we need to transform them by hand to the, respectively, Xi_c and Xi_b.
752 #ifdef debug_heavyHadrons
753 G4int initialNewProjCode = NewProjCode, initialNewTargCode = NewTargCode;
754 #endif
755 if ( NewProjCode == 4322 ) NewProjCode = 4232; // Xi_c'+ -> Xi_c+
756 else if ( NewProjCode == 4312 ) NewProjCode = 4132; // Xi_c'0 -> Xi_c0
757 else if ( NewProjCode == 5312 ) NewProjCode = 5132; // Xi_b'- -> Xi_b-
758 else if ( NewProjCode == 5322 ) NewProjCode = 5232; // Xi_b'0 -> Xi_b0
759 if ( NewTargCode == 4322 ) NewTargCode = 4232; // Xi_c'+ -> Xi_c+
760 else if ( NewTargCode == 4312 ) NewTargCode = 4132; // Xi_c'0 -> Xi_c0
761 else if ( NewTargCode == 5312 ) NewTargCode = 5132; // Xi_b'- -> Xi_b-
762 else if ( NewTargCode == 5322 ) NewTargCode = 5232; // Xi_b'0 -> Xi_b0
763 #ifdef debug_heavyHadrons
764 if ( NewProjCode != initialNewProjCode || NewTargCode != initialNewTargCode ) {
765 G4cout << "G4DiffractiveExcitation::ExciteParticipants_doChargeExchange : forcing (inexisting in G4)" << G4endl
766 << "\t heavy baryon into an existing one:" << G4endl;
767 if ( NewProjCode != initialNewProjCode ) {
768 G4cout << "\t \t projectile baryon with pdgCode=" << initialNewProjCode
769 << " into pdgCode=" << NewProjCode << G4endl;
770 }
771 if ( NewTargCode != initialNewTargCode ) {
772 G4cout << "\t \t target baryon with pdgCode=" << initialNewTargCode
773 << " into pdgCode=" << NewTargCode << G4endl;
774 }
775 }
776 #endif
777
778 // Sampling of the masses of the projectile and target nucleons.
779 // Because of energy conservation, the ordering of the sampling matters:
780 // randomly, half of the time we sample first the target nucleon mass and
781 // then the projectile nucleon mass, and the other half of the time we
782 // sample first the projectile nucleon mass and then the target nucleon mass.
783 G4VSplitableHadron* firstHadron = target;
784 G4VSplitableHadron* secondHadron = projectile;
785 G4int firstHadronCode = NewTargCode, secondHadronCode = NewProjCode;
786 G4double massConstraint = common.M0projectile;
787 G4bool isFirstTarget = true;
788 if ( G4UniformRand() < 0.5 ) {
789 // Sample first the projectile nucleon mass, then the target nucleon mass.
790 firstHadron = projectile; secondHadron = target;
791 firstHadronCode = NewProjCode; secondHadronCode = NewTargCode;
792 massConstraint = common.M0target;
793 isFirstTarget = false;
794 }
795 G4double Mtest1st = 0.0, Mtest2nd = 0.0, Mmin1st = 0.0, Mmin2nd = 0.0;
796 for ( int iSamplingCase = 0; iSamplingCase < 2; iSamplingCase++ ) {
797 G4VSplitableHadron* aHadron = firstHadron;
798 G4int aHadronCode = firstHadronCode;
799 if ( iSamplingCase == 1 ) { // Second nucleon mass sampling
800 aHadron = secondHadron; aHadronCode = secondHadronCode; massConstraint = Mtest1st;
801 }
802 G4double MtestHadron = 0.0, MminHadron = 0.0;
803 if ( aHadron->GetStatus() == 1 || aHadron->GetStatus() == 2 ) {
804 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( aHadronCode );
805 if ( ! TestParticle ) return returnCode; // Not possible to find such a hadron: unsuccessfully ended, nothing else can be done
806 MminHadron = common.BrW.GetMinimumMass( TestParticle );
807 if ( common.SqrtS - massConstraint < MminHadron ) return returnCode; // Kinematically impossible: unsuccessfully ended, nothing else can be done
808 if ( TestParticle->GetPDGWidth() == 0.0 ) {
809 MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass() );
810 } else {
811 const G4int maxNumberOfAttempts = 50;
812 G4int attempts = 0;
813 while ( attempts < maxNumberOfAttempts ) {
814 attempts++;
815 MtestHadron = common.BrW.SampleMass( TestParticle, TestParticle->GetPDGMass()
816 + 5.0*TestParticle->GetPDGWidth() );
817 if ( common.SqrtS < MtestHadron + massConstraint ) {
818 continue; // Kinematically unacceptable: try again
819 } else {
820 break; // Kinematically acceptable: the mass sampling is successful
821 }
822 }
823 if ( attempts >= maxNumberOfAttempts ) return returnCode; // All attempts failed: unsuccessfully ended, nothing else can be done
824 }
825 }
826 if ( iSamplingCase == 0 ) {
827 Mtest1st = MtestHadron; Mmin1st = MminHadron;
828 } else {
829 Mtest2nd = MtestHadron; Mmin2nd = MminHadron;
830 }
831 } // End for loop on the two sampling cases (1st and 2nd)
832 if ( isFirstTarget ) {
833 MtestTr = Mtest1st; MtestPr = Mtest2nd;
834 common.MminTarget = Mmin1st; common.MminProjectile = Mmin2nd;
835 } else {
836 MtestTr = Mtest2nd; MtestPr = Mtest1st;
837 common.MminTarget = Mmin2nd; common.MminProjectile = Mmin1st;
838 }
839
840 if ( MtestPr != 0.0 ) {
841 common.M0projectile = MtestPr;
842 common.M0projectile2 = common.M0projectile * common.M0projectile;
843 common.ProjectileDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
844 common.ProjectileNonDiffStateMinMass = common.M0projectile + 220.0*MeV; // 220 MeV=m_pi+80 MeV
845 }
846 if ( MtestTr != 0.0 ) {
847 common.M0target = MtestTr;
848 common.M0target2 = common.M0target * common.M0target;
849 common.TargetDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
850 common.TargetNonDiffStateMinMass = common.M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
851 }
852
853 } // End of if ( common.absProjectilePDGcode < 1000 )
854
855 // If we assume that final state hadrons after the charge exchange will be
856 // in the ground states, we have to put
857 if ( common.SqrtS < common.M0projectile + common.M0target ) return returnCode; // unsuccessfully ended, nothing else can be done
858
859 common.PZcms2 = ( sqr( common.S ) + sqr( common.M0projectile2 ) + sqr( common.M0target2 )
860 - 2.0 * ( common.S * ( common.M0projectile2 + common.M0target2 )
861 + common.M0projectile2 * common.M0target2 ) ) / 4.0 / common.S;
862 #ifdef debugFTFexictation
863 G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
864 << "At the end// NewTargCode " << NewTargCode << G4endl
865 << "M0pr M0tr SqS " << common.M0projectile << " " << common.M0target << " "
866 << common.SqrtS << G4endl
867 << "M0pr2 M0tr2 SqS " << common.M0projectile2 << " " << common.M0target2 << " "
868 << common.SqrtS << G4endl
869 << "PZcms2 after the change " << common.PZcms2 << G4endl << G4endl;
870 #endif
871 if ( common.PZcms2 < 0.0 ) return returnCode; // It can be if energy is not sufficient for Delta;
872 // unsuccessfully ended, nothing else can be done
873 projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
874 target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
875 common.PZcms = std::sqrt( common.PZcms2 );
876 common.Pprojectile.setPz( common.PZcms );
877 common.Pprojectile.setE( std::sqrt( common.M0projectile2 + common.PZcms2 ) );
878 common.Ptarget.setPz( -common.PZcms );
879 common.Ptarget.setE( std::sqrt( common.M0target2 + common.PZcms2 ) );
880 #ifdef debugFTFexictation
881 G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << common.Pprojectile << G4endl
882 << common.Ptarget << G4endl << common.Pprojectile + common.Ptarget << G4endl;
883 #endif
884
885 if ( projectile->GetStatus() != 0 ) projectile->SetStatus( 2 );
886 if ( target->GetStatus() != 0 ) target->SetStatus( 2 );
887
888 // Check for possible excitation of the participants
889 if ( common.SqrtS < common.M0projectile + common.TargetDiffStateMinMass ||
890 common.SqrtS < common.ProjectileDiffStateMinMass + common.M0target ||
891 common.ProbOfDiffraction == 0.0 ) common.ProbExc = 0.0;
892
893 if ( G4UniformRand() > common.ProbExc ) { // Make elastic scattering
894 #ifdef debugFTFexictation
895 G4cout << "Make elastic scattering of new hadrons" << G4endl;
896 #endif
897 common.Pprojectile.transform( common.toLab );
898 common.Ptarget.transform( common.toLab );
899 projectile->Set4Momentum( common.Pprojectile );
900 target->Set4Momentum( common.Ptarget );
901 G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
902 #ifdef debugFTFexictation
903 G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
904 << G4endl << projectile->Get4Momentum() << G4endl << target->Get4Momentum() << G4endl
905 << projectile->Get4Momentum() + target->Get4Momentum() << " "
906 << (projectile->Get4Momentum() + target->Get4Momentum()).mag() << G4endl;
907 #endif
908 if ( Result ) returnCode = 0; // successfully ended and nothing else needs to be done
909 return returnCode;
910 }
911
912 #ifdef debugFTFexictation
913 G4cout << "Make excitation of new hadrons" << G4endl;
914 #endif
915
916 // Redefinition of ProbOfDiffraction because the probabilities are changed due to quark exchange
917 common.ProbOfDiffraction = common.ProbProjectileDiffraction + common.ProbTargetDiffraction;
918 if ( common.ProbOfDiffraction != 0.0 ) {
919 common.ProbProjectileDiffraction /= common.ProbOfDiffraction;
920 common.ProbTargetDiffraction /= common.ProbOfDiffraction;
921 }
922
923 return returnCode = 1; // successfully completed, but the work needs to be continued
924}
G4long
long G4long
Definition: G4Types.hh:87
G4DiffractiveExcitation::NewNucleonId
G4int NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const
Definition: G4DiffractiveExcitation.cc:1537
G4DiffractiveExcitation::UnpackMeson
void UnpackMeson(G4int IdPDG, G4int &Q1, G4int &Q2) const
Definition: G4DiffractiveExcitation.cc:1496
G4DiffractiveExcitation::UnpackBaryon
void UnpackBaryon(G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const
Definition: G4DiffractiveExcitation.cc:1526
G4ElasticHNScattering::ElasticScattering
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
Definition: G4ElasticHNScattering.cc:64
G4FTFParameters::GetDeltaProbAtQuarkExchange
G4double GetDeltaProbAtQuarkExchange()
Definition: G4FTFParameters.hh:694
G4FTFParameters::GetProbOfSameQuarkExchange
G4double GetProbOfSameQuarkExchange()
Definition: G4FTFParameters.hh:698
G4ParticleDefinition::GetPDGWidth
G4double GetPDGWidth() const
Definition: G4ParticleDefinition.hh:112
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:586
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4VSplitableHadron::SetStatus
void SetStatus(const G4int aStatus)
Definition: G4VSplitableHadron.hh:167
G4VSplitableHadron::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4VSplitableHadron.hh:137
G4VSplitableHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4VSplitableHadron.hh:122
G4VSplitableHadron::SetDefinition
void SetDefinition(const G4ParticleDefinition *aDefinition)
Definition: G4VSplitableHadron.hh:132
G4VSplitableHadron::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4VSplitableHadron.hh:127

References common(), G4ElasticHNScattering::ElasticScattering(), G4ParticleTable::FindParticle(), G4cout, G4endl, G4UniformRand, G4VSplitableHadron::Get4Momentum(), G4VSplitableHadron::GetDefinition(), G4FTFParameters::GetDeltaProbAtQuarkExchange(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGiIsospin(), G4ParticleDefinition::GetPDGMass(), G4ParticleDefinition::GetPDGWidth(), G4FTFParameters::GetProbOfSameQuarkExchange(), G4VSplitableHadron::GetStatus(), MeV, NewNucleonId(), G4VSplitableHadron::Set4Momentum(), G4VSplitableHadron::SetDefinition(), G4VSplitableHadron::SetStatus(), sqr(), UnpackBaryon(), and UnpackMeson().

Referenced by ExciteParticipants().

◆ ExciteParticipants_doDiffraction()

G4bool G4DiffractiveExcitation::ExciteParticipants_doDiffraction ( G4VSplitableHadron *  projectile,
G4VSplitableHadron *  target,
G4FTFParameters *  theParameters,
G4DiffractiveExcitation::CommonVariables &  common 
) const
private

Definition at line 928 of file G4DiffractiveExcitation.cc.

931 {
932 // Second of the three utility methods used only by ExciteParticipants:
933 // it does the sampling for the diffraction case, either projectile or target diffraction.
934
935 G4bool isProjectileDiffraction = false;
936 if ( G4UniformRand() < common.ProbProjectileDiffraction ) { // projectile diffraction
937 isProjectileDiffraction = true;
938 #ifdef debugFTFexictation
939 G4cout << "projectile diffraction" << G4endl;
940 #endif
941 common.ProjMassT2 = common.ProjectileDiffStateMinMass2;
942 common.ProjMassT = common.ProjectileDiffStateMinMass;
943 common.TargMassT2 = common.M0target2;
944 common.TargMassT = common.M0target;
945 } else { // target diffraction
946 #ifdef debugFTFexictation
947 G4cout << "Target diffraction" << G4endl;
948 #endif
949 common.ProjMassT2 = common.M0projectile2;
950 common.ProjMassT = common.M0projectile;
951 common.TargMassT2 = common.TargetDiffStateMinMass2;
952 common.TargMassT = common.TargetDiffStateMinMass;
953 }
954
955 // Check whether the interaction is possible
956 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
957
958 common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
959 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
960 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
961 if ( common.PZcms2 < 0.0 ) return false;
962 common.maxPtSquare = common.PZcms2;
963
964 G4double DiffrAveragePt2 = theParameters->GetAvaragePt2ofElasticScattering()*1.2;
965 G4bool loopCondition = true;
966 G4int whilecount = 0;
967 do { // Generate pt and mass of projectile
968
969 whilecount++;
970 if ( whilecount > 1000 ) {
971 common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
972 return false; // Ignore this interaction
973 };
974
975 common.Qmomentum = G4LorentzVector( GaussianPt( DiffrAveragePt2, common.maxPtSquare ), 0 );
976 common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
977 if ( isProjectileDiffraction ) { // projectile diffraction
978 common.ProjMassT2 = common.ProjectileDiffStateMinMass2 + common.Pt2;
979 common.TargMassT2 = common.M0target2 + common.Pt2;
980 } else { // target diffraction
981 common.ProjMassT2 = common.M0projectile2 + common.Pt2;
982 common.TargMassT2 = common.TargetDiffStateMinMass2 + common.Pt2;
983 }
984 common.ProjMassT = std::sqrt( common.ProjMassT2 );
985 common.TargMassT = std::sqrt( common.TargMassT2 );
986 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
987
988 common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
989 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
990 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
991 if ( common.PZcms2 < 0.0 ) continue;
992
993 common.PZcms = std::sqrt( common.PZcms2 );
994 if ( isProjectileDiffraction ) { // projectile diffraction
995 common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
996 common.PMinusMax = common.SqrtS - common.TargMassT;
997 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
998 common.TMinusNew = common.SqrtS - common.PMinusNew;
999 common.Qminus = common.Ptarget.minus() - common.TMinusNew;
1000 common.TPlusNew = common.TargMassT2 / common.TMinusNew;
1001 common.Qplus = common.Ptarget.plus() - common.TPlusNew;
1002 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1003 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1004 loopCondition = ( ( common.Pprojectile + common.Qmomentum ).mag2() <
1005 common.ProjectileDiffStateMinMass2 );
1006 } else { // target diffraction
1007 common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
1008 common.TPlusMax = common.SqrtS - common.ProjMassT;
1009 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1010 common.PPlusNew = common.SqrtS - common.TPlusNew;
1011 common.Qplus = common.PPlusNew - common.Pprojectile.plus();
1012 common.PMinusNew = common.ProjMassT2 / common.PPlusNew;
1013 common.Qminus = common.PMinusNew - common.Pprojectile.minus();
1014 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1015 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1016 loopCondition = ( ( common.Ptarget - common.Qmomentum ).mag2() <
1017 common.TargetDiffStateMinMass2 );
1018 }
1019
1020 } while ( loopCondition ); /* Loop checking, 10.08.2015, A.Ribon */
1021 // Repeat the sampling because there was not any excitation
1022
1023 if ( isProjectileDiffraction ) { // projectile diffraction
1024 projectile->SetStatus( 0 );
1025 if ( projectile->GetStatus() == 2 ) projectile->SetStatus( 1 );
1026 if ( target->GetStatus() == 1 && target->GetSoftCollisionCount() == 0 ) target->SetStatus( 2 );
1027 } else { // target diffraction
1028 target->SetStatus( 0 );
1029 }
1030
1031 return true;
1032}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
CLHEP::Hep3Vector::mag2
double mag2() const
G4DiffractiveExcitation::ChooseP
G4double ChooseP(G4double Pmin, G4double Pmax) const
Definition: G4DiffractiveExcitation.cc:1444
G4DiffractiveExcitation::GaussianPt
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
Definition: G4DiffractiveExcitation.cc:1461
G4FTFParameters::GetAvaragePt2ofElasticScattering
G4double GetAvaragePt2ofElasticScattering()
Definition: G4FTFParameters.hh:688

References ChooseP(), common(), G4cout, G4endl, G4UniformRand, GaussianPt(), G4FTFParameters::GetAvaragePt2ofElasticScattering(), G4VSplitableHadron::GetSoftCollisionCount(), G4VSplitableHadron::GetStatus(), CLHEP::Hep3Vector::mag2(), G4VSplitableHadron::SetStatus(), and sqr().

Referenced by ExciteParticipants().

◆ ExciteParticipants_doNonDiffraction()

G4bool G4DiffractiveExcitation::ExciteParticipants_doNonDiffraction ( G4VSplitableHadron *  projectile,
G4VSplitableHadron *  target,
G4FTFParameters *  theParameters,
G4DiffractiveExcitation::CommonVariables &  common 
) const
private

Definition at line 1036 of file G4DiffractiveExcitation.cc.

1040 {
1041 // Third of the three utility methods used only by ExciteParticipants:
1042 // it does the sampling for the non-diffraction case.
1043
1044 #ifdef debugFTFexictation
1045 G4cout << "Non-diffraction process" << G4endl;
1046 #endif
1047
1048 // Check whether the interaction is possible
1049 common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2;
1050 common.ProjMassT = common.ProjectileNonDiffStateMinMass;
1051 common.TargMassT2 = common.TargetNonDiffStateMinMass2;
1052 common.TargMassT = common.TargetNonDiffStateMinMass;
1053 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) return false;
1054
1055 common.PZcms2 = ( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
1056 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1057 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1058 common.maxPtSquare = common.PZcms2;
1059
1060 G4int whilecount = 0;
1061 do { // Generate pt and masses
1062
1063 whilecount++;
1064 if ( whilecount > 1000 ) {
1065 common.Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1066 return false; // Ignore this interaction
1067 };
1068
1069 common.Qmomentum = G4LorentzVector( GaussianPt( theParameters->GetAveragePt2(),
1070 common.maxPtSquare ), 0 );
1071 common.Pt2 = G4ThreeVector( common.Qmomentum.vect() ).mag2();
1072 common.ProjMassT2 = common.ProjectileNonDiffStateMinMass2 + common.Pt2;
1073 common.ProjMassT = std::sqrt( common.ProjMassT2 );
1074 common.TargMassT2 = common.TargetNonDiffStateMinMass2 + common.Pt2;
1075 common.TargMassT = std::sqrt( common.TargMassT2 );
1076 if ( common.SqrtS < common.ProjMassT + common.TargMassT ) continue;
1077
1078 common.PZcms2 =( sqr( common.S ) + sqr( common.ProjMassT2 ) + sqr( common.TargMassT2 )
1079 - 2.0 * ( common.S * ( common.ProjMassT2 + common.TargMassT2 )
1080 + common.ProjMassT2 * common.TargMassT2 ) ) / 4.0 / common.S;
1081 if ( common.PZcms2 < 0.0 ) continue;
1082
1083 common.PZcms = std::sqrt( common.PZcms2 );
1084 common.PMinusMin = std::sqrt( common.ProjMassT2 + common.PZcms2 ) - common.PZcms;
1085 common.PMinusMax = common.SqrtS - common.TargMassT;
1086 common.TPlusMin = std::sqrt( common.TargMassT2 + common.PZcms2 ) - common.PZcms;
1087 common.TPlusMax = common.SqrtS - common.ProjMassT;
1088
1089 if ( G4UniformRand() <= 0.5 ) { // Random choice projectile or target sampling
1090 if ( G4UniformRand() < theParameters->GetProbLogDistrPrD() ) {
1091 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1092 } else {
1093 common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*G4UniformRand() + common.PMinusMin;
1094 }
1095 if ( G4UniformRand() < theParameters->GetProbLogDistr() ) {
1096 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1097 } else {
1098 common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*G4UniformRand() + common.TPlusMin;
1099 }
1100 } else {
1101 if ( G4UniformRand() < theParameters->GetProbLogDistr() ) {
1102 common.TPlusNew = ChooseP( common.TPlusMin, common.TPlusMax );
1103 } else {
1104 common.TPlusNew = ( common.TPlusMax - common.TPlusMin )*G4UniformRand() + common.TPlusMin;
1105 }
1106 if ( G4UniformRand() < theParameters->GetProbLogDistrPrD() ) {
1107 common.PMinusNew = ChooseP( common.PMinusMin, common.PMinusMax );
1108 } else {
1109 common.PMinusNew = ( common.PMinusMax - common.PMinusMin )*G4UniformRand() + common.PMinusMin;
1110 }
1111 }
1112
1113 common.Qminus = common.PMinusNew - common.Pprojectile.minus();
1114 common.Qplus = -( common.TPlusNew - common.Ptarget.plus() );
1115 common.Qmomentum.setPz( (common.Qplus - common.Qminus)/2.0 );
1116 common.Qmomentum.setE( (common.Qplus + common.Qminus)/2.0 );
1117 #ifdef debugFTFexictation
1118 G4cout <<"Sampled: Mpr, MdifPr, Mtr, MdifTr "<<G4endl
1119 << ( common.Pprojectile + common.Qmomentum ).mag() << " "
1120 << common.ProjectileNonDiffStateMinMass << G4endl
1121 << ( common.Ptarget - common.Qmomentum ).mag() << " "
1122 << common.TargetNonDiffStateMinMass << G4endl;
1123 #endif
1124
1125 } while ( ( common.Pprojectile + common.Qmomentum ).mag2() <
1126 common.ProjectileNonDiffStateMinMass2 || // No double Diffraction
1127 ( common.Ptarget - common.Qmomentum ).mag2() <
1128 common.TargetNonDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
1129
1130 projectile->SetStatus( 0 );
1131 target->SetStatus( 0 );
1132 return true;
1133}
G4FTFParameters::GetProbLogDistrPrD
G4double GetProbLogDistrPrD()
Definition: G4FTFParameters.hh:722
G4FTFParameters::GetAveragePt2
G4double GetAveragePt2()
Definition: G4FTFParameters.hh:718
G4FTFParameters::GetProbLogDistr
G4double GetProbLogDistr()
Definition: G4FTFParameters.hh:726

References ChooseP(), common(), G4cout, G4endl, G4UniformRand, GaussianPt(), G4FTFParameters::GetAveragePt2(), G4FTFParameters::GetProbLogDistr(), G4FTFParameters::GetProbLogDistrPrD(), CLHEP::Hep3Vector::mag2(), and sqr().

Referenced by ExciteParticipants().

◆ GaussianPt()

G4ThreeVector G4DiffractiveExcitation::GaussianPt ( G4double  AveragePt2,
G4double  maxPtSquare 
) const
private

Definition at line 1461 of file G4DiffractiveExcitation.cc.

1461 {
1462 // @@ this method is used in FTFModel as well. Should go somewhere common!
1463 G4double Pt2( 0.0 );
1464 if ( AveragePt2 <= 0.0 ) {
1465 Pt2 = 0.0;
1466 } else {
1467 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1468 ( G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1469 }
1470 G4double Pt = std::sqrt( Pt2 );
1471 G4double phi = G4UniformRand() * twopi;
1472 return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1473}

References G4Exp(), G4Log(), G4UniformRand, and twopi.

Referenced by ExciteParticipants_doDiffraction(), and ExciteParticipants_doNonDiffraction().

◆ GetQuarkFractionOfKink()

G4double G4DiffractiveExcitation::GetQuarkFractionOfKink ( G4double  zmin,
G4double  zmax 
) const
private

Definition at line 1478 of file G4DiffractiveExcitation.cc.

1478 {
1479 G4double z, yf;
1480 const G4int maxNumberOfLoops = 10000;
1481 G4int loopCounter = 0;
1482 do {
1483 z = zmin + G4UniformRand() * (zmax - zmin);
1484 yf = z*z + sqr(1.0 - z);
1485 } while ( ( G4UniformRand() > yf ) &&
1486 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1487 if ( loopCounter >= maxNumberOfLoops ) {
1488 z = 0.5*(zmin + zmax); // Just something acceptable, without any physics consideration.
1489 }
1490 return z;
1491}

References G4UniformRand, and sqr().

Referenced by CreateStrings().

◆ LambdaF()

G4double G4DiffractiveExcitation::LambdaF ( G4double  sqrM,
G4double  sqrM1,
G4double  sqrM2 
) const
private

◆ NewNucleonId()

G4int G4DiffractiveExcitation::NewNucleonId ( G4int  Q1,
G4int  Q2,
G4int  Q3 
) const
private

Definition at line 1537 of file G4DiffractiveExcitation.cc.

1537 {
1538 // Order the three integers in such a way that Q1 >= Q2 >= Q3
1539 G4int TmpQ( 0 );
1540 if ( Q3 > Q2 ) {
1541 TmpQ = Q2;
1542 Q2 = Q3;
1543 Q3 = TmpQ;
1544 } else if ( Q3 > Q1 ) {
1545 TmpQ = Q1;
1546 Q1 = Q3;
1547 Q3 = TmpQ;
1548 }
1549 if ( Q2 > Q1 ) {
1550 TmpQ = Q1;
1551 Q1 = Q2;
1552 Q2 = TmpQ;
1553 }
1554 // By now Q1 >= Q2 >= Q3
1555 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1556 return NewCode;
1557}

Referenced by ExciteParticipants_doChargeExchange().

◆ operator!=()

G4bool G4DiffractiveExcitation::operator!= ( const G4DiffractiveExcitation &  right) const
private

Definition at line 1587 of file G4DiffractiveExcitation.cc.

1587 {
1588 throw G4HadronicException( __FILE__, __LINE__,
1589 "G4DiffractiveExcitation != operator not meant to be called" );
1590}

◆ operator=()

const G4DiffractiveExcitation & G4DiffractiveExcitation::operator= ( const G4DiffractiveExcitation &  right)
private

Definition at line 1570 of file G4DiffractiveExcitation.cc.

1570 {
1571 throw G4HadronicException( __FILE__, __LINE__,
1572 "G4DiffractiveExcitation = operator not meant to be called" );
1573 return *this;
1574}

◆ operator==()

G4bool G4DiffractiveExcitation::operator== ( const G4DiffractiveExcitation &  right) const
private

Definition at line 1579 of file G4DiffractiveExcitation.cc.

1579 {
1580 throw G4HadronicException( __FILE__, __LINE__,
1581 "G4DiffractiveExcitation == operator not meant to be called" );
1582}

◆ UnpackBaryon()

void G4DiffractiveExcitation::UnpackBaryon ( G4int  IdPDG,
G4int &  Q1,
G4int &  Q2,
G4int &  Q3 
) const
private

Definition at line 1526 of file G4DiffractiveExcitation.cc.

1527 {
1528 Q1 = IdPDG / 1000;
1529 Q2 = (IdPDG % 1000) / 100;
1530 Q3 = (IdPDG % 100) / 10;
1531 return;
1532}

Referenced by ExciteParticipants_doChargeExchange().

◆ UnpackMeson()

void G4DiffractiveExcitation::UnpackMeson ( G4int  IdPDG,
G4int &  Q1,
G4int &  Q2 
) const
private

Definition at line 1496 of file G4DiffractiveExcitation.cc.

1496 {
1497 G4int absIdPDG = std::abs( IdPDG );
1498 if ( ! ( absIdPDG == 111 || absIdPDG == 221 || absIdPDG == 331 || // Pi0 , Eta , Eta'
1499 absIdPDG == 441 || absIdPDG == 443 || absIdPDG == 553 ) ) { // Etac , J/psi , Upsilon
1500 // All other projectile mesons (including charmed and bottom ones)
1501 Q1 = absIdPDG / 100;
1502 Q2 = (absIdPDG % 100) / 10;
1503 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1504 if ( IdPDG < 0 ) anti *= -1;
1505 Q1 *= anti;
1506 Q2 *= -1 * anti;
1507 } else {
1508 if ( absIdPDG == 441 || absIdPDG == 443 ) { // Etac , J/psi
1509 Q1 = 4; Q2 = -4;
1510 } else if ( absIdPDG == 553 ) { // Upsilon
1511 Q1 = 5; Q2 = -5;
1512 } else { // Pi0 , Eta , Eta'
1513 if ( G4UniformRand() < 0.5 ) {
1514 Q1 = 1; Q2 = -1;
1515 } else {
1516 Q1 = 2; Q2 = -2;
1517 }
1518 }
1519 }
1520 return;
1521}
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112

References G4UniformRand, and G4INCL::Math::max().

Referenced by ExciteParticipants_doChargeExchange().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.3