- p1[index][0] += px;
- p1[index][1] += py;
- p1[index][2] += pz;
- p1[index][3] += e;
- continue;
- }
- } // massive particles
-
- //
- // Change light-cone kinematics rel. to initial parton
- //
- Double_t eppzOld = e + pl;
- Double_t empzOld = e - pl;
-
- Double_t eppzNew = (1. - z) * eppzOld;
- Double_t empzNew = empzOld - mt2 * z / eppzOld;
- Double_t eNew0 = 0.5 * (eppzNew + empzNew);
- Double_t pzNew0 = 0.5 * (eppzNew - empzNew);
-
- Double_t ptNew;
- //
- // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
- Double_t mt2New = eppzNew * empzNew;
- if (mt2New < 1.e-8) mt2New = 0.;
-
- if (m * m > mt2New) {
+ //
+ //
+ // Transform into frame in which initial parton is along z-axis
+ //
+ TVector3 v(px, py, pz);
+ v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
+ Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
+
+ Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
+ Double_t mt2 = jt * jt + m * m;
+ Double_t zmax = 1.;
+ //
+ // Kinematic limit on z
+ //
+ if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
+ //
+ // Change light-cone kinematics rel. to initial parton
+ //
+ Double_t eppzOld = e + pl;
+ Double_t empzOld = e - pl;
+
+ Double_t eppzNew = (1. - z) * eppzOld;
+ Double_t empzNew = empzOld - mt2 * z / eppzOld;
+ Double_t eNew = 0.5 * (eppzNew + empzNew);
+ Double_t plNew = 0.5 * (eppzNew - empzNew);
+
+ Double_t jtNew;
+ //
+ // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
+ Double_t mt2New = eppzNew * empzNew;
+ if (mt2New < 1.e-8) mt2New = 0.;
+ if (z < zmax) {
+ if (m * m > mt2New) {
+ //
+ // This should not happen
+ //
+ Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
+ jtNew = 0;
+ } else {
+ jtNew = TMath::Sqrt(mt2New - m * m);
+ }
+ } else {
+ // If pT is to small (probably a leading massive particle) we scale only the energy
+ // This can cause negative masses of the radiated gluon
+ // Let's hope for the best ...
+ jtNew = jt;
+ eNew = TMath::Sqrt(plNew * plNew + mt2);
+
+ }
+ //
+ // Calculate new px, py
+ //
+ Double_t pxNew = 0;
+ Double_t pyNew = 0;
+
+ if (jt>0) {
+ pxNew = jtNew / jt * pxs;
+ pyNew = jtNew / jt * pys;
+ }
+// Double_t dpx = pxs - pxNew;
+// Double_t dpy = pys - pyNew;
+// Double_t dpz = pl - plNew;
+// Double_t de = e - eNew;
+// Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
+// printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
+// printf("New mass (2) %e %e \n", pxNew, pyNew);
+ //
+ // Rotate back
+ //
+ TVector3 w(pxNew, pyNew, plNew);
+ w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
+ pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
+
+ p1[index][0] += pxNew;
+ p1[index][1] += pyNew;
+ p1[index][2] += plNew;
+ p1[index][3] += eNew;
+ //
+ // Updated 4-momentum vectors
+ //
+ pNew[icount][0] = pxNew;
+ pNew[icount][1] = pyNew;
+ pNew[icount][2] = plNew;
+ pNew[icount][3] = eNew;
+ kNew[icount] = i;
+ icount++;
+ } // parton loop