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