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