033aa86303b1ff3df72386a9d5c1f1f0bcfa0d60
[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.10 2004/02/13 11:08:16 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/02/13 11:08:16 $ 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 << endl;
321  cout << " *AliCollider::Init* Standard Pythia initialisation." << endl;
322  cout << " Beam particle : " << beam << " Target particle : " << target
323       << " Frame = " << frame << " Energy = " << win
324       << endl;
325 }
326 ///////////////////////////////////////////////////////////////////////////
327 void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win)
328 {
329 // Initialisation of the underlying Pythia generator package for the generation
330 // of nucleus-nucleus interactions.
331 // In addition to the Pythia standard arguments 'frame' and 'win', the user
332 // can specify here (Z,A) values of the projectile and target nuclei.
333 //
334 // Note : The 'win' value denotes either the cms energy per nucleon-nucleon collision
335 //        (i.e. frame="cms") or the momentum per nucleon in all other cases.
336 //
337 // The event number is reset to 0.
338  fEventnum=0;
339  fNucl=1;
340  fFrame=frame;
341  fWin=win;
342  fZproj=0;
343  fAproj=0;
344  fZtarg=0;
345  fAtarg=0;
346  fFracpp=0;
347  fFracnp=0;
348  fFracpn=0;
349  fFracnn=0;
350
351  if (ap<1 || at<1 || zp>ap || zt>at)
352  {
353   cout << endl;
354   cout << " *AliCollider::Init* Invalid input value(s). Zproj = " << zp
355        << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
356   return;
357  }
358
359  fZproj=zp;
360  fAproj=ap;
361  fZtarg=zt;
362  fAtarg=at;
363
364  cout << endl;
365  cout << " *AliCollider::Init* Nucleus-Nucleus generator initialisation." << endl;
366  cout << " Zproj = " << zp << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at
367       << " Frame = " << frame << " Energy = " << win
368       << endl;
369 }
370 ///////////////////////////////////////////////////////////////////////////
371 void AliCollider::GetFractions(Float_t zp,Float_t ap,Float_t zt,Float_t at)
372 {
373 // Determine the fractions for the various N-N collision processes.
374 // The various processes are : p+p, n+p, p+n and n+n.
375  if (zp<0) zp=0;
376  if (zt<0) zt=0;
377
378  fFracpp=0;
379  fFracnp=0;
380  fFracpn=0;
381  fFracnn=0;
382
383  if (ap>0 && at>0)
384  {
385   fFracpp=(zp/ap)*(zt/at);
386   fFracnp=(1.-zp/ap)*(zt/at);
387   fFracpn=(zp/ap)*(1.-zt/at);
388   fFracnn=(1.-zp/ap)*(1.-zt/at);
389  }
390 }
391 ///////////////////////////////////////////////////////////////////////////
392 void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
393 {
394 // Generate one event.
395 // In case of a nucleus-nucleus interaction, the argument 'npt' denotes
396 // the number of participant nucleons.
397 // Normally also the spectator tracks will be stored into the event structure.
398 // The spectator tracks have a negative user Id to distinguish them from the
399 // ordinary generated tracks.
400 // In case the user has selected the creation of vertex structures, the spectator
401 // tracks will be linked to the primary vertex.
402 // However, specification of npt<0 will suppress the storage of spectator tracks.
403 // In the latter case abs(npt) will be taken as the number of participants.  
404 // In case of a standard Pythia run for 'elementary' particle interactions,
405 // the value of npt is totally irrelevant.
406 //
407 // The argument 'mlist' denotes the list mode used for Pylist().
408 // Note : mlist<0 suppresses the invokation of Pylist().
409 // By default, no listing is produced (i.e. mlist=-1).
410 //
411 // The argument 'medit' denotes the edit mode used for Pyedit().
412 // Note : medit<0 suppresses the invokation of Pyedit().
413 // By default, only 'stable' final particles are kept (i.e. medit=1). 
414 //
415 // In the case of a standard Pythia run concerning 'elementary' particle
416 // interactions, the projectile and target particle ID's for the created
417 // event structure are set to the corresponding Pythia KF codes.
418 // All the A and Z values are in that case set to zero.
419 // In case of a nucleus-nucleus interaction, the proper A and Z values for 
420 // the projectile and target particles are set in the event structure.
421 // However, in this case both particle ID's are set to zero.
422 //
423 // Note : Only in case an event passed the selection criteria as specified
424 //        via SelectEvent(), the event will appear on the output file.
425
426  fEventnum++; 
427
428  Int_t specmode=1;
429  if (npt<0)
430  {
431   specmode=0;
432   npt=abs(npt);
433  }
434
435  // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
436  Int_t ncols[4]={0,0,0,0};
437
438  Int_t zp=0;
439  Int_t ap=0;
440  Int_t zt=0;
441  Int_t at=0;
442
443  Int_t ncol=1;
444  if (fNucl)
445  {
446   if (npt<1 || npt>(fAproj+fAtarg))
447   {
448    cout << " *AliCollider::MakeEvent* Invalid input value. npt = " << npt
449         << " Aproj = " << fAproj << " Atarg = " << fAtarg << endl;
450    return;
451   }
452
453   // Determine the number of nucleon-nucleon collisions
454   ncol=npt/2;
455   if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
456
457   // Determine the number of the various types of N+N interactions
458   zp=fZproj;
459   ap=fAproj;
460   zt=fZtarg;
461   at=fAtarg;
462   Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
463   if (ap>at) maxa=1;
464   Float_t* rans=new Float_t[ncol];
465   fRan.Uniform(rans,ncol);
466   Float_t rndm=0;
467   for (Int_t i=0; i<ncol; i++)
468   {
469    GetFractions(zp,ap,zt,at);
470    rndm=rans[i];
471    if (rndm<=fFracpp) // p+p interaction
472    {
473     ncols[0]++;
474     if (maxa==2)
475     {
476      at--;
477      zt--;
478     } 
479     else
480     {
481      ap--;
482      zp--;
483     }
484    }
485    if (rndm>fFracpp && rndm<=(fFracpp+fFracnp)) // n+p interaction
486    {
487     ncols[1]++;
488     if (maxa==2)
489     {
490      at--;
491      zt--;
492     } 
493     else
494     {
495      ap--;
496     }
497    }
498    if (rndm>(fFracpp+fFracnp) && rndm<=(fFracpp+fFracnp+fFracpn)) // p+n interaction
499    {
500     ncols[2]++;
501     if (maxa==2)
502     {
503      at--;
504     } 
505     else
506     {
507      ap--;
508      zp--;
509     }
510    }
511    if (rndm>(fFracpp+fFracnp+fFracpn)) // n+n interaction
512    {
513     ncols[3]++; 
514     if (maxa==2)
515     {
516      at--;
517     } 
518     else
519     {
520      ap--;
521     }
522    }
523   }
524   delete [] rans;
525  }
526
527  if (!(fEventnum%fPrintfreq))
528  {
529   cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
530        << endl;
531   if (fNucl)
532   {
533    cout << " npart = " << npt << " ncol = " << ncol 
534         << " ncolpp = " << ncols[0] << " ncolnp = " << ncols[1]
535         << " ncolpn = " << ncols[2] << " ncolnn = " << ncols[3] << endl;
536   }
537  }
538
539  if (!fEvent)
540  {
541   fEvent=new AliEvent();
542   fEvent->SetOwner();
543   fEvent->SetName(GetName());
544   fEvent->SetTitle(GetTitle());
545  }
546
547  fEvent->Reset();
548  fEvent->SetRunNumber(fRunnum);
549  fEvent->SetEventNumber(fEventnum);
550
551  AliTrack t;
552  Ali3Vector p;
553  AliPosition r,rx;
554  Float_t v[3];
555  AliVertex vert;
556  Ali3Vector pproj,ptarg;
557
558  if (fVertexmode)
559  {
560   // Make sure the primary vertex gets correct location and Id=1
561   v[0]=0;
562   v[1]=0;
563   v[2]=0;
564   r.SetPosition(v,"car");
565   v[0]=fResolution;
566   v[1]=fResolution;
567   v[2]=fResolution;
568   r.SetPositionErrors(v,"car");
569
570   vert.SetId(1);
571   vert.SetTrackCopy(0);
572   vert.SetVertexCopy(0);
573   vert.SetPosition(r);
574   fEvent->AddVertex(vert,0);
575  }
576
577  Int_t kf=0;
578  Float_t charge=0,mass=0;
579  TString name;
580
581  Int_t ntypes=4;
582
583  // Singular settings for a normal Pythia elementary particle interation 
584  if (!fNucl)
585  {
586   ntypes=1;
587   ncols[0]=1;
588  }
589
590  // Generate all the various collisions
591  fSelect=0;      // Flag to indicate whether the total event is selected or not 
592  Int_t select=0; // Flag to indicate whether the sub-event is selected or not 
593  Int_t first=1;  // Flag to indicate the first collision process
594  Double_t pnucl;
595  Int_t npart=0,ntk=0;
596  Double_t dist=0;
597  for (Int_t itype=0; itype<ntypes; itype++)
598  {
599   if (fNucl)
600   {
601    if (itype==0 && ncols[itype]) Initialize(fFrame,"p","p",fWin);
602    if (itype==1 && ncols[itype]) Initialize(fFrame,"n","p",fWin);
603    if (itype==2 && ncols[itype]) Initialize(fFrame,"p","n",fWin);
604    if (itype==3 && ncols[itype]) Initialize(fFrame,"n","n",fWin);
605   }
606   for (Int_t jcol=0; jcol<ncols[itype]; jcol++)
607   {
608    GenerateEvent();
609
610    select=IsSelected();
611    if (select) fSelect=1;
612
613    if (first) // Store projectile and target information in the event structure
614    {
615     if (fNucl)
616     {
617      v[0]=GetP(1,1);
618      v[1]=GetP(1,2);
619      v[2]=GetP(1,3);
620      pproj.SetVector(v,"car");
621      pnucl=pproj.GetNorm();
622      fEvent->SetProjectile(fAproj,fZproj,pnucl);
623      v[0]=GetP(2,1);
624      v[1]=GetP(2,2);
625      v[2]=GetP(2,3);
626      ptarg.SetVector(v,"car");
627      pnucl=ptarg.GetNorm();
628      fEvent->SetTarget(fAtarg,fZtarg,pnucl);
629     }
630     else
631     {
632      v[0]=GetP(1,1);
633      v[1]=GetP(1,2);
634      v[2]=GetP(1,3);
635      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
636      kf=GetK(1,2);
637      fEvent->SetProjectile(0,0,pnucl,kf);
638      v[0]=GetP(2,1);
639      v[1]=GetP(2,2);
640      v[2]=GetP(2,3);
641      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
642      kf=GetK(2,2);
643      fEvent->SetTarget(0,0,pnucl,kf);
644     }
645     first=0;
646    }
647
648    if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
649
650    if (mlist>=0 && select) Pylist(mlist);
651
652    npart=GetN();
653    for (Int_t jpart=1; jpart<=npart; jpart++)
654    {
655     kf=GetK(jpart,2);
656     charge=Pychge(kf)/3.;
657     mass=GetP(jpart,5);
658     name=GetPyname(kf);
659
660     // 3-momentum in GeV/c
661     v[0]=GetP(jpart,1);
662     v[1]=GetP(jpart,2);
663     v[2]=GetP(jpart,3);
664     p.SetVector(v,"car");
665
666     // Production location in cm.
667     v[0]=GetV(jpart,1)/10;
668     v[1]=GetV(jpart,2)/10;
669     v[2]=GetV(jpart,3)/10;
670     r.SetPosition(v,"car");
671
672     ntk++;
673
674     t.Reset();
675     t.SetId(ntk);
676     t.SetParticleCode(kf);
677     t.SetName(name.Data());
678     t.SetCharge(charge);
679     t.SetMass(mass);
680     t.Set3Momentum(p);
681     t.SetBeginPoint(r);
682
683     fEvent->AddTrack(t);
684
685     // Build the vertex structures if requested
686     if (fVertexmode)
687     {
688      // Check if track belongs within the resolution to an existing vertex
689      Int_t add=0;  
690      for (Int_t jv=1; jv<=fEvent->GetNvertices(); jv++)
691      {
692       AliVertex* vx=fEvent->GetVertex(jv);
693       if (vx)
694       {
695        rx=vx->GetPosition();
696        dist=rx.GetDistance(r);
697        if (dist < fResolution)
698        {
699         AliTrack* tx=fEvent->GetIdTrack(ntk);
700         if (tx)
701         { 
702          vx->AddTrack(tx);
703          add=1;
704         }
705        }
706       }
707       if (add) break; // No need to look further for vertex candidates
708      }
709
710      // If track was not close enough to an existing vertex
711      // a new secondary vertex is created      
712      if (!add && fVertexmode>1)
713      {
714       AliTrack* tx=fEvent->GetIdTrack(ntk);
715       if (tx)
716       {
717        v[0]=fResolution;
718        v[1]=fResolution;
719        v[2]=fResolution;
720        r.SetPositionErrors(v,"car");
721        vert.Reset();
722        vert.SetTrackCopy(0);
723        vert.SetVertexCopy(0);
724        vert.SetId((fEvent->GetNvertices())+1);
725        vert.SetPosition(r);
726        vert.AddTrack(tx);
727        fEvent->AddVertex(vert,0);
728       } 
729      }
730     }
731    } // End of loop over the produced particles for each collision
732   } // End of loop over number of collisions for each type
733  } // End of loop over collision types
734
735  // Link sec. vertices to primary if requested
736  // Note that also the connecting tracks are automatically created
737  if (fVertexmode>2)
738  {
739   AliVertex* vp=fEvent->GetIdVertex(1); // Primary vertex
740   if (vp)
741   {
742    for (Int_t i=2; i<=fEvent->GetNvertices(); i++)
743    {
744     AliVertex* vx=fEvent->GetVertex(i);
745     if (vx)
746     {
747      if (vx->GetId() != 1) vp->AddVertex(vx);
748     }
749    }
750   }
751  }
752
753  // Include the spectator tracks in the event structure.
754  if (fNucl && specmode)
755  {
756   v[0]=0;
757   v[1]=0;
758   v[2]=0;
759   r.SetPosition(v,"car");
760
761   zp=fZproj-(ncols[0]+ncols[2]);
762   if (zp<0) zp=0;
763   ap=fAproj-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
764   if (ap<0) ap=0;
765   zt=fZtarg-(ncols[0]+ncols[1]);
766   if (zt<0) zt=0;
767   at=fAtarg-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
768   if (at<0) at=0;
769
770   Int_t nspec=0;
771
772   if (pproj.GetNorm() > fSpecpmin)
773   {
774    kf=2212; // Projectile spectator protons
775    charge=Pychge(kf)/3.;
776    mass=GetPMAS(Pycomp(kf),1);
777    name=GetPyname(kf);
778    for (Int_t iprojp=1; iprojp<=zp; iprojp++)
779    {
780     nspec++;
781     t.Reset();
782     t.SetId(-nspec);
783     t.SetParticleCode(kf);
784     t.SetName(name.Data());
785     t.SetTitle("Projectile spectator proton");
786     t.SetCharge(charge);
787     t.SetMass(mass);
788     t.Set3Momentum(pproj);
789     t.SetBeginPoint(r);
790
791     fEvent->AddTrack(t);
792    }
793
794    kf=2112; // Projectile spectator neutrons
795    charge=Pychge(kf)/3.;
796    mass=GetPMAS(Pycomp(kf),1);
797    name=GetPyname(kf);
798    for (Int_t iprojn=1; iprojn<=(ap-zp); iprojn++)
799    {
800     nspec++;
801     t.Reset();
802     t.SetId(-nspec);
803     t.SetParticleCode(kf);
804     t.SetName(name.Data());
805     t.SetTitle("Projectile spectator neutron");
806     t.SetCharge(charge);
807     t.SetMass(mass);
808     t.Set3Momentum(pproj);
809     t.SetBeginPoint(r);
810
811     fEvent->AddTrack(t);
812    }
813   }
814
815   if (ptarg.GetNorm() > fSpecpmin)
816   {
817    kf=2212; // Target spectator protons
818    charge=Pychge(kf)/3.;
819    mass=GetPMAS(Pycomp(kf),1);
820    name=GetPyname(kf);
821    for (Int_t itargp=1; itargp<=zt; itargp++)
822    {
823     nspec++;
824     t.Reset();
825     t.SetId(-nspec);
826     t.SetParticleCode(kf);
827     t.SetName(name.Data());
828     t.SetTitle("Target spectator proton");
829     t.SetCharge(charge);
830     t.SetMass(mass);
831     t.Set3Momentum(ptarg);
832     t.SetBeginPoint(r);
833
834     fEvent->AddTrack(t);
835    }
836
837    kf=2112; // Target spectator neutrons
838    charge=Pychge(kf)/3.;
839    mass=GetPMAS(Pycomp(kf),1);
840    name=GetPyname(kf);
841    for (Int_t itargn=1; itargn<=(at-zt); itargn++)
842    {
843     nspec++;
844     t.Reset();
845     t.SetId(-nspec);
846     t.SetParticleCode(kf);
847     t.SetName(name.Data());
848     t.SetTitle("Target spectator neutron");
849     t.SetCharge(charge);
850     t.SetMass(mass);
851     t.Set3Momentum(ptarg);
852     t.SetBeginPoint(r);
853
854     fEvent->AddTrack(t);
855    }
856   }
857
858  // Link the spectator tracks to the primary vertex.
859  if (fVertexmode)
860  {
861   AliVertex* vp=fEvent->GetIdVertex(1);
862   if (vp)
863   {
864    for (Int_t ispec=1; ispec<=nspec; ispec++)
865    {
866     AliTrack* tx=fEvent->GetIdTrack(-ispec);
867     if (tx) vp->AddTrack(tx);
868    }
869   }
870  }
871 }
872
873  if (mlist && !(fEventnum%fPrintfreq)) cout << endl; // Create empty output line after the event
874
875  if (fOutTree && fSelect) fOutTree->Fill();
876 }
877 ///////////////////////////////////////////////////////////////////////////
878 AliEvent* AliCollider::GetEvent(Int_t select)
879 {
880 // Provide pointer to the generated event structure.
881 //
882 // select = 0 : Always return the pointer to the generated event.
883 //          1 : Only return the pointer to the generated event in case
884 //              the event passed the selection criteria as specified via
885 //              SelectEvent(). Otherwise the value 0 will be returned.
886 //
887 // By invoking GetEvent() the default of select=0 will be used.
888
889  if (!select || fSelect)
890  {
891   return fEvent;
892  }
893  else
894  {
895   return 0;
896  }
897 }
898 ///////////////////////////////////////////////////////////////////////////
899 void AliCollider::EndRun()
900 {
901 // Properly close the output file (if needed).
902  if (fOutFile)
903  {
904   fOutFile->Write();
905   fOutFile->Close();
906   cout << " *AliCollider::EndRun* Output file correctly closed." << endl;
907  }
908 }
909 ///////////////////////////////////////////////////////////////////////////
910 void AliCollider::SetStable(Int_t id,Int_t mode)
911 {
912 // Declare whether a particle must be regarded as stable or not.
913 // The parameter "id" indicates the Pythia KF particle code, which
914 // basically is the PDG particle identifier code.
915 // The parameter "mode" indicates the action to be taken.
916 //
917 // mode = 0 : Particle will be able to decay
918 //        1 : Particle will be regarded as stable.
919 //
920 // In case the user does NOT explicitly invoke this function, the standard
921 // Pythia settings for the decay tables are used.
922 //
923 // When this function is invoked without the "mode" argument, then the
924 // default of mode=1 will be used for the specified particle.
925 //
926 // Notes :
927 // -------
928 // 1) This function should be invoked after the initialisation call
929 //    to AliCollider::Init.
930 // 2) Due to the internals of Pythia, there is no need to specify particles
931 //    and their corresponding anti-particles separately as (un)stable.
932 //    Once a particle has been declared (un)stable, the corresponding 
933 //    anti-particle will be treated in the same way.
934
935  if (mode==0 || mode==1)
936  {
937   Int_t kc=Pycomp(id);
938   Int_t decay=1-mode;
939   if (kc>0)
940   {
941    SetMDCY(kc,1,decay);
942   }
943   else 
944   {
945    cout << " *AliCollider::SetStable* Unknown particle code. id = " << id << endl;
946   }
947  }
948  else
949  {
950   cout << " *AliCollider::SetStable* Invalid parameter. mode = " << mode << endl;
951  }
952 }
953 ///////////////////////////////////////////////////////////////////////////
954 void AliCollider::SelectEvent(Int_t id)
955 {
956 // Add a particle to the event selection list.
957 // The parameter "id" indicates the Pythia KF particle code, which
958 // basically is the PDG particle identifier code.
959 // In case the user has built a selection list via this procedure, only the
960 // events in which one of the particles specified in the list was generated
961 // will be kept. 
962 // The investigation of the generated particles takes place when the complete
963 // event is in memory, including all (shortlived) mother particles and resonances.
964 // So, the settings of the various particle decay modes have no influence on
965 // the event selection described here.
966 //
967 // If no list has been specified, all events will be accepted.  
968 //
969 // Note : id=0 will delete the selection list.
970 //
971 // Be aware of the fact that severe selection criteria (i.e. selecting only
972 // rare events) may result in long runtimes before an event sample has been
973 // obtained.
974 //
975  if (!id)
976  {
977   if (fSelections)
978   {
979    delete fSelections;
980    fSelections=0;
981   }
982  }
983  else
984  {
985   Int_t kc=Pycomp(id);
986   if (!fSelections)
987   {
988    fSelections=new TArrayI(1);
989    fSelections->AddAt(kc,0);
990   }
991   else
992   {
993    Int_t exist=0;
994    Int_t size=fSelections->GetSize();
995    for (Int_t i=0; i<size; i++)
996    {
997     if (kc==fSelections->At(i))
998     {
999      exist=1;
1000      break;
1001     }
1002    }
1003   
1004    if (!exist)
1005    {
1006     fSelections->Set(size+1);
1007     fSelections->AddAt(kc,size);
1008    }
1009   }
1010  }
1011 }
1012 ///////////////////////////////////////////////////////////////////////////
1013 Int_t AliCollider::GetSelectionFlag()
1014 {
1015 // Return the value of the selection flag for the total event.
1016 // When the event passed the selection criteria as specified via
1017 // SelectEvent() the value 1 is returned, otherwise the value 0 is returned.
1018  return fSelect;
1019 }
1020 ///////////////////////////////////////////////////////////////////////////
1021 Int_t AliCollider::IsSelected()
1022 {
1023 // Check whether the generated (sub)event contains one of the particles
1024 // specified in the selection list via SelectEvent().
1025 // If this is the case or when no selection list is present, the value 1
1026 // will be returned, indicating the event is selected to be kept.
1027 // Otherwise the value 0 will be returned.
1028
1029  if (!fSelections) return 1; 
1030
1031  Int_t nsel=fSelections->GetSize();
1032  Int_t npart=GetN();
1033  Int_t kf,kc;
1034
1035  Int_t select=0;
1036  for (Int_t jpart=1; jpart<=npart; jpart++)
1037  {
1038   kf=GetK(jpart,2);
1039   kc=Pycomp(kf);
1040   for (Int_t i=0; i<nsel; i++)
1041   {
1042    if (kc==fSelections->At(i))
1043    {
1044     select=1;
1045     break;
1046    }
1047   }
1048   if (select) break;
1049  }
1050  return select;
1051 }
1052 ///////////////////////////////////////////////////////////////////////////
1053 void AliCollider::SetSpectatorPmin(Float_t pmin)
1054 {
1055 // Set minimal momentum in GeV/c for spectator tracks to be stored.
1056 // Spectator tracks with a momentum below this threshold will not be stored
1057 // in the (output) event structure.
1058 // This facility allows to minimise the output file size.
1059 // Note that when the user wants to boost the event into another reference
1060 // frame these spectator tracks might have got momenta above the threshold.
1061 // However, when the spectator tracks were not stored in the event structure
1062 // in the original frame, there is no way to retreive them anymore. 
1063  fSpecpmin=pmin;
1064 }
1065 ///////////////////////////////////////////////////////////////////////////
1066 Float_t AliCollider::GetSpectatorPmin()
1067 {
1068 // Provide the minimal spectator momentum in GeV/c.
1069  return fSpecpmin;
1070 }
1071 ///////////////////////////////////////////////////////////////////////////
1072 TString AliCollider::GetPyname(Int_t kf)
1073 {
1074 // Provide the correctly truncated Pythia particle name for PGD code kf
1075 //
1076 // The TPythia6::Pyname returned name is copied into a TString and truncated
1077 // at the first blank to prevent funny trailing characters due to incorrect
1078 // stripping of empty characters in TPythia6::Pyname.
1079 // The truncation at the first blank is allowed due to the Pythia convention
1080 // that particle names never contain blanks.
1081  char name[16];
1082  TString sname;
1083  Pyname(kf,name);
1084  sname=name[0];
1085  for (Int_t i=1; i<16; i++)
1086  {
1087   if (name[i]==' ') break;
1088   sname=sname+name[i];
1089  }
1090  return sname;
1091 }
1092 ///////////////////////////////////////////////////////////////////////////