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