Step size region by region.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythiaJets.cxx
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 //
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
34 ClassImp(AliGenPythiaJets)
35
36 AliGenPythiaJets::AliGenPythiaJets()
37                  :AliGenPythia()
38 {
39 // Default Constructor
40 }
41
42 AliGenPythiaJets::AliGenPythiaJets(Int_t npart)
43                  :AliGenPythia(npart)
44 {
45     fName = "PythiaJets";
46     fTitle= "Jet Generator using PYTHIA";
47 }
48
49 AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia):
50     AliGenPythia(Pythia)
51 {
52 // copy constructor
53     Pythia.Copy(*this);
54 }
55
56 AliGenPythiaJets::~AliGenPythiaJets()
57 {
58 // Destructor
59 }
60
61 void 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
81 void 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 //
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);
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();
168                 PushTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
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
193 Bool_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           
219 AliGenPythiaJets& AliGenPythiaJets::operator=(const  AliGenPythiaJets& rhs)
220 {
221 // Assignment operator
222     rhs.Copy(*this);
223     return *this;
224 }
225
226 void  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