]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - RALICE/AliCollider.cxx
Use of appropriate sensor depending response objects in SPD simulation
[u/mrichter/AliRoot.git] / RALICE / AliCollider.cxx
index 6882859ed513c1814ce4430eaa287c4ee50f002b..a494d06f44e4bc35713b574c477641d30c7eff33 100644 (file)
@@ -1,4 +1,19 @@
-// $Id: AliCollider.cxx,v 1.1 2002/11/27 21:25:52 nick Exp $
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// $Id: AliCollider.cxx,v 1.12 2004/05/04 15:33:04 nick Exp $
 
 ///////////////////////////////////////////////////////////////////////////
 // Class AliCollider
@@ -38,6 +53,8 @@
 //
 //  gen->Init("fixt",zp,ap,zt,at,158);
 //
+//  gen->SetTitle("SPS Pb-Pb collision at 158A GeV/c beam energy");
+//
 //  Int_t nevents=5;
 //
 //  AliRandom rndm;
@@ -73,6 +90,8 @@
 //
 //  gen->Init("fixt","nu_mu","p",1e11);
 //
+//  gen->SetTitle("Atmospheric nu_mu-p interaction at 1e20 eV");
+//
 //  Int_t nevents=10;
 //
 //  for (Int_t i=0; i<nevents; i++)
 //
 //   AliEvent* evt=gen->GetEvent();
 //  
-//   evt->Info();
+//   evt->Data();
 //  }
 //
 //  gen->EndRun();
 //
 //
 //--- Author: Nick van Eijndhoven 22-nov-2002 Utrecht University
-//- Modified: NvE $Date: 2002/11/27 21:25:52 $ Utrecht University
+//- Modified: NvE $Date: 2004/05/04 15:33:04 $ Utrecht University
 ///////////////////////////////////////////////////////////////////////////
 
 #include "AliCollider.h"
+#include "Riostream.h"
  
 ClassImp(AliCollider) // Class implementation to enable ROOT I/O
  
-AliCollider::AliCollider()
+AliCollider::AliCollider() : TPythia6()
 {
 // Default constructor.
 // All variables initialised to default values.
+//
+// Some Pythia default MC parameters are automatically modified to provide
+// more suitable running conditions for soft processes in view of
+// nucleus-nucleus interactions and astrophysical processes.
+// The user may initialise the generator with all the default Pythia
+// parameters and obtain full user control to modify the settings by means
+// of the SetUserControl memberfunction.
+//
+// Refer to the SetElastic memberfunction for the inclusion of elastic
+// and diffractive processes.
+// By default these processes are not included.
+
  fVertexmode=0;    // No vertex structure creation
  fResolution=1e-5; // Standard resolution is 0.1 micron
  fRunnum=0;
  fEventnum=0;
  fPrintfreq=1;
+ fUserctrl=0; // Automatic optimisation of some MC parameters 
+ fElastic=0;  // No elastic and diffractive processes
 
  fEvent=0;
 
+ fSpecpmin=0;
+
  fFrame="none";
  fWin=0;
 
@@ -123,6 +159,13 @@ AliCollider::AliCollider()
 
  fOutFile=0;
  fOutTree=0;
+
+ fSelections=0;
+ fSelect=0;
+
+ TString s=GetName();
+ s+=" (AliCollider)";
+ SetName(s.Data());
 }
 ///////////////////////////////////////////////////////////////////////////
 AliCollider::~AliCollider()
@@ -143,6 +186,11 @@ AliCollider::~AliCollider()
   delete fOutTree;
   fOutTree=0;
  }
+ if (fSelections)
+ {
+  delete fSelections;
+  fSelections=0;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliCollider::SetOutputFile(TString s)
@@ -224,7 +272,7 @@ void AliCollider::SetVertexMode(Int_t mode)
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCollider::GetVertexMode()
+Int_t AliCollider::GetVertexMode() const
 {
 // Provide the current mode for vertex structure creation.
  return fVertexmode;
@@ -239,7 +287,7 @@ void AliCollider::SetResolution(Double_t res)
  fResolution=fabs(res);
 }
 ///////////////////////////////////////////////////////////////////////////
-Double_t AliCollider::GetResolution()
+Double_t AliCollider::GetResolution() const
 {
 // Provide the current resolution (in cm) for resolving (sec.) vertices.
  return fResolution;
@@ -252,7 +300,7 @@ void AliCollider::SetRunNumber(Int_t run)
  fRunnum=run;
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCollider::GetRunNumber()
+Int_t AliCollider::GetRunNumber() const
 {
 // Provide the user defined run number.
  return fRunnum;
@@ -265,33 +313,121 @@ void AliCollider::SetPrintFreq(Int_t n)
  fPrintfreq=n;
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCollider::GetPrintFreq()
+Int_t AliCollider::GetPrintFreq() const
 {
 // Provide the user selected print frequency.
  return fPrintfreq;
 }
 ///////////////////////////////////////////////////////////////////////////
+void AliCollider::SetUserControl(Int_t flag)
+{
+// Set the user control flag w.r.t. disabling automatic optimisation
+// of some Pythia default MC parameters for soft interactions in view of
+// nucleus-nucleus collisions and astrophysical processes.
+// Flag = 0 : Limited user control (automatic optimisation enabled)
+//        1 : Full user control (automatic optimisation disabled)
+// By default the user control is set to 0 (i.e. automatic optimisation).
+// See the Init() memberfunctions for further details w.r.t. the optimisations.
+ fUserctrl=flag;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCollider::GetUserControl() const
+{
+// Provide the value of the user control flag.
+ return fUserctrl;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCollider::SetElastic(Int_t flag)
+{
+// Set the flag w.r.t. inclusion of elastic and diffractive processes.
+// By default these processes are not included.
+// Flag = 0 : Do not include elastic and diffractive processes
+//        1 : Elastic and diffractive processes will be included
+ fElastic=flag;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCollider::GetElastic() const
+{
+// Provide the value of the control flag for elastic and diffractive processes.
+ return fElastic;
+}
+///////////////////////////////////////////////////////////////////////////
 void AliCollider::Init(char* frame,char* beam,char* target,Float_t win)
 {
 // Initialisation of the underlying Pythia generator package.
+// The event number is reset to 0.
 // This routine just invokes TPythia6::Initialize(...) and the arguments
 // have the corresponding meaning.
-// The event number is reset to 0.
+// Some Pythia default MC parameters are automatically modified to provide
+// more suitable running conditions for soft processes in view of
+// astrophysical processes.
+// The optimisations consist of : 
+// * Usage of real photons for photon beams of targets
+// * Minimum CMS energy of 3 GeV for the event
+// * Activation of the default K factor values
+//   with separate settings for ordinary and color annihilation graphs.
+// The user may initialise the generator with all the default Pythia
+// parameters and obtain full user control to modify the settings by means
+// of invoking the SetUserControl memberfunction before this initialisation.
+// Note that the inclusion of elastic and diffractive processes is controlled
+// by invokation of the SetElastic memberfunction before this initialisation,
+// irrespective of the UserControl selection.
+
+ if (!fUserctrl) // Optimisation of some MC parameters
+ {
+  SetMSTP(14,10); // Real photons for photon beams or targets
+  SetPARP(2,3.);  // Minimum CMS energy for the event
+  SetMSTP(33,2);  // Activate K factor. Separate for ordinary and color annih. graphs
+ }
+
+ if (fElastic) SetMSEL(2); // Include low-Pt, elastic and diffractive events
+
  fEventnum=0;
  fNucl=0;
  fFrame=frame;
  fWin=win;
  Initialize(frame,beam,target,win);
+
+ cout << endl;
+ cout << " *AliCollider::Init* Standard Pythia initialisation." << endl;
+ cout << " Beam particle : " << beam << " Target particle : " << target
+      << " Frame = " << frame << " Energy = " << win
+      << endl;
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win)
 {
 // Initialisation of the underlying Pythia generator package for the generation
 // of nucleus-nucleus interactions.
-// In addition to the Pythia standard arguments 'frame' and 'win', the user
-// can specify here (Z,A) values of the projectile and target nuclei and the number
-// 'npart' of the participant nucleons for this collision.
 // The event number is reset to 0.
+// In addition to the Pythia standard arguments 'frame' and 'win', the user
+// can specify here (Z,A) values of the projectile and target nuclei.
+//
+// Note : The 'win' value denotes either the cms energy per nucleon-nucleon collision
+//        (i.e. frame="cms") or the momentum per nucleon in all other cases.
+//
+// Some Pythia default MC parameters are automatically modified to provide
+// more suitable running conditions for soft processes in view of
+// nucleus-nucleus interactions and astrophysical processes.
+// The optimisations consist of : 
+// * Minimum CMS energy of 3 GeV for the event
+// * Activation of the default K factor values
+//   with separate settings for ordinary and color annihilation graphs.
+// The user may initialise the generator with all the default Pythia
+// parameters and obtain full user control to modify the settings by means
+// of invoking the SetUserControl memberfunction before this initialisation.
+// Note that the inclusion of elastic and diffractive processes is controlled
+// by invokation of the SetElastic memberfunction before this initialisation,
+// irrespective of the UserControl selection.
+
+ if (!fUserctrl) // Optimisation of some MC parameters
+ {
+  SetPARP(2,3.);  // Minimum CMS energy for the event
+  SetMSTP(33,2);  // Activate K factor. Separate for ordinary and color annih. graphs
+ }
+
+ if (fElastic) SetMSEL(2); // Include low-Pt, elastic and diffractive events
+
  fEventnum=0;
  fNucl=1;
  fFrame=frame;
@@ -307,6 +443,7 @@ void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t w
 
  if (ap<1 || at<1 || zp>ap || zt>at)
  {
+  cout << endl;
   cout << " *AliCollider::Init* Invalid input value(s). Zproj = " << zp
        << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
   return;
@@ -317,6 +454,7 @@ void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t w
  fZtarg=zt;
  fAtarg=at;
 
+ cout << endl;
  cout << " *AliCollider::Init* Nucleus-Nucleus generator initialisation." << endl;
  cout << " Zproj = " << zp << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at
       << " Frame = " << frame << " Energy = " << win
@@ -349,18 +487,53 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
 // Generate one event.
 // In case of a nucleus-nucleus interaction, the argument 'npt' denotes
 // the number of participant nucleons.
+// Normally also the spectator tracks will be stored into the event structure.
+// The spectator tracks have a negative user Id to distinguish them from the
+// ordinary generated tracks.
+// In case the user has selected the creation of vertex structures, the spectator
+// tracks will be linked to the primary vertex.
+// However, specification of npt<0 will suppress the storage of spectator tracks.
+// In the latter case abs(npt) will be taken as the number of participants.  
 // In case of a standard Pythia run for 'elementary' particle interactions,
 // the value of npt is totally irrelevant.
+//
+// The argument 'mlist' denotes the list mode used for Pylist().
+// Note : mlist<0 suppresses the invokation of Pylist().
+// By default, no listing is produced (i.e. mlist=-1).
+//
 // The argument 'medit' denotes the edit mode used for Pyedit().
+// Note : medit<0 suppresses the invokation of Pyedit().
 // By default, only 'stable' final particles are kept (i.e. medit=1). 
-// The argument 'mlist' denotes the list mode used for Pylist().
-// By default, no listing is produced (i.e. mlist=0).
+//
+// In the case of a standard Pythia run concerning 'elementary' particle
+// interactions, the projectile and target particle ID's for the created
+// event structure are set to the corresponding Pythia KF codes.
+// All the A and Z values are in that case set to zero.
+// In case of a nucleus-nucleus interaction, the proper A and Z values for 
+// the projectile and target particles are set in the event structure.
+// However, in this case both particle ID's are set to zero.
+//
+// Note : Only in case an event passed the selection criteria as specified
+//        via SelectEvent(), the event will appear on the output file.
 
  fEventnum++; 
 
+ Int_t specmode=1;
+ if (npt<0)
+ {
+  specmode=0;
+  npt=abs(npt);
+ }
+
  // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
  Int_t ncols[4]={0,0,0,0};
 
+ Int_t zp=0;
+ Int_t ap=0;
+ Int_t zt=0;
+ Int_t at=0;
+
+ Int_t ncol=1;
  if (fNucl)
  {
   if (npt<1 || npt>(fAproj+fAtarg))
@@ -371,14 +544,14 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
   }
 
   // Determine the number of nucleon-nucleon collisions
-  Int_t ncol=npt/2.;
+  ncol=npt/2;
   if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
 
   // Determine the number of the various types of N+N interactions
-  Int_t zp=fZproj;
-  Int_t ap=fAproj;
-  Int_t zt=fZtarg;
-  Int_t at=fAtarg;
+  zp=fZproj;
+  ap=fAproj;
+  zt=fZtarg;
+  at=fAtarg;
   Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
   if (ap>at) maxa=1;
   Float_t* rans=new Float_t[ncol];
@@ -442,47 +615,61 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
    }
   }
   delete [] rans;
+ }
 
-  if (!(fEventnum%fPrintfreq))
+ if (!(fEventnum%fPrintfreq))
+ {
+  cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
+       << endl;
+  if (fNucl)
   {
-   cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
-        << endl;
    cout << " npart = " << npt << " ncol = " << ncol 
         << " ncolpp = " << ncols[0] << " ncolnp = " << ncols[1]
         << " ncolpn = " << ncols[2] << " ncolnn = " << ncols[3] << endl;
   }
-
  }
 
  if (!fEvent)
  {
   fEvent=new AliEvent();
   fEvent->SetOwner();
+  fEvent->SetName(GetName());
+  fEvent->SetTitle(GetTitle());
  }
 
  fEvent->Reset();
  fEvent->SetRunNumber(fRunnum);
  fEvent->SetEventNumber(fEventnum);
 
+ AliTrack t;
+ Ali3Vector p;
+ AliPosition r,rx;
+ Float_t v[3];
  AliVertex vert;
+ Ali3Vector pproj,ptarg;
+
  if (fVertexmode)
  {
   // Make sure the primary vertex gets correct location and Id=1
+  v[0]=0;
+  v[1]=0;
+  v[2]=0;
+  r.SetPosition(v,"car");
+  v[0]=fResolution;
+  v[1]=fResolution;
+  v[2]=fResolution;
+  r.SetPositionErrors(v,"car");
+
   vert.SetId(1);
   vert.SetTrackCopy(0);
   vert.SetVertexCopy(0);
+  vert.SetPosition(r);
   fEvent->AddVertex(vert,0);
  }
 
- AliTrack t;
- Ali3Vector p;
- AliPosition r,rx;
- Float_t v[3];
-
- Int_t kf=0,kc=0;
+ Int_t kf=0;
  Float_t charge=0,mass=0;
-
- TMCParticle* part=0;
+ TString name;
 
  Int_t ntypes=4;
 
@@ -494,6 +681,10 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
  }
 
  // Generate all the various collisions
+ fSelect=0;      // Flag to indicate whether the total event is selected or not 
+ Int_t select=0; // Flag to indicate whether the sub-event is selected or not 
+ Int_t first=1;  // Flag to indicate the first collision process
+ Double_t pnucl;
  Int_t npart=0,ntk=0;
  Double_t dist=0;
  for (Int_t itype=0; itype<ntypes; itype++)
@@ -509,43 +700,78 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
   {
    GenerateEvent();
 
-   Pyedit(medit); // Define which particles are to be kept
+   select=IsSelected();
+   if (select) fSelect=1;
 
-   if (mlist) Pylist(mlist);
+   if (first) // Store projectile and target information in the event structure
+   {
+    if (fNucl)
+    {
+     v[0]=GetP(1,1);
+     v[1]=GetP(1,2);
+     v[2]=GetP(1,3);
+     pproj.SetVector(v,"car");
+     pnucl=pproj.GetNorm();
+     fEvent->SetProjectile(fAproj,fZproj,pnucl);
+     v[0]=GetP(2,1);
+     v[1]=GetP(2,2);
+     v[2]=GetP(2,3);
+     ptarg.SetVector(v,"car");
+     pnucl=ptarg.GetNorm();
+     fEvent->SetTarget(fAtarg,fZtarg,pnucl);
+    }
+    else
+    {
+     v[0]=GetP(1,1);
+     v[1]=GetP(1,2);
+     v[2]=GetP(1,3);
+     pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+     kf=GetK(1,2);
+     fEvent->SetProjectile(0,0,pnucl,kf);
+     v[0]=GetP(2,1);
+     v[1]=GetP(2,2);
+     v[2]=GetP(2,3);
+     pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+     kf=GetK(2,2);
+     fEvent->SetTarget(0,0,pnucl,kf);
+    }
+    first=0;
+   }
 
-   ImportParticles();
-   npart=0;
-   if (fParticles) npart=fParticles->GetEntries();
+   if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
 
-   for (Int_t jpart=0; jpart<npart; jpart++)
+   if (mlist>=0 && select)
    {
-    part=(TMCParticle*)fParticles->At(jpart);
-    if (!part) continue;
-
-    kf=part->GetKF();
-    kc=Pycomp(kf);
+    Pylist(mlist);
+    cout << endl;
+   }
 
-    charge=GetKCHG(kc,1)/3.;
-    if (kf<0) charge*=-1;
-    mass=GetPMAS(kc,1);
+   npart=GetN();
+   for (Int_t jpart=1; jpart<=npart; jpart++)
+   {
+    kf=GetK(jpart,2);
+    charge=Pychge(kf)/3.;
+    mass=GetP(jpart,5);
+    name=GetPyname(kf);
 
     // 3-momentum in GeV/c
-    v[0]=part->GetPx();
-    v[1]=part->GetPy();
-    v[2]=part->GetPz();
+    v[0]=GetP(jpart,1);
+    v[1]=GetP(jpart,2);
+    v[2]=GetP(jpart,3);
     p.SetVector(v,"car");
 
     // Production location in cm.
-    v[0]=(part->GetVx())/10;
-    v[1]=(part->GetVy())/10;
-    v[2]=(part->GetVz())/10;
-    r.SetVector(v,"car");
+    v[0]=GetV(jpart,1)/10;
+    v[1]=GetV(jpart,2)/10;
+    v[2]=GetV(jpart,3)/10;
+    r.SetPosition(v,"car");
 
     ntk++;
 
     t.Reset();
     t.SetId(ntk);
     t.SetParticleCode(kf);
+    t.SetName(name.Data());
     t.SetCharge(charge);
     t.SetMass(mass);
     t.Set3Momentum(p);
@@ -585,6 +811,10 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
       AliTrack* tx=fEvent->GetIdTrack(ntk);
       if (tx)
       {
+       v[0]=fResolution;
+       v[1]=fResolution;
+       v[2]=fResolution;
+       r.SetPositionErrors(v,"car");
        vert.Reset();
        vert.SetTrackCopy(0);
        vert.SetVertexCopy(0);
@@ -617,14 +847,158 @@ void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
   }
  }
 
- if (mlist) cout << endl; // Create empty output line after the event
- if (fOutTree) fOutTree->Fill();
+ // Include the spectator tracks in the event structure.
+ if (fNucl && specmode)
+ {
+  v[0]=0;
+  v[1]=0;
+  v[2]=0;
+  r.SetPosition(v,"car");
+
+  zp=fZproj-(ncols[0]+ncols[2]);
+  if (zp<0) zp=0;
+  ap=fAproj-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
+  if (ap<0) ap=0;
+  zt=fZtarg-(ncols[0]+ncols[1]);
+  if (zt<0) zt=0;
+  at=fAtarg-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
+  if (at<0) at=0;
+
+  Int_t nspec=0;
+
+  if (pproj.GetNorm() > fSpecpmin)
+  {
+   kf=2212; // Projectile spectator protons
+   charge=Pychge(kf)/3.;
+   mass=GetPMAS(Pycomp(kf),1);
+   name=GetPyname(kf);
+   for (Int_t iprojp=1; iprojp<=zp; iprojp++)
+   {
+    nspec++;
+    t.Reset();
+    t.SetId(-nspec);
+    t.SetParticleCode(kf);
+    t.SetName(name.Data());
+    t.SetTitle("Projectile spectator proton");
+    t.SetCharge(charge);
+    t.SetMass(mass);
+    t.Set3Momentum(pproj);
+    t.SetBeginPoint(r);
+
+    fEvent->AddTrack(t);
+   }
+
+   kf=2112; // Projectile spectator neutrons
+   charge=Pychge(kf)/3.;
+   mass=GetPMAS(Pycomp(kf),1);
+   name=GetPyname(kf);
+   for (Int_t iprojn=1; iprojn<=(ap-zp); iprojn++)
+   {
+    nspec++;
+    t.Reset();
+    t.SetId(-nspec);
+    t.SetParticleCode(kf);
+    t.SetName(name.Data());
+    t.SetTitle("Projectile spectator neutron");
+    t.SetCharge(charge);
+    t.SetMass(mass);
+    t.Set3Momentum(pproj);
+    t.SetBeginPoint(r);
+
+    fEvent->AddTrack(t);
+   }
+  }
+
+  if (ptarg.GetNorm() > fSpecpmin)
+  {
+   kf=2212; // Target spectator protons
+   charge=Pychge(kf)/3.;
+   mass=GetPMAS(Pycomp(kf),1);
+   name=GetPyname(kf);
+   for (Int_t itargp=1; itargp<=zt; itargp++)
+   {
+    nspec++;
+    t.Reset();
+    t.SetId(-nspec);
+    t.SetParticleCode(kf);
+    t.SetName(name.Data());
+    t.SetTitle("Target spectator proton");
+    t.SetCharge(charge);
+    t.SetMass(mass);
+    t.Set3Momentum(ptarg);
+    t.SetBeginPoint(r);
+
+    fEvent->AddTrack(t);
+   }
+
+   kf=2112; // Target spectator neutrons
+   charge=Pychge(kf)/3.;
+   mass=GetPMAS(Pycomp(kf),1);
+   name=GetPyname(kf);
+   for (Int_t itargn=1; itargn<=(at-zt); itargn++)
+   {
+    nspec++;
+    t.Reset();
+    t.SetId(-nspec);
+    t.SetParticleCode(kf);
+    t.SetName(name.Data());
+    t.SetTitle("Target spectator neutron");
+    t.SetCharge(charge);
+    t.SetMass(mass);
+    t.Set3Momentum(ptarg);
+    t.SetBeginPoint(r);
+
+    fEvent->AddTrack(t);
+   }
+  }
+
+ // Link the spectator tracks to the primary vertex.
+ if (fVertexmode)
+ {
+  AliVertex* vp=fEvent->GetIdVertex(1);
+  if (vp)
+  {
+   for (Int_t ispec=1; ispec<=nspec; ispec++)
+   {
+    AliTrack* tx=fEvent->GetIdTrack(-ispec);
+    if (tx) vp->AddTrack(tx);
+   }
+  }
+ }
+}
+
+ if (!(fEventnum%fPrintfreq) && (mlist || fEvent))
+ {
+  if (fEvent)
+  {
+   cout << " Number of tracks in the event structure : "
+        << fEvent->GetNtracks() << endl;
+  }
+  cout << endl; // Create empty output line after the event
+ }
+
+ if (fOutTree && fSelect) fOutTree->Fill();
 }
 ///////////////////////////////////////////////////////////////////////////
-AliEvent* AliCollider::GetEvent()
+AliEvent* AliCollider::GetEvent(Int_t select) const
 {
 // Provide pointer to the generated event structure.
- return fEvent;
+//
+// select = 0 : Always return the pointer to the generated event.
+//          1 : Only return the pointer to the generated event in case
+//              the event passed the selection criteria as specified via
+//              SelectEvent(). Otherwise the value 0 will be returned.
+//
+// By invoking GetEvent() the default of select=0 will be used.
+
+ if (!select || fSelect)
+ {
+  return fEvent;
+ }
+ else
+ {
+  return 0;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliCollider::EndRun()
@@ -638,3 +1012,186 @@ void AliCollider::EndRun()
  }
 }
 ///////////////////////////////////////////////////////////////////////////
+void AliCollider::SetStable(Int_t id,Int_t mode)
+{
+// Declare whether a particle must be regarded as stable or not.
+// The parameter "id" indicates the Pythia KF particle code, which
+// basically is the PDG particle identifier code.
+// The parameter "mode" indicates the action to be taken.
+//
+// mode = 0 : Particle will be able to decay
+//        1 : Particle will be regarded as stable.
+//
+// In case the user does NOT explicitly invoke this function, the standard
+// Pythia settings for the decay tables are used.
+//
+// When this function is invoked without the "mode" argument, then the
+// default of mode=1 will be used for the specified particle.
+//
+// Notes :
+// -------
+// 1) This function should be invoked after the initialisation call
+//    to AliCollider::Init.
+// 2) Due to the internals of Pythia, there is no need to specify particles
+//    and their corresponding anti-particles separately as (un)stable.
+//    Once a particle has been declared (un)stable, the corresponding 
+//    anti-particle will be treated in the same way.
+
+ if (mode==0 || mode==1)
+ {
+  Int_t kc=Pycomp(id);
+  Int_t decay=1-mode;
+  if (kc>0)
+  {
+   SetMDCY(kc,1,decay);
+  }
+  else 
+  {
+   cout << " *AliCollider::SetStable* Unknown particle code. id = " << id << endl;
+  }
+ }
+ else
+ {
+  cout << " *AliCollider::SetStable* Invalid parameter. mode = " << mode << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCollider::SelectEvent(Int_t id)
+{
+// Add a particle to the event selection list.
+// The parameter "id" indicates the Pythia KF particle code, which
+// basically is the PDG particle identifier code.
+// In case the user has built a selection list via this procedure, only the
+// events in which one of the particles specified in the list was generated
+// will be kept. 
+// The investigation of the generated particles takes place when the complete
+// event is in memory, including all (shortlived) mother particles and resonances.
+// So, the settings of the various particle decay modes have no influence on
+// the event selection described here.
+//
+// If no list has been specified, all events will be accepted.  
+//
+// Note : id=0 will delete the selection list.
+//
+// Be aware of the fact that severe selection criteria (i.e. selecting only
+// rare events) may result in long runtimes before an event sample has been
+// obtained.
+//
+ if (!id)
+ {
+  if (fSelections)
+  {
+   delete fSelections;
+   fSelections=0;
+  }
+ }
+ else
+ {
+  Int_t kc=Pycomp(id);
+  if (!fSelections)
+  {
+   fSelections=new TArrayI(1);
+   fSelections->AddAt(kc,0);
+  }
+  else
+  {
+   Int_t exist=0;
+   Int_t size=fSelections->GetSize();
+   for (Int_t i=0; i<size; i++)
+   {
+    if (kc==fSelections->At(i))
+    {
+     exist=1;
+     break;
+    }
+   }
+  
+   if (!exist)
+   {
+    fSelections->Set(size+1);
+    fSelections->AddAt(kc,size);
+   }
+  }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCollider::GetSelectionFlag() const
+{
+// Return the value of the selection flag for the total event.
+// When the event passed the selection criteria as specified via
+// SelectEvent() the value 1 is returned, otherwise the value 0 is returned.
+ return fSelect;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCollider::IsSelected()
+{
+// Check whether the generated (sub)event contains one of the particles
+// specified in the selection list via SelectEvent().
+// If this is the case or when no selection list is present, the value 1
+// will be returned, indicating the event is selected to be kept.
+// Otherwise the value 0 will be returned.
+
+ if (!fSelections) return 1; 
+
+ Int_t nsel=fSelections->GetSize();
+ Int_t npart=GetN();
+ Int_t kf,kc;
+
+ Int_t select=0;
+ for (Int_t jpart=1; jpart<=npart; jpart++)
+ {
+  kf=GetK(jpart,2);
+  kc=Pycomp(kf);
+  for (Int_t i=0; i<nsel; i++)
+  {
+   if (kc==fSelections->At(i))
+   {
+    select=1;
+    break;
+   }
+  }
+  if (select) break;
+ }
+ return select;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCollider::SetSpectatorPmin(Float_t pmin)
+{
+// Set minimal momentum in GeV/c for spectator tracks to be stored.
+// Spectator tracks with a momentum below this threshold will not be stored
+// in the (output) event structure.
+// This facility allows to minimise the output file size.
+// Note that when the user wants to boost the event into another reference
+// frame these spectator tracks might have got momenta above the threshold.
+// However, when the spectator tracks were not stored in the event structure
+// in the original frame, there is no way to retreive them anymore. 
+ fSpecpmin=pmin;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCollider::GetSpectatorPmin() const
+{
+// Provide the minimal spectator momentum in GeV/c.
+ return fSpecpmin;
+}
+///////////////////////////////////////////////////////////////////////////
+TString AliCollider::GetPyname(Int_t kf)
+{
+// Provide the correctly truncated Pythia particle name for PGD code kf
+//
+// The TPythia6::Pyname returned name is copied into a TString and truncated
+// at the first blank to prevent funny trailing characters due to incorrect
+// stripping of empty characters in TPythia6::Pyname.
+// The truncation at the first blank is allowed due to the Pythia convention
+// that particle names never contain blanks.
+ char name[16];
+ TString sname;
+ Pyname(kf,name);
+ sname=name[0];
+ for (Int_t i=1; i<16; i++)
+ {
+  if (name[i]==' ') break;
+  sname=sname+name[i];
+ }
+ return sname;
+}
+///////////////////////////////////////////////////////////////////////////