1 /**************************************************************************
2 * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////
17 // Afterburner to generate (anti)deuterons simulating coalescence of
18 // (anti)nucleons as in A. J. Baltz et al., Phys. lett B 325(1994)7.
19 // If the nucleon generator does not provide a spatial description of
20 // the source the afterburner can provide one.
22 // There two options for the source: a thermal picture where all nucleons
23 // are placed randomly and homogeneously in an spherical volume, or
24 // an expansion one where they are projected onto its surface.
26 // An (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron
27 // with momentum difference less than ~ 300MeV and relative distance less than
28 // ~ 2.1fm. Only 3/4 of these clusters are accepted due to the deuteron spin.
30 // When there are more than one pair fulfilling the coalescence conditions,
31 // the selected pair can be the one with the first partner, with
32 // the lowest relative momentum, the lowest relative distance or both.
33 //////////////////////////////////////////////////////////////////////////
35 // Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
36 // Arturo Menchaca <menchaca@fisica.unam.mx>
39 #include "Riostream.h"
40 #include "TParticle.h"
44 #include "TMCProcess.h"
49 #include "AliCollisionGeometry.h"
50 #include "AliGenCocktailEventHeader.h"
51 #include "AliGenCocktailEntry.h"
52 #include "AliGenCocktailAfterBurner.h"
53 #include "AliGenDeuteron.h"
56 ClassImp(AliGenDeuteron)
58 AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
64 ,fClusterType(cluster)
76 fSign = sign > 0 ? 1:-1;
80 AliGenDeuteron::~AliGenDeuteron()
87 void AliGenDeuteron::Init()
90 // Standard AliGenerator initializer
94 void AliGenDeuteron::Generate()
97 // Look for coalescence of (anti)nucleons in the nucleon source
98 // provided by the generator cocktail
100 AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
101 if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
103 AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
105 // find nuclear radius, only for first generator and projectile=target
106 Bool_t collisionGeometry=0;
107 if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
110 Int_t ap, zp, at, zt;
111 gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
112 gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
113 if(ap != at) AliWarning("projectile and target have different size");
114 fR = 1.29*TMath::Power(at, 1./3.);
115 collisionGeometry = 1;
118 for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
120 gener->SetActiveEventNumber(ns);
122 if(fModel != kNone && collisionGeometry)
124 fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
125 fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
127 if(fB >= 2.*fR) continue; // no collision
132 gener->GetActiveEventHeader()->PrimaryVertex(primVtx);
133 TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);
135 // emerging nucleons from the collision
136 fCurStack = gener->GetStack(ns);
139 AliWarning("no event stack");
143 TList* protons = new TList();
144 protons->SetOwner(kFALSE);
145 TList* neutrons = new TList();
146 neutrons->SetOwner(kFALSE);
148 for (Int_t i=0; i < fCurStack->GetNtrack(); ++i)
150 TParticle* iParticle = fCurStack->Particle(i);
152 if(iParticle->TestBit(kDoneBit)) continue;
154 // select only nucleons within the freeze-out volume
155 TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());
156 if((r-r0).Mag() > fMaxRadius*1.e-13) continue;
158 Int_t pdgCode = iParticle->GetPdgCode();
159 if(pdgCode == fSign*2212)// (anti)proton
161 FixProductionVertex(iParticle);
162 protons->Add(iParticle);
164 else if(pdgCode == fSign*2112) // (anti)neutron
166 FixProductionVertex(iParticle);
167 neutrons->Add(iParticle);
172 if(fClusterType==kFirstPartner)
174 FirstPartner(protons, neutrons);
178 WeightMatrix(protons, neutrons);
181 protons->Clear("nodelete");
182 neutrons->Clear("nodelete");
191 Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
194 // Coalescence conditions as in
195 // A. J. Baltz et al., Phys. lett B 325(1994)7
197 // returns < 0 if coalescence is not possible
198 // otherwise returns a coalescence probability
200 TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
201 TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
203 TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
204 TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
206 Double_t deltaP = (p2-p1).Mag(); // relative momentum
207 if( deltaP >= fPmax) return -1.;
209 Double_t deltaR = (v2-v1).Mag(); // relative distance (cm)
210 if(deltaR >= fRmax*1.e-13) return -1.;
212 if(Rndm() > fSpinProb) return -1.; // spin
214 if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
215 if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
217 return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
220 void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
223 // Clusters are made with the first nucleon pair that fulfill
224 // the coalescence conditions, starting with the protons
226 TIter p_next(protons);
227 while(TParticle* n0 = (TParticle*) p_next())
229 TParticle* partner = 0;
230 TIter n_next(neutrons);
231 while(TParticle* n1 = (TParticle*) n_next() )
233 if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
238 if(partner == 0) continue; // with next proton
240 PushDeuteron(n0, partner);
242 // Remove from the list for the next iteration
243 neutrons->Remove(partner);
247 void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
250 // Build all possible nucleon pairs with their own probability
251 // and select only those with the highest coalescence probability
253 Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
255 TParticle** cProton = new TParticle*[nMaxPairs];
256 TParticle** cNeutron = new TParticle*[nMaxPairs];
257 Double_t* cWeight = new Double_t[nMaxPairs];
259 // build all pairs with probability > 0
261 TIter p_next(protons);
262 while(TParticle* n1 = (TParticle*) p_next())
264 TIter n_next(neutrons);
265 while(TParticle* n2 = (TParticle*) n_next() )
267 Double_t weight = this->GetCoalescenceProbability(n1,n2);
268 if(weight < 0) continue;
272 cWeight[cIdx] = weight;
278 Int_t nPairs = cIdx + 1;
280 // find the interacting pairs:
281 // remove repeated pairs and select only
282 // the pair with the highest coalescence probability
284 Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
286 TParticle** iProton = new TParticle*[nMaxIntPair];
287 TParticle** iNeutron = new TParticle*[nMaxIntPair];
288 Double_t* iWeight = new Double_t[nMaxIntPair];
295 for(Int_t i=0; i < nPairs; ++i)
297 if(cWeight[i] > wMax)
304 if(j == -1 ) break; // end
306 // Save the interacting pair
308 iProton[iIdx] = cProton[j];
309 iNeutron[iIdx] = cNeutron[j];
310 iWeight[iIdx] = cWeight[j];
312 // invalidate all combinations with these pairs for the next iteration
313 for(Int_t i=0; i < nPairs; ++i)
315 if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
316 if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
320 Int_t nIntPairs = iIdx + 1;
326 // Add the (anti)deuterons to the current event stack
327 for(Int_t i=0; i<nIntPairs; ++i)
329 TParticle* n1 = iProton[i];
330 TParticle* n2 = iNeutron[i];
339 void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
342 // Create an (anti)deuteron from parent1 and parent2,
343 // add to the current stack and set kDoneBit for the parents
345 const Double_t kDeuteronMass = 1.87561282;
346 const Int_t kDeuteronPdg = 1000010020;
349 TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
350 TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
353 // production vertex same as the parent1's
354 TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
357 Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
359 // Add a new (anti)deuteron to current event stack
360 fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
361 pN.X(), pN.Y(), pN.Z(), energy,
362 vN.X(), vN.Y(), vN.Z(), parent1->T(),
363 0., 0., 0., kPNCapture, fNtrk, 1., 0);
365 // Set kDoneBit for the parents
366 parent1->SetBit(kDoneBit);
367 parent2->SetBit(kDoneBit);
370 void AliGenDeuteron::FixProductionVertex(TParticle* i)
373 // Replace current generator nucleon spatial distribution
374 // with a custom distribution according to the selected model
376 if(fModel == kNone || fModel > kExpansion) return;
378 // semi-axis from collision geometry (fm)
379 Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.);
380 Double_t b = fTimeLength + fR - fB/2.;
381 Double_t c = fTimeLength;
387 if(fModel == kThermal)
389 // uniformly ditributed in the volume on an ellipsoid
390 // random (r,theta,phi) unit sphere
391 Double_t r = TMath::Power(Rndm(),1./3.);
392 Double_t theta = TMath::ACos(2.*Rndm()-1.);
393 Double_t phi = 2.*TMath::Pi()*Rndm();
395 // transform coordenates
396 xx = a*r*TMath::Sin(theta)*TMath::Cos(phi);
397 yy = b*r*TMath::Sin(theta)*TMath::Sin(phi);
398 zz = c*r*TMath::Cos(theta);
400 else if(fModel == kExpansion)
402 // project into the surface of an ellipsoid
403 xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());
404 yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());
405 zz = c*TMath::Cos(i->Theta());
408 // rotate by the reaction plane angle
409 Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR);
410 Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR);
413 // translate by the production vertex (cm)
414 i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());