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