New class for deuteron generation

Eulogio Serradilla <eulogio.serradilla@ciemat.es>

diff --git a/EVGEN/AliGenDeuteron.cxx b/EVGEN/AliGenDeuteron.cxx
new file mode 100644 (file)
index 0000000..2dd4eb5
--- /dev/null
+++ b/EVGEN/AliGenDeuteron.cxx
@@ -0,0 +1,214 @@
+#include "Riostream.h"\r
+#include "TParticle.h"\r
+#include "AliStack.h"\r
+#include "AliGenDeuteron.h"\r
+#include "AliGenCocktailAfterBurner.h"\r
+#include "TMath.h"\r
+#include "TMCProcess.h"\r
+#include "TList.h"\r
+#include "TVector3.h"\r
+#include "TRandom3.h"\r
+#include "AliMC.h"\r
+/*\r
+A very simple model based on the article "Physics Letters B325 (1994) 7-12". \r
+The only addition is to model the freeze-out at which the\r
+light nuclei can form for the energies of the LHC, since PYTHIA does not give\r
+any detailed spatial description of the generated particles at the level of a\r
+few fermis.\r
+\r
+There are two options to place the nucleons: a thermal picture where all\r
+nucleons are placed randomly and homogeneously in an spherical volume of a\r
+few fermis, and expansion one where they are projected on its surface.\r
+\r
+A (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron\r
+with momentum difference less than ~ 300MeV and relative distance less than a\r
+~ 2.1fm. Only 3/4 of this clusters are accepted due to the deuteron spin.\r
+*/\r
+\r
+ClassImp(AliGenDeuteron)\r
+\r
+\r
+AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
+fDeuteronMass(1.87561282), //pdg\r
+fPmax(pmax), // GeV/c\r
+fRmax(rmax*1.e-13), //cm\r
+fRsrc(rsrc*1.e-13), // cm\r
+fSpinProb(0.75),\r
+fModel(model),\r
+fSign(1)\r
+{\r
+       fSign = sign > 0 ? 1:-1;\r
+}\r
+\r
+AliGenDeuteron::~AliGenDeuteron()\r
+{\r
+}\r
+\r
+void AliGenDeuteron::Init()\r
+{\r
+       // Standard AliGenerator Initializer\r
+}\r
+\r
+void AliGenDeuteron::Generate()\r
+{\r
+// Modify the stack of each event, adding the new particles at the end and\r
+// not transporting the parents.\r
+\r
+       Info("Generate","freeze-out model : %d (0 expansion, 1 thermal)",fModel);\r
+       Info("Generate","relative momentum: %g GeV/c",fPmax);\r
+       Info("Generate","relative distance: %g cm",fRmax);\r
+       Info("Generate","source radius    : %g cm",fRsrc);\r
+       Info("Generate","spin probability : %g ",fSpinProb);\r
+       Info("Generate","sign             : %d ",fSign);\r
+       \r
+       TRandom3 rnd(0);\r
+       TRandom3 rnd2(0);\r
+       \r
+       Int_t ntr;\r
+       \r
+       // Get the cocktail generator\r
+       AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();\r
+       \r
+       // Loop over events\r
+       for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)\r
+       {\r
+               gener->SetActiveEventNumber(ns);\r
+               \r
+               AliStack* stack = gener->GetStack(ns); // Stack of event ns\r
+               if(!stack)\r
+               {\r
+                       Info("Generate", "no stack, exiting");\r
+                       return;\r
+               }\r
+               \r
+               // find emerging nucleons from the collision of current event\r
+               \r
+               TList* protons = new TList();\r
+               protons->SetOwner(kFALSE);\r
+               TList* neutrons = new TList();\r
+               neutrons->SetOwner(kFALSE);\r
+               \r
+               // FIXME: primary vertex\r
+               //TVector3 r0(AliGenerator::fVertex[0],AliGenerator::fVertex[1],AliGenerator::fVertex[2]);\r
+               \r
+               // workaround for primary vertex\r
+               TParticle* prim = stack->Particle(0);\r
+               TVector3 r0(prim->Vx(),prim->Vy(),prim->Vz()); // primary vertex\r
+               \r
+               for (Int_t i=0; i<stack->GetNtrack(); ++i)\r
+               {\r
+                       TParticle* iParticle = stack->Particle(i);\r
+                       \r
+                       if(iParticle->TestBit(kDoneBit)) continue;\r
+                       // select only nucleons within the freeze-out volume\r
+                       TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());\r
+                       if((r-r0).Mag()>fRsrc) continue;\r
+                       \r
+                       Int_t pdgCode = iParticle->GetPdgCode();\r
+                       if(pdgCode == fSign*2212)// (anti)proton\r
+                       {\r
+                               FixProductionVertex(iParticle,rnd2);\r
+                               protons->Add(iParticle);\r
+                       }\r
+                       else if(pdgCode == fSign*2112) // (anti)neutron\r
+                       {\r
+                               FixProductionVertex(iParticle,rnd2);\r
+                               neutrons->Add(iParticle);\r
+                       }\r
+               }\r
+               \r
+               // coalescence conditions\r
+               // A.J. Baltz et al., Phys. lett B 325(1994)7\r
+               \r
+               TIter p_next(protons);\r
+               while(TParticle* n0 = (TParticle*) p_next())\r
+               {\r
+                       TParticle* iProton = 0;\r
+                       TParticle* jNeutron = 0;\r
+                       \r
+                       // Select the first nucleon\r
+                       iProton = n0;\r
+                       \r
+                       TVector3 v0(n0->Vx(), n0->Vy(), n0->Vz());\r
+                       TVector3 p0(n0->Px(), n0->Py(), n0->Pz()), p1(0,0,0);\r
+                       \r
+                       // See if there is a neibourgh...\r
+                       TIter n_next(neutrons);\r
+                       while(TParticle* n1 = (TParticle*) n_next() )\r
+                       {\r
+                               TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());\r
+                               TVector3 p(n1->Px(), n1->Py(), n1->Pz());\r
+                               \r
+                               if((p-p0).Mag() > fPmax) continue;\r
+                               if((v-v0).Mag() > fRmax) continue;\r
+                               \r
+                               jNeutron = n1;\r
+                               p1 = p;\r
+                               break;\r
+                       }\r
+                       \r
+                       if(jNeutron == 0) continue; // with next proton\r
+                       \r
+                       if(rnd.Rndm() > fSpinProb) continue;\r
+                       \r
+                       // neutron captured!\r
+                       \r
+                       TVector3 pN = p0+p1;\r
+                       TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex\r
+                       \r
+                       // E^2 = p^2 + m^2\r
+                       Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);\r
+                       \r
+                       // Add a new (anti)deuteron to the current event stack\r
+                       stack->PushTrack(1, stack->Particles()->IndexOf(iProton), fSign*1000010020,\r
+                                        pN.X(), pN.Y(), pN.Z(), EN,\r
+                                        vN.X(), vN.Y(), vN.Z(), iProton->T(),\r
+                                        0., 0., 0., kPNCapture, ntr,1.,0);\r
+                       \r
+                       //Info("Generate","neutron capture (NEW DEUTERON)");\r
+                       \r
+                       // Set bit not to transport for the nucleons\r
+                       iProton->SetBit(kDoneBit);\r
+                       jNeutron->SetBit(kDoneBit);\r
+                       \r
+                       // Remove from the list\r
+                       neutrons->Remove(jNeutron);\r
+               }\r
+               \r
+               protons->Clear("nodelete");\r
+               neutrons->Clear("nodelete");\r
+               \r
+               delete protons;\r
+               delete neutrons;\r
+       }\r
+       \r
+       Info("Generate","Coalescence afterburner: DONE");\r
+}\r
+\r
+// create the freeze-out nucleon distribution around the collision vertex\r
+void AliGenDeuteron::FixProductionVertex(TParticle* i, TRandom3& rnd)\r
+{\r
+       Double_t x,y,z;\r
+       \r
+       if(fModel == kThermal) // homogeneous volume\r
+       {\r
+               // random (r,theta,phi)\r
+               Double_t r = fRsrc*TMath::Power(rnd.Rndm(),1./3.);\r
+               Double_t theta = TMath::ACos(2.*rnd.Rndm()-1.);\r
+               Double_t phi = 2.*TMath::Pi()*rnd.Rndm();\r
+       \r
+               // transform coordenates\r
+               x = r*TMath::Sin(theta)*TMath::Cos(phi);\r
+               y = r*TMath::Sin(theta)*TMath::Sin(phi);\r
+               z = r*TMath::Cos(theta);\r
+       }\r
+       else // projection into the surface of an sphere\r
+       {\r
+               x = fRsrc*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());\r
+               y = fRsrc*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());\r
+               z = fRsrc*TMath::Cos(i->Theta());\r
+       }\r
+       \r
+       // assume nucleons in the freeze-out have the same production vertex\r
+       i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
+}\r

diff --git a/EVGEN/AliGenDeuteron.h b/EVGEN/AliGenDeuteron.h
new file mode 100644 (file)
index 0000000..1c56f56
--- /dev/null
+++ b/EVGEN/AliGenDeuteron.h
@@ -0,0 +1,57 @@
+#ifndef _ALIGENDEUTERON_H_\r
+#define _ALIGENDEUTERON_H_\r
+\r
+// Afterburner to simulate coalescence of nucleons\r
+//\r
+\r
+#include "AliGenerator.h"\r
+\r
+class AliGenDeuteron: public AliGenerator\r
+{\r
+\r
+ public:\r
+\r
+       AliGenDeuteron(Int_t sign=1, Double_t pmax=0.3, Double_t rmax=2.1, Double_t rsrc=1.5, Int_t model=1);\r
+       virtual ~AliGenDeuteron();\r
+\r
+       virtual void Init();\r
+       virtual void Generate();\r
+       \r
+       Double_t GetCoalescenceDistance() const { return fRmax*1.e+13; }\r
+       Double_t GetCoalescenceMomentum() const { return fPmax; }\r
+       Double_t GetSourceRadius() const { return fRsrc*1.e+13; }\r
+       Double_t GetSpinProbability() const { return fSpinProb; }\r
+       Int_t GetFreezeOutModel() const { return fModel; }\r
+       Int_t GetSign() const { return fSign;}\r
+       \r
+       void SetCoalescenceDistance(Double_t r=2.1) { fRmax = r*1.e-13; }\r
+       void SetCoalescenceMomentum(Double_t p=0.3) { fPmax = p; }\r
+       void SetSourceRadius(Double_t r=1.5) { fRsrc = r*1.e-13; }\r
+       void SetSpinProbability(Double_t s=0.75) { fSpinProb = s; }\r
+       void SetFreezeOutModel(Int_t model = kThermal) { fModel = model; }\r
+       void SetSign(Int_t sign) {  fSign = sign > 0 ? 1 : -1;}\r
+       \r
+ public:\r
+       enum { kExpansion, kThermal };\r
+       \r
+ private:\r
+ \r
+       AliGenDeuteron(const AliGenDeuteron &other);\r
+       AliGenDeuteron& operator=(const AliGenDeuteron &other);\r
+       \r
+       void FixProductionVertex(class TParticle* i, class TRandom3& rnd);\r
+       \r
+ private:\r
+       \r
+       const Double_t fDeuteronMass;\r
+       Double_t fPmax; // Maximum p-n momentum difference (GeV/c)\r
+       Double_t fRmax; // Maximum p-n distance (cm)\r
+       Double_t fRsrc; // Emitting source radius (cm)\r
+       Double_t fSpinProb; // cluster formation probability due to spin\r
+       Int_t fModel;   // How to place the nucleons at freeze-out stage\r
+       Int_t fSign; // +1 for deuteron, -1 for antideuterons\r
+       \r
+       ClassDef(AliGenDeuteron,1)\r
+};\r
+\r
+#endif // _ALIGENDEUTERON_H_\r

diff --git a/EVGEN/EVGENLinkDef.h b/EVGEN/EVGENLinkDef.h
index 90b07c152fe64a7eba33911d909d078dbcfa6f85..9bf193d0484fa6ded113c202ab3122593c36f77c 100644 (file)
--- a/EVGEN/EVGENLinkDef.h
+++ b/EVGEN/EVGENLinkDef.h
@@ -62,4 +62,5 @@
 #pragma link C++ class  AliGenFunction+;
 #pragma link C++ class  AliGenTHnSparse+;
 #pragma link C++ class  AliOmegaDalitz+;
+#pragma link C++ class  AliGenDeuteron+;
 #endif

diff --git a/EVGEN/libEVGEN.pkg b/EVGEN/libEVGEN.pkg
index 240416abb159dc624032b306c51f2400533b128e..4ce1be35c0d41426223e9bf83a8784be0f3ecda0 100644 (file)
--- a/EVGEN/libEVGEN.pkg
+++ b/EVGEN/libEVGEN.pkg
@@ -21,7 +21,7 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliGenReaderEMD.cxx AliDecayerPolarized.cxx AliGenCorrHF.cxx AliGenCosmicsParam.cxx \
                AliGenKrypton.cxx AliGenThermalPhotons.cxx AliGenPromptPhotons.cxx\
                AliGenPileup.cxx AliGenFunction.cxx AliGenTHnSparse.cxx\
-               AliOmegaDalitz.cxx
+               AliOmegaDalitz.cxx AliGenDeuteron.cxx
 
 
 # Headerfiles for this particular package (Path respect to own directory)