Line data Source code
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"
41 : #include "AliStack.h"
42 : #include "AliRun.h"
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"
54 : #include "AliLog.h"
55 :
56 6 : ClassImp(AliGenDeuteron)
57 :
58 0 : AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
59 0 : :fSign(1)
60 0 : ,fPmax(pmax)
61 0 : ,fRmax(rmax)
62 0 : ,fSpinProb(1)
63 0 : ,fClusterType(cluster)
64 0 : ,fModel(0)
65 0 : ,fTimeLength(2.5)
66 0 : ,fB(0)
67 0 : ,fR(0)
68 0 : ,fPsiR(0)
69 0 : ,fCurStack(0)
70 0 : {
71 : //
72 : // constructor
73 : //
74 0 : fSign = sign > 0 ? 1:-1;
75 :
76 0 : }
77 :
78 0 : AliGenDeuteron::~AliGenDeuteron()
79 0 : {
80 : //
81 : // Default destructor
82 : //
83 0 : }
84 :
85 : void AliGenDeuteron::Init()
86 : {
87 : //
88 : // Standard AliGenerator initializer
89 : //
90 0 : }
91 :
92 : void AliGenDeuteron::Generate()
93 : {
94 : //
95 : // Look for coalescence of (anti)nucleons in the nucleon source
96 : // provided by the generator cocktail
97 : //
98 0 : AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
99 0 : if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
100 :
101 0 : AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
102 :
103 : // find nuclear radius, only for first generator and projectile=target
104 : Bool_t collisionGeometry=0;
105 0 : if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
106 : {
107 0 : TString name;
108 0 : Int_t ap, zp, at, zt;
109 0 : gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
110 0 : gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
111 0 : if(ap != at) AliWarning("projectile and target have different size");
112 0 : fR = 1.29*TMath::Power(at, 1./3.);
113 : collisionGeometry = 1;
114 0 : }
115 :
116 0 : for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
117 : {
118 0 : gener->SetActiveEventNumber(ns);
119 :
120 0 : if(fModel != kNone && collisionGeometry)
121 : {
122 0 : fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
123 0 : fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
124 :
125 0 : if(fB >= 2.*fR) continue; // no collision
126 : }
127 :
128 0 : fCurStack = gener->GetStack(ns);
129 0 : if(!fCurStack)
130 : {
131 0 : AliWarning("no event stack");
132 0 : return;
133 : }
134 :
135 0 : TList* protons = new TList();
136 0 : protons->SetOwner(kFALSE);
137 0 : TList* neutrons = new TList();
138 0 : neutrons->SetOwner(kFALSE);
139 :
140 : // particles produced by the generator
141 0 : for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
142 : {
143 0 : TParticle* iParticle = fCurStack->Particle(i);
144 :
145 0 : if(iParticle->GetStatusCode() != 1) continue;
146 :
147 0 : Int_t pdgCode = iParticle->GetPdgCode();
148 0 : if(pdgCode == fSign*2212)// (anti)proton
149 : {
150 0 : FixProductionVertex(iParticle);
151 0 : protons->Add(iParticle);
152 0 : }
153 0 : else if(pdgCode == fSign*2112) // (anti)neutron
154 : {
155 0 : FixProductionVertex(iParticle);
156 0 : neutrons->Add(iParticle);
157 0 : }
158 0 : }
159 :
160 : // look for clusters
161 0 : if(fClusterType==kFirstPartner)
162 : {
163 0 : FirstPartner(protons, neutrons);
164 0 : }
165 : else
166 : {
167 0 : WeightMatrix(protons, neutrons);
168 : }
169 :
170 0 : protons->Clear("nodelete");
171 0 : neutrons->Clear("nodelete");
172 :
173 0 : delete protons;
174 0 : delete neutrons;
175 0 : }
176 :
177 0 : AliInfo("DONE");
178 0 : }
179 :
180 : Double_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 : //
189 : const Double_t kProtonMass = 0.938272013;
190 : const Double_t kNeutronMass = 0.939565378;
191 :
192 0 : TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
193 0 : TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
194 :
195 0 : TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
196 0 : TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
197 :
198 0 : Double_t deltaP = this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
199 0 : if( deltaP >= fPmax) return -1.;
200 :
201 0 : Double_t deltaR = (v2-v1).Mag(); // relative distance (cm)
202 0 : if(deltaR >= fRmax*1.e-13) return -1.;
203 :
204 0 : if(Rndm() > fSpinProb) return -1.; // spin
205 :
206 0 : if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
207 0 : if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
208 :
209 0 : return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
210 0 : }
211 :
212 : void 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 0 : TIter p_next(protons);
219 0 : while(TParticle* n0 = (TParticle*) p_next())
220 : {
221 : TParticle* partner = 0;
222 0 : TIter n_next(neutrons);
223 0 : while(TParticle* n1 = (TParticle*) n_next() )
224 : {
225 0 : if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
226 : partner = n1;
227 0 : break;
228 : }
229 :
230 0 : if(partner == 0) continue; // with next proton
231 :
232 0 : PushDeuteron(n0, partner);
233 :
234 : // Remove from the list for the next iteration
235 0 : neutrons->Remove(partner);
236 0 : }
237 0 : }
238 :
239 : void 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 0 : Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
246 :
247 0 : TParticle** cProton = new TParticle*[nMaxPairs];
248 0 : TParticle** cNeutron = new TParticle*[nMaxPairs];
249 0 : Double_t* cWeight = new Double_t[nMaxPairs];
250 :
251 : // build all pairs with probability > 0
252 : Int_t cIdx = -1;
253 0 : TIter p_next(protons);
254 0 : while(TParticle* n1 = (TParticle*) p_next())
255 : {
256 0 : TIter n_next(neutrons);
257 0 : while(TParticle* n2 = (TParticle*) n_next() )
258 : {
259 0 : Double_t weight = this->GetCoalescenceProbability(n1,n2);
260 0 : if(weight < 0) continue;
261 0 : ++cIdx;
262 0 : cProton[cIdx] = n1;
263 0 : cNeutron[cIdx] = n2;
264 0 : cWeight[cIdx] = weight;
265 0 : }
266 0 : n_next.Reset();
267 0 : }
268 0 : p_next.Reset();
269 :
270 0 : 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 0 : Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
277 :
278 0 : TParticle** iProton = new TParticle*[nMaxIntPair];
279 0 : TParticle** iNeutron = new TParticle*[nMaxIntPair];
280 0 : Double_t* iWeight = new Double_t[nMaxIntPair];
281 :
282 : Int_t iIdx = -1;
283 0 : while(kTRUE)
284 : {
285 : Int_t j = -1;
286 : Double_t wMax = 0;
287 0 : for(Int_t i=0; i < nPairs; ++i)
288 : {
289 0 : if(cWeight[i] > wMax)
290 : {
291 : wMax=cWeight[i];
292 : j = i;
293 0 : }
294 : }
295 :
296 0 : if(j == -1 ) break; // end
297 :
298 : // Save the interacting pair
299 0 : ++iIdx;
300 0 : iProton[iIdx] = cProton[j];
301 0 : iNeutron[iIdx] = cNeutron[j];
302 0 : iWeight[iIdx] = cWeight[j];
303 :
304 : // invalidate all combinations with these pairs for the next iteration
305 0 : for(Int_t i=0; i < nPairs; ++i)
306 : {
307 0 : if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
308 0 : if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
309 : }
310 0 : }
311 :
312 0 : Int_t nIntPairs = iIdx + 1;
313 :
314 0 : delete[] cProton;
315 0 : delete[] cNeutron;
316 0 : delete[] cWeight;
317 :
318 : // Add the (anti)deuterons to the current event stack
319 0 : for(Int_t i=0; i<nIntPairs; ++i)
320 : {
321 0 : TParticle* n1 = iProton[i];
322 0 : TParticle* n2 = iNeutron[i];
323 0 : PushDeuteron(n1,n2);
324 : }
325 :
326 0 : delete[] iProton;
327 0 : delete[] iNeutron;
328 0 : delete[] iWeight;
329 0 : }
330 :
331 : void 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 0 : TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
342 0 : TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
343 0 : TVector3 pN = p1+p2;
344 :
345 : // production vertex same as the parent1's
346 0 : TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
347 :
348 : // E^2 = p^2 + m^2
349 0 : Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
350 :
351 0 : Int_t ntrk = 0;
352 : Double_t weight = 1;
353 : Int_t is = 1; // final state particle
354 :
355 : // Add a new (anti)deuteron to current event stack
356 0 : fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
357 0 : pN.X(), pN.Y(), pN.Z(), energy,
358 0 : vN.X(), vN.Y(), vN.Z(), parent1->T(),
359 : 0., 0., 0., kPNCapture, ntrk, weight, is);
360 :
361 : // change the status code of the parents
362 0 : parent1->SetStatusCode(kCluster);
363 0 : parent2->SetStatusCode(kCluster);
364 :
365 : // Set kDoneBit for the parents
366 0 : parent1->SetBit(kDoneBit);
367 0 : parent2->SetBit(kDoneBit);
368 0 : }
369 :
370 : void 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 0 : if(fModel == kNone || fModel > kExpansion) return;
377 :
378 : // semi-axis from collision geometry (fm)
379 0 : Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.);
380 0 : 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 0 : if(fModel == kThermal)
388 : {
389 : // uniformly ditributed in the volume on an ellipsoid
390 : // random (r,theta,phi) unit sphere
391 0 : Double_t r = TMath::Power(Rndm(),1./3.);
392 0 : Double_t theta = TMath::ACos(2.*Rndm()-1.);
393 0 : Double_t phi = 2.*TMath::Pi()*Rndm();
394 :
395 : // transform coordenates
396 0 : xx = a*r*TMath::Sin(theta)*TMath::Cos(phi);
397 0 : yy = b*r*TMath::Sin(theta)*TMath::Sin(phi);
398 0 : zz = c*r*TMath::Cos(theta);
399 0 : }
400 0 : else if(fModel == kExpansion)
401 : {
402 : // project into the surface of an ellipsoid
403 0 : xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());
404 0 : yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());
405 0 : zz = c*TMath::Cos(i->Theta());
406 0 : }
407 :
408 : // rotate by the reaction plane angle
409 0 : Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR);
410 0 : 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 0 : i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
415 0 : }
416 :
417 : 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
418 : {
419 : //
420 : // square of the center of mass energy
421 : //
422 0 : Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
423 0 : Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
424 :
425 0 : return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
426 : }
427 :
428 : 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
429 : {
430 : //
431 : // momentum in the CM frame
432 : //
433 0 : Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
434 :
435 0 : return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
436 : }
437 :
438 : Double_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 0 : return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
444 : }
445 :
|
