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