Change libdummypythia to libdummypythia6
[u/mrichter/AliRoot.git] / EVGEN / AliGenPythia.cxx
CommitLineData
4c039060 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$
00d6ce7d 18Revision 1.20 2000/09/06 14:29:33 morsch
19Use AliPythia for event generation an AliDecayPythia for decays.
20Correct handling of "nodecay" option
21
18edb254 22Revision 1.19 2000/07/11 18:24:56 fca
23Coding convention corrections + few minor bug fixes
24
aee8290b 25Revision 1.18 2000/06/30 12:40:34 morsch
26Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
27Geant units (cm).
28
1512e357 29Revision 1.17 2000/06/09 20:34:07 morsch
30All coding rule violations except RS3 corrected
31
f87cfe57 32Revision 1.16 2000/05/15 15:04:20 morsch
33The full event is written for fNtrack = -1
34Coding rule violations corrected.
35
5ddeb374 36Revision 1.15 2000/04/26 10:14:24 morsch
37Particles array has one entry more than pythia particle list. Upper bound of
38particle loop changed to np-1 (R. Guernane, AM)
39
2a0e6f5b 40Revision 1.14 2000/04/05 08:36:13 morsch
41Check status code of particles in Pythia event
42to avoid double counting as partonic state and final state particle.
43
23211d3e 44Revision 1.13 1999/11/09 07:38:48 fca
45Changes for compatibility with version 2.23 of ROOT
46
084c1b4a 47Revision 1.12 1999/11/03 17:43:20 fca
48New version from G.Martinez & A.Morsch
49
886b6f73 50Revision 1.11 1999/09/29 09:24:14 fca
51Introduction of the Copyright and cvs Log
4c039060 52*/
53
fe4da5cc 54#include "AliGenPythia.h"
18edb254 55#include "AliDecayerPythia.h"
fe4da5cc 56#include "AliRun.h"
57#include "AliPythia.h"
18edb254 58#include "AliPDG.h"
1578254f 59#include <TParticle.h>
18edb254 60#include <TSystem.h>
f87cfe57 61
fe4da5cc 62 ClassImp(AliGenPythia)
63
64AliGenPythia::AliGenPythia()
65 :AliGenerator()
66{
18edb254 67// Default Constructor
68 fDecayer = new AliDecayerPythia();
fe4da5cc 69}
70
71AliGenPythia::AliGenPythia(Int_t npart)
72 :AliGenerator(npart)
73{
74// default charm production at 5. 5 TeV
75// semimuonic decay
76// structure function GRVHO
77//
78 fXsection = 0.;
79 fParentSelect.Set(5);
80 fChildSelect.Set(5);
81 for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
82 SetProcess();
83 SetStrucFunc();
886b6f73 84 SetForceDecay();
fe4da5cc 85 SetPtHard();
86 SetEnergyCMS();
18edb254 87 fDecayer = new AliDecayerPythia();
fe4da5cc 88}
89
f87cfe57 90AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
91{
92// copy constructor
93}
94
fe4da5cc 95AliGenPythia::~AliGenPythia()
96{
5ddeb374 97// Destructor
fe4da5cc 98}
99
100void AliGenPythia::Init()
101{
5ddeb374 102// Initialisation
18edb254 103 SetMC(AliPythia::Instance());
fe4da5cc 104 fPythia=(AliPythia*) fgMCEvGen;
105//
106 fParentWeight=1./Float_t(fNpart);
107//
108// Forward Paramters to the AliPythia object
18edb254 109 // gSystem->Exec("ln -s $ALICE_ROOT/data/Decay.table fort.1");
110 // fPythia->Pyupda(2,1);
111 // gSystem->Exec("rm fort.1");
00d6ce7d 112
113 fDecayer->SetForceDecay(fForceDecay);
18edb254 114 fDecayer->Init();
00d6ce7d 115
18edb254 116
fe4da5cc 117 fPythia->SetCKIN(3,fPtHardMin);
118 fPythia->SetCKIN(4,fPtHardMax);
119 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
18edb254 120
121 // fPythia->Pylist(0);
122 // fPythia->Pystat(2);
fe4da5cc 123// Parent and Children Selection
124 switch (fProcess)
125 {
126 case charm:
127
128 fParentSelect[0]=411;
129 fParentSelect[1]=421;
130 fParentSelect[2]=431;
131 fParentSelect[3]=4122;
132 break;
133 case charm_unforced:
134
135 fParentSelect[0]=411;
136 fParentSelect[1]=421;
137 fParentSelect[2]=431;
138 fParentSelect[3]=4122;
139 break;
140 case beauty:
141 fParentSelect[0]=511;
142 fParentSelect[1]=521;
143 fParentSelect[2]=531;
144 fParentSelect[3]=5122;
145 break;
146 case beauty_unforced:
147 fParentSelect[0]=511;
148 fParentSelect[1]=521;
149 fParentSelect[2]=531;
150 fParentSelect[3]=5122;
151 break;
152 case jpsi_chi:
153 case jpsi:
154 fParentSelect[0]=443;
155 break;
5be1fe76 156 case mb:
157 break;
fe4da5cc 158 }
159
160 switch (fForceDecay)
161 {
162 case semielectronic:
163 case dielectron:
164 case b_jpsi_dielectron:
165 case b_psip_dielectron:
18edb254 166 fChildSelect[0]=kElectron;
fe4da5cc 167 break;
168 case semimuonic:
169 case dimuon:
170 case b_jpsi_dimuon:
171 case b_psip_dimuon:
5be1fe76 172 case pitomu:
173 case katomu:
18edb254 174 fChildSelect[0]=kMuonMinus;
fe4da5cc 175 break;
18edb254 176 case hadronicD:
177 fChildSelect[0]=kPiPlus;
178 fChildSelect[1]=kKPlus;
179 break;
886b6f73 180 case all:
181 case nodecay:
18edb254 182 break;
fe4da5cc 183 }
184}
185
186void AliGenPythia::Generate()
187{
5ddeb374 188// Generate one event
00d6ce7d 189 fDecayer->ForceDecay();
fe4da5cc 190
5be1fe76 191 Float_t polar[3] = {0,0,0};
192 Float_t origin[3]= {0,0,0};
5ddeb374 193 Float_t originP[3]= {0,0,0};
5be1fe76 194 Float_t origin0[3]= {0,0,0};
1512e357 195 Float_t p[3], pP[4];
196// Float_t random[6];
09fd3ea2 197 static TClonesArray *particles;
5be1fe76 198// converts from mm/c to s
199 const Float_t kconv=0.001/2.999792458e8;
fe4da5cc 200
201
202//
fe4da5cc 203 Int_t nt=0;
5ddeb374 204 Int_t ntP=0;
fe4da5cc 205 Int_t jev=0;
7921755b 206 Int_t j, kf;
09fd3ea2 207
208 if(!particles) particles=new TClonesArray("TParticle",1000);
fe4da5cc 209
210 fTrials=0;
211 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
aee8290b 212 if(fVertexSmear==kPerEvent) {
1512e357 213 fPythia->SetMSTP(151,1);
fe4da5cc 214 for (j=0;j<3;j++) {
1512e357 215 fPythia->SetPARP(151+j, fOsigma[j]/10.);
fe4da5cc 216 }
aee8290b 217 } else if (fVertexSmear==kPerTrack) {
fe4da5cc 218 fPythia->SetMSTP(151,0);
fe4da5cc 219 }
1512e357 220
fe4da5cc 221 while(1)
222 {
084c1b4a 223 fPythia->Pyevnt();
18edb254 224// fPythia->Pylist(1);
fe4da5cc 225 fTrials++;
f1092809 226 fPythia->ImportParticles(particles,"All");
fe4da5cc 227 Int_t np = particles->GetEntriesFast();
228 printf("\n **************************************************%d\n",np);
229 Int_t nc=0;
230 if (np == 0 ) continue;
5be1fe76 231 if (fProcess != mb) {
2a0e6f5b 232 for (Int_t i = 0; i<np-1; i++) {
1578254f 233 TParticle * iparticle = (TParticle *) particles->At(i);
23211d3e 234 Int_t ks = iparticle->GetStatusCode();
a8228d85 235 kf = CheckPDGCode(iparticle->GetPdgCode());
2a0e6f5b 236 if (ks==21) continue;
237
18edb254 238 fChildWeight=(fDecayer->GetPartialBranchingRatio(kf))*fParentWeight;
fe4da5cc 239//
240// Parent
5be1fe76 241 if (ParentSelected(TMath::Abs(kf))) {
242 if (KinematicSelection(iparticle)) {
243 if (nc==0) {
244//
245// Store information concerning the hard scattering process
246//
5ddeb374 247 Float_t massP = fPythia->GetPARI(13);
248 Float_t ptP = fPythia->GetPARI(17);
249 Float_t yP = fPythia->GetPARI(37);
250 Float_t xmtP = sqrt(ptP*ptP+massP*massP);
251 Float_t ty = Float_t(TMath::TanH(yP));
252 pP[0] = ptP;
253 pP[1] = 0;
254 pP[2] = xmtP*ty/sqrt(1.-ty*ty);
255 pP[3] = massP;
5be1fe76 256 gAlice->SetTrack(0,-1,-1,
5ddeb374 257 pP,originP,polar,
258 0,"Hard Scat.",ntP,fParentWeight);
259 gAlice->KeepTrack(ntP);
5be1fe76 260 }
261 nc++;
fe4da5cc 262//
263// store parent track information
1578254f 264 p[0]=iparticle->Px();
265 p[1]=iparticle->Py();
266 p[2]=iparticle->Pz();
1512e357 267 origin[0]=origin0[0]+iparticle->Vx()*10.;
268 origin[1]=origin0[1]+iparticle->Vy()*10.;
269 origin[2]=origin0[2]+iparticle->Vz()*10.;
fe4da5cc 270
1578254f 271 Int_t ifch=iparticle->GetFirstDaughter();
272 Int_t ilch=iparticle->GetLastDaughter();
18edb254 273
274 if ((ifch !=0 && ilch !=0) || fForceDecay == nodecay) {
275 Int_t trackit=0;
276 if (fForceDecay == nodecay) trackit = 1;
277 gAlice->SetTrack(trackit,ntP,kf,
5be1fe76 278 p,origin,polar,
279 0,"Primary",nt,fParentWeight);
280 gAlice->KeepTrack(nt);
281 Int_t iparent = nt;
fe4da5cc 282//
283// Children
18edb254 284 if (fForceDecay != nodecay) {
285 for (j=ifch; j<=ilch; j++)
286 {
287 TParticle * ichild =
1578254f 288 (TParticle *) particles->At(j-1);
18edb254 289 kf = CheckPDGCode(ichild->GetPdgCode());
fe4da5cc 290//
291//
18edb254 292 if (ChildSelected(TMath::Abs(kf))) {
1512e357 293 origin[0]=origin0[0]+ichild->Vx()*10.;
294 origin[1]=origin0[1]+ichild->Vy()*10.;
295 origin[2]=origin0[2]+ichild->Vz()*10.;
1578254f 296 p[0]=ichild->Px();
297 p[1]=ichild->Py();
298 p[2]=ichild->Pz();
299 Float_t tof=kconv*ichild->T();
300 gAlice->SetTrack(fTrackIt, iparent, kf,
5be1fe76 301 p,origin,polar,
302 tof,"Decay",nt,fChildWeight);
303 gAlice->KeepTrack(nt);
18edb254 304 } // select child
305 } // child loop
306 }
5be1fe76 307 }
308 } // kinematic selection
309 } // select particle
310 } // particle loop
311 } else {
2a0e6f5b 312 for (Int_t i = 0; i<np-1; i++) {
1578254f 313 TParticle * iparticle = (TParticle *) particles->At(i);
a8228d85 314 kf = CheckPDGCode(iparticle->GetPdgCode());
1578254f 315 Int_t ks = iparticle->GetStatusCode();
2a0e6f5b 316
1578254f 317 if (ks==1 && kf!=0 && KinematicSelection(iparticle)) {
5be1fe76 318 nc++;
319//
320// store track information
1578254f 321 p[0]=iparticle->Px();
322 p[1]=iparticle->Py();
323 p[2]=iparticle->Pz();
1512e357 324 origin[0]=origin0[0]+iparticle->Vx()*10;
325 origin[1]=origin0[1]+iparticle->Vy()*10;
326 origin[2]=origin0[2]+iparticle->Vz()*10;
1578254f 327 Float_t tof=kconv*iparticle->T();
328 gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
5be1fe76 329 tof,"Primary",nt);
330 gAlice->KeepTrack(nt);
331 } // select particle
332 } // particle loop
333 printf("\n I've put %i particles on the stack \n",nc);
334 } // mb ?
fe4da5cc 335 if (nc > 0) {
5be1fe76 336 jev+=nc;
5ddeb374 337 if (jev >= fNpart || fNpart == -1) {
fe4da5cc 338 fKineBias=Float_t(fNpart)/Float_t(fTrials);
23211d3e 339 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
fe4da5cc 340 break;
341 }
342 }
343 } // event loop
344// adjust weight due to kinematic selection
345 AdjustWeights();
346// get cross-section
347 fXsection=fPythia->GetPARI(1);
348}
349
350Bool_t AliGenPythia::ParentSelected(Int_t ip)
351{
5ddeb374 352// True if particle is in list of parent particles to be selected
fe4da5cc 353 for (Int_t i=0; i<5; i++)
354 {
355 if (fParentSelect[i]==ip) return kTRUE;
356 }
357 return kFALSE;
358}
359
360Bool_t AliGenPythia::ChildSelected(Int_t ip)
361{
5ddeb374 362// True if particle is in list of decay products to be selected
fe4da5cc 363 for (Int_t i=0; i<5; i++)
364 {
365 if (fChildSelect[i]==ip) return kTRUE;
366 }
367 return kFALSE;
368}
369
1578254f 370Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
fe4da5cc 371{
5ddeb374 372// Perform kinematic selection
1578254f 373 Float_t px=particle->Px();
374 Float_t py=particle->Py();
375 Float_t pz=particle->Pz();
376 Float_t e=particle->Energy();
fe4da5cc 377
378//
379// transverse momentum cut
380 Float_t pt=TMath::Sqrt(px*px+py*py);
381 if (pt > fPtMax || pt < fPtMin)
382 {
383// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
384 return kFALSE;
385 }
386//
387// momentum cut
388 Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
389 if (p > fPMax || p < fPMin)
390 {
391// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
392 return kFALSE;
393 }
394
395//
396// theta cut
5be1fe76 397 Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
fe4da5cc 398 if (theta > fThetaMax || theta < fThetaMin)
399 {
400// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
401 return kFALSE;
402 }
403
404//
405// rapidity cut
406 Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
407 if (y > fYMax || y < fYMin)
408 {
409// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
410 return kFALSE;
411 }
412
413//
414// phi cut
f87cfe57 415 Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
fe4da5cc 416 if (phi > fPhiMax || phi < fPhiMin)
417 {
418// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
419 return kFALSE;
420 }
421
422 return kTRUE;
423}
424void AliGenPythia::AdjustWeights()
425{
5ddeb374 426// Adjust the weights after generation of all events
427//
428 TClonesArray *partArray = gAlice->Particles();
429 TParticle *part;
fe4da5cc 430 Int_t ntrack=gAlice->GetNtrack();
431 for (Int_t i=0; i<ntrack; i++) {
5ddeb374 432 part= (TParticle*) partArray->UncheckedAt(i);
433 part->SetWeight(part->GetWeight()*fKineBias);
fe4da5cc 434 }
435}
436
a8228d85 437Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
438{
439//
23211d3e 440// If the particle is in a diffractive state, then take action accordingly
a8228d85 441 switch (pdgcode) {
18edb254 442 case 91:
443 return 92;
a8228d85 444 case 110:
445 //rho_diff0 -- difficult to translate, return rho0
446 return 113;
447 case 210:
448 //pi_diffr+ -- change to pi+
449 return 211;
450 case 220:
451 //omega_di0 -- change to omega0
452 return 223;
453 case 330:
454 //phi_diff0 -- return phi0
455 return 333;
456 case 440:
457 //J/psi_di0 -- return J/psi
458 return 443;
459 case 2110:
460 //n_diffr -- return neutron
461 return 2112;
462 case 2210:
463 //p_diffr+ -- return proton
464 return 2212;
465 }
466 //non diffractive state -- return code unchanged
467 return pdgcode;
468}
469
f87cfe57 470AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
471{
472// Assignment operator
473 return *this;
474}
fe4da5cc 475
476
18edb254 477void AliGenPythia::Streamer(TBuffer &R__b)
478{
479 // Stream an object of class AliGenPythia.
480
481 if (R__b.IsReading()) {
482 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
483 AliGenerator::Streamer(R__b);
484 R__b >> (Int_t&)fProcess;
485 R__b >> (Int_t&)fStrucFunc;
486 R__b >> (Int_t&)fForceDecay;
487 R__b >> fEnergyCMS;
488 R__b >> fKineBias;
489 R__b >> fTrials;
490 fParentSelect.Streamer(R__b);
491 fChildSelect.Streamer(R__b);
492 R__b >> fXsection;
493// (AliPythia::Instance())->Streamer(R__b);
494 R__b >> fPtHardMin;
495 R__b >> fPtHardMax;
496// if (fDecayer) fDecayer->Streamer(R__b);
497 } else {
498 R__b.WriteVersion(AliGenPythia::IsA());
499 AliGenerator::Streamer(R__b);
500 R__b << (Int_t)fProcess;
501 R__b << (Int_t)fStrucFunc;
502 R__b << (Int_t)fForceDecay;
503 R__b << fEnergyCMS;
504 R__b << fKineBias;
505 R__b << fTrials;
506 fParentSelect.Streamer(R__b);
507 fChildSelect.Streamer(R__b);
508 R__b << fXsection;
509// R__b << fPythia;
510 R__b << fPtHardMin;
511 R__b << fPtHardMax;
512 // fDecayer->Streamer(R__b);
513 }
514}
515
516
517
fe4da5cc 518
5be1fe76 519
520
521
522