#include <Riostream.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TF2.h>
#include <TH1.h>
#include <TH2.h>
#include <TObjArray.h>
#include <TPDGCode.h>
#include <TParticle.h>
#include <TDatabasePDG.h>
#include <TROOT.h>
#include "AliGeVSimParticle.h"
#include "AliGenGeVSim.h"
#include "AliGenGeVSimEventHeader.h"
#include "AliGenerator.h"
#include "AliRun.h"
using std::cout;
using std::endl;
ClassImp(AliGenGeVSim)
AliGenGeVSim::AliGenGeVSim() :
AliGenerator(-1),
fModel(0),
fPsi(0),
fIsMultTotal(kTRUE),
fPtFormula(0),
fYFormula(0),
fPhiFormula(0),
fCurrentForm(0),
fPtYHist(0),
fPartTypes(0)
{
for (Int_t i=0; i<4; i++)
fPtYFormula[i] = 0;
for (Int_t i=0; i<2; i++)
fHist[i] = 0;
}
AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal)
: AliGenerator(-1),
fModel(0),
fPsi(psi),
fIsMultTotal(isMultTotal),
fPtFormula(0),
fYFormula(0),
fPhiFormula(0),
fCurrentForm(0),
fPtYHist(0),
fPartTypes(0)
{
if (psi < 0 || psi > 360 )
Error ("AliGenGeVSim", "Reaction plane angle ( %13.3f )out of range [0..360]", psi);
fPsi = psi * TMath::Pi() / 180. ;
fIsMultTotal = isMultTotal;
fPartTypes = new TObjArray();
InitFormula();
}
AliGenGeVSim::~AliGenGeVSim() {
if (fPartTypes != NULL) delete fPartTypes;
}
Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) const {
if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
return kTRUE;
}
Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
if ( TestBit(kThetaRange) ) {
Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]);
if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
}
if ( TestBit(kMomentumRange) ) {
Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
if ( momentum < fPMin || momentum > fPMax) return kFALSE;
}
return kTRUE;
}
static Double_t aPtForm(Double_t * x, Double_t * par) {
return x[0] * TMath::Exp( -sqrt(par[0]*par[0] + x[0]*x[0]) / par[1]);
}
static Double_t aYForm(Double_t * x, Double_t * par) {
return TMath::Exp ( - x[0]*x[0] / (2 * par[0]*par[0] ) );
}
static Double_t aPtYFormula0(Double_t *x, Double_t * par) {
Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
return x[0] * aFormE * TMath::Exp(-aFormE/par[1]);
}
static Double_t aPtYFormula1(Double_t *x, Double_t * par) {
Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
return x[0] * aFormE / ( TMath::Exp( aFormE / par[1]) - 1 );
}
static Double_t aPtYFormula2(Double_t *x, Double_t * par) {
Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
Double_t aFormG = 1 / TMath::Sqrt((1.-par[2])*(1.+par[2]));
Double_t aFormYp = par[2]*TMath::Sqrt( (par[0]*par[0] + x[0]*x[0])
* (TMath::CosH(x[1])-par[0])*(TMath::CosH(x[1])+par[0]))
/( par[1]*TMath::Sqrt((1.-par[2])*(1.+par[2])));
return x[0] * aFormE * TMath::Exp( - aFormG * aFormE / par[1])
*( TMath::SinH(aFormYp)/aFormYp
+ par[1]/(aFormG*aFormE)
* ( TMath::SinH(aFormYp)/aFormYp-TMath::CosH(aFormYp) ) );
}
static Double_t aPhiForm(Double_t * x, Double_t * par) {
return 1 + 2*par[1]*TMath::Cos(x[0]-par[0])
+ 2*par[2]*TMath::Cos(2*(x[0]-par[0]));
}
void AliGenGeVSim::InitFormula() {
fPtFormula = new TF1("gevsimPt", &aPtForm, 0, 3, 2);
fYFormula = new TF1("gevsimRapidity", &aYForm, -3, 3,1);
fPtFormula->SetParNames("mass", "temperature");
fPtFormula->SetParameters(1., 1.);
fYFormula->SetParName(0, "sigmaY");
fYFormula->SetParameter(0, 1.);
fPtFormula->SetNpx(100);
fYFormula->SetNpx(100);
fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2);
fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2);
fPtYFormula[2] = new TF2("gevsimPtY_4", &aPtYFormula2, 0, 3, -2, 2, 3);
fPtYFormula[3] = 0;
const char* kParamNames[3] = {"mass", "temperature", "expVel"};
Int_t i, j;
for (i=0; i<3; i++) {
fPtYFormula[i]->SetNpx(100);
fPtYFormula[i]->SetNpy(100);
for (j=0; j<3; j++) {
if ( i != 2 && j == 2 ) continue;
fPtYFormula[i]->SetParName(j, kParamNames[j]);
fPtYFormula[i]->SetParameter(j, 0.5);
}
}
fPhiFormula = new TF1("gevsimPhi", &aPhiForm, 0, 2*TMath::Pi(), 3);
fPhiFormula->SetParNames("psi", "directed", "elliptic");
fPhiFormula->SetParameters(0., 0., 0.);
fPhiFormula->SetNpx(180);
}
void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
if (fPartTypes == NULL)
fPartTypes = new TObjArray();
fPartTypes->Add(part);
}
void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
fIsMultTotal = isTotal;
}
Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
static const char* ending[] = {"", "Rndm"};
static const char* patt1 = "gevsim%s%s";
static const char* patt2 = "gevsim%d%s%s";
char buffer[80];
TF1 *form;
Float_t scaler = 1.;
Int_t i, j;
for (i=0; i<2; i++) {
for (j=0; j<2; j++) {
form = 0;
if (i == 0) snprintf(buffer, 80, patt1, params[paramId], ending[j]);
else snprintf(buffer, 80, patt2, pdg, params[paramId], ending[j]);
form = (TF1 *)gROOT->GetFunction(buffer);
if (form != 0) {
if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber());
if (j == 1) {
form->SetParameter(0, gAlice->GetEvNumber());
scaler *= form->GetRandom();
}
}
}
}
return scaler;
}
void AliGenGeVSim::DetermineReactionPlane() {
TF1 *form;
form = 0;
form = (TF1 *)gROOT->GetFunction("gevsimPsi");
if (form) fPsi = form->Eval(gAlice->GetEvNumber()) * TMath::Pi() / 180;
form = 0;
form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
if (form) fPsi = form->GetRandom() * TMath::Pi() / 180;
cout << "Psi = " << fPsi << "\t" << (Int_t)(fPsi*180./TMath::Pi()) << endl;
fPhiFormula->SetParameter(0, fPsi);
}
Float_t AliGenGeVSim::GetdNdYToTotal() {
Float_t integ, mag;
const Double_t kMaxPt = 3.0, kMaxY = 2.;
if (fModel == 1) {
integ = fYFormula->Integral(-kMaxY, kMaxY);
mag = fYFormula->Eval(0);
return integ/mag;
}
if (fModel > 1 && fModel < 6) {
integ = ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY);
mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -0.1, 0.1) / 0.2;
return integ/mag;
}
if (fModel == 6) {
integ = fHist[1]->Integral();
mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
return integ/mag;
}
if (fModel == 7) {
Int_t yBins = fPtYHist->GetNbinsY();
Int_t ptBins = fPtYHist->GetNbinsX();
integ = fPtYHist->Integral(0, ptBins, 0, yBins);
mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
return integ/mag;
}
return 1;
}
void AliGenGeVSim::SetFormula(Int_t pdg) {
char buff[40];
const char* msg[4] = {
"Custom Formula for Pt Y distribution not found [pdg = %d]",
"Histogram for Pt distribution not found [pdg = %d]",
"Histogram for Y distribution not found [pdg = %d]",
"HIstogram for Pt Y dostribution not found [pdg = %d]"
};
const char* pattern[8] = {
"gevsimDistPtY", "gevsimDist%dPtY",
"gevsimHistPt", "gevsimHist%dPt",
"gevsimHistY", "gevsimHist%dY",
"gevsimHistPtY", "gevsimHist%dPtY"
};
const char *where = "SetFormula";
if (fModel < 1 || fModel > 7)
Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
if (fModel == 1) fCurrentForm = fPtFormula;
if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
if (fModel == 5) {
fCurrentForm = 0;
fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
if (!fCurrentForm) {
snprintf(buff, 40, pattern[1], pdg);
fCurrentForm = (TF2*)gROOT->GetFunction(buff);
if (!fCurrentForm) Error(where, msg[0], pdg);
}
}
if (fModel == 6) {
for (Int_t i=0; i<2; i++) {
fHist[i] = 0;
fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
if (!fHist[i]) {
snprintf(buff, 40, pattern[3+2*i], pdg);
fHist[i] = (TH1D*)gROOT->FindObject(buff);
if (!fHist[i]) Error(where, msg[1+i], pdg);
}
}
}
if (fModel == 7) {
fPtYHist = 0;
fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
if (!fPtYHist) {
snprintf(buff, 40, pattern[7], pdg);
fPtYHist = (TH2D*)gROOT->FindObject(buff);
}
if (!fPtYHist) Error(where, msg[3], pdg);
}
}
void AliGenGeVSim:: AdjustFormula() {
const Double_t kMaxPt = 3.0;
const Double_t kMaxY = 2.0;
Double_t minPt, maxPt, minY, maxY;
if (fModel > 4) return;
if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
else maxPt = kMaxPt;
if (TestBit(kPtRange)) minPt = fPtMin;
else minPt = 0;
if (TestBit(kPtRange) && fPtMin > kMaxPt )
Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
if (TestBit(kYRange)) {
minY = fYMin;
maxY = fYMax;
} else {
minY = -kMaxY;
maxY = kMaxY;
}
if (fModel == 1) {
fPtFormula->SetRange(fPtMin, maxPt);
fYFormula->SetRange(fYMin, fYMax);
}
if (fModel > 1)
((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
if (TestBit(kPhiRange))
fPhiFormula->SetRange(fPhiMin, fPhiMax);
}
void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
if (fModel == 1) {
pt = fPtFormula->GetRandom();
y = fYFormula->GetRandom();
return;
}
if (fModel > 1 && fModel < 6) {
((TF2*)fCurrentForm)->GetRandom2(pt, y);
return;
}
if (fModel == 6) {
pt = fHist[0]->GetRandom();
y = fHist[1]->GetRandom();
}
if (fModel == 7) {
fPtYHist->GetRandom2(pt, y);
return;
}
}
void AliGenGeVSim::Init() {
}
void AliGenGeVSim::Generate() {
PDG_t pdg;
Float_t mass;
Float_t orgin[3] = {0,0,0};
Float_t polar[3] = {0,0,0};
Float_t time = 0;
Float_t multiplicity = 0;
Bool_t isMultTotal = kTRUE;
Float_t paramScaler;
Float_t directedScaller = 1., ellipticScaller = 1.;
TLorentzVector *v = new TLorentzVector(0,0,0,0);
const Int_t kParent = -1;
Int_t id;
VertexInternal();
orgin[0] = fVertex[0];
orgin[1] = fVertex[1];
orgin[2] = fVertex[2];
time = fTime;
TDatabasePDG *db = TDatabasePDG::Instance();
TParticlePDG *type;
AliGeVSimParticle *partType;
Int_t nType, nParticle, nParam;
const Int_t kNParams = 6;
DetermineReactionPlane();
for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
partType = (AliGeVSimParticle *)fPartTypes->At(nType);
pdg = (PDG_t)partType->GetPdgCode();
type = db->GetParticle(pdg);
mass = type->Mass();
fModel = partType->GetModel();
SetFormula(pdg);
fCurrentForm->SetParameter("mass", mass);
for (nParam = 0; nParam < kNParams; nParam++) {
paramScaler = FindScaler(nParam, pdg);
if (nParam == 0)
fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
if (nParam == 1 && fModel == 1)
fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
if (nParam == 2 && fModel == 4) {
Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity();
if (totalExpVal == 0.0) totalExpVal = 0.0001;
if (totalExpVal == 1.0) totalExpVal = 9.9999;
fCurrentForm->SetParameter("expVel", totalExpVal);
}
if (nParam == 3) directedScaller = paramScaler;
if (nParam == 4) ellipticScaller = paramScaler;
if (nParam == 5) {
if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
else isMultTotal = fIsMultTotal;
multiplicity = paramScaler * partType->GetMultiplicity();
multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
}
}
if (partType->IsFlowSimple()) {
fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
}
AdjustFormula();
Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
nParticle = 0;
while (nParticle < multiplicity) {
Double_t pt, y, phi;
Float_t p[3] = {0,0,0};
GetRandomPtY(pt, y);
if (!partType->IsFlowSimple()) {
fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
}
phi = fPhiFormula->GetRandom();
if (!isMultTotal) nParticle++;
if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
v->SetPtEtaPhiM(pt, y, phi, mass);
p[0] = v->Px();
p[1] = v->Py();
p[2] = v->Pz();
if ( !CheckAcceptance(p) ) continue;
PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);
if (isMultTotal) nParticle++;
}
}
AliGenGeVSimEventHeader *header = new AliGenGeVSimEventHeader("GeVSim header");
TArrayF eventVertex(3,orgin);
header->SetPrimaryVertex(eventVertex);
header->SetInteractionTime(time);
header->SetEventPlane(fPsi);
header->SetEllipticFlow(fPhiFormula->GetParameter(2));
gAlice->SetGenEventHeader(header);
delete v;
}