//
// ATLAS Tuning
//
- SetMSTP(51,7); // CTEQ5L pdf
+
+ SetMSTP(51, kCTEQ5L); // CTEQ5L pdf
SetMSTP(81,1); // Multiple Interactions ON
SetMSTP(82,4); // Double Gaussian Model
-void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t qTransport, Float_t maxLength, Int_t iECMethod)
+void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod)
{
// Initializes
// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
fQuenchingWeights = new AliQuenchingWeights();
fQuenchingWeights->InitMult();
- fQuenchingWeights->SetQTransport(qTransport);
+ fQuenchingWeights->SetK(k);
fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
- fQuenchingWeights->SetLengthMax(Int_t(maxLength));
- fQuenchingWeights->SampleEnergyLoss();
-
}
static Float_t eMean = 0.;
static Int_t icall = 0;
- Double_t p0[2][5];
- Double_t p1[2][5];
- Double_t p2[2][5];
- Int_t klast[2] = {-1, -1};
+ Double_t p0[4][5];
+ Double_t p1[4][5];
+ Double_t p2[4][5];
+ Int_t klast[4] = {-1, -1, -1, -1};
Int_t numpart = fPyjets->N;
- Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0.;
- Double_t pxq[2], pyq[2], pzq[2], eq[2], yq[2], mq[2], pq[2], phiq[2], thetaq[2], ptq[2];
- Bool_t quenched[2];
- Double_t phi;
- Double_t zInitial[2], wjtKick[2];
- Int_t nGluon[2];
-
+ Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
+ Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
+ Bool_t quenched[4];
+ Double_t zInitial[4], wjtKick[4];
+ Int_t nGluon[4];
+ Int_t qPdg[4];
Int_t imo, kst, pdg;
//
-// Primary partons
+// Sore information about Primary partons
+//
+// j =
+// 0, 1 partons from hard scattering
+// 2, 3 partons from initial state radiation
+//
+ for (Int_t i = 2; i <= 7; i++) {
+ Int_t j = 0;
+ // Skip gluons that participate in hard scattering
+ if (i == 4 || i == 5) continue;
+ // Gluons from hard Scattering
+ if (i == 6 || i == 7) {
+ j = i - 6;
+ pxq[j] = fPyjets->P[0][i];
+ pyq[j] = fPyjets->P[1][i];
+ pzq[j] = fPyjets->P[2][i];
+ eq[j] = fPyjets->P[3][i];
+ mq[j] = fPyjets->P[4][i];
+ } else {
+ // Gluons from initial state radiation
+ //
+ // Obtain 4-momentum vector from difference between original parton and parton after gluon
+ // radiation. Energy is calculated independently because initial state radition does not
+ // conserve strictly momentum and energy for each partonic system independently.
+ //
+ // Not very clean. Should be improved !
+ //
+ //
+ j = i;
+ pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
+ pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
+ pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
+ mq[j] = fPyjets->P[4][i];
+ eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
+ }
+//
+// Calculate some kinematic variables
//
-
-
-
- for (Int_t i = 6; i <= 7; i++) {
- Int_t j = i - 6;
-
- pxq[j] = fPyjets->P[0][i];
- pyq[j] = fPyjets->P[1][i];
- pzq[j] = fPyjets->P[2][i];
- eq[j] = fPyjets->P[3][i];
- mq[j] = fPyjets->P[4][i];
yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
- phi = phiq[j];
-
- // Quench only central jets
- if (TMath::Abs(yq[j]) > 2.5) {
+ qPdg[j] = fPyjets->K[1][i];
+ }
+
+ Double_t int0[4];
+ Double_t int1[4];
+
+ fGlauber->GetI0I1ForPythia(4, phiq, int0, int1, 15.);
+
+ for (Int_t j = 0; j < 4; j++) {
+ //
+ // Quench only central jets and with E > 10.
+ //
+
+
+ Int_t itype = (qPdg[j] == 21) ? 2 : 1;
+ Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
+
+ if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
zInitial[j] = 0.;
} else {
- pdg = fPyjets->K[1][i];
-
- // Get length in nucleus
- Double_t l;
- fGlauber->GetLengthsForPythia(1, &phi, &l, -1.);
- //
- // Energy loss for given length and parton typr
- Int_t itype = (pdg == 21) ? 2 : 1;
-
- Double_t eloss = fQuenchingWeights->GetELossRandom(itype, l, eq[j]);
- if (eq[j] > 80. && TMath::Abs(yq[j]) < 0.5) {
+ if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
icall ++;
eMean += eloss;
}
-
//
// Extra pt
-
- wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->GetQTransport());
+ Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
+ wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
//
// Fractional energy loss
zInitial[j] = eloss / eq[j];
if (zInitial[j] == 1.) zInitial[j] = 0.95;
//
// Some debug printing
- printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n",
- j, itype, eq[j], phi, l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
+
- zInitial[j] = 1.;
- while (zInitial[j] >= 0.95) zInitial[j] = gRandom->Exp(0.2);
+// printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n",
+// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
+
+// zInitial[j] = 0.8;
+// while (zInitial[j] >= 0.95) zInitial[j] = gRandom->Exp(0.2);
}
quenched[j] = (zInitial[j] > 0.01);
} // primary partons
-
+
Double_t pNew[1000][4];
Int_t kNew[1000];
Int_t icount = 0;
//
// System Loop
- for (Int_t isys = 0; isys < 2; isys++) {
+ for (Int_t isys = 0; isys < 4; isys++) {
// Skip to next system if not quenched.
if (!quenched[isys]) continue;
// Loop on radiation events
for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
- Double_t zHeavy = zInitial[isys];
-//
-
while (1) {
icount = 0;
for (Int_t k = 0; k < 4; k++)
// Quarks and gluons only
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 != (isys + 7) && (imo % 1000) != (isys + 7)) continue;
+ Int_t imom = imo % 1000;
+ if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
+ if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
+
// Skip comment lines
if (kst != 1 && kst != 2) continue;
theta = TMath::ATan2(pt, pz);
//
-// Save 4-momentum sum for balancing
- Int_t index = imo - 7;
- if (index >= 1000) index = imo % 1000 - 7;
+// Save 4-momentum sum for balancing
+ Int_t index = isys;
p0[index][0] += px;
p0[index][1] += py;
// Fractional energy loss
Double_t z = zInitial[index];
+
// Don't fully quench radiated gluons
//
if (imo > 1000) {
// This small factor makes sure that the gluons are not too close in phase space to avoid recombination
//
- z = 0.05;
+ z = 0.02;
}
+// printf("z: %d %f\n", imo, z);
+
//
-
- if (m > 0.) z = zHeavy;
//
//
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff