bug fixes
[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 "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  ,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                 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                 
141                 // particles produced by the generator
142                 for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
143                 {
144                         TParticle* iParticle = fCurStack->Particle(i);
145                         
146                         if(iParticle->GetStatusCode() != 1) continue;
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 //
190         const Double_t kProtonMass  = 0.938272013;
191         const Double_t kNeutronMass = 0.939565378;
192         
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         
199         Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
200         if( deltaP >= fPmax) return -1.;
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;
284         while(kTRUE)
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
353         fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
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
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