Access to the number of associated clusters (M.Ivanov)
[u/mrichter/AliRoot.git] / RALICE / AliCollider.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // $Id: AliCollider.cxx,v 1.9 2004/01/12 08:23:22 nick Exp $
17
18 ///////////////////////////////////////////////////////////////////////////
19 // Class AliCollider
20 // Pythia based universal physics event generator.
21 // This class is derived from TPythia6 and has some extensions to
22 // support also generation of nucleus-nucleus interactions and to allow
23 // investigation of the effect of detector resolving power.
24 // Furthermore, the produced event information is provided in a format
25 // using the AliEvent structure.
26 // For the produced AliTrack objects, the particle ID code is set to the
27 // Pythia KF value, which is compatible with the PDG identifier.
28 // This will allow a direct analysis of the produced data using the
29 // Ralice physics analysis tools.
30 //
31 // For further details concerning the produced output structure,
32 // see the docs of the memberfunctions SetVertexMode and SetResolution.
33 //
34 // Example job of minimum biased Pb+Pb interactions :
35 // --------------------------------------------------
36 // {
37 //  gSystem->Load("libEG");
38 //  gSystem->Load("libEGPythia6");
39 //  gSystem->Load("ralice");
40 //
41 //  AliCollider* gen=new AliCollider();
42 //
43 //  gen->SetOutputFile("test.root");
44 //  gen->SetVertexMode(3);    
45 //  gen->SetResolution(1e-4); // 1 micron vertex resolution
46 //
47 //  gen->SetRunNumber(1);
48 //
49 //  Int_t zp=82;
50 //  Int_t ap=208;
51 //  Int_t zt=82;
52 //  Int_t at=208;
53 //
54 //  gen->Init("fixt",zp,ap,zt,at,158);
55 //
56 //  gen->SetTitle("SPS Pb-Pb collision at 158A GeV/c beam energy");
57 //
58 //  Int_t nevents=5;
59 //
60 //  AliRandom rndm;
61 //  Float_t* rans=new Float_t[nevents];
62 //  rndm.Uniform(rans,nevents,2,ap+at);
63 //  Int_t npart;
64 //  for (Int_t i=0; i<nevents; i++)
65 //  {
66 //   npart=rans[i];
67 //   gen->MakeEvent(npart);
68 //
69 //   AliEvent* evt=gen->GetEvent();
70 //  
71 //   evt->List();
72 //  }
73 //
74 //  gen->EndRun();
75 // }
76 //
77 //
78 // Example job of a cosmic nu+p atmospheric interaction.
79 // -----------------------------------------------------
80 // {
81 //  gSystem->Load("libEG");
82 //  gSystem->Load("libEGPythia6");
83 //  gSystem->Load("ralice");
84 //
85 //  AliCollider* gen=new AliCollider();
86 //
87 //  gen->SetOutputFile("test.root");
88 //
89 //  gen->SetRunNumber(1);
90 //
91 //  gen->Init("fixt","nu_mu","p",1e11);
92 //
93 //  gen->SetTitle("Atmospheric nu_mu-p interaction at 1e20 eV");
94 //
95 //  Int_t nevents=10;
96 //
97 //  for (Int_t i=0; i<nevents; i++)
98 //  {
99 //   gen->MakeEvent(0,1);
100 //
101 //   AliEvent* evt=gen->GetEvent();
102 //  
103 //   evt->Data();
104 //  }
105 //
106 //  gen->EndRun();
107 // }
108 //
109 //
110 //--- Author: Nick van Eijndhoven 22-nov-2002 Utrecht University
111 //- Modified: NvE $Date: 2004/01/12 08:23:22 $ Utrecht University
112 ///////////////////////////////////////////////////////////////////////////
113
114 #include "AliCollider.h"
115 #include "Riostream.h"
116  
117 ClassImp(AliCollider) // Class implementation to enable ROOT I/O
118  
119 AliCollider::AliCollider() : TPythia6()
120 {
121 // Default constructor.
122 // All variables initialised to default values.
123  fVertexmode=0;    // No vertex structure creation
124  fResolution=1e-5; // Standard resolution is 0.1 micron
125  fRunnum=0;
126  fEventnum=0;
127  fPrintfreq=1;
128
129  fEvent=0;
130
131  fSpecpmin=0;
132
133  fFrame="none";
134  fWin=0;
135
136  fNucl=0;
137  fZproj=0;
138  fAproj=0;
139  fZtarg=0;
140  fAtarg=0;
141  fFracpp=0;
142  fFracnp=0;
143  fFracpn=0;
144  fFracnn=0;
145
146  fOutFile=0;
147  fOutTree=0;
148
149  fSelections=0;
150  fSelect=0;
151
152  TString s=GetName();
153  s+=" (AliCollider)";
154  SetName(s.Data());
155 }
156 ///////////////////////////////////////////////////////////////////////////
157 AliCollider::~AliCollider()
158 {
159 // Default destructor
160  if (fEvent)
161  {
162   delete fEvent;
163   fEvent=0;
164  }
165  if (fOutFile)
166  {
167   delete fOutFile;
168   fOutFile=0;
169  }
170  if (fOutTree)
171  {
172   delete fOutTree;
173   fOutTree=0;
174  }
175  if (fSelections)
176  {
177   delete fSelections;
178   fSelections=0;
179  }
180 }
181 ///////////////////////////////////////////////////////////////////////////
182 void AliCollider::SetOutputFile(TString s)
183 {
184 // Create the output file containing all the data in ROOT output format.
185  if (fOutFile)
186  {
187   delete fOutFile;
188   fOutFile=0;
189  }
190  fOutFile=new TFile(s.Data(),"RECREATE","AliCollider data");
191
192  if (fOutTree)
193  {
194   delete fOutTree;
195   fOutTree=0;
196  }
197  fOutTree=new TTree("T","AliCollider event data");
198
199  Int_t bsize=32000;
200  Int_t split=0;
201  fOutTree->Branch("Events","AliEvent",&fEvent,bsize,split);
202 }
203 ///////////////////////////////////////////////////////////////////////////
204 void AliCollider::SetVertexMode(Int_t mode)
205 {
206 // Set the mode of the vertex structure creation.
207 //
208 // By default all generated tracks will only appear in the AliEvent
209 // structure without any primary (and secondary) vertex structure.
210 // The user can build the vertex structure if he/she wants by means
211 // of the beginpoint location of each AliTrack.
212 //
213 // However, one can also let AliCollider automatically create
214 // the primary (and secondary) vertex structure(s).
215 // In this case the primary vertex is given Id=1 and all sec. vertices
216 // are given Id's 2,3,4,....
217 // All vertices are created as standalone entities in the AliEvent structure
218 // without any linking between the various vertices.
219 // For this automated process, the user-selected resolution
220 // (see SetResolution) is used to decide whether or not certain vertex
221 // locations can be resolved.
222 // In case no vertex creation is selected (i.e. the default mode=0),
223 // the value of the resolution is totally irrelevant.
224 //
225 // The user can also let AliCollider automatically connect the sec. vertices
226 // to the primary vertex (i.e. mode=3). This process will also automatically
227 // generate the tracks connecting the vertices.
228 // Note that the result of the mode=3 operation may be very sensitive to
229 // the resolution parameter. Therefore, no attempt is made to distinguish
230 // between secondary, tertiary etc... vertices. All sec. vertices are
231 // linked to the primary one.
232 //  
233 // Irrespective of the selected mode, all generated tracks can be obtained
234 // directly from the AliEvent structure.
235 // In case (sec.) vertex creation is selected, all generated vertices can
236 // also be obtained directly from the AliEvent structure. 
237 // These (sec.) vertices contain only the corresponding pointers to the various
238 // tracks which are stored in the AliEvent structure.
239 //
240 // Overview of vertex creation modes :
241 // -----------------------------------
242 // mode = 0 ==> No vertex structure will be created
243 //        1 ==> Only primary vertex structure will be created
244 //        2 ==> Unconnected primary and secondary vertices will be created
245 //        3 ==> Primary and secondary vertices will be created where all the
246 //              sec. vertices will be connected to the primary vertex.
247 //              Also the vertex connecting tracks will be automatically
248 //              generated. 
249 //
250  if (mode<0 || mode >3)
251  {
252   cout << " *AliCollider::SetVertexMode* Invalid argument mode : " << mode << endl;
253   fVertexmode=0;
254  }
255  else
256  {
257   fVertexmode=mode;
258  }
259 }
260 ///////////////////////////////////////////////////////////////////////////
261 Int_t AliCollider::GetVertexMode()
262 {
263 // Provide the current mode for vertex structure creation.
264  return fVertexmode;
265 }
266 ///////////////////////////////////////////////////////////////////////////
267 void AliCollider::SetResolution(Double_t res)
268 {
269 // Set the resolution (in cm) for resolving (sec.) vertices.
270 // By default this resolution is set to 0.1 micron.
271 // Note : In case no vertex creation has been selected, the value of
272 //        the resolution is totally irrelevant.
273  fResolution=fabs(res);
274 }
275 ///////////////////////////////////////////////////////////////////////////
276 Double_t AliCollider::GetResolution()
277 {
278 // Provide the current resolution (in cm) for resolving (sec.) vertices.
279  return fResolution;
280 }
281 ///////////////////////////////////////////////////////////////////////////
282 void AliCollider::SetRunNumber(Int_t run)
283 {
284 // Set the user defined run number.
285 // By default the run number is set to 0.
286  fRunnum=run;
287 }
288 ///////////////////////////////////////////////////////////////////////////
289 Int_t AliCollider::GetRunNumber()
290 {
291 // Provide the user defined run number.
292  return fRunnum;
293 }
294 ///////////////////////////////////////////////////////////////////////////
295 void AliCollider::SetPrintFreq(Int_t n)
296 {
297 // Set the print frequency for every 'n' events.
298 // By default the printfrequency is set to 1 (i.e. every event).
299  fPrintfreq=n;
300 }
301 ///////////////////////////////////////////////////////////////////////////
302 Int_t AliCollider::GetPrintFreq()
303 {
304 // Provide the user selected print frequency.
305  return fPrintfreq;
306 }
307 ///////////////////////////////////////////////////////////////////////////
308 void AliCollider::Init(char* frame,char* beam,char* target,Float_t win)
309 {
310 // Initialisation of the underlying Pythia generator package.
311 // This routine just invokes TPythia6::Initialize(...) and the arguments
312 // have the corresponding meaning.
313 // The event number is reset to 0.
314  fEventnum=0;
315  fNucl=0;
316  fFrame=frame;
317  fWin=win;
318  Initialize(frame,beam,target,win);
319
320  cout << " *AliCollider::Init* Standard Pythia initialisation." << endl;
321  cout << " Beam particle : " << beam << " Target particle : " << target
322       << " Frame = " << frame << " Energy = " << win
323       << endl;
324 }
325 ///////////////////////////////////////////////////////////////////////////
326 void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win)
327 {
328 // Initialisation of the underlying Pythia generator package for the generation
329 // of nucleus-nucleus interactions.
330 // In addition to the Pythia standard arguments 'frame' and 'win', the user
331 // can specify here (Z,A) values of the projectile and target nuclei.
332 //
333 // Note : The 'win' value denotes either the cms energy per nucleon-nucleon collision
334 //        (i.e. frame="cms") or the momentum per nucleon in all other cases.
335 //
336 // The event number is reset to 0.
337  fEventnum=0;
338  fNucl=1;
339  fFrame=frame;
340  fWin=win;
341  fZproj=0;
342  fAproj=0;
343  fZtarg=0;
344  fAtarg=0;
345  fFracpp=0;
346  fFracnp=0;
347  fFracpn=0;
348  fFracnn=0;
349
350  if (ap<1 || at<1 || zp>ap || zt>at)
351  {
352   cout << " *AliCollider::Init* Invalid input value(s). Zproj = " << zp
353        << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
354   return;
355  }
356
357  fZproj=zp;
358  fAproj=ap;
359  fZtarg=zt;
360  fAtarg=at;
361
362  cout << " *AliCollider::Init* Nucleus-Nucleus generator initialisation." << endl;
363  cout << " Zproj = " << zp << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at
364       << " Frame = " << frame << " Energy = " << win
365       << endl;
366 }
367 ///////////////////////////////////////////////////////////////////////////
368 void AliCollider::GetFractions(Float_t zp,Float_t ap,Float_t zt,Float_t at)
369 {
370 // Determine the fractions for the various N-N collision processes.
371 // The various processes are : p+p, n+p, p+n and n+n.
372  if (zp<0) zp=0;
373  if (zt<0) zt=0;
374
375  fFracpp=0;
376  fFracnp=0;
377  fFracpn=0;
378  fFracnn=0;
379
380  if (ap>0 && at>0)
381  {
382   fFracpp=(zp/ap)*(zt/at);
383   fFracnp=(1.-zp/ap)*(zt/at);
384   fFracpn=(zp/ap)*(1.-zt/at);
385   fFracnn=(1.-zp/ap)*(1.-zt/at);
386  }
387 }
388 ///////////////////////////////////////////////////////////////////////////
389 void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
390 {
391 // Generate one event.
392 // In case of a nucleus-nucleus interaction, the argument 'npt' denotes
393 // the number of participant nucleons.
394 // Normally also the spectator tracks will be stored into the event structure.
395 // The spectator tracks have a negative user Id to distinguish them from the
396 // ordinary generated tracks.
397 // In case the user has selected the creation of vertex structures, the spectator
398 // tracks will be linked to the primary vertex.
399 // However, specification of npt<0 will suppress the storage of spectator tracks.
400 // In the latter case abs(npt) will be taken as the number of participants.  
401 // In case of a standard Pythia run for 'elementary' particle interactions,
402 // the value of npt is totally irrelevant.
403 //
404 // The argument 'mlist' denotes the list mode used for Pylist().
405 // Note : mlist<0 suppresses the invokation of Pylist().
406 // By default, no listing is produced (i.e. mlist=-1).
407 //
408 // The argument 'medit' denotes the edit mode used for Pyedit().
409 // Note : medit<0 suppresses the invokation of Pyedit().
410 // By default, only 'stable' final particles are kept (i.e. medit=1). 
411 //
412 // In the case of a standard Pythia run concerning 'elementary' particle
413 // interactions, the projectile and target particle ID's for the created
414 // event structure are set to the corresponding Pythia KF codes.
415 // All the A and Z values are in that case set to zero.
416 // In case of a nucleus-nucleus interaction, the proper A and Z values for 
417 // the projectile and target particles are set in the event structure.
418 // However, in this case both particle ID's are set to zero.
419 //
420 // Note : Only in case an event passed the selection criteria as specified
421 //        via SelectEvent(), the event will appear on the output file.
422
423  fEventnum++; 
424
425  Int_t specmode=1;
426  if (npt<0)
427  {
428   specmode=0;
429   npt=abs(npt);
430  }
431
432  // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
433  Int_t ncols[4]={0,0,0,0};
434
435  Int_t zp=0;
436  Int_t ap=0;
437  Int_t zt=0;
438  Int_t at=0;
439
440  Int_t ncol=1;
441  if (fNucl)
442  {
443   if (npt<1 || npt>(fAproj+fAtarg))
444   {
445    cout << " *AliCollider::MakeEvent* Invalid input value. npt = " << npt
446         << " Aproj = " << fAproj << " Atarg = " << fAtarg << endl;
447    return;
448   }
449
450   // Determine the number of nucleon-nucleon collisions
451   ncol=npt/2;
452   if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
453
454   // Determine the number of the various types of N+N interactions
455   zp=fZproj;
456   ap=fAproj;
457   zt=fZtarg;
458   at=fAtarg;
459   Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
460   if (ap>at) maxa=1;
461   Float_t* rans=new Float_t[ncol];
462   fRan.Uniform(rans,ncol);
463   Float_t rndm=0;
464   for (Int_t i=0; i<ncol; i++)
465   {
466    GetFractions(zp,ap,zt,at);
467    rndm=rans[i];
468    if (rndm<=fFracpp) // p+p interaction
469    {
470     ncols[0]++;
471     if (maxa==2)
472     {
473      at--;
474      zt--;
475     } 
476     else
477     {
478      ap--;
479      zp--;
480     }
481    }
482    if (rndm>fFracpp && rndm<=(fFracpp+fFracnp)) // n+p interaction
483    {
484     ncols[1]++;
485     if (maxa==2)
486     {
487      at--;
488      zt--;
489     } 
490     else
491     {
492      ap--;
493     }
494    }
495    if (rndm>(fFracpp+fFracnp) && rndm<=(fFracpp+fFracnp+fFracpn)) // p+n interaction
496    {
497     ncols[2]++;
498     if (maxa==2)
499     {
500      at--;
501     } 
502     else
503     {
504      ap--;
505      zp--;
506     }
507    }
508    if (rndm>(fFracpp+fFracnp+fFracpn)) // n+n interaction
509    {
510     ncols[3]++; 
511     if (maxa==2)
512     {
513      at--;
514     } 
515     else
516     {
517      ap--;
518     }
519    }
520   }
521   delete [] rans;
522  }
523
524  if (!(fEventnum%fPrintfreq))
525  {
526   cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
527        << endl;
528   if (fNucl)
529   {
530    cout << " npart = " << npt << " ncol = " << ncol 
531         << " ncolpp = " << ncols[0] << " ncolnp = " << ncols[1]
532         << " ncolpn = " << ncols[2] << " ncolnn = " << ncols[3] << endl;
533   }
534  }
535
536  if (!fEvent)
537  {
538   fEvent=new AliEvent();
539   fEvent->SetOwner();
540   fEvent->SetName(GetName());
541   fEvent->SetTitle(GetTitle());
542  }
543
544  fEvent->Reset();
545  fEvent->SetRunNumber(fRunnum);
546  fEvent->SetEventNumber(fEventnum);
547
548  AliTrack t;
549  Ali3Vector p;
550  AliPosition r,rx;
551  Float_t v[3];
552  AliVertex vert;
553  Ali3Vector pproj,ptarg;
554
555  if (fVertexmode)
556  {
557   // Make sure the primary vertex gets correct location and Id=1
558   v[0]=0;
559   v[1]=0;
560   v[2]=0;
561   r.SetPosition(v,"car");
562   v[0]=fResolution;
563   v[1]=fResolution;
564   v[2]=fResolution;
565   r.SetPositionErrors(v,"car");
566
567   vert.SetId(1);
568   vert.SetTrackCopy(0);
569   vert.SetVertexCopy(0);
570   vert.SetPosition(r);
571   fEvent->AddVertex(vert,0);
572  }
573
574  Int_t kf=0;
575  Float_t charge=0,mass=0;
576  TString name;
577
578  Int_t ntypes=4;
579
580  // Singular settings for a normal Pythia elementary particle interation 
581  if (!fNucl)
582  {
583   ntypes=1;
584   ncols[0]=1;
585  }
586
587  // Generate all the various collisions
588  fSelect=0;      // Flag to indicate whether the total event is selected or not 
589  Int_t select=0; // Flag to indicate whether the sub-event is selected or not 
590  Int_t first=1;  // Flag to indicate the first collision process
591  Double_t pnucl;
592  Int_t npart=0,ntk=0;
593  Double_t dist=0;
594  for (Int_t itype=0; itype<ntypes; itype++)
595  {
596   if (fNucl)
597   {
598    if (itype==0 && ncols[itype]) Initialize(fFrame,"p","p",fWin);
599    if (itype==1 && ncols[itype]) Initialize(fFrame,"n","p",fWin);
600    if (itype==2 && ncols[itype]) Initialize(fFrame,"p","n",fWin);
601    if (itype==3 && ncols[itype]) Initialize(fFrame,"n","n",fWin);
602   }
603   for (Int_t jcol=0; jcol<ncols[itype]; jcol++)
604   {
605    GenerateEvent();
606
607    select=IsSelected();
608    if (select) fSelect=1;
609
610    if (first) // Store projectile and target information in the event structure
611    {
612     if (fNucl)
613     {
614      v[0]=GetP(1,1);
615      v[1]=GetP(1,2);
616      v[2]=GetP(1,3);
617      pproj.SetVector(v,"car");
618      pnucl=pproj.GetNorm();
619      fEvent->SetProjectile(fAproj,fZproj,pnucl);
620      v[0]=GetP(2,1);
621      v[1]=GetP(2,2);
622      v[2]=GetP(2,3);
623      ptarg.SetVector(v,"car");
624      pnucl=ptarg.GetNorm();
625      fEvent->SetTarget(fAtarg,fZtarg,pnucl);
626     }
627     else
628     {
629      v[0]=GetP(1,1);
630      v[1]=GetP(1,2);
631      v[2]=GetP(1,3);
632      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
633      kf=GetK(1,2);
634      fEvent->SetProjectile(0,0,pnucl,kf);
635      v[0]=GetP(2,1);
636      v[1]=GetP(2,2);
637      v[2]=GetP(2,3);
638      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
639      kf=GetK(2,2);
640      fEvent->SetTarget(0,0,pnucl,kf);
641     }
642     first=0;
643    }
644
645    if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
646
647    if (mlist>=0 && select) Pylist(mlist);
648
649    npart=GetN();
650    for (Int_t jpart=1; jpart<=npart; jpart++)
651    {
652     kf=GetK(jpart,2);
653     charge=Pychge(kf)/3.;
654     mass=GetP(jpart,5);
655     name=GetPyname(kf);
656
657     // 3-momentum in GeV/c
658     v[0]=GetP(jpart,1);
659     v[1]=GetP(jpart,2);
660     v[2]=GetP(jpart,3);
661     p.SetVector(v,"car");
662
663     // Production location in cm.
664     v[0]=GetV(jpart,1)/10;
665     v[1]=GetV(jpart,2)/10;
666     v[2]=GetV(jpart,3)/10;
667     r.SetPosition(v,"car");
668
669     ntk++;
670
671     t.Reset();
672     t.SetId(ntk);
673     t.SetParticleCode(kf);
674     t.SetName(name.Data());
675     t.SetCharge(charge);
676     t.SetMass(mass);
677     t.Set3Momentum(p);
678     t.SetBeginPoint(r);
679
680     fEvent->AddTrack(t);
681
682     // Build the vertex structures if requested
683     if (fVertexmode)
684     {
685      // Check if track belongs within the resolution to an existing vertex
686      Int_t add=0;  
687      for (Int_t jv=1; jv<=fEvent->GetNvertices(); jv++)
688      {
689       AliVertex* vx=fEvent->GetVertex(jv);
690       if (vx)
691       {
692        rx=vx->GetPosition();
693        dist=rx.GetDistance(r);
694        if (dist < fResolution)
695        {
696         AliTrack* tx=fEvent->GetIdTrack(ntk);
697         if (tx)
698         { 
699          vx->AddTrack(tx);
700          add=1;
701         }
702        }
703       }
704       if (add) break; // No need to look further for vertex candidates
705      }
706
707      // If track was not close enough to an existing vertex
708      // a new secondary vertex is created      
709      if (!add && fVertexmode>1)
710      {
711       AliTrack* tx=fEvent->GetIdTrack(ntk);
712       if (tx)
713       {
714        v[0]=fResolution;
715        v[1]=fResolution;
716        v[2]=fResolution;
717        r.SetPositionErrors(v,"car");
718        vert.Reset();
719        vert.SetTrackCopy(0);
720        vert.SetVertexCopy(0);
721        vert.SetId((fEvent->GetNvertices())+1);
722        vert.SetPosition(r);
723        vert.AddTrack(tx);
724        fEvent->AddVertex(vert,0);
725       } 
726      }
727     }
728    } // End of loop over the produced particles for each collision
729   } // End of loop over number of collisions for each type
730  } // End of loop over collision types
731
732  // Link sec. vertices to primary if requested
733  // Note that also the connecting tracks are automatically created
734  if (fVertexmode>2)
735  {
736   AliVertex* vp=fEvent->GetIdVertex(1); // Primary vertex
737   if (vp)
738   {
739    for (Int_t i=2; i<=fEvent->GetNvertices(); i++)
740    {
741     AliVertex* vx=fEvent->GetVertex(i);
742     if (vx)
743     {
744      if (vx->GetId() != 1) vp->AddVertex(vx);
745     }
746    }
747   }
748  }
749
750  // Include the spectator tracks in the event structure.
751  if (fNucl && specmode)
752  {
753   v[0]=0;
754   v[1]=0;
755   v[2]=0;
756   r.SetPosition(v,"car");
757
758   zp=fZproj-(ncols[0]+ncols[2]);
759   if (zp<0) zp=0;
760   ap=fAproj-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
761   if (ap<0) ap=0;
762   zt=fZtarg-(ncols[0]+ncols[1]);
763   if (zt<0) zt=0;
764   at=fAtarg-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
765   if (at<0) at=0;
766
767   Int_t nspec=0;
768
769   if (pproj.GetNorm() > fSpecpmin)
770   {
771    kf=2212; // Projectile spectator protons
772    charge=Pychge(kf)/3.;
773    mass=GetPMAS(Pycomp(kf),1);
774    name=GetPyname(kf);
775    for (Int_t iprojp=1; iprojp<=zp; iprojp++)
776    {
777     nspec++;
778     t.Reset();
779     t.SetId(-nspec);
780     t.SetParticleCode(kf);
781     t.SetName(name.Data());
782     t.SetTitle("Projectile spectator proton");
783     t.SetCharge(charge);
784     t.SetMass(mass);
785     t.Set3Momentum(pproj);
786     t.SetBeginPoint(r);
787
788     fEvent->AddTrack(t);
789    }
790
791    kf=2112; // Projectile spectator neutrons
792    charge=Pychge(kf)/3.;
793    mass=GetPMAS(Pycomp(kf),1);
794    name=GetPyname(kf);
795    for (Int_t iprojn=1; iprojn<=(ap-zp); iprojn++)
796    {
797     nspec++;
798     t.Reset();
799     t.SetId(-nspec);
800     t.SetParticleCode(kf);
801     t.SetName(name.Data());
802     t.SetTitle("Projectile spectator neutron");
803     t.SetCharge(charge);
804     t.SetMass(mass);
805     t.Set3Momentum(pproj);
806     t.SetBeginPoint(r);
807
808     fEvent->AddTrack(t);
809    }
810   }
811
812   if (ptarg.GetNorm() > fSpecpmin)
813   {
814    kf=2212; // Target spectator protons
815    charge=Pychge(kf)/3.;
816    mass=GetPMAS(Pycomp(kf),1);
817    name=GetPyname(kf);
818    for (Int_t itargp=1; itargp<=zt; itargp++)
819    {
820     nspec++;
821     t.Reset();
822     t.SetId(-nspec);
823     t.SetParticleCode(kf);
824     t.SetName(name.Data());
825     t.SetTitle("Target spectator proton");
826     t.SetCharge(charge);
827     t.SetMass(mass);
828     t.Set3Momentum(ptarg);
829     t.SetBeginPoint(r);
830
831     fEvent->AddTrack(t);
832    }
833
834    kf=2112; // Target spectator neutrons
835    charge=Pychge(kf)/3.;
836    mass=GetPMAS(Pycomp(kf),1);
837    name=GetPyname(kf);
838    for (Int_t itargn=1; itargn<=(at-zt); itargn++)
839    {
840     nspec++;
841     t.Reset();
842     t.SetId(-nspec);
843     t.SetParticleCode(kf);
844     t.SetName(name.Data());
845     t.SetTitle("Target spectator neutron");
846     t.SetCharge(charge);
847     t.SetMass(mass);
848     t.Set3Momentum(ptarg);
849     t.SetBeginPoint(r);
850
851     fEvent->AddTrack(t);
852    }
853   }
854
855  // Link the spectator tracks to the primary vertex.
856  if (fVertexmode)
857  {
858   AliVertex* vp=fEvent->GetIdVertex(1);
859   if (vp)
860   {
861    for (Int_t ispec=1; ispec<=nspec; ispec++)
862    {
863     AliTrack* tx=fEvent->GetIdTrack(-ispec);
864     if (tx) vp->AddTrack(tx);
865    }
866   }
867  }
868 }
869
870  if (mlist && !(fEventnum%fPrintfreq)) cout << endl; // Create empty output line after the event
871
872  if (fOutTree && fSelect) fOutTree->Fill();
873 }
874 ///////////////////////////////////////////////////////////////////////////
875 AliEvent* AliCollider::GetEvent(Int_t select)
876 {
877 // Provide pointer to the generated event structure.
878 //
879 // select = 0 : Always return the pointer to the generated event.
880 //          1 : Only return the pointer to the generated event in case
881 //              the event passed the selection criteria as specified via
882 //              SelectEvent(). Otherwise the value 0 will be returned.
883 //
884 // By invoking GetEvent() the default of select=0 will be used.
885
886  if (!select || fSelect)
887  {
888   return fEvent;
889  }
890  else
891  {
892   return 0;
893  }
894 }
895 ///////////////////////////////////////////////////////////////////////////
896 void AliCollider::EndRun()
897 {
898 // Properly close the output file (if needed).
899  if (fOutFile)
900  {
901   fOutFile->Write();
902   fOutFile->Close();
903   cout << " *AliCollider::EndRun* Output file correctly closed." << endl;
904  }
905 }
906 ///////////////////////////////////////////////////////////////////////////
907 void AliCollider::SetStable(Int_t id,Int_t mode)
908 {
909 // Declare whether a particle must be regarded as stable or not.
910 // The parameter "id" indicates the Pythia KF particle code, which
911 // basically is the PDG particle identifier code.
912 // The parameter "mode" indicates the action to be taken.
913 //
914 // mode = 0 : Particle will be able to decay
915 //        1 : Particle will be regarded as stable.
916 //
917 // In case the user does NOT explicitly invoke this function, the standard
918 // Pythia settings for the decay tables are used.
919 //
920 // When this function is invoked without the "mode" argument, then the
921 // default of mode=1 will be used for the specified particle.
922 //
923 // Notes :
924 // -------
925 // 1) This function should be invoked after the initialisation call
926 //    to AliCollider::Init.
927 // 2) Due to the internals of Pythia, there is no need to specify particles
928 //    and their corresponding anti-particles separately as (un)stable.
929 //    Once a particle has been declared (un)stable, the corresponding 
930 //    anti-particle will be treated in the same way.
931
932  if (mode==0 || mode==1)
933  {
934   Int_t kc=Pycomp(id);
935   Int_t decay=1-mode;
936   if (kc>0)
937   {
938    SetMDCY(kc,1,decay);
939   }
940   else 
941   {
942    cout << " *AliCollider::SetStable* Unknown particle code. id = " << id << endl;
943   }
944  }
945  else
946  {
947   cout << " *AliCollider::SetStable* Invalid parameter. mode = " << mode << endl;
948  }
949 }
950 ///////////////////////////////////////////////////////////////////////////
951 void AliCollider::SelectEvent(Int_t id)
952 {
953 // Add a particle to the event selection list.
954 // The parameter "id" indicates the Pythia KF particle code, which
955 // basically is the PDG particle identifier code.
956 // In case the user has built a selection list via this procedure, only the
957 // events in which one of the particles specified in the list was generated
958 // will be kept. 
959 // The investigation of the generated particles takes place when the complete
960 // event is in memory, including all (shortlived) mother particles and resonances.
961 // So, the settings of the various particle decay modes have no influence on
962 // the event selection described here.
963 //
964 // If no list has been specified, all events will be accepted.  
965 //
966 // Note : id=0 will delete the selection list.
967 //
968 // Be aware of the fact that severe selection criteria (i.e. selecting only
969 // rare events) may result in long runtimes before an event sample has been
970 // obtained.
971 //
972  if (!id)
973  {
974   if (fSelections)
975   {
976    delete fSelections;
977    fSelections=0;
978   }
979  }
980  else
981  {
982   Int_t kc=Pycomp(id);
983   if (!fSelections)
984   {
985    fSelections=new TArrayI(1);
986    fSelections->AddAt(kc,0);
987   }
988   else
989   {
990    Int_t exist=0;
991    Int_t size=fSelections->GetSize();
992    for (Int_t i=0; i<size; i++)
993    {
994     if (kc==fSelections->At(i))
995     {
996      exist=1;
997      break;
998     }
999    }
1000   
1001    if (!exist)
1002    {
1003     fSelections->Set(size+1);
1004     fSelections->AddAt(kc,size);
1005    }
1006   }
1007  }
1008 }
1009 ///////////////////////////////////////////////////////////////////////////
1010 Int_t AliCollider::GetSelectionFlag()
1011 {
1012 // Return the value of the selection flag for the total event.
1013 // When the event passed the selection criteria as specified via
1014 // SelectEvent() the value 1 is returned, otherwise the value 0 is returned.
1015  return fSelect;
1016 }
1017 ///////////////////////////////////////////////////////////////////////////
1018 Int_t AliCollider::IsSelected()
1019 {
1020 // Check whether the generated (sub)event contains one of the particles
1021 // specified in the selection list via SelectEvent().
1022 // If this is the case or when no selection list is present, the value 1
1023 // will be returned, indicating the event is selected to be kept.
1024 // Otherwise the value 0 will be returned.
1025
1026  if (!fSelections) return 1; 
1027
1028  Int_t nsel=fSelections->GetSize();
1029  Int_t npart=GetN();
1030  Int_t kf,kc;
1031
1032  Int_t select=0;
1033  for (Int_t jpart=1; jpart<=npart; jpart++)
1034  {
1035   kf=GetK(jpart,2);
1036   kc=Pycomp(kf);
1037   for (Int_t i=0; i<nsel; i++)
1038   {
1039    if (kc==fSelections->At(i))
1040    {
1041     select=1;
1042     break;
1043    }
1044   }
1045   if (select) break;
1046  }
1047  return select;
1048 }
1049 ///////////////////////////////////////////////////////////////////////////
1050 void AliCollider::SetSpectatorPmin(Float_t pmin)
1051 {
1052 // Set minimal momentum in GeV/c for spectator tracks to be stored.
1053 // Spectator tracks with a momentum below this threshold will not be stored
1054 // in the (output) event structure.
1055 // This facility allows to minimise the output file size.
1056 // Note that when the user wants to boost the event into another reference
1057 // frame these spectator tracks might have got momenta above the threshold.
1058 // However, when the spectator tracks were not stored in the event structure
1059 // in the original frame, there is no way to retreive them anymore. 
1060  fSpecpmin=pmin;
1061 }
1062 ///////////////////////////////////////////////////////////////////////////
1063 Float_t AliCollider::GetSpectatorPmin()
1064 {
1065 // Provide the minimal spectator momentum in GeV/c.
1066  return fSpecpmin;
1067 }
1068 ///////////////////////////////////////////////////////////////////////////
1069 TString AliCollider::GetPyname(Int_t kf)
1070 {
1071 // Provide the correctly truncated Pythia particle name for PGD code kf
1072 //
1073 // The TPythia6::Pyname returned name is copied into a TString and truncated
1074 // at the first blank to prevent funny trailing characters due to incorrect
1075 // stripping of empty characters in TPythia6::Pyname.
1076 // The truncation at the first blank is allowed due to the Pythia convention
1077 // that particle names never contain blanks.
1078  char name[16];
1079  TString sname;
1080  Pyname(kf,name);
1081  sname=name[0];
1082  for (Int_t i=1; i<16; i++)
1083  {
1084   if (name[i]==' ') break;
1085   sname=sname+name[i];
1086  }
1087  return sname;
1088 }
1089 ///////////////////////////////////////////////////////////////////////////