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