3980c287e1ba3a08135951cb087f54fdf5cd4121
[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 /*
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