01-dev-2004 NvE Installation scripts modified to make use of environment variable...
[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.12 2004/05/04 15:33:04 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/05/04 15:33:04 $ 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 //
124 // Some Pythia default MC parameters are automatically modified to provide
125 // more suitable running conditions for soft processes in view of
126 // nucleus-nucleus interactions and astrophysical processes.
127 // The user may initialise the generator with all the default Pythia
128 // parameters and obtain full user control to modify the settings by means
129 // of the SetUserControl memberfunction.
130 //
131 // Refer to the SetElastic memberfunction for the inclusion of elastic
132 // and diffractive processes.
133 // By default these processes are not included.
134
135  fVertexmode=0;    // No vertex structure creation
136  fResolution=1e-5; // Standard resolution is 0.1 micron
137  fRunnum=0;
138  fEventnum=0;
139  fPrintfreq=1;
140  fUserctrl=0; // Automatic optimisation of some MC parameters 
141  fElastic=0;  // No elastic and diffractive processes
142
143  fEvent=0;
144
145  fSpecpmin=0;
146
147  fFrame="none";
148  fWin=0;
149
150  fNucl=0;
151  fZproj=0;
152  fAproj=0;
153  fZtarg=0;
154  fAtarg=0;
155  fFracpp=0;
156  fFracnp=0;
157  fFracpn=0;
158  fFracnn=0;
159
160  fOutFile=0;
161  fOutTree=0;
162
163  fSelections=0;
164  fSelect=0;
165
166  TString s=GetName();
167  s+=" (AliCollider)";
168  SetName(s.Data());
169 }
170 ///////////////////////////////////////////////////////////////////////////
171 AliCollider::~AliCollider()
172 {
173 // Default destructor
174  if (fEvent)
175  {
176   delete fEvent;
177   fEvent=0;
178  }
179  if (fOutFile)
180  {
181   delete fOutFile;
182   fOutFile=0;
183  }
184  if (fOutTree)
185  {
186   delete fOutTree;
187   fOutTree=0;
188  }
189  if (fSelections)
190  {
191   delete fSelections;
192   fSelections=0;
193  }
194 }
195 ///////////////////////////////////////////////////////////////////////////
196 void AliCollider::SetOutputFile(TString s)
197 {
198 // Create the output file containing all the data in ROOT output format.
199  if (fOutFile)
200  {
201   delete fOutFile;
202   fOutFile=0;
203  }
204  fOutFile=new TFile(s.Data(),"RECREATE","AliCollider data");
205
206  if (fOutTree)
207  {
208   delete fOutTree;
209   fOutTree=0;
210  }
211  fOutTree=new TTree("T","AliCollider event data");
212
213  Int_t bsize=32000;
214  Int_t split=0;
215  fOutTree->Branch("Events","AliEvent",&fEvent,bsize,split);
216 }
217 ///////////////////////////////////////////////////////////////////////////
218 void AliCollider::SetVertexMode(Int_t mode)
219 {
220 // Set the mode of the vertex structure creation.
221 //
222 // By default all generated tracks will only appear in the AliEvent
223 // structure without any primary (and secondary) vertex structure.
224 // The user can build the vertex structure if he/she wants by means
225 // of the beginpoint location of each AliTrack.
226 //
227 // However, one can also let AliCollider automatically create
228 // the primary (and secondary) vertex structure(s).
229 // In this case the primary vertex is given Id=1 and all sec. vertices
230 // are given Id's 2,3,4,....
231 // All vertices are created as standalone entities in the AliEvent structure
232 // without any linking between the various vertices.
233 // For this automated process, the user-selected resolution
234 // (see SetResolution) is used to decide whether or not certain vertex
235 // locations can be resolved.
236 // In case no vertex creation is selected (i.e. the default mode=0),
237 // the value of the resolution is totally irrelevant.
238 //
239 // The user can also let AliCollider automatically connect the sec. vertices
240 // to the primary vertex (i.e. mode=3). This process will also automatically
241 // generate the tracks connecting the vertices.
242 // Note that the result of the mode=3 operation may be very sensitive to
243 // the resolution parameter. Therefore, no attempt is made to distinguish
244 // between secondary, tertiary etc... vertices. All sec. vertices are
245 // linked to the primary one.
246 //  
247 // Irrespective of the selected mode, all generated tracks can be obtained
248 // directly from the AliEvent structure.
249 // In case (sec.) vertex creation is selected, all generated vertices can
250 // also be obtained directly from the AliEvent structure. 
251 // These (sec.) vertices contain only the corresponding pointers to the various
252 // tracks which are stored in the AliEvent structure.
253 //
254 // Overview of vertex creation modes :
255 // -----------------------------------
256 // mode = 0 ==> No vertex structure will be created
257 //        1 ==> Only primary vertex structure will be created
258 //        2 ==> Unconnected primary and secondary vertices will be created
259 //        3 ==> Primary and secondary vertices will be created where all the
260 //              sec. vertices will be connected to the primary vertex.
261 //              Also the vertex connecting tracks will be automatically
262 //              generated. 
263 //
264  if (mode<0 || mode >3)
265  {
266   cout << " *AliCollider::SetVertexMode* Invalid argument mode : " << mode << endl;
267   fVertexmode=0;
268  }
269  else
270  {
271   fVertexmode=mode;
272  }
273 }
274 ///////////////////////////////////////////////////////////////////////////
275 Int_t AliCollider::GetVertexMode() const
276 {
277 // Provide the current mode for vertex structure creation.
278  return fVertexmode;
279 }
280 ///////////////////////////////////////////////////////////////////////////
281 void AliCollider::SetResolution(Double_t res)
282 {
283 // Set the resolution (in cm) for resolving (sec.) vertices.
284 // By default this resolution is set to 0.1 micron.
285 // Note : In case no vertex creation has been selected, the value of
286 //        the resolution is totally irrelevant.
287  fResolution=fabs(res);
288 }
289 ///////////////////////////////////////////////////////////////////////////
290 Double_t AliCollider::GetResolution() const
291 {
292 // Provide the current resolution (in cm) for resolving (sec.) vertices.
293  return fResolution;
294 }
295 ///////////////////////////////////////////////////////////////////////////
296 void AliCollider::SetRunNumber(Int_t run)
297 {
298 // Set the user defined run number.
299 // By default the run number is set to 0.
300  fRunnum=run;
301 }
302 ///////////////////////////////////////////////////////////////////////////
303 Int_t AliCollider::GetRunNumber() const
304 {
305 // Provide the user defined run number.
306  return fRunnum;
307 }
308 ///////////////////////////////////////////////////////////////////////////
309 void AliCollider::SetPrintFreq(Int_t n)
310 {
311 // Set the print frequency for every 'n' events.
312 // By default the printfrequency is set to 1 (i.e. every event).
313  fPrintfreq=n;
314 }
315 ///////////////////////////////////////////////////////////////////////////
316 Int_t AliCollider::GetPrintFreq() const
317 {
318 // Provide the user selected print frequency.
319  return fPrintfreq;
320 }
321 ///////////////////////////////////////////////////////////////////////////
322 void AliCollider::SetUserControl(Int_t flag)
323 {
324 // Set the user control flag w.r.t. disabling automatic optimisation
325 // of some Pythia default MC parameters for soft interactions in view of
326 // nucleus-nucleus collisions and astrophysical processes.
327 // Flag = 0 : Limited user control (automatic optimisation enabled)
328 //        1 : Full user control (automatic optimisation disabled)
329 // By default the user control is set to 0 (i.e. automatic optimisation).
330 // See the Init() memberfunctions for further details w.r.t. the optimisations.
331  fUserctrl=flag;
332 }
333 ///////////////////////////////////////////////////////////////////////////
334 Int_t AliCollider::GetUserControl() const
335 {
336 // Provide the value of the user control flag.
337  return fUserctrl;
338 }
339 ///////////////////////////////////////////////////////////////////////////
340 void AliCollider::SetElastic(Int_t flag)
341 {
342 // Set the flag w.r.t. inclusion of elastic and diffractive processes.
343 // By default these processes are not included.
344 // Flag = 0 : Do not include elastic and diffractive processes
345 //        1 : Elastic and diffractive processes will be included
346  fElastic=flag;
347 }
348 ///////////////////////////////////////////////////////////////////////////
349 Int_t AliCollider::GetElastic() const
350 {
351 // Provide the value of the control flag for elastic and diffractive processes.
352  return fElastic;
353 }
354 ///////////////////////////////////////////////////////////////////////////
355 void AliCollider::Init(char* frame,char* beam,char* target,Float_t win)
356 {
357 // Initialisation of the underlying Pythia generator package.
358 // The event number is reset to 0.
359 // This routine just invokes TPythia6::Initialize(...) and the arguments
360 // have the corresponding meaning.
361 // Some Pythia default MC parameters are automatically modified to provide
362 // more suitable running conditions for soft processes in view of
363 // astrophysical processes.
364 // The optimisations consist of : 
365 // * Usage of real photons for photon beams of targets
366 // * Minimum CMS energy of 3 GeV for the event
367 // * Activation of the default K factor values
368 //   with separate settings for ordinary and color annihilation graphs.
369 // The user may initialise the generator with all the default Pythia
370 // parameters and obtain full user control to modify the settings by means
371 // of invoking the SetUserControl memberfunction before this initialisation.
372 // Note that the inclusion of elastic and diffractive processes is controlled
373 // by invokation of the SetElastic memberfunction before this initialisation,
374 // irrespective of the UserControl selection.
375
376  if (!fUserctrl) // Optimisation of some MC parameters
377  {
378   SetMSTP(14,10); // Real photons for photon beams or targets
379   SetPARP(2,3.);  // Minimum CMS energy for the event
380   SetMSTP(33,2);  // Activate K factor. Separate for ordinary and color annih. graphs
381  }
382
383  if (fElastic) SetMSEL(2); // Include low-Pt, elastic and diffractive events
384
385  fEventnum=0;
386  fNucl=0;
387  fFrame=frame;
388  fWin=win;
389  Initialize(frame,beam,target,win);
390
391  cout << endl;
392  cout << " *AliCollider::Init* Standard Pythia initialisation." << endl;
393  cout << " Beam particle : " << beam << " Target particle : " << target
394       << " Frame = " << frame << " Energy = " << win
395       << endl;
396 }
397 ///////////////////////////////////////////////////////////////////////////
398 void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win)
399 {
400 // Initialisation of the underlying Pythia generator package for the generation
401 // of nucleus-nucleus interactions.
402 // The event number is reset to 0.
403 // In addition to the Pythia standard arguments 'frame' and 'win', the user
404 // can specify here (Z,A) values of the projectile and target nuclei.
405 //
406 // Note : The 'win' value denotes either the cms energy per nucleon-nucleon collision
407 //        (i.e. frame="cms") or the momentum per nucleon in all other cases.
408 //
409 // Some Pythia default MC parameters are automatically modified to provide
410 // more suitable running conditions for soft processes in view of
411 // nucleus-nucleus interactions and astrophysical processes.
412 // The optimisations consist of : 
413 // * Minimum CMS energy of 3 GeV for the event
414 // * Activation of the default K factor values
415 //   with separate settings for ordinary and color annihilation graphs.
416 // The user may initialise the generator with all the default Pythia
417 // parameters and obtain full user control to modify the settings by means
418 // of invoking the SetUserControl memberfunction before this initialisation.
419 // Note that the inclusion of elastic and diffractive processes is controlled
420 // by invokation of the SetElastic memberfunction before this initialisation,
421 // irrespective of the UserControl selection.
422
423  if (!fUserctrl) // Optimisation of some MC parameters
424  {
425   SetPARP(2,3.);  // Minimum CMS energy for the event
426   SetMSTP(33,2);  // Activate K factor. Separate for ordinary and color annih. graphs
427  }
428
429  if (fElastic) SetMSEL(2); // Include low-Pt, elastic and diffractive events
430
431  fEventnum=0;
432  fNucl=1;
433  fFrame=frame;
434  fWin=win;
435  fZproj=0;
436  fAproj=0;
437  fZtarg=0;
438  fAtarg=0;
439  fFracpp=0;
440  fFracnp=0;
441  fFracpn=0;
442  fFracnn=0;
443
444  if (ap<1 || at<1 || zp>ap || zt>at)
445  {
446   cout << endl;
447   cout << " *AliCollider::Init* Invalid input value(s). Zproj = " << zp
448        << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
449   return;
450  }
451
452  fZproj=zp;
453  fAproj=ap;
454  fZtarg=zt;
455  fAtarg=at;
456
457  cout << endl;
458  cout << " *AliCollider::Init* Nucleus-Nucleus generator initialisation." << endl;
459  cout << " Zproj = " << zp << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at
460       << " Frame = " << frame << " Energy = " << win
461       << endl;
462 }
463 ///////////////////////////////////////////////////////////////////////////
464 void AliCollider::GetFractions(Float_t zp,Float_t ap,Float_t zt,Float_t at)
465 {
466 // Determine the fractions for the various N-N collision processes.
467 // The various processes are : p+p, n+p, p+n and n+n.
468  if (zp<0) zp=0;
469  if (zt<0) zt=0;
470
471  fFracpp=0;
472  fFracnp=0;
473  fFracpn=0;
474  fFracnn=0;
475
476  if (ap>0 && at>0)
477  {
478   fFracpp=(zp/ap)*(zt/at);
479   fFracnp=(1.-zp/ap)*(zt/at);
480   fFracpn=(zp/ap)*(1.-zt/at);
481   fFracnn=(1.-zp/ap)*(1.-zt/at);
482  }
483 }
484 ///////////////////////////////////////////////////////////////////////////
485 void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
486 {
487 // Generate one event.
488 // In case of a nucleus-nucleus interaction, the argument 'npt' denotes
489 // the number of participant nucleons.
490 // Normally also the spectator tracks will be stored into the event structure.
491 // The spectator tracks have a negative user Id to distinguish them from the
492 // ordinary generated tracks.
493 // In case the user has selected the creation of vertex structures, the spectator
494 // tracks will be linked to the primary vertex.
495 // However, specification of npt<0 will suppress the storage of spectator tracks.
496 // In the latter case abs(npt) will be taken as the number of participants.  
497 // In case of a standard Pythia run for 'elementary' particle interactions,
498 // the value of npt is totally irrelevant.
499 //
500 // The argument 'mlist' denotes the list mode used for Pylist().
501 // Note : mlist<0 suppresses the invokation of Pylist().
502 // By default, no listing is produced (i.e. mlist=-1).
503 //
504 // The argument 'medit' denotes the edit mode used for Pyedit().
505 // Note : medit<0 suppresses the invokation of Pyedit().
506 // By default, only 'stable' final particles are kept (i.e. medit=1). 
507 //
508 // In the case of a standard Pythia run concerning 'elementary' particle
509 // interactions, the projectile and target particle ID's for the created
510 // event structure are set to the corresponding Pythia KF codes.
511 // All the A and Z values are in that case set to zero.
512 // In case of a nucleus-nucleus interaction, the proper A and Z values for 
513 // the projectile and target particles are set in the event structure.
514 // However, in this case both particle ID's are set to zero.
515 //
516 // Note : Only in case an event passed the selection criteria as specified
517 //        via SelectEvent(), the event will appear on the output file.
518
519  fEventnum++; 
520
521  Int_t specmode=1;
522  if (npt<0)
523  {
524   specmode=0;
525   npt=abs(npt);
526  }
527
528  // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
529  Int_t ncols[4]={0,0,0,0};
530
531  Int_t zp=0;
532  Int_t ap=0;
533  Int_t zt=0;
534  Int_t at=0;
535
536  Int_t ncol=1;
537  if (fNucl)
538  {
539   if (npt<1 || npt>(fAproj+fAtarg))
540   {
541    cout << " *AliCollider::MakeEvent* Invalid input value. npt = " << npt
542         << " Aproj = " << fAproj << " Atarg = " << fAtarg << endl;
543    return;
544   }
545
546   // Determine the number of nucleon-nucleon collisions
547   ncol=npt/2;
548   if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
549
550   // Determine the number of the various types of N+N interactions
551   zp=fZproj;
552   ap=fAproj;
553   zt=fZtarg;
554   at=fAtarg;
555   Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
556   if (ap>at) maxa=1;
557   Float_t* rans=new Float_t[ncol];
558   fRan.Uniform(rans,ncol);
559   Float_t rndm=0;
560   for (Int_t i=0; i<ncol; i++)
561   {
562    GetFractions(zp,ap,zt,at);
563    rndm=rans[i];
564    if (rndm<=fFracpp) // p+p interaction
565    {
566     ncols[0]++;
567     if (maxa==2)
568     {
569      at--;
570      zt--;
571     } 
572     else
573     {
574      ap--;
575      zp--;
576     }
577    }
578    if (rndm>fFracpp && rndm<=(fFracpp+fFracnp)) // n+p interaction
579    {
580     ncols[1]++;
581     if (maxa==2)
582     {
583      at--;
584      zt--;
585     } 
586     else
587     {
588      ap--;
589     }
590    }
591    if (rndm>(fFracpp+fFracnp) && rndm<=(fFracpp+fFracnp+fFracpn)) // p+n interaction
592    {
593     ncols[2]++;
594     if (maxa==2)
595     {
596      at--;
597     } 
598     else
599     {
600      ap--;
601      zp--;
602     }
603    }
604    if (rndm>(fFracpp+fFracnp+fFracpn)) // n+n interaction
605    {
606     ncols[3]++; 
607     if (maxa==2)
608     {
609      at--;
610     } 
611     else
612     {
613      ap--;
614     }
615    }
616   }
617   delete [] rans;
618  }
619
620  if (!(fEventnum%fPrintfreq))
621  {
622   cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
623        << endl;
624   if (fNucl)
625   {
626    cout << " npart = " << npt << " ncol = " << ncol 
627         << " ncolpp = " << ncols[0] << " ncolnp = " << ncols[1]
628         << " ncolpn = " << ncols[2] << " ncolnn = " << ncols[3] << endl;
629   }
630  }
631
632  if (!fEvent)
633  {
634   fEvent=new AliEvent();
635   fEvent->SetOwner();
636   fEvent->SetName(GetName());
637   fEvent->SetTitle(GetTitle());
638  }
639
640  fEvent->Reset();
641  fEvent->SetRunNumber(fRunnum);
642  fEvent->SetEventNumber(fEventnum);
643
644  AliTrack t;
645  Ali3Vector p;
646  AliPosition r,rx;
647  Float_t v[3];
648  AliVertex vert;
649  Ali3Vector pproj,ptarg;
650
651  if (fVertexmode)
652  {
653   // Make sure the primary vertex gets correct location and Id=1
654   v[0]=0;
655   v[1]=0;
656   v[2]=0;
657   r.SetPosition(v,"car");
658   v[0]=fResolution;
659   v[1]=fResolution;
660   v[2]=fResolution;
661   r.SetPositionErrors(v,"car");
662
663   vert.SetId(1);
664   vert.SetTrackCopy(0);
665   vert.SetVertexCopy(0);
666   vert.SetPosition(r);
667   fEvent->AddVertex(vert,0);
668  }
669
670  Int_t kf=0;
671  Float_t charge=0,mass=0;
672  TString name;
673
674  Int_t ntypes=4;
675
676  // Singular settings for a normal Pythia elementary particle interation 
677  if (!fNucl)
678  {
679   ntypes=1;
680   ncols[0]=1;
681  }
682
683  // Generate all the various collisions
684  fSelect=0;      // Flag to indicate whether the total event is selected or not 
685  Int_t select=0; // Flag to indicate whether the sub-event is selected or not 
686  Int_t first=1;  // Flag to indicate the first collision process
687  Double_t pnucl;
688  Int_t npart=0,ntk=0;
689  Double_t dist=0;
690  for (Int_t itype=0; itype<ntypes; itype++)
691  {
692   if (fNucl)
693   {
694    if (itype==0 && ncols[itype]) Initialize(fFrame,"p","p",fWin);
695    if (itype==1 && ncols[itype]) Initialize(fFrame,"n","p",fWin);
696    if (itype==2 && ncols[itype]) Initialize(fFrame,"p","n",fWin);
697    if (itype==3 && ncols[itype]) Initialize(fFrame,"n","n",fWin);
698   }
699   for (Int_t jcol=0; jcol<ncols[itype]; jcol++)
700   {
701    GenerateEvent();
702
703    select=IsSelected();
704    if (select) fSelect=1;
705
706    if (first) // Store projectile and target information in the event structure
707    {
708     if (fNucl)
709     {
710      v[0]=GetP(1,1);
711      v[1]=GetP(1,2);
712      v[2]=GetP(1,3);
713      pproj.SetVector(v,"car");
714      pnucl=pproj.GetNorm();
715      fEvent->SetProjectile(fAproj,fZproj,pnucl);
716      v[0]=GetP(2,1);
717      v[1]=GetP(2,2);
718      v[2]=GetP(2,3);
719      ptarg.SetVector(v,"car");
720      pnucl=ptarg.GetNorm();
721      fEvent->SetTarget(fAtarg,fZtarg,pnucl);
722     }
723     else
724     {
725      v[0]=GetP(1,1);
726      v[1]=GetP(1,2);
727      v[2]=GetP(1,3);
728      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
729      kf=GetK(1,2);
730      fEvent->SetProjectile(0,0,pnucl,kf);
731      v[0]=GetP(2,1);
732      v[1]=GetP(2,2);
733      v[2]=GetP(2,3);
734      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
735      kf=GetK(2,2);
736      fEvent->SetTarget(0,0,pnucl,kf);
737     }
738     first=0;
739    }
740
741    if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
742
743    if (mlist>=0 && select) Pylist(mlist);
744
745    npart=GetN();
746    for (Int_t jpart=1; jpart<=npart; jpart++)
747    {
748     kf=GetK(jpart,2);
749     charge=Pychge(kf)/3.;
750     mass=GetP(jpart,5);
751     name=GetPyname(kf);
752
753     // 3-momentum in GeV/c
754     v[0]=GetP(jpart,1);
755     v[1]=GetP(jpart,2);
756     v[2]=GetP(jpart,3);
757     p.SetVector(v,"car");
758
759     // Production location in cm.
760     v[0]=GetV(jpart,1)/10;
761     v[1]=GetV(jpart,2)/10;
762     v[2]=GetV(jpart,3)/10;
763     r.SetPosition(v,"car");
764
765     ntk++;
766
767     t.Reset();
768     t.SetId(ntk);
769     t.SetParticleCode(kf);
770     t.SetName(name.Data());
771     t.SetCharge(charge);
772     t.SetMass(mass);
773     t.Set3Momentum(p);
774     t.SetBeginPoint(r);
775
776     fEvent->AddTrack(t);
777
778     // Build the vertex structures if requested
779     if (fVertexmode)
780     {
781      // Check if track belongs within the resolution to an existing vertex
782      Int_t add=0;  
783      for (Int_t jv=1; jv<=fEvent->GetNvertices(); jv++)
784      {
785       AliVertex* vx=fEvent->GetVertex(jv);
786       if (vx)
787       {
788        rx=vx->GetPosition();
789        dist=rx.GetDistance(r);
790        if (dist < fResolution)
791        {
792         AliTrack* tx=fEvent->GetIdTrack(ntk);
793         if (tx)
794         { 
795          vx->AddTrack(tx);
796          add=1;
797         }
798        }
799       }
800       if (add) break; // No need to look further for vertex candidates
801      }
802
803      // If track was not close enough to an existing vertex
804      // a new secondary vertex is created      
805      if (!add && fVertexmode>1)
806      {
807       AliTrack* tx=fEvent->GetIdTrack(ntk);
808       if (tx)
809       {
810        v[0]=fResolution;
811        v[1]=fResolution;
812        v[2]=fResolution;
813        r.SetPositionErrors(v,"car");
814        vert.Reset();
815        vert.SetTrackCopy(0);
816        vert.SetVertexCopy(0);
817        vert.SetId((fEvent->GetNvertices())+1);
818        vert.SetPosition(r);
819        vert.AddTrack(tx);
820        fEvent->AddVertex(vert,0);
821       } 
822      }
823     }
824    } // End of loop over the produced particles for each collision
825   } // End of loop over number of collisions for each type
826  } // End of loop over collision types
827
828  // Link sec. vertices to primary if requested
829  // Note that also the connecting tracks are automatically created
830  if (fVertexmode>2)
831  {
832   AliVertex* vp=fEvent->GetIdVertex(1); // Primary vertex
833   if (vp)
834   {
835    for (Int_t i=2; i<=fEvent->GetNvertices(); i++)
836    {
837     AliVertex* vx=fEvent->GetVertex(i);
838     if (vx)
839     {
840      if (vx->GetId() != 1) vp->AddVertex(vx);
841     }
842    }
843   }
844  }
845
846  // Include the spectator tracks in the event structure.
847  if (fNucl && specmode)
848  {
849   v[0]=0;
850   v[1]=0;
851   v[2]=0;
852   r.SetPosition(v,"car");
853
854   zp=fZproj-(ncols[0]+ncols[2]);
855   if (zp<0) zp=0;
856   ap=fAproj-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
857   if (ap<0) ap=0;
858   zt=fZtarg-(ncols[0]+ncols[1]);
859   if (zt<0) zt=0;
860   at=fAtarg-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
861   if (at<0) at=0;
862
863   Int_t nspec=0;
864
865   if (pproj.GetNorm() > fSpecpmin)
866   {
867    kf=2212; // Projectile spectator protons
868    charge=Pychge(kf)/3.;
869    mass=GetPMAS(Pycomp(kf),1);
870    name=GetPyname(kf);
871    for (Int_t iprojp=1; iprojp<=zp; iprojp++)
872    {
873     nspec++;
874     t.Reset();
875     t.SetId(-nspec);
876     t.SetParticleCode(kf);
877     t.SetName(name.Data());
878     t.SetTitle("Projectile spectator proton");
879     t.SetCharge(charge);
880     t.SetMass(mass);
881     t.Set3Momentum(pproj);
882     t.SetBeginPoint(r);
883
884     fEvent->AddTrack(t);
885    }
886
887    kf=2112; // Projectile spectator neutrons
888    charge=Pychge(kf)/3.;
889    mass=GetPMAS(Pycomp(kf),1);
890    name=GetPyname(kf);
891    for (Int_t iprojn=1; iprojn<=(ap-zp); iprojn++)
892    {
893     nspec++;
894     t.Reset();
895     t.SetId(-nspec);
896     t.SetParticleCode(kf);
897     t.SetName(name.Data());
898     t.SetTitle("Projectile spectator neutron");
899     t.SetCharge(charge);
900     t.SetMass(mass);
901     t.Set3Momentum(pproj);
902     t.SetBeginPoint(r);
903
904     fEvent->AddTrack(t);
905    }
906   }
907
908   if (ptarg.GetNorm() > fSpecpmin)
909   {
910    kf=2212; // Target spectator protons
911    charge=Pychge(kf)/3.;
912    mass=GetPMAS(Pycomp(kf),1);
913    name=GetPyname(kf);
914    for (Int_t itargp=1; itargp<=zt; itargp++)
915    {
916     nspec++;
917     t.Reset();
918     t.SetId(-nspec);
919     t.SetParticleCode(kf);
920     t.SetName(name.Data());
921     t.SetTitle("Target spectator proton");
922     t.SetCharge(charge);
923     t.SetMass(mass);
924     t.Set3Momentum(ptarg);
925     t.SetBeginPoint(r);
926
927     fEvent->AddTrack(t);
928    }
929
930    kf=2112; // Target spectator neutrons
931    charge=Pychge(kf)/3.;
932    mass=GetPMAS(Pycomp(kf),1);
933    name=GetPyname(kf);
934    for (Int_t itargn=1; itargn<=(at-zt); itargn++)
935    {
936     nspec++;
937     t.Reset();
938     t.SetId(-nspec);
939     t.SetParticleCode(kf);
940     t.SetName(name.Data());
941     t.SetTitle("Target spectator neutron");
942     t.SetCharge(charge);
943     t.SetMass(mass);
944     t.Set3Momentum(ptarg);
945     t.SetBeginPoint(r);
946
947     fEvent->AddTrack(t);
948    }
949   }
950
951  // Link the spectator tracks to the primary vertex.
952  if (fVertexmode)
953  {
954   AliVertex* vp=fEvent->GetIdVertex(1);
955   if (vp)
956   {
957    for (Int_t ispec=1; ispec<=nspec; ispec++)
958    {
959     AliTrack* tx=fEvent->GetIdTrack(-ispec);
960     if (tx) vp->AddTrack(tx);
961    }
962   }
963  }
964 }
965
966  if (!(fEventnum%fPrintfreq) && (mlist || fEvent))
967  {
968   if (fEvent)
969   {
970    cout << " Number of tracks in the event structure : "
971         << fEvent->GetNtracks() << endl;
972   }
973   cout << endl; // Create empty output line after the event
974  }
975
976  if (fOutTree && fSelect) fOutTree->Fill();
977 }
978 ///////////////////////////////////////////////////////////////////////////
979 AliEvent* AliCollider::GetEvent(Int_t select) const
980 {
981 // Provide pointer to the generated event structure.
982 //
983 // select = 0 : Always return the pointer to the generated event.
984 //          1 : Only return the pointer to the generated event in case
985 //              the event passed the selection criteria as specified via
986 //              SelectEvent(). Otherwise the value 0 will be returned.
987 //
988 // By invoking GetEvent() the default of select=0 will be used.
989
990  if (!select || fSelect)
991  {
992   return fEvent;
993  }
994  else
995  {
996   return 0;
997  }
998 }
999 ///////////////////////////////////////////////////////////////////////////
1000 void AliCollider::EndRun()
1001 {
1002 // Properly close the output file (if needed).
1003  if (fOutFile)
1004  {
1005   fOutFile->Write();
1006   fOutFile->Close();
1007   cout << " *AliCollider::EndRun* Output file correctly closed." << endl;
1008  }
1009 }
1010 ///////////////////////////////////////////////////////////////////////////
1011 void AliCollider::SetStable(Int_t id,Int_t mode)
1012 {
1013 // Declare whether a particle must be regarded as stable or not.
1014 // The parameter "id" indicates the Pythia KF particle code, which
1015 // basically is the PDG particle identifier code.
1016 // The parameter "mode" indicates the action to be taken.
1017 //
1018 // mode = 0 : Particle will be able to decay
1019 //        1 : Particle will be regarded as stable.
1020 //
1021 // In case the user does NOT explicitly invoke this function, the standard
1022 // Pythia settings for the decay tables are used.
1023 //
1024 // When this function is invoked without the "mode" argument, then the
1025 // default of mode=1 will be used for the specified particle.
1026 //
1027 // Notes :
1028 // -------
1029 // 1) This function should be invoked after the initialisation call
1030 //    to AliCollider::Init.
1031 // 2) Due to the internals of Pythia, there is no need to specify particles
1032 //    and their corresponding anti-particles separately as (un)stable.
1033 //    Once a particle has been declared (un)stable, the corresponding 
1034 //    anti-particle will be treated in the same way.
1035
1036  if (mode==0 || mode==1)
1037  {
1038   Int_t kc=Pycomp(id);
1039   Int_t decay=1-mode;
1040   if (kc>0)
1041   {
1042    SetMDCY(kc,1,decay);
1043   }
1044   else 
1045   {
1046    cout << " *AliCollider::SetStable* Unknown particle code. id = " << id << endl;
1047   }
1048  }
1049  else
1050  {
1051   cout << " *AliCollider::SetStable* Invalid parameter. mode = " << mode << endl;
1052  }
1053 }
1054 ///////////////////////////////////////////////////////////////////////////
1055 void AliCollider::SelectEvent(Int_t id)
1056 {
1057 // Add a particle to the event selection list.
1058 // The parameter "id" indicates the Pythia KF particle code, which
1059 // basically is the PDG particle identifier code.
1060 // In case the user has built a selection list via this procedure, only the
1061 // events in which one of the particles specified in the list was generated
1062 // will be kept. 
1063 // The investigation of the generated particles takes place when the complete
1064 // event is in memory, including all (shortlived) mother particles and resonances.
1065 // So, the settings of the various particle decay modes have no influence on
1066 // the event selection described here.
1067 //
1068 // If no list has been specified, all events will be accepted.  
1069 //
1070 // Note : id=0 will delete the selection list.
1071 //
1072 // Be aware of the fact that severe selection criteria (i.e. selecting only
1073 // rare events) may result in long runtimes before an event sample has been
1074 // obtained.
1075 //
1076  if (!id)
1077  {
1078   if (fSelections)
1079   {
1080    delete fSelections;
1081    fSelections=0;
1082   }
1083  }
1084  else
1085  {
1086   Int_t kc=Pycomp(id);
1087   if (!fSelections)
1088   {
1089    fSelections=new TArrayI(1);
1090    fSelections->AddAt(kc,0);
1091   }
1092   else
1093   {
1094    Int_t exist=0;
1095    Int_t size=fSelections->GetSize();
1096    for (Int_t i=0; i<size; i++)
1097    {
1098     if (kc==fSelections->At(i))
1099     {
1100      exist=1;
1101      break;
1102     }
1103    }
1104   
1105    if (!exist)
1106    {
1107     fSelections->Set(size+1);
1108     fSelections->AddAt(kc,size);
1109    }
1110   }
1111  }
1112 }
1113 ///////////////////////////////////////////////////////////////////////////
1114 Int_t AliCollider::GetSelectionFlag() const
1115 {
1116 // Return the value of the selection flag for the total event.
1117 // When the event passed the selection criteria as specified via
1118 // SelectEvent() the value 1 is returned, otherwise the value 0 is returned.
1119  return fSelect;
1120 }
1121 ///////////////////////////////////////////////////////////////////////////
1122 Int_t AliCollider::IsSelected()
1123 {
1124 // Check whether the generated (sub)event contains one of the particles
1125 // specified in the selection list via SelectEvent().
1126 // If this is the case or when no selection list is present, the value 1
1127 // will be returned, indicating the event is selected to be kept.
1128 // Otherwise the value 0 will be returned.
1129
1130  if (!fSelections) return 1; 
1131
1132  Int_t nsel=fSelections->GetSize();
1133  Int_t npart=GetN();
1134  Int_t kf,kc;
1135
1136  Int_t select=0;
1137  for (Int_t jpart=1; jpart<=npart; jpart++)
1138  {
1139   kf=GetK(jpart,2);
1140   kc=Pycomp(kf);
1141   for (Int_t i=0; i<nsel; i++)
1142   {
1143    if (kc==fSelections->At(i))
1144    {
1145     select=1;
1146     break;
1147    }
1148   }
1149   if (select) break;
1150  }
1151  return select;
1152 }
1153 ///////////////////////////////////////////////////////////////////////////
1154 void AliCollider::SetSpectatorPmin(Float_t pmin)
1155 {
1156 // Set minimal momentum in GeV/c for spectator tracks to be stored.
1157 // Spectator tracks with a momentum below this threshold will not be stored
1158 // in the (output) event structure.
1159 // This facility allows to minimise the output file size.
1160 // Note that when the user wants to boost the event into another reference
1161 // frame these spectator tracks might have got momenta above the threshold.
1162 // However, when the spectator tracks were not stored in the event structure
1163 // in the original frame, there is no way to retreive them anymore. 
1164  fSpecpmin=pmin;
1165 }
1166 ///////////////////////////////////////////////////////////////////////////
1167 Float_t AliCollider::GetSpectatorPmin() const
1168 {
1169 // Provide the minimal spectator momentum in GeV/c.
1170  return fSpecpmin;
1171 }
1172 ///////////////////////////////////////////////////////////////////////////
1173 TString AliCollider::GetPyname(Int_t kf)
1174 {
1175 // Provide the correctly truncated Pythia particle name for PGD code kf
1176 //
1177 // The TPythia6::Pyname returned name is copied into a TString and truncated
1178 // at the first blank to prevent funny trailing characters due to incorrect
1179 // stripping of empty characters in TPythia6::Pyname.
1180 // The truncation at the first blank is allowed due to the Pythia convention
1181 // that particle names never contain blanks.
1182  char name[16];
1183  TString sname;
1184  Pyname(kf,name);
1185  sname=name[0];
1186  for (Int_t i=1; i<16; i++)
1187  {
1188   if (name[i]==' ') break;
1189   sname=sname+name[i];
1190  }
1191  return sname;
1192 }
1193 ///////////////////////////////////////////////////////////////////////////