]>
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" | |
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 | |
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) | |
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) | |
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 | ||
be66cf70 | 129 | fCurStack = gener->GetStack(ns); |
130 | if(!fCurStack) | |
131 | { | |
132 | AliWarning("no event stack"); | |
133 | return; | |
134 | } | |
135 | ||
136 | TList* protons = new TList(); | |
137 | protons->SetOwner(kFALSE); | |
138 | TList* neutrons = new TList(); | |
139 | neutrons->SetOwner(kFALSE); | |
140 | ||
95145d47 | 141 | // particles produced by the generator |
142 | for (Int_t i=0; i < fCurStack->GetNprimary(); ++i) | |
be66cf70 | 143 | { |
144 | TParticle* iParticle = fCurStack->Particle(i); | |
145 | ||
95145d47 | 146 | if(iParticle->GetStatusCode() != 1) continue; |
be66cf70 | 147 | |
148 | Int_t pdgCode = iParticle->GetPdgCode(); | |
149 | if(pdgCode == fSign*2212)// (anti)proton | |
150 | { | |
151 | FixProductionVertex(iParticle); | |
152 | protons->Add(iParticle); | |
153 | } | |
154 | else if(pdgCode == fSign*2112) // (anti)neutron | |
155 | { | |
156 | FixProductionVertex(iParticle); | |
157 | neutrons->Add(iParticle); | |
158 | } | |
159 | } | |
160 | ||
161 | // look for clusters | |
162 | if(fClusterType==kFirstPartner) | |
163 | { | |
164 | FirstPartner(protons, neutrons); | |
165 | } | |
166 | else | |
167 | { | |
168 | WeightMatrix(protons, neutrons); | |
169 | } | |
170 | ||
171 | protons->Clear("nodelete"); | |
172 | neutrons->Clear("nodelete"); | |
173 | ||
174 | delete protons; | |
175 | delete neutrons; | |
176 | } | |
177 | ||
178 | AliInfo("DONE"); | |
179 | } | |
180 | ||
181 | Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const | |
182 | { | |
183 | // | |
184 | // Coalescence conditions as in | |
185 | // A. J. Baltz et al., Phys. lett B 325(1994)7 | |
186 | // | |
187 | // returns < 0 if coalescence is not possible | |
188 | // otherwise returns a coalescence probability | |
189 | // | |
95145d47 | 190 | const Double_t kProtonMass = 0.938272013; |
191 | const Double_t kNeutronMass = 0.939565378; | |
192 | ||
be66cf70 | 193 | TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz()); |
194 | TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz()); | |
195 | ||
196 | TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz()); | |
197 | TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz()); | |
198 | ||
95145d47 | 199 | Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame |
200 | if( deltaP >= fPmax) return -1.; | |
be66cf70 | 201 | |
202 | Double_t deltaR = (v2-v1).Mag(); // relative distance (cm) | |
203 | if(deltaR >= fRmax*1.e-13) return -1.; | |
204 | ||
205 | if(Rndm() > fSpinProb) return -1.; // spin | |
206 | ||
207 | if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax; | |
208 | if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax; | |
209 | ||
210 | return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax); | |
211 | } | |
212 | ||
213 | void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons) | |
214 | { | |
215 | // | |
216 | // Clusters are made with the first nucleon pair that fulfill | |
217 | // the coalescence conditions, starting with the protons | |
218 | // | |
219 | TIter p_next(protons); | |
220 | while(TParticle* n0 = (TParticle*) p_next()) | |
221 | { | |
222 | TParticle* partner = 0; | |
223 | TIter n_next(neutrons); | |
224 | while(TParticle* n1 = (TParticle*) n_next() ) | |
225 | { | |
226 | if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron | |
227 | partner = n1; | |
228 | break; | |
229 | } | |
230 | ||
231 | if(partner == 0) continue; // with next proton | |
232 | ||
233 | PushDeuteron(n0, partner); | |
234 | ||
235 | // Remove from the list for the next iteration | |
236 | neutrons->Remove(partner); | |
237 | } | |
238 | } | |
239 | ||
240 | void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons) | |
241 | { | |
242 | // | |
243 | // Build all possible nucleon pairs with their own probability | |
244 | // and select only those with the highest coalescence probability | |
245 | // | |
246 | Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize(); | |
247 | ||
248 | TParticle** cProton = new TParticle*[nMaxPairs]; | |
249 | TParticle** cNeutron = new TParticle*[nMaxPairs]; | |
250 | Double_t* cWeight = new Double_t[nMaxPairs]; | |
251 | ||
252 | // build all pairs with probability > 0 | |
253 | Int_t cIdx = -1; | |
254 | TIter p_next(protons); | |
255 | while(TParticle* n1 = (TParticle*) p_next()) | |
256 | { | |
257 | TIter n_next(neutrons); | |
258 | while(TParticle* n2 = (TParticle*) n_next() ) | |
259 | { | |
260 | Double_t weight = this->GetCoalescenceProbability(n1,n2); | |
261 | if(weight < 0) continue; | |
262 | ++cIdx; | |
263 | cProton[cIdx] = n1; | |
264 | cNeutron[cIdx] = n2; | |
265 | cWeight[cIdx] = weight; | |
266 | } | |
267 | n_next.Reset(); | |
268 | } | |
269 | p_next.Reset(); | |
270 | ||
271 | Int_t nPairs = cIdx + 1; | |
272 | ||
273 | // find the interacting pairs: | |
274 | // remove repeated pairs and select only | |
275 | // the pair with the highest coalescence probability | |
276 | ||
277 | Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize()); | |
278 | ||
279 | TParticle** iProton = new TParticle*[nMaxIntPair]; | |
280 | TParticle** iNeutron = new TParticle*[nMaxIntPair]; | |
281 | Double_t* iWeight = new Double_t[nMaxIntPair]; | |
282 | ||
283 | Int_t iIdx = -1; | |
95145d47 | 284 | while(kTRUE) |
be66cf70 | 285 | { |
286 | Int_t j = -1; | |
287 | Double_t wMax = 0; | |
288 | for(Int_t i=0; i < nPairs; ++i) | |
289 | { | |
290 | if(cWeight[i] > wMax) | |
291 | { | |
292 | wMax=cWeight[i]; | |
293 | j = i; | |
294 | } | |
295 | } | |
296 | ||
297 | if(j == -1 ) break; // end | |
298 | ||
299 | // Save the interacting pair | |
300 | ++iIdx; | |
301 | iProton[iIdx] = cProton[j]; | |
302 | iNeutron[iIdx] = cNeutron[j]; | |
303 | iWeight[iIdx] = cWeight[j]; | |
304 | ||
305 | // invalidate all combinations with these pairs for the next iteration | |
306 | for(Int_t i=0; i < nPairs; ++i) | |
307 | { | |
308 | if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.; | |
309 | if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.; | |
310 | } | |
311 | } | |
312 | ||
313 | Int_t nIntPairs = iIdx + 1; | |
314 | ||
315 | delete[] cProton; | |
316 | delete[] cNeutron; | |
317 | delete[] cWeight; | |
318 | ||
319 | // Add the (anti)deuterons to the current event stack | |
320 | for(Int_t i=0; i<nIntPairs; ++i) | |
321 | { | |
322 | TParticle* n1 = iProton[i]; | |
323 | TParticle* n2 = iNeutron[i]; | |
324 | PushDeuteron(n1,n2); | |
325 | } | |
326 | ||
327 | delete[] iProton; | |
328 | delete[] iNeutron; | |
329 | delete[] iWeight; | |
330 | } | |
331 | ||
332 | void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2) | |
333 | { | |
334 | // | |
335 | // Create an (anti)deuteron from parent1 and parent2, | |
336 | // add to the current stack and set kDoneBit for the parents | |
337 | // | |
338 | const Double_t kDeuteronMass = 1.87561282; | |
339 | const Int_t kDeuteronPdg = 1000010020; | |
340 | ||
341 | // momentum | |
342 | TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz()); | |
343 | TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz()); | |
344 | TVector3 pN = p1+p2; | |
345 | ||
346 | // production vertex same as the parent1's | |
347 | TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz()); | |
348 | ||
349 | // E^2 = p^2 + m^2 | |
350 | Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass); | |
351 | ||
352 | // Add a new (anti)deuteron to current event stack | |
0f52fba4 | 353 | fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg, |
be66cf70 | 354 | pN.X(), pN.Y(), pN.Z(), energy, |
355 | vN.X(), vN.Y(), vN.Z(), parent1->T(), | |
356 | 0., 0., 0., kPNCapture, fNtrk, 1., 0); | |
357 | ||
358 | // Set kDoneBit for the parents | |
359 | parent1->SetBit(kDoneBit); | |
360 | parent2->SetBit(kDoneBit); | |
361 | } | |
362 | ||
363 | void AliGenDeuteron::FixProductionVertex(TParticle* i) | |
364 | { | |
365 | // | |
366 | // Replace current generator nucleon spatial distribution | |
367 | // with a custom distribution according to the selected model | |
368 | // | |
369 | if(fModel == kNone || fModel > kExpansion) return; | |
370 | ||
371 | // semi-axis from collision geometry (fm) | |
372 | Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.); | |
373 | Double_t b = fTimeLength + fR - fB/2.; | |
374 | Double_t c = fTimeLength; | |
375 | ||
376 | Double_t xx = 0; | |
377 | Double_t yy = 0; | |
378 | Double_t zz = 0; | |
379 | ||
380 | if(fModel == kThermal) | |
381 | { | |
382 | // uniformly ditributed in the volume on an ellipsoid | |
383 | // random (r,theta,phi) unit sphere | |
384 | Double_t r = TMath::Power(Rndm(),1./3.); | |
385 | Double_t theta = TMath::ACos(2.*Rndm()-1.); | |
386 | Double_t phi = 2.*TMath::Pi()*Rndm(); | |
387 | ||
388 | // transform coordenates | |
389 | xx = a*r*TMath::Sin(theta)*TMath::Cos(phi); | |
390 | yy = b*r*TMath::Sin(theta)*TMath::Sin(phi); | |
391 | zz = c*r*TMath::Cos(theta); | |
392 | } | |
393 | else if(fModel == kExpansion) | |
394 | { | |
395 | // project into the surface of an ellipsoid | |
396 | xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi()); | |
397 | yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi()); | |
398 | zz = c*TMath::Cos(i->Theta()); | |
399 | } | |
400 | ||
401 | // rotate by the reaction plane angle | |
402 | Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR); | |
403 | Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR); | |
404 | Double_t z = zz; | |
405 | ||
406 | // translate by the production vertex (cm) | |
407 | i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T()); | |
408 | } | |
409 | ||
95145d47 | 410 | 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 |
411 | { | |
412 | // | |
413 | // square of the center of mass energy | |
414 | // | |
415 | Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1); | |
416 | Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2); | |
417 | ||
418 | return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z)); | |
419 | } | |
420 | ||
421 | 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 | |
422 | { | |
423 | // | |
424 | // momentum in the CM frame | |
425 | // | |
426 | Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2); | |
427 | ||
428 | return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s)); | |
429 | } | |
430 | ||
431 | Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const | |
432 | { | |
433 | // | |
434 | // momentum in the CM frame | |
435 | // | |
436 | return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2); | |
437 | } | |
438 |