We've been changing the afterburner code to be able to use it with other
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 11 Mar 2011 12:38:15 +0000 (12:38 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 11 Mar 2011 12:38:15 +0000 (12:38 +0000)
generators and to add other optional ways of finding clusters. It works
reasonably well for pp although for PbPb still needs more testing to find a set
of good parameters.
Eulogio Serradilla <eulogio.serradilla@cern.ch>

EVGEN/AliGenDeuteron.cxx
EVGEN/AliGenDeuteron.h

index 1995115..7d0d570 100644 (file)
-/**************************************************************************\r
- * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-\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
-// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,\r
-//         Arturo Menchaca <menchaca@fisica.unam.mx>\r
-//\r
-\r
-#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 "AliMC.h"\r
-#include "AliRun.h"\r
-#include "TArrayF.h"\r
-#include "AliGenCocktailEventHeader.h"\r
-\r
-ClassImp(AliGenDeuteron)\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
-       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
-               // primary vertex of current event\r
-               TArrayF primVtx;\r
-               gener->GetActiveEventHeader()->PrimaryVertex(primVtx);\r
-               TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);\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);\r
-                               protons->Add(iParticle);\r
-                       }\r
-                       else if(pdgCode == fSign*2112) // (anti)neutron\r
-                       {\r
-                               FixProductionVertex(iParticle);\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(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)\r
-{\r
-    // Fix for the production vertex\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(Rndm(),1./3.);\r
-               Double_t theta = TMath::ACos(2.*Rndm()-1.);\r
-               Double_t phi = 2.*TMath::Pi()*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 at the freeze-out have the same production vertex\r
-       i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
-}\r
+/**************************************************************************
+ * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////
+// Afterburner to generate (anti)deuterons simulating coalescence of
+// (anti)nucleons as in A. J. Baltz et al., Phys. lett B 325(1994)7.
+// If the nucleon generator does not provide a spatial description of 
+// the source the afterburner can provide one.
+//
+// There two options for the source: a thermal picture where all nucleons
+// are placed randomly and homogeneously in an spherical volume, or
+// an expansion one where they are projected onto its surface.
+//
+// An (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron
+// with momentum difference less than ~ 300MeV and relative distance less than
+// ~ 2.1fm. Only 3/4 of these clusters are accepted due to the deuteron spin.
+//
+// When there are more than one pair fulfilling the coalescence conditions,
+// the selected pair can be the one with the first partner, with
+// the lowest relative momentum, the lowest relative distance or both.
+//////////////////////////////////////////////////////////////////////////
+
+// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
+//         Arturo Menchaca <menchaca@fisica.unam.mx>
+//
+
+#include "Riostream.h"
+#include "TParticle.h"
+#include "AliRun.h"
+#include "AliStack.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"
+
+ClassImp(AliGenDeuteron)
+
+AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
+ :fSign(1)
+ ,fPmax(pmax)
+ ,fRmax(rmax)
+ ,fSpinProb(0.75)
+ ,fMaxRadius(1000.)
+ ,fClusterType(cluster)
+ ,fModel(0)
+ ,fTimeLength(2.5)
+ ,fB(0)
+ ,fR(0)
+ ,fPsiR(0)
+ ,fCurStack(0)
+ ,fNtrk(0)
+{
+//
+// constructor
+//
+       fSign = sign > 0 ? 1:-1;
+       
+}
+
+AliGenDeuteron::~AliGenDeuteron()
+{
+//
+// Default destructor
+//
+}
+
+void AliGenDeuteron::Init()
+{
+//
+// Standard AliGenerator initializer
+//
+}
+
+void AliGenDeuteron::Generate()
+{
+//
+// Look for coalescence of (anti)nucleons in the nucleon source
+// provided by the generator cocktail
+//
+       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();
+       
+       // find nuclear radius, only for first generator and projectile=target
+       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; // no collision
+               }
+               
+               // primary vertex
+               TArrayF primVtx;
+               gener->GetActiveEventHeader()->PrimaryVertex(primVtx);
+               TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);
+               
+               // emerging nucleons from the collision
+               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->GetNtrack(); ++i)
+               {
+                       TParticle* iParticle = fCurStack->Particle(i);
+                       
+                       if(iParticle->TestBit(kDoneBit)) continue;
+                       
+                       // select only nucleons within the freeze-out volume
+                       TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());
+                       if((r-r0).Mag() > fMaxRadius*1.e-13) continue;
+                       
+                       Int_t pdgCode = iParticle->GetPdgCode();
+                       if(pdgCode == fSign*2212)// (anti)proton
+                       {
+                               FixProductionVertex(iParticle);
+                               protons->Add(iParticle);
+                       }
+                       else if(pdgCode == fSign*2112) // (anti)neutron
+                       {
+                               FixProductionVertex(iParticle);
+                               neutrons->Add(iParticle);
+                       }
+               }
+               
+               // look for clusters
+               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
+{
+//
+// Coalescence conditions as in
+// A. J. Baltz et al., Phys. lett B 325(1994)7
+//
+// returns < 0 if coalescence is not possible
+// otherwise returns a coalescence probability
+//
+       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 = (p2-p1).Mag();       // relative momentum
+       if( deltaP >= fPmax)       return -1.;
+       
+       Double_t deltaR = (v2-v1).Mag();       // relative distance (cm)
+       if(deltaR >= fRmax*1.e-13) return -1.;
+       
+       if(Rndm() > fSpinProb)    return -1.;  // spin
+       
+       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)
+{
+//
+// Clusters are made with the first nucleon pair that fulfill
+// the coalescence conditions, starting with the protons
+//
+       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; // with next neutron
+                       partner = n1;
+                       break;
+               }
+               
+               if(partner == 0) continue; // with next proton
+               
+               PushDeuteron(n0, partner);
+               
+               // Remove from the list for the next iteration
+               neutrons->Remove(partner);
+       }
+}
+
+void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
+{
+//
+// Build all possible nucleon pairs with their own probability
+// and select only those with the highest coalescence probability
+//
+       Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
+       
+       TParticle** cProton = new TParticle*[nMaxPairs];
+       TParticle** cNeutron = new TParticle*[nMaxPairs];
+       Double_t*  cWeight = new Double_t[nMaxPairs];
+       
+       // build all pairs with probability > 0
+       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;
+       
+       // find the interacting pairs:
+       // remove repeated pairs and select only
+       // the pair with the highest coalescence probability
+       
+       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(true)
+       {
+               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; // end
+               
+               // Save the interacting pair
+               ++iIdx;
+               iProton[iIdx]  = cProton[j];
+               iNeutron[iIdx] = cNeutron[j];
+               iWeight[iIdx]  = cWeight[j];
+               
+               // invalidate all combinations with these pairs for the next iteration
+               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;
+       
+       // Add the (anti)deuterons to the current event stack
+       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)
+{
+//
+// Create an (anti)deuteron from parent1 and parent2,
+// add to the current stack and set kDoneBit for the parents
+//
+       const Double_t kDeuteronMass = 1.87561282;
+       const Int_t kDeuteronPdg = 1000010020;
+       
+       // momentum
+       TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
+       TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
+       TVector3 pN = p1+p2;
+       
+       // production vertex same as the parent1's
+       TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
+       
+       // E^2 = p^2 + m^2
+       Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
+       
+       // Add a new (anti)deuteron to current event stack
+       fCurStack->PushTrack(1, fCurStack->Particles()->IndexOf(parent1), fSign*kDeuteronPdg,
+                        pN.X(), pN.Y(), pN.Z(), energy,
+                        vN.X(), vN.Y(), vN.Z(), parent1->T(),
+                        0., 0., 0., kPNCapture, fNtrk, 1., 0);
+       
+       // Set kDoneBit for the parents
+       parent1->SetBit(kDoneBit);
+       parent2->SetBit(kDoneBit);
+}
+
+void AliGenDeuteron::FixProductionVertex(TParticle* i)
+{
+//
+// Replace current generator nucleon spatial distribution
+// with a custom distribution according to the selected model
+//
+       if(fModel == kNone || fModel > kExpansion) return;
+       
+       // semi-axis from collision geometry (fm)
+       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)
+       {
+               // uniformly ditributed in the volume on an ellipsoid
+               // random (r,theta,phi) unit sphere
+               Double_t r = TMath::Power(Rndm(),1./3.);
+               Double_t theta = TMath::ACos(2.*Rndm()-1.);
+               Double_t phi = 2.*TMath::Pi()*Rndm();
+               
+               // transform coordenates
+               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)
+       {
+               // project into the surface of an ellipsoid
+               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());
+       }
+       
+       // rotate by the reaction plane angle
+       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;
+       
+       // translate by the production vertex (cm)
+       i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
+}
+
index 4266ac2..9d974bf 100644 (file)
@@ -1,61 +1,78 @@
-#ifndef _ALIGENDEUTERON_H_\r
-#define _ALIGENDEUTERON_H_\r
-\r
-/* Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice                               */\r
-\r
-// Afterburner to simulate coalescence of nucleons\r
-// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>\r
-//         Arturo Menchaca <menchaca@fisica.unam.mx>\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);\r
-       \r
- private:\r
-       \r
-       const Double_t fDeuteronMass; // Deuteron mass\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
+#ifndef ALIGENDEUTERON_H
+#define ALIGENDEUTERON_H
+
+
+/* Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// Afterburner to simulate coalescence of (anti)nucleons
+// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>
+//         Arturo Menchaca <menchaca@fisica.unam.mx>
+
+#include "AliGenerator.h"
+
+
+class AliGenDeuteron: public AliGenerator
+{
+
+ public:
+
+       AliGenDeuteron(Int_t sign=1, Double_t pmax=0.3, Double_t rmax=2.1, Int_t cluster=0 );
+       virtual ~AliGenDeuteron();
+
+       virtual void Init();
+       virtual void Generate();
+       
+       Int_t GetSign() const { return fSign;}
+       Double_t GetCoalescenceMomentum() const { return fPmax; }
+       Double_t GetCoalescenceDistance() const { return fRmax; }
+       Double_t GetSpinProbability() const { return fSpinProb; }
+       Double_t GetFreezeOutRadius() const { return fMaxRadius; }
+       Double_t GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const;
+       Int_t GetClusterType() const { return fClusterType; }
+       Int_t GetFreezeOutModel() const { return fModel; }
+       Double_t GetFreezeOutTime() const { return fTimeLength; }
+       
+       void SetSign(Int_t sign) {  fSign = sign > 0 ? 1 : -1;}
+       void SetCoalescenceMomentum(Double_t p) { fPmax = p; }
+       void SetCoalescenceDistance(Double_t r) { fRmax = r; }
+       void SetSpinProbability(Double_t s) { fSpinProb = s; }
+       void SetFreezeOutRadius(Double_t r) { fMaxRadius = r; }
+       void SetClusterType(Int_t cluster) { fClusterType = cluster; }
+       void SetFreezeOutModel(Int_t model, Double_t timeLength=2.5) { fModel = model; fTimeLength=timeLength;}
+       
+ public:
+
+       enum { kFirstPartner, kLowestMomentum, kLowestDistance, kBoth };
+       enum { kNone, kThermal, kExpansion };
+       
+ private:
+       AliGenDeuteron(const AliGenDeuteron &other);
+       AliGenDeuteron& operator=(const AliGenDeuteron &other);
+       
+       void FixProductionVertex(class TParticle* i);
+       void FirstPartner(const TList* protons, TList* neutrons);
+       void WeightMatrix(const TList* protons, const TList* neutrons);
+       void PushDeuteron(TParticle* parent1, TParticle* parent2);
+       
+ private:
+       
+       Int_t fSign;          // +1 for deuterons, -1 for antideuterons
+       Double_t fPmax;       // Maximum p-n momentum difference (GeV/c)
+       Double_t fRmax;       // Maximum p-n distance (fm)
+       Double_t fSpinProb;   // cluster formation probability due to spin
+       Double_t fMaxRadius;  // Maximum freeze-out radius (fm)
+       Int_t fClusterType;   // Probability criteria to find clusters
+       Int_t fModel;         // Model to override generator source
+       Double_t fTimeLength; // Thermal and chemical freeze-out time (fm/c)
+       Double_t fB;          // Impact parameter (fm)
+       Double_t fR;          // Projectile/Target nuclear radius (fm)
+       Double_t fPsiR;       // Reaction plane angle
+       AliStack* fCurStack;  //! current event stack
+       Int_t fNtrk;          //! number of the stored track
+       
+       ClassDef(AliGenDeuteron,1)
+};
+
+#endif // ALIGENDEUTERON_H