# define pycell pycell_
# define pyshow pyshow_
# define pyrobo pyrobo_
+# define pyquen pyquen_
# define type_of_call
#else
# define pyclus PYCLUS
# define pycell PYCELL
# define pyrobo PYROBO
+# define pyquen PYQUEN
# define type_of_call _stdcall
#endif
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 &);
+extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
//_____________________________________________________________________________
fProcess = process;
fEcms = energy;
fStrucFunc = strucfunc;
-// don't decay p0
- SetMDCY(Pycomp(111),1,0);
-// select structure function
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+ SetMDCY(Pycomp(111) ,1,0);
+ SetMDCY(Pycomp(310) ,1,0);
+ SetMDCY(Pycomp(3122),1,0);
+ SetMDCY(Pycomp(3112),1,0);
+ SetMDCY(Pycomp(3212),1,0);
+ SetMDCY(Pycomp(3222),1,0);
+ SetMDCY(Pycomp(3312),1,0);
+ SetMDCY(Pycomp(3322),1,0);
+ SetMDCY(Pycomp(3334),1,0);
+ // select structure function
SetMSTP(52,2);
SetMSTP(51,strucfunc);
//
//
// ATLAS Tuning
//
- SetMSTP(51,7); // CTEQ5L pdf
+
+ SetMSTP(51, kCTEQ5L); // CTEQ5L pdf
SetMSTP(81,1); // Multiple Interactions ON
SetMSTP(82,4); // Double Gaussian Model
// ATLAS Tuning
//
- SetMSTP(51,7); // CTEQ5L pdf
+ SetMSTP(51,kCTEQ5L); // CTEQ5L pdf
SetMSTP(81,1); // Multiple Interactions ON
SetMSTP(82,4); // Double Gaussian Model
// Set c-quark mass
SetPMAS(4,1,1.2);
+ break;
+ case kPyDPlusPbPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with Pb-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+
+ // All QCD processes
+ SetMSEL(1);
+
+ // No multiple interactions
+ SetMSTP(81,0);
+ SetPARP(81,0.0);
+ SetPARP(82,0.0);
+
+ // Initial/final parton shower on (Pythia default)
+ SetMSTP(61,1);
+ SetMSTP(71,1);
+
+ // 2nd order alpha_s
+ SetMSTP(2,2);
+
+ // QCD scales
+ SetMSTP(32,2);
+ SetPARP(34,1.0);
+
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.304);
+ SetPARP(93,6.52);
+
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+
break;
case kPyCharmpPbMNR:
case kPyD0pPbMNR:
// Set c-quark mass
SetPMAS(4,1,1.2);
+ break;
+ case kPyDPluspPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with p-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+
+ // All QCD processes
+ SetMSEL(1);
+
+ // No multiple interactions
+ SetMSTP(81,0);
+ SetPARP(81,0.0);
+ SetPARP(82,0.0);
+
+ // Initial/final parton shower on (Pythia default)
+ SetMSTP(61,1);
+ SetMSTP(71,1);
+
+ // 2nd order alpha_s
+ SetMSTP(2,2);
+
+ // QCD scales
+ SetMSTP(32,2);
+ SetPARP(34,1.0);
+
+ // Intrinsic <kT>
+ SetMSTP(91,1);
+ SetPARP(91,1.16);
+ SetPARP(93,5.8);
+
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+
break;
case kPyCharmppMNR:
case kPyD0ppMNR:
// Set c-quark mass
SetPMAS(4,1,1.2);
+ break;
+ case kPyDPlusppMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+
+ // All QCD processes
+ SetMSEL(1);
+
+ // No multiple interactions
+ SetMSTP(81,0);
+ SetPARP(81,0.0);
+ SetPARP(82,0.0);
+
+ // Initial/final parton shower on (Pythia default)
+ SetMSTP(61,1);
+ SetMSTP(71,1);
+
+ // 2nd order alpha_s
+ SetMSTP(2,2);
+
+ // QCD scales
+ SetMSTP(32,2);
+ SetPARP(34,1.0);
+
+ // Intrinsic <kT^2>
+ SetMSTP(91,1);
+ SetPARP(91,1.);
+ SetPARP(93,5.);
+
+ // Set c-quark mass
+ SetPMAS(4,1,1.2);
+
break;
case kPyBeautyPbPbMNR:
// Tuning of Pythia parameters aimed to get a resonable agreement
-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};
- Int_t kglu[2];
+ 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 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) {
- zInitial[j] = 0.;
+ qPdg[j] = fPyjets->K[1][i];
+ }
+
+ Double_t int0[4];
+ Double_t int1[4];
+
+ fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 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.) {
+ fZQuench[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];
+ fZQuench[j] = eloss / eq[j];
//
// Avoid complete loss
//
- if (zInitial[j] == 1.) zInitial[j] = 0.95;
+ if (fZQuench[j] == 1.) fZQuench[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]);
+
+// fZQuench[j] = 0.8;
+// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2);
}
- quenched[j] = (zInitial[j] > 0.01);
+ quenched[j] = (fZQuench[j] > 0.01);
} // primary partons
-
+
+
+
Double_t pNew[1000][4];
Int_t kNew[1000];
Int_t icount = 0;
+ Double_t zquench[4];
+
//
// 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;
- nGluon[isys] = 1 + Int_t(zInitial[isys] / (1. - zInitial[isys]));
+ nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
if (nGluon[isys] > 6) nGluon[isys] = 6;
- zInitial[isys] = 1. - TMath::Power(1. - zInitial[isys], 1./Double_t(nGluon[isys]));
+ zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
// 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];
+ Double_t z = zquench[index];
+
// Don't fully quench radiated gluons
//
// 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;
//
//
p2[isys][4] = TMath::Sqrt(p2[isys][4]);
break;
} else {
- printf("Warning negative mass squared in system %d %f ! \n", isys, zInitial[isys]);
+ printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
printf("4-Momentum: %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 !\n");
// Add the gluons
//
Int_t ish = 0;
- Int_t iGlu, jmin, jmax, iNew;
+ Int_t iGlu;
if (!quenched[isys]) continue;
//
// Last parton from shower i
}
// this->Pylist(1);
} // end quench
+
+
+void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
+{
+ // Igor Lokthine's quenching routine
+ pyquen(a, ibf, b);
+}
+
+void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
+{
+ // Return event specific quenching parameters
+ xp = fXJet;
+ yp = fYJet;
+ for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
+
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff