]>
Commit | Line | Data |
---|---|---|
1 | /************************************************************************** | |
2 | * Copyright(c) 1998-1999, 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 | /* $Id$ */ | |
17 | ||
18 | // Generator using Hydjet as an external generator | |
19 | // The main Hydjet options are accessable for the user through this interface. | |
20 | // Uses the THydjet implementation of TGenerator. | |
21 | // Author: | |
22 | // Rafael Diaz Valdes (rafael.diaz.valdes@cern.ch) | |
23 | // | |
24 | #include <Riostream.h> | |
25 | #include <THydjet.h> | |
26 | #include <TParticle.h> | |
27 | #include <TClonesArray.h> | |
28 | ||
29 | #include "AliGenHydjet.h" | |
30 | #include "AliLog.h" | |
31 | #include "AliGenHydjetEventHeader.h" | |
32 | #include "AliRun.h" | |
33 | #include "AliPythiaRndm.h" | |
34 | ||
35 | using std::cout; | |
36 | using std::endl; | |
37 | ||
38 | ClassImp(AliGenHydjet) | |
39 | ||
40 | AliGenHydjet::AliGenHydjet(Int_t npart) : | |
41 | AliGenMC(npart), | |
42 | //initial parameters | |
43 | fFrame("CMS"), // Reference frame | |
44 | fAtomicWeigth(207), // Projectile-Target atomic weight | |
45 | fIfbtype(0), // centrality type | |
46 | fFixImpactParam(0), // fixed impact parameter | |
47 | fMinImpactParam(0.), // minimum impact parameter | |
48 | fMaxImpactParam(1.), // maximum impact parameter | |
49 | fSoftMult(20000), // mean soft multiplicity | |
50 | //hydro parameters | |
51 | fJetProd(2), // flag Jet production (nhsel) | |
52 | fYflow(5.), // max longitudinal flow | |
53 | fTflow(1.), // max transverse flow | |
54 | fSoftFract(1.), // Soft multiplicity fraction | |
55 | // PYTHIA parameters | |
56 | fDijetProd(1), // flag dijet production | |
57 | fMinPtHard(10.), // min pt hard | |
58 | fStructFunction(7), // Structure Function (default CTEQ5M) | |
59 | fMultipleInt(0), // flag multiple interaction | |
60 | ||
61 | fHydjet(0), | |
62 | ||
63 | fSelectAll(0), | |
64 | fFlavor(0) | |
65 | ||
66 | { | |
67 | // Default PbPb collisions at 5. 5 TeV | |
68 | // | |
69 | fEnergyCMS = 5500.; //Energy cms | |
70 | fName = "Hydjet"; | |
71 | fTitle = "Particle Generator using Hydjet"; | |
72 | // Set random number generator | |
73 | if (!AliPythiaRndm::GetPythiaRandom()) | |
74 | AliPythiaRndm::SetPythiaRandom(GetRandom()); | |
75 | } | |
76 | ||
77 | AliGenHydjet::~AliGenHydjet() | |
78 | { | |
79 | // Destructor | |
80 | } | |
81 | ||
82 | void AliGenHydjet::Init() | |
83 | { | |
84 | // Initialisation | |
85 | ||
86 | SetMC(new THydjet(fEnergyCMS, fFrame,fAtomicWeigth, | |
87 | fIfbtype,fMinImpactParam, fMaxImpactParam, | |
88 | fFixImpactParam,fSoftMult)); | |
89 | ||
90 | fHydjet=(THydjet*) fMCEvGen; | |
91 | ||
92 | // init hydro common blocks | |
93 | fHydjet->SetNHSEL(fJetProd); | |
94 | fHydjet->SetYLFL(fYflow); | |
95 | fHydjet->SetYTFL(fTflow); | |
96 | fHydjet->SetFPART(fSoftFract); | |
97 | // init PYTHIA common block | |
98 | fHydjet->SetMSEL(fDijetProd); | |
99 | fHydjet->SetPTMIN(fMinPtHard); | |
100 | fHydjet->SetCKIN(3,fMinPtHard); | |
101 | fHydjet->SetMSTP(51,fStructFunction); | |
102 | fHydjet->SetMSTP(81,0); | |
103 | fHydjet->SetMSTU(21,1); | |
104 | fHydjet->SetPARU(14,1.); | |
105 | ||
106 | // Initialize Hydjet | |
107 | fHydjet->Initialize(); | |
108 | ||
109 | } | |
110 | ||
111 | void AliGenHydjet::Generate() | |
112 | { | |
113 | ||
114 | /////////////////////////////////////////////////////////////////////////////////////////// | |
115 | Float_t polar[3] = {0,0,0}; | |
116 | Float_t origin[3] = {0,0,0}; | |
117 | Float_t origin0[3] = {0,0,0}; | |
118 | Float_t time0 = 0.; | |
119 | Float_t p[3]; | |
120 | Float_t tof; | |
121 | ||
122 | // converts from mm/c to s | |
123 | const Float_t kconv = 0.001/2.999792458e8; | |
124 | // | |
125 | Int_t nt = 0; | |
126 | Int_t j, kf, ks, imo; | |
127 | kf = 0; | |
128 | ||
129 | for (j = 0;j < 3; j++) origin0[j] = fOrigin[j]; | |
130 | time0 = fTimeOrigin; | |
131 | if(fVertexSmear == kPerEvent) { | |
132 | Vertex(); | |
133 | for (j=0; j < 3; j++) origin0[j] = fVertex[j]; | |
134 | time0 = fTime; | |
135 | } | |
136 | ||
137 | while(1) | |
138 | { | |
139 | // Generate one event | |
140 | fHydjet->GenerateEvent(); | |
141 | fHydjet->ImportParticles(&fParticles,"All"); | |
142 | ||
143 | Int_t np = fParticles.GetEntriesFast(); | |
144 | Int_t nc = 0; | |
145 | if (np == 0 ) continue; | |
146 | Int_t i; | |
147 | Int_t* pSelected = new Int_t[np]; | |
148 | ||
149 | for (i = 0; i < np; i++) { | |
150 | pSelected[i] = 0; | |
151 | } | |
152 | ||
153 | // Get event vertex | |
154 | fVertex[0] = origin0[0]; | |
155 | fVertex[1] = origin0[1]; | |
156 | fVertex[2] = origin0[2]; | |
157 | fTime = time0; | |
158 | ||
159 | // | |
160 | // Now select the final state particles | |
161 | // | |
162 | ||
163 | for(i = 0; i<np; i++){ | |
164 | TParticle * iparticle = (TParticle *) fParticles.At(i); | |
165 | // Is this a final state particle ? | |
166 | if (!Stable(iparticle)) continue; | |
167 | Bool_t selected = kTRUE; | |
168 | kf = iparticle->GetPdgCode(); | |
169 | ks = iparticle->GetStatusCode(); | |
170 | if(!fSelectAll) selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf); | |
171 | // Put particle on the stack if selected | |
172 | if(selected) { | |
173 | nc++; | |
174 | pSelected[i] = 1; | |
175 | } // selected | |
176 | } // particle loop final state | |
177 | ||
178 | ||
179 | // | |
180 | // Write particles to stack | |
181 | // | |
182 | ||
183 | for(i = 0; i<np; i++) { | |
184 | TParticle * iparticle = (TParticle *) fParticles.At(i); | |
185 | Bool_t hasDaughter = (iparticle->GetFirstDaughter() > 0); | |
186 | if(pSelected[i]){ | |
187 | kf = iparticle->GetPdgCode(); | |
188 | ks = iparticle->GetStatusCode(); | |
189 | p[0] = iparticle->Px(); | |
190 | p[1] = iparticle->Py(); | |
191 | p[2] = iparticle->Pz(); | |
192 | origin[0] = fVertex[0]+iparticle->Vx()/10; | |
193 | origin[1] = fVertex[1]+iparticle->Vy()/10; | |
194 | origin[2] = fVertex[2]+iparticle->Vz()/10; | |
195 | tof = fTime + kconv*iparticle->T(); | |
196 | ||
197 | imo = -1; | |
198 | //imo = iparticle->GetFirstMother(); | |
199 | //if(imo == 0) imo =-1; | |
200 | Bool_t tFlag = (fTrackIt && !hasDaughter); | |
201 | PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks); | |
202 | KeepTrack(nt); | |
203 | } // if selected | |
204 | } // particle loop | |
205 | delete[] pSelected; | |
206 | ||
207 | AliInfo(Form("\n I've put %i particles on the stack \n",nc)); | |
208 | if (nc > 0) break; | |
209 | } // event loop | |
210 | MakeHeader(); | |
211 | } | |
212 | ||
213 | ||
214 | Bool_t AliGenHydjet::SelectFlavor(Int_t pid) | |
215 | { | |
216 | // Select flavor of particle | |
217 | // 0: all | |
218 | // 4: charm and beauty | |
219 | // 5: beauty | |
220 | Bool_t res = 0; | |
221 | if (fFlavor == 0) { | |
222 | res = kTRUE; | |
223 | } else { | |
224 | Int_t ifl = TMath::Abs(pid/100); | |
225 | if (ifl > 10) ifl/=10; | |
226 | res = (fFlavor == ifl); | |
227 | } | |
228 | return res; | |
229 | } | |
230 | ||
231 | Bool_t AliGenHydjet::Stable(TParticle* particle) const | |
232 | { | |
233 | // Return true for a stable particle | |
234 | if (particle->GetFirstDaughter() == 0 ) | |
235 | { | |
236 | return kTRUE; | |
237 | } else { | |
238 | return kFALSE; | |
239 | } | |
240 | } | |
241 | ||
242 | ||
243 | /////////////////////////////////// | |
244 | void AliGenHydjet::MakeHeader() | |
245 | { | |
246 | // Builds the event header, to be called after each event | |
247 | AliGenEventHeader* header = new AliGenHydjetEventHeader("Hydjet"); | |
248 | ((AliGenHydjetEventHeader*) header)->SetNProduced(fHydjet->GetN()); | |
249 | ((AliGenHydjetEventHeader*) header)->SetImpactParameter(fHydjet->GetBGEN()); | |
250 | ((AliGenHydjetEventHeader*) header)->SetNbcol(fHydjet->GetNBCOL()); | |
251 | ((AliGenHydjetEventHeader*) header)->SetNpart(fHydjet->GetNPART()); | |
252 | ((AliGenHydjetEventHeader*) header)->SetNpyt(fHydjet->GetNPYT()); | |
253 | ((AliGenHydjetEventHeader*) header)->SetNhyd(fHydjet->GetNHYD()); | |
254 | // draw generation values | |
255 | cout << "*********Hyflow ********** " << endl; | |
256 | cout << " YTFL " << fHydjet->GetYTFL() << endl; | |
257 | cout << " YLFL " << fHydjet->GetYLFL() << endl; | |
258 | cout << " FPART " << fHydjet->GetFPART() << endl; | |
259 | cout << "*********Hyjpar ********* " << endl; | |
260 | cout << " GetNHSEL() " << fHydjet->GetNHSEL() << endl; | |
261 | cout << " GetPTMIN() " << fHydjet->GetPTMIN() << endl; | |
262 | cout << " GetNJET() " << fHydjet->GetNJET() << endl; | |
263 | cout << "*********Hyfpar ********* " << endl; | |
264 | cout << " GetBGEN() " << fHydjet->GetBGEN() << endl; | |
265 | cout << " GetNBCOL() " << fHydjet->GetNBCOL() << endl; | |
266 | cout << " GetNPART() " << fHydjet->GetNPART() << endl; | |
267 | cout << " GetNPYT() " << fHydjet->GetNPYT() << endl; | |
268 | cout << " GetNHYD() " << fHydjet->GetNHYD() << endl; | |
269 | ||
270 | ||
271 | // Event Vertex | |
272 | header->SetPrimaryVertex(fVertex); | |
273 | header->SetInteractionTime(fTime); | |
274 | AddHeader(header); | |
275 | fCollisionGeometry = (AliGenHydjetEventHeader*) header; | |
276 | } | |
277 | ||
278 | void AliGenHydjet::AddHeader(AliGenEventHeader* header) | |
279 | { | |
280 | // Passes header either to the container or to gAlice | |
281 | if (fContainer) { | |
282 | fContainer->AddHeader(header); | |
283 | } else { | |
284 | gAlice->SetGenEventHeader(header); | |
285 | } | |
286 | } | |
287 | ||
288 | ||
289 | void AliGenHydjet::Copy(TObject &) const | |
290 | { | |
291 | Fatal("Copy","Not implemented!\n"); | |
292 | } | |
293 | ||
294 | AliGenHydjet& AliGenHydjet::operator=(const AliGenHydjet& rhs) | |
295 | { | |
296 | rhs.Copy(*this); | |
297 | return (*this); | |
298 | } | |
299 |