New class added.
[u/mrichter/AliRoot.git] / EVGEN / AliGenDeuteron.cxx
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 #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(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
360         fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
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