+ numpart += ish;
+ (fPyjets->N) += ish;
+
+ if (ish == 1) {
+ fPyjets->P[0][iGlu] = p2[i][0];
+ fPyjets->P[1][iGlu] = p2[i][1];
+ fPyjets->P[2][iGlu] = p2[i][2];
+ fPyjets->P[3][iGlu] = p2[i][3];
+ fPyjets->P[4][iGlu] = p2[i][4];
+
+ fPyjets->K[0][iGlu] = 2;
+ fPyjets->K[1][iGlu] = 21;
+ fPyjets->K[2][iGlu] = fPyjets->K[2][iNew] + 1000;
+ fPyjets->K[3][iGlu] = -1;
+ fPyjets->K[4][iGlu] = -1;
+ } else {
+ //
+ // Split gluon in rest frame.
+ //
+ Double_t bx = p2[i][0] / p2[i][3];
+ Double_t by = p2[i][1] / p2[i][3];
+ Double_t bz = p2[i][2] / p2[i][3];
+ Double_t pst = p2[i][4] / 2.;
+ //
+ // Isotropic decay ????
+ Double_t cost = 2. * gRandom->Rndm() - 1.;
+ Double_t sint = TMath::Sqrt(1. - cost * cost);
+ Double_t phi = 2. * TMath::Pi() * gRandom->Rndm();
+
+ Double_t pz1 = pst * cost;
+ Double_t pz2 = -pst * cost;
+ Double_t pt1 = pst * sint;
+ Double_t pt2 = -pst * sint;
+ Double_t px1 = pt1 * TMath::Cos(phi);
+ Double_t py1 = pt1 * TMath::Sin(phi);
+ Double_t px2 = pt2 * TMath::Cos(phi);
+ Double_t py2 = pt2 * TMath::Sin(phi);
+
+ fPyjets->P[0][iGlu] = px1;
+ fPyjets->P[1][iGlu] = py1;
+ fPyjets->P[2][iGlu] = pz1;
+ fPyjets->P[3][iGlu] = pst;
+ fPyjets->P[4][iGlu] = 0.;
+
+ fPyjets->K[0][iGlu] = 2;
+ fPyjets->K[1][iGlu] = 21;
+ fPyjets->K[2][iGlu] = fPyjets->K[2][iNew] + 1000;
+ fPyjets->K[3][iGlu] = -1;
+ fPyjets->K[4][iGlu] = -1;
+
+ fPyjets->P[0][iGlu+1] = px2;
+ fPyjets->P[1][iGlu+1] = py2;
+ fPyjets->P[2][iGlu+1] = pz2;
+ fPyjets->P[3][iGlu+1] = pst;
+ fPyjets->P[4][iGlu+1] = 0.;
+
+ fPyjets->K[0][iGlu+1] = 2;
+ fPyjets->K[1][iGlu+1] = 21;
+ fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew] + 1000;
+ fPyjets->K[3][iGlu+1] = -1;
+ fPyjets->K[4][iGlu+1] = -1;
+ SetMSTU(1,0);
+ SetMSTU(2,0);
+ //
+ // Boost back
+ //
+ Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
+ }
+ } // end adding gluons
+ //
+ // Check energy conservation
+ Double_t pxs = 0.;
+ Double_t pys = 0.;
+ Double_t pzs = 0.;
+ Double_t es = 14000.;
+
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ kst = fPyjets->K[0][i];
+ if (kst != 1 && kst != 2) continue;
+ pxs += fPyjets->P[0][i];
+ pys += fPyjets->P[1][i];
+ pzs += fPyjets->P[2][i];
+ es -= fPyjets->P[3][i];
+ }
+ if (TMath::Abs(pxs) > 1.e-2 ||
+ TMath::Abs(pys) > 1.e-2 ||
+ 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");