]>
Commit | Line | Data |
---|---|---|
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 | 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(1) | |
63 | ,fClusterType(cluster) | |
64 | ,fModel(0) | |
65 | ,fTimeLength(2.5) | |
66 | ,fB(0) | |
67 | ,fR(0) | |
68 | ,fPsiR(0) | |
69 | ,fCurStack(0) | |
70 | { | |
71 | // | |
72 | // constructor | |
73 | // | |
74 | fSign = sign > 0 ? 1:-1; | |
75 | ||
76 | } | |
77 | ||
78 | AliGenDeuteron::~AliGenDeuteron() | |
79 | { | |
80 | // | |
81 | // Default destructor | |
82 | // | |
83 | } | |
84 | ||
85 | void AliGenDeuteron::Init() | |
86 | { | |
87 | // | |
88 | // Standard AliGenerator initializer | |
89 | // | |
90 | } | |
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 | 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 | ||
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 | ||
140 | // particles produced by the generator | |
141 | for (Int_t i=0; i < fCurStack->GetNprimary(); ++i) | |
142 | { | |
143 | TParticle* iParticle = fCurStack->Particle(i); | |
144 | ||
145 | if(iParticle->GetStatusCode() != 1) continue; | |
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 | ||
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 | 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 | ||
198 | Double_t deltaP = this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame | |
199 | if( deltaP >= fPmax) return -1.; | |
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 | ||
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 | 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 | ||
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 | 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; | |
283 | while(kTRUE) | |
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 | ||
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 | 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 | ||
351 | 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 | fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg, | |
357 | pN.X(), pN.Y(), pN.Z(), energy, | |
358 | 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 | parent1->SetStatusCode(kCluster); | |
363 | parent2->SetStatusCode(kCluster); | |
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 | ||
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 | 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 | ||
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 | 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 | ||
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 | return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2); | |
444 | } | |
445 |