printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f\n",
j, itype, eq[j], phi, l, eloss, wjtKick[j]);
}
-
quenched[j] = (zInitial[j] > 0.01);
-
}
//
wjtKick[1] = wjtKick[1] / TMath::Sqrt(Double_t(kGluons));
// this->Pylist(1);
+// Arrays to store particle 4-momenta to be changed
+//
+/*
+ Double_t** pNew = new Double_t* [numpart];
+ for (Int_t i = 0; i < numpart; i++) pNew[i] = new Double_t [4];
+ Int_t* kNew = new Int_t [numpart];
+*/
+
+ Double_t pNew[1000][4];
+ Int_t kNew[1000];
+ Int_t icount = 0;
+//
+// Radiation Loop
for (Int_t iglu = 0; iglu < kGluons; iglu++) {
- for (Int_t k = 0; k < 4; k++)
- {
- p0[0][k] = 0.; p0[1][k] = 0.;
- p1[0][k] = 0.; p1[1][k] = 0.;
- p2[0][k] = 0.; p2[1][k] = 0.;
- }
-
- Int_t nq[2] = {0, 0};
-
- for (Int_t i = 0; i < numpart; i++)
- {
- imo = fPyjets->K[2][i];
- kst = fPyjets->K[0][i];
- pdg = fPyjets->K[1][i];
+ for (Int_t isys = 0; isys < 2; isys++) {
+// Skip to next system if not quenched.
+ Double_t zHeavy = zInitial[isys];
+ if (!quenched[isys]) continue;
+//
+ while (1) {
+ icount = 0;
+ for (Int_t k = 0; k < 4; k++)
+ {
+ p0[isys][k] = 0.;
+ p1[isys][k] = 0.;
+ p2[isys][k] = 0.;
+ }
+// Loop over partons
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ imo = fPyjets->K[2][i];
+ kst = fPyjets->K[0][i];
+ pdg = fPyjets->K[1][i];
+
+
+
// Quarks and gluons only
- if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
+ if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
// Particles from hard scattering only
- if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
- if (imo != 7 && imo != 8 && imo != 1007 && imo != 1008) continue;
-
+ if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
+ if (imo != isys + 7 && imo != 1000 + isys + 7) continue;
+
// Skip comment lines
- if (kst != 1 && kst != 2) continue;
+ if (kst != 1 && kst != 2) continue;
//
// Parton kinematic
- px = fPyjets->P[0][i];
- py = fPyjets->P[1][i];
- pz = fPyjets->P[2][i];
- e = fPyjets->P[3][i];
- m = fPyjets->P[4][i];
- pt = TMath::Sqrt(px * px + py * py);
- p = TMath::Sqrt(px * px + py * py + pz * pz);
- phi = TMath::Pi() + TMath::ATan2(-py, -px);
- theta = TMath::ATan2(pt, pz);
-
+ px = fPyjets->P[0][i];
+ py = fPyjets->P[1][i];
+ pz = fPyjets->P[2][i];
+ e = fPyjets->P[3][i];
+ m = fPyjets->P[4][i];
+ pt = TMath::Sqrt(px * px + py * py);
+ p = TMath::Sqrt(px * px + py * py + pz * pz);
+ phi = TMath::Pi() + TMath::ATan2(-py, -px);
+ theta = TMath::ATan2(pt, pz);
+
//
// Save 4-momentum sum for balancing
- Int_t index = imo - 7;
- if (index >= 1000) index -= 1000;
-
- p0[index][0] += px;
- p0[index][1] += py;
- p0[index][2] += pz;
- p0[index][3] += e;
-
+ Int_t index = imo - 7;
+ if (index >= 1000) index -= 1000;
+
+ p0[index][0] += px;
+ p0[index][1] += py;
+ p0[index][2] += pz;
+ p0[index][3] += e;
+
// Don't quench radiated gluons
//
- if (imo == 1007 || imo == 1008) {
- p1[index][0] += px;
- p1[index][1] += py;
- p1[index][2] += pz;
- p1[index][3] += e;
- continue;
- }
-
+ if (imo == 1000 + isys + 7) {
+ p1[index][0] += px;
+ p1[index][1] += py;
+ p1[index][2] += pz;
+ p1[index][3] += e;
+ continue;
+ }
+
//
-
- klast[index] = i;
+
+ klast[index] = i;
+
//
// Fractional energy loss
- Double_t z = zInitial[index];
- if (!quenched[index]) continue;
- //
- //
- // 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 mt = TMath::Sqrt(mt2);
- Double_t zmin = 0.;
- Double_t zmax = 1.;
- //
- // z**2 mt**2 - pt**2 + pt'**2 > 0
- // z**2 mt**2 + mt'**2 + m**2 - pt**2
- // z**2 mt**2 + (1-z)**2 mt**2 + m**2 - pt**2
- // mt**2(z**2 + 1 + z**2 - 2z) + m**2 - pt**2
- // mt**2(2z**2 + 1 - 2z) + m**2 - pt**2 > 0
- // mt**2(2z**2 + 1 - 2z) + 2 m**2 - mt**2 > 0
- // mt**2(2z**2 - 2z) + 2 m**2 > 0
- // z mt**2 (1 - z) - m**2 < 0
- // z**2 - z + 1/4 > 1/4 - m**2/mt**2
- // (z-1/2)**2 > 1/4 - m**2/mt**2
- // |z-1/2| > sqrt(1/4 - m**2/mt**2)
- //
- // m/mt < 1/2
- // mt > 2m
- //
- if (mt < 2. * m) {
- printf("No phase space for quenching !: mt (%e) < 2 m (%e) \n", mt, m);
- p1[index][0] += px;
- p1[index][1] += py;
- p1[index][2] += pz;
- p1[index][3] += e;
- continue;
- } else {
- zmin = 0.5 - TMath::Sqrt(0.25 - m * m / mt2);
- if (z < zmin) {
- printf("No phase space for quenching ??: z (%e) < zmin (%e) \n", z, zmin);
-// z = zmin * 1.01;
-
- p1[index][0] += px;
- p1[index][1] += py;
- p1[index][2] += pz;
- p1[index][3] += e;
- continue;
-
- }
- }
- //
- // Kinematic limit on z
- //
-
- if (m > 0.) {
- zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
- if (z > zmax) {
- printf("We have to put z to the kinematic limit %e %e \n", z, zmax);
- z = 0.9999 * zmax;
- } // z > zmax
- if (z < 0.01) {
+ Double_t z = zInitial[index];
+ if (m > 0.) z = zHeavy;
+
+ //
+ //
+ // 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);
+ if (z > zmax) {
+ printf("We have to put z to the kinematic limit %e %e \n", z, zmax);
+ z = 0.9999 * zmax;
+ } // z > zmax
+
+ if (z < 0.01) {
//
// If z is too small, there is no phase space for quenching
//
- printf("No phase space for quenching ! %e \n", z);
+ printf("No phase space for quenching ! %e \n", z);
+ p1[index][0] += px;
+ p1[index][1] += py;
+ p1[index][2] += pz;
+ p1[index][3] += e;
+ continue;
+ }
+ } // massive particles
- 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.;
+ //
+ // 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 (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);
+ }
+
- if (m * m > mt2New) {
+ //
+ // Calculate new px, py
+ //
+ Double_t pxNew = jtNew / jt * pxs;
+ Double_t 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
//
- // This should not happen
+ // Check if there was phase-space for quenching
//
- Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
- ptNew = 0;
- } else {
- ptNew = TMath::Sqrt(mt2New - m * m);
- }
+ if (icount == 0) quenched[isys] = kFALSE;
+ if (!quenched[isys]) break;
+
+ for (Int_t j = 0; j < 4; j++)
+ {
+ p2[isys][j] = p0[isys][j] - p1[isys][j];
+ }
+ p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
+
+ if (p2[isys][4] > 0.) {
+ p2[isys][4] = TMath::Sqrt(p2[isys][4]);
+ break;
+ } else {
+ printf("Warning negative mass squared in system %d %f ! \n", isys, zInitial[isys]);
+
+ printf("Kinematics %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
+ if (p2[isys][4] < -0.01) {
+ printf("Negative mass squared ! Let's try to fix this by decreasing z\n");
+// this->Pylist(1);
+
+ } else {
+ p2[isys][4] = 0.;
+ break;
+ }
+ }
+ //
+ // jt-kick
+ //
+ /*
+ TVector3 v(p2[k][0], p2[k][1], p2[k][2]);
+ v.RotateZ(-phiq[k]);
+ v.RotateY(-thetaq[k]);
+ Double_t px = v.X(); Double_t py = v.Y(); Double_t pz = v.Z();
+ Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
+ Double_t jtKick = wjtKick[k] * TMath::Sqrt(-TMath::Log(r));
+ Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
+ px += jtKick * TMath::Cos(phiKick);
+ py += jtKick * TMath::Sin(phiKick);
+ TVector3 w(px, py, pz);
+ w.RotateY(thetaq[k]);
+ w.RotateZ(phiq[k]);
+ p2[k][0] = w.X(); p2[k][1] = w.Y(); p2[k][2] = w.Z();
+ p2[k][3] = TMath::Sqrt(p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2] + p2[k][4] * p2[k][4]);
+ */
+ zHeavy *= 0.98;
+ printf("zHeavy lowered to %f\n", zHeavy);
+ if (zHeavy < 0.01) {
+ printf("No success ! \n");
+ icount = 0;
+ quenched[isys] = kFALSE;
+ break;
+ }
+ } // iteration on z
- //
- // Calculate new px, py
- //
- Double_t pxNew0 = ptNew / jt * pxs;
- Double_t pyNew0 = ptNew / jt * pys;
-/*
- Double_t dpx = pxs - pxNew0;
- Double_t dpy = pys - pyNew0;
- Double_t dpz = pl - pzNew0;
- Double_t de = e - eNew0;
- Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
-*/
- //
- // Rotate back
- //
- TVector3 w(pxNew0, pyNew0, pzNew0);
- w.RotateY(thetaq[index]);
- w.RotateZ(phiq[index]);
- pxNew0 = w.X(); pyNew0 = w.Y(); pzNew0 = w.Z();
-
-
- p1[index][0] += pxNew0;
- p1[index][1] += pyNew0;
- p1[index][2] += pzNew0;
- p1[index][3] += eNew0;
- //
- // Update event record
- //
- fPyjets->P[0][i] = pxNew0;
- fPyjets->P[1][i] = pyNew0;
- fPyjets->P[2][i] = pzNew0;
- fPyjets->P[3][i] = eNew0;
- nq[index]++;
-
- }
-
- //
- // Gluons
- //
-
- for (Int_t k = 0; k < 2; k++)
- {
- //
- // Check if there was phase-space for quenching
- //
- if (nq[k] == 0) quenched[k] = kFALSE;
-
- if (!quenched[k]) continue;
-
- for (Int_t j = 0; j < 4; j++)
- {
- p2[k][j] = p0[k][j] - p1[k][j];
+// Update event record
+ for (Int_t k = 0; k < icount; k++) {
+// printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
+ fPyjets->P[0][kNew[k]] = pNew[k][0];
+ fPyjets->P[1][kNew[k]] = pNew[k][1];
+ fPyjets->P[2][kNew[k]] = pNew[k][2];
+ fPyjets->P[3][kNew[k]] = pNew[k][3];
}
- p2[k][4] = p2[k][3] * p2[k][3] - p2[k][0] * p2[k][0] - p2[k][1] * p2[k][1] - p2[k][2] * p2[k][2];
-
- if (p2[k][4] > 0.) {
- p2[k][4] = TMath::Sqrt(p2[k][4]);
- } else {
- printf("Warning negative mass squared in system %d %f ! \n", k, zInitial[k]);
- printf("Kinematics %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[k][0], p2[k][1], p2[k][2], p2[k][3], p2[k][4]);
- if (p2[k][4] < -0.1) Fatal("Boost", "Negative mass squared !");
- p2[k][4] = 0.;
- }
- //
- // jt-kick
- //
- /*
- TVector3 v(p2[k][0], p2[k][1], p2[k][2]);
- v.RotateZ(-phiq[k]);
- v.RotateY(-thetaq[k]);
- Double_t px = v.X(); Double_t py = v.Y(); Double_t pz = v.Z();
- Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
- Double_t jtKick = wjtKick[k] * TMath::Sqrt(-TMath::Log(r));
- Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
- px += jtKick * TMath::Cos(phiKick);
- py += jtKick * TMath::Sin(phiKick);
- TVector3 w(px, py, pz);
- w.RotateY(thetaq[k]);
- w.RotateZ(phiq[k]);
- p2[k][0] = w.X(); p2[k][1] = w.Y(); p2[k][2] = w.Z();
- p2[k][3] = TMath::Sqrt(p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2] + p2[k][4] * p2[k][4]);
- */
- }
-
+ } // System
//
// Add the gluons
//
-
Int_t ish = 0;
for (Int_t i = 0; i < 2; i++) {
Int_t jmin, jmax, iGlu, iNew;
numpart += ish;
(fPyjets->N) += ish;
-
+
if (ish == 1) {
fPyjets->P[0][iGlu] = p2[i][0];
fPyjets->P[1][iGlu] = p2[i][1];
//
Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
}
- } // end adding gluons
- //
+ }
+
// Check energy conservation
Double_t pxs = 0.;
Double_t pys = 0.;
TMath::Abs(pzs) > 1.e-1) {
printf("%e %e %e %e\n", pxs, pys, pzs, es);
this->Pylist(1);
- Fatal("Quench()", "4-Momentum non-conservation");
+ Fatal("Quench()", "4-Momentum non-conservation");
}
-
- } // end quenchin loop
- // Clean-up
+// this->Pylist(1);
+
+ } // end quenching loop
+// Clean-up
for (Int_t i = 0; i < numpart; i++)
{
- imo = fPyjets->K[2][i];
- if (imo > 1000) fPyjets->K[2][i] -= 1000;
+ imo = fPyjets->K[2][i];
+ if (imo > 1000) fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
}
-
+
+
+
+// delete[] kNew;
+// delete[] pNew;
+
} // end quench
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff