36 {
37#ifndef QUIRKS_STANDALONE
40#endif
41
43 const G4StepPoint*
ps =
step->GetPostStepPoint();
44 const G4ParticleDefinition*
particle =
track->GetParticleDefinition();
45 const Quirk* quirk =
dynamic_cast<const Quirk*
>(
particle);
46 if (quirk == 0) return;
47
48 G4int
i = (
particle->GetPDGEncoding() > 0) ? 0 : 1;
50 m_x[
i] = G4LorentzVector(
ps->GetPosition(), CLHEP::c_light *
ps->GetGlobalTime());
51 m_p[
i] =
track->GetDynamicParticle()->Get4Momentum();
52
53 G4bool doPrint = (m_numSteps != 0 &&
track->GetCurrentStepNumber() % m_numSteps == 0);
55 G4double progress =
m_x[
i].rho();
56 if (progress >=
m_iStep[i] * m_step) {
58 doPrint = true;
59 }
60
61 if (doPrint) {
62 const InfracolorForce* string[2];
64 string[1-
i] =
string[
i]->GetReactionForce();
65
66 G4cout.precision(14);
67 G4cout <<
"quirk " <<
i <<
" step " <<
track->GetCurrentStepNumber() <<
" status " <<
track->GetTrackStatus() << G4endl;
68 G4cout <<
"x0: " <<
m_x[0] << G4endl;
69 G4cout <<
"x1: " <<
m_x[1] << G4endl;
70
71
72 G4cout <<
"dx: " <<
m_x[1] -
m_x[0] << G4endl;
73
74 G4double
force =
string[0]->GetStringForce();
75 G4LorentzVector
ss[2];
76 for(G4int i = 0;
i < 2;
i++)
ss[i] =
string[i]->GetSumStrings();
77 G4LorentzVector
dx =
ss[0] -
ss[1];
79 G4LorentzVector ptot =
ps +
m_p[0] +
m_p[1];
80 G4cout <<
"dx: " <<
dx << G4endl;
81 G4cout <<
"s0: " <<
ss[0] << G4endl;
82 G4cout <<
"s1: " <<
ss[1] << G4endl;
83 G4cout <<
"p0: " <<
m_p[0] << G4endl;
84 G4cout <<
"p1: " <<
m_p[1] << G4endl;
85 G4cout <<
"ps: " <<
ps << G4endl;
86 G4cout << "p: " << ptot << G4endl;
87
89 G4ThreeVector
L =
dx.vect().cross(p1s);
90 G4ThreeVector Excm = p1s.t() *
dx.vect() - p1s.vect() *
dx.t();
91 for(G4int i = 0;
i < 2;
i++) {
92 L +=
string[
i]->GetAngMomentum();
93 Excm +=
string[
i]->GetMomentOfE();
94 }
95 L -= Excm.cross(ptot.vect() / ptot.t());
96 Excm += ptot.t() *
m_x[0].vect() - ptot.vect() *
m_x[0].t();
97 G4cout <<
"L: " <<
L << G4endl;
98 G4cout << "cm: " << Excm / ptot.t() << G4endl;
99
100 G4LorentzVector
u = ptot / ptot.m();
101 G4LorentzVector L4(L *
u.t(), L *
u.vect());
102 G4double
p0 = sqrt(
sqr(ptot.m()/2) -
sqr(quirk->GetPDGMass()));
103 G4cout <<
"p0: " <<
p0 << G4endl;
104 G4cout <<
"min: " << - L4.m() /
p0 << G4endl;
105 G4cout <<
"max: " << (ptot.m() - 2*quirk->GetPDGMass()) /
force << G4endl;
106
107 G4cout << "n: " << string[0]->GetNStrings() << "\t" << string[1]->GetNStrings() << "\t" << string[0]->GetNStrings() + string[1]->GetNStrings() << G4endl;
108 G4cout << G4endl;
109 }
110
111 if (false) {
114 G4double t2A =
m_x[
i].t();
115 G4double t2B =
m_x[1-
i].t();
116 G4double dtA = t2A - t1A;
117 G4double dtB = t2B - t1B;
118 if (t2B > t1A && t2A > t1B && dtA > 0 && dtB > 0) {
119 G4ThreeVector vA = (
m_x[
i].vect() -
m_xprev[i].vect()) / dtA;
120 G4ThreeVector vB = (
m_x[1-
i].vect() -
m_xprev[1-i].vect()) / dtB;
121 G4ThreeVector xA =
m_xprev[
i].vect();
122 G4ThreeVector xB =
m_xprev[1-
i].vect() + vB * (t1A - t1B);
123 G4double dvsq = (vB - vA).
mag2();
124 if (dvsq > 0) {
125 G4double
dt = - (vB - vA) * (xB - xA) / dvsq;
126 G4double
t = t1A +
dt;
127 if (t1A <= t && t <= t2A && t1B <= t && t <= t2B) {
128 G4double dist = ((xB - xA) + (vB - vA) *
dt).
mag();
129 G4cout <<
"mindist: " << dist <<
" " <<
t << G4endl;
130 G4cout << G4endl;
131 }
132 }
133 }
134 }
135 }
Scalar mag() const
mag method
Scalar mag2() const
mag2 method - forward to squaredNorm()
G4LorentzVector m_xprev[2]
const InfracolorForce & GetStringIn() const
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
@ u
Enums for curvilinear frames.