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