Updates for mixing.
[u/mrichter/AliRoot.git] / EVGEN / AliGenDeuteron.cxx
1 /**************************************************************************\r
2  * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 //\r
17 // A very simple model based on the article "Physics Letters B325 (1994) 7-12".\r
18 // The only addition is to model the freeze-out at which the\r
19 // light nuclei can form for the energies of the LHC, since PYTHIA does not give\r
20 // any detailed spatial description of the generated particles at the level of a\r
21 // few fermis.\r
22 //\r
23 // There are two options to place the nucleons: a thermal picture where all\r
24 // nucleons are placed randomly and homogeneously in an spherical volume of a\r
25 // few fermis, and expansion one where they are projected on its surface.\r
26 //\r
27 // A (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron\r
28 // with momentum difference less than ~ 300MeV and relative distance less than a\r
29 // ~ 2.1fm. Only 3/4 of this clusters are accepted due to the deuteron spin.\r
30 //\r
31 // Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,\r
32 //         Arturo Menchaca <menchaca@fisica.unam.mx>\r
33 //\r
34 \r
35 #include "Riostream.h"\r
36 #include "TParticle.h"\r
37 #include "AliStack.h"\r
38 #include "AliGenDeuteron.h"\r
39 #include "AliGenCocktailAfterBurner.h"\r
40 #include "TMath.h"\r
41 #include "TMCProcess.h"\r
42 #include "TList.h"\r
43 #include "TVector3.h"\r
44 #include "AliMC.h"\r
45 #include "TArrayF.h"\r
46 #include "AliGenCocktailEventHeader.h"\r
47 \r
48 ClassImp(AliGenDeuteron)\r
49 \r
50 AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
51 fDeuteronMass(1.87561282), //pdg\r
52 fPmax(pmax), // GeV/c\r
53 fRmax(rmax*1.e-13), //cm\r
54 fRsrc(rsrc*1.e-13), // cm\r
55 fSpinProb(0.75),\r
56 fModel(model),\r
57 fSign(1)\r
58 {\r
59         fSign = sign > 0 ? 1:-1;\r
60 }\r
61 \r
62 AliGenDeuteron::~AliGenDeuteron()\r
63 {\r
64 }\r
65 \r
66 void AliGenDeuteron::Init()\r
67 {\r
68         // Standard AliGenerator Initializer\r
69 }\r
70 \r
71 void AliGenDeuteron::Generate()\r
72 {\r
73 // Modify the stack of each event, adding the new particles at the end and\r
74 // not transporting the parents.\r
75 \r
76         Info("Generate","freeze-out model : %d (0 expansion, 1 thermal)",fModel);\r
77         Info("Generate","relative momentum: %g GeV/c",fPmax);\r
78         Info("Generate","relative distance: %g cm",fRmax);\r
79         Info("Generate","source radius    : %g cm",fRsrc);\r
80         Info("Generate","spin probability : %g ",fSpinProb);\r
81         Info("Generate","sign             : %d ",fSign);\r
82         \r
83         Int_t ntr;\r
84         \r
85         // Get the cocktail generator\r
86         AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();\r
87         \r
88         // Loop over events\r
89         for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)\r
90         {\r
91                 gener->SetActiveEventNumber(ns);\r
92                 \r
93                 AliStack* stack = gener->GetStack(ns); // Stack of event ns\r
94                 if(!stack)\r
95                 {\r
96                         Info("Generate", "no stack, exiting");\r
97                         return;\r
98                 }\r
99                 \r
100                 // find emerging nucleons from the collision of current event\r
101                 \r
102                 TList* protons = new TList();\r
103                 protons->SetOwner(kFALSE);\r
104                 TList* neutrons = new TList();\r
105                 neutrons->SetOwner(kFALSE);\r
106                 \r
107                 // primary vertex of current event\r
108                 TArrayF primVtx;\r
109                 gener->GetActiveEventHeader()->PrimaryVertex(primVtx);\r
110                 TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);\r
111                 \r
112                 for (Int_t i=0; i<stack->GetNtrack(); ++i)\r
113                 {\r
114                         TParticle* iParticle = stack->Particle(i);\r
115                         \r
116                         if(iParticle->TestBit(kDoneBit)) continue;\r
117                         // select only nucleons within the freeze-out volume\r
118                         TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());\r
119                         if((r-r0).Mag()>fRsrc) continue;\r
120                         \r
121                         Int_t pdgCode = iParticle->GetPdgCode();\r
122                         if(pdgCode == fSign*2212)// (anti)proton\r
123                         {\r
124                                 FixProductionVertex(iParticle);\r
125                                 protons->Add(iParticle);\r
126                         }\r
127                         else if(pdgCode == fSign*2112) // (anti)neutron\r
128                         {\r
129                                 FixProductionVertex(iParticle);\r
130                                 neutrons->Add(iParticle);\r
131                         }\r
132                 }\r
133                 \r
134                 // coalescence conditions\r
135                 // A.J. Baltz et al., Phys. lett B 325(1994)7\r
136                 \r
137                 TIter p_next(protons);\r
138                 while(TParticle* n0 = (TParticle*) p_next())\r
139                 {\r
140                         TParticle* iProton = 0;\r
141                         TParticle* jNeutron = 0;\r
142                         \r
143                         // Select the first nucleon\r
144                         iProton = n0;\r
145                         \r
146                         TVector3 v0(n0->Vx(), n0->Vy(), n0->Vz());\r
147                         TVector3 p0(n0->Px(), n0->Py(), n0->Pz()), p1(0,0,0);\r
148                         \r
149                         // See if there is a neibourgh...\r
150                         TIter n_next(neutrons);\r
151                         while(TParticle* n1 = (TParticle*) n_next() )\r
152                         {\r
153                                 TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());\r
154                                 TVector3 p(n1->Px(), n1->Py(), n1->Pz());\r
155                                 \r
156                                 if((p-p0).Mag() > fPmax) continue;\r
157                                 if((v-v0).Mag() > fRmax) continue;\r
158                                 \r
159                                 jNeutron = n1;\r
160                                 p1 = p;\r
161                                 break;\r
162                         }\r
163                         \r
164                         if(jNeutron == 0) continue; // with next proton\r
165                         \r
166                         if(Rndm() > fSpinProb) continue;\r
167                         \r
168                         // neutron captured!\r
169                         \r
170                         TVector3 pN = p0+p1;\r
171                         TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex\r
172                         \r
173                         // E^2 = p^2 + m^2\r
174                         Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);\r
175                         \r
176                         // Add a new (anti)deuteron to the current event stack\r
177                         stack->PushTrack(1, stack->Particles()->IndexOf(iProton), fSign*1000010020,\r
178                                          pN.X(), pN.Y(), pN.Z(), EN,\r
179                                          vN.X(), vN.Y(), vN.Z(), iProton->T(),\r
180                                          0., 0., 0., kPNCapture, ntr,1.,0);\r
181                         \r
182                         //Info("Generate","neutron capture (NEW DEUTERON)");\r
183                         \r
184                         // Set bit not to transport for the nucleons\r
185                         iProton->SetBit(kDoneBit);\r
186                         jNeutron->SetBit(kDoneBit);\r
187                         \r
188                         // Remove from the list\r
189                         neutrons->Remove(jNeutron);\r
190                 }\r
191                 \r
192                 protons->Clear("nodelete");\r
193                 neutrons->Clear("nodelete");\r
194                 \r
195                 delete protons;\r
196                 delete neutrons;\r
197         }\r
198         \r
199         Info("Generate","Coalescence afterburner: DONE");\r
200 }\r
201 \r
202 // create the freeze-out nucleon distribution around the collision vertex\r
203 void AliGenDeuteron::FixProductionVertex(TParticle* i)\r
204 {\r
205         Double_t x,y,z;\r
206         \r
207         if(fModel == kThermal) // homogeneous volume\r
208         {\r
209                 // random (r,theta,phi)\r
210                 Double_t r = fRsrc*TMath::Power(Rndm(),1./3.);\r
211                 Double_t theta = TMath::ACos(2.*Rndm()-1.);\r
212                 Double_t phi = 2.*TMath::Pi()*Rndm();\r
213         \r
214                 // transform coordenates\r
215                 x = r*TMath::Sin(theta)*TMath::Cos(phi);\r
216                 y = r*TMath::Sin(theta)*TMath::Sin(phi);\r
217                 z = r*TMath::Cos(theta);\r
218         }\r
219         else // projection into the surface of an sphere\r
220         {\r
221                 x = fRsrc*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());\r
222                 y = fRsrc*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());\r
223                 z = fRsrc*TMath::Cos(i->Theta());\r
224         }\r
225         \r
226         // assume nucleons at the freeze-out have the same production vertex\r
227         i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
228 }\r