1 /**************************************************************************
\r
2 * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\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
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
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
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
31 // Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
\r
32 // Arturo Menchaca <menchaca@fisica.unam.mx>
\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
41 #include "TMCProcess.h"
\r
43 #include "TVector3.h"
\r
45 #include "TArrayF.h"
\r
46 #include "AliGenCocktailEventHeader.h"
\r
48 ClassImp(AliGenDeuteron)
\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
59 fSign = sign > 0 ? 1:-1;
\r
62 AliGenDeuteron::~AliGenDeuteron()
\r
66 void AliGenDeuteron::Init()
\r
68 // Standard AliGenerator Initializer
\r
71 void AliGenDeuteron::Generate()
\r
73 // Modify the stack of each event, adding the new particles at the end and
\r
74 // not transporting the parents.
\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
85 // Get the cocktail generator
\r
86 AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
\r
89 for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
\r
91 gener->SetActiveEventNumber(ns);
\r
93 AliStack* stack = gener->GetStack(ns); // Stack of event ns
\r
96 Info("Generate", "no stack, exiting");
\r
100 // find emerging nucleons from the collision of current event
\r
102 TList* protons = new TList();
\r
103 protons->SetOwner(kFALSE);
\r
104 TList* neutrons = new TList();
\r
105 neutrons->SetOwner(kFALSE);
\r
107 // primary vertex of current event
\r
109 gener->GetActiveEventHeader()->PrimaryVertex(primVtx);
\r
110 TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);
\r
112 for (Int_t i=0; i<stack->GetNtrack(); ++i)
\r
114 TParticle* iParticle = stack->Particle(i);
\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
121 Int_t pdgCode = iParticle->GetPdgCode();
\r
122 if(pdgCode == fSign*2212)// (anti)proton
\r
124 FixProductionVertex(iParticle);
\r
125 protons->Add(iParticle);
\r
127 else if(pdgCode == fSign*2112) // (anti)neutron
\r
129 FixProductionVertex(iParticle);
\r
130 neutrons->Add(iParticle);
\r
134 // coalescence conditions
\r
135 // A.J. Baltz et al., Phys. lett B 325(1994)7
\r
137 TIter p_next(protons);
\r
138 while(TParticle* n0 = (TParticle*) p_next())
\r
140 TParticle* iProton = 0;
\r
141 TParticle* jNeutron = 0;
\r
143 // Select the first nucleon
\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
149 // See if there is a neibourgh...
\r
150 TIter n_next(neutrons);
\r
151 while(TParticle* n1 = (TParticle*) n_next() )
\r
153 TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());
\r
154 TVector3 p(n1->Px(), n1->Py(), n1->Pz());
\r
156 if((p-p0).Mag() > fPmax) continue;
\r
157 if((v-v0).Mag() > fRmax) continue;
\r
164 if(jNeutron == 0) continue; // with next proton
\r
166 if(Rndm() > fSpinProb) continue;
\r
168 // neutron captured!
\r
170 TVector3 pN = p0+p1;
\r
171 TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex
\r
174 Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);
\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
182 //Info("Generate","neutron capture (NEW DEUTERON)");
\r
184 // Set bit not to transport for the nucleons
\r
185 iProton->SetBit(kDoneBit);
\r
186 jNeutron->SetBit(kDoneBit);
\r
188 // Remove from the list
\r
189 neutrons->Remove(jNeutron);
\r
192 protons->Clear("nodelete");
\r
193 neutrons->Clear("nodelete");
\r
199 Info("Generate","Coalescence afterburner: DONE");
\r
202 // create the freeze-out nucleon distribution around the collision vertex
\r
203 void AliGenDeuteron::FixProductionVertex(TParticle* i)
\r
207 if(fModel == kThermal) // homogeneous volume
\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
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
219 else // projection into the surface of an sphere
\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
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