#ifndef WIN32
# define pyclus pyclus_
# define pycell pycell_
+# define pyshow pyshow_
+# define pyrobo pyrobo_
# define type_of_call
#else
# define pyclus PYCLUS
# define pycell PYCELL
+# define pyrobo PYROBO
# define type_of_call _stdcall
#endif
extern "C" void type_of_call pyclus(Int_t & );
extern "C" void type_of_call pycell(Int_t & );
+extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
//_____________________________________________________________________________
pycell(njet);
}
+void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
+{
+// Call Pythia jet reconstruction algorithm
+//
+ Int_t numpart = fPyjets->N;
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ if (fPyjets->K[2][i] == 7) ip1 = i+1;
+ if (fPyjets->K[2][i] == 8) ip2 = i+1;
+ }
+
+
+ qmax = 2. * GetVINT(51);
+ printf("Pyshow %d %d %f", ip1, ip2, qmax);
+
+ pyshow(ip1, ip2, qmax);
+}
+
+void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
+{
+ pyrobo(imi, ima, the, phi, bex, bey, bez);
+}
+
+
+
+void AliPythia::Quench()
+{
+//
+//
+// Simple Jet Quenching routine:
+// =============================
+// The jet formed by all final state partons radiated by the parton created
+// in the hard collisions is quenched by a factor z using:
+// (E + p_z)new = (1-z) (E + p_z)old
+//
+// The lost momentum is first balanced by one gluon with virtuality > 0.
+// Subsequently the gluon splits to yield two gluons with E = p.
+//
+ Float_t p0[2][5];
+ Float_t p1[2][5];
+ Float_t p2[2][5];
+ Int_t klast[2] = {-1, -1};
+ Int_t kglu[2];
+ for (Int_t i = 0; i < 4; i++)
+ {
+ p0[0][i] = 0.;
+ p0[1][i] = 0.;
+ p1[0][i] = 0.;
+ p1[1][i] = 0.;
+ p2[0][i] = 0.;
+ p2[1][i] = 0.;
+ }
+
+ Int_t numpart = fPyjets->N;
+
+ for (Int_t i = 0; i < numpart; i++)
+ {
+ Int_t imo = fPyjets->K[2][i];
+ Int_t kst = fPyjets->K[0][i];
+ Int_t pdg = fPyjets->K[1][i];
+
+// Quarks and gluons only
+ if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
+
+// Particles from hard scattering only
+
+
+ Float_t px = fPyjets->P[0][i];
+ Float_t py = fPyjets->P[1][i];
+ Float_t pz = fPyjets->P[2][i];
+ Float_t e = fPyjets->P[3][i];
+ Float_t m = fPyjets->P[4][i];
+ Float_t pt = TMath::Sqrt(px * px + py * py);
+// Skip comment lines
+ if (kst != 1 && kst != 2) continue;
+
+ Float_t mt = TMath::Sqrt(px * px + py * py + m * m);
+
+ //
+ // Some cuts to be in a save kinematic region
+ //
+ if (imo != 7 && imo != 8) continue;
+ Int_t index = imo - 7;
+ klast[index] = i;
+
+ p0[index][0] += px;
+ p0[index][1] += py;
+ p0[index][2] += pz;
+ p0[index][3] += e;
+
+ //
+ // Fix z
+ //
+
+ Float_t z = 0.2;
+ Float_t eppzOld = e + pz;
+ Float_t empzOld = e - pz;
+
+
+ //
+ // Kinematics of the original parton
+ //
+
+ Float_t eppzNew = (1. - z) * eppzOld;
+ Float_t empzNew = empzOld - mt * mt * z / eppzOld;
+ Float_t eNew0 = 0.5 * (eppzNew + empzNew);
+ Float_t pzNew0 = 0.5 * (eppzNew - empzNew);
+ //
+ // Skip if pt too small
+ //
+ if (m * m > eppzNew * empzNew) continue;
+ Float_t ptNew = TMath::Sqrt(eppzNew * empzNew - m * m);
+ Float_t pxNew0 = ptNew / pt * px;
+ Float_t pyNew0 = ptNew / pt * py;
+
+ 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;
+
+ }
+
+ //
+ // Gluons
+ //
+
+ for (Int_t k = 0; k < 2; k++)
+ {
+ for (Int_t j = 0; j < 4; j++)
+ {
+ p2[k][j] = p0[k][j] - p1[k][j];
+ }
+ 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.)
+ {
+
+ //
+ // Bring gluon back to mass shell via momentum scaling
+ // (momentum will not be conserved, but energy)
+ //
+ // not used anymore
+/*
+ Float_t psq = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
+ Float_t fact = TMath::Sqrt(1. + p2[k][4] / psq);
+ p2[k][0] *= fact;
+ p2[k][1] *= fact;
+ p2[k][2] *= fact;
+ p2[k][3] = TMath::Sqrt(psq) * fact;
+ p2[k][4] = 0.;
+*/
+ }
+ }
+
+ if (p2[0][4] > 0.) {
+ p2[0][4] = TMath::Sqrt(p2[0][4]);
+ } else {
+ printf("Warning negative mass squared ! \n");
+ }
+
+ if (p2[1][4] > 0.) {
+ p2[1][4] = TMath::Sqrt(p2[1][4]);
+ } else {
+ printf("Warning negative mass squared ! \n");
+ }
+
+ //
+ // Add the gluons
+ //
+
+
+ for (Int_t i = 0; i < 2; i++) {
+ Int_t ish, jmin, jmax, iGlu, iNew;
+ Int_t in = klast[i];
+ ish = 0;
+
+ if (in == -1) continue;
+ if (i == 1 && klast[1] > klast[0]) in += ish;
+
+ jmin = in - 1;
+ ish = 1;
+
+ if (p2[i][4] > 0) ish = 2;
+
+ iGlu = in;
+ iNew = in + ish;
+ jmax = numpart + ish - 1;
+
+ if (fPyjets->K[0][in-1] == 1 || fPyjets->K[0][in-1] == 21 || fPyjets->K[0][in-1] == 11) {
+ jmin = in;
+ iGlu = in + 1;
+ iNew = in;
+ }
+
+ kglu[i] = iGlu;
+
+ for (Int_t j = jmax; j > jmin; j--)
+ {
+ for (Int_t k = 0; k < 5; k++) {
+ fPyjets->K[k][j] = fPyjets->K[k][j-ish];
+ fPyjets->P[k][j] = fPyjets->P[k][j-ish];
+ fPyjets->V[k][j] = fPyjets->V[k][j-ish];
+ }
+ } // end shifting
+ 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];
+ 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];
+
+ Float_t pst = p2[i][4] / 2.;
+
+ Float_t cost = 2. * gRandom->Rndm() - 1.;
+ Float_t sint = TMath::Sqrt(1. - cost * cost);
+ Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
+
+ Float_t pz1 = pst * cost;
+ Float_t pz2 = -pst * cost;
+ Float_t pt1 = pst * sint;
+ Float_t pt2 = -pst * sint;
+ Float_t px1 = pt1 * TMath::Cos(phi);
+ Float_t py1 = pt1 * TMath::Sin(phi);
+ Float_t px2 = pt2 * TMath::Cos(phi);
+ Float_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];
+ 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];
+ 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
+} // end quench
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff