1 // $Id: AliCollider.cxx,v 1.2 2002/11/29 13:54:52 nick Exp $
3 ///////////////////////////////////////////////////////////////////////////
5 // Pythia based universal physics event generator.
6 // This class is derived from TPythia6 and has some extensions to
7 // support also generation of nucleus-nucleus interactions and to allow
8 // investigation of the effect of detector resolving power.
9 // Furthermore, the produced event information is provided in a format
10 // using the AliEvent structure.
11 // For the produced AliTrack objects, the particle ID code is set to the
12 // Pythia KF value, which is compatible with the PDG identifier.
13 // This will allow a direct analysis of the produced data using the
14 // Ralice physics analysis tools.
16 // For further details concerning the produced output structure,
17 // see the docs of the memberfunctions SetVertexMode and SetResolution.
19 // Example job of minimum biased Pb+Pb interactions :
20 // --------------------------------------------------
22 // gSystem->Load("libEG");
23 // gSystem->Load("libEGPythia6");
24 // gSystem->Load("ralice");
26 // AliCollider* gen=new AliCollider();
28 // gen->SetOutputFile("test.root");
29 // gen->SetVertexMode(3);
30 // gen->SetResolution(1e-4); // 1 micron vertex resolution
32 // gen->SetRunNumber(1);
39 // gen->Init("fixt",zp,ap,zt,at,158);
44 // Float_t* rans=new Float_t[nevents];
45 // rndm.Uniform(rans,nevents,2,ap+at);
47 // for (Int_t i=0; i<nevents; i++)
50 // gen->MakeEvent(npart);
52 // AliEvent* evt=gen->GetEvent();
61 // Example job of a cosmic nu+p atmospheric interaction.
62 // -----------------------------------------------------
64 // gSystem->Load("libEG");
65 // gSystem->Load("libEGPythia6");
66 // gSystem->Load("ralice");
68 // AliCollider* gen=new AliCollider();
70 // gen->SetOutputFile("test.root");
72 // gen->SetRunNumber(1);
74 // gen->Init("fixt","nu_mu","p",1e11);
78 // for (Int_t i=0; i<nevents; i++)
80 // gen->MakeEvent(0,1);
82 // AliEvent* evt=gen->GetEvent();
91 //--- Author: Nick van Eijndhoven 22-nov-2002 Utrecht University
92 //- Modified: NvE $Date: 2002/11/29 13:54:52 $ Utrecht University
93 ///////////////////////////////////////////////////////////////////////////
95 #include "AliCollider.h"
97 ClassImp(AliCollider) // Class implementation to enable ROOT I/O
99 AliCollider::AliCollider()
101 // Default constructor.
102 // All variables initialised to default values.
103 fVertexmode=0; // No vertex structure creation
104 fResolution=1e-5; // Standard resolution is 0.1 micron
127 ///////////////////////////////////////////////////////////////////////////
128 AliCollider::~AliCollider()
130 // Default destructor
147 ///////////////////////////////////////////////////////////////////////////
148 void AliCollider::SetOutputFile(TString s)
150 // Create the output file containing all the data in ROOT output format.
156 fOutFile=new TFile(s.Data(),"RECREATE","AliCollider data");
163 fOutTree=new TTree("T","AliCollider event data");
167 fOutTree->Branch("Events","AliEvent",&fEvent,bsize,split);
169 ///////////////////////////////////////////////////////////////////////////
170 void AliCollider::SetVertexMode(Int_t mode)
172 // Set the mode of the vertex structure creation.
174 // By default all generated tracks will only appear in the AliEvent
175 // structure without any primary (and secondary) vertex structure.
176 // The user can build the vertex structure if he/she wants by means
177 // of the beginpoint location of each AliTrack.
179 // However, one can also let AliCollider automatically create
180 // the primary (and secondary) vertex structure(s).
181 // In this case the primary vertex is given Id=1 and all sec. vertices
182 // are given Id's 2,3,4,....
183 // All vertices are created as standalone entities in the AliEvent structure
184 // without any linking between the various vertices.
185 // For this automated process, the user-selected resolution
186 // (see SetResolution) is used to decide whether or not certain vertex
187 // locations can be resolved.
188 // In case no vertex creation is selected (i.e. the default mode=0),
189 // the value of the resolution is totally irrelevant.
191 // The user can also let AliCollider automatically connect the sec. vertices
192 // to the primary vertex (i.e. mode=3). This process will also automatically
193 // generate the tracks connecting the vertices.
194 // Note that the result of the mode=3 operation may be very sensitive to
195 // the resolution parameter. Therefore, no attempt is made to distinguish
196 // between secondary, tertiary etc... vertices. All sec. vertices are
197 // linked to the primary one.
199 // Irrespective of the selected mode, all generated tracks can be obtained
200 // directly from the AliEvent structure.
201 // In case (sec.) vertex creation is selected, all generated vertices can
202 // also be obtained directly from the AliEvent structure.
203 // These (sec.) vertices contain only the corresponding pointers to the various
204 // tracks which are stored in the AliEvent structure.
206 // Overview of vertex creation modes :
207 // -----------------------------------
208 // mode = 0 ==> No vertex structure will be created
209 // 1 ==> Only primary vertex structure will be created
210 // 2 ==> Unconnected primary and secondary vertices will be created
211 // 3 ==> Primary and secondary vertices will be created where all the
212 // sec. vertices will be connected to the primary vertex.
213 // Also the vertex connecting tracks will be automatically
216 if (mode<0 || mode >3)
218 cout << " *AliCollider::SetVertexMode* Invalid argument mode : " << mode << endl;
226 ///////////////////////////////////////////////////////////////////////////
227 Int_t AliCollider::GetVertexMode()
229 // Provide the current mode for vertex structure creation.
232 ///////////////////////////////////////////////////////////////////////////
233 void AliCollider::SetResolution(Double_t res)
235 // Set the resolution (in cm) for resolving (sec.) vertices.
236 // By default this resolution is set to 0.1 micron.
237 // Note : In case no vertex creation has been selected, the value of
238 // the resolution is totally irrelevant.
239 fResolution=fabs(res);
241 ///////////////////////////////////////////////////////////////////////////
242 Double_t AliCollider::GetResolution()
244 // Provide the current resolution (in cm) for resolving (sec.) vertices.
247 ///////////////////////////////////////////////////////////////////////////
248 void AliCollider::SetRunNumber(Int_t run)
250 // Set the user defined run number.
251 // By default the run number is set to 0.
254 ///////////////////////////////////////////////////////////////////////////
255 Int_t AliCollider::GetRunNumber()
257 // Provide the user defined run number.
260 ///////////////////////////////////////////////////////////////////////////
261 void AliCollider::SetPrintFreq(Int_t n)
263 // Set the print frequency for every 'n' events.
264 // By default the printfrequency is set to 1 (i.e. every event).
267 ///////////////////////////////////////////////////////////////////////////
268 Int_t AliCollider::GetPrintFreq()
270 // Provide the user selected print frequency.
273 ///////////////////////////////////////////////////////////////////////////
274 void AliCollider::Init(char* frame,char* beam,char* target,Float_t win)
276 // Initialisation of the underlying Pythia generator package.
277 // This routine just invokes TPythia6::Initialize(...) and the arguments
278 // have the corresponding meaning.
279 // The event number is reset to 0.
284 Initialize(frame,beam,target,win);
286 ///////////////////////////////////////////////////////////////////////////
287 void AliCollider::Init(char* frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win)
289 // Initialisation of the underlying Pythia generator package for the generation
290 // of nucleus-nucleus interactions.
291 // In addition to the Pythia standard arguments 'frame' and 'win', the user
292 // can specify here (Z,A) values of the projectile and target nuclei.
294 // Note : The 'win' value denotes either the cms energy per nucleon-nucleon collision
295 // (i.e. frame="cms") or the momentum per nucleon in all other cases.
297 // The event number is reset to 0.
311 if (ap<1 || at<1 || zp>ap || zt>at)
313 cout << " *AliCollider::Init* Invalid input value(s). Zproj = " << zp
314 << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
323 cout << " *AliCollider::Init* Nucleus-Nucleus generator initialisation." << endl;
324 cout << " Zproj = " << zp << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at
325 << " Frame = " << frame << " Energy = " << win
328 ///////////////////////////////////////////////////////////////////////////
329 void AliCollider::GetFractions(Float_t zp,Float_t ap,Float_t zt,Float_t at)
331 // Determine the fractions for the various N-N collision processes.
332 // The various processes are : p+p, n+p, p+n and n+n.
343 fFracpp=(zp/ap)*(zt/at);
344 fFracnp=(1.-zp/ap)*(zt/at);
345 fFracpn=(zp/ap)*(1.-zt/at);
346 fFracnn=(1.-zp/ap)*(1.-zt/at);
349 ///////////////////////////////////////////////////////////////////////////
350 void AliCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
352 // Generate one event.
353 // In case of a nucleus-nucleus interaction, the argument 'npt' denotes
354 // the number of participant nucleons.
355 // In case of a standard Pythia run for 'elementary' particle interactions,
356 // the value of npt is totally irrelevant.
358 // The argument 'medit' denotes the edit mode used for Pyedit().
359 // Note : medit<0 suppresses the invokation of Pyedit().
360 // By default, only 'stable' final particles are kept (i.e. medit=1).
362 // The argument 'mlist' denotes the list mode used for Pylist().
363 // Note : mlist<0 suppresses the invokation of Pylist().
364 // By default, no listing is produced (i.e. mlist=-1).
366 // In the case of a standard Pythia run concerning 'elementary' particle
367 // interactions, the projectile and target particle ID's for the created
368 // event structure are set to the corresponding Pythia KF codes.
369 // All the A and Z values are in that case set to zero.
370 // In case of a nucleus-nucleus interaction, the proper A and Z values for
371 // the projectile and target particles are set in the event structure.
372 // However, in this case both particle ID's are set to zero.
376 // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
377 Int_t ncols[4]={0,0,0,0};
381 if (npt<1 || npt>(fAproj+fAtarg))
383 cout << " *AliCollider::MakeEvent* Invalid input value. npt = " << npt
384 << " Aproj = " << fAproj << " Atarg = " << fAtarg << endl;
388 // Determine the number of nucleon-nucleon collisions
390 if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
392 // Determine the number of the various types of N+N interactions
397 Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
399 Float_t* rans=new Float_t[ncol];
400 fRan.Uniform(rans,ncol);
402 for (Int_t i=0; i<ncol; i++)
404 GetFractions(zp,ap,zt,at);
406 if (rndm<=fFracpp) // p+p interaction
420 if (rndm>fFracpp && rndm<=(fFracpp+fFracnp)) // n+p interaction
433 if (rndm>(fFracpp+fFracnp) && rndm<=(fFracpp+fFracnp+fFracpn)) // p+n interaction
446 if (rndm>(fFracpp+fFracnp+fFracpn)) // n+n interaction
461 if (!(fEventnum%fPrintfreq))
463 cout << " *AliCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum
465 cout << " npart = " << npt << " ncol = " << ncol
466 << " ncolpp = " << ncols[0] << " ncolnp = " << ncols[1]
467 << " ncolpn = " << ncols[2] << " ncolnn = " << ncols[3] << endl;
474 fEvent=new AliEvent();
479 fEvent->SetRunNumber(fRunnum);
480 fEvent->SetEventNumber(fEventnum);
490 // Make sure the primary vertex gets correct location and Id=1
494 r.SetPosition(v,"car");
498 r.SetPositionErrors(v,"car");
501 vert.SetTrackCopy(0);
502 vert.SetVertexCopy(0);
504 fEvent->AddVertex(vert,0);
508 Float_t charge=0,mass=0;
514 // Singular settings for a normal Pythia elementary particle interation
521 // Generate all the various collisions
522 Int_t first=1; // Flag to indicate the first collision process
526 for (Int_t itype=0; itype<ntypes; itype++)
530 if (itype==0 && ncols[itype]) Initialize(fFrame,"p","p",fWin);
531 if (itype==1 && ncols[itype]) Initialize(fFrame,"n","p",fWin);
532 if (itype==2 && ncols[itype]) Initialize(fFrame,"p","n",fWin);
533 if (itype==3 && ncols[itype]) Initialize(fFrame,"n","n",fWin);
535 for (Int_t jcol=0; jcol<ncols[itype]; jcol++)
539 if (first) // Store projectile and target information in the event structure
546 pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
547 fEvent->SetProjectile(fAproj,fZproj,pnucl);
551 pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
552 fEvent->SetTarget(fAtarg,fZtarg,pnucl);
559 pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
561 fEvent->SetProjectile(0,0,pnucl,kf);
565 pnucl=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
567 fEvent->SetTarget(0,0,pnucl,kf);
572 if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
574 if (mlist >= 0) Pylist(mlist);
578 if (fParticles) npart=fParticles->GetEntries();
580 for (Int_t jpart=0; jpart<npart; jpart++)
582 part=(TMCParticle*)fParticles->At(jpart);
588 charge=GetKCHG(kc,1)/3.;
589 if (kf<0) charge*=-1;
592 // 3-momentum in GeV/c
596 p.SetVector(v,"car");
598 // Production location in cm.
599 v[0]=(part->GetVx())/10;
600 v[1]=(part->GetVy())/10;
601 v[2]=(part->GetVz())/10;
602 r.SetPosition(v,"car");
608 t.SetParticleCode(kf);
616 // Build the vertex structures if requested
619 // Check if track belongs within the resolution to an existing vertex
621 for (Int_t jv=1; jv<=fEvent->GetNvertices(); jv++)
623 AliVertex* vx=fEvent->GetVertex(jv);
626 rx=vx->GetPosition();
627 dist=rx.GetDistance(r);
628 if (dist < fResolution)
630 AliTrack* tx=fEvent->GetIdTrack(ntk);
638 if (add) break; // No need to look further for vertex candidates
641 // If track was not close enough to an existing vertex
642 // a new secondary vertex is created
643 if (!add && fVertexmode>1)
645 AliTrack* tx=fEvent->GetIdTrack(ntk);
651 r.SetPositionErrors(v,"car");
653 vert.SetTrackCopy(0);
654 vert.SetVertexCopy(0);
655 vert.SetId((fEvent->GetNvertices())+1);
658 fEvent->AddVertex(vert,0);
662 } // End of loop over the produced particles for each collision
663 } // End of loop over number of collisions for each type
664 } // End of loop over collision types
666 // Link sec. vertices to primary if requested
667 // Note that also the connecting tracks are automatically created
670 AliVertex* vp=fEvent->GetIdVertex(1); // Primary vertex
673 for (Int_t i=2; i<=fEvent->GetNvertices(); i++)
675 AliVertex* vx=fEvent->GetVertex(i);
678 if (vx->GetId() != 1) vp->AddVertex(vx);
684 if (mlist) cout << endl; // Create empty output line after the event
685 if (fOutTree) fOutTree->Fill();
687 ///////////////////////////////////////////////////////////////////////////
688 AliEvent* AliCollider::GetEvent()
690 // Provide pointer to the generated event structure.
693 ///////////////////////////////////////////////////////////////////////////
694 void AliCollider::EndRun()
696 // Properly close the output file (if needed).
701 cout << " *AliCollider::EndRun* Output file correctly closed." << endl;
704 ///////////////////////////////////////////////////////////////////////////