]>
Commit | Line | Data |
---|---|---|
e3b21c34 | 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 | |
539d3694 | 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 | |
539d3694 | 44 | #include "AliMC.h"\r |
e3b21c34 | 45 | #include "TArrayF.h"\r |
46 | #include "AliGenCocktailEventHeader.h"\r | |
539d3694 | 47 | \r |
48 | ClassImp(AliGenDeuteron)\r | |
49 | \r | |
539d3694 | 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 | |
539d3694 | 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 | |
e3b21c34 | 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 | |
539d3694 | 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 | |
e3b21c34 | 124 | FixProductionVertex(iParticle);\r |
539d3694 | 125 | protons->Add(iParticle);\r |
126 | }\r | |
127 | else if(pdgCode == fSign*2112) // (anti)neutron\r | |
128 | {\r | |
e3b21c34 | 129 | FixProductionVertex(iParticle);\r |
539d3694 | 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 | |
e3b21c34 | 166 | if(Rndm() > fSpinProb) continue;\r |
539d3694 | 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 | |
e3b21c34 | 203 | void AliGenDeuteron::FixProductionVertex(TParticle* i)\r |
539d3694 | 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 | |
e3b21c34 | 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 | |
539d3694 | 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 | |
e3b21c34 | 226 | // assume nucleons at the freeze-out have the same production vertex\r |
539d3694 | 227 | i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r |
228 | }\r |