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.6 2003/08/29 09:05:11 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 //  Int_t nevents=5;
  57 //
  58 //  AliRandom rndm;
  59 //  Float_t* rans=new Float_t[nevents];
  60 //  rndm.Uniform(rans,nevents,2,ap+at);
  61 //  Int_t npart;
  62 //  for (Int_t i=0; i<nevents; i++)
  63 //  {
  64 //   npart=rans[i];
  65 //   gen->MakeEvent(npart);
  66 //
  67 //   AliEvent* evt=gen->GetEvent();
  68 //
  69 //   evt->List();
  70 //  }
  71 //
  72 //  gen->EndRun();
  73 // }
  74 //
  75 //
  76 // Example job of a cosmic nu+p atmospheric interaction.
  77 // -----------------------------------------------------
  78 // {
  79 //  gSystem->Load("libEG");
  80 //  gSystem->Load("libEGPythia6");
  81 //  gSystem->Load("ralice");
  82 //
  83 //  AliCollider* gen=new AliCollider();
  84 //
  85 //  gen->SetOutputFile("test.root");
  86 //
  87 //  gen->SetRunNumber(1);
  88 //
  89 //  gen->Init("fixt","nu_mu","p",1e11);
  90 //
  91 //  Int_t nevents=10;
  92 //
  93 //  for (Int_t i=0; i<nevents; i++)
  94 //  {
  95 //   gen->MakeEvent(0,1);
  96 //
  97 //   AliEvent* evt=gen->GetEvent();
  98 //
  99 //   evt->Data();
 100 //  }
 101 //
 102 //  gen->EndRun();
 103 // }
 104 //
 105 //
 106 //--- Author: Nick van Eijndhoven 22-nov-2002 Utrecht University
 107 //- Modified: NvE $Date: 2003/08/29 09:05:11 $ Utrecht University
 108 ///////////////////////////////////////////////////////////////////////////
 109
 110 #include "AliCollider.h"
 111 #include "Riostream.h"
 112
 113 ClassImp(AliCollider) // Class implementation to enable ROOT I/O
 114
 115 AliCollider::AliCollider() : TPythia6()
 116 {
 117 // Default constructor.
 118 // All variables initialised to default values.
 119  fVertexmode=0;    // No vertex structure creation
 120  fResolution=1e-5; // Standard resolution is 0.1 micron
 121  fRunnum=0;
 122  fEventnum=0;
 123  fPrintfreq=1;
 124
 125  fEvent=0;
 126
 127  fFrame="none";
 128  fWin=0;
 129
 130  fNucl=0;
 131  fZproj=0;
 132  fAproj=0;
 133  fZtarg=0;
 134  fAtarg=0;
 135  fFracpp=0;
 136  fFracnp=0;
 137  fFracpn=0;
 138  fFracnn=0;
 139
 140  fOutFile=0;
 141  fOutTree=0;
 142 }
 143 ///////////////////////////////////////////////////////////////////////////
 144 AliCollider::~AliCollider()
 145 {
 146 // Default destructor
 147  if (fEvent)
 148  {
 149   delete fEvent;
 150   fEvent=0;
 151  }
 152  if (fOutFile)
 153  {
 154   delete fOutFile;
 155   fOutFile=0;
 156  }
 157  if (fOutTree)
 158  {
 159   delete fOutTree;
 160   fOutTree=0;
 161  }
 162 }
 163 ///////////////////////////////////////////////////////////////////////////
 164 void AliCollider::SetOutputFile(TString s)
 165 {
 166 // Create the output file containing all the data in ROOT output format.
 167  if (fOutFile)
 168  {
 169   delete fOutFile;
 170   fOutFile=0;
 171  }
 172  fOutFile=new TFile(s.Data(),"RECREATE","AliCollider data");
 173
 174  if (fOutTree)
 175  {
 176   delete fOutTree;
 177   fOutTree=0;
 178  }
 179  fOutTree=new TTree("T","AliCollider event data");
 180
 181  Int_t bsize=32000;
 182  Int_t split=0;
 183  fOutTree->Branch("Events","AliEvent",&fEvent,bsize,split);
 184 }
 185 ///////////////////////////////////////////////////////////////////////////
 186 void AliCollider::SetVertexMode(Int_t mode)
 187 {
 188 // Set the mode of the vertex structure creation.
 189 //
 190 // By default all generated tracks will only appear in the AliEvent
 191 // structure without any primary (and secondary) vertex structure.
 192 // The user can build the vertex structure if he/she wants by means
 193 // of the beginpoint location of each AliTrack.
 194 //
 195 // However, one can also let AliCollider automatically create
 196 // the primary (and secondary) vertex structure(s).
 197 // In this case the primary vertex is given Id=1 and all sec. vertices
 198 // are given Id's 2,3,4,....
 199 // All vertices are created as standalone entities in the AliEvent structure
 200 // without any linking between the various vertices.
 201 // For this automated process, the user-selected resolution
 202 // (see SetResolution) is used to decide whether or not certain vertex
 203 // locations can be resolved.
 204 // In case no vertex creation is selected (i.e. the default mode=0),
 205 // the value of the resolution is totally irrelevant.
 206 //
 207 // The user can also let AliCollider automatically connect the sec. vertices
 208 // to the primary vertex (i.e. mode=3). This process will also automatically
 209 // generate the tracks connecting the vertices.
 210 // Note that the result of the mode=3 operation may be very sensitive to
 211 // the resolution parameter. Therefore, no attempt is made to distinguish
 212 // between secondary, tertiary etc... vertices. All sec. vertices are
 213 // linked to the primary one.
 214 //
 215 // Irrespective of the selected mode, all generated tracks can be obtained
 216 // directly from the AliEvent structure.
 217 // In case (sec.) vertex creation is selected, all generated vertices can
 218 // also be obtained directly from the AliEvent structure.
 219 // These (sec.) vertices contain only the corresponding pointers to the various
 220 // tracks which are stored in the AliEvent structure.
 221 //
 222 // Overview of vertex creation modes :
 223 // -----------------------------------
 224 // mode = 0 ==> No vertex structure will be created
 225 //        1 ==> Only primary vertex structure will be created
 226 //        2 ==> Unconnected primary and secondary vertices will be created
 227 //        3 ==> Primary and secondary vertices will be created where all the
 228 //              sec. vertices will be connected to the primary vertex.
 229 //              Also the vertex connecting tracks will be automatically
 230 //              generated.
 231 //
 232  if (mode<0 || mode >3)
 233  {
 234   cout << " *AliCollider::SetVertexMode* Invalid argument mode : " << mode << endl;
 235   fVertexmode=0;
 236  }
 237  else
 238  {
 239   fVertexmode=mode;
 240  }
 241 }
 242 ///////////////////////////////////////////////////////////////////////////
 243 Int_t AliCollider::GetVertexMode()
 244 {
 245 // Provide the current mode for vertex structure creation.
 246  return fVertexmode;
 247 }
 248 ///////////////////////////////////////////////////////////////////////////
 249 void AliCollider::SetResolution(Double_t res)
 250 {
 251 // Set the resolution (in cm) for resolving (sec.) vertices.
 252 // By default this resolution is set to 0.1 micron.
 253 // Note : In case no vertex creation has been selected, the value of
 254 //        the resolution is totally irrelevant.
 255  fResolution=fabs(res);
 256 }
 257 ///////////////////////////////////////////////////////////////////////////
 258 Double_t AliCollider::GetResolution()
 259 {
 260 // Provide the current resolution (in cm) for resolving (sec.) vertices.
 261  return fResolution;
 262 }
 263 ///////////////////////////////////////////////////////////////////////////
 264 void AliCollider::SetRunNumber(Int_t run)
 265 {
 266 // Set the user defined run number.
 267 // By default the run number is set to 0.
 268  fRunnum=run;
 269 }
 270 ///////////////////////////////////////////////////////////////////////////
 271 Int_t AliCollider::GetRunNumber()
 272 {
 273 // Provide the user defined run number.
 274  return fRunnum;
 275 }
 276 ///////////////////////////////////////////////////////////////////////////
 277 void AliCollider::SetPrintFreq(Int_t n)
 278 {
 279 // Set the print frequency for every 'n' events.
 280 // By default the printfrequency is set to 1 (i.e. every event).
 281  fPrintfreq=n;
 282 }
 283 ///////////////////////////////////////////////////////////////////////////
 284 Int_t AliCollider::GetPrintFreq()
 285 {
 286 // Provide the user selected print frequency.
 287  return fPrintfreq;
 288 }
 289 ///////////////////////////////////////////////////////////////////////////
 290 void AliCollider::Init(char* frame,char* beam,char* target,Float_t win)
 291 {
 292 // Initialisation of the underlying Pythia generator package.
 293 // This routine just invokes TPythia6::Initialize(...) and the arguments
 294 // have the corresponding meaning.
 295 // The event number is reset to 0.
 296  fEventnum=0;
 297  fNucl=0;
 298  fFrame=frame;
 299  fWin=win;
 300  Initialize(frame,beam,target,win);
 301
 302  cout << " *AliCollider::Init* Standard Pythia initialisation." << endl;
 303  cout << " Beam particle : " << beam << " Target particle : " << target
 304       << " Frame = " << frame << " Energy = " << win
 305       << endl;
 306 }
 307 ///////////////////////////////////////////////////////////////////////////
 308 void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win)
 309 {
 310 // Initialisation of the underlying Pythia generator package for the generation
 311 // of nucleus-nucleus interactions.
 312 // In addition to the Pythia standard arguments 'frame' and 'win', the user
 313 // can specify here (Z,A) values of the projectile and target nuclei.
 314 //
 315 // Note : The 'win' value denotes either the cms energy per nucleon-nucleon collision
 316 //        (i.e. frame="cms") or the momentum per nucleon in all other cases.
 317 //
 318 // The event number is reset to 0.
 319  fEventnum=0;
 320  fNucl=1;
 321  fFrame=frame;
 322  fWin=win;
 323  fZproj=0;
 324  fAproj=0;
 325  fZtarg=0;
 326  fAtarg=0;
 327  fFracpp=0;
 328  fFracnp=0;
 329  fFracpn=0;
 330  fFracnn=0;
 331
 332  if (ap<1 || at<1 || zp>ap || zt>at)
 333  {
 334   cout << " *AliCollider::Init* Invalid input value(s). Zproj = " << zp
 335        << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
 336   return;
 337  }
 338
 339  fZproj=zp;
 340  fAproj=ap;
 341  fZtarg=zt;
 342  fAtarg=at;
 343
 344  cout << " *AliCollider::Init* Nucleus-Nucleus generator initialisation." << endl;
 345  cout << " Zproj = " << zp << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at
 346       << " Frame = " << frame << " Energy = " << win
 347       << endl;
 348 }
 349 ///////////////////////////////////////////////////////////////////////////
 350 void AliCollider::GetFractions(Float_t zp,Float_t ap,Float_t zt,Float_t at)
 351 {
 352 // Determine the fractions for the various N-N collision processes.
 353 // The various processes are : p+p, n+p, p+n and n+n.
 354  if (zp<0) zp=0;
 355  if (zt<0) zt=0;
 356
 357  fFracpp=0;
 358  fFracnp=0;
 359  fFracpn=0;
 360  fFracnn=0;
 361
 362  if (ap>0 && at>0)
 363  {
 364   fFracpp=(zp/ap)*(zt/at);
 365   fFracnp=(1.-zp/ap)*(zt/at);
 366   fFracpn=(zp/ap)*(1.-zt/at);
 367   fFracnn=(1.-zp/ap)*(1.-zt/at);
 368  }
 369 }
 370 ///////////////////////////////////////////////////////////////////////////
 371 void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
 372 {
 373 // Generate one event.
 374 // In case of a nucleus-nucleus interaction, the argument 'npt' denotes
 375 // the number of participant nucleons.
 376 // In case of a standard Pythia run for 'elementary' particle interactions,
 377 // the value of npt is totally irrelevant.
 378 //
 379 // The argument 'mlist' denotes the list mode used for Pylist().
 380 // Note : mlist<0 suppresses the invokation of Pylist().
 381 // By default, no listing is produced (i.e. mlist=-1).
 382 //
 383 // The argument 'medit' denotes the edit mode used for Pyedit().
 384 // Note : medit<0 suppresses the invokation of Pyedit().
 385 // By default, only 'stable' final particles are kept (i.e. medit=1).
 386 //
 387 // In the case of a standard Pythia run concerning 'elementary' particle
 388 // interactions, the projectile and target particle ID's for the created
 389 // event structure are set to the corresponding Pythia KF codes.
 390 // All the A and Z values are in that case set to zero.
 391 // In case of a nucleus-nucleus interaction, the proper A and Z values for
 392 // the projectile and target particles are set in the event structure.
 393 // However, in this case both particle ID's are set to zero.
 394
 395  fEventnum++;
 396
 397  // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
 398  Int_t ncols[4]={0,0,0,0};
 399
 400  Int_t ncol=1;
 401  if (fNucl)
 402  {
 403   if (npt<1 || npt>(fAproj+fAtarg))
 404   {
 405    cout << " *AliCollider::MakeEvent* Invalid input value. npt = " << npt
 406         << " Aproj = " << fAproj << " Atarg = " << fAtarg << endl;
 407    return;
 408   }
 409
 410   // Determine the number of nucleon-nucleon collisions
 411   ncol=npt/2;
 412   if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
 413
 414   // Determine the number of the various types of N+N interactions
 415   Int_t zp=fZproj;
 416   Int_t ap=fAproj;
 417   Int_t zt=fZtarg;
 418   Int_t at=fAtarg;
 419   Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
 420   if (ap>at) maxa=1;
 421   Float_t* rans=new Float_t[ncol];
 422   fRan.Uniform(rans,ncol);
 423   Float_t rndm=0;
 424   for (Int_t i=0; i<ncol; i++)
 425   {
 426    GetFractions(zp,ap,zt,at);
 427    rndm=rans[i];
 428    if (rndm<=fFracpp) // p+p interaction
 429    {
 430     ncols[0]++;
 431     if (maxa==2)
 432     {
 433      at--;
 434      zt--;
 435     }
 436     else
 437     {
 438      ap--;
 439      zp--;
 440     }
 441    }
 442    if (rndm>fFracpp && rndm<=(fFracpp+fFracnp)) // n+p interaction
 443    {
 444     ncols[1]++;
 445     if (maxa==2)
 446     {
 447      at--;
 448      zt--;
 449     }
 450     else
 451     {
 452      ap--;
 453     }
 454    }
 455    if (rndm>(fFracpp+fFracnp) && rndm<=(fFracpp+fFracnp+fFracpn)) // p+n interaction
 456    {
 457     ncols[2]++;
 458     if (maxa==2)
 459     {
 460      at--;
 461     }
 462     else
 463     {
 464      ap--;
 465      zp--;
 466     }
 467    }
 468    if (rndm>(fFracpp+fFracnp+fFracpn)) // n+n interaction
 469    {
 470     ncols[3]++;
 471     if (maxa==2)
 472     {
 473      at--;
 474     }
 475     else
 476     {
 477      ap--;
 478     }
 479    }
 480   }
 481   delete [] rans;
 482  }
 483
 484  if (!(fEventnum%fPrintfreq))
 485  {
 486   cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
 487        << endl;
 488   if (fNucl)
 489   {
 490    cout << " npart = " << npt << " ncol = " << ncol
 491         << " ncolpp = " << ncols[0] << " ncolnp = " << ncols[1]
 492         << " ncolpn = " << ncols[2] << " ncolnn = " << ncols[3] << endl;
 493   }
 494  }
 495
 496  if (!fEvent)
 497  {
 498   fEvent=new AliEvent();
 499   fEvent->SetOwner();
 500  }
 501
 502  fEvent->Reset();
 503  fEvent->SetRunNumber(fRunnum);
 504  fEvent->SetEventNumber(fEventnum);
 505
 506  AliTrack t;
 507  Ali3Vector p;
 508  AliPosition r,rx;
 509  Float_t v[3];
 510  AliVertex vert;
 511
 512  if (fVertexmode)
 513  {
 514   // Make sure the primary vertex gets correct location and Id=1
 515   v[0]=0;
 516   v[1]=0;
 517   v[2]=0;
 518   r.SetPosition(v,"car");
 519   v[0]=fResolution;
 520   v[1]=fResolution;
 521   v[2]=fResolution;
 522   r.SetPositionErrors(v,"car");
 523
 524   vert.SetId(1);
 525   vert.SetTrackCopy(0);
 526   vert.SetVertexCopy(0);
 527   vert.SetPosition(r);
 528   fEvent->AddVertex(vert,0);
 529  }
 530
 531  Int_t kf=0;
 532  Float_t charge=0,mass=0;
 533
 534  Int_t ntypes=4;
 535
 536  // Singular settings for a normal Pythia elementary particle interation
 537  if (!fNucl)
 538  {
 539   ntypes=1;
 540   ncols[0]=1;
 541  }
 542
 543  // Generate all the various collisions
 544  Int_t first=1; // Flag to indicate the first collision process
 545  Double_t pnucl;
 546  Int_t npart=0,ntk=0;
 547  Double_t dist=0;
 548  for (Int_t itype=0; itype<ntypes; itype++)
 549  {
 550   if (fNucl)
 551   {
 552    if (itype==0 && ncols[itype]) Initialize(fFrame,"p","p",fWin);
 553    if (itype==1 && ncols[itype]) Initialize(fFrame,"n","p",fWin);
 554    if (itype==2 && ncols[itype]) Initialize(fFrame,"p","n",fWin);
 555    if (itype==3 && ncols[itype]) Initialize(fFrame,"n","n",fWin);
 556   }
 557   for (Int_t jcol=0; jcol<ncols[itype]; jcol++)
 558   {
 559    GenerateEvent();
 560
 561    if (first) // Store projectile and target information in the event structure
 562    {
 563     if (fNucl)
 564     {
 565      v[0]=GetP(1,1);
 566      v[1]=GetP(1,2);
 567      v[2]=GetP(1,3);
 568      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
 569      fEvent->SetProjectile(fAproj,fZproj,pnucl);
 570      v[0]=GetP(2,1);
 571      v[1]=GetP(2,2);
 572      v[2]=GetP(2,3);
 573      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
 574      fEvent->SetTarget(fAtarg,fZtarg,pnucl);
 575     }
 576     else
 577     {
 578      v[0]=GetP(1,1);
 579      v[1]=GetP(1,2);
 580      v[2]=GetP(1,3);
 581      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
 582      kf=GetK(1,2);
 583      fEvent->SetProjectile(0,0,pnucl,kf);
 584      v[0]=GetP(2,1);
 585      v[1]=GetP(2,2);
 586      v[2]=GetP(2,3);
 587      pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
 588      kf=GetK(2,2);
 589      fEvent->SetTarget(0,0,pnucl,kf);
 590     }
 591     first=0;
 592    }
 593
 594    if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
 595
 596    if (mlist >= 0) Pylist(mlist);
 597
 598    npart=GetN();
 599    for (Int_t jpart=1; jpart<=npart; jpart++)
 600    {
 601     kf=GetK(jpart,2);
 602     charge=Pychge(kf)/3.;
 603     mass=GetP(jpart,5);
 604
 605     // 3-momentum in GeV/c
 606     v[0]=GetP(jpart,1);
 607     v[1]=GetP(jpart,2);
 608     v[2]=GetP(jpart,3);
 609     p.SetVector(v,"car");
 610
 611     // Production location in cm.
 612     v[0]=GetV(jpart,1)/10;
 613     v[1]=GetV(jpart,2)/10;
 614     v[2]=GetV(jpart,3)/10;
 615     r.SetPosition(v,"car");
 616
 617     ntk++;
 618
 619     t.Reset();
 620     t.SetId(ntk);
 621     t.SetParticleCode(kf);
 622     t.SetCharge(charge);
 623     t.SetMass(mass);
 624     t.Set3Momentum(p);
 625     t.SetBeginPoint(r);
 626
 627     fEvent->AddTrack(t);
 628
 629     // Build the vertex structures if requested
 630     if (fVertexmode)
 631     {
 632      // Check if track belongs within the resolution to an existing vertex
 633      Int_t add=0;
 634      for (Int_t jv=1; jv<=fEvent->GetNvertices(); jv++)
 635      {
 636       AliVertex* vx=fEvent->GetVertex(jv);
 637       if (vx)
 638       {
 639        rx=vx->GetPosition();
 640        dist=rx.GetDistance(r);
 641        if (dist < fResolution)
 642        {
 643         AliTrack* tx=fEvent->GetIdTrack(ntk);
 644         if (tx)
 645         {
 646          vx->AddTrack(tx);
 647          add=1;
 648         }
 649        }
 650       }
 651       if (add) break; // No need to look further for vertex candidates
 652      }
 653
 654      // If track was not close enough to an existing vertex
 655      // a new secondary vertex is created
 656      if (!add && fVertexmode>1)
 657      {
 658       AliTrack* tx=fEvent->GetIdTrack(ntk);
 659       if (tx)
 660       {
 661        v[0]=fResolution;
 662        v[1]=fResolution;
 663        v[2]=fResolution;
 664        r.SetPositionErrors(v,"car");
 665        vert.Reset();
 666        vert.SetTrackCopy(0);
 667        vert.SetVertexCopy(0);
 668        vert.SetId((fEvent->GetNvertices())+1);
 669        vert.SetPosition(r);
 670        vert.AddTrack(tx);
 671        fEvent->AddVertex(vert,0);
 672       }
 673      }
 674     }
 675    } // End of loop over the produced particles for each collision
 676   } // End of loop over number of collisions for each type
 677  } // End of loop over collision types
 678
 679  // Link sec. vertices to primary if requested
 680  // Note that also the connecting tracks are automatically created
 681  if (fVertexmode>2)
 682  {
 683   AliVertex* vp=fEvent->GetIdVertex(1); // Primary vertex
 684   if (vp)
 685   {
 686    for (Int_t i=2; i<=fEvent->GetNvertices(); i++)
 687    {
 688     AliVertex* vx=fEvent->GetVertex(i);
 689     if (vx)
 690     {
 691      if (vx->GetId() != 1) vp->AddVertex(vx);
 692     }
 693    }
 694   }
 695  }
 696
 697  if (mlist && !(fEventnum%fPrintfreq)) cout << endl; // Create empty output line after the event
 698  if (fOutTree) fOutTree->Fill();
 699 }
 700 ///////////////////////////////////////////////////////////////////////////
 701 AliEvent* AliCollider::GetEvent()
 702 {
 703 // Provide pointer to the generated event structure.
 704  return fEvent;
 705 }
 706 ///////////////////////////////////////////////////////////////////////////
 707 void AliCollider::EndRun()
 708 {
 709 // Properly close the output file (if needed).
 710  if (fOutFile)
 711  {
 712   fOutFile->Write();
 713   fOutFile->Close();
 714   cout << " *AliCollider::EndRun* Output file correctly closed." << endl;
 715  }
 716 }
 717 ///////////////////////////////////////////////////////////////////////////