EPOS/AliGenEpos.cxx

   1 /*
   2  * AliGenEpos.cpp
   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 "AliLog.h"
  14 #include "AliGenEventHeader.h"
  15 #include "AliGenEposEventHeader.h"
  16
  17 ClassImp(AliGenEpos)
  18
  19 AliGenEpos::AliGenEpos() : AliGenMC(),
  20                 fBmin(0),
  21                 fBmax(10000),
  22                 fPhiMin(0),
  23                 fPhiMax(2*3.1415926),
  24                 fFilterModelOutput(kFALSE) {
  25         SetMC(new TEpos());
  26 }
  27
  28 AliGenEpos::AliGenEpos(Int_t npart) : AliGenMC(npart),
  29                 fBmin(0),
  30                 fBmax(10000),
  31                 fPhiMin(0),
  32                 fPhiMax(2*3.1415926),
  33                 fFilterModelOutput(kFALSE) {
  34         SetMC(new TEpos());
  35 }
  36
  37 void AliGenEpos::Init() {
  38         AliGenMC::Init();
  39         TEpos *epos = GetTEpos();
  40         epos->SetLaproj(this->fZProjectile);
  41         epos->SetMaproj(this->fAProjectile);
  42         epos->SetLatarg(this->fZTarget);
  43         epos->SetMatarg(this->fATarget);
  44         epos->SetPhimin(this->fPhiMin);
  45         epos->SetPhimax(this->fPhiMax);
  46         epos->SetBminim(this->fBmin);
  47         epos->SetBmaxim(this->fBmax);
  48         epos->SetEcms(this->fEnergyCMS);
  49         GetTEpos()->Initialize();
  50 }
  51
  52 void AliGenEpos::Generate() {
  53           Float_t polar[3]   =   {0,0,0};
  54           Float_t origin0[3]  =   {0,0,0};
  55           Float_t origin[3]   =   {0,0,0};
  56
  57           Int_t nt  = 0; //output parameter for PushTrack
  58
  59           Vertex();
  60           for (int j=0; j < 3; j++) origin0[j] = fVertex[j];
  61
  62           // Generate one event
  63
  64           GetTEpos()->GenerateEvent();
  65           AliWarning("Generated");
  66           GetTEpos()->ImportParticles(&fParticles);
  67
  68           Int_t np = fParticles.GetEntriesFast();
  69           AliWarning(Form("Imported %d particles", np));
  70
  71           Int_t *idsOnStack = NULL;
  72           idsOnStack = new Int_t[np];
  73           TParticle *iparticle;
  74
  75           for (int i = 0; i < np; i++) {
  76                   iparticle = (TParticle *) fParticles.At(i);
  77                   //Bool_t isNullEntry = iparticle->GetStatusCode() == 0;
  78                   //Bool_t isCommentOrUnknown = iparticle->GetStatusCode() > 2;
  79                   Bool_t hasDecayed = iparticle->GetStatusCode() == 2;
  80                   Bool_t isFinalState = iparticle->GetStatusCode() == 1;
  81                   Int_t imo = iparticle->GetFirstMother();
  82                   Bool_t  hasMother = (imo >=0);
  83
  84
  85                   if (isFinalState) {
  86                           origin[0] = iparticle->Vx();
  87                           origin[1] = iparticle->Vy();
  88                           origin[2] = iparticle->Vz();
  89                           //doubled track with freeze out coordinates for femtoscopy
  90                           PushTrack(0,
  91                                           imo>=0?idsOnStack[imo]:-1,
  92                                           iparticle->GetPdgCode(),
  93                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
  94                                 origin[0], origin[1], origin[2],
  95                                 iparticle->T(),
  96                                 polar[0],polar[1],polar[2],
  97                                 hasMother ? kPDecay:kPNoProcess,nt);
  98
  99                       idsOnStack[i] = nt;
 100                       fNprimaries++;
 101                       KeepTrack(nt);
 102
 103                       //real track with smeared vertex
 104                       origin[0] += origin0[0];
 105                       origin[1] += origin0[1];
 106                       origin[2] += origin0[2];
 107                           PushTrack(1,
 108                                           nt,   //doubled track as mother
 109                                           iparticle->GetPdgCode(),
 110                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
 111                                 origin[0], origin[1], origin[2],
 112                                 iparticle->T(),
 113                                 polar[0],polar[1],polar[2],
 114                                 kPDecay,nt);
 115                       fNprimaries++;
 116                       KeepTrack(nt);
 117                   } else if(hasDecayed || !fFilterModelOutput) {
 118                           // don't track it and don't smear vertex
 119                           origin[0] = iparticle->Vx();
 120                           origin[1] = iparticle->Vy();
 121                           origin[2] = iparticle->Vz();
 122                           PushTrack(0,
 123                                           imo>=0?idsOnStack[imo]:-1,
 124                                           iparticle->GetPdgCode(),
 125                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
 126                                 origin[0], origin[1], origin[2],
 127                                 iparticle->T(),
 128                                 polar[0],polar[1],polar[2],
 129                                 hasMother ? kPDecay:kPNoProcess,nt);
 130                       idsOnStack[i] = nt;
 131                       fNprimaries++;
 132                       KeepTrack(nt);
 133                   } else if(fFilterModelOutput) {
 134                           //filtered internal model objects
 135                           idsOnStack[i] = -1; // to erase mother field in its dauthers
 136                   }
 137           }
 138           SetHighWaterMark(fNprimaries);
 139           TArrayF eventVertex;
 140           eventVertex.Set(3);
 141           eventVertex[0] = origin0[0];
 142           eventVertex[1] = origin0[1];
 143           eventVertex[2] = origin0[2];
 144         // Builds the event header, to be called after each event
 145           AliGenEposEventHeader* header = new AliGenEposEventHeader("EPOS");
 146
 147           header->SetNProduced(fNprimaries);
 148           header->SetPrimaryVertex(eventVertex);
 149
 150           header->SetImpactParameter(GetTEpos()->GetBimevt());
 151           header->SetReactionPlaneAngle(GetTEpos()->GetPhievt());
 152
 153           header->SetHardScatters(0);
 154           header->SetParticipants(GetTEpos()->GetNpjevt(), GetTEpos()->GetNtgevt());
 155           header->SetCollisions(GetTEpos()->GetKolevt(), 0, 0, GetTEpos()->GetNpjevt());
 156           header->SetSpectators(GetTEpos()->GetJpnevt(), GetTEpos()->GetJppevt(), GetTEpos()->GetJtnevt(), GetTEpos()->GetJtpevt());
 157
 158         // Event Vertex
 159           header->SetPrimaryVertex(fVertex);
 160           AddHeader(header);
 161           fCollisionGeometry = (AliGenEposEventHeader*)  header;
 162
 163           delete[] idsOnStack;
 164 }
 165
 166 AliGenEpos::~AliGenEpos() {
 167 }
 168