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