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