From: morsch Date: Wed, 24 Sep 2003 12:34:38 +0000 (+0000) Subject: Simple jet quenching algorythm added. X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=452af8c788be92fcf9c3d726d0453b34867512ab Simple jet quenching algorythm added. --- diff --git a/PYTHIA6/AliPythia.cxx b/PYTHIA6/AliPythia.cxx index 5613d528c7b..798e76a9e82 100644 --- a/PYTHIA6/AliPythia.cxx +++ b/PYTHIA6/AliPythia.cxx @@ -23,15 +23,20 @@ ClassImp(AliPythia) #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 &); //_____________________________________________________________________________ @@ -497,3 +502,288 @@ void AliPythia::Pycell(Int_t& njet) 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 + + diff --git a/PYTHIA6/AliPythia.h b/PYTHIA6/AliPythia.h index b0d7375b017..c91fa9872ef 100644 --- a/PYTHIA6/AliPythia.h +++ b/PYTHIA6/AliPythia.h @@ -51,9 +51,12 @@ class AliPythia : public TPythia6, public AliRndm virtual void SetDecayTable(); virtual void Pycell(Int_t& nclus); virtual void Pyclus(Int_t& nclus); + virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax); + virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez); + // return instance of the singleton static AliPythia* Instance(); - + virtual void Quench(); protected: Process_t fProcess; // Process type Float_t fEcms; // Centre of mass energy