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