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