Classed imported from EVGEN.
[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$
18Revision 1.3 2003/02/26 13:08:22 morsch
19TParticle included.
20
21Revision 1.2 2003/02/26 10:26:32 morsch
22AliGenPythiaJets: jets centered on EMCAL
23
24Revision 1.1 2003/01/17 04:10:31 morsch
25First 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
44ClassImp(AliGenPythiaJets)
45
46AliGenPythiaJets::AliGenPythiaJets()
47 :AliGenPythia()
48{
49// Default Constructor
50}
51
52AliGenPythiaJets::AliGenPythiaJets(Int_t npart)
53 :AliGenPythia(npart)
54{
55 fName = "PythiaJets";
56 fTitle= "Jet Generator using PYTHIA";
57}
58
59AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia)
60{
61// copy constructor
62 Pythia.Copy(*this);
63}
64
65AliGenPythiaJets::~AliGenPythiaJets()
66{
67// Destructor
68}
69
70void 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
90void 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
202Bool_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
228AliGenPythiaJets& AliGenPythiaJets::operator=(const AliGenPythiaJets& rhs)
229{
230// Assignment operator
231 return *this;
232}
233
234void 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