/*
$Log$
+Revision 1.60 2002/07/19 14:49:03 morsch
+Typo corrected.
+
+Revision 1.59 2002/07/19 14:35:36 morsch
+Count total number of trials. Print mean Q, x1, x2.
+
+Revision 1.58 2002/07/17 10:04:09 morsch
+SetYHard method added.
+
+Revision 1.57 2002/05/22 13:22:53 morsch
+Process kPyMbNonDiffr added.
+
+Revision 1.56 2002/04/26 10:30:01 morsch
+Option kPyBeautyPbMNR added. (N. Carrer).
+
+Revision 1.55 2002/04/17 10:23:56 morsch
+Coding Rule violations corrected.
+
+Revision 1.54 2002/03/28 11:49:10 morsch
+Pass status code in SetTrack.
+
+Revision 1.53 2002/03/25 14:51:13 morsch
+New stack-fill and count options introduced (N. Carrer).
+
+Revision 1.51 2002/03/06 08:46:57 morsch
+- Loop until np-1
+- delete dyn. alloc. arrays (N. Carrer)
+
+Revision 1.50 2002/03/03 13:48:50 morsch
+Option kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
+NLO calculations (Nicola Carrer).
+
+Revision 1.49 2002/02/08 16:50:50 morsch
+Add name and title in constructor.
+
+Revision 1.48 2001/12/20 11:44:28 morsch
+Add kinematic bias for direct gamma production.
+
+Revision 1.47 2001/12/19 14:45:00 morsch
+Store number of trials in header.
+
+Revision 1.46 2001/12/19 10:36:19 morsch
+Add possibility if jet kinematic biasing.
+
+Revision 1.45 2001/11/28 08:06:52 morsch
+Use fMaxLifeTime parameter.
+
+Revision 1.44 2001/11/27 13:13:07 morsch
+Maximum lifetime for long-lived particles to be put on the stack is parameter.
+It can be set via SetMaximumLifetime(..).
+
+Revision 1.43 2001/10/21 18:35:56 hristov
+Several pointers were set to zero in the default constructors to avoid memory management problems
+
+Revision 1.42 2001/10/15 08:21:55 morsch
+Vertex truncation settings moved to AliGenMC.
+
+Revision 1.41 2001/10/08 08:45:42 morsch
+Possibility of vertex cut added.
+
+Revision 1.40 2001/09/25 11:30:23 morsch
+Pass event vertex to header.
+
+Revision 1.39 2001/07/27 17:09:36 morsch
+Use local SetTrack, KeepTrack and SetHighWaterMark methods
+to delegate either to local stack or to stack owned by AliRun.
+(Piotr Skowronski, A.M.)
+
+Revision 1.38 2001/07/13 10:58:54 morsch
+- Some coded moved to AliGenMC
+- Improved handling of secondary vertices.
+
+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).
Introduction of the Copyright and cvs Log
*/
+//
+// Generator using the TPythia interface (via AliPythia)
+// to generate pp collisions.
+// Using SetNuclei() also nuclear modifications to the structure functions
+// can be taken into account. This makes, of course, only sense for the
+// generation of the products of hard processes (heavy flavor, jets ...)
+//
+// andreas.morsch@cern.ch
+//
+
+
#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();
+ fParticles = 0;
+ fPythia = 0;
+ fDecayer = new AliDecayerPythia();
SetEventListRange();
+ SetJetPhiRange();
+ SetJetEtaRange();
+ SetGammaPhiRange();
+ SetGammaEtaRange();
}
AliGenPythia::AliGenPythia(Int_t npart)
- :AliGenerator(npart)
+ :AliGenMC(npart)
{
// default charm production at 5. 5 TeV
// semimuonic decay
// structure function GRVHO
//
+ fName = "Pythia";
+ fTitle= "Particle Generator using PYTHIA";
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();
SetPtHard();
+ SetYHard();
SetEnergyCMS();
fDecayer = new AliDecayerPythia();
// Set random number generator
sRandom=fRandom;
+ fFlavorSelect = 0;
+ // Produced particles
+ fParticles = new TClonesArray("TParticle",1000);
+ fEventVertex.Set(3);
SetEventListRange();
+ SetJetPhiRange();
+ SetJetEtaRange();
+ SetGammaPhiRange();
+ SetGammaEtaRange();
+ // Options determining what to keep in the stack (Heavy flavour generation)
+ fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
+ fFeedDownOpt = kTRUE; // allow feed down from higher family
+ // Fragmentation on/off
+ fFragmentation = kTRUE;
+ // Default counting mode
+ fCountMode = kCountAll;
}
AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
{
// copy constructor
+ Pythia.Copy(*this);
}
AliGenPythia::~AliGenPythia()
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();
fPythia->SetCKIN(3,fPtHardMin);
- fPythia->SetCKIN(4,fPtHardMax);
+ fPythia->SetCKIN(4,fPtHardMax);
+ fPythia->SetCKIN(7,fYHardMin);
+ fPythia->SetCKIN(8,fYHardMax);
+
if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
+ // Fragmentation?
+ if (fFragmentation) {
+ fPythia->SetMSTP(111,1);
+ } else {
+ fPythia->SetMSTP(111,0);
+ }
+
fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
// fPythia->Pylist(0);
switch (fProcess)
{
case kPyCharm:
- fParentSelect[0]=411;
- fParentSelect[1]=421;
- fParentSelect[2]=431;
- fParentSelect[3]=4122;
- break;
case kPyCharmUnforced:
- fParentSelect[0]=411;
- fParentSelect[1]=421;
- fParentSelect[2]=431;
- fParentSelect[3]=4122;
+ case kPyCharmPbMNR:
+ fParentSelect[0] = 411;
+ fParentSelect[1] = 421;
+ fParentSelect[2] = 431;
+ fParentSelect[3] = 4122;
+ fFlavorSelect = 4;
+ break;
+ case kPyD0PbMNR:
+ fParentSelect[0] = 421;
+ fFlavorSelect = 4;
break;
case kPyBeauty:
- fParentSelect[0]=511;
- fParentSelect[1]=521;
- fParentSelect[2]=531;
- fParentSelect[3]=5122;
+ case kPyBeautyPbMNR:
+ 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 kPyMbNonDiffr:
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;
- }
+//
+// This counts the total number of calls to Pyevnt() per run.
+ fTrialsRun = 0;
+ fQ = 0.;
+ fX1 = 0.;
+ fX2 = 0.;
+ fNev = 0 ;
+//
+ 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;
- for (j=0;j<3;j++) origin0[j]=fOrigin[j];
+
+ // Set collision vertex position
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();
- printf("\n **************************************************%d\n",np);
- Int_t nc=0;
+
+ fPythia->ImportParticles(fParticles,"All");
+
+//
+//
+//
+ Int_t i;
+
+ Int_t np = fParticles->GetEntriesFast();
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);
+// 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];
+ Int_t* trackIt = new Int_t[np];
+ for (i=0; i< np; i++) {
+ pParent[i] = -1;
+ pSelected[i] = 0;
+ trackIt[i] = 0;
+ }
+
+ Int_t nc = 0; // Total n. of selected particles
+ Int_t nParents = 0; // Selected parents
+ Int_t nTkbles = 0; // Trackable particles
+ if (fProcess != kPyMb && fProcess != kPyJets &&
+ fProcess != kPyDirectGamma &&
+ fProcess != kPyMbNonDiffr) {
+
+ for (i = 0; i<np; i++) {
+ 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());
+ }
+ // What to keep in Stack?
+ Bool_t flavorOK = kFALSE;
+ Bool_t selectOK = kFALSE;
+ if (fFeedDownOpt) {
+ if (kfl >= fFlavorSelect) flavorOK = kTRUE;
+ } else {
+ if (kfl > fFlavorSelect) {
+ nc = -1;
+ break;
+ }
+ if (kfl == fFlavorSelect) flavorOK = kTRUE;
+ }
+ switch (fStackFillOpt) {
+ case kFlavorSelection:
+ selectOK = kTRUE;
+ break;
+ case kParentSelection:
+ if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
+ break;
+ }
+ if (flavorOK && selectOK) {
//
- 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))
+ {
+ continue;
+ }
+ pSelected[i] = 1;
+ if (ParentSelected(kf)) ++nParents; // Update parent count
+// 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))
+ {
+ continue;
+ }
//
-// 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
- 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,-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);
+ if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
+ {
+ 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;
+ // Count quarks only if you did not include fragmentation
+ if (fFragmentation && kf <= 10) 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) > fMaxLifeTime)) 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;
+ }
+ if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
+//
+//
+
+ } // particle selection loop
+ if (nc > 0) {
+ for (i = 0; i<np; i++) {
+ if (!pSelected[i]) continue;
+ TParticle * iparticle = (TParticle *) fParticles->At(i);
+ kf = CheckPDGCode(iparticle->GetPdgCode());
+ Int_t ks = iparticle->GetStatusCode();
+ 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.;
+ 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);
+ pParent[i] = nt;
+ KeepTrack(nt);
+ } // SetTrack loop
+ }
+ } else {
+ nc = GenerateMB();
} // mb ?
+
+ if (pParent) delete[] pParent;
+ if (pSelected) delete[] pSelected;
+ if (trackIt) delete[] trackIt;
+
if (nc > 0) {
- jev+=nc;
+ switch (fCountMode) {
+ case kCountAll:
+ // printf(" Count all \n");
+ jev += nc;
+ break;
+ case kCountParents:
+ // printf(" Count parents \n");
+ jev += nParents;
+ break;
+ case kCountTrackables:
+ // printf(" Count trackable \n");
+ jev += nTkbles;
+ break;
+ }
if (jev >= fNpart || fNpart == -1) {
fKineBias=Float_t(fNpart)/Float_t(fTrials);
printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+
+ fQ += fPythia->GetVINT(51);
+ fX1 += fPythia->GetVINT(41);
+ fX2 += fPythia->GetVINT(42);
+ fTrialsRun += fTrials;
+ fNev++;
+ MakeHeader();
break;
}
}
} // event loop
- gAlice->SetHighWaterMark(nt);
+ SetHighWaterMark(nt);
// adjust weight due to kinematic selection
AdjustWeights();
// get cross-section
fXsection=fPythia->GetPARI(1);
}
-void AliGenPythia::FinishRun()
+Int_t AliGenPythia::GenerateMB()
{
-// 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)
-{
-// True if particle is in list of decay products to be selected
- if (fForceDecay == kAll) return kTRUE;
-
- for (Int_t i=0; i<5; i++)
- {
- if (fChildSelect[i]==ip) return kTRUE;
- }
- return kFALSE;
-}
-
-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();
-
//
-// 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;
- }
+// Min Bias selection and other global selections
//
-// 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;
+ 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};
+// converts from mm/c to s
+ const Float_t kconv=0.001/2.999792458e8;
+
+ Int_t np = fParticles->GetEntriesFast();
+ Int_t* pParent = new Int_t[np];
+ for (i=0; i< np; i++) pParent[i] = -1;
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+ TParticle* jet1 = (TParticle *) fParticles->At(6);
+ TParticle* jet2 = (TParticle *) fParticles->At(7);
+ if (!CheckTrigger(jet1, jet2)) return 0;
}
+ for (i = 0; i<np; 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();
+ 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;
+
+ iparent = (ipa > -1) ? pParent[ipa] : -1;
+
//
-// 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;
- }
+// store track information
+ 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.;
+ Float_t tof=kconv*iparticle->T();
+ SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
+ tof, kPPrimary, nt, 1., ks);
+ KeepTrack(nt);
+ pParent[i] = nt;
+ } // select particle
+ } // particle loop
+
+ if (pParent) delete[] pParent;
+
+ printf("\n I've put %i particles on the stack \n",nc);
+ return nc;
+}
-//
-// 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;
- }
- }
-//
-// 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;
- }
+void AliGenPythia::FinishRun()
+{
+// Print x-section summary
+ fPythia->Pystat(1);
+ fQ /= fNev;
+ fX1 /= fNev;
+ fX2 /= fNev;
+ printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
+ printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
+
- return kTRUE;
}
+
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() const
+{
+// Builds the event header, to be called after each event
+ AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
+ ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+ ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
+ header->SetPrimaryVertex(fEventVertex);
+ gAlice->SetGenEventHeader(header);
+}
+
+Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) const
+{
+// Check the kinematic trigger condition
+//
+ Double_t eta[2];
+ eta[0] = jet1->Eta();
+ eta[1] = jet2->Eta();
+ Double_t phi[2];
+ phi[0] = jet1->Phi();
+ phi[1] = jet2->Phi();
+ Int_t pdg[2];
+ pdg[0] = jet1->GetPdgCode();
+ pdg[1] = jet2->GetPdgCode();
+ Bool_t triggered = kFALSE;
+
+ if (fProcess == kPyJets) {
+ //Check eta range first...
+ if ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet) ||
+ (eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet))
+ {
+ //Eta is okay, now check phi range
+ if ((phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet) ||
+ (phi[1] < fPhiMaxJet && phi[1] > fPhiMinJet))
+ {
+ triggered = kTRUE;
+ }
+ }
+ } else {
+ Int_t ij = 0;
+ Int_t ig = 1;
+ if (pdg[0] == kGamma) {
+ ij = 1;
+ ig = 0;
+ }
+ //Check eta range first...
+ if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
+ (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
+ {
+ //Eta is okay, now check phi range
+ if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
+ (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
+ {
+ triggered = kTRUE;
+ }
+ }
+ }
+ return triggered;
+}
AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
{