]>
Commit | Line | Data |
---|---|---|
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" | |
41 | #include "AliRun.h" | |
42 | #include "AliStack.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 | ||
55 | ClassImp(AliGenDeuteron) | |
56 | ||
57 | AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster) | |
58 | :fSign(1) | |
59 | ,fPmax(pmax) | |
60 | ,fRmax(rmax) | |
61 | ,fSpinProb(0.75) | |
62 | ,fMaxRadius(1000.) | |
63 | ,fClusterType(cluster) | |
64 | ,fModel(0) | |
65 | ,fTimeLength(2.5) | |
66 | ,fB(0) | |
67 | ,fR(0) | |
68 | ,fPsiR(0) | |
69 | ,fCurStack(0) | |
70 | ,fNtrk(0) | |
71 | { | |
72 | // | |
73 | // constructor | |
74 | // | |
75 | fSign = sign > 0 ? 1:-1; | |
76 | ||
77 | } | |
78 | ||
79 | AliGenDeuteron::~AliGenDeuteron() | |
80 | { | |
81 | // | |
82 | // Default destructor | |
83 | // | |
84 | } | |
85 | ||
86 | void AliGenDeuteron::Init() | |
87 | { | |
88 | // | |
89 | // Standard AliGenerator initializer | |
90 | // | |
91 | } | |
92 | ||
93 | void AliGenDeuteron::Generate() | |
94 | { | |
95 | // | |
96 | // Look for coalescence of (anti)nucleons in the nucleon source | |
97 | // provided by the generator cocktail | |
98 | // | |
99 | AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType)); | |
100 | if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength)); | |
101 | ||
102 | AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator(); | |
103 | ||
104 | // find nuclear radius, only for first generator and projectile=target | |
105 | Bool_t collisionGeometry=0; | |
106 | if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry()) | |
107 | { | |
108 | TString name; | |
109 | Int_t ap, zp, at, zt; | |
110 | gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp); | |
111 | gener->FirstGenerator()->Generator()->GetTarget(name,at,zt); | |
112 | if(ap != at) AliWarning("projectile and target have different size"); | |
113 | fR = 1.29*TMath::Power(at, 1./3.); | |
114 | collisionGeometry = 1; | |
115 | } | |
116 | ||
117 | for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns) | |
118 | { | |
119 | gener->SetActiveEventNumber(ns); | |
120 | ||
121 | if(fModel != kNone && collisionGeometry) | |
122 | { | |
123 | fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle(); | |
124 | fB = (gener->GetCollisionGeometry(ns))->ImpactParameter(); | |
125 | ||
126 | if(fB >= 2.*fR) continue; // no collision | |
127 | } | |
128 | ||
129 | // primary vertex | |
130 | TArrayF primVtx; | |
131 | gener->GetActiveEventHeader()->PrimaryVertex(primVtx); | |
132 | TVector3 r0(primVtx[0],primVtx[1],primVtx[2]); | |
133 | ||
134 | // emerging nucleons from the collision | |
135 | fCurStack = gener->GetStack(ns); | |
136 | if(!fCurStack) | |
137 | { | |
138 | AliWarning("no event stack"); | |
139 | return; | |
140 | } | |
141 | ||
142 | TList* protons = new TList(); | |
143 | protons->SetOwner(kFALSE); | |
144 | TList* neutrons = new TList(); | |
145 | neutrons->SetOwner(kFALSE); | |
146 | ||
147 | for (Int_t i=0; i < fCurStack->GetNtrack(); ++i) | |
148 | { | |
149 | TParticle* iParticle = fCurStack->Particle(i); | |
150 | ||
151 | if(iParticle->TestBit(kDoneBit)) continue; | |
152 | ||
153 | // select only nucleons within the freeze-out volume | |
154 | TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz()); | |
155 | if((r-r0).Mag() > fMaxRadius*1.e-13) continue; | |
156 | ||
157 | Int_t pdgCode = iParticle->GetPdgCode(); | |
158 | if(pdgCode == fSign*2212)// (anti)proton | |
159 | { | |
160 | FixProductionVertex(iParticle); | |
161 | protons->Add(iParticle); | |
162 | } | |
163 | else if(pdgCode == fSign*2112) // (anti)neutron | |
164 | { | |
165 | FixProductionVertex(iParticle); | |
166 | neutrons->Add(iParticle); | |
167 | } | |
168 | } | |
169 | ||
170 | // look for clusters | |
171 | if(fClusterType==kFirstPartner) | |
172 | { | |
173 | FirstPartner(protons, neutrons); | |
174 | } | |
175 | else | |
176 | { | |
177 | WeightMatrix(protons, neutrons); | |
178 | } | |
179 | ||
180 | protons->Clear("nodelete"); | |
181 | neutrons->Clear("nodelete"); | |
182 | ||
183 | delete protons; | |
184 | delete neutrons; | |
185 | } | |
186 | ||
187 | AliInfo("DONE"); | |
188 | } | |
189 | ||
190 | Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const | |
191 | { | |
192 | // | |
193 | // Coalescence conditions as in | |
194 | // A. J. Baltz et al., Phys. lett B 325(1994)7 | |
195 | // | |
196 | // returns < 0 if coalescence is not possible | |
197 | // otherwise returns a coalescence probability | |
198 | // | |
199 | TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz()); | |
200 | TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz()); | |
201 | ||
202 | TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz()); | |
203 | TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz()); | |
204 | ||
205 | Double_t deltaP = (p2-p1).Mag(); // relative momentum | |
206 | if( deltaP >= fPmax) return -1.; | |
207 | ||
208 | Double_t deltaR = (v2-v1).Mag(); // relative distance (cm) | |
209 | if(deltaR >= fRmax*1.e-13) return -1.; | |
210 | ||
211 | if(Rndm() > fSpinProb) return -1.; // spin | |
212 | ||
213 | if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax; | |
214 | if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax; | |
215 | ||
216 | return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax); | |
217 | } | |
218 | ||
219 | void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons) | |
220 | { | |
221 | // | |
222 | // Clusters are made with the first nucleon pair that fulfill | |
223 | // the coalescence conditions, starting with the protons | |
224 | // | |
225 | TIter p_next(protons); | |
226 | while(TParticle* n0 = (TParticle*) p_next()) | |
227 | { | |
228 | TParticle* partner = 0; | |
229 | TIter n_next(neutrons); | |
230 | while(TParticle* n1 = (TParticle*) n_next() ) | |
231 | { | |
232 | if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron | |
233 | partner = n1; | |
234 | break; | |
235 | } | |
236 | ||
237 | if(partner == 0) continue; // with next proton | |
238 | ||
239 | PushDeuteron(n0, partner); | |
240 | ||
241 | // Remove from the list for the next iteration | |
242 | neutrons->Remove(partner); | |
243 | } | |
244 | } | |
245 | ||
246 | void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons) | |
247 | { | |
248 | // | |
249 | // Build all possible nucleon pairs with their own probability | |
250 | // and select only those with the highest coalescence probability | |
251 | // | |
252 | Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize(); | |
253 | ||
254 | TParticle** cProton = new TParticle*[nMaxPairs]; | |
255 | TParticle** cNeutron = new TParticle*[nMaxPairs]; | |
256 | Double_t* cWeight = new Double_t[nMaxPairs]; | |
257 | ||
258 | // build all pairs with probability > 0 | |
259 | Int_t cIdx = -1; | |
260 | TIter p_next(protons); | |
261 | while(TParticle* n1 = (TParticle*) p_next()) | |
262 | { | |
263 | TIter n_next(neutrons); | |
264 | while(TParticle* n2 = (TParticle*) n_next() ) | |
265 | { | |
266 | Double_t weight = this->GetCoalescenceProbability(n1,n2); | |
267 | if(weight < 0) continue; | |
268 | ++cIdx; | |
269 | cProton[cIdx] = n1; | |
270 | cNeutron[cIdx] = n2; | |
271 | cWeight[cIdx] = weight; | |
272 | } | |
273 | n_next.Reset(); | |
274 | } | |
275 | p_next.Reset(); | |
276 | ||
277 | Int_t nPairs = cIdx + 1; | |
278 | ||
279 | // find the interacting pairs: | |
280 | // remove repeated pairs and select only | |
281 | // the pair with the highest coalescence probability | |
282 | ||
283 | Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize()); | |
284 | ||
285 | TParticle** iProton = new TParticle*[nMaxIntPair]; | |
286 | TParticle** iNeutron = new TParticle*[nMaxIntPair]; | |
287 | Double_t* iWeight = new Double_t[nMaxIntPair]; | |
288 | ||
289 | Int_t iIdx = -1; | |
290 | while(true) | |
291 | { | |
292 | Int_t j = -1; | |
293 | Double_t wMax = 0; | |
294 | for(Int_t i=0; i < nPairs; ++i) | |
295 | { | |
296 | if(cWeight[i] > wMax) | |
297 | { | |
298 | wMax=cWeight[i]; | |
299 | j = i; | |
300 | } | |
301 | } | |
302 | ||
303 | if(j == -1 ) break; // end | |
304 | ||
305 | // Save the interacting pair | |
306 | ++iIdx; | |
307 | iProton[iIdx] = cProton[j]; | |
308 | iNeutron[iIdx] = cNeutron[j]; | |
309 | iWeight[iIdx] = cWeight[j]; | |
310 | ||
311 | // invalidate all combinations with these pairs for the next iteration | |
312 | for(Int_t i=0; i < nPairs; ++i) | |
313 | { | |
314 | if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.; | |
315 | if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.; | |
316 | } | |
317 | } | |
318 | ||
319 | Int_t nIntPairs = iIdx + 1; | |
320 | ||
321 | delete[] cProton; | |
322 | delete[] cNeutron; | |
323 | delete[] cWeight; | |
324 | ||
325 | // Add the (anti)deuterons to the current event stack | |
326 | for(Int_t i=0; i<nIntPairs; ++i) | |
327 | { | |
328 | TParticle* n1 = iProton[i]; | |
329 | TParticle* n2 = iNeutron[i]; | |
330 | PushDeuteron(n1,n2); | |
331 | } | |
332 | ||
333 | delete[] iProton; | |
334 | delete[] iNeutron; | |
335 | delete[] iWeight; | |
336 | } | |
337 | ||
338 | void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2) | |
339 | { | |
340 | // | |
341 | // Create an (anti)deuteron from parent1 and parent2, | |
342 | // add to the current stack and set kDoneBit for the parents | |
343 | // | |
344 | const Double_t kDeuteronMass = 1.87561282; | |
345 | const Int_t kDeuteronPdg = 1000010020; | |
346 | ||
347 | // momentum | |
348 | TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz()); | |
349 | TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz()); | |
350 | TVector3 pN = p1+p2; | |
351 | ||
352 | // production vertex same as the parent1's | |
353 | TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz()); | |
354 | ||
355 | // E^2 = p^2 + m^2 | |
356 | Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass); | |
357 | ||
358 | // Add a new (anti)deuteron to current event stack | |
0f52fba4 | 359 | fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg, |
be66cf70 | 360 | pN.X(), pN.Y(), pN.Z(), energy, |
361 | vN.X(), vN.Y(), vN.Z(), parent1->T(), | |
362 | 0., 0., 0., kPNCapture, fNtrk, 1., 0); | |
363 | ||
364 | // Set kDoneBit for the parents | |
365 | parent1->SetBit(kDoneBit); | |
366 | parent2->SetBit(kDoneBit); | |
367 | } | |
368 | ||
369 | void AliGenDeuteron::FixProductionVertex(TParticle* i) | |
370 | { | |
371 | // | |
372 | // Replace current generator nucleon spatial distribution | |
373 | // with a custom distribution according to the selected model | |
374 | // | |
375 | if(fModel == kNone || fModel > kExpansion) return; | |
376 | ||
377 | // semi-axis from collision geometry (fm) | |
378 | Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.); | |
379 | Double_t b = fTimeLength + fR - fB/2.; | |
380 | Double_t c = fTimeLength; | |
381 | ||
382 | Double_t xx = 0; | |
383 | Double_t yy = 0; | |
384 | Double_t zz = 0; | |
385 | ||
386 | if(fModel == kThermal) | |
387 | { | |
388 | // uniformly ditributed in the volume on an ellipsoid | |
389 | // random (r,theta,phi) unit sphere | |
390 | Double_t r = TMath::Power(Rndm(),1./3.); | |
391 | Double_t theta = TMath::ACos(2.*Rndm()-1.); | |
392 | Double_t phi = 2.*TMath::Pi()*Rndm(); | |
393 | ||
394 | // transform coordenates | |
395 | xx = a*r*TMath::Sin(theta)*TMath::Cos(phi); | |
396 | yy = b*r*TMath::Sin(theta)*TMath::Sin(phi); | |
397 | zz = c*r*TMath::Cos(theta); | |
398 | } | |
399 | else if(fModel == kExpansion) | |
400 | { | |
401 | // project into the surface of an ellipsoid | |
402 | xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi()); | |
403 | yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi()); | |
404 | zz = c*TMath::Cos(i->Theta()); | |
405 | } | |
406 | ||
407 | // rotate by the reaction plane angle | |
408 | Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR); | |
409 | Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR); | |
410 | Double_t z = zz; | |
411 | ||
412 | // translate by the production vertex (cm) | |
413 | i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T()); | |
414 | } | |
415 |