#include "AliPythia.h"
#include "AliPythiaRndm.h"
#include "AliRun.h"
+#include "AliStack.h"
+#include "AliRunLoader.h"
+#include "AliMC.h"
- ClassImp(AliGenPythia)
+ClassImp(AliGenPythia)
AliGenPythia::AliGenPythia()
:AliGenMC()
SetGammaPhiRange();
SetGammaEtaRange();
SetPtKick();
+ SetQuench();
+
fSetNuclei = kFALSE;
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
// Produced particles
fParticles = new TClonesArray("TParticle",1000);
fHeader = 0;
- fEventVertex.Set(3);
SetEventListRange();
SetJetPhiRange();
SetJetEtaRange();
SetGammaPhiRange();
SetGammaEtaRange();
SetJetReconstructionMode();
+ SetQuench();
SetPtKick();
// Options determining what to keep in the stack (Heavy flavour generation)
fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
}
AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
+ :AliGenMC(Pythia)
{
// copy constructor
Pythia.Copy(*this);
// Initialisation
SetMC(AliPythia::Instance());
- fPythia=(AliPythia*) fgMCEvGen;
+ fPythia=(AliPythia*) fMCEvGen;
+
//
fParentWeight=1./Float_t(fNpart);
//
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
- Float_t p[3];
+ Float_t p[4];
// converts from mm/c to s
const Float_t kconv=0.001/2.999792458e8;
//
Int_t jev=0;
Int_t j, kf;
fTrials=0;
+
// Set collision vertex position
- if(fVertexSmear==kPerEvent) {
- fPythia->SetMSTP(151,1);
- for (j=0;j<3;j++) {
- fPythia->SetPARP(151+j, fOsigma[j]*10.);
- }
- } else if (fVertexSmear==kPerTrack) {
- fPythia->SetMSTP(151,0);
- }
+ if (fVertexSmear == kPerEvent) Vertex();
+
// event loop
while(1)
{
+ if (fQuench) {
+ fPythia->SetMSTJ(1, 0);
+ }
fPythia->Pyevnt();
- if (gAlice->GetEvNumber()>=fDebugEventFirst &&
- gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
- fTrials++;
+
+ if (fQuench) {
+ fPythia->Quench();
+ fPythia->SetMSTJ(1, 1);
+ fPythia->Pyexec();
+ }
+ if (gAlice) {
+ if (gAlice->GetEvNumber()>=fDebugEventFirst &&
+ gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
+ }
+
+ fTrials++;
fPythia->ImportParticles(fParticles,"All");
Boost();
Int_t np = fParticles->GetEntriesFast();
if (np == 0 ) continue;
-// Get event vertex and discard the event if the Z coord. is too big
- TParticle *iparticle = (TParticle *) fParticles->At(0);
- Float_t distz = iparticle->Vz()/10.;
- if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
//
- fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
- fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
- fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
+
//
Int_t* pParent = new Int_t[np];
Int_t* pSelected = new Int_t[np];
fProcess != kPyMbNonDiffr) {
for (i = 0; i<np; i++) {
- iparticle = (TParticle *) fParticles->At(i);
+ TParticle* iparticle = (TParticle *) fParticles->At(i);
Int_t ks = iparticle->GetStatusCode();
kf = CheckPDGCode(iparticle->GetPdgCode());
// No initial state partons
// Track final state particle
if (ks == 1) trackIt[i] = 1;
// Track semi-stable particles
- if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
+ if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
// Track particles selected by process if undecayed.
if (fForceDecay == kNoDecay) {
if (ParentSelected(kf)) trackIt[i] = 1;
p[0] = iparticle->Px();
p[1] = iparticle->Py();
p[2] = iparticle->Pz();
- origin[0] = fOrigin[0]+iparticle->Vx()/10.;
- origin[1] = fOrigin[1]+iparticle->Vy()/10.;
- origin[2] = fOrigin[2]+iparticle->Vz()/10.;
+ p[3] = iparticle->Energy();
+
+ origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
+ origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
+ origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
+
Float_t tof = kconv*iparticle->T();
Int_t ipa = iparticle->GetFirstMother()-1;
Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
- SetTrack(fTrackIt*trackIt[i] ,
- iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
+
+ PushTrack(fTrackIt*trackIt[i], iparent, kf,
+ p[0], p[1], p[2], p[3],
+ origin[0], origin[1], origin[2], tof,
+ polar[0], polar[1], polar[2],
+ kPPrimary, nt, 1., ks);
pParent[i] = nt;
KeepTrack(nt);
- } // SetTrack loop
+ } // PushTrack loop
}
} else {
nc = GenerateMB();
} // event loop
SetHighWaterMark(nt);
// adjust weight due to kinematic selection
- AdjustWeights();
+// AdjustWeights();
// get cross-section
fXsection=fPythia->GetPARI(1);
}
//
Int_t i, kf, nt, iparent;
Int_t nc = 0;
- Float_t p[3];
+ Float_t p[4];
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
// converts from mm/c to s
p[0] = iparticle->Px();
p[1] = iparticle->Py();
p[2] = iparticle->Pz();
- origin[0] = fOrigin[0]+iparticle->Vx()/10.;
- origin[1] = fOrigin[1]+iparticle->Vy()/10.;
- origin[2] = fOrigin[2]+iparticle->Vz()/10.;
+ p[3] = iparticle->Energy();
+
+ origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
+ origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
+ origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
+
Float_t tof=kconv*iparticle->T();
- SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
- tof, kPPrimary, nt, 1., ks);
+
+ PushTrack(fTrackIt*trackIt, iparent, kf,
+ p[0], p[1], p[2], p[3],
+ origin[0], origin[1], origin[2], tof,
+ polar[0], polar[1], polar[2],
+ kPPrimary, nt, 1., ks);
KeepTrack(nt);
pParent[i] = nt;
} // select particle
{
// Adjust the weights after generation of all events
//
- TParticle *part;
- Int_t ntrack=gAlice->GetNtrack();
- for (Int_t i=0; i<ntrack; i++) {
- part= gAlice->Particle(i);
- part->SetWeight(part->GetWeight()*fKineBias);
+ if (gAlice) {
+ TParticle *part;
+ Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
+ for (Int_t i=0; i<ntrack; i++) {
+ part= gAlice->GetMCApp()->Particle(i);
+ part->SetWeight(part->GetWeight()*fKineBias);
+ }
}
}
((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
//
// Event Vertex
- fHeader->SetPrimaryVertex(fEventVertex);
+ fHeader->SetPrimaryVertex(fVertex);
//
// Jets that have triggered
if (fProcess == kPyJets)
jets[3][i]);
}
}
- gAlice->SetGenEventHeader(fHeader);
+ if (gAlice) gAlice->SetGenEventHeader(fHeader);
}
AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
{
// Assignment operator
+ rhs.Copy(*this);
return *this;
}
// Load event into Pythia Common Block
//
-
- Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries();
+ AliRunLoader* rl = AliRunLoader::GetRunLoader();
+ Int_t npart = (rl->Stack())-> GetNprimary();
(fPythia->GetPyjets())->N = npart;
for (Int_t part = 0; part < npart; part++) {
- TParticle *MPart = gAlice->Particle(part);
+ TParticle *MPart = (rl->Stack())->Particle(part);
Int_t kf = MPart->GetPdgCode();
Int_t ks = MPart->GetStatusCode();
Float_t px = MPart->Px();
Float_t py = MPart->Py();
Float_t pz = MPart->Pz();
Float_t e = MPart->Energy();
- Float_t p = TMath::Sqrt(px * px + py * py + pz * pz);
- Float_t m = TMath::Sqrt(e * e - p * p);
+ Float_t m = MPart->GetCalcMass();
(fPythia->GetPyjets())->P[0][part] = px;
}
}
-void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax,
- Int_t& njets, Float_t jets [4][50])
+void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
{
//
// Calls the Pythia jet finding algorithm to find jets in the current event
Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
Float_t e = (fPythia->GetPyjets())->P[3][n+i];
Float_t pt = TMath::Sqrt(px * px + py * py);
- Float_t phi = TMath::ATan2(py,px);
+ Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
Float_t theta = TMath::ATan2(pt,pz);
Float_t et = e * TMath::Sin(theta);
Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
if (
eta > fEtaMinJet && eta < fEtaMaxJet &&
- phi > fPhiMinJet && eta < fPhiMaxJet &&
+ phi > fPhiMinJet && phi < fPhiMaxJet &&
et > fEtMinJet && et < fEtMaxJet
)
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff