doxy: do not show whitespace diffs on bulk edit
[u/mrichter/AliRoot.git] / EVGEN / AliGenDeuteron.cxx
CommitLineData
be66cf70 1/**************************************************************************
2 * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
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.
21//
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.
25//
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.
29//
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//////////////////////////////////////////////////////////////////////////
34
35// Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
36// Arturo Menchaca <menchaca@fisica.unam.mx>
37//
38
39#include "Riostream.h"
40#include "TParticle.h"
be66cf70 41#include "AliStack.h"
95145d47 42#include "AliRun.h"
be66cf70 43#include "TMath.h"
44#include "TMCProcess.h"
45#include "TList.h"
46#include "TVector3.h"
47#include "AliMC.h"
48#include "TArrayF.h"
49#include "AliCollisionGeometry.h"
50#include "AliGenCocktailEventHeader.h"
51#include "AliGenCocktailEntry.h"
52#include "AliGenCocktailAfterBurner.h"
53#include "AliGenDeuteron.h"
c93255fe 54#include "AliLog.h"
be66cf70 55
56ClassImp(AliGenDeuteron)
57
58AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
59 :fSign(1)
60 ,fPmax(pmax)
61 ,fRmax(rmax)
95145d47 62 ,fSpinProb(1)
be66cf70 63 ,fClusterType(cluster)
64 ,fModel(0)
65 ,fTimeLength(2.5)
66 ,fB(0)
67 ,fR(0)
68 ,fPsiR(0)
69 ,fCurStack(0)
be66cf70 70{
71//
72// constructor
73//
74 fSign = sign > 0 ? 1:-1;
75
76}
77
78AliGenDeuteron::~AliGenDeuteron()
79{
80//
81// Default destructor
82//
83}
84
85void AliGenDeuteron::Init()
86{
87//
88// Standard AliGenerator initializer
89//
90}
91
92void AliGenDeuteron::Generate()
93{
94//
95// Look for coalescence of (anti)nucleons in the nucleon source
96// provided by the generator cocktail
97//
98 AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
99 if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
100
101 AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
102
103 // find nuclear radius, only for first generator and projectile=target
104 Bool_t collisionGeometry=0;
105 if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
106 {
107 TString name;
108 Int_t ap, zp, at, zt;
109 gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
110 gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
111 if(ap != at) AliWarning("projectile and target have different size");
112 fR = 1.29*TMath::Power(at, 1./3.);
113 collisionGeometry = 1;
114 }
115
116 for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
117 {
118 gener->SetActiveEventNumber(ns);
119
120 if(fModel != kNone && collisionGeometry)
121 {
122 fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
123 fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
124
125 if(fB >= 2.*fR) continue; // no collision
126 }
127
be66cf70 128 fCurStack = gener->GetStack(ns);
129 if(!fCurStack)
130 {
131 AliWarning("no event stack");
132 return;
133 }
134
135 TList* protons = new TList();
136 protons->SetOwner(kFALSE);
137 TList* neutrons = new TList();
138 neutrons->SetOwner(kFALSE);
139
95145d47 140 // particles produced by the generator
141 for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
be66cf70 142 {
143 TParticle* iParticle = fCurStack->Particle(i);
144
95145d47 145 if(iParticle->GetStatusCode() != 1) continue;
be66cf70 146
147 Int_t pdgCode = iParticle->GetPdgCode();
148 if(pdgCode == fSign*2212)// (anti)proton
149 {
150 FixProductionVertex(iParticle);
151 protons->Add(iParticle);
152 }
153 else if(pdgCode == fSign*2112) // (anti)neutron
154 {
155 FixProductionVertex(iParticle);
156 neutrons->Add(iParticle);
157 }
158 }
159
160 // look for clusters
161 if(fClusterType==kFirstPartner)
162 {
163 FirstPartner(protons, neutrons);
164 }
165 else
166 {
167 WeightMatrix(protons, neutrons);
168 }
169
170 protons->Clear("nodelete");
171 neutrons->Clear("nodelete");
172
173 delete protons;
174 delete neutrons;
175 }
176
177 AliInfo("DONE");
178}
179
180Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
181{
182//
183// Coalescence conditions as in
184// A. J. Baltz et al., Phys. lett B 325(1994)7
185//
186// returns < 0 if coalescence is not possible
187// otherwise returns a coalescence probability
188//
95145d47 189 const Double_t kProtonMass = 0.938272013;
190 const Double_t kNeutronMass = 0.939565378;
191
be66cf70 192 TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
193 TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
194
195 TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
196 TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
197
3fdac9ed 198 Double_t deltaP = this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
95145d47 199 if( deltaP >= fPmax) return -1.;
be66cf70 200
201 Double_t deltaR = (v2-v1).Mag(); // relative distance (cm)
202 if(deltaR >= fRmax*1.e-13) return -1.;
203
204 if(Rndm() > fSpinProb) return -1.; // spin
205
206 if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
207 if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
208
209 return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
210}
211
212void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
213{
214//
215// Clusters are made with the first nucleon pair that fulfill
216// the coalescence conditions, starting with the protons
217//
218 TIter p_next(protons);
219 while(TParticle* n0 = (TParticle*) p_next())
220 {
221 TParticle* partner = 0;
222 TIter n_next(neutrons);
223 while(TParticle* n1 = (TParticle*) n_next() )
224 {
225 if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
226 partner = n1;
227 break;
228 }
229
230 if(partner == 0) continue; // with next proton
231
232 PushDeuteron(n0, partner);
233
234 // Remove from the list for the next iteration
235 neutrons->Remove(partner);
236 }
237}
238
239void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
240{
241//
242// Build all possible nucleon pairs with their own probability
243// and select only those with the highest coalescence probability
244//
245 Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
246
247 TParticle** cProton = new TParticle*[nMaxPairs];
248 TParticle** cNeutron = new TParticle*[nMaxPairs];
249 Double_t* cWeight = new Double_t[nMaxPairs];
250
251 // build all pairs with probability > 0
252 Int_t cIdx = -1;
253 TIter p_next(protons);
254 while(TParticle* n1 = (TParticle*) p_next())
255 {
256 TIter n_next(neutrons);
257 while(TParticle* n2 = (TParticle*) n_next() )
258 {
259 Double_t weight = this->GetCoalescenceProbability(n1,n2);
260 if(weight < 0) continue;
261 ++cIdx;
262 cProton[cIdx] = n1;
263 cNeutron[cIdx] = n2;
264 cWeight[cIdx] = weight;
265 }
266 n_next.Reset();
267 }
268 p_next.Reset();
269
270 Int_t nPairs = cIdx + 1;
271
272 // find the interacting pairs:
273 // remove repeated pairs and select only
274 // the pair with the highest coalescence probability
275
276 Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
277
278 TParticle** iProton = new TParticle*[nMaxIntPair];
279 TParticle** iNeutron = new TParticle*[nMaxIntPair];
280 Double_t* iWeight = new Double_t[nMaxIntPair];
281
282 Int_t iIdx = -1;
95145d47 283 while(kTRUE)
be66cf70 284 {
285 Int_t j = -1;
286 Double_t wMax = 0;
287 for(Int_t i=0; i < nPairs; ++i)
288 {
289 if(cWeight[i] > wMax)
290 {
291 wMax=cWeight[i];
292 j = i;
293 }
294 }
295
296 if(j == -1 ) break; // end
297
298 // Save the interacting pair
299 ++iIdx;
300 iProton[iIdx] = cProton[j];
301 iNeutron[iIdx] = cNeutron[j];
302 iWeight[iIdx] = cWeight[j];
303
304 // invalidate all combinations with these pairs for the next iteration
305 for(Int_t i=0; i < nPairs; ++i)
306 {
307 if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
308 if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
309 }
310 }
311
312 Int_t nIntPairs = iIdx + 1;
313
314 delete[] cProton;
315 delete[] cNeutron;
316 delete[] cWeight;
317
318 // Add the (anti)deuterons to the current event stack
319 for(Int_t i=0; i<nIntPairs; ++i)
320 {
321 TParticle* n1 = iProton[i];
322 TParticle* n2 = iNeutron[i];
323 PushDeuteron(n1,n2);
324 }
325
326 delete[] iProton;
327 delete[] iNeutron;
328 delete[] iWeight;
329}
330
331void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
332{
333//
334// Create an (anti)deuteron from parent1 and parent2,
335// add to the current stack and set kDoneBit for the parents
336//
337 const Double_t kDeuteronMass = 1.87561282;
338 const Int_t kDeuteronPdg = 1000010020;
339
340 // momentum
341 TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
342 TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
343 TVector3 pN = p1+p2;
344
345 // production vertex same as the parent1's
346 TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
347
348 // E^2 = p^2 + m^2
349 Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
350
3d87465c 351 Int_t ntrk = 0;
352 Double_t weight = 1;
353 Int_t is = 1; // final state particle
354
be66cf70 355 // Add a new (anti)deuteron to current event stack
0f52fba4 356 fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
be66cf70 357 pN.X(), pN.Y(), pN.Z(), energy,
358 vN.X(), vN.Y(), vN.Z(), parent1->T(),
3d87465c 359 0., 0., 0., kPNCapture, ntrk, weight, is);
360
361 // change the status code of the parents
362 parent1->SetStatusCode(kCluster);
363 parent2->SetStatusCode(kCluster);
be66cf70 364
365 // Set kDoneBit for the parents
366 parent1->SetBit(kDoneBit);
367 parent2->SetBit(kDoneBit);
368}
369
370void AliGenDeuteron::FixProductionVertex(TParticle* i)
371{
372//
373// Replace current generator nucleon spatial distribution
374// with a custom distribution according to the selected model
375//
376 if(fModel == kNone || fModel > kExpansion) return;
377
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;
382
383 Double_t xx = 0;
384 Double_t yy = 0;
385 Double_t zz = 0;
386
387 if(fModel == kThermal)
388 {
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();
394
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);
399 }
400 else if(fModel == kExpansion)
401 {
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());
406 }
407
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);
411 Double_t z = zz;
412
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());
415}
416
95145d47 417Double_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
418{
419//
420// square of the center of mass energy
421//
422 Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
423 Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
424
425 return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
426}
427
428Double_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
429{
430//
431// momentum in the CM frame
432//
433 Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
434
435 return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
436}
437
438Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
439{
440//
441// momentum in the CM frame
442//
443 return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
444}
445