]>
Commit | Line | Data |
---|---|---|
8d2cd130 | 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 | /* | |
17 | $Log$ | |
18 | Revision 1.3 2003/02/26 13:08:22 morsch | |
19 | TParticle included. | |
20 | ||
21 | Revision 1.2 2003/02/26 10:26:32 morsch | |
22 | AliGenPythiaJets: jets centered on EMCAL | |
23 | ||
24 | Revision 1.1 2003/01/17 04:10:31 morsch | |
25 | First commit. | |
26 | */ | |
27 | ||
28 | // | |
29 | // Generator using the TPythia interface (via AliPythia) | |
30 | // to generate jets in pp collisions. | |
31 | // Using SetNuclei() also nuclear modifications to the structure functions | |
32 | // can be taken into account. | |
33 | // Using SetQuenchingFactor(f) quenched jets can be modelled by superimposing | |
34 | // two jets with energies e * f and e * (1-f) | |
35 | // | |
36 | // andreas.morsch@cern.ch | |
37 | // | |
38 | ||
39 | ||
40 | #include "AliGenPythiaJets.h" | |
41 | #include "AliRun.h" | |
42 | #include <TParticle.h> | |
43 | ||
44 | ClassImp(AliGenPythiaJets) | |
45 | ||
46 | AliGenPythiaJets::AliGenPythiaJets() | |
47 | :AliGenPythia() | |
48 | { | |
49 | // Default Constructor | |
50 | } | |
51 | ||
52 | AliGenPythiaJets::AliGenPythiaJets(Int_t npart) | |
53 | :AliGenPythia(npart) | |
54 | { | |
55 | fName = "PythiaJets"; | |
56 | fTitle= "Jet Generator using PYTHIA"; | |
57 | } | |
58 | ||
59 | AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia) | |
60 | { | |
61 | // copy constructor | |
62 | Pythia.Copy(*this); | |
63 | } | |
64 | ||
65 | AliGenPythiaJets::~AliGenPythiaJets() | |
66 | { | |
67 | // Destructor | |
68 | } | |
69 | ||
70 | void AliGenPythiaJets::Init() | |
71 | { | |
72 | // Initialization | |
73 | // | |
74 | printf("AliGenPythiaJets::Init() \n"); | |
75 | ||
76 | AliGenPythia::Init(); | |
77 | ||
78 | if (fQuench > 0.) { | |
79 | fEtMinJetQ[0] = fEtMinJet * fQuench; | |
80 | fEtMaxJetQ[0] = fEtMaxJet * fQuench; | |
81 | fEtMinJetQ[1] = fEtMinJet * (1. - fQuench); | |
82 | fEtMaxJetQ[1] = fEtMaxJet * (1. - fQuench); | |
83 | fPtHardMinQ[0] = fPtHardMin * fQuench; | |
84 | fPtHardMaxQ[0] = fPtHardMax * fQuench; | |
85 | fPtHardMinQ[1] = fPtHardMin * (1. - fQuench); | |
86 | fPtHardMaxQ[1] = fPtHardMax * (1. - fQuench); | |
87 | } | |
88 | } | |
89 | ||
90 | void AliGenPythiaJets::Generate() | |
91 | { | |
92 | // Generate one event | |
93 | fDecayer->ForceDecay(); | |
94 | ||
95 | Float_t polar[3] = {0,0,0}; | |
96 | Float_t origin[3] = {0,0,0}; | |
97 | Float_t p[3]; | |
98 | // converts from mm/c to s | |
99 | const Float_t kconv=0.001/2.999792458e8; | |
100 | // | |
101 | Int_t nt = 0; | |
102 | Int_t nc = 0; | |
103 | Int_t jev = 0; | |
104 | Int_t j, kf, iparent; | |
105 | fTrials=0; | |
106 | // | |
107 | // Set collision vertex position | |
108 | if(fVertexSmear==kPerEvent) { | |
109 | fPythia->SetMSTP(151,1); | |
110 | for (j=0;j<3;j++) { | |
111 | fPythia->SetPARP(151+j, fOsigma[j]*10.); | |
112 | } | |
113 | } else if (fVertexSmear==kPerTrack) { | |
114 | fPythia->SetMSTP(151,0); | |
115 | } | |
116 | // Event loop | |
117 | while(1) | |
118 | { | |
119 | if (fQuench > 0.) { | |
120 | fPythia->SetCKIN(3,fPtHardMinQ[jev]); | |
121 | fPythia->SetCKIN(4,fPtHardMaxQ[jev]); | |
122 | fEtMinJet = fEtMinJetQ[jev]; | |
123 | fEtMaxJet = fEtMaxJetQ[jev]; | |
124 | } | |
125 | ||
126 | fPythia->Pyevnt(); | |
127 | if (gAlice->GetEvNumber()>=fDebugEventFirst && | |
128 | gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1); | |
129 | fTrials++; | |
130 | // | |
131 | // Has this jet triggered | |
132 | // | |
133 | if ((fEtMinJet != -1) && ! CheckTrigger()) continue; | |
134 | // | |
135 | fPythia->ImportParticles(fParticles,"All"); | |
136 | Int_t i; | |
137 | Int_t np = fParticles->GetEntriesFast(); | |
138 | if (np == 0 ) continue; | |
139 | // Get event vertex and discard the event if the z coord. is too big | |
140 | TParticle *iparticle = (TParticle *) fParticles->At(0); | |
141 | Float_t distz = iparticle->Vz()/10.; | |
142 | if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue; | |
143 | // | |
144 | // | |
145 | fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0); | |
146 | fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1); | |
147 | fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2); | |
148 | ||
149 | Int_t* pParent = new Int_t[np]; | |
150 | for (i=0; i< np; i++) pParent[i] = -1; | |
151 | ||
152 | ||
153 | // | |
154 | for (i = 0; i<np; i++) { | |
155 | Int_t trackIt = 0; | |
156 | TParticle * iparticle = (TParticle *) fParticles->At(i); | |
157 | kf = CheckPDGCode(iparticle->GetPdgCode()); | |
158 | Int_t ks = iparticle->GetStatusCode(); | |
159 | Int_t km = iparticle->GetFirstMother(); | |
160 | if ((ks == 1 && kf !=0 && KinematicSelection(iparticle, 0)) || | |
161 | (ks != 1) || | |
162 | (fProcess == kPyJets && ks == 21 && km == 0 && i > 1)) { | |
163 | nc++; | |
164 | if (ks == 1) trackIt = 1; | |
165 | Int_t ipa = iparticle->GetFirstMother() - 1; | |
166 | iparent = (ipa > -1) ? pParent[ipa] : -1; | |
167 | // | |
168 | // Store track information | |
169 | // | |
170 | p[0] = iparticle->Px(); | |
171 | p[1] = iparticle->Py(); | |
172 | p[2] = iparticle->Pz(); | |
173 | origin[0] = fOrigin[0]+iparticle->Vx()/10.; | |
174 | origin[1] = fOrigin[1]+iparticle->Vy()/10.; | |
175 | origin[2] = fOrigin[2]+iparticle->Vz()/10.; | |
176 | Float_t tof=kconv*iparticle->T(); | |
177 | SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar, | |
178 | tof, kPPrimary, nt, 1., ks); | |
179 | KeepTrack(nt); | |
180 | pParent[i] = nt; | |
181 | } // select particle | |
182 | } // particle loop | |
183 | ||
184 | if (pParent) delete[] pParent; | |
185 | printf("\n AliGenPythiaJets: I've put %i particles on the stack \n",nc); | |
186 | if (nc > 0) { | |
187 | jev += 1; | |
188 | if ((fQuench <= 0.) || (fQuench > 0. && jev == 2)) { | |
189 | fKineBias=Float_t(fNpart)/Float_t(fTrials); | |
190 | printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev); | |
191 | fNev++; | |
192 | MakeHeader(); | |
193 | break; | |
194 | } | |
195 | } | |
196 | } | |
197 | SetHighWaterMark(nt); | |
198 | // Get cross-section | |
199 | fXsection=fPythia->GetPARI(1); | |
200 | } | |
201 | ||
202 | Bool_t AliGenPythiaJets::CheckTrigger() | |
203 | { | |
204 | // Check the kinematic trigger condition | |
205 | // | |
206 | Bool_t triggered = kFALSE; | |
207 | Int_t njets = 0; | |
208 | Int_t ntrig = 0; | |
209 | Float_t jets[4][10]; | |
210 | // | |
211 | // Use Pythia clustering on parton level to determine jet axis | |
212 | // | |
213 | GetJets(njets, ntrig, jets); | |
214 | ||
215 | if (ntrig) { | |
216 | triggered = kTRUE; | |
217 | Float_t px = jets[0][0]; | |
218 | Float_t py = jets[1][0]; | |
219 | Float_t pz = jets[2][0]; | |
220 | Float_t e = jets[3][0]; | |
221 | Float_t beta = pz/e; | |
222 | Float_t phi = TMath::Pi()+TMath::ATan2(-py,-px); | |
223 | TransformEvent(beta, -2. * TMath::Pi() / 3. + phi); | |
224 | } | |
225 | return triggered; | |
226 | } | |
227 | ||
228 | AliGenPythiaJets& AliGenPythiaJets::operator=(const AliGenPythiaJets& rhs) | |
229 | { | |
230 | // Assignment operator | |
231 | return *this; | |
232 | } | |
233 | ||
234 | void AliGenPythiaJets::TransformEvent(Float_t beta, Float_t phi) | |
235 | { | |
236 | // | |
237 | // Perform Lorentz Transformation and Rotation | |
238 | // | |
239 | Float_t gamma = 1./TMath::Sqrt(1. - beta * beta); | |
240 | Int_t npart = (fPythia->GetPyjets())->N; | |
241 | ||
242 | for (Int_t part = 0; part < npart; part++) { | |
243 | Float_t px = (fPythia->GetPyjets())->P[0][part]; | |
244 | Float_t py = (fPythia->GetPyjets())->P[1][part]; | |
245 | Float_t pz = (fPythia->GetPyjets())->P[2][part]; | |
246 | Float_t e = (fPythia->GetPyjets())->P[3][part]; | |
247 | // | |
248 | // Lorentz Transform | |
249 | // | |
250 | Float_t pzt = gamma * pz - gamma * beta * e; | |
251 | Float_t et = -gamma * beta * pz + gamma * e; | |
252 | // | |
253 | // Rotation | |
254 | // | |
255 | Float_t pxt = TMath::Cos(phi) * px + TMath::Sin(phi) * py; | |
256 | Float_t pyt = - TMath::Sin(phi) * px + TMath::Cos(phi) * py; | |
257 | // | |
258 | // | |
259 | (fPythia->GetPyjets())->P[0][part] = pxt; | |
260 | (fPythia->GetPyjets())->P[1][part] = pyt; | |
261 | (fPythia->GetPyjets())->P[2][part] = pzt; | |
262 | (fPythia->GetPyjets())->P[3][part] = et; | |
263 | } | |
264 | } | |
265 |