EPOS/AliGenEpos.cxx

   1 //
   2 // AliGenEpos.cxx
   3 //
   4 //  ALICE event generator based on EPOS model from Klaus Werner
   5 //
   6 //  Created on: Feb 28, 2009
   7 //      Author: Piotr Ostrowski, postrow@if.pw.edu.pl
   8 //
   9
  10 #include "AliGenEpos.h"
  11 #include "TEpos.h"
  12 #include "TParticle.h"
  13 #include "TMath.h"
  14 #include "AliLog.h"
  15 #include "AliGenEventHeader.h"
  16 #include "AliGenEposEventHeader.h"
  17
  18 ClassImp(AliGenEpos)
  19
  20 AliGenEpos::AliGenEpos() : AliGenMC(),
  21                 fBmin(0),
  22                 fBmax(10000),
  23                 fPhiMin(0),
  24                 fPhiMax(TMath::TwoPi()),
  25                 fFilterModelOutput(kFALSE) {
  26         SetMC(new TEpos());
  27 }
  28
  29 AliGenEpos::AliGenEpos(Int_t npart) : AliGenMC(npart),
  30                 fBmin(0),
  31                 fBmax(10000),
  32                 fPhiMin(0),
  33                 fPhiMax(TMath::TwoPi()),
  34                 fFilterModelOutput(kFALSE) {
  35         SetMC(new TEpos());
  36 }
  37
  38 void AliGenEpos::Init() {
  39   // Sets up TEpos
  40         AliGenMC::Init();
  41         TEpos *epos = GetTEpos();
  42         epos->SetLaproj(this->fZProjectile);
  43         epos->SetMaproj(this->fAProjectile);
  44         epos->SetLatarg(this->fZTarget);
  45         epos->SetMatarg(this->fATarget);
  46         epos->SetPhimin(this->fPhiMin);
  47         epos->SetPhimax(this->fPhiMax);
  48         epos->SetBminim(this->fBmin);
  49         epos->SetBmaxim(this->fBmax);
  50         epos->SetEcms(this->fEnergyCMS);
  51         epos->SetSplitting(kTRUE);
  52         GetTEpos()->Initialize();
  53 }
  54
  55 void AliGenEpos::Generate() {
  56   // Does actual generation and output conversion
  57           Float_t polar[3]    =   {0,0,0};
  58           Float_t origin0[3]  =   {0,0,0};
  59           Float_t origin[3]   =   {0,0,0};
  60           Float_t time0 = 0.;
  61           Float_t time  = 0.;
  62           fNprimaries = 0;
  63           Int_t nt  = 0; //output parameter for PushTrack
  64
  65           Vertex();
  66           for (int j=0; j < 3; j++) origin0[j] = fVertex[j];
  67           time0 = fTime;
  68
  69           // Generate one event
  70
  71           GetTEpos()->GenerateEvent();
  72           //AliWarning("Generated");
  73           GetTEpos()->ImportParticles(&fParticles);
  74
  75           Int_t np = fParticles.GetEntriesFast();
  76           AliWarning(Form("Imported %d particles", np));
  77
  78           Int_t *idsOnStack = NULL;
  79           idsOnStack = new Int_t[np];
  80           for (int i = 0; i < np; i++) idsOnStack[i] = 0;
  81           TParticle *iparticle;
  82
  83           for (int i = 0; i < np; i++) {
  84                   iparticle = (TParticle *) fParticles.At(i);
  85                   //Bool_t isNullEntry = iparticle->GetStatusCode() == 0;
  86                   //Bool_t isCommentOrUnknown = iparticle->GetStatusCode() > 2;
  87                   Bool_t hasDecayed = iparticle->GetStatusCode() >= 2;
  88                   Bool_t isFinalState = iparticle->GetStatusCode() == 1;
  89                   Int_t imo = iparticle->GetFirstMother();
  90                   Bool_t  hasMother = (imo >=0);
  91
  92
  93                   if (isFinalState) {
  94                           origin[0] = iparticle->Vx();
  95                           origin[1] = iparticle->Vy();
  96                           origin[2] = iparticle->Vz();
  97                           time      = iparticle->T();
  98                           //doubled track with freeze out coordinates for femtoscopy
  99                           PushTrack(0, imo>=0?idsOnStack[imo]:-1,
 100                                 iparticle->GetPdgCode(),
 101                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
 102                                 origin[0], origin[1], origin[2], time,
 103                                 polar[0],polar[1],polar[2],
 104                                 hasMother ? kPDecay:kPNoProcess, nt, 1.,
 105                                 iparticle->GetStatusCode());
 106
 107                       idsOnStack[i] = nt;
 108                       fNprimaries++;
 109                       KeepTrack(nt);
 110
 111                       //real track with smeared vertex
 112                       origin[0] += origin0[0];
 113                       origin[1] += origin0[1];
 114                       origin[2] += origin0[2];
 115                       time      += time0;
 116                           PushTrack(1,
 117                                           nt,   //doubled track as mother
 118                                           iparticle->GetPdgCode(),
 119                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
 120                                 origin[0], origin[1], origin[2],
 121                                 time,
 122                                 polar[0],polar[1],polar[2],
 123                                 kPDecay,nt);
 124                       fNprimaries++;
 125                       KeepTrack(nt);
 126                   } else if(hasDecayed || !fFilterModelOutput) {
 127                           // don't track it and don't smear vertex
 128                           origin[0] = iparticle->Vx();
 129                           origin[1] = iparticle->Vy();
 130                           origin[2] = iparticle->Vz();
 131                           time      = iparticle->T();
 132                           PushTrack(0, imo>=0?idsOnStack[imo]:-1,
 133                                 iparticle->GetPdgCode(),
 134                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
 135                                 origin[0], origin[1], origin[2], time,
 136                                 polar[0],polar[1],polar[2],
 137                                 hasMother ? kPDecay:kPNoProcess, nt, 1.,
 138                                 iparticle->GetStatusCode());
 139                       idsOnStack[i] = nt;
 140                       fNprimaries++;
 141                       KeepTrack(nt);
 142                   } else if(fFilterModelOutput) {
 143                           //filtered internal model objects
 144                           idsOnStack[i] = -1; // to erase mother field in its dauthers
 145                   }
 146           }
 147           SetHighWaterMark(fNprimaries);
 148           TArrayF eventVertex;
 149           eventVertex.Set(3);
 150           eventVertex[0] = origin0[0];
 151           eventVertex[1] = origin0[1];
 152           eventVertex[2] = origin0[2];
 153         // Builds the event header, to be called after each event
 154           AliGenEposEventHeader* header = new AliGenEposEventHeader("EPOS");
 155
 156           header->SetNProduced(fNprimaries);
 157           header->SetPrimaryVertex(eventVertex);
 158           header->SetInteractionTime(time0);
 159
 160           header->SetImpactParameter(GetTEpos()->GetBimevt());
 161           header->SetReactionPlaneAngle(GetTEpos()->GetPhievt());
 162
 163           header->SetHardScatters(0);
 164           header->SetParticipants(GetTEpos()->GetNpjevt(), GetTEpos()->GetNtgevt());
 165           header->SetCollisions(GetTEpos()->GetKolevt(), 0, 0, GetTEpos()->GetNpjevt());
 166           header->SetSpectators(GetTEpos()->GetJpnevt(), GetTEpos()->GetJppevt(), GetTEpos()->GetJtnevt(), GetTEpos()->GetJtpevt());
 167
 168         // Event Vertex
 169           header->SetPrimaryVertex(fVertex);
 170
 171           header->SetBimevt(GetTEpos()->GetBimevt());
 172           header->SetPhievt(GetTEpos()->GetPhievt());
 173           header->SetKolevt(GetTEpos()->GetKolevt());
 174           header->SetKoievt(GetTEpos()->GetKoievt());
 175           header->SetPmxevt(GetTEpos()->GetPmxevt());
 176           header->SetEgyevt(GetTEpos()->GetEgyevt());
 177           header->SetNpjevt(GetTEpos()->GetNpjevt());
 178           header->SetNtgevt(GetTEpos()->GetNtgevt());
 179           header->SetNpnevt(GetTEpos()->GetNpnevt());
 180           header->SetNppevt(GetTEpos()->GetNppevt());
 181           header->SetNtnevt(GetTEpos()->GetNtnevt());
 182           header->SetNtpevt(GetTEpos()->GetNtpevt());
 183           header->SetJpnevt(GetTEpos()->GetJpnevt());
 184           header->SetJppevt(GetTEpos()->GetJppevt());
 185           header->SetJtnevt(GetTEpos()->GetJtnevt());
 186           header->SetJtpevt(GetTEpos()->GetJtpevt());
 187           header->SetXbjevt(GetTEpos()->GetXbjevt());
 188           header->SetQsqevt(GetTEpos()->GetQsqevt());
 189           header->SetNglevt(GetTEpos()->GetNglevt());
 190           header->SetZppevt(GetTEpos()->GetZppevt());
 191           header->SetZptevt(GetTEpos()->GetZptevt());
 192     AddHeader(header);
 193     fCollisionGeometry = (AliGenEposEventHeader*)  header;
 194
 195     delete[] idsOnStack;
 196 }
 197
 198 AliGenEpos::~AliGenEpos() {
 199 }
 200