[u/mrichter/AliRoot.git] / EPOS / AliGenEpos.cxx

/*
 * AliGenEpos.cpp
 *
 *  ALICE event generator based on EPOS model from Klaus Werner
 *
 *  Created on: Feb 28, 2009
 *      Author: Piotr Ostrowski, postrow@if.pw.edu.pl
 */

#include "AliGenEpos.h"
#include "TEpos.h"
#include "TParticle.h"
#include "AliLog.h"
#include "AliGenEventHeader.h"
#include "AliGenEposEventHeader.h"

ClassImp(AliGenEpos)

AliGenEpos::AliGenEpos() : AliGenMC(),
		fBmin(0),
		fBmax(10000),
		fPhiMin(0),
		fPhiMax(2*3.1415926),
		fFilterModelOutput(kFALSE) {
	SetMC(new TEpos());
}

AliGenEpos::AliGenEpos(Int_t npart) : AliGenMC(npart),
		fBmin(0),
		fBmax(10000),
		fPhiMin(0),
		fPhiMax(2*3.1415926),
		fFilterModelOutput(kFALSE) {
	SetMC(new TEpos());
}

void AliGenEpos::Init() {
	AliGenMC::Init();
	TEpos *epos = GetTEpos();
	epos->SetLaproj(this->fZProjectile);
	epos->SetMaproj(this->fAProjectile);
	epos->SetLatarg(this->fZTarget);
	epos->SetMatarg(this->fATarget);
	epos->SetPhimin(this->fPhiMin);
	epos->SetPhimax(this->fPhiMax);
	epos->SetBminim(this->fBmin);
	epos->SetBmaxim(this->fBmax);
	epos->SetEcms(this->fEnergyCMS);
	GetTEpos()->Initialize();
}

void AliGenEpos::Generate() {
	  Float_t polar[3]   =   {0,0,0};
	  Float_t origin0[3]  =   {0,0,0};
	  Float_t origin[3]   =   {0,0,0};

	  Int_t nt  = 0; //output parameter for PushTrack

	  Vertex();
	  for (int j=0; j < 3; j++) origin0[j] = fVertex[j];

	  // Generate one event

	  GetTEpos()->GenerateEvent();
	  AliWarning("Generated");
	  GetTEpos()->ImportParticles(&fParticles);

	  Int_t np = fParticles.GetEntriesFast();
	  AliWarning(Form("Imported %d particles", np));

	  Int_t *idsOnStack = NULL;
	  idsOnStack = new Int_t[np];
	  TParticle *iparticle;

	  for (int i = 0; i < np; i++) {
		  iparticle = (TParticle *) fParticles.At(i);
		  //Bool_t isNullEntry = iparticle->GetStatusCode() == 0;
		  //Bool_t isCommentOrUnknown = iparticle->GetStatusCode() > 2;
		  Bool_t hasDecayed = iparticle->GetStatusCode() == 2;
		  Bool_t isFinalState = iparticle->GetStatusCode() == 1;
		  Int_t imo = iparticle->GetFirstMother();
		  Bool_t  hasMother = (imo >=0);


		  if (isFinalState) {
			  origin[0] = iparticle->Vx();
			  origin[1] = iparticle->Vy();
			  origin[2] = iparticle->Vz();
			  //doubled track with freeze out coordinates for femtoscopy
			  PushTrack(0,
					  imo>=0?idsOnStack[imo]:-1,
					  iparticle->GetPdgCode(),
				iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
				origin[0], origin[1], origin[2],
				iparticle->T(),
				polar[0],polar[1],polar[2],
				hasMother ? kPDecay:kPNoProcess,nt);

		      idsOnStack[i] = nt;
		      fNprimaries++;
		      KeepTrack(nt);

		      //real track with smeared vertex
		      origin[0] += origin0[0];
		      origin[1] += origin0[1];
		      origin[2] += origin0[2];
			  PushTrack(1,
					  nt,   //doubled track as mother
					  iparticle->GetPdgCode(),
				iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
				origin[0], origin[1], origin[2],
				iparticle->T(),
				polar[0],polar[1],polar[2],
				kPDecay,nt);
		      fNprimaries++;
		      KeepTrack(nt);
		  } else if(hasDecayed || !fFilterModelOutput) {
			  // don't track it and don't smear vertex
			  origin[0] = iparticle->Vx();
			  origin[1] = iparticle->Vy();
			  origin[2] = iparticle->Vz();
			  PushTrack(0,
					  imo>=0?idsOnStack[imo]:-1,
					  iparticle->GetPdgCode(),
				iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
				origin[0], origin[1], origin[2],
				iparticle->T(),
				polar[0],polar[1],polar[2],
				hasMother ? kPDecay:kPNoProcess,nt);
		      idsOnStack[i] = nt;
		      fNprimaries++;
		      KeepTrack(nt);
		  } else if(fFilterModelOutput) {
			  //filtered internal model objects
			  idsOnStack[i] = -1; // to erase mother field in its dauthers
		  }
	  }
	  SetHighWaterMark(fNprimaries);
	  TArrayF eventVertex;
	  eventVertex.Set(3);
	  eventVertex[0] = origin0[0];
	  eventVertex[1] = origin0[1];
	  eventVertex[2] = origin0[2];
	// Builds the event header, to be called after each event
	  AliGenEposEventHeader* header = new AliGenEposEventHeader("EPOS");

	  header->SetNProduced(fNprimaries);
	  header->SetPrimaryVertex(eventVertex);

	  header->SetImpactParameter(GetTEpos()->GetBimevt());
	  header->SetReactionPlaneAngle(GetTEpos()->GetPhievt());

	  header->SetHardScatters(0);
	  header->SetParticipants(GetTEpos()->GetNpjevt(), GetTEpos()->GetNtgevt());
	  header->SetCollisions(GetTEpos()->GetKolevt(), 0, 0, GetTEpos()->GetNpjevt());
	  header->SetSpectators(GetTEpos()->GetJpnevt(), GetTEpos()->GetJppevt(), GetTEpos()->GetJtnevt(), GetTEpos()->GetJtpevt());

	// Event Vertex
	  header->SetPrimaryVertex(fVertex);
	  AddHeader(header);
	  fCollisionGeometry = (AliGenEposEventHeader*)  header;

	  delete[] idsOnStack;
}

AliGenEpos::~AliGenEpos() {
}
Commit	Line	Data
9ef1c2d9	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