#include "Riostream.h"
#include "TParticle.h"
#include "AliStack.h"
#include "AliRun.h"
#include "TMath.h"
#include "TMCProcess.h"
#include "TList.h"
#include "TVector3.h"
#include "AliMC.h"
#include "TArrayF.h"
#include "AliCollisionGeometry.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenCocktailEntry.h"
#include "AliGenCocktailAfterBurner.h"
#include "AliGenDeuteron.h"
#include "AliLog.h"
ClassImp(AliGenDeuteron)
AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
:fSign(1)
,fPmax(pmax)
,fRmax(rmax)
,fSpinProb(1)
,fClusterType(cluster)
,fModel(0)
,fTimeLength(2.5)
,fB(0)
,fR(0)
,fPsiR(0)
,fCurStack(0)
{
fSign = sign > 0 ? 1:-1;
}
AliGenDeuteron::~AliGenDeuteron()
{
}
void AliGenDeuteron::Init()
{
}
void AliGenDeuteron::Generate()
{
AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
Bool_t collisionGeometry=0;
if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
{
TString name;
Int_t ap, zp, at, zt;
gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
if(ap != at) AliWarning("projectile and target have different size");
fR = 1.29*TMath::Power(at, 1./3.);
collisionGeometry = 1;
}
for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
{
gener->SetActiveEventNumber(ns);
if(fModel != kNone && collisionGeometry)
{
fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
if(fB >= 2.*fR) continue;
}
fCurStack = gener->GetStack(ns);
if(!fCurStack)
{
AliWarning("no event stack");
return;
}
TList* protons = new TList();
protons->SetOwner(kFALSE);
TList* neutrons = new TList();
neutrons->SetOwner(kFALSE);
for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
{
TParticle* iParticle = fCurStack->Particle(i);
if(iParticle->GetStatusCode() != 1) continue;
Int_t pdgCode = iParticle->GetPdgCode();
if(pdgCode == fSign*2212)
{
FixProductionVertex(iParticle);
protons->Add(iParticle);
}
else if(pdgCode == fSign*2112)
{
FixProductionVertex(iParticle);
neutrons->Add(iParticle);
}
}
if(fClusterType==kFirstPartner)
{
FirstPartner(protons, neutrons);
}
else
{
WeightMatrix(protons, neutrons);
}
protons->Clear("nodelete");
neutrons->Clear("nodelete");
delete protons;
delete neutrons;
}
AliInfo("DONE");
}
Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
{
const Double_t kProtonMass = 0.938272013;
const Double_t kNeutronMass = 0.939565378;
TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass);
if( deltaP >= fPmax) return -1.;
Double_t deltaR = (v2-v1).Mag();
if(deltaR >= fRmax*1.e-13) return -1.;
if(Rndm() > fSpinProb) return -1.;
if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
}
void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
{
TIter p_next(protons);
while(TParticle* n0 = (TParticle*) p_next())
{
TParticle* partner = 0;
TIter n_next(neutrons);
while(TParticle* n1 = (TParticle*) n_next() )
{
if(GetCoalescenceProbability(n0, n1) < 0 ) continue;
partner = n1;
break;
}
if(partner == 0) continue;
PushDeuteron(n0, partner);
neutrons->Remove(partner);
}
}
void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
{
Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
TParticle** cProton = new TParticle*[nMaxPairs];
TParticle** cNeutron = new TParticle*[nMaxPairs];
Double_t* cWeight = new Double_t[nMaxPairs];
Int_t cIdx = -1;
TIter p_next(protons);
while(TParticle* n1 = (TParticle*) p_next())
{
TIter n_next(neutrons);
while(TParticle* n2 = (TParticle*) n_next() )
{
Double_t weight = this->GetCoalescenceProbability(n1,n2);
if(weight < 0) continue;
++cIdx;
cProton[cIdx] = n1;
cNeutron[cIdx] = n2;
cWeight[cIdx] = weight;
}
n_next.Reset();
}
p_next.Reset();
Int_t nPairs = cIdx + 1;
Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
TParticle** iProton = new TParticle*[nMaxIntPair];
TParticle** iNeutron = new TParticle*[nMaxIntPair];
Double_t* iWeight = new Double_t[nMaxIntPair];
Int_t iIdx = -1;
while(kTRUE)
{
Int_t j = -1;
Double_t wMax = 0;
for(Int_t i=0; i < nPairs; ++i)
{
if(cWeight[i] > wMax)
{
wMax=cWeight[i];
j = i;
}
}
if(j == -1 ) break;
++iIdx;
iProton[iIdx] = cProton[j];
iNeutron[iIdx] = cNeutron[j];
iWeight[iIdx] = cWeight[j];
for(Int_t i=0; i < nPairs; ++i)
{
if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
}
}
Int_t nIntPairs = iIdx + 1;
delete[] cProton;
delete[] cNeutron;
delete[] cWeight;
for(Int_t i=0; i<nIntPairs; ++i)
{
TParticle* n1 = iProton[i];
TParticle* n2 = iNeutron[i];
PushDeuteron(n1,n2);
}
delete[] iProton;
delete[] iNeutron;
delete[] iWeight;
}
void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
{
const Double_t kDeuteronMass = 1.87561282;
const Int_t kDeuteronPdg = 1000010020;
TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
TVector3 pN = p1+p2;
TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
Int_t ntrk = 0;
Double_t weight = 1;
Int_t is = 1;
fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
pN.X(), pN.Y(), pN.Z(), energy,
vN.X(), vN.Y(), vN.Z(), parent1->T(),
0., 0., 0., kPNCapture, ntrk, weight, is);
parent1->SetStatusCode(kCluster);
parent2->SetStatusCode(kCluster);
parent1->SetBit(kDoneBit);
parent2->SetBit(kDoneBit);
}
void AliGenDeuteron::FixProductionVertex(TParticle* i)
{
if(fModel == kNone || fModel > kExpansion) return;
Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.);
Double_t b = fTimeLength + fR - fB/2.;
Double_t c = fTimeLength;
Double_t xx = 0;
Double_t yy = 0;
Double_t zz = 0;
if(fModel == kThermal)
{
Double_t r = TMath::Power(Rndm(),1./3.);
Double_t theta = TMath::ACos(2.*Rndm()-1.);
Double_t phi = 2.*TMath::Pi()*Rndm();
xx = a*r*TMath::Sin(theta)*TMath::Cos(phi);
yy = b*r*TMath::Sin(theta)*TMath::Sin(phi);
zz = c*r*TMath::Cos(theta);
}
else if(fModel == kExpansion)
{
xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());
yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());
zz = c*TMath::Cos(i->Theta());
}
Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR);
Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR);
Double_t z = zz;
i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
}
Double_t AliGenDeuteron::GetS(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
{
Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
}
Double_t AliGenDeuteron::GetPcm(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
{
Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
}
Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
{
return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
}