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