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

//
// AliGenEpos.cxx
//
//  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 "TMath.h"
#include "AliLog.h"
#include "AliGenEventHeader.h"
#include "AliGenEposEventHeader.h"

ClassImp(AliGenEpos)

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

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

void AliGenEpos::Init() {
  // Sets up TEpos
	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);
	epos->SetSplitting(kTRUE);
	GetTEpos()->Initialize();
}

void AliGenEpos::Generate() {
  // Does actual generation and output conversion
	  Float_t polar[3]    =   {0,0,0};
	  Float_t origin0[3]  =   {0,0,0};
	  Float_t origin[3]   =   {0,0,0};
	  Float_t time0 = 0.;
	  Float_t time  = 0.;
	  fNprimaries = 0;
	  Int_t nt  = 0; //output parameter for PushTrack

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

	  // 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];
	  for (int i = 0; i < np; i++) idsOnStack[i] = 0;
	  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();
			  time      = iparticle->T();
			  //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], time,
				polar[0],polar[1],polar[2],
				hasMother ? kPDecay:kPNoProcess, nt, 1., 
				iparticle->GetStatusCode());

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

		      //real track with smeared vertex
		      origin[0] += origin0[0];
		      origin[1] += origin0[1];
		      origin[2] += origin0[2];
		      time      += time0;
			  PushTrack(1,
					  nt,   //doubled track as mother
					  iparticle->GetPdgCode(),
				iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
				origin[0], origin[1], origin[2],
				time,
				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();
			  time      = iparticle->T();
			  PushTrack(0, imo>=0?idsOnStack[imo]:-1,
				iparticle->GetPdgCode(),
				iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
				origin[0], origin[1], origin[2], time,
				polar[0],polar[1],polar[2],
				hasMother ? kPDecay:kPNoProcess, nt, 1., 
				iparticle->GetStatusCode());
		      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->SetInteractionTime(time0);

	  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);
	  
	  header->SetBimevt(GetTEpos()->GetBimevt());
	  header->SetPhievt(GetTEpos()->GetPhievt());
	  header->SetKolevt(GetTEpos()->GetKolevt());
	  header->SetKoievt(GetTEpos()->GetKoievt());
	  header->SetPmxevt(GetTEpos()->GetPmxevt());
	  header->SetEgyevt(GetTEpos()->GetEgyevt());
	  header->SetNpjevt(GetTEpos()->GetNpjevt());
	  header->SetNtgevt(GetTEpos()->GetNtgevt());
	  header->SetNpnevt(GetTEpos()->GetNpnevt());
	  header->SetNppevt(GetTEpos()->GetNppevt());
	  header->SetNtnevt(GetTEpos()->GetNtnevt());
	  header->SetNtpevt(GetTEpos()->GetNtpevt());
	  header->SetJpnevt(GetTEpos()->GetJpnevt());
	  header->SetJppevt(GetTEpos()->GetJppevt());
	  header->SetJtnevt(GetTEpos()->GetJtnevt());
	  header->SetJtpevt(GetTEpos()->GetJtpevt());
	  header->SetXbjevt(GetTEpos()->GetXbjevt());
	  header->SetQsqevt(GetTEpos()->GetQsqevt());
	  header->SetNglevt(GetTEpos()->GetNglevt());
	  header->SetZppevt(GetTEpos()->GetZppevt());
	  header->SetZptevt(GetTEpos()->GetZptevt());
    AddHeader(header);
    fCollisionGeometry = (AliGenEposEventHeader*)  header;
    
    delete[] idsOnStack;
}

AliGenEpos::~AliGenEpos() {
}
Commit	Line	Data
9033ee49	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	//
9ef1c2d9	9
	10	#include "AliGenEpos.h"
	11	#include "TEpos.h"
	12	#include "TParticle.h"
879c2333	13	#include "TMath.h"
9ef1c2d9	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),
879c2333	24	fPhiMax(TMath::TwoPi()),
9ef1c2d9	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),
879c2333	33	fPhiMax(TMath::TwoPi()),
9ef1c2d9	34	fFilterModelOutput(kFALSE) {
	35	SetMC(new TEpos());
	36	}
	37
	38	void AliGenEpos::Init() {
9033ee49	39	// Sets up TEpos
9ef1c2d9	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);
06c8cc00	51	epos->SetSplitting(kTRUE);
9ef1c2d9	52	GetTEpos()->Initialize();
	53	}
	54
	55	void AliGenEpos::Generate() {
9033ee49	56	// Does actual generation and output conversion
e6348652	57	Float_t polar[3] = {0,0,0};
9ef1c2d9	58	Float_t origin0[3] = {0,0,0};
9ef1c2d9	59	Float_t origin[3] = {0,0,0};
21391258	60	Float_t time0 = 0.;
21391258	61	Float_t time = 0.;
879c2333	62	fNprimaries = 0;
9ef1c2d9	63	Int_t nt = 0; //output parameter for PushTrack
	64
	65	Vertex();
	66	for (int j=0; j < 3; j++) origin0[j] = fVertex[j];
21391258	67	time0 = fTime;
9ef1c2d9	68
	69	// Generate one event
	70
	71	GetTEpos()->GenerateEvent();
079fb251	72	//AliWarning("Generated");
9ef1c2d9	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];
e6348652	80	for (int i = 0; i < np; i++) idsOnStack[i] = 0;
9ef1c2d9	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;
879c2333	87	Bool_t hasDecayed = iparticle->GetStatusCode() >= 2;
9ef1c2d9	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();
21391258	97	time = iparticle->T();
9ef1c2d9	98	//doubled track with freeze out coordinates for femtoscopy
079fb251	99	PushTrack(0, imo>=0?idsOnStack[imo]:-1,
079fb251	100	iparticle->GetPdgCode(),
9ef1c2d9	101	iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
079fb251	102	origin[0], origin[1], origin[2], time,
9ef1c2d9	103	polar[0],polar[1],polar[2],
079fb251	104	hasMother ? kPDecay:kPNoProcess, nt, 1.,
079fb251	105	iparticle->GetStatusCode());
9ef1c2d9	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];
21391258	115	time += time0;
9ef1c2d9	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],
21391258	121	time,
9ef1c2d9	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();
21391258	131	time = iparticle->T();
079fb251	132	PushTrack(0, imo>=0?idsOnStack[imo]:-1,
079fb251	133	iparticle->GetPdgCode(),
9ef1c2d9	134	iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
079fb251	135	origin[0], origin[1], origin[2], time,
9ef1c2d9	136	polar[0],polar[1],polar[2],
079fb251	137	hasMother ? kPDecay:kPNoProcess, nt, 1.,
079fb251	138	iparticle->GetStatusCode());
9ef1c2d9	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);
21391258	158	header->SetInteractionTime(time0);
9ef1c2d9	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);
33652044	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;
9ef1c2d9	196	}
	197
	198	AliGenEpos::~AliGenEpos() {
	199	}
	200