/*
$Log$
+Revision 1.30 2001/06/15 07:55:04 morsch
+Put only first generation decay products on the stack.
+
Revision 1.29 2001/03/27 10:58:41 morsch
Initialize decayer before generation. Important if run inside cocktail.
fYPara = 0;
fParam = 0;
fAnalog = kAnalog;
- SetCutOnChild();
- SetChildMomentumRange();
- SetChildPtRange();
- SetChildPhiRange();
- SetChildThetaRange();
SetDeltaPt();
//
// Set random number generator
sRandom = fRandom;
}
-AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Int_t param, char* tname):AliGenerator(npart)
+AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Int_t param, char* tname):AliGenMC(npart)
{
// Constructor using number of particles parameterisation id and library
fYPara = 0;
fParam = param;
fAnalog = kAnalog;
- fChildSelect.Set(5);
- for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
SetForceDecay();
- SetCutOnChild();
- SetChildMomentumRange();
- SetChildPtRange();
- SetChildPhiRange();
- SetChildThetaRange();
SetDeltaPt();
//
// Set random number generator
//____________________________________________________________
-AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenerator(npart)
+AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenMC(npart)
{
// Constructor using parameterisation id and number of particles
//
Double_t (*PtPara) (Double_t*, Double_t*),
Double_t (*YPara ) (Double_t* ,Double_t*),
Int_t (*IpPara) (TRandom *))
- :AliGenerator(npart)
+ :AliGenMC(npart)
{
// Constructor
// Gines Martinez 1/10/99
fDecayer->Init();
//
- switch (fForceDecay)
- {
- case kSemiElectronic:
- case kDiElectron:
- case kBJpsiDiElectron:
- case kBPsiPrimeDiElectron:
- fChildSelect[0]=11;
- break;
- case kSemiMuonic:
- case kDiMuon:
- case kBJpsiDiMuon:
- case kBPsiPrimeDiMuon:
- fChildSelect[0]=13;
- break;
- case kPiToMu:
- fChildSelect[0]=13;
- break;
- case kKaToMu:
- fChildSelect[0]=13;
- break;
- case kHadronicD:
-// Implement me !!
- break;
- case kNoDecay:
- break;
- case kAll:
- break;
- }
+ AliGenMC::Init();
}
//____________________________________________________________
Float_t pt, pl, ptot; // Transverse, logitudinal and total momenta of the parent particle
Float_t phi, theta; // Phi and theta spherical angles of the parent particle momentum
Float_t p[3], pc[3],
- och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
+ och[3]; // Momentum, polarisation and origin of the children particles from lujet
Float_t ty, xmt;
- Int_t nt, i, j, kfch[10];
+ Int_t nt, i, j;
Float_t wgtp, wgtch;
Double_t dummy;
static TClonesArray *particles;
//
- if(!particles) particles=new TClonesArray("TParticle",1000);
+ if(!particles) particles = new TClonesArray("TParticle",1000);
static TDatabasePDG *pDataBase = new TDatabasePDG();
if(!pDataBase) pDataBase = new TDatabasePDG();
// Calculating vertex position per event
for (j=0;j<3;j++) origin0[j]=fOrigin[j];
if(fVertexSmear==kPerEvent) {
-// Rndm(random,6);
-// for (j=0;j<3;j++) {
-// origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
-// TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-// TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-// }
-// }
Vertex();
for (j=0;j<3;j++) origin0[j]=fVertex[j];
}
Int_t ipa=0;
+
// Generating fNpart particles
while (ipa<fNpart) {
while(1) {
ptot=TMath::Sqrt(pt*pt+pl*pl);
// Cut on momentum
if(ptot<fPMin || ptot>fPMax) continue;
+//
p[0]=pt*TMath::Cos(phi);
p[1]=pt*TMath::Sin(phi);
p[2]=pl;
// if fForceDecay == none Primary particle is tracked by GEANT
// (In the latest, make sure that GEANT actually does all the decays you want)
//
+
if (fForceDecay != kNoDecay) {
// Using lujet to decay particle
Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
// select decay particles
Int_t np=fDecayer->ImportParticles(particles);
Int_t ncsel=0;
-
+ Int_t* pFlag = new Int_t[np];
+ Int_t* pParent = new Int_t[np];
+ Int_t* pSelected = new Int_t[np];
+ Int_t* trackIt = new Int_t[np];
+
+ for (i=0; i<np; i++) {
+ pFlag[i] = 0;
+ pSelected[i] = 0;
+ pParent[i] = -1;
+ }
+
if (np >1) {
TParticle* iparticle = (TParticle *) particles->At(0);
- Int_t ipF = iparticle->GetFirstDaughter();
- Int_t ipL = iparticle->GetLastDaughter();
- for (i = ipF-1; i<ipL; i++) {
+ Int_t ipF, ipL;
+ for (i = 1; i<np ; i++) {
+ trackIt[i] = 1;
iparticle = (TParticle *) particles->At(i);
Int_t kf = iparticle->GetPdgCode();
+ Int_t ks = iparticle->GetStatusCode();
+// flagged particle
+
+ if (pFlag[i] == 1) {
+ printf("\n deselected %d", i);
+ ipF = iparticle->GetFirstDaughter();
+ ipL = iparticle->GetLastDaughter();
+ if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
+ continue;
+ }
+
+// flag decay products of particles with long life-time (c tau > .3 mum)
+
+ if (ks != 1) {
+ TParticlePDG *particle = pDataBase->GetParticle(kf);
+
+ Double_t lifeTime = fDecayer->GetLifetime(kf);
+ Double_t mass = particle->Mass();
+ Double_t width = particle->Width();
+ printf("\n lifetime %d %e %e %e ", i, lifeTime, mass, width);
+ if (lifeTime > (Double_t) 1.e-15) {
+ ipF = iparticle->GetFirstDaughter();
+ ipL = iparticle->GetLastDaughter();
+ if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
+ } else{
+ trackIt[i] = 0;
+ pSelected[i] = 1;
+ }
+ } // ks==1 ?
//
// children
- if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll)
+ if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
{
- pc[0]=iparticle->Px();
- pc[1]=iparticle->Py();
- pc[2]=iparticle->Pz();
- och[0]=origin0[0]+iparticle->Vx()/10;
- och[1]=origin0[1]+iparticle->Vy()/10;
- och[2]=origin0[2]+iparticle->Vz()/10;
if (fCutOnChild) {
- Float_t ptChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
- Float_t pChild=TMath::Sqrt(ptChild*ptChild+pc[2]*pc[2]);
- Float_t thetaChild=TMath::ATan2(ptChild,pc[2]);
- Float_t phiChild=TMath::ATan2(pc[1],pc[0]);
- Bool_t childok =
- ((ptChild > fChildPtMin && ptChild <fChildPtMax) &&
- (pChild > fChildPMin && pChild <fChildPMax) &&
- (thetaChild > fChildThetaMin && thetaChild <fChildThetaMax) &&
- (phiChild > fChildPhiMin && phiChild <fChildPhiMax));
- if(childok)
- {
- pch[ncsel][0]=pc[0];
- pch[ncsel][1]=pc[1];
- pch[ncsel][2]=pc[2];
- kfch[ncsel]=kf;
+ pc[0]=iparticle->Px();
+ pc[1]=iparticle->Py();
+ pc[2]=iparticle->Pz();
+ Bool_t childok = KinematicSelection(iparticle, 1);
+ if(childok) {
+ pSelected[i] = 1;
ncsel++;
} else {
ncsel=-1;
break;
} // child kine cuts
} else {
- pch[ncsel][0]=pc[0];
- pch[ncsel][1]=pc[1];
- pch[ncsel][2]=pc[2];
- kfch[ncsel]=kf;
+ pSelected[i] = 1;
ncsel++;
} // if child selection
} // select muon
if ((fCutOnChild && ncsel >0) || !fCutOnChild){
ipa++;
//
-// parent
- gAlice->
- SetTrack(0,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
- iparent=nt;
+// Parent
+ gAlice->SetTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
+ pParent[0] = nt;
gAlice->KeepTrack(nt);
- for (i=0; i< ncsel; i++) {
- gAlice->SetTrack(fTrackIt,iparent,kfch[i],
- &pch[i][0],och,polar,
- 0,kPDecay,nt,wgtch);
- gAlice->KeepTrack(nt);
+//
+// Decay Products
+//
+ for (i = 1; i < np; i++) {
+ if (pSelected[i]) {
+ TParticle* iparticle = (TParticle *) particles->At(i);
+ Int_t kf = iparticle->GetPdgCode();
+ Int_t ipa = iparticle->GetFirstMother()-1;
+
+ och[0] = origin0[0]+iparticle->Vx()/10;
+ och[1] = origin0[1]+iparticle->Vy()/10;
+ och[2] = origin0[2]+iparticle->Vz()/10;
+ pc[0] = iparticle->Px();
+ pc[1] = iparticle->Py();
+ pc[2] = iparticle->Pz();
+
+ if (ipa > -1) {
+ iparent = pParent[ipa];
+ } else {
+ iparent = -1;
+ }
+ printf("\n SetTrack %d %d %d %d", i, kf, iparent, trackIt[i]);
+ gAlice->SetTrack(fTrackIt*trackIt[i], iparent, kf,
+ pc, och, polar,
+ 0, kPDecay, nt, wgtch);
+ pParent[i] = nt;
+ gAlice->KeepTrack(nt);
+ }
}
} // Decays by Lujet
+
+ if (pFlag) delete[] pFlag;
+ if (pParent) delete[] pParent;
+ if (pSelected) delete[] pSelected;
+ if (trackIt) delete[] trackIt;
} // kinematic selection
else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
{
} // while
} // event loop
gAlice->SetHighWaterMark(nt);
-
}
-Bool_t AliGenParam::ChildSelected(Int_t ip)
-{
-// True if particle is in list of selected children
- for (Int_t i=0; i<5; i++)
- {
- if (fChildSelect[i]==ip) return kTRUE;
- }
- return kFALSE;
-}
-
-Bool_t AliGenParam::KinematicSelection(TParticle *particle)
-{
-// Perform kinematic cuts
- Float_t px=particle->Px();
- Float_t py=particle->Py();
- Float_t pz=particle->Pz();
-//
-// momentum cut
- Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
- if (p > fPMax || p < fPMin)
- {
-// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
- return kFALSE;
- }
- Float_t pt=TMath::Sqrt(px*px+py*py);
-
-//
-// theta cut
- Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
- if (theta > fThetaMax || theta < fThetaMin)
- {
-// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
- return kFALSE;
- }
-
- return kTRUE;
-}
-
-
AliGenParam& AliGenParam::operator=(const AliGenParam& rhs)
{
// Assignment operator
/*
$Log$
+Revision 1.37 2001/06/28 11:17:28 morsch
+SetEventListRange setter added. Events in specified range are listed for
+debugging. (Yuri Kharlov)
+
Revision 1.36 2001/03/30 07:05:49 morsch
Final print-out in finish run.
Write parton system for jet-production (preliminary solution).
*/
#include "AliGenPythia.h"
+#include "AliGenPythiaEventHeader.h"
#include "AliDecayerPythia.h"
#include "AliRun.h"
#include "AliPythia.h"
ClassImp(AliGenPythia)
AliGenPythia::AliGenPythia()
- :AliGenerator()
+ :AliGenMC()
{
// Default Constructor
fDecayer = new AliDecayerPythia();
}
AliGenPythia::AliGenPythia(Int_t npart)
- :AliGenerator(npart)
+ :AliGenMC(npart)
{
// default charm production at 5. 5 TeV
// semimuonic decay
fXsection = 0.;
fNucA1=0;
fNucA2=0;
- fParentSelect.Set(5);
- fChildSelect.Set(5);
- for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
SetProcess();
SetStrucFunc();
SetForceDecay();
// Set random number generator
sRandom=fRandom;
SetEventListRange();
+ fFlavorSelect = 0;
+ // Produced particles
+ fParticles = new TClonesArray("TParticle",1000);
}
AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
fParentWeight=1./Float_t(fNpart);
//
// Forward Paramters to the AliPythia object
- // gSystem->Exec("ln -s $ALICE_ROOT/data/Decay.table fort.1");
- // fPythia->Pyupda(2,1);
- // gSystem->Exec("rm fort.1");
-
fDecayer->SetForceDecay(fForceDecay);
fDecayer->Init();
switch (fProcess)
{
case kPyCharm:
- fParentSelect[0]=411;
- fParentSelect[1]=421;
- fParentSelect[2]=431;
- fParentSelect[3]=4122;
+ fParentSelect[0] = 411;
+ fParentSelect[1] = 421;
+ fParentSelect[2] = 431;
+ fParentSelect[3] = 4122;
+ fFlavorSelect = 4;
break;
case kPyCharmUnforced:
- fParentSelect[0]=411;
- fParentSelect[1]=421;
- fParentSelect[2]=431;
- fParentSelect[3]=4122;
+ fParentSelect[0] = 411;
+ fParentSelect[1] = 421;
+ fParentSelect[2] = 431;
+ fParentSelect[3]= 4122;
+ fFlavorSelect = 4;
break;
case kPyBeauty:
- fParentSelect[0]=511;
- fParentSelect[1]=521;
- fParentSelect[2]=531;
- fParentSelect[3]=5122;
+ fParentSelect[0]= 511;
+ fParentSelect[1]= 521;
+ fParentSelect[2]= 531;
+ fParentSelect[3]= 5122;
+ fParentSelect[4]= 5132;
+ fParentSelect[5]= 5232;
+ fParentSelect[6]= 5332;
+ fFlavorSelect = 5;
break;
case kPyBeautyUnforced:
- fParentSelect[0]=511;
- fParentSelect[1]=521;
- fParentSelect[2]=531;
- fParentSelect[3]=5122;
+ fParentSelect[0] = 511;
+ fParentSelect[1] = 521;
+ fParentSelect[2] = 531;
+ fParentSelect[3] = 5122;
+ fParentSelect[4] = 5132;
+ fParentSelect[5] = 5232;
+ fParentSelect[6] = 5332;
+ fFlavorSelect = 5;
break;
case kPyJpsiChi:
case kPyJpsi:
- fParentSelect[0]=443;
+ fParentSelect[0] = 443;
break;
case kPyMb:
case kPyJets:
case kPyDirectGamma:
break;
}
-
- switch (fForceDecay)
- {
- case kSemiElectronic:
- case kDiElectron:
- case kBJpsiDiElectron:
- case kBPsiPrimeDiElectron:
- fChildSelect[0]=kElectron;
- break;
- case kSemiMuonic:
- case kDiMuon:
- case kBJpsiDiMuon:
- case kBPsiPrimeDiMuon:
- case kPiToMu:
- case kKaToMu:
- fChildSelect[0]=kMuonMinus;
- break;
- case kHadronicD:
- fChildSelect[0]=kPiPlus;
- fChildSelect[1]=kKPlus;
- break;
- case kAll:
- case kNoDecay:
- break;
- }
+ AliGenMC::Init();
}
void AliGenPythia::Generate()
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
- Float_t originP[3] = {0,0,0};
Float_t origin0[3] = {0,0,0};
- Float_t p[3], pP[4];
-// Float_t random[6];
- static TClonesArray *particles;
+ Float_t p[3];
// converts from mm/c to s
const Float_t kconv=0.001/2.999792458e8;
-
-
//
Int_t nt=0;
- Int_t ntP=0;
Int_t jev=0;
Int_t j, kf;
-
- if(!particles) particles=new TClonesArray("TParticle",1000);
-
fTrials=0;
+
+// Set collision vertex position
for (j=0;j<3;j++) origin0[j]=fOrigin[j];
if(fVertexSmear==kPerEvent) {
fPythia->SetMSTP(151,1);
for (j=0;j<3;j++) {
- fPythia->SetPARP(151+j, fOsigma[j]/10.);
+ fPythia->SetPARP(151+j, fOsigma[j]*10.);
}
} else if (fVertexSmear==kPerTrack) {
fPythia->SetMSTP(151,0);
}
-
+// event loop
while(1)
{
fPythia->Pyevnt();
if (gAlice->GetEvNumber()>=fDebugEventFirst &&
gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
fTrials++;
- fPythia->ImportParticles(particles,"All");
- Int_t np = particles->GetEntriesFast();
+ fPythia->ImportParticles(fParticles,"All");
+
+//
+//
+//
+ Int_t i;
+
+ Int_t np = fParticles->GetEntriesFast();
+ Int_t* pParent = new Int_t[np];
+ Int_t* pSelected = new Int_t[np];
+ Int_t* trackIt = new Int_t[np];
+ for (i=0; i< np-1; i++) {
+ pParent[i] = -1;
+ pSelected[i] = 0;
+ }
printf("\n **************************************************%d\n",np);
- Int_t nc=0;
+ Int_t nc = 0;
if (np == 0 ) continue;
if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
- for (Int_t i = 0; i<np-1; i++) {
- TParticle * iparticle = (TParticle *) particles->At(i);
+ for (i = 0; i<np-1; i++) {
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
Int_t ks = iparticle->GetStatusCode();
kf = CheckPDGCode(iparticle->GetPdgCode());
+// No initial state partons
if (ks==21) continue;
- fChildWeight=(fDecayer->GetPartialBranchingRatio(kf))*fParentWeight;
//
-// Parent
- if (ParentSelected(TMath::Abs(kf))) {
- if (KinematicSelection(iparticle)) {
- if (nc==0) {
+// Heavy Flavor Selection
//
-// Store information concerning the hard scattering process
+ // quark ?
+ kf = TMath::Abs(kf);
+ Int_t kfl = kf;
+ // meson ?
+ if (kfl > 10) kfl/=100;
+ // baryon
+ if (kfl > 10) kfl/=10;
+ if (kfl > 10) kfl/=10;
+
+ Int_t ipa = iparticle->GetFirstMother()-1;
+ Int_t kfMo = 0;
+
+ if (ipa > -1) {
+ TParticle * mother = (TParticle *) fParticles->At(ipa);
+ kfMo = TMath::Abs(mother->GetPdgCode());
+ }
+// printf("\n particle (all) %d %d %d", i, pSelected[i], kf);
+ if (kfl >= fFlavorSelect) {
//
- Float_t massP = fPythia->GetPARI(13);
- Float_t ptP = fPythia->GetPARI(17);
- Float_t yP = fPythia->GetPARI(37);
- Float_t xmtP = sqrt(ptP*ptP+massP*massP);
- Float_t ty = Float_t(TMath::TanH(yP));
- pP[0] = ptP;
- pP[1] = 0;
- pP[2] = xmtP*ty/sqrt(1.-ty*ty);
- pP[3] = massP;
- gAlice->SetTrack(0,-1,-1,
- pP,originP,polar,
- 0,kPPrimary,ntP,fParentWeight);
-// 0,"Hard Scat.",ntP,fParentWeight);
- gAlice->KeepTrack(ntP);
- }
- nc++;
+// Heavy flavor hadron or quark
//
-// store parent track information
- p[0]=iparticle->Px();
- p[1]=iparticle->Py();
- p[2]=iparticle->Pz();
- origin[0]=origin0[0]+iparticle->Vx()/10.;
- origin[1]=origin0[1]+iparticle->Vy()/10.;
- origin[2]=origin0[2]+iparticle->Vz()/10.;
-
- Int_t ifch=iparticle->GetFirstDaughter();
- Int_t ilch=iparticle->GetLastDaughter();
-
- if ((ifch !=0 && ilch !=0) || fForceDecay == kNoDecay) {
- Int_t trackit=0;
- if (fForceDecay == kNoDecay) trackit = 1;
- gAlice->SetTrack(trackit,ntP,kf,
- p,origin,polar,
- 0,kPPrimary,nt,fParentWeight);
- gAlice->KeepTrack(nt);
- Int_t iparent = nt;
+// Kinematic seletion on final state heavy flavor mesons
+ if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
+ {
+ nc = -1;
+ break;
+ }
+ pSelected[i] = 1;
+// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
+ } else {
+// Kinematic seletion on decay products
+ if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
+ && !KinematicSelection(iparticle, 1))
+ {
+ nc = -1;
+ break;
+ }
//
-// Children
- if (fForceDecay != kNoDecay) {
- for (j=ifch; j<=ilch; j++)
- {
- TParticle * ichild =
- (TParticle *) particles->At(j-1);
- kf = CheckPDGCode(ichild->GetPdgCode());
+// Decay products
+// Select if mother was selected and is not tracked
+
+ if (pSelected[ipa] &&
+ !trackIt[ipa] && // mother will be tracked ?
+ kfMo != 5 && // mother is b-quark, don't store fragments
+ kfMo != 4 && // mother is c-quark, don't store fragments
+ kf != 92) // don't store string
+ {
//
+// Semi-stable or de-selected: diselect decay products:
//
- if (ChildSelected(TMath::Abs(kf))) {
- origin[0]=origin0[0]+ichild->Vx()/10.;
- origin[1]=origin0[1]+ichild->Vy()/10.;
- origin[2]=origin0[2]+ichild->Vz()/10.;
- p[0]=ichild->Px();
- p[1]=ichild->Py();
- p[2]=ichild->Pz();
- Float_t tof=kconv*ichild->T();
- gAlice->SetTrack(fTrackIt, iparent, kf,
- p,origin,polar,
- tof,kPDecay,nt,fChildWeight);
- gAlice->KeepTrack(nt);
- } // select child
- } // child loop
- }
- }
- } // kinematic selection
- } // select particle
- } // particle loop
- } else {
- for (Int_t i = 0; i<np-1; i++) {
- TParticle * iparticle = (TParticle *) particles->At(i);
- kf = CheckPDGCode(iparticle->GetPdgCode());
- Int_t ks = iparticle->GetStatusCode();
- Int_t km = iparticle->GetFirstMother();
- // printf("\n process %d %d\n", ks,km);
-
- if ((ks==1 && kf!=0 && KinematicSelection(iparticle)) ||
- (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
- nc++;
//
-// store track information
+ if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > 10e-15)
+ {
+ Int_t ipF = iparticle->GetFirstDaughter();
+ Int_t ipL = iparticle->GetLastDaughter();
+ if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
+ }
+// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
+ pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
+ }
+ }
+ if (pSelected[i] == -1) pSelected[i] = 0;
+ if (!pSelected[i]) continue;
+ nc++;
+// Decision on tracking
+ trackIt[i] = 0;
+//
+// Track final state particle
+ if (ks == 1) trackIt[i] = 1;
+// Track semi-stable particles
+ if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15)) trackIt[i] = 1;
+// Track particles selected by process if undecayed.
+ if (fForceDecay == kNoDecay) {
+ if (ParentSelected(kf)) trackIt[i] = 1;
+ } else {
+ if (ParentSelected(kf)) trackIt[i] = 0;
+ }
+//
+//
+
+ } // particle selection loop
+ if (nc > -1) {
+ for (i = 0; i<np-1; i++) {
+ if (!pSelected[i]) continue;
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
+ kf = CheckPDGCode(iparticle->GetPdgCode());
p[0]=iparticle->Px();
p[1]=iparticle->Py();
p[2]=iparticle->Pz();
origin[1]=origin0[1]+iparticle->Vy()/10.;
origin[2]=origin0[2]+iparticle->Vz()/10.;
Float_t tof=kconv*iparticle->T();
- gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
- tof,kPPrimary,nt);
- gAlice->KeepTrack(nt);
- } // select particle
- } // particle loop
- printf("\n I've put %i particles on the stack \n",nc);
+ Int_t ipa = iparticle->GetFirstMother()-1;
+ Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
+// printf("\n SetTrack %d %d %d %d", i, iparent, kf, trackIt[i]);
+ gAlice->SetTrack(fTrackIt*trackIt[i] ,
+ iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
+ pParent[i] = nt;
+ gAlice->KeepTrack(nt);
+ } // SetTrack loop
+ }
+ } else {
+ nc = GenerateMB();
} // mb ?
+
if (nc > 0) {
jev+=nc;
if (jev >= fNpart || fNpart == -1) {
fXsection=fPythia->GetPARI(1);
}
-void AliGenPythia::FinishRun()
-{
-// Print x-section summary
- fPythia->Pystat(1);
-}
-
-Bool_t AliGenPythia::ParentSelected(Int_t ip)
-{
-// True if particle is in list of parent particles to be selected
- for (Int_t i=0; i<5; i++)
- {
- if (fParentSelect[i]==ip) return kTRUE;
- }
- return kFALSE;
-}
-
-Bool_t AliGenPythia::ChildSelected(Int_t ip)
+Int_t AliGenPythia::GenerateMB()
{
-// True if particle is in list of decay products to be selected
- if (fForceDecay == kAll) return kTRUE;
+ Int_t i, kf, nt, iparent;
+ Int_t nc = 0;
+ Float_t p[3];
+ Float_t polar[3] = {0,0,0};
+ Float_t origin[3] = {0,0,0};
+ Float_t origin0[3] = {0,0,0};
+// converts from mm/c to s
+ const Float_t kconv=0.001/2.999792458e8;
- for (Int_t i=0; i<5; i++)
- {
- if (fChildSelect[i]==ip) return kTRUE;
- }
- return kFALSE;
-}
+ Int_t np = fParticles->GetEntriesFast();
+ Int_t* pParent = new Int_t[np];
+ for (i=0; i< np-1; i++) pParent[i] = -1;
+
+ for (i = 0; i<np-1; i++) {
+ Int_t trackIt = 0;
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
+ kf = CheckPDGCode(iparticle->GetPdgCode());
+ Int_t ks = iparticle->GetStatusCode();
+ Int_t km = iparticle->GetFirstMother();
+// printf("\n Particle: %d %d %d", i, kf, ks);
+
+ if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
+ (ks != 1) ||
+ (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+ nc++;
+ if (ks == 1) trackIt = 1;
+ Int_t ipa = iparticle->GetFirstMother()-1;
-Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
-{
-// Perform kinematic selection
- Float_t px=particle->Px();
- Float_t py=particle->Py();
- Float_t pz=particle->Pz();
- Float_t e=particle->Energy();
+ iparent = (ipa > -1) ? pParent[ipa] : -1;
//
-// transverse momentum cut
- Float_t pt=TMath::Sqrt(px*px+py*py);
- if (pt > fPtMax || pt < fPtMin)
- {
-// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
- return kFALSE;
- }
-//
-// momentum cut
- Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
- if (p > fPMax || p < fPMin)
- {
-// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
- return kFALSE;
- }
+// store track information
+ p[0]=iparticle->Px();
+ p[1]=iparticle->Py();
+ p[2]=iparticle->Pz();
+ origin[0]=origin0[0]+iparticle->Vx()/10.;
+ origin[1]=origin0[1]+iparticle->Vy()/10.;
+ origin[2]=origin0[2]+iparticle->Vz()/10.;
+ Float_t tof=kconv*iparticle->T();
+ gAlice->SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
+ tof, kPPrimary, nt);
+ gAlice->KeepTrack(nt);
+ pParent[i] = nt;
+ } // select particle
+ } // particle loop
+
+ if (pParent) delete[] pParent;
-//
-// theta cut
- Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
- if (theta > fThetaMax || theta < fThetaMin)
- {
-// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
- return kFALSE;
- }
-
-//
-// rapidity cut
- if ( (e-pz)<=0 || (e+pz)<=0 ) {
- return kFALSE;
- }
- else {
- Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
- if (y > fYMax || y < fYMin)
- {
-// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
- return kFALSE;
- }
- }
+ printf("\n I've put %i particles on the stack \n",nc);
+ MakeHeader();
+ return nc;
+}
-//
-// phi cut
- Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
- if (phi > fPhiMax || phi < fPhiMin)
- {
-// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
- return kFALSE;
- }
- return kTRUE;
+void AliGenPythia::FinishRun()
+{
+// Print x-section summary
+ fPythia->Pystat(1);
}
+
void AliGenPythia::AdjustWeights()
{
// Adjust the weights after generation of all events
part->SetWeight(part->GetWeight()*fKineBias);
}
}
-
-Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
-{
-//
-// If the particle is in a diffractive state, then take action accordingly
- switch (pdgcode) {
- case 91:
- return 92;
- case 110:
- //rho_diff0 -- difficult to translate, return rho0
- return 113;
- case 210:
- //pi_diffr+ -- change to pi+
- return 211;
- case 220:
- //omega_di0 -- change to omega0
- return 223;
- case 330:
- //phi_diff0 -- return phi0
- return 333;
- case 440:
- //J/psi_di0 -- return J/psi
- return 443;
- case 2110:
- //n_diffr -- return neutron
- return 2112;
- case 2210:
- //p_diffr+ -- return proton
- return 2212;
- }
- //non diffractive state -- return code unchanged
- return pdgcode;
-}
-
void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
{
fNucA1 = a1;
fNucA2 = a2;
}
+
+
+void AliGenPythia::MakeHeader()
+{
+// Builds the event header, to be called after each event
+ AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
+ ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+ gAlice->SetGenEventHeader(header);
+}
AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff