#include "AliPythia6.h"
#include "AliStack.h"
#include "AliPythiaRndm.h"
#include "AliFastGlauber.h"
#include "AliQuenchingWeights.h"
#include "TVector3.h"
#include "TParticle.h"
#include "PyquenCommon.h"
ClassImp(AliPythia6)
#ifndef WIN32
# define pyclus pyclus_
# define pycell pycell_
# define pyshow pyshow_
# define pyshowq pyshowq_
# define pyrobo pyrobo_
# define pyquen pyquen_
# define pyevnw pyevnw_
# define pyjoin pyjoin_
# define qpygin0 qpygin0_
# define setpowwght setpowwght_
# define type_of_call
#else
# define pyclus PYCLUS
# define pycell PYCELL
# define pyshow PYSHOW
# define pyshowq PYSHOWQ
# define pyrobo PYROBO
# define pyquen PYQUEN
# define pyevnw PYEVNW
# define pyjoin PYJOIN
# define qpygin0 QPYGIN0
# define setpowwght SETPOWWGHT
# define type_of_call _stdcall
#endif
extern "C" void type_of_call pyjoin(Int_t &, Int_t * );
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 pyshowq(Int_t &, Int_t &, Double_t &);
extern "C" void type_of_call qpygin0();
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 &);
extern "C" void type_of_call pyevnw();
extern "C" void type_of_call setpowwght(Double_t &);
AliPythia6* AliPythia6::fgAliPythia=NULL;
AliPythia6::AliPythia6():
TPythia6(),
AliPythiaBase(),
fProcess(kPyMb),
fEcms(0.),
fStrucFunc(kCTEQ5L),
fProjectile("p"),
fTarget("p"),
fXJet(0.),
fYJet(0.),
fNGmax(30),
fZmax(0.97),
fGlauber(0),
fQuenchingWeights(0)
{
Int_t i;
for (i = 0; i < 501; i++) fDefMDCY[i] = 0;
for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
for (i = 0; i < 4; i++) fZQuench[i] = 0;
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
fGlauber = 0;
fQuenchingWeights = 0;
}
AliPythia6::AliPythia6(const AliPythia6& pythia):
TPythia6(),
AliPythiaBase(),
fProcess(kPyMb),
fEcms(0.),
fStrucFunc(kCTEQ5L),
fProjectile("p"),
fTarget("p"),
fXJet(0.),
fYJet(0.),
fNGmax(30),
fZmax(0.97),
fGlauber(0),
fQuenchingWeights(0)
{
Int_t i;
for (i = 0; i < 501; i++) fDefMDCY[i] = 0;
for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
for (i = 0; i < 4; i++) fZQuench[i] = 0;
pythia.Copy(*this);
}
void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t )
{
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
fProcess = process;
fEcms = energy;
fStrucFunc = strucfunc;
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);
SetMSTP(52,2);
SetMSTP(51,AliStructFuncType::PDFsetIndex(strucfunc));
SetMSTU(16,2);
for (Int_t i=1; i<= 200; i++) {
SetMSUB(i,0);
}
switch (process)
{
case kPyOldUEQ2ordered:
SetMSTP(81,1);
SetMSTP(82,4);
SetPARP(83,0.5);
SetPARP(84,0.4);
SetPARP(82,2.0);
SetPARP(89,1800);
SetPARP(90,0.25);
SetPARP(85,0.9);
SetPARP(86,0.95);
SetPARP(67,4);
SetPARP(71,4);
SetPARJ(81,0.29);
break;
case kPyOldUEQ2ordered2:
SetMSTP(81,1);
SetMSTP(82,4);
SetPARP(83,0.5);
SetPARP(84,0.4);
SetPARP(82,2.0);
SetPARP(89,1800);
SetPARP(90,0.16);
SetPARP(85,0.9);
SetPARP(86,0.95);
SetPARP(67,4);
SetPARP(71,4);
SetPARJ(81,0.29);
break;
case kPyOldPopcorn:
SetMSEL(1);
SetMSTJ(12,3);
SetMSTP(88,2);
SetMSTJ(1,1);
AtlasTuning();
break;
case kPyCharm:
SetMSEL(4);
SetPMAS(4,1,1.2);
SetMSTP(91,1);
SetPARP(91,1.);
SetPARP(93,5.);
break;
case kPyBeauty:
SetMSEL(5);
SetPMAS(5,1,4.75);
break;
case kPyJpsi:
SetMSEL(0);
SetMSUB(86,1);
break;
case kPyJpsiChi:
SetMSEL(0);
SetMSUB(86,1);
SetMSUB(87,1);
SetMSUB(88,1);
SetMSUB(89,1);
break;
case kPyCharmUnforced:
SetMSEL(0);
SetMSUB(28,1);
SetMSUB(53,1);
SetMSUB(68,1);
break;
case kPyBeautyUnforced:
SetMSEL(0);
SetMSUB(28,1);
SetMSUB(53,1);
SetMSUB(68,1);
break;
case kPyMb:
SetMSEL(0);
SetMSUB(92,1);
SetMSUB(93,1);
SetMSUB(94,1);
SetMSUB(95,1);
AtlasTuning();
break;
case kPyMbAtlasTuneMC09:
SetMSEL(0);
SetMSUB(92,1);
SetMSUB(93,1);
SetMSUB(94,1);
SetMSUB(95,1);
AtlasTuningMC09();
break;
case kPyMbWithDirectPhoton:
SetMSEL(0);
SetMSUB(92,1);
SetMSUB(93,1);
SetMSUB(94,1);
SetMSUB(95,1);
SetMSUB(14,1);
SetMSUB(18,1);
SetMSUB(29,1);
SetMSUB(114,1);
SetMSUB(115,1);
AtlasTuning();
break;
case kPyMbDefault:
SetMSEL(0);
SetMSUB(92,1);
SetMSUB(93,1);
SetMSUB(94,1);
SetMSUB(95,1);
break;
case kPyLhwgMb:
SetMSEL(0);
SetMSUB(92,1);
SetMSUB(93,1);
SetMSUB(94,1);
SetMSUB(95,1);
SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));
SetMSTP(52,2);
SetMSTP(68,1);
SetMSTP(70,2);
SetMSTP(81,1);
SetMSTP(82,4);
SetMSTP(88,1);
SetPARP(82,2.3);
SetPARP(83,0.5);
SetPARP(84,0.5);
SetPARP(85,0.9);
SetPARP(90,0.2);
break;
case kPyMbNonDiffr:
SetMSEL(0);
SetMSUB(95,1);
AtlasTuning();
break;
case kPyMbMSEL1:
ConfigHeavyFlavor();
SetMSTP(91,1);
SetPARP(91,1.);
SetPARP(93,5.);
SetPMAS(4,1,1.2);
SetPMAS(5,1,4.78);
SetPARP(71,4.);
AtlasTuning();
break;
case kPyJets:
SetMSEL(1);
SetPARP(67,2.5);
SetMSTP(82,4);
SetPARP(82,2.0);
SetPARP(84,0.4);
SetPARP(85,0.90) ;
SetPARP(86,0.95);
SetPARP(89,1800.);
SetPARP(90,0.25);
break;
case kPyDirectGamma:
SetMSEL(10);
break;
case kPyCharmPbPbMNR:
case kPyD0PbPbMNR:
case kPyDPlusPbPbMNR:
case kPyDPlusStrangePbPbMNR:
ConfigHeavyFlavor();
SetMSTP(91,1);
SetPARP(91,1.304);
SetPARP(93,6.52);
SetPMAS(4,1,1.2);
break;
case kPyCharmpPbMNR:
case kPyD0pPbMNR:
case kPyDPluspPbMNR:
case kPyDPlusStrangepPbMNR:
ConfigHeavyFlavor();
SetMSTP(91,1);
SetPARP(91,1.16);
SetPARP(93,5.8);
SetPMAS(4,1,1.2);
break;
case kPyCharmppMNR:
case kPyD0ppMNR:
case kPyDPlusppMNR:
case kPyDPlusStrangeppMNR:
case kPyLambdacppMNR:
ConfigHeavyFlavor();
SetMSTP(91,1);
SetPARP(91,1.);
SetPARP(93,5.);
SetPMAS(4,1,1.2);
break;
case kPyCharmppMNRwmi:
ConfigHeavyFlavor();
SetMSTP(91,1);
SetPARP(91,1.);
SetPARP(93,5.);
SetPMAS(4,1,1.2);
AtlasTuning();
break;
case kPyBeautyPbPbMNR:
ConfigHeavyFlavor();
SetPARP(67,1.0);
SetPARP(71,1.0);
SetMSTP(91,1);
SetPARP(91,2.035);
SetPARP(93,10.17);
SetPMAS(5,1,4.75);
break;
case kPyBeautypPbMNR:
ConfigHeavyFlavor();
SetPARP(67,1.0);
SetPARP(71,1.0);
SetMSTP(91,1);
SetPARP(91,1.60);
SetPARP(93,8.00);
SetPMAS(5,1,4.75);
break;
case kPyBeautyppMNR:
ConfigHeavyFlavor();
SetPARP(67,1.0);
SetPARP(71,1.0);
SetMSTP(91,1);
SetPARP(91,1.);
SetPARP(93,5.);
SetPMAS(5,1,4.75);
break;
case kPyBeautyJets:
case kPyBeautyppMNRwmi:
ConfigHeavyFlavor();
SetPARP(67,1.0);
SetPARP(71,1.0);
SetMSTP(91,1);
SetPARP(91,1.);
SetPARP(93,5.);
SetPMAS(5,1,4.75);
AtlasTuning();
break;
case kPyW:
SetMSEL(0);
SetMSUB(2,1);
SetMSTP(61,1);
SetMSTP(71,1);
break;
case kPyZ:
SetMSEL(0);
SetMSUB(1,1);
SetMSTP(43,2);
SetMSTP(61,1);
SetMSTP(71,1);
break;
case kPyZgamma:
SetMSEL(0);
SetMSUB(1,1);
SetMSTP(61,1);
SetMSTP(71,1);
break;
case kPyMBRSingleDiffraction:
case kPyMBRDoubleDiffraction:
case kPyMBRCentralDiffraction:
break;
case kPyJetsPWHG:
SetMSTP(86,1);
SetMSTU(22,10);
SetMSTU(26,20);
break;
case kPyCharmPWHG:
case kPyBeautyPWHG:
case kPyWPWHG:
SetMSTU(26,20);
}
SetMSTP(41,1);
if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG || process == kPyWPWHG) {
Initialize("USER","","",0.);
} else {
Initialize("CMS",fProjectile,fTarget,fEcms);
}
}
Int_t AliPythia6::CheckedLuComp(Int_t kf)
{
Int_t kc=Pycomp(kf);
return kc;
}
void AliPythia6::SetNuclei(Int_t a1, Int_t a2)
{
SetMSTP(52,2);
SetMSTP(192, a1);
SetMSTP(193, a2);
}
AliPythia6* AliPythia6::Instance()
{
if (fgAliPythia) {
return fgAliPythia;
} else {
fgAliPythia = new AliPythia6();
return fgAliPythia;
}
}
void AliPythia6::PrintParticles()
{
Int_t np = 0;
char* name = new char[16];
for (Int_t kf=0; kf<1000000; kf++) {
for (Int_t c = 1; c > -2; c-=2) {
Int_t kc = Pycomp(c*kf);
if (kc) {
Float_t mass = GetPMAS(kc,1);
Float_t width = GetPMAS(kc,2);
Float_t tau = GetPMAS(kc,4);
Pyname(kf,name);
np++;
printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
c*kf, name, mass, width, tau);
}
}
}
printf("\n Number of particles %d \n \n", np);
}
void AliPythia6::ResetDecayTable()
{
Int_t i;
for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
}
void AliPythia6::SetDecayTable()
{
Int_t i;
for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
}
void AliPythia6::Pyjoin(Int_t& npart, Int_t *ipart)
{
pyjoin(npart, ipart);
}
void AliPythia6::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
{
pyshowq(ip1, ip2, qmax);
}
void AliPythia6::Qpygin0()
{
qpygin0();
}
void AliPythia6::Pyclus(Int_t& njet)
{
pyclus(njet);
}
void AliPythia6::Pycell(Int_t& njet)
{
pycell(njet);
}
void AliPythia6::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
{
Int_t n = GetN();
px = GetPyjets()->P[0][n+i];
py = GetPyjets()->P[1][n+i];
pz = GetPyjets()->P[2][n+i];
e = GetPyjets()->P[3][n+i];
}
void AliPythia6::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
{
pyshow(ip1, ip2, qmax);
}
void AliPythia6::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 AliPythia6::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
{
fGlauber = AliFastGlauber::Instance();
fGlauber->Init(2);
fGlauber->SetCentralityClass(cMin, cMax);
fQuenchingWeights = new AliQuenchingWeights();
fQuenchingWeights->InitMult();
fQuenchingWeights->SetK(k);
fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
fNGmax = ngmax;
fZmax = zmax;
}
void AliPythia6::Quench()
{
static Float_t eMean = 0.;
static Int_t icall = 0;
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., 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] = {0., 0., 0., 0.};
Int_t nGluon[4];
Int_t qPdg[4];
Int_t imo, kst, pdg;
for (Int_t i = 2; i <= 7; i++) {
Int_t j = 0;
if (i == 4 || i == 5) continue;
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 {
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]);
}
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]);
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++) {
Int_t itype = (qPdg[j] == 21) ? 2 : 1;
Double_t eloss = fQuenchingWeights->GetELossRandomK(itype, int0[j], int1[j], eq[j]);
if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
fZQuench[j] = 0.;
} else {
if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
icall ++;
eMean += eloss;
}
Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
fZQuench[j] = eloss / eq[j];
if (fZQuench[j] > fZmax) fZQuench[j] = fZmax;
}
quenched[j] = (fZQuench[j] > 0.01);
}
Double_t pNew[1000][4];
Int_t kNew[1000];
Int_t icount = 0;
Double_t zquench[4];
for (Int_t isys = 0; isys < 4; isys++) {
if (!quenched[isys]) continue;
nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
Int_t igMin = -1;
Int_t igMax = -1;
Double_t pg[4] = {0., 0., 0., 0.};
for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
while (1) {
icount = 0;
for (Int_t k = 0; k < 4; k++)
{
p0[isys][k] = 0.;
p1[isys][k] = 0.;
p2[isys][k] = 0.;
}
for (Int_t i = 0; i < numpart; i++)
{
imo = fPyjets->K[2][i];
kst = fPyjets->K[0][i];
pdg = fPyjets->K[1][i];
if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
Int_t imom = imo % 1000;
if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
if (kst != 1 && kst != 2) continue;
px = fPyjets->P[0][i];
py = fPyjets->P[1][i];
pz = fPyjets->P[2][i];
e = fPyjets->P[3][i];
m = fPyjets->P[4][i];
pt = TMath::Sqrt(px * px + py * py);
p = TMath::Sqrt(px * px + py * py + pz * pz);
phi = TMath::Pi() + TMath::ATan2(-py, -px);
theta = TMath::ATan2(pt, pz);
Int_t index = isys;
p0[index][0] += px;
p0[index][1] += py;
p0[index][2] += pz;
p0[index][3] += e;
klast[index] = i;
Double_t z = zquench[index];
if (imo > 1000) {
z = 0.02;
}
TVector3 v(px, py, pz);
v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
Double_t mt2 = jt * jt + m * m;
Double_t zmax = 1.;
if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
Double_t eppzOld = e + pl;
Double_t empzOld = e - pl;
Double_t eppzNew = (1. - z) * eppzOld;
Double_t empzNew = empzOld - mt2 * z / eppzOld;
Double_t eNew = 0.5 * (eppzNew + empzNew);
Double_t plNew = 0.5 * (eppzNew - empzNew);
Double_t jtNew;
Double_t mt2New = eppzNew * empzNew;
if (mt2New < 1.e-8) mt2New = 0.;
if (z < zmax) {
if (m * m > mt2New) {
Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
jtNew = 0;
} else {
jtNew = TMath::Sqrt(mt2New - m * m);
}
} else {
jtNew = jt;
eNew = TMath::Sqrt(plNew * plNew + mt2);
}
Double_t pxNew = 0;
Double_t pyNew = 0;
if (jt > 0.) {
pxNew = jtNew / jt * pxs;
pyNew = jtNew / jt * pys;
}
TVector3 w(pxNew, pyNew, plNew);
w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
p1[index][0] += pxNew;
p1[index][1] += pyNew;
p1[index][2] += plNew;
p1[index][3] += eNew;
pNew[icount][0] = pxNew;
pNew[icount][1] = pyNew;
pNew[icount][2] = plNew;
pNew[icount][3] = eNew;
kNew[icount] = i;
icount++;
}
if (icount == 0) quenched[isys] = kFALSE;
if (!quenched[isys]) break;
for (Int_t j = 0; j < 4; j++)
{
p2[isys][j] = p0[isys][j] - p1[isys][j];
}
p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
if (p2[isys][4] > 0.) {
p2[isys][4] = TMath::Sqrt(p2[isys][4]);
break;
} else {
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");
p2[isys][4] = 0.;
p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
break;
} else {
p2[isys][4] = 0.;
break;
}
}
}
for (Int_t k = 0; k < icount; k++) {
fPyjets->P[0][kNew[k]] = pNew[k][0];
fPyjets->P[1][kNew[k]] = pNew[k][1];
fPyjets->P[2][kNew[k]] = pNew[k][2];
fPyjets->P[3][kNew[k]] = pNew[k][3];
}
Int_t ish = 0;
Int_t iGlu;
if (!quenched[isys]) continue;
Int_t in = klast[isys];
if (in == -1) continue;
if (isys == 1 && klast[1] > klast[0]) in += ish;
ish = 1;
if (p2[isys][4] > 0.05) ish = 2;
iGlu = numpart;
if (iglu == 0) igMin = iGlu;
igMax = iGlu;
numpart += ish;
(fPyjets->N) += ish;
if (ish == 1) {
fPyjets->P[0][iGlu] = p2[isys][0];
fPyjets->P[1][iGlu] = p2[isys][1];
fPyjets->P[2][iGlu] = p2[isys][2];
fPyjets->P[3][iGlu] = p2[isys][3];
fPyjets->P[4][iGlu] = p2[isys][4];
fPyjets->K[0][iGlu] = 1;
if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
fPyjets->K[1][iGlu] = 21;
fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
fPyjets->K[3][iGlu] = -1;
fPyjets->K[4][iGlu] = -1;
pg[0] += p2[isys][0];
pg[1] += p2[isys][1];
pg[2] += p2[isys][2];
pg[3] += p2[isys][3];
} else {
Double_t bx = p2[isys][0] / p2[isys][3];
Double_t by = p2[isys][1] / p2[isys][3];
Double_t bz = p2[isys][2] / p2[isys][3];
Double_t pst = p2[isys][4] / 2.;
Double_t cost = 2. * gRandom->Rndm() - 1.;
Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
Double_t phis = 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(phis);
Double_t py1 = pt1 * TMath::Sin(phis);
Double_t px2 = pt2 * TMath::Cos(phis);
Double_t py2 = pt2 * TMath::Sin(phis);
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] = 1 ;
fPyjets->K[1][iGlu] = 21;
fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 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] = 1;
if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
fPyjets->K[1][iGlu+1] = 21;
fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
fPyjets->K[3][iGlu+1] = -1;
fPyjets->K[4][iGlu+1] = -1;
SetMSTU(1,0);
SetMSTU(2,0);
Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
}
}
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);
}
}
for (Int_t i = 0; i < numpart; i++)
{
imo = fPyjets->K[2][i];
if (imo > 1000) {
fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
}
}
}
void AliPythia6::Pyquen(Double_t a, Int_t ibf, Double_t b)
{
pyquen(a, ibf, b);
}
void AliPythia6::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
{
PYQPAR.t0 = t0;
PYQPAR.tau0 = tau0;
PYQPAR.nf = nf;
PYQPAR.iengl = iengl;
PYQPAR.iangl = iangl;
}
void AliPythia6::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
{
Int_t npart = stack -> GetNprimary();
Int_t n0 = 0;
if (!flag) {
GetPyjets()->N = npart;
} else {
n0 = GetPyjets()->N;
GetPyjets()->N = n0 + npart;
}
for (Int_t part = 0; part < npart; part++) {
TParticle *mPart = stack->Particle(part);
Int_t kf = mPart->GetPdgCode();
Int_t ks = mPart->GetStatusCode();
Int_t idf = mPart->GetFirstDaughter();
Int_t idl = mPart->GetLastDaughter();
if (reHadr) {
if (ks == 11 || ks == 12) {
ks -= 10;
idf = -1;
idl = -1;
}
}
Float_t px = mPart->Px();
Float_t py = mPart->Py();
Float_t pz = mPart->Pz();
Float_t e = mPart->Energy();
Float_t m = mPart->GetCalcMass();
(GetPyjets())->P[0][part+n0] = px;
(GetPyjets())->P[1][part+n0] = py;
(GetPyjets())->P[2][part+n0] = pz;
(GetPyjets())->P[3][part+n0] = e;
(GetPyjets())->P[4][part+n0] = m;
(GetPyjets())->K[1][part+n0] = kf;
(GetPyjets())->K[0][part+n0] = ks;
(GetPyjets())->K[3][part+n0] = idf + 1;
(GetPyjets())->K[4][part+n0] = idl + 1;
(GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
}
}
void AliPythia6::Pyevnw()
{
pyevnw();
}
void AliPythia6::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
{
xp = fXJet;
yp = fYJet;
for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
}
void AliPythia6::ConfigHeavyFlavor()
{
SetMSEL(1);
SetMSTP(81,0);
SetPARP(81, 0.);
SetPARP(82, 0.);
SetMSTP(61,1);
SetMSTP(71,1);
SetMSTP(2,2);
SetMSTP(32,2);
SetPARP(34,1.0);
}
void AliPythia6::AtlasTuning()
{
SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ5L));
SetMSTP(81,1);
SetMSTP(82,4);
SetPARP(81,1.9);
SetPARP(82,1.8);
SetPARP(89,1000.);
SetPARP(90,0.16);
SetPARP(83,0.5);
SetPARP(84,0.5);
SetPARP(85,0.33);
SetPARP(86,0.66);
SetPARP(67,1);
}
void AliPythia6::AtlasTuningMC09()
{
printf("ATLAS New TUNE MC09\n");
SetMSTP(81,21);
SetMSTP(82, 4);
SetMSTP(52, 2);
SetMSTP(51, 20650);
SetMSTP(70, 0);
SetMSTP(72, 1);
SetMSTP(88, 1);
SetMSTP(90, 0);
SetPARP(78, 0.3);
SetPARP(80, 0.1);
SetPARP(82, 2.3);
SetPARP(83, 0.8);
SetPARP(84, 0.7);
SetPARP(90, 0.25);
SetPARJ(81, 0.29);
SetMSTP(95, 6);
SetPARJ(41, 0.3);
SetPARJ(42, 0.58);
SetPARJ(46, 0.75);
SetPARP(89,1800.);
}
void AliPythia6::SetWeightPower(Double_t pow)
{
setpowwght(pow);
SetMSTP(142, 1);
}
void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
{
SetCKIN(3, ptmin);
SetCKIN(4, ptmax);
}
void AliPythia6::SetYHardRange(Float_t ymin, Float_t ymax)
{
SetCKIN(7, ymin);
SetCKIN(8, ymax);
}
void AliPythia6::SetFragmentation(Int_t flag)
{
SetMSTP(111, flag);
}
void AliPythia6::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
{
SetMSTP(61, flag1);
SetMSTP(71, flag2);
}
void AliPythia6::SetIntrinsicKt(Float_t kt)
{
if (kt > 0.) {
SetMSTP(91,1);
SetPARP(91,kt);
SetPARP(93, 4. * kt);
} else {
SetMSTP(91,0);
}
}
void AliPythia6::SwitchHFOff()
{
SetMSTP(58, 3);
SetMSTJ(45, 3);
}
void AliPythia6::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
{
SetPARU(51, etamax);
SetMSTU(51, neta);
SetMSTU(52, nphi);
SetPARU(58, thresh);
SetPARU(52, etseed);
SetPARU(53, minet);
SetPARU(54, r);
SetMSTU(54, 2);
}
void AliPythia6::ModifiedSplitting()
{
SetPARJ(200, 0.8);
SetMSTJ(41, 1);
SetMSTJ(42, 2);
SetMSTJ(44, 2);
SetMSTJ(47, 0);
SetMSTJ(50, 0);
SetPARJ(82, 1.);
}
void AliPythia6::SwitchHadronisationOff()
{
SetMSTJ(1, 0);
}
void AliPythia6::SwitchHadronisationOn()
{
SetMSTJ(1, 1);
}
void AliPythia6::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
{
q = GetVINT(51);
x1 = GetVINT(41);
x2 = GetVINT(42);
}
Float_t AliPythia6::GetXSection()
{
return (GetPARI(1));
}
Float_t AliPythia6::GetPtHard()
{
return GetVINT(47);
}
Int_t AliPythia6::ProcessCode()
{
return GetMSTI(1);
}
void AliPythia6::PrintStatistics()
{
Pystat(1);
}
void AliPythia6::EventListing()
{
Pylist(2);
}
AliPythia6& AliPythia6::operator=(const AliPythia6& rhs)
{
rhs.Copy(*this);
return *this;
}
void AliPythia6::Copy(TObject&) const
{
Fatal("Copy","Not implemented!\n");
}