/*
$Log$
+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.
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"
:AliGenMC()
{
// Default Constructor
- fDecayer = new AliDecayerPythia();
+ fParticles = 0;
+ fPythia = 0;
+ fDecayer = new AliDecayerPythia();
SetEventListRange();
+ SetJetPhiRange();
+ SetJetEtaRange();
}
AliGenPythia::AliGenPythia(Int_t npart)
// semimuonic decay
// structure function GRVHO
//
+ fName = "Pythia";
+ fTitle= "Particle Generator using PYTHIA";
fXsection = 0.;
fNucA1=0;
fNucA2=0;
SetStrucFunc();
SetForceDecay();
SetPtHard();
+ SetYHard();
SetEnergyCMS();
fDecayer = new AliDecayerPythia();
// Set random number generator
sRandom=fRandom;
- SetEventListRange();
fFlavorSelect = 0;
// Produced particles
fParticles = new TClonesArray("TParticle",1000);
fEventVertex.Set(3);
+ SetEventListRange();
+ SetJetPhiRange();
+ SetJetEtaRange();
+ // 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()
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;
- fFlavorSelect = 4;
- break;
case kPyCharmUnforced:
- fParentSelect[0] = 411;
- fParentSelect[1] = 421;
- fParentSelect[2] = 431;
- fParentSelect[3]= 4122;
- fFlavorSelect = 4;
+ 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:
+ case kPyBeautyPbMNR:
fParentSelect[0]= 511;
fParentSelect[1]= 521;
fParentSelect[2]= 531;
fParentSelect[0] = 443;
break;
case kPyMb:
+ case kPyMbNonDiffr:
case kPyJets:
case kPyDirectGamma:
break;
}
+//
+// This counts the total number of calls to Pyevnt() per run.
+ fTrialsRun = 0;
+ fQ = 0.;
+ fX1 = 0.;
+ fX2 = 0.;
+ fNev = 0.;
+//
AliGenMC::Init();
}
{
// Generate one event
fDecayer->ForceDecay();
+
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
- Float_t origin0[3] = {0,0,0};
Float_t p[3];
// converts from mm/c to s
const Float_t kconv=0.001/2.999792458e8;
Int_t j, kf;
fTrials=0;
-// Set collision vertex position
- 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++) {
if (gAlice->GetEvNumber()>=fDebugEventFirst &&
gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
fTrials++;
+
fPythia->ImportParticles(fParticles,"All");
//
Int_t i;
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];
Int_t* trackIt = new Int_t[np];
- for (i=0; i< np-1; i++) {
+ for (i=0; i< np; i++) {
pParent[i] = -1;
pSelected[i] = 0;
+ trackIt[i] = 0;
}
- printf("\n **************************************************%d\n",np);
- Int_t nc = 0;
- if (np == 0 ) continue;
-// Get event vertex
- TParticle* iparticle;
- iparticle = (TParticle *) fParticles->At(0);
- fEventVertex[0] = iparticle->Vx()/10.;
- fEventVertex[1] = iparticle->Vy()/10.;
- fEventVertex[2] = iparticle->Vz()/10.;
-//
- if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
- for (i = 0; i<np-1; i++) {
+
+ 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());
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) {
+ // 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) {
//
// Heavy flavor hadron or quark
//
// Kinematic seletion on final state heavy flavor mesons
if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
{
- nc = -1;
- break;
+ 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))
{
- nc = -1;
- break;
+ continue;
}
//
// Decay products
// Semi-stable or de-selected: diselect decay products:
//
//
- if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > 10e-15)
+ if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
{
Int_t ipF = iparticle->GetFirstDaughter();
Int_t ipL = iparticle->GetLastDaughter();
}
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)> 10e-15)) 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;
} else {
if (ParentSelected(kf)) trackIt[i] = 0;
}
+ if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
//
//
} // particle selection loop
- if (nc > -1) {
- for (i = 0; i<np-1; i++) {
+ if (nc > 0) {
+ for (i = 0; i<np; 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[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();
- Int_t ipa = iparticle->GetFirstMother()-1;
+ 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;
-// printf("\n SetTrack %d %d %d %d", i, iparent, kf, trackIt[i]);
SetTrack(fTrackIt*trackIt[i] ,
- iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
+ iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
pParent[i] = nt;
KeepTrack(nt);
} // SetTrack loop
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);
+ fPythia->GetVINT(41),
+ fPythia->GetVINT(42),
+ fPythia->GetVINT(51));
+
+ fQ += fPythia->GetVINT(51);
+ fX1 += fPythia->GetVINT(41);
+ fX2 += fPythia->GetVINT(42);
+ fTrialsRun += fTrials;
+ fNev++;
+ MakeHeader();
break;
}
}
Int_t AliGenPythia::GenerateMB()
{
+//
+// Min Bias selection and other global selections
+//
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;
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++) {
+ 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();
-// 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;
-
+
iparent = (ipa > -1) ? pParent[ipa] : -1;
-
+
//
// 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.;
+ 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);
+ tof, kPPrimary, nt, 1., ks);
KeepTrack(nt);
pParent[i] = nt;
} // select particle
if (pParent) delete[] pParent;
printf("\n I've put %i particles on the stack \n",nc);
- MakeHeader();
return nc;
}
{
// 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);
+
+
}
void AliGenPythia::AdjustWeights()
}
-void AliGenPythia::MakeHeader()
+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)
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff