]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenPythia.cxx
A pointer was set to zero in the default constructor to avoid memory management problems
[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$
4c07486d 18Revision 1.41 2001/10/08 08:45:42 morsch
19Possibility of vertex cut added.
20
a0b2b296 21Revision 1.40 2001/09/25 11:30:23 morsch
22Pass event vertex to header.
23
c5d36e84 24Revision 1.39 2001/07/27 17:09:36 morsch
25Use local SetTrack, KeepTrack and SetHighWaterMark methods
26to delegate either to local stack or to stack owned by AliRun.
27(Piotr Skowronski, A.M.)
28
a99cf51f 29Revision 1.38 2001/07/13 10:58:54 morsch
30- Some coded moved to AliGenMC
31- Improved handling of secondary vertices.
32
fff02fee 33Revision 1.37 2001/06/28 11:17:28 morsch
34SetEventListRange setter added. Events in specified range are listed for
35debugging. (Yuri Kharlov)
36
11dfaf8d 37Revision 1.36 2001/03/30 07:05:49 morsch
38Final print-out in finish run.
39Write parton system for jet-production (preliminary solution).
40
9be6226f 41Revision 1.35 2001/03/09 13:03:40 morsch
42Process_t and Struc_Func_t moved to AliPythia.h
43
f1a48a38 44Revision 1.34 2001/02/14 15:50:40 hristov
45The last particle in event marked using SetHighWaterMark
46
28adac24 47Revision 1.33 2001/01/30 09:23:12 hristov
48Streamers removed (R.Brun)
49
a8a6107b 50Revision 1.32 2001/01/26 19:55:51 hristov
51Major upgrade of AliRoot code
52
2ab0c725 53Revision 1.31 2001/01/17 10:54:31 hristov
54Better protection against FPE
55
0f9e52a3 56Revision 1.30 2000/12/18 08:55:35 morsch
57Make AliPythia dependent generartors work with new scheme of random number generation
58
3356c022 59Revision 1.29 2000/12/04 11:22:03 morsch
60Init of sRandom as in 1.15
61
097ec7c1 62Revision 1.28 2000/12/02 11:41:39 morsch
63Use SetRandom() to initialize random number generator in constructor.
64
1ed25e11 65Revision 1.27 2000/11/30 20:29:02 morsch
66Initialise static variable sRandom in constructor: sRandom = fRandom;
67
ad1df918 68Revision 1.26 2000/11/30 07:12:50 alibrary
69Introducing new Rndm and QA classes
70
65fb704d 71Revision 1.25 2000/10/18 19:11:27 hristov
72Division by zero fixed
73
919c2096 74Revision 1.24 2000/09/18 10:41:35 morsch
75Add possibility to use nuclear structure functions from PDF library V8.
76
811826d8 77Revision 1.23 2000/09/14 14:05:40 morsch
78dito
79
819f84ad 80Revision 1.22 2000/09/14 14:02:22 morsch
81- Correct conversion from mm to cm when passing particle vertex to MC.
82- Correct handling of fForceDecay == all.
83
4451ef92 84Revision 1.21 2000/09/12 14:14:55 morsch
85Call fDecayer->ForceDecay() at the beginning of Generate().
86
00d6ce7d 87Revision 1.20 2000/09/06 14:29:33 morsch
88Use AliPythia for event generation an AliDecayPythia for decays.
89Correct handling of "nodecay" option
90
18edb254 91Revision 1.19 2000/07/11 18:24:56 fca
92Coding convention corrections + few minor bug fixes
93
aee8290b 94Revision 1.18 2000/06/30 12:40:34 morsch
95Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
96Geant units (cm).
97
1512e357 98Revision 1.17 2000/06/09 20:34:07 morsch
99All coding rule violations except RS3 corrected
100
f87cfe57 101Revision 1.16 2000/05/15 15:04:20 morsch
102The full event is written for fNtrack = -1
103Coding rule violations corrected.
104
5ddeb374 105Revision 1.15 2000/04/26 10:14:24 morsch
106Particles array has one entry more than pythia particle list. Upper bound of
107particle loop changed to np-1 (R. Guernane, AM)
108
2a0e6f5b 109Revision 1.14 2000/04/05 08:36:13 morsch
110Check status code of particles in Pythia event
111to avoid double counting as partonic state and final state particle.
112
23211d3e 113Revision 1.13 1999/11/09 07:38:48 fca
114Changes for compatibility with version 2.23 of ROOT
115
084c1b4a 116Revision 1.12 1999/11/03 17:43:20 fca
117New version from G.Martinez & A.Morsch
118
886b6f73 119Revision 1.11 1999/09/29 09:24:14 fca
120Introduction of the Copyright and cvs Log
4c039060 121*/
122
fe4da5cc 123#include "AliGenPythia.h"
fff02fee 124#include "AliGenPythiaEventHeader.h"
18edb254 125#include "AliDecayerPythia.h"
fe4da5cc 126#include "AliRun.h"
127#include "AliPythia.h"
18edb254 128#include "AliPDG.h"
1578254f 129#include <TParticle.h>
18edb254 130#include <TSystem.h>
f87cfe57 131
fe4da5cc 132 ClassImp(AliGenPythia)
133
134AliGenPythia::AliGenPythia()
fff02fee 135 :AliGenMC()
fe4da5cc 136{
18edb254 137// Default Constructor
138 fDecayer = new AliDecayerPythia();
11dfaf8d 139 SetEventListRange();
fe4da5cc 140}
141
142AliGenPythia::AliGenPythia(Int_t npart)
fff02fee 143 :AliGenMC(npart)
fe4da5cc 144{
145// default charm production at 5. 5 TeV
146// semimuonic decay
147// structure function GRVHO
148//
149 fXsection = 0.;
811826d8 150 fNucA1=0;
151 fNucA2=0;
fe4da5cc 152 SetProcess();
153 SetStrucFunc();
886b6f73 154 SetForceDecay();
fe4da5cc 155 SetPtHard();
156 SetEnergyCMS();
18edb254 157 fDecayer = new AliDecayerPythia();
1ed25e11 158 // Set random number generator
097ec7c1 159 sRandom=fRandom;
11dfaf8d 160 SetEventListRange();
fff02fee 161 fFlavorSelect = 0;
162 // Produced particles
163 fParticles = new TClonesArray("TParticle",1000);
c5d36e84 164 fEventVertex.Set(3);
fe4da5cc 165}
166
f87cfe57 167AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
168{
169// copy constructor
170}
171
fe4da5cc 172AliGenPythia::~AliGenPythia()
173{
5ddeb374 174// Destructor
fe4da5cc 175}
176
11dfaf8d 177void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
178{
179 // Set a range of event numbers, for which a table
180 // of generated particle will be printed
181 fDebugEventFirst = eventFirst;
182 fDebugEventLast = eventLast;
183 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
184}
185
fe4da5cc 186void AliGenPythia::Init()
187{
5ddeb374 188// Initialisation
18edb254 189 SetMC(AliPythia::Instance());
fe4da5cc 190 fPythia=(AliPythia*) fgMCEvGen;
191//
192 fParentWeight=1./Float_t(fNpart);
193//
194// Forward Paramters to the AliPythia object
00d6ce7d 195 fDecayer->SetForceDecay(fForceDecay);
18edb254 196 fDecayer->Init();
00d6ce7d 197
18edb254 198
fe4da5cc 199 fPythia->SetCKIN(3,fPtHardMin);
200 fPythia->SetCKIN(4,fPtHardMax);
811826d8 201 if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
fe4da5cc 202 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
18edb254 203
204 // fPythia->Pylist(0);
205 // fPythia->Pystat(2);
fe4da5cc 206// Parent and Children Selection
207 switch (fProcess)
208 {
f1a48a38 209 case kPyCharm:
fff02fee 210 fParentSelect[0] = 411;
211 fParentSelect[1] = 421;
212 fParentSelect[2] = 431;
213 fParentSelect[3] = 4122;
214 fFlavorSelect = 4;
fe4da5cc 215 break;
f1a48a38 216 case kPyCharmUnforced:
fff02fee 217 fParentSelect[0] = 411;
218 fParentSelect[1] = 421;
219 fParentSelect[2] = 431;
220 fParentSelect[3]= 4122;
221 fFlavorSelect = 4;
fe4da5cc 222 break;
f1a48a38 223 case kPyBeauty:
fff02fee 224 fParentSelect[0]= 511;
225 fParentSelect[1]= 521;
226 fParentSelect[2]= 531;
227 fParentSelect[3]= 5122;
228 fParentSelect[4]= 5132;
229 fParentSelect[5]= 5232;
230 fParentSelect[6]= 5332;
231 fFlavorSelect = 5;
fe4da5cc 232 break;
f1a48a38 233 case kPyBeautyUnforced:
fff02fee 234 fParentSelect[0] = 511;
235 fParentSelect[1] = 521;
236 fParentSelect[2] = 531;
237 fParentSelect[3] = 5122;
238 fParentSelect[4] = 5132;
239 fParentSelect[5] = 5232;
240 fParentSelect[6] = 5332;
241 fFlavorSelect = 5;
fe4da5cc 242 break;
f1a48a38 243 case kPyJpsiChi:
244 case kPyJpsi:
fff02fee 245 fParentSelect[0] = 443;
fe4da5cc 246 break;
f1a48a38 247 case kPyMb:
248 case kPyJets:
249 case kPyDirectGamma:
5be1fe76 250 break;
fe4da5cc 251 }
fff02fee 252 AliGenMC::Init();
fe4da5cc 253}
254
255void AliGenPythia::Generate()
256{
5ddeb374 257// Generate one event
00d6ce7d 258 fDecayer->ForceDecay();
a0b2b296 259
9be6226f 260 Float_t polar[3] = {0,0,0};
261 Float_t origin[3] = {0,0,0};
fff02fee 262 Float_t p[3];
5be1fe76 263// converts from mm/c to s
264 const Float_t kconv=0.001/2.999792458e8;
fe4da5cc 265//
fe4da5cc 266 Int_t nt=0;
267 Int_t jev=0;
7921755b 268 Int_t j, kf;
fe4da5cc 269 fTrials=0;
fff02fee 270
271// Set collision vertex position
aee8290b 272 if(fVertexSmear==kPerEvent) {
1512e357 273 fPythia->SetMSTP(151,1);
fe4da5cc 274 for (j=0;j<3;j++) {
fff02fee 275 fPythia->SetPARP(151+j, fOsigma[j]*10.);
fe4da5cc 276 }
aee8290b 277 } else if (fVertexSmear==kPerTrack) {
fe4da5cc 278 fPythia->SetMSTP(151,0);
fe4da5cc 279 }
fff02fee 280// event loop
fe4da5cc 281 while(1)
282 {
084c1b4a 283 fPythia->Pyevnt();
11dfaf8d 284 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
285 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
fe4da5cc 286 fTrials++;
fff02fee 287 fPythia->ImportParticles(fParticles,"All");
288
289//
290//
291//
292 Int_t i;
293
294 Int_t np = fParticles->GetEntriesFast();
a0b2b296 295 if (np == 0 ) continue;
296// Get event vertex and discard the event if the Z coord. is too big
4c07486d 297 TParticle *iparticle = (TParticle *) fParticles->At(0);
298 Float_t distz = iparticle->Vz()/10.;
299 if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
a0b2b296 300//
4c07486d 301 fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
302 fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
303 fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
a0b2b296 304//
fff02fee 305 Int_t* pParent = new Int_t[np];
306 Int_t* pSelected = new Int_t[np];
307 Int_t* trackIt = new Int_t[np];
308 for (i=0; i< np-1; i++) {
309 pParent[i] = -1;
310 pSelected[i] = 0;
311 }
fe4da5cc 312 printf("\n **************************************************%d\n",np);
fff02fee 313 Int_t nc = 0;
f1a48a38 314 if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
fff02fee 315 for (i = 0; i<np-1; i++) {
c5d36e84 316 iparticle = (TParticle *) fParticles->At(i);
23211d3e 317 Int_t ks = iparticle->GetStatusCode();
a8228d85 318 kf = CheckPDGCode(iparticle->GetPdgCode());
fff02fee 319// No initial state partons
2a0e6f5b 320 if (ks==21) continue;
fe4da5cc 321//
fff02fee 322// Heavy Flavor Selection
5be1fe76 323//
fff02fee 324 // quark ?
325 kf = TMath::Abs(kf);
326 Int_t kfl = kf;
327 // meson ?
328 if (kfl > 10) kfl/=100;
329 // baryon
330 if (kfl > 10) kfl/=10;
331 if (kfl > 10) kfl/=10;
332
333 Int_t ipa = iparticle->GetFirstMother()-1;
334 Int_t kfMo = 0;
335
336 if (ipa > -1) {
337 TParticle * mother = (TParticle *) fParticles->At(ipa);
338 kfMo = TMath::Abs(mother->GetPdgCode());
339 }
340// printf("\n particle (all) %d %d %d", i, pSelected[i], kf);
341 if (kfl >= fFlavorSelect) {
5be1fe76 342//
fff02fee 343// Heavy flavor hadron or quark
fe4da5cc 344//
fff02fee 345// Kinematic seletion on final state heavy flavor mesons
346 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
347 {
348 nc = -1;
349 break;
350 }
351 pSelected[i] = 1;
352// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
353 } else {
354// Kinematic seletion on decay products
355 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
356 && !KinematicSelection(iparticle, 1))
357 {
358 nc = -1;
359 break;
360 }
fe4da5cc 361//
fff02fee 362// Decay products
363// Select if mother was selected and is not tracked
364
365 if (pSelected[ipa] &&
366 !trackIt[ipa] && // mother will be tracked ?
367 kfMo != 5 && // mother is b-quark, don't store fragments
368 kfMo != 4 && // mother is c-quark, don't store fragments
369 kf != 92) // don't store string
370 {
fe4da5cc 371//
fff02fee 372// Semi-stable or de-selected: diselect decay products:
fe4da5cc 373//
5be1fe76 374//
fff02fee 375 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > 10e-15)
376 {
377 Int_t ipF = iparticle->GetFirstDaughter();
378 Int_t ipL = iparticle->GetLastDaughter();
379 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
380 }
381// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
382 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
383 }
384 }
385 if (pSelected[i] == -1) pSelected[i] = 0;
386 if (!pSelected[i]) continue;
387 nc++;
388// Decision on tracking
389 trackIt[i] = 0;
390//
391// Track final state particle
392 if (ks == 1) trackIt[i] = 1;
393// Track semi-stable particles
394 if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15)) trackIt[i] = 1;
395// Track particles selected by process if undecayed.
396 if (fForceDecay == kNoDecay) {
397 if (ParentSelected(kf)) trackIt[i] = 1;
398 } else {
399 if (ParentSelected(kf)) trackIt[i] = 0;
400 }
401//
402//
403
404 } // particle selection loop
405 if (nc > -1) {
406 for (i = 0; i<np-1; i++) {
407 if (!pSelected[i]) continue;
408 TParticle * iparticle = (TParticle *) fParticles->At(i);
409 kf = CheckPDGCode(iparticle->GetPdgCode());
a0b2b296 410 p[0] = iparticle->Px();
411 p[1] = iparticle->Py();
412 p[2] = iparticle->Pz();
413 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
414 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
415 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
416 Float_t tof = kconv*iparticle->T();
417 Int_t ipa = iparticle->GetFirstMother()-1;
fff02fee 418 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a99cf51f 419 SetTrack(fTrackIt*trackIt[i] ,
fff02fee 420 iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
421 pParent[i] = nt;
a99cf51f 422 KeepTrack(nt);
fff02fee 423 } // SetTrack loop
424 }
425 } else {
426 nc = GenerateMB();
5be1fe76 427 } // mb ?
fff02fee 428
fe4da5cc 429 if (nc > 0) {
5be1fe76 430 jev+=nc;
5ddeb374 431 if (jev >= fNpart || fNpart == -1) {
fe4da5cc 432 fKineBias=Float_t(fNpart)/Float_t(fTrials);
23211d3e 433 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
fe4da5cc 434 break;
435 }
436 }
437 } // event loop
a99cf51f 438 SetHighWaterMark(nt);
fe4da5cc 439// adjust weight due to kinematic selection
440 AdjustWeights();
441// get cross-section
442 fXsection=fPythia->GetPARI(1);
443}
444
fff02fee 445Int_t AliGenPythia::GenerateMB()
fe4da5cc 446{
fff02fee 447 Int_t i, kf, nt, iparent;
448 Int_t nc = 0;
449 Float_t p[3];
450 Float_t polar[3] = {0,0,0};
451 Float_t origin[3] = {0,0,0};
fff02fee 452// converts from mm/c to s
453 const Float_t kconv=0.001/2.999792458e8;
4451ef92 454
fff02fee 455 Int_t np = fParticles->GetEntriesFast();
456 Int_t* pParent = new Int_t[np];
457 for (i=0; i< np-1; i++) pParent[i] = -1;
458
459 for (i = 0; i<np-1; i++) {
460 Int_t trackIt = 0;
461 TParticle * iparticle = (TParticle *) fParticles->At(i);
462 kf = CheckPDGCode(iparticle->GetPdgCode());
463 Int_t ks = iparticle->GetStatusCode();
464 Int_t km = iparticle->GetFirstMother();
465// printf("\n Particle: %d %d %d", i, kf, ks);
466
467 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
468 (ks != 1) ||
469 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
470 nc++;
471 if (ks == 1) trackIt = 1;
472 Int_t ipa = iparticle->GetFirstMother()-1;
fe4da5cc 473
fff02fee 474 iparent = (ipa > -1) ? pParent[ipa] : -1;
fe4da5cc 475
476//
fff02fee 477// store track information
a0b2b296 478 p[0] = iparticle->Px();
479 p[1] = iparticle->Py();
480 p[2] = iparticle->Pz();
481 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
482 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
483 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
fff02fee 484 Float_t tof=kconv*iparticle->T();
a99cf51f 485 SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
a0b2b296 486 tof, kPPrimary, nt);
a99cf51f 487 KeepTrack(nt);
fff02fee 488 pParent[i] = nt;
489 } // select particle
490 } // particle loop
491
492 if (pParent) delete[] pParent;
fe4da5cc 493
fff02fee 494 printf("\n I've put %i particles on the stack \n",nc);
495 MakeHeader();
496 return nc;
497}
fe4da5cc 498
fe4da5cc 499
fff02fee 500void AliGenPythia::FinishRun()
501{
502// Print x-section summary
503 fPythia->Pystat(1);
fe4da5cc 504}
fff02fee 505
fe4da5cc 506void AliGenPythia::AdjustWeights()
507{
5ddeb374 508// Adjust the weights after generation of all events
509//
5ddeb374 510 TParticle *part;
fe4da5cc 511 Int_t ntrack=gAlice->GetNtrack();
512 for (Int_t i=0; i<ntrack; i++) {
2ab0c725 513 part= gAlice->Particle(i);
5ddeb374 514 part->SetWeight(part->GetWeight()*fKineBias);
fe4da5cc 515 }
516}
811826d8 517
518void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
519{
520// Treat protons as inside nuclei with mass numbers a1 and a2
521 fNucA1 = a1;
522 fNucA2 = a2;
523}
fff02fee 524
525
526void AliGenPythia::MakeHeader()
527{
528// Builds the event header, to be called after each event
529 AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
c5d36e84 530 ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
531 header->SetPrimaryVertex(fEventVertex);
532 gAlice->SetGenEventHeader(header);
fff02fee 533}
811826d8 534
535
f87cfe57 536AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
537{
538// Assignment operator
539 return *this;
540}
fe4da5cc 541
542
811826d8 543
a8a6107b 544#ifdef never
18edb254 545void AliGenPythia::Streamer(TBuffer &R__b)
546{
547 // Stream an object of class AliGenPythia.
548
549 if (R__b.IsReading()) {
550 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
551 AliGenerator::Streamer(R__b);
552 R__b >> (Int_t&)fProcess;
553 R__b >> (Int_t&)fStrucFunc;
554 R__b >> (Int_t&)fForceDecay;
555 R__b >> fEnergyCMS;
556 R__b >> fKineBias;
557 R__b >> fTrials;
558 fParentSelect.Streamer(R__b);
559 fChildSelect.Streamer(R__b);
560 R__b >> fXsection;
561// (AliPythia::Instance())->Streamer(R__b);
562 R__b >> fPtHardMin;
563 R__b >> fPtHardMax;
564// if (fDecayer) fDecayer->Streamer(R__b);
565 } else {
566 R__b.WriteVersion(AliGenPythia::IsA());
567 AliGenerator::Streamer(R__b);
568 R__b << (Int_t)fProcess;
569 R__b << (Int_t)fStrucFunc;
570 R__b << (Int_t)fForceDecay;
571 R__b << fEnergyCMS;
572 R__b << fKineBias;
573 R__b << fTrials;
574 fParentSelect.Streamer(R__b);
575 fChildSelect.Streamer(R__b);
576 R__b << fXsection;
577// R__b << fPythia;
578 R__b << fPtHardMin;
579 R__b << fPtHardMax;
580 // fDecayer->Streamer(R__b);
581 }
582}
a8a6107b 583#endif
18edb254 584