// //
////////////////////////////////////////////////////////////////////////////
+#include <iostream>
+#include <string>
+
#include "TUHKMgen.h"
#ifndef DATABASE_PDG
#include "DatabasePDG.h"
#include <TParticle.h>
#include <TClonesArray.h>
#include <TMCProcess.h>
-#include "TDatabasePDG.h"
-#include "TSystem.h"
+#include <TDatabasePDG.h>
+#include <TSystem.h>
#include "AliGenUHKM.h"
#include "AliRun.h"
#include "AliGenHijingEventHeader.h"
#include "AliLog.h"
-#include <iostream>
-#include <string>
using namespace std;
-//class TClonesArray;
-//enum TMCProcess;
-//class AliRun;
-//class TSystem;
-//class AliDecayer;
ClassImp(AliGenUHKM)
{
// Default constructor setting up default reasonable parameter values
// for central Pb+Pb collisions at 5.5TeV
- // cout << "AliGenUHKM::AliGenUHKM() IN" << endl;
// LHC
fHydjetParams.fSqrtS=5500; //LHC
fHydjetParams.fIenglu=0;
fHydjetParams.fIanglu=0;
*/
- strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
- strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
+ strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255);
+ strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255);
for(Int_t i=0; i<500; i++) {
fStableFlagPDG[i] = 0;
fStableFlagStatus[i] = kFALSE;
}
fStableFlagged = 0;
- // cout << "AliGenUHKM::AliGenUHKM() OUT" << endl;
}
//_______________________________________
// Constructor specifying the size of the particle table
// and setting up default reasonable parameter values
// for central Pb+Pb collisions at 5.5TeV
- // cout << "AliGenUHKM::AliGenUHKM(Int_t) IN" << endl;
+
fName = "UHKM";
fTitle= "Particle Generator using UHKM 3.0";
fHydjetParams.fIanglu=0;
*/
- strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
- strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
+ strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255);
+ strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255);
for(Int_t i=0; i<500; i++) {
fStableFlagPDG[i] = 0;
fStableFlagStatus[i] = kFALSE;
}
fStableFlagged = 0;
- // cout << "AliGenUHKM::AliGenUHKM(Int_t) OUT" << endl;
}
//__________________________________________
AliGenUHKM::~AliGenUHKM()
{
// Destructor, do nothing
+ // delete fParticles;
}
void AliGenUHKM::SetAllParametersRHIC()
// further to the model.
// HYDJET++ is initialized (average multiplicities are calculated, particle species definitions and decay
// channels are loaded, etc.)
- // cout << "AliGenUHKM::Init() IN" << endl;
SetMC(new TUHKMgen());
fUHKMgen = (TUHKMgen*) fMCEvGen;
SetAllParameters();
AliGenMC::Init();
-
fUHKMgen->Initialize();
- CheckPDGTable(); //
-
- fUHKMgen->Print(); //
- // cout << "AliGenUHKM::Init() OUT" << endl;
+ CheckPDGTable();
+ fUHKMgen->Print();
}
//________________________________________
{
// Generate one HYDJET++ event, get the output and push particles further
// to AliRoot's stack
- // cout << "AliGenUHKM::Generate() IN" << endl;
+
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
Float_t origin0[3] = {0,0,0};
+ Float_t time0 = 0.;
+ Float_t p[3];
Float_t v[3];
- Float_t mass, energy;
+ Float_t mass=0.0, energy=0.0;
Vertex();
for(Int_t j=0; j<3; j++) origin0[j] = fVertex[j];
+ time0 = fTime;
- fTrials = 0;
- // cout << "AliGenUHKM::Generate() fTrials = " << fTrials << endl;
-
- Int_t nt = 0;
-
+ // Generate the event and import particles
fUHKMgen->GenerateEvent();
- fTrials++;
-
fUHKMgen->ImportParticles(&fParticles,"All");
-
Int_t np = fParticles.GetEntriesFast();
- // cout << "AliGenUHKM::Generate() GetEntries " <<np<< endl;
+ Int_t nt = 0;
+ // Handle the IDs of particles on the stack
Int_t* idsOnStack = new Int_t[np];
Int_t* newPos = new Int_t[np];
- for(Int_t i=0; i<np; i++) newPos[i] = i;
-
- //_________ Loop for particle selection
- // for(Int_t i=1; i<np; i++) {
- for(Int_t i=1; i<np; i++) {
- // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- // the particle indexes are 0 based but fParticles is a 1 based array
- // -1 is the trivial code (when it does not exist)
- // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- TParticle *iparticle = (TParticle*)fParticles.At(i);
- // cout << "AliGenUHKM::Generate() particle #" << i << " in fParticles *********************"<< endl;
+ for(Int_t i=0; i<np; i++) {
+ newPos[i] = i;
+ idsOnStack[i] = -1;
+ }
+
- Int_t kf = iparticle->GetPdgCode();
- // if(kf==130 || kf==-130)
- // cout << "AliGenUHKM::Generate() PDG = " << kf << endl;
+ // Generate a random phi used to rotate the whole event
+ Double_t eventRotation = gRandom->Rndm()*TMath::Pi();
+ TParticle *iparticle;
+ Double_t partMomPhi=0.0, partPt=0.0;
+ Double_t partVtxPhi=0.0, partVtxR=0.0;
+ //_________ Loop for particles in the stack
+ for(Int_t i=0; i<np; i++) {
+ iparticle = (TParticle*)fParticles.At(i);
+ Int_t kf = iparticle->GetPdgCode();
Bool_t hasMother = (iparticle->GetFirstMother() >= 0);
-
- // cout << "AliGenUHKM::Generate() mother index in fParticles = "
- // << (iparticle->GetFirstMother()==-1 ? -1 : iparticle->GetFirstMother()+1)
- // << " ; hasMother = " << hasMother << endl;
- // cout << "AliGenUHKM::Generate() type (0-hydro;1-jet) = " << iparticle->GetStatusCode() << endl;
Bool_t hasDaughter = (iparticle->GetNDaughters() > 0);
-
- // cout << "AliGenUHKM::Generate() n.daughters = " << iparticle->GetNDaughters()
- //<< " ; hasDaughter = " << hasDaughter << endl;
-
-
+
if(hasDaughter) {
- // cout << "AliGenUHKM::Generate() decayed particle (not trackable)" << endl;
// This particle has decayed
// It will not be tracked
// Add it only once with coordinates not
// smeared with primary vertex position
- Float_t p[3] = {p[0] = iparticle->Px(),
- p[1] = iparticle->Py(),
- p[2] = iparticle->Pz()};
+
+ // rotate the direction of the particle
+ partMomPhi = TMath::ATan2(iparticle->Py(), iparticle->Px());
+ partPt = TMath::Hypot(iparticle->Px(), iparticle->Py());
+ p[0] = partPt*TMath::Cos(partMomPhi+eventRotation);
+ p[1] = partPt*TMath::Sin(partMomPhi+eventRotation);
+ p[2] = iparticle->Pz();
+
mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
- v[0] = iparticle->Vx();
- v[1] = iparticle->Vy();
- v[2] = iparticle->Vz();
- Float_t time = iparticle->T();
- // if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0
+ // rotate the freezeout point
+ partVtxPhi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
+ partVtxR = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
+ v[0] = partVtxR*TMath::Cos(partVtxPhi + eventRotation);
+ v[1] = partVtxR*TMath::Cos(partVtxPhi + eventRotation);
+ v[2] = iparticle->Vz();
+ Float_t time = iparticle->T();
+
Int_t imo = -1;
if(hasMother) {
- imo = iparticle->GetFirstMother(); //index of mother particle in fParticles
+ imo = iparticle->GetFirstMother(); //index of mother particle in fParticles
} // if has mother
Bool_t trackFlag = kFALSE; // tFlag is kFALSE --> do not track the particle
-
- // printf("Pushing Track %d with status %d mother %d\n", kf, trackFlag, imo>=0?idsOnStack[imo]:imo);
- PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
- p[0], p[1], p[2], energy,
- v[0], v[1], v[2], time,
- polar[0], polar[1], polar[2],
- (hasMother ? kPDecay : kPNoProcess), nt);
- // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
- // << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
+
+ PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf,
+ p[0], p[1], p[2], energy,
+ v[0], v[1], v[2], time,
+ polar[0], polar[1], polar[2],
+ (hasMother ? kPDecay : kPNoProcess), nt);
idsOnStack[i] = nt;
+
fNprimaries++;
KeepTrack(nt);
}
else {
- // cout << "AliGenUHKM::Generate() final particle --> push it twice on the stack" << endl;
// This is a final state particle
// It will be tracked
// Add it TWICE to the stack !!!
// this one will not be tracked
// Second time with event-wide c0ordinates and vertex smearing
// this one will be tracked
- Float_t p[3] = {p[0] = iparticle->Px(),
- p[1] = iparticle->Py(),
- p[2] = iparticle->Pz()};
- mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
+
+ // rotate the direction of the particle
+ partMomPhi = TMath::ATan2(iparticle->Py(), iparticle->Px());
+ partPt = TMath::Hypot(iparticle->Px(), iparticle->Py());
+ p[0] = partPt*TMath::Cos(partMomPhi+eventRotation);
+ p[1] = partPt*TMath::Sin(partMomPhi+eventRotation);
+ p[2] = iparticle->Pz();
+
energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
- v[0] = iparticle->Vx();
- v[1] = iparticle->Vy();
- v[2] = iparticle->Vz();
+ // rotate the freezeout point
+ partVtxPhi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
+ partVtxR = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
+ v[0] = partVtxR*TMath::Cos(partVtxPhi + eventRotation);
+ v[1] = partVtxR*TMath::Cos(partVtxPhi + eventRotation);
+ v[2] = iparticle->Vz();
+
Int_t type = iparticle->GetStatusCode(); // 1-from jet / 0-from hydro
- //if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0
Int_t coeffT=1;
if(type==1) coeffT=-1; //to separate particles from jets
-
+
+
Int_t imo = -1;
if(hasMother) {
- imo = iparticle->GetFirstMother();
+ imo = iparticle->GetFirstMother();
} // if has mother
+
Bool_t trackFlag = kFALSE; // tFlag = kFALSE --> do not track this one, its for femtoscopy
- PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf,
- p[0], p[1], p[2], energy,
- v[0], v[1], v[2], (iparticle->T())*coeffT,
- polar[0], polar[1], polar[2],
- hasMother ? kPDecay:kPNoProcess, nt);
- // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
- // << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl;
-
+ PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf,
+ p[0], p[1], p[2], energy,
+ v[0], v[1], v[2], (iparticle->T())*coeffT, // freeze-out time is negative if the particle comes from jet
+ polar[0], polar[1], polar[2],
+ hasMother ? kPDecay:kPNoProcess, nt);
+
idsOnStack[i] = nt;
fNprimaries++;
KeepTrack(nt);
-
+
origin[0] = origin0[0]+v[0];
origin[1] = origin0[1]+v[1];
origin[2] = origin0[2]+v[2];
+ Float_t time = time0+iparticle->T();
imo = nt;
- trackFlag = kTRUE; // tFlag = kTRUE --> track this one
- // cout << "AliGenUHKM::Generate() trackFlag = " << trackFlag << endl;
-
+ trackFlag = fTrackIt; // Track this particle, unless otherwise specified by fTrackIt
+
PushTrack(trackFlag, imo, kf,
- p[0], p[1], p[2], energy,
- origin[0], origin[1], origin[2], iparticle->T(),
- polar[0], polar[1], polar[2],
- hasMother ? kPDecay:kPNoProcess, nt);
- // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt
- // << "; mother index on stack = " << imo << endl;
+ p[0], p[1], p[2], energy,
+ origin[0], origin[1], origin[2], time,
+ polar[0], polar[1], polar[2],
+ hasMother ? kPDecay:kPNoProcess, nt);
+
fNprimaries++;
KeepTrack(nt);
}
}
-
+
SetHighWaterMark(fNprimaries);
TArrayF eventVertex;
eventVertex[0] = origin0[0];
eventVertex[1] = origin0[1];
eventVertex[2] = origin0[2];
+ Float_t eventTime = time0;
// Builds the event header, to be called after each event
AliGenEventHeader* header = new AliGenHijingEventHeader("UHKM");
+ Double_t b = 0.;
+ Double_t npart = 0;
+ Double_t nbin = 0;
+ fUHKMgen->GetCentrality(b, npart, nbin);
+ printf("********** %13.3f %13.3f %13.3f \n", b, npart, nbin);
+
+
((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
- ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
+ ((AliGenHijingEventHeader*) header)->SetInteractionTime(eventTime);
+ ((AliGenHijingEventHeader*) header)->SetImpactParameter(b);
((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
((AliGenHijingEventHeader*) header)->SetHardScatters(0);
- ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
- ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
+ ((AliGenHijingEventHeader*) header)->SetParticipants(Int_t(npart), 0);
+ ((AliGenHijingEventHeader*) header)->SetCollisions(Int_t(nbin), 0, 0, 0);
((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(0);//evrot);
header->SetPrimaryVertex(fVertex);
+ header->SetInteractionTime(fTime);
AddHeader(header);
fCollisionGeometry = (AliGenHijingEventHeader*) header;
- delete idsOnStack;
-
- // gAlice->SetGenEventHeader(header);
-
- // printf(" Finish Generate .. %d ..\n",nt);
- // cout << "AliGenUHKM::Generate() OUT" << endl;
+ delete [] idsOnStack;
+ delete [] newPos;
}
void AliGenUHKM::Copy(TObject &) const
void AliGenUHKM::SetAllParameters() {
// Forward all input parameters to the TGenerator::TUHKMgen object
- // cout << "AliGenUHKM::SetAllParameters() IN" << endl;
-
fUHKMgen->SetEcms(fHydjetParams.fSqrtS);
fUHKMgen->SetBmin(fHydjetParams.fBmin);
fUHKMgen->SetBmax(fHydjetParams.fBmax);
fUHKMgen->SetCoordAsymmPar(fHydjetParams.fEpsilon);
fUHKMgen->SetGammaS(fHydjetParams.fCorrS);
- // fUHKMgen->SetHdec(fHydjetParams.fTime);
fUHKMgen->SetEtaType(fHydjetParams.fEtaType);
fUHKMgen->SetFlagWeakDecay(fHydjetParams.fWeakDecay);
// fUHKMgen->SetMinimumMass(fMinMass);
// fUHKMgen->SetMaximumMass(fMaxMass);
- cout << "AliGenUHKM::Init() no. stable flagged particles = " << fStableFlagged << endl;
+ // cout << "AliGenUHKM::Init() no. stable flagged particles = " << fStableFlagged << endl;
for(Int_t i=0; i<fStableFlagged; i++) {
- cout << "AliGenUHKM::Init() flag no. " << i
- << " PDG = " << fStableFlagPDG[i]
- << " flag = " << fStableFlagStatus[i] << endl;
+ // cout << "AliGenUHKM::Init() flag no. " << i
+ // << " PDG = " << fStableFlagPDG[i]
+ // << " flag = " << fStableFlagStatus[i] << endl;
fUHKMgen->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
}
- cout<<" Print all parameters "<<endl;
- cout<<" SqrtS = "<<fHydjetParams.fSqrtS<<endl;
- cout<<" Bmin = "<< fHydjetParams.fBmin<<endl;
+ cout<<" Print all parameters "<<endl;
+ cout<<" SqrtS = "<<fHydjetParams.fSqrtS<<endl;
+ cout<<" Bmin = "<< fHydjetParams.fBmin<<endl;
cout<<" Bmax= "<<fHydjetParams.fBmax<<endl;
cout<<" Aw= "<<fHydjetParams.fAw<<endl;
cout<<" Seed= "<<fHydjetParams.fSeed<<endl;
cout<<" --------Flags------ "<<endl;
cout<<" GammaS= "<<fHydjetParams.fCorrS<<endl;
- // fUHKMgen->SetHdec(fHydjetParams.fTime<<endl;
cout<<" EtaType= "<<fHydjetParams.fEtaType<<endl;
cout<<" FlagWeakDecay= "<<fHydjetParams.fWeakDecay<<endl;
// Add temporarely all particle definitions from HYDJET++ which miss in the ROOT's PDG tables
// to the TDatabasePDG table.
- // cout << "AliGenUHKM::CheckPDGTable() IN" << endl;
- //TDabasePDG *rootPDG = TDatabasePDG::Instance(); // ROOT's PDG table
DatabasePDG *uhkmPDG = fUHKMgen->PDGInfo(); // UHKM's PDG table
TParticlePDG *rootTestParticle;
ParticlePDG *uhkmTestParticle;
- // cout << "particles with good status in UHKM table = " << uhkmPDG->GetNParticles() << endl;
// loop over all particles in the SHARE table
for(Int_t i=0; i<uhkmPDG->GetNParticles(); i++) {
- // cout << "particle #" << i << " ================" << endl;
// get a particle specie
uhkmTestParticle = uhkmPDG->GetPDGParticleByIndex(i);
- // cout << "PDG = " << uhkmTestParticle->GetPDG() << endl;
// check if this code exists in ROOT's table
rootTestParticle = TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG());
if(!rootTestParticle) { // if not then add it to the ROOT's PDG database
uhkmTestParticle->GetWidth(),
(Int_t(uhkmTestParticle->GetBaryonNumber())==0 ? "meson" : "baryon"),
uhkmTestParticle->GetPDG());
- // cout << "Not found in ROOT's PDG database --> added now" << endl;
if(uhkmTestParticle->GetWidth()<1e-10)
- cout << uhkmTestParticle->GetPDG() << "its mass "
+ cout << uhkmTestParticle->GetPDG() << " with mass "
<< TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Mass()
<< TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Width() << endl;
}
- // else
- // cout << "Found in ROOT's PDG database --> not added" << endl;
} // end for
- // cout << "AliGenUHKM::CheckPDGTable() OUT" << endl;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff