Geant4  10.02.p03
Public Member Functions | Private Member Functions | List of all members
G4DiffractiveExcitation Class Reference

#include <G4DiffractiveExcitation.hh>

Collaboration diagram for G4DiffractiveExcitation:

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
 

Private Member Functions

 G4DiffractiveExcitation (const G4DiffractiveExcitation &right)
 
const G4DiffractiveExcitation & operator= (const G4DiffractiveExcitation &right)
 
int operator== (const G4DiffractiveExcitation &right) const
 
int operator!= (const G4DiffractiveExcitation &right) const
 
G4double LambdaF (G4double sqrM, G4double sqrM1, G4double sqrM2) const
 
G4ThreeVector GaussianPt (G4double AveragePt2, G4double maxPtSquare) const
 
G4double ChooseP (G4double Pmin, G4double Pmax) const
 
G4double GetQuarkFractionOfKink (G4double zmin, G4double zmax) const
 
void UnpackMeson (G4int IdPDG, G4int &Q1, G4int &Q2) const
 
void UnpackBaryon (G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const
 
G4int NewNucleonId (G4int Q1, G4int Q2, G4int Q3) const
 

Detailed Description

Definition at line 51 of file G4DiffractiveExcitation.hh.

Constructor & Destructor Documentation

◆ G4DiffractiveExcitation() [1/2]

G4DiffractiveExcitation::G4DiffractiveExcitation ( )

Definition at line 84 of file G4DiffractiveExcitation.cc.

84 {}

◆ ~G4DiffractiveExcitation()

G4DiffractiveExcitation::~G4DiffractiveExcitation ( )
virtual

Definition at line 89 of file G4DiffractiveExcitation.cc.

89 {}

◆ G4DiffractiveExcitation() [2/2]

G4DiffractiveExcitation::G4DiffractiveExcitation ( const G4DiffractiveExcitation &  right)
private

Definition at line 1640 of file G4DiffractiveExcitation.cc.

1640  {
1641  throw G4HadronicException( __FILE__, __LINE__,
1642  "G4DiffractiveExcitation copy contructor not meant to be called" );
1643 }

Member Function Documentation

◆ ChooseP()

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

Definition at line 1532 of file G4DiffractiveExcitation.cc.

1532  {
1533  // Choose an x between Xmin and Xmax with P(x) ~ 1/x .
1534  // To be improved...
1535  G4double range = Pmax - Pmin;
1536  if ( Pmin <= 0.0 || range <= 0.0 ) {
1537  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1538  throw G4HadronicException( __FILE__, __LINE__,
1539  "G4DiffractiveExcitation::ChooseP : Invalid arguments " );
1540  }
1541  G4double P = Pmin * G4Pow::GetInstance()->powA( Pmax/Pmin, G4UniformRand() );
1542  //G4double P = (Pmax - Pmin) * G4UniformRand() + Pmin;
1543  return P;
1544 }
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
P
static double P[]
Definition: nf_gammaFunctions.cc:89
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CreateStrings()

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

Definition at line 1182 of file G4DiffractiveExcitation.cc.

1186  {
1187 
1188  //G4cout << "Create Strings SplitUp " << hadron << G4endl
1189  // << "Defin " << hadron->GetDefinition() << G4endl
1190  // << "Defin " << hadron->GetDefinition()->GetPDGEncoding() << G4endl;
1191 
1192  hadron->SplitUp();
1193 
1194  G4Parton* start = hadron->GetNextParton();
1195  if ( start == NULL ) {
1196  G4cout << " G4FTFModel::String() Error: No start parton found" << G4endl;
1197  FirstString = 0; SecondString = 0;
1198  return;
1199  }
1200 
1201  G4Parton* end = hadron->GetNextParton();
1202  if ( end == NULL ) {
1203  G4cout << " G4FTFModel::String() Error: No end parton found" << G4endl;
1204  FirstString = 0; SecondString = 0;
1205  return;
1206  }
1207 
1208  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1209  // << G4endl
1210  // << "Create string " << start->GetPDGcode() << " " << end->GetPDGcode() << G4endl;
1211 
1212  G4LorentzVector Phadron = hadron->Get4Momentum();
1213  //G4cout << "String mom " << Phadron << G4endl;
1214  G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 );
1215  G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 );
1216  G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 );
1217  G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 );
1218  G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 );
1219 
1220  G4int PDGcode_startQ = std::abs( start->GetDefinition()->GetPDGEncoding() );
1221  G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding() );
1222  //G4cout << "PDGcode_startQ " << PDGcode_startQ << " PDGcode_endQ " << PDGcode_endQ << G4endl;
1223 
1224  G4double Wmin( 0.0 );
1225  if ( isProjectile ) {
1226  Wmin = theParameters->GetProjMinDiffMass();
1227  } else {
1228  Wmin = theParameters->GetTarMinDiffMass();
1229  }
1230 
1231  G4double W = hadron->Get4Momentum().mag();
1232  //G4cout << "Wmin W " << Wmin << " " << W << G4endl;
1233  //G4int Uzhi; G4cin >> Uzhi;
1234  G4double W2 = W*W;
1235  G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); // x2( 0.0 )
1236  G4bool Kink = false;
1237 
1238  if ( ! ( ( start->GetDefinition()->GetParticleSubType() == "di_quark" &&
1239  end->GetDefinition()->GetParticleSubType() == "di_quark" ) ||
1240  ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1241  end->GetDefinition()->GetParticleSubType() == "quark" ) ) ) {
1242  // Kinky strings are allowed only for qq-q strings;
1243  // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
1244  // according to the analysis of Pbar P interactions
1245 
1246  if ( W > Wmin ) { // Kink is possible
1247  if ( hadron->GetStatus() == 0 ) { // VU 10.04.2012
1248  G4double Pt2kink = theParameters->GetPt2Kink(); // For non-diffractive
1249 // Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) ); // Uzhi 18 Sept. 2014
1250  if(Pt2kink) // Uzhi 18 Sept. 2014
1251  {Pt = std::sqrt( Pt2kink * ( G4Pow::GetInstance()->powA( W2/16.0/Pt2kink + 1.0, G4UniformRand() ) - 1.0 ) );} // Uzhi 18 Sept. 2014
1252  else {Pt=0.;} // Uzhi 18 Sept. 2014
1253  } else {
1254  Pt = 0.0;
1255  }
1256 
1257  if ( Pt > 500.0*MeV ) {
1258  G4double Ymax = G4Log( W/2.0/Pt + std::sqrt( W2/4.0/Pt/Pt - 1.0 ) );
1259  G4double Y = Ymax*( 1.0 - 2.0*G4UniformRand() );
1260  x1 = 1.0 - Pt/W * G4Exp( Y );
1261  x3 = 1.0 - Pt/W * G4Exp(-Y );
1262  //x2 = 2.0 - x1 - x3;
1263 
1264  G4double Mass_startQ = 650.0*MeV;
1265  if ( PDGcode_startQ < 3 ) Mass_startQ = 325.0*MeV;
1266  if ( PDGcode_startQ == 3 ) Mass_startQ = 500.0*MeV;
1267  if ( PDGcode_startQ == 4 ) Mass_startQ = 1600.0*MeV;
1268  G4double Mass_endQ = 650.0*MeV;
1269  if ( PDGcode_endQ < 3 ) Mass_endQ = 325.0*MeV;
1270  if ( PDGcode_endQ == 3 ) Mass_endQ = 500.0*MeV;
1271  if ( PDGcode_endQ == 4 ) Mass_endQ = 1600.0*MeV;
1272 
1273  G4double P2_1 = W2*x1*x1/4.0 - Mass_endQ*Mass_endQ;
1274  G4double P2_3 = W2*x3*x3/4.0 - Mass_startQ*Mass_startQ;
1275  G4double P2_2 = sqr( (2.0 - x1 - x3)*W/2.0 );
1276  if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) {
1277  Kink = false;
1278  } else {
1279  G4double P_1 = std::sqrt( P2_1 );
1280  G4double P_2 = std::sqrt( P2_2 );
1281  G4double P_3 = std::sqrt( P2_3 );
1282  G4double CosT12 = ( P2_3 - P2_1 - P2_2 ) / (2.0*P_1*P_2);
1283  G4double CosT13 = ( P2_2 - P2_1 - P2_3 ) / (2.0*P_1*P_3);
1284  //Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated 11.12.09
1285 
1286  if ( std::abs( CosT12 ) > 1.0 || std::abs( CosT13 ) > 1.0 ) {
1287  Kink = false;
1288  } else {
1289  Kink = true;
1290  Pt = P_2 * std::sqrt( 1.0 - CosT12*CosT12 ); // because system was rotated 11.12.09
1291  Pstart.setPx( -Pt ); Pstart.setPy( 0.0 ); Pstart.setPz( P_3*CosT13 );
1292  Pend.setPx( 0.0 ); Pend.setPy( 0.0 ); Pend.setPz( P_1 );
1293  Pkink.setPx( Pt ); Pkink.setPy( 0.0 ); Pkink.setPz( P_2*CosT12 );
1294  Pstart.setE( x3*W/2.0 );
1295  Pkink.setE( Pkink.vect().mag() );
1296  Pend.setE( x1*W/2.0 );
1297 
1298  G4double XkQ = GetQuarkFractionOfKink( 0.0, 1.0 );
1299  if ( Pkink.getZ() > 0.0 ) {
1300  if ( XkQ > 0.5 ) {
1301  PkinkQ1 = XkQ*Pkink;
1302  } else {
1303  PkinkQ1 = (1.0 - XkQ)*Pkink;
1304  }
1305  } else {
1306  if ( XkQ > 0.5 ) {
1307  PkinkQ1 = (1.0 - XkQ)*Pkink;
1308  } else {
1309  PkinkQ1 = XkQ*Pkink;
1310  }
1311  }
1312 
1313  PkinkQ2 = Pkink - PkinkQ1;
1314  // Minimizing Pt1^2+Pt3^2
1315  G4double Cos2Psi = ( sqr(x1) - sqr(x3) + 2.0*sqr( x3*CosT13 ) ) /
1316  std::sqrt( sqr( sqr(x1) - sqr(x3) ) + sqr( 2.0*x1*x3*CosT13 ) );
1317  G4double Psi = std::acos( Cos2Psi );
1318 
1319  G4LorentzRotation Rotate;
1320  if ( isProjectile ) {
1321  Rotate.rotateY( Psi );
1322  } else {
1323  Rotate.rotateY( pi + Psi ); // Uzhi 18 Sept. 2014
1324 // Rotate.rotateY( pi - Psi ); // Uzhi 18 Sept. 2014
1325  }
1326  Rotate.rotateZ( twopi * G4UniformRand() );
1327  Pstart *= Rotate;
1328  Pkink *= Rotate;
1329  PkinkQ1 *= Rotate;
1330  PkinkQ2 *= Rotate;
1331  Pend *= Rotate;
1332  }
1333  } // End of if ( P2_1 <= 0.0 || P2_3 <= 0.0 )
1334  } // End of if ( Pt > 500.0*MeV )
1335  } // End of if ( W > Wmin ) : check for a kink
1336  } // end of qq-q string selection
1337 
1338  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetParticleSubType() << " "
1339  // << end->GetDefinition()->GetParticleSubType() << G4endl;
1340  //G4cout << "Kink " << Kink << " " << start->GetDefinition()->GetPDGEncoding() << " "
1341  // << end->GetDefinition()->GetPDGEncoding() << G4endl;
1342  //G4int Uzhi; G4cin >> Uzhi;
1343 
1344  if ( Kink ) { // Kink is possible
1345 
1346  //G4cout << "Kink is sampled!" << G4endl;
1347  std::vector< G4double > QuarkProbabilitiesAtGluonSplitUp =
1348  theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1349 
1350  G4int QuarkInGluon( 1 ); G4double Ksi = G4UniformRand();
1351  for ( unsigned int Iq = 0; Iq < 3; Iq++ ) {
1352  //G4cout << "Iq " << Iq << G4endl;
1353  if ( Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq] ) QuarkInGluon++;
1354  }
1355  //G4cout << "Last Iq " << QuarkInGluon << G4endl;
1356  G4Parton* Gquark = new G4Parton( QuarkInGluon );
1357  G4Parton* Ganti_quark = new G4Parton( -QuarkInGluon );
1358  //G4cout << "Lorentz " << G4endl;
1359 
1360  G4LorentzRotation toCMS( -1 * Phadron.boostVector() );
1361  G4LorentzRotation toLab( toCMS.inverse() );
1362  //G4cout << "Pstart " << Pstart << G4endl;
1363 //G4cout << "Pend " << Pend << G4endl;
1364 //G4cout << "Kink1 " <<PkinkQ1<<G4endl;
1365 //G4cout << "Kink2 " <<PkinkQ2<<G4endl;
1366 //G4cout << "Pstart " << Pstart << G4endl<<G4endl;
1367 
1368  Pstart.transform( toLab ); start->Set4Momentum( Pstart );
1369  PkinkQ1.transform( toLab );
1370  PkinkQ2.transform( toLab );
1371  Pend.transform( toLab ); end->Set4Momentum( Pend );
1372  //G4cout << "Pstart " << Pstart << G4endl;
1373  //G4cout << "Pend " << Pend << G4endl;
1374  //G4cout << "Defin " << hadron->GetDefinition()<< G4endl;
1375  //G4cout << "Defin " << hadron->GetDefinition()->GetPDGEncoding()<< G4endl;
1376 
1377  //G4int absPDGcode = std::abs( hadron->GetDefinition()->GetPDGEncoding() );
1378  G4int absPDGcode = 1500; // 23 Dec
1379  if ( start->GetDefinition()->GetParticleSubType() == "quark" &&
1380  end->GetDefinition()->GetParticleSubType() == "quark" ) {
1381  absPDGcode = 110;
1382  }
1383  //G4cout << "absPDGcode " << absPDGcode << G4endl;
1384 
1385  if ( absPDGcode < 1000 ) { // meson
1386  if ( isProjectile ) { // Projectile
1387  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1388 
1389 FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1390  SecondString = new G4ExcitedString( Gquark, start , +1 );
1391  Ganti_quark->Set4Momentum( PkinkQ1 );
1392  Gquark->Set4Momentum( PkinkQ2 );
1393  } else { // Anti_Quark on the end
1394  FirstString = new G4ExcitedString( end , Gquark, +1 );
1395  SecondString = new G4ExcitedString( Ganti_quark, start , +1 );
1396  Gquark->Set4Momentum( PkinkQ1 );
1397  Ganti_quark->Set4Momentum( PkinkQ2 );
1398  }
1399  } else { // Target
1400  if ( end->GetDefinition()->GetPDGEncoding() > 0 ) { // A quark on the end
1401  FirstString = new G4ExcitedString( Ganti_quark, end , -1 );
1402  SecondString = new G4ExcitedString( start , Gquark, -1 );
1403  Ganti_quark->Set4Momentum( PkinkQ2 );
1404  Gquark->Set4Momentum( PkinkQ1 );
1405  } else { // Anti_Quark on the end
1406  FirstString = new G4ExcitedString( Gquark, end , -1 );
1407  SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1408  Gquark->Set4Momentum( PkinkQ2 );
1409  Ganti_quark->Set4Momentum( PkinkQ1 );
1410  }
1411  }
1412  } else { // Baryon/AntiBaryon
1413 
1414 //G4cout<<"isProjectile "<<isProjectile<<G4endl;
1415 //G4cout<<"end "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1416 //G4cout<<PkinkQ1<<G4endl;
1417 //G4cout<<PkinkQ2<<G4endl;
1418 //G4cout<<"start "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1419 //G4int Uzhi; G4cin>>Uzhi;
1420  if ( isProjectile ) { // Projectile
1421  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1422  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1423  FirstString = new G4ExcitedString( end , Gquark, +1 ); // Open Uzhi
1424  SecondString = new G4ExcitedString( Ganti_quark, start , +1 ); // Open Uzhi
1425  Gquark->Set4Momentum( PkinkQ1 );
1426  Ganti_quark->Set4Momentum( PkinkQ2 );
1427 
1428 // FirstString = new G4ExcitedString( Gquark, end, +1 ); // Uzhi 18 Sept. 2014
1429 // SecondString = new G4ExcitedString( start, Ganti_quark, +1 ); // Uzhi 18 Sept. 2014
1430 
1431  } else { // Anti_DiQuark on the end or quark
1432  FirstString = new G4ExcitedString( end , Ganti_quark, +1 );
1433  SecondString = new G4ExcitedString( Gquark, start , +1 );
1434  Ganti_quark->Set4Momentum( PkinkQ1 );
1435  Gquark->Set4Momentum( PkinkQ2 );
1436  }
1437  } else { // Target
1438  if ( end->GetDefinition()->GetParticleType() == "diquarks" &&
1439  end->GetDefinition()->GetPDGEncoding() > 0 ) { // DiQuark on the end
1440 // FirstString = new G4ExcitedString( Gquark, end , -1 );
1441 // SecondString = new G4ExcitedString( start , Ganti_quark, -1 );
1442  Gquark->Set4Momentum( PkinkQ1 );
1443  Ganti_quark->Set4Momentum( PkinkQ2 );
1444 
1445  FirstString = new G4ExcitedString( end, Gquark, -1 );
1446  SecondString = new G4ExcitedString( Ganti_quark, start, -1 );
1447 
1448  } else { // Anti_DiQuark on the end or Q
1449  FirstString = new G4ExcitedString( Ganti_quark, end , -1 ); // Uzhi ?????
1450  SecondString = new G4ExcitedString( start , Gquark, -1 );
1451  Gquark->Set4Momentum( PkinkQ2 );
1452  Ganti_quark->Set4Momentum( PkinkQ1 );
1453  }
1454  }
1455  }
1456 
1457  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1458  FirstString->SetPosition( hadron->GetPosition() );
1459  SecondString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1460  SecondString->SetPosition( hadron->GetPosition() );
1461 
1462  } else { // End of kink is possible: Kink is impossible
1463 
1464  //G4cout << start << " " << start->GetPDGcode() << " " << end << " " << end->GetPDGcode()
1465  // << G4endl;
1466 /* // Uzhi 18 Sept. 2014
1467  if ( isProjectile ) {
1468  FirstString = new G4ExcitedString( end, start, +1 );
1469  } else {
1470  FirstString = new G4ExcitedString( start, end, -1 );
1471  }
1472 */ // Uzhi 18 Sept. 2014
1473 
1474  FirstString = new G4ExcitedString( end, start, +1 ); // Uzhi 18 Sept. 2014
1475 
1476  FirstString->SetTimeOfCreation( hadron->GetTimeOfCreation() );
1477  FirstString->SetPosition( hadron->GetPosition() );
1478  SecondString = 0;
1479 
1480  // momenta of string ends
1481  G4double Momentum = hadron->Get4Momentum().vect().mag();
1482  G4double Plus = hadron->Get4Momentum().e() + Momentum;
1483  G4double Minus = hadron->Get4Momentum().e() - Momentum;
1484  G4ThreeVector tmp;
1485  if ( Momentum > 0.0 ) {
1486  tmp.set( hadron->Get4Momentum().px(),
1487  hadron->Get4Momentum().py(),
1488  hadron->Get4Momentum().pz() );
1489  tmp /= Momentum;
1490  } else {
1491  tmp.set( 0.0, 0.0, 1.0 );
1492  }
1493  G4LorentzVector Pstart1( tmp, 0.0 );
1494  G4LorentzVector Pend1( tmp, 0.0 );
1495  if ( isProjectile ) {
1496  Pstart1 *= (-1.0)*Minus/2.0;
1497  Pend1 *= (+1.0)*Plus /2.0;
1498  } else {
1499  Pstart1 *= (+1.0)*Plus/ 2.0;
1500  Pend1 *= (-1.0)*Minus/2.0;
1501  }
1502  Momentum = -Pstart1.mag();
1503  Pstart1.setT( Momentum ); // It is assumed that quark has m=0.
1504  Momentum = -Pend1.mag();
1505  Pend1.setT( Momentum ); // It is assumed that di-quark has m=0.
1506  start->Set4Momentum( Pstart1 );
1507  end->Set4Momentum( Pend1 );
1508  SecondString = 0;
1509 
1510  } // End of kink is impossible
1511 
1512  //G4cout << "Quarks in the string at creation" << FirstString->GetRightParton()->GetPDGcode()
1513  // << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl
1514  // << FirstString << " " << SecondString << G4endl;
1515 
1516  #ifdef G4_FTFDEBUG
1517  G4cout << " generated string flavors " << start->GetPDGcode() << " / "
1518  << end->GetPDGcode() << G4endl << " generated string momenta: quark "
1519  << start->Get4Momentum() << "mass : " << start->Get4Momentum().mag() << G4endl
1520  << " generated string momenta: Diquark " << end->Get4Momentum() << "mass : "
1521  << end->Get4Momentum().mag() << G4endl << " sum of ends "
1522  << Pstart + Pend << G4endl << " Original "
1523  << hadron->Get4Momentum() << G4endl;
1524  #endif
1525 
1526  return;
1527 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
MeV
static const double MeV
Definition: G4SIunits.hh:211
tmp
Float_t tmp
Definition: advanced/microbeam/plot.C:37
G4Parton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:140
G4FTFParameters::GetProjMinDiffMass
G4double GetProjMinDiffMass()
Definition: G4FTFParameters.hh:431
G4FTFParameters::GetTarMinDiffMass
G4double GetTarMinDiffMass()
Definition: G4FTFParameters.hh:439
Y
Float_t Y
Definition: advanced/microbeam/plot.C:39
G4Parton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:145
G4Parton::GetDefinition
G4ParticleDefinition * GetDefinition()
Definition: G4Parton.hh:158
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:141
G4FTFParameters::GetPt2Kink
G4double GetPt2Kink()
Definition: G4FTFParameters.hh:461
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:162
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4ExcitedString::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4ExcitedString.hh:161
G4ExcitedString::SetPosition
void SetPosition(const G4ThreeVector &aPosition)
Definition: G4ExcitedString.hh:155
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4ExcitedString::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition: G4ExcitedString.hh:143
G4bool
bool G4bool
Definition: G4Types.hh:79
CLHEP::Hep3Vector::mag
double mag() const
twopi
static const double twopi
Definition: G4SIunits.hh:75
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
CLHEP::HepLorentzRotation::transform
HepLorentzRotation & transform(const HepBoost &b)
G4DiffractiveExcitation::GetQuarkFractionOfKink
G4double GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
Definition: G4DiffractiveExcitation.cc:1566
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:174
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
pi
static const double pi
Definition: G4SIunits.hh:74
G4Parton::GetPDGcode
G4int GetPDGcode() const
Definition: G4Parton.hh:124
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
G4endl
#define G4endl
Definition: G4ios.hh:61
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition: G4ParticleDefinition.hh:142
Ymax
int Ymax
Definition: LoadPlotNtuple.C:18
G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp
std::vector< G4double > GetQuarkProbabilitiesAtGluonSplitUp()
Definition: G4FTFParameters.hh:465
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ExciteParticipants()

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

Definition at line 94 of file G4DiffractiveExcitation.cc.

97  {
98 
99  #ifdef debugFTFexictation
100  G4cout << G4endl << "FTF ExciteParticipants --------------" << G4endl;
101  #endif
102 
103  // Projectile parameters
104  G4LorentzVector Pprojectile = projectile->Get4Momentum();
105  if ( Pprojectile.z() < 0.0 ) return false;
106 
107  G4int ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding();
108  G4int absProjectilePDGcode = std::abs( ProjectilePDGcode );
109  G4double M0projectile = Pprojectile.mag();
110  G4double ProjectileRapidity = Pprojectile.rapidity();
111 
112  // Target parameters
113  G4LorentzVector Ptarget = target->Get4Momentum();
114  G4int TargetPDGcode = target->GetDefinition()->GetPDGEncoding();
115  G4int absTargetPDGcode = std::abs( TargetPDGcode );
116  G4double M0target = Ptarget.mag();
117  //G4double TargetRapidity = Ptarget.rapidity();
118 
119  // Kinematical properties of the interactions
120  G4LorentzVector Psum = Pprojectile + Ptarget; // 4-momentum in CMS
121  G4double S = Psum.mag2();
122  G4double SqrtS = std::sqrt( S );
123  //Uzhi_SqrtS = std::sqrt( S );
124 
125  // Check off-shellness of the participants
126  G4SampleResonance BrW; // Uzhi Oct. 2014
127 
128  G4bool PutOnMassShell( false );
129 
130  G4double MminProjectile(0.); // Uzhi Oct. 2014
131  MminProjectile = BrW.GetMinimumMass(projectile->GetDefinition()); // Uzhi Oct. 2014
132 //G4double M0projectile = MminProjectile; // With de-excitation Uzhi Oct. 2014
133 //G4double M0projectile = Pprojectile.mag(); // Without de-excitation, see above
134 
135  if ( M0projectile < MminProjectile )
136  {
137  PutOnMassShell = true;
138  M0projectile = BrW.SampleMass(projectile->GetDefinition(),projectile->GetDefinition()->GetPDGMass() + 5.0*projectile->GetDefinition()->GetPDGWidth());
139  }
140 
141  G4double M0projectile2 = M0projectile * M0projectile;
142  G4double ProjectileDiffStateMinMass = theParameters->GetProjMinDiffMass(); // Uzhi Oct 2014
143  G4double ProjectileNonDiffStateMinMass = theParameters->GetProjMinNonDiffMass(); // Uzhi Oct 2014
144  if ( M0projectile > ProjectileDiffStateMinMass ) { // Uzhi Oct 2014
145  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV;
146  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV;
147  if(absProjectilePDGcode > 3000) { // Strange baryon // Uzhi Nov. 2014
148  ProjectileDiffStateMinMass += 140.0*MeV; // Uzhi Nov. 2014
149  ProjectileNonDiffStateMinMass += 140.0*MeV; // Uzhi Nov. 2014
150  } // Uzhi Nov. 2014
151  }
152 
153  G4double MminTarget(0.); // Uzhi Oct. 2014
154  MminTarget = BrW.GetMinimumMass(target->GetDefinition()); // Uzhi Oct. 2014
155  if ( M0target < MminTarget )
156  {
157  PutOnMassShell = true;
158  M0target = BrW.SampleMass(target->GetDefinition(),target->GetDefinition()->GetPDGMass() + 5.0*target->GetDefinition()->GetPDGWidth());
159  }
160 
161  G4double M0target2 = M0target * M0target;
162  G4double TargetDiffStateMinMass = theParameters->GetTarMinDiffMass(); // Uzhi Oct 2014
163  G4double TargetNonDiffStateMinMass = theParameters->GetTarMinNonDiffMass(); // Uzhi Oct 2014
164  if ( M0target > TargetDiffStateMinMass ) { // Uzhi Oct 2014
165  TargetDiffStateMinMass = M0target + 220.0*MeV;
166  TargetNonDiffStateMinMass = M0target + 220.0*MeV;
167  if(absTargetPDGcode > 3000) { // Strange baryon // Uzhi Nov. 2014
168  TargetDiffStateMinMass += 140.0*MeV; // Uzhi Nov. 2014
169  TargetNonDiffStateMinMass += 140.0*MeV; // Uzhi Nov. 2014
170  } // Uzhi Nov. 2014
171  };
172 
173  #ifdef debugFTFexictation
174  G4cout << "Proj Targ PDGcodes " << ProjectilePDGcode << " " << TargetPDGcode << G4endl
175  << "M0projectile Y " << M0projectile << " " << ProjectileRapidity << G4endl;
176  //G4cout << "M0target Y " << M0target << " " << TargetRapidity << G4endl;
177  G4cout << "Pproj " << Pprojectile << G4endl << "Ptarget " << Ptarget << G4endl;
178  #endif
179 
180  G4double AveragePt2 = theParameters->GetAveragePt2();
181 // G4double ProbLogDistrPrD = theParameters->GetProbLogDistrPrD(); // Uzhi Oct 2014 ***
182  G4double ProbLogDistr = theParameters->GetProbLogDistr();
183  G4double SumMasses = M0projectile + M0target; // + 220.0*MeV; // Uzhi Nov. 2014
184 
185  // Transform momenta to cms and then rotate parallel to z axis;
186  G4LorentzRotation toCms( -1 * Psum.boostVector() );
187 
188  G4LorentzVector Ptmp = toCms * Pprojectile;
189  if ( Ptmp.pz() <= 0.0 ) return false; // "String" moving backwards in CMS, abort collision!
190 
191  toCms.rotateZ( -1*Ptmp.phi() );
192  toCms.rotateY( -1*Ptmp.theta() );
193  G4LorentzRotation toLab(toCms.inverse());
194  Pprojectile.transform( toCms );
195  Ptarget.transform( toCms );
196 
197  G4double PZcms2, PZcms;
198 
199  #ifdef debugFTFexictation
200  G4cout << "SqrtS " << SqrtS << G4endl << "M0pr M0tr SumM+220 " << M0projectile << " "
201  << M0target << " " << SumMasses << G4endl;
202  #endif
203 
204  if ( SqrtS < M0projectile + M0target ) return false;
205  if ( SqrtS < SumMasses ) return false;
206  // The model cannot work at low energy
207 
208  PZcms2 = ( S*S + M0projectile2*M0projectile2 + M0target2*M0target2
209  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
210 
211  #ifdef debugFTFexictation
212  G4cout << "PZcms2 after PutOnMassShell " << PZcms2 << G4endl;
213  #endif
214 
215  if ( PZcms2 < 0 ) return false;
216  // It can be in an interaction with off-shell nuclear nucleon
217 
218  PZcms = std::sqrt( PZcms2 );
219  if ( PutOnMassShell ) {
220  if ( Pprojectile.z() > 0.0 ) {
221  Pprojectile.setPz( PZcms );
222  Ptarget.setPz( -PZcms );
223  } else {
224  Pprojectile.setPz( -PZcms );
225  Ptarget.setPz( PZcms );
226  };
227  Pprojectile.setE( std::sqrt( M0projectile2 +
228  Pprojectile.x()*Pprojectile.x() +
229  Pprojectile.y()*Pprojectile.y() +
230  PZcms2 ) );
231  Ptarget.setE( std::sqrt( M0target2 +
232  Ptarget.x()*Ptarget.x() +
233  Ptarget.y()*Ptarget.y() +
234  PZcms2 ) );
235  }
236 
237  G4double maxPtSquare(0.); // = PZcms2;
238 
239  //Uzhi_QEnex = 0;
240  //Uzhi_QEexc = 0;
241  //Uzhi_targetdiffraction = 0;
242  //Uzhi_projectilediffraction = 0;
243  //Uzhi_nondiffraction = 0;
244  //G4int UzhiPrD( 0 ), UzhiTrD( 0 ), UzhiND( 0 );
245 
246  #ifdef debugFTFexictation
247  G4cout << "Start --------------------" << G4endl << "Proj M0 Mdif Mndif " << M0projectile
248  << " " << ProjectileDiffStateMinMass << " " << ProjectileNonDiffStateMinMass << G4endl
249  << "Targ M0 Mdif Mndif " << M0target << " " << TargetDiffStateMinMass << " "
250  << TargetNonDiffStateMinMass << G4endl << "SqrtS " << SqrtS << G4endl
251  << "Proj CMS " << Pprojectile << G4endl << "Targ CMS " << Ptarget << G4endl;
252  #endif
253 
254  // Charge exchange can be possible
255  // Getting the values needed for exchange
256  // Check for possible quark exchange
257  G4double QeNoExc = theParameters->GetProcProb( 0, ProjectileRapidity );
258  G4double QeExc = theParameters->GetProcProb( 1, ProjectileRapidity )*theParameters->GetProcProb( 4, ProjectileRapidity );
259  G4double ProbProjectileDiffraction = theParameters->GetProcProb( 2, ProjectileRapidity );
260  G4double ProbTargetDiffraction = theParameters->GetProcProb( 3, ProjectileRapidity );
261 
262 // =================================== Apr. 19 ========================================= Vova
263 //QeNoExc = 1.;
264 //QeExc = 0.;
265 //ProbProjectileDiffraction = 0.0;
266 //ProbTargetDiffraction = 0.0;
267 //AveragePt2 = 0.;
268 // =================================== Apr. 19 =========================================
269 
270  if(QeNoExc+QeExc+ProbProjectileDiffraction+ProbTargetDiffraction > 1.) // Uzhi Nov. 2014
271  {QeNoExc=1.0-QeExc-ProbProjectileDiffraction-ProbTargetDiffraction;}
272 //QeNoExc=0.;
273  G4double ProbExc( 0.0 );
274  if ( QeExc + QeNoExc != 0.0 ) ProbExc = QeExc/(QeExc + QeNoExc);
275  G4double DeltaProbAtQuarkExchange = theParameters->GetDeltaProbAtQuarkExchange();
276  G4double DeltaMass = G4ParticleTable::GetParticleTable()->FindParticle( 2224 )->GetPDGMass();
277 
278 //ProbProjectileDiffraction = 0.5; // Vova
279 //ProbTargetDiffraction = 0.5; // Vova
280  G4double ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction;
281 
282  #ifdef debugFTFexictation
283  G4cout << "Proc Probs " << QeNoExc << " " << QeExc << " " << ProbProjectileDiffraction
284  << " " << ProbTargetDiffraction << G4endl
285  << "ProjectileRapidity " << ProjectileRapidity << G4endl;
286 // G4int Uzhi; G4cin >> Uzhi;
287  #endif
288 
289  G4ParticleDefinition* TestParticle(0); // Uzhi Oct. 2014
290  G4double MtestPr(0.), MtestTr(0.); // Uzhi Oct. 2014
291 
292  if ( 1.0 - QeExc - QeNoExc > 0.0 ) {
293  ProbProjectileDiffraction /= ( 1.0 - QeExc - QeNoExc );
294  ProbTargetDiffraction /= ( 1.0 - QeExc - QeNoExc );
295  }
296 
297  if ( G4UniformRand() < QeExc + QeNoExc ) {
298 
299  #ifdef debugFTFexictation
300  G4cout << "Q exchange --------------------------" << G4endl;
301  #endif
302 
303  G4int NewProjCode( 0 ), NewTargCode( 0 );
304  G4int ProjQ1( 0 ), ProjQ2( 0 ), ProjQ3( 0 );
305 
306  // Projectile unpacking
307  if ( absProjectilePDGcode < 1000 ) { // projectile is meson
308  UnpackMeson( ProjectilePDGcode, ProjQ1, ProjQ2 );
309  } else { // projectile is baryon
310  UnpackBaryon( ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3 );
311  }
312 
313  // Target unpacking
314  G4int TargQ1( 0 ), TargQ2( 0 ), TargQ3( 0 );
315  UnpackBaryon( TargetPDGcode, TargQ1, TargQ2, TargQ3 );
316 
317  #ifdef debugFTFexictation
318  G4cout << "Proj Quarks " << ProjQ1 << " " << ProjQ2 << " " << ProjQ3 << G4endl
319  << "Targ Quarks " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl;
320  #endif
321 
322  // Sampling of exchanged quarks
323  G4int ProjExchangeQ( 0 );
324  G4int TargExchangeQ( 0 );
325 
326  if ( absProjectilePDGcode < 1000 )
327  { // projectile is meson
328 
329  if ( ProjQ1 > 0 )
330  { // ProjQ1 is quark
331  ProjExchangeQ = ProjQ1;
332  //------------------------------- Exchange of non-identical quarks is allowed
333  G4int NpossibleStates=3; // =====================================================
334 //
335  NpossibleStates=0; // =====================================================
336  if(ProjQ1 != TargQ1) NpossibleStates++;
337  if(ProjQ1 != TargQ2) NpossibleStates++;
338  if(ProjQ1 != TargQ3) NpossibleStates++;
339 //
340  G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
341 //G4cout<<"NpossibleStates Nsampled "<<NpossibleStates<<" "<<Nsampled<<G4endl;
342 //
343  NpossibleStates=0;
344  if(ProjQ1 != TargQ1)
345  {
346  NpossibleStates++;
347  if(NpossibleStates == Nsampled) {TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ1 = TargExchangeQ;}
348  }
349  if(ProjQ1 != TargQ2)
350  {
351  NpossibleStates++;
352  if(NpossibleStates == Nsampled) {TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ1 = TargExchangeQ;}
353  }
354  if(ProjQ1 != TargQ3)
355  {
356  NpossibleStates++;
357  if(NpossibleStates == Nsampled) {TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ1 = TargExchangeQ;}
358  }
359 //
360 
361 //if(Nsampled == 1) {TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ1 = TargExchangeQ;}
362 //else if(Nsampled == 2) {TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ1 = TargExchangeQ;}
363 //else {TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ1 = TargExchangeQ;}
364  }
365  else
366  { // ProjQ2 is quark
367  ProjExchangeQ = ProjQ2;
368  //------------------------------- Exchange of non-identical quarks is allowed
369  G4int NpossibleStates=3;
370 //
371  NpossibleStates=0;
372  if(ProjQ2 != TargQ1) NpossibleStates++;
373  if(ProjQ2 != TargQ2) NpossibleStates++;
374  if(ProjQ2 != TargQ3) NpossibleStates++;
375 //
376  G4int Nsampled = G4RandFlat::shootInt( G4long( NpossibleStates ) ) + 1;
377 //
378  NpossibleStates=0;
379  if(ProjQ2 != TargQ1)
380  {
381  NpossibleStates++;
382  if(NpossibleStates == Nsampled) {TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ2 = TargExchangeQ;}
383  }
384  if(ProjQ2 != TargQ2)
385  {
386  NpossibleStates++;
387  if(NpossibleStates == Nsampled) {TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ2 = TargExchangeQ;}
388  }
389  if(ProjQ2 != TargQ3)
390  {
391  NpossibleStates++;
392  if(NpossibleStates == Nsampled) {TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ2 = TargExchangeQ;}
393  }
394 //
395 
396 //if(Nsampled == 1) {TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjQ2 = TargExchangeQ;}
397 //else if(Nsampled == 2) {TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjQ2 = TargExchangeQ;}
398 //else {TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjQ2 = TargExchangeQ;}
399  } // End of if ( ProjQ1 > 0 )
400 
401  #ifdef debugFTFexictation
402  G4cout << "Exchanged Qs in Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
403  #endif
404 
405  G4int aProjQ1 = std::abs( ProjQ1 );
406  G4int aProjQ2 = std::abs( ProjQ2 );
407 
408  G4bool ProjExcited = false; // Uzhi Oct 2014
409 
410  G4int attempts=0; // Uzhi Oct 2014 start
411  while(attempts < 50) /* Loop checking, 10.08.2015, A.Ribon */
412  {// Determination of a new projectile ID which garanty energy-momentum conservation
413  attempts++;
414 
415  G4double Ksi = G4UniformRand();
416 
417  if ( aProjQ1 == aProjQ2 )
418  {
419  if ( aProjQ1 != 3 )
420  {
421  NewProjCode = 111; // Pi0-meson
422  if ( Ksi < 0.5 )
423  {
424  NewProjCode = 221; // Eta -meson
425  if ( Ksi < 0.25 ) {NewProjCode = 331;} // Eta'-meson
426  }
427  } else
428  {
429  NewProjCode = 221; // Eta -meson
430  if( Ksi < 0.5 ) {NewProjCode = 331;} // Eta'-meson
431  }
432  } else
433  {
434  if ( aProjQ1 > aProjQ2 )
435  {
436  NewProjCode = aProjQ1*100 + aProjQ2*10 + 1;
437  } else
438  {
439  NewProjCode = aProjQ2*100 + aProjQ1*10 + 1;
440  }
441  }
442 
443  #ifdef debugFTFexictation
444  G4cout << "NewProjCode " << NewProjCode << G4endl;
445  #endif
446 
447  ProjExcited = false;
448  if ( G4UniformRand() < 0.5 )
449  {
450  NewProjCode += 2; // Excited meson
451  ProjExcited = true;
452  }
453 // if ( aProjQ1 != aProjQ2 ) NewProjCode *= ( ProjectilePDGcode / absProjectilePDGcode ); // Uzhi 27 Nov. 2014
454 // Uzhi 27 Nov. 2014
455  G4int Qquarks=0;
456  if ( aProjQ1 == 1 ) {Qquarks -= ProjQ1;}
457  else if( aProjQ1 == 2 ) {Qquarks += ProjQ1;}
458  else {Qquarks -= ProjQ1/aProjQ1;}
459 
460  if ( aProjQ2 == 1 ) {Qquarks -= ProjQ2;}
461  else if( aProjQ2 == 2 ) {Qquarks += ProjQ2;}
462  else {Qquarks -= ProjQ2/aProjQ2;}
463 
464  if( Qquarks < 0 ) NewProjCode *=(-1);
465 // Uzhi 27 Nov. 2014
466 
467  #ifdef debugFTFexictation
468  G4cout << "NewProjCode +2 or 0 " << NewProjCode << G4endl;
469  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
470  G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquarks<<G4endl;
471  G4cout<<"+++++++++++++++++++++++++++++++++++++++"<<G4endl;
472  #endif
473 
474 // --------------------------------------------------------------------------------- Proj
475  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
476  if(!TestParticle) continue;
477 
478 //MminProjectile=TestParticle->GetPDGMass(); // ??????????????????????
479  MminProjectile=BrW.GetMinimumMass(TestParticle);
480 
481  if(SqrtS-M0target < MminProjectile) continue;
482 
483  MtestPr = BrW.SampleMass(TestParticle, TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth());
484 // G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode )->GetPDGMass(); // Uzhi 2014
485 
486  #ifdef debugFTFexictation
487  G4cout << "TestParticle Name " << NewProjCode << " " << TestParticle->GetParticleName()<< G4endl;
488  G4cout << "MtestPart MtestPart0 "<<MtestPr<<" "<<TestParticle->GetPDGMass()<<G4endl;
489  G4cout << "M0projectile projectile PDGMass " << M0projectile << " "
490  << projectile->GetDefinition()->GetPDGMass() << G4endl;
491  #endif
492 
493 // --------------------------------------------------------------------------------- Targ
494  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
495 
496  #ifdef debugFTFexictation
497  G4cout << "New TrQ " << TargQ1 << " " << TargQ2 << " " << TargQ3 << G4endl
498  << "NewTargCode " << NewTargCode << G4endl;
499  #endif
500 
501  if( TargQ1 != TargQ2 && TargQ1 != TargQ3 && TargQ2 != TargQ3 )
502  { // Lambda or Sigma0 ???
503  if ( G4UniformRand() < 0.5 ) {NewTargCode+=2;}
504  else { if ( G4UniformRand() < 0.75 ) NewTargCode=3122;}
505  }
506  else if( TargQ1 == TargQ2 && TargQ1 == TargQ3 )
507  {
508  NewTargCode += 2; ProjExcited = true; //Create Delta isobar
509  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Delta was the target
510  if ( G4UniformRand() > DeltaProbAtQuarkExchange )
511  {NewTargCode += 2; ProjExcited = true;} // Save Delta isobar
512  else {} // De-excite initial Delta isobar
513  } else if ( ! ProjExcited &&
514  G4UniformRand() < DeltaProbAtQuarkExchange && // Nucleon was the target
515  SqrtS > M0projectile + DeltaMass ) { // Create Delta isobar
516  NewTargCode +=2; // Save initial nucleon
517  } else {}
518 
519  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
520 
521  if(!TestParticle) continue;
522 
523  #ifdef debugFTFexictation
524  G4cout << "New targ " << NewTargCode << " " << TestParticle->GetParticleName() << G4endl;
525  #endif
526 
527  MminTarget=BrW.GetMinimumMass(TestParticle);
528 
529  if(SqrtS-MtestPr < MminTarget) continue;
530 
531  MtestTr = BrW.SampleMass(TestParticle,TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth());
532 
533  if(SqrtS > MtestPr+MtestTr) break;
534  } // End of while(attempts < 50)//===============================
535 
536  if(attempts >= 50) return false; // ==============================
537 /*
538  if ( MtestPr > Pprojectile.mag() ) {M0projectile = MtestPr;}
539  else
540  {
541  if ( std::abs( MtestPr - M0projectile ) //projectile->GetDefinition()->GetPDGMass() ) // Uzhi Oct. 2014
542  < 140.0*MeV ) {
543  M0projectile = MtestPr;
544  }
545  }
546 */
547  if ( MtestPr >= Pprojectile.mag() ) {M0projectile = MtestPr;} // Uzhi 18 Nov. 2014
548  else if (projectile->GetStatus() != 0 ) {M0projectile = MtestPr;} // Uzhi 18 Nov. 2014
549 
550 
551  #ifdef debugFTFexictation
552  G4cout << "M0projectile After check " << M0projectile << G4endl;
553  #endif
554 
555  M0projectile2 = M0projectile * M0projectile;
556  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
557  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
558 
559  if ( MtestTr >= Ptarget.mag() ) {M0target = MtestTr;} // Uzhi 18 Nov. 2014
560  else if (target->GetStatus() != 0 ) {M0target = MtestTr;} // Uzhi 18 Nov. 2014
561 /*
562 M0target = MtestTr; // Uzhi 18 Nov. 2014
563  if ( MtestTr > Ptarget.mag() ) {M0target = MtestTr;}
564  else
565  {
566  if ( std::abs( MtestTr - M0target ) // target->GetDefinition()->GetPDGMass() ) Uzhi Oct. 2014
567  < 140.0*MeV ) {
568  M0target = MtestTr;
569  }
570  }
571 */
572  M0target2 = M0target * M0target;
573 
574  #ifdef debugFTFexictation
575  G4cout << "New targ M0 M0^2 " << M0target << " " << M0target2 << G4endl;
576  #endif
577 
578  TargetDiffStateMinMass = M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
579  TargetNonDiffStateMinMass = M0target + 220.0*MeV; // 220 MeV=m_pi+80 MeV;
580 
581  } else { // of the if ( absProjectilePDGcode < 1000 ) ;
582  // The projectile is baryon now
583 
584  G4double Same = theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0.
585 // G4bool ProjDeltaHasCreated( false ); // Uzhi Oct. 2014
586 // G4bool TargDeltaHasCreated( false ); // Uzhi Oct. 2014
587 
588  G4double Ksi = G4UniformRand();
589  if ( G4UniformRand() < 0.5 ) { // Sampling exchange quark from proj. or targ.
590  // Sampling exchanged quark from the projectile
591 
592  if ( Ksi < 0.333333 ) {
593  ProjExchangeQ = ProjQ1;
594  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
595  ProjExchangeQ = ProjQ2;
596  } else {
597  ProjExchangeQ = ProjQ3;
598  }
599 
600  if ( ProjExchangeQ != TargQ1 || G4UniformRand() < Same ) {
601  TargExchangeQ = TargQ1; TargQ1 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
602  } else {
603  if ( ProjExchangeQ != TargQ2 || G4UniformRand() < Same ) {
604  TargExchangeQ = TargQ2; TargQ2 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
605  } else {
606  TargExchangeQ = TargQ3; TargQ3 = ProjExchangeQ; ProjExchangeQ = TargExchangeQ;
607  }
608  }
609 
610  #ifdef debugFTFexictation
611  G4cout << "Exchange Qs Pr Tr " << ProjExchangeQ << " " << TargExchangeQ << G4endl;
612  #endif
613 
614  if ( Ksi < 0.333333 ) {
615  ProjQ1 = ProjExchangeQ;
616  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
617  ProjQ2 = ProjExchangeQ;
618  } else {
619  ProjQ3 = ProjExchangeQ;
620  }
621 
622  } else { // Sampling exchanged quark from the target
623 
624  if ( Ksi < 0.333333 ) {
625  TargExchangeQ = TargQ1;
626  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
627  TargExchangeQ = TargQ2;
628  } else {
629  TargExchangeQ = TargQ3;
630  }
631  if ( TargExchangeQ != ProjQ1 || G4UniformRand() < Same ) {
632  ProjExchangeQ = ProjQ1; ProjQ1 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
633  } else {
634  if ( TargExchangeQ != ProjQ2 || G4UniformRand() < Same ) {
635  ProjExchangeQ = ProjQ2; ProjQ2 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
636  } else {
637  ProjExchangeQ = ProjQ3; ProjQ3 = TargExchangeQ; TargExchangeQ = ProjExchangeQ;
638  }
639  }
640 
641  if ( Ksi < 0.333333 ) {
642  TargQ1 = TargExchangeQ;
643  } else if ( 0.333333 <= Ksi && Ksi < 0.666667 ) {
644  TargQ2 = TargExchangeQ;
645  } else {
646  TargQ3 = TargExchangeQ;
647  }
648 
649  } // End of quark sampling for the baryons
650 
651  NewProjCode = NewNucleonId( ProjQ1, ProjQ2, ProjQ3 );
652  NewTargCode = NewNucleonId( TargQ1, TargQ2, TargQ3 );
653 
654  G4int attempts=0; // Uzhi Oct 2014 start
655  while(attempts < 50) /* Loop checking, 10.08.2015, A.Ribon */
656  {// Determination of a new projectile ID which garanty energy-momentum conservation
657  attempts++;
658 
659  if ( ProjQ1 == ProjQ2 && ProjQ1 == ProjQ3 ) {
660  NewProjCode += 2; // ProjDeltaHasCreated = true; // Uzhi Oct. 2014
661  } else if ( projectile->GetDefinition()->GetPDGiIsospin() == 3 ) { // Projectile was Delta
662  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
663  NewProjCode += 2; //ProjDeltaHasCreated = true; // Uzhi Oct. 2014
664  } else {
665  NewProjCode += 0; //ProjDeltaHasCreated = false; // Uzhi Oct. 2014
666  }
667  } else { // Projectile was Nucleon
668  if ( G4UniformRand() < DeltaProbAtQuarkExchange && SqrtS > DeltaMass + M0target ) {
669  NewProjCode += 2; //ProjDeltaHasCreated = true; // Uzhi Oct. 2014
670  } else {
671  NewProjCode += 0; //ProjDeltaHasCreated = false; // Uzhi Oct. 2014
672  }
673  }
674 
675  if ( TargQ1 == TargQ2 && TargQ1 == TargQ3 ) {
676  NewTargCode += 2; //TargDeltaHasCreated = true; // Uzhi Oct. 2014
677  } else if ( target->GetDefinition()->GetPDGiIsospin() == 3 ) { // Target was Delta
678  if ( G4UniformRand() > DeltaProbAtQuarkExchange ) {
679  NewTargCode += 2; //TargDeltaHasCreated = true; // Uzhi Oct. 2014
680  } else {
681  NewTargCode += 0; //TargDeltaHasCreated = false; // Uzhi Oct. 2014
682  }
683  } else { // Target was Nucleon
684  if ( G4UniformRand() < DeltaProbAtQuarkExchange && SqrtS > M0projectile + DeltaMass ) {
685  NewTargCode += 2; //TargDeltaHasCreated = true; // Uzhi Oct. 2014
686  } else {
687  NewTargCode += 0; //TargDeltaHasCreated = false; // Uzhi Oct. 2014
688  }
689  }
690 
691  #ifdef debugFTFexictation
692  G4cout << "NewProjCode NewTargCode " << NewProjCode << " " << NewTargCode << G4endl;
693 // G4int Uzhi; G4cin >> Uzhi;
694  #endif
695 
696  if ( absProjectilePDGcode == NewProjCode && absTargetPDGcode == NewTargCode ) {
697  } // Nothing was changed! It is not right!?
698 
699  // Forming baryons
700 /*
701  if ( G4UniformRand() > 0.5 ) { // Uzhi Oct. 2014
702  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
703  } else {
704  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
705  }
706 */
707 /*
708  if ( ProjDeltaHasCreated ) {
709  if ( G4UniformRand() > 0.5 ) {
710  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
711  } else {
712  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
713  }
714  }
715 
716  if ( TargDeltaHasCreated ) {
717  if ( G4UniformRand() > 0.5 ) {
718  ProbProjectileDiffraction = 1.0; ProbTargetDiffraction = 0.0;
719  } else {
720  ProbProjectileDiffraction = 0.0; ProbTargetDiffraction = 1.0;
721  }
722  }
723 */
724 // --------------------------------------------------------------------------------- Proj
725  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode );
726  if(!TestParticle) continue;
727 
728  MminProjectile=BrW.GetMinimumMass(TestParticle);
729 
730  if(SqrtS-M0target < MminProjectile) continue;
731 
732  MtestPr = BrW.SampleMass(TestParticle,TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth());
733 
734 // --------------------------------------------------------------------------------- Targ
735  TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode );
736  if(!TestParticle) continue;
737 
738  MminTarget=BrW.GetMinimumMass(TestParticle);
739 
740  if(SqrtS-MtestPr < MminTarget) continue;
741 
742  MtestTr = BrW.SampleMass(TestParticle,TestParticle->GetPDGMass() + 5.0*TestParticle->GetPDGWidth());
743 
744  if(SqrtS > MtestPr+MtestTr) break;
745  } // End of while(attempts < 50)//===============================
746 
747  if(attempts >= 50) return false; // ==============================
748 
749  if ( MtestPr >= Pprojectile.mag() ) {M0projectile = MtestPr;}
750  else if (projectile->GetStatus() != 0 ) {M0projectile = MtestPr;} // Uzhi 18 Nov. 2014
751  M0projectile2 = M0projectile * M0projectile;
752  ProjectileDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
753  ProjectileNonDiffStateMinMass = M0projectile + 220.0*MeV; //220 MeV=m_pi+80 MeV
754 
755  if ( MtestTr >= Ptarget.mag() ) {M0target = MtestTr;}
756  else if (target->GetStatus() != 0 ) {M0target = MtestTr;} // Uzhi 18 Nov. 2014
757  M0target2 = M0target * M0target;
758  TargetDiffStateMinMass = M0target + 220.0*MeV; //220 MeV=m_pi+80 MeV;
759  TargetNonDiffStateMinMass = M0target + 220.0*MeV; //220 MeV=m_pi+80 MeV;
760 
761  } // End of if ( absProjectilePDGcode < 1000 )
762 //--------------------------------------------------------------------------------------
763 
764  // If we assume that final state hadrons after the charge exchange will be
765  // in the ground states, we have to put
766  if ( SqrtS < M0projectile + M0target ) return false;
767 
768  PZcms2 = ( S*S + M0projectile2*M0projectile2 + M0target2*M0target2
769  - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S;
770 
771  #ifdef debugFTFexictation
772  G4cout << "At the end// NewProjCode " << NewProjCode << G4endl
773  << "At the end// NewTargCode " << NewTargCode << G4endl
774  << "M0pr M0tr SqS " << M0projectile << " " << M0target << " " << SqrtS << G4endl
775  << "M0pr2 M0tr2 SqS " << M0projectile2 << " " << M0target2 << " " << SqrtS << G4endl
776  << "PZcms2 after the change " << PZcms2 << G4endl << G4endl;
777  #endif
778 
779  if ( PZcms2 < 0 ) return false; // It can be if energy is not sufficient for Delta
780 
781  projectile->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewProjCode ) );
782  target->SetDefinition( G4ParticleTable::GetParticleTable()->FindParticle( NewTargCode ) );
783 
784  PZcms = std::sqrt( PZcms2 );
785  Pprojectile.setPz( PZcms );
786  Pprojectile.setE( std::sqrt( M0projectile2 + PZcms2 ) );
787  Ptarget.setPz( -PZcms );
788  Ptarget.setE( std::sqrt( M0target2 + PZcms2 ) );
789 
790  if(projectile->GetStatus() != 0 ) projectile->SetStatus(2); // Uzhi 18 Nov. 2014
791  if(target->GetStatus() != 0 ) target->SetStatus(2); // Uzhi 18 Nov. 2014
792 
793  #ifdef debugFTFexictation
794  G4cout << "Proj Targ and Proj+Targ in CMS" << G4endl << Pprojectile << G4endl << Ptarget
795  << G4endl << Pprojectile + Ptarget << G4endl;
796  #endif
797 
798 //--------------------- Check for possible excitation of the participants -------------------
799  if((SqrtS < M0projectile + TargetDiffStateMinMass) || // Uzhi Oct 2014
800  (SqrtS < ProjectileDiffStateMinMass + M0target) ||
801  (ProbOfDiffraction == 0.) ) ProbExc=0.;// Uzhi Oct 2014
802 
803  if ( G4UniformRand() > ProbExc ) { // Make elastic scattering
804 
805  #ifdef debugFTFexictation
806  G4cout << "Make elastic scattering of new hadrons" << G4endl;
807  #endif
808 
809  Pprojectile.transform( toLab );
810  Ptarget.transform( toLab );
811 
812  projectile->Set4Momentum( Pprojectile );
813  target->Set4Momentum( Ptarget );
814 
815  G4bool Result = theElastic->ElasticScattering( projectile, target, theParameters );
816 
817  #ifdef debugFTFexictation
818  G4cout << "Result of el. scatt " << Result << G4endl << "Proj Targ and Proj+Targ in Lab"
819  << G4endl << projectile->Get4Momentum() << G4endl << target->Get4Momentum() << G4endl
820  << projectile->Get4Momentum() + target->Get4Momentum() << " " << (projectile->Get4Momentum() + target->Get4Momentum()).mag() << G4endl;
821  #endif
822 
823 //Uzhi_QEnex++;
824  return Result;
825  }
826 //Uzhi_QEexc++;
827 
828  #ifdef debugFTFexictation
829  G4cout << "Make excitation of new hadrons" << G4endl;
830  #endif
831 
832 // Redefinition of ProbOfDiffraction because the probabilities are changed due to quark exchange
833 
834  ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction; // Uzhi Oct. 2014
835  if ( ProbOfDiffraction != 0.0 ) {
836  ProbProjectileDiffraction /= ProbOfDiffraction;
837  ProbTargetDiffraction /= ProbOfDiffraction;
838  }
839 //Uzhi_QEnex++;
840  } // End of if ( G4UniformRand() < QeExc + QeNoExc ) , i.e. of the charge exchange part
841 
842  ProbOfDiffraction = ProbProjectileDiffraction + ProbTargetDiffraction;
843 
844  #ifdef debugFTFexictation
845  G4cout << "Excitation --------------------" << G4endl
846  << "Proj M0 MdMin MndMin " << M0projectile << " " << ProjectileDiffStateMinMass << " "
847  << ProjectileNonDiffStateMinMass << G4endl
848  << "Targ M0 MdMin MndMin " << M0target << " " << TargetDiffStateMinMass << " "
849  << TargetNonDiffStateMinMass << G4endl << "SqrtS " << SqrtS << G4endl
850  << "Prob: ProjDiff TargDiff + Sum " << ProbProjectileDiffraction << " "
851  << ProbTargetDiffraction << " " << ProbOfDiffraction << G4endl;
852  #endif
853 
854  if ( ProbOfDiffraction != 0.0 ) {
855  ProbProjectileDiffraction /= ProbOfDiffraction;
856  } else {
857  ProbProjectileDiffraction = 0.0;
858  }
859 
860  #ifdef debugFTFexictation
861  G4cout << "Prob: ProjDiff TargDiff + Sum " << ProbProjectileDiffraction << " "
862  << ProbTargetDiffraction << " " << ProbOfDiffraction << G4endl;
863  #endif
864 
865  G4double ProjectileDiffStateMinMass2 = sqr( ProjectileDiffStateMinMass );
866  G4double ProjectileNonDiffStateMinMass2 = sqr( ProjectileNonDiffStateMinMass );
867  G4double TargetDiffStateMinMass2 = sqr( TargetDiffStateMinMass );
868  G4double TargetNonDiffStateMinMass2 = sqr( TargetNonDiffStateMinMass );
869 
870  G4double Pt2;
871  G4double ProjMassT2, ProjMassT;
872  G4double TargMassT2, TargMassT;
873  G4double PMinusMin, PMinusMax;
874  //G4double PPlusMin , PPlusMax;
875  G4double TPlusMin, TPlusMax;
876  G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
877  G4LorentzVector Qmomentum;
878  G4double Qminus, Qplus;
879  G4int whilecount = 0;
880 
881  // Choose a process
882  if ( G4UniformRand() < ProbOfDiffraction ) {
883 
884  if ( G4UniformRand() < ProbProjectileDiffraction ) { // projectile diffraction
885 
886  #ifdef debugFTFexictation
887  G4cout << "projectile diffraction" << G4endl;
888  #endif
889 
890  //UzhiPrD++;
891 
892  do { // while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileDiffStateMinMass2 )
893 
894  //Uzhi_projectilediffraction = 1;
895  //Uzhi_targetdiffraction = 0;
896  //Uzhi_Mx2 = 1.0;
897 
898  // Generate pt and mass of projectile
899 
900  whilecount++;
901  if ( whilecount > 1000 ) {
902  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
903  return false; // Ignore this interaction
904  };
905 
906  // Check that the interaction is possible
907  ProjMassT2 = ProjectileDiffStateMinMass2;
908  ProjMassT = ProjectileDiffStateMinMass;
909  TargMassT2 = M0target2;
910  TargMassT = M0target;
911  if ( SqrtS < ProjMassT + TargMassT ) return false;
912 
913  PZcms2 =( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
914  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
915 
916  if ( PZcms2 < 0 ) return false;
917 
918  maxPtSquare = PZcms2;
919 
920  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
921 
922  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
923  ProjMassT2 = ProjectileDiffStateMinMass2 + Pt2;
924  ProjMassT = std::sqrt( ProjMassT2 );
925  TargMassT2 = M0target2 + Pt2;
926  TargMassT = std::sqrt( TargMassT2 );
927  if ( SqrtS < ProjMassT + TargMassT ) continue;
928 
929  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
930  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
931 
932  if ( PZcms2 < 0 ) continue;
933 
934  PZcms = std::sqrt( PZcms2 );
935  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
936  PMinusMax = SqrtS - TargMassT;
937 
938 // PMinusNew = PMinusMax; //ChooseP( PMinusMin, PMinusMax ); ==== Apr. 19 ==================== Vova
939  PMinusNew = ChooseP( PMinusMin, PMinusMax );
940 
941  TMinusNew = SqrtS - PMinusNew;
942  Qminus = Ptarget.minus() - TMinusNew;
943  TPlusNew = TargMassT2 / TMinusNew;
944  Qplus = Ptarget.plus() - TPlusNew;
945  Qmomentum.setPz( (Qplus - Qminus)/2 );
946  Qmomentum.setE( (Qplus + Qminus)/2 );
947 
948  } while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
949  // Repeat the sampling because there was not any excitation
950 
951 // projectile->SetStatus( 1*projectile->GetStatus() ); // Uzhi Oct 2014
952 
953  if(projectile->GetStatus() == 2) projectile->SetStatus(1); // Uzhi Oct 2014
954  if((target->GetStatus() == 1) && (target->GetSoftCollisionCount() == 0)) // Uzhi Oct 2014
955  target->SetStatus(2); // Uzhi Oct 2014
956 
957  } else { // Target diffraction
958 
959  #ifdef debugFTFexictation
960  G4cout << "Target diffraction" << G4endl;
961  #endif
962 
963  //UzhiTrD++;
964 
965  do { // while ( ( Ptarget - Qmomentum ).mag2() < TargetDiffStateMinMass2 )
966 
967  //Uzhi_projectilediffraction = 0;
968  //Uzhi_targetdiffraction = 1;
969  //Uzhi_Mx2 = 1.0;
970 
971  // Generate pt and target mass
972 
973  whilecount++;
974  if ( whilecount > 1000 ) {
975  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
976  return false; // Ignore this interaction
977  };
978 
979  // Check that the interaction is possible
980  ProjMassT2 = M0projectile2;
981  ProjMassT = M0projectile;
982 
983  TargMassT2 = TargetDiffStateMinMass2;
984  TargMassT = TargetDiffStateMinMass;
985 
986  if ( SqrtS < ProjMassT + TargMassT ) return false;
987 
988  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
989  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
990 
991  if ( PZcms2 < 0 ) return false;
992 
993  maxPtSquare = PZcms2;
994 
995  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
996 
997  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
998  ProjMassT2 = M0projectile2 + Pt2;
999  ProjMassT = std::sqrt( ProjMassT2 );
1000  TargMassT2 = TargetDiffStateMinMass2 + Pt2;
1001  TargMassT = std::sqrt( TargMassT2 );
1002  if ( SqrtS < ProjMassT + TargMassT ) continue;
1003 
1004  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1005  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1006 
1007  if ( PZcms2 < 0 ) continue;
1008 
1009  PZcms = std::sqrt( PZcms2 );
1010  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
1011 //TPlusMax = std::sqrt( TargMassT2 + PZcms2 ) + PZcms;
1012  TPlusMax = SqrtS - ProjMassT;
1013 
1014  TPlusNew = ChooseP( TPlusMin, TPlusMax );
1015 //TPlusNew = TPlusMin;
1016 
1017  PPlusNew = SqrtS - TPlusNew;
1018  Qplus = PPlusNew - Pprojectile.plus();
1019  PMinusNew = ProjMassT2 / PPlusNew;
1020  Qminus = PMinusNew - Pprojectile.minus();
1021  Qmomentum.setPz( (Qplus - Qminus)/2 );
1022  Qmomentum.setE( (Qplus + Qminus)/2 );
1023 
1024  } while ( ( Ptarget - Qmomentum ).mag2() < TargetDiffStateMinMass2 ); /* Loop checking, 10.08.2015, A.Ribon */
1025  // Repeat the sampling because there was not any excitation
1026 
1027 // target->SetStatus( 1*target->GetStatus() ); // Uzhi Oct 2014
1028 
1029  if((projectile->GetStatus() == 1) && (projectile->GetSoftCollisionCount() == 0)) // Uzhi Oct 2014
1030  projectile->SetStatus(2); // Uzhi Oct 2014
1031  if(target->GetStatus() == 2) target->SetStatus(1); // Uzhi Oct 2014
1032 
1033  } // End of if ( G4UniformRand() < ProbProjectileDiffraction )
1034 
1035  } else { // Non-diffraction process
1036 
1037  #ifdef debugFTFexictation
1038  G4cout << "Non-diffraction process" << G4endl;
1039  #endif
1040 
1041 //UzhiND++;
1042 //Uzhi_QEnex++;
1043  do { // while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileNonDiffStateMinMass2 || ...
1044 
1045  //Uzhi_projectilediffraction = 0;
1046  //Uzhi_targetdiffraction = 0;
1047  //Uzhi_Mx2 = 1.0;
1048 
1049  // Generate pt and masses
1050 
1051  whilecount++;
1052  if ( whilecount > 1000 ) {
1053  Qmomentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1054  return false; // Ignore this interaction
1055  };
1056 
1057  // Check that the interaction is possible
1058  ProjMassT2 = ProjectileNonDiffStateMinMass2;
1059  ProjMassT = ProjectileNonDiffStateMinMass;
1060  TargMassT2 = TargetNonDiffStateMinMass2;
1061  TargMassT = TargetNonDiffStateMinMass;
1062  if ( SqrtS < ProjMassT + TargMassT ) return false;
1063 
1064  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1065  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1066 
1067  if ( PZcms2 < 0 ) return false;
1068 
1069  maxPtSquare = PZcms2;
1070 
1071  Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0 );
1072 
1073  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
1074  ProjMassT2 = ProjectileNonDiffStateMinMass2 + Pt2;
1075  ProjMassT = std::sqrt( ProjMassT2 );
1076  TargMassT2 = TargetNonDiffStateMinMass2 + Pt2;
1077  TargMassT = std::sqrt( TargMassT2 );
1078  if ( SqrtS < ProjMassT + TargMassT ) continue;
1079 
1080  PZcms2 =( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
1081  -2.0*S*ProjMassT2 - 2.0*S*TargMassT2 -2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
1082 
1083  if ( PZcms2 < 0 ) continue;
1084 
1085  PZcms = std::sqrt( PZcms2 );
1086  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
1087 //PMinusMax = std::sqrt( ProjMassT2 + PZcms2 ) + PZcms;
1088  PMinusMax = SqrtS - TargMassT;
1089 
1090  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
1091 //TPlusMax = std::sqrt( TargMassT2 + PZcms2 ) + PZcms;
1092  TPlusMax = SqrtS - ProjMassT; // Uzhi 18 Sept. 2014
1093 
1094 /*
1095  if ( G4UniformRand() < ProbLogDistrPrD ) { // Uzhi Oct 2014
1096  PMinusNew = ChooseP( PMinusMin, PMinusMax );
1097  } else {
1098  PMinusNew = ( PMinusMax - PMinusMin )*G4UniformRand() + PMinusMin;
1099  }
1100 // Qminus = PMinusNew - Pprojectile.minus();
1101 */
1102 
1103 // TPlusMax = SqrtS - PMinusNew;
1104 
1105  if ( G4UniformRand() < ProbLogDistr ) {
1106  PMinusNew = ChooseP( PMinusMin, PMinusMax );
1107  TPlusNew = ChooseP( TPlusMin, TPlusMax );
1108  } else {
1109  PMinusNew = ( PMinusMax - PMinusMin )*G4UniformRand() + PMinusMin;
1110  TPlusNew = ( TPlusMax - TPlusMin )*G4UniformRand() + TPlusMin;
1111  }
1112 
1113  Qminus = PMinusNew - Pprojectile.minus();
1114 
1115  Qplus = -( TPlusNew - Ptarget.plus() );
1116  Qmomentum.setPz( (Qplus - Qminus)/2 );
1117  Qmomentum.setE( (Qplus + Qminus)/2 );
1118 
1119  #ifdef debugFTFexictation
1120  G4cout << ( Pprojectile + Qmomentum ).mag2() << " " << ProjectileNonDiffStateMinMass2
1121  << G4endl << ( Ptarget - Qmomentum ).mag2() << " "
1122  << TargetNonDiffStateMinMass2 << G4endl;
1123  G4cout<<"To continue - enter any integer"<<G4endl;
1124 // G4int Uzhi; G4cin >> Uzhi;
1125  #endif
1126 
1127  } while ( ( Pprojectile + Qmomentum ).mag2() < ProjectileNonDiffStateMinMass2 || //No double Diffraction
1128  ( Ptarget - Qmomentum ).mag2() < TargetNonDiffStateMinMass2 || // ); //
1129  ( Pprojectile + Qmomentum ).pz() < 0.); /* Loop checking, 10.08.2015, A.Ribon */
1130 
1131  projectile->SetStatus( 0*projectile->GetStatus() );
1132  target->SetStatus( 0*target->GetStatus() );
1133 
1134  } // End of if ( G4UniformRand() < ProbOfDiffraction )
1135 
1136  Pprojectile += Qmomentum;
1137  Ptarget -= Qmomentum;
1138 
1139  // Transform back and update SplitableHadron Participant.
1140  Pprojectile.transform( toLab );
1141  Ptarget.transform( toLab );
1142 
1143  // Calculation of the creation time
1144  //Uzhi 9.11 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
1145  //Uzhi 9.11 projectile->SetPosition( target->GetPosition() );
1146  // Creation time and position of target nucleon were determined in
1147  // ReggeonCascade() of G4FTFModel
1148  //
1149  //if ( Uzhi_projectilediffraction != 0 ) {
1150  // Uzhi_Mx2 = Pprojectile.mag2(); Uzhi_modT = ( target->Get4Momentum() - Ptarget ).mag2();
1151  //}
1152  //if ( Uzhi_targetdiffraction != 0 ) {
1153  // Uzhi_Mx2 = Ptarget.mag2(); Uzhi_modT = ( projectile->Get4Momentum() - Pprojectile ).mag2();
1154  //}
1155  //if ( Uzhi_QE != 0 ) {
1156  // Uzhi_projectilediffraction = 0;
1157  // Uzhi_targetdiffraction = 0;
1158  // Uzhi_Mx2 = 1.0;
1159  //}
1160 
1161  #ifdef debugFTFexictation
1162  G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl;
1163  #endif
1164 
1165  projectile->Set4Momentum( Pprojectile );
1166  target->Set4Momentum( Ptarget );
1167  projectile->IncrementCollisionCount( 1 );
1168  target->IncrementCollisionCount( 1 );
1169 
1170  //Uzhi_projectilediffraction = UzhiPrD;
1171  //Uzhi_targetdiffraction = UzhiTrD;
1172  //Uzhi_nondiffraction = UzhiND;
1173  //G4cout << Uzhi_projectilediffraction << " " << Uzhi_targetdiffraction << " "
1174  // << Uzhi_nondiffraction << G4endl;
1175 
1176  return true;
1177 }
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:566
G4DiffractiveExcitation::ChooseP
G4double ChooseP(G4double Pmin, G4double Pmax) const
Definition: G4DiffractiveExcitation.cc:1532
G4FTFParameters::GetProjMinDiffMass
G4double GetProjMinDiffMass()
Definition: G4FTFParameters.hh:431
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:76
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4DiffractiveExcitation::UnpackBaryon
void UnpackBaryon(G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const
Definition: G4DiffractiveExcitation.cc:1606
G4FTFParameters::GetProbLogDistr
G4double GetProbLogDistr()
Definition: G4FTFParameters.hh:455
G4FTFParameters::GetTarMinDiffMass
G4double GetTarMinDiffMass()
Definition: G4FTFParameters.hh:439
G4long
long G4long
Definition: G4Types.hh:80
CLHEP::HepLorentzVector::rotateZ
HepLorentzVector & rotateZ(double)
Definition: LorentzVector.cc:80
G4int
int G4int
Definition: G4Types.hh:78
G4SampleResonance::SampleMass
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
Definition: G4SampleResonance.cc:110
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4DiffractiveExcitation::GaussianPt
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
Definition: G4DiffractiveExcitation.cc:1549
G4bool
bool G4bool
Definition: G4Types.hh:79
G4SampleResonance::GetMinimumMass
G4double GetMinimumMass(const G4ParticleDefinition *p) const
Definition: G4SampleResonance.cc:46
CLHEP::HepLorentzRotation::transform
HepLorentzRotation & transform(const HepBoost &b)
G4FTFParameters::GetProcProb
G4double GetProcProb(const G4int ProcN, const G4double y)
Definition: G4FTFParameters.cc:754
G4FTFParameters::GetAveragePt2
G4double GetAveragePt2()
Definition: G4FTFParameters.hh:447
G4FTFParameters::GetTarMinNonDiffMass
G4double GetTarMinNonDiffMass()
Definition: G4FTFParameters.hh:443
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
target
cout<< "-> Edep in the target
Definition: analysis.C:54
CLHEP::HepLorentzVector::rotateY
HepLorentzVector & rotateY(double)
Definition: LorentzVector.cc:76
G4FTFParameters::GetDeltaProbAtQuarkExchange
G4double GetDeltaProbAtQuarkExchange()
Definition: G4FTFParameters.hh:423
G4FTFParameters::GetProbOfSameQuarkExchange
G4double GetProbOfSameQuarkExchange()
Definition: G4FTFParameters.hh:427
G4endl
#define G4endl
Definition: G4ios.hh:61
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4DiffractiveExcitation::NewNucleonId
G4int NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const
Definition: G4DiffractiveExcitation.cc:1617
G4FTFParameters::GetProjMinNonDiffMass
G4double GetProjMinNonDiffMass()
Definition: G4FTFParameters.hh:435
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4ElasticHNScattering::ElasticScattering
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
Definition: G4ElasticHNScattering.cc:64
G4DiffractiveExcitation::UnpackMeson
void UnpackMeson(G4int IdPDG, G4int &Q1, G4int &Q2) const
Definition: G4DiffractiveExcitation.cc:1584
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GaussianPt()

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

Definition at line 1549 of file G4DiffractiveExcitation.cc.

1549  {
1550  // @@ this method is used in FTFModel as well. Should go somewhere common!
1551  G4double Pt2( 0.0 );
1552  if ( AveragePt2 <= 0.0 ) {
1553  Pt2 = 0.0;
1554  } else {
1555  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
1556  ( G4Exp( -maxPtSquare/AveragePt2 ) - 1.0 ) );
1557  }
1558  G4double Pt = std::sqrt( Pt2 );
1559  G4double phi = G4UniformRand() * twopi;
1560  return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 );
1561 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
twopi
static const double twopi
Definition: G4SIunits.hh:75
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetQuarkFractionOfKink()

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

Definition at line 1566 of file G4DiffractiveExcitation.cc.

1566  {
1567  G4double z, yf;
1568  const G4int maxNumberOfLoops = 10000;
1569  G4int loopCounter = 0;
1570  do {
1571  z = zmin + G4UniformRand() * (zmax - zmin);
1572  yf = z*z + sqr(1.0 - z);
1573  } while ( ( G4UniformRand() > yf ) &&
1574  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1575  if ( loopCounter >= maxNumberOfLoops ) {
1576  z = 0.5*(zmin + zmax); // Just something acceptable, without any physics consideration.
1577  }
1578  return z;
1579 }
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

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

1617  {
1618  G4int TmpQ( 0 );
1619  if ( Q3 > Q2 ) {
1620  TmpQ = Q2;
1621  Q2 = Q3;
1622  Q3 = TmpQ;
1623  } else if ( Q3 > Q1 ) {
1624  TmpQ = Q1;
1625  Q1 = Q3;
1626  Q3 = TmpQ;
1627  }
1628  if ( Q2 > Q1 ) {
1629  TmpQ = Q1;
1630  Q1 = Q2;
1631  Q2 = TmpQ;
1632  }
1633  G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + 2;
1634  return NewCode;
1635 }
G4int
int G4int
Definition: G4Types.hh:78
Here is the caller graph for this function:

◆ operator!=()

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

Definition at line 1665 of file G4DiffractiveExcitation.cc.

1665  {
1666  throw G4HadronicException( __FILE__, __LINE__,
1667  "G4DiffractiveExcitation != operator not meant to be called" );
1668 }

◆ operator=()

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

Definition at line 1648 of file G4DiffractiveExcitation.cc.

1648  {
1649  throw G4HadronicException( __FILE__, __LINE__,
1650  "G4DiffractiveExcitation = operator not meant to be called" );
1651  return *this;
1652 }

◆ operator==()

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

Definition at line 1657 of file G4DiffractiveExcitation.cc.

1657  {
1658  throw G4HadronicException( __FILE__, __LINE__,
1659  "G4DiffractiveExcitation == operator not meant to be called" );
1660 }

◆ UnpackBaryon()

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

Definition at line 1606 of file G4DiffractiveExcitation.cc.

1607  {
1608  Q1 = IdPDG / 1000;
1609  Q2 = (IdPDG % 1000) / 100;
1610  Q3 = (IdPDG % 100) / 10;
1611  return;
1612 }
Here is the caller graph for this function:

◆ UnpackMeson()

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

Definition at line 1584 of file G4DiffractiveExcitation.cc.

1584  {
1585  G4int absIdPDG = std::abs( IdPDG );
1586  if(!((absIdPDG == 111)||(absIdPDG == 221)||(absIdPDG == 331))) // Uzhi Oct. 2014
1587  { // Ordinary mesons =======================
1588  Q1 = absIdPDG / 100;
1589  Q2 = (absIdPDG % 100) / 10;
1590  G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) % 2 );
1591  if ( IdPDG < 0 ) anti *= -1;
1592  Q1 *= anti;
1593  Q2 *= -1 * anti;
1594  }
1595  else
1596  { // Pi0, Eta, Eta' =======================
1597  if( G4UniformRand() < 0.5 ) {Q1 = 1; Q2 = -1;}
1598  else {Q1 = 2; Q2 = -2;}
1599  }
1600  return;
1601 }
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
Here is the caller graph for this function:
The documentation for this class was generated from the following files: