1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.38 2001/07/13 10:58:54 morsch
19 - Some coded moved to AliGenMC
20 - Improved handling of secondary vertices.
22 Revision 1.37 2001/06/28 11:17:28 morsch
23 SetEventListRange setter added. Events in specified range are listed for
24 debugging. (Yuri Kharlov)
26 Revision 1.36 2001/03/30 07:05:49 morsch
27 Final print-out in finish run.
28 Write parton system for jet-production (preliminary solution).
30 Revision 1.35 2001/03/09 13:03:40 morsch
31 Process_t and Struc_Func_t moved to AliPythia.h
33 Revision 1.34 2001/02/14 15:50:40 hristov
34 The last particle in event marked using SetHighWaterMark
36 Revision 1.33 2001/01/30 09:23:12 hristov
37 Streamers removed (R.Brun)
39 Revision 1.32 2001/01/26 19:55:51 hristov
40 Major upgrade of AliRoot code
42 Revision 1.31 2001/01/17 10:54:31 hristov
43 Better protection against FPE
45 Revision 1.30 2000/12/18 08:55:35 morsch
46 Make AliPythia dependent generartors work with new scheme of random number generation
48 Revision 1.29 2000/12/04 11:22:03 morsch
49 Init of sRandom as in 1.15
51 Revision 1.28 2000/12/02 11:41:39 morsch
52 Use SetRandom() to initialize random number generator in constructor.
54 Revision 1.27 2000/11/30 20:29:02 morsch
55 Initialise static variable sRandom in constructor: sRandom = fRandom;
57 Revision 1.26 2000/11/30 07:12:50 alibrary
58 Introducing new Rndm and QA classes
60 Revision 1.25 2000/10/18 19:11:27 hristov
61 Division by zero fixed
63 Revision 1.24 2000/09/18 10:41:35 morsch
64 Add possibility to use nuclear structure functions from PDF library V8.
66 Revision 1.23 2000/09/14 14:05:40 morsch
69 Revision 1.22 2000/09/14 14:02:22 morsch
70 - Correct conversion from mm to cm when passing particle vertex to MC.
71 - Correct handling of fForceDecay == all.
73 Revision 1.21 2000/09/12 14:14:55 morsch
74 Call fDecayer->ForceDecay() at the beginning of Generate().
76 Revision 1.20 2000/09/06 14:29:33 morsch
77 Use AliPythia for event generation an AliDecayPythia for decays.
78 Correct handling of "nodecay" option
80 Revision 1.19 2000/07/11 18:24:56 fca
81 Coding convention corrections + few minor bug fixes
83 Revision 1.18 2000/06/30 12:40:34 morsch
84 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
87 Revision 1.17 2000/06/09 20:34:07 morsch
88 All coding rule violations except RS3 corrected
90 Revision 1.16 2000/05/15 15:04:20 morsch
91 The full event is written for fNtrack = -1
92 Coding rule violations corrected.
94 Revision 1.15 2000/04/26 10:14:24 morsch
95 Particles array has one entry more than pythia particle list. Upper bound of
96 particle loop changed to np-1 (R. Guernane, AM)
98 Revision 1.14 2000/04/05 08:36:13 morsch
99 Check status code of particles in Pythia event
100 to avoid double counting as partonic state and final state particle.
102 Revision 1.13 1999/11/09 07:38:48 fca
103 Changes for compatibility with version 2.23 of ROOT
105 Revision 1.12 1999/11/03 17:43:20 fca
106 New version from G.Martinez & A.Morsch
108 Revision 1.11 1999/09/29 09:24:14 fca
109 Introduction of the Copyright and cvs Log
112 #include "AliGenPythia.h"
113 #include "AliGenPythiaEventHeader.h"
114 #include "AliDecayerPythia.h"
116 #include "AliPythia.h"
118 #include <TParticle.h>
121 ClassImp(AliGenPythia)
123 AliGenPythia::AliGenPythia()
126 // Default Constructor
127 fDecayer = new AliDecayerPythia();
131 AliGenPythia::AliGenPythia(Int_t npart)
134 // default charm production at 5. 5 TeV
136 // structure function GRVHO
146 fDecayer = new AliDecayerPythia();
147 // Set random number generator
151 // Produced particles
152 fParticles = new TClonesArray("TParticle",1000);
155 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
160 AliGenPythia::~AliGenPythia()
165 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
167 // Set a range of event numbers, for which a table
168 // of generated particle will be printed
169 fDebugEventFirst = eventFirst;
170 fDebugEventLast = eventLast;
171 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
174 void AliGenPythia::Init()
177 SetMC(AliPythia::Instance());
178 fPythia=(AliPythia*) fgMCEvGen;
180 fParentWeight=1./Float_t(fNpart);
182 // Forward Paramters to the AliPythia object
183 fDecayer->SetForceDecay(fForceDecay);
187 fPythia->SetCKIN(3,fPtHardMin);
188 fPythia->SetCKIN(4,fPtHardMax);
189 if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
190 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
192 // fPythia->Pylist(0);
193 // fPythia->Pystat(2);
194 // Parent and Children Selection
198 fParentSelect[0] = 411;
199 fParentSelect[1] = 421;
200 fParentSelect[2] = 431;
201 fParentSelect[3] = 4122;
204 case kPyCharmUnforced:
205 fParentSelect[0] = 411;
206 fParentSelect[1] = 421;
207 fParentSelect[2] = 431;
208 fParentSelect[3]= 4122;
212 fParentSelect[0]= 511;
213 fParentSelect[1]= 521;
214 fParentSelect[2]= 531;
215 fParentSelect[3]= 5122;
216 fParentSelect[4]= 5132;
217 fParentSelect[5]= 5232;
218 fParentSelect[6]= 5332;
221 case kPyBeautyUnforced:
222 fParentSelect[0] = 511;
223 fParentSelect[1] = 521;
224 fParentSelect[2] = 531;
225 fParentSelect[3] = 5122;
226 fParentSelect[4] = 5132;
227 fParentSelect[5] = 5232;
228 fParentSelect[6] = 5332;
233 fParentSelect[0] = 443;
243 void AliGenPythia::Generate()
245 // Generate one event
246 fDecayer->ForceDecay();
248 Float_t polar[3] = {0,0,0};
249 Float_t origin[3] = {0,0,0};
250 Float_t origin0[3] = {0,0,0};
252 // converts from mm/c to s
253 const Float_t kconv=0.001/2.999792458e8;
260 // Set collision vertex position
261 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
262 if(fVertexSmear==kPerEvent) {
263 fPythia->SetMSTP(151,1);
265 fPythia->SetPARP(151+j, fOsigma[j]*10.);
267 } else if (fVertexSmear==kPerTrack) {
268 fPythia->SetMSTP(151,0);
274 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
275 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
277 fPythia->ImportParticles(fParticles,"All");
284 Int_t np = fParticles->GetEntriesFast();
285 Int_t* pParent = new Int_t[np];
286 Int_t* pSelected = new Int_t[np];
287 Int_t* trackIt = new Int_t[np];
288 for (i=0; i< np-1; i++) {
292 printf("\n **************************************************%d\n",np);
294 if (np == 0 ) continue;
295 if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
296 for (i = 0; i<np-1; i++) {
297 TParticle * iparticle = (TParticle *) fParticles->At(i);
298 Int_t ks = iparticle->GetStatusCode();
299 kf = CheckPDGCode(iparticle->GetPdgCode());
300 // No initial state partons
301 if (ks==21) continue;
303 // Heavy Flavor Selection
309 if (kfl > 10) kfl/=100;
311 if (kfl > 10) kfl/=10;
312 if (kfl > 10) kfl/=10;
314 Int_t ipa = iparticle->GetFirstMother()-1;
318 TParticle * mother = (TParticle *) fParticles->At(ipa);
319 kfMo = TMath::Abs(mother->GetPdgCode());
321 // printf("\n particle (all) %d %d %d", i, pSelected[i], kf);
322 if (kfl >= fFlavorSelect) {
324 // Heavy flavor hadron or quark
326 // Kinematic seletion on final state heavy flavor mesons
327 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
333 // printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
335 // Kinematic seletion on decay products
336 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
337 && !KinematicSelection(iparticle, 1))
344 // Select if mother was selected and is not tracked
346 if (pSelected[ipa] &&
347 !trackIt[ipa] && // mother will be tracked ?
348 kfMo != 5 && // mother is b-quark, don't store fragments
349 kfMo != 4 && // mother is c-quark, don't store fragments
350 kf != 92) // don't store string
353 // Semi-stable or de-selected: diselect decay products:
356 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > 10e-15)
358 Int_t ipF = iparticle->GetFirstDaughter();
359 Int_t ipL = iparticle->GetLastDaughter();
360 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
362 // printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
363 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
366 if (pSelected[i] == -1) pSelected[i] = 0;
367 if (!pSelected[i]) continue;
369 // Decision on tracking
372 // Track final state particle
373 if (ks == 1) trackIt[i] = 1;
374 // Track semi-stable particles
375 if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15)) trackIt[i] = 1;
376 // Track particles selected by process if undecayed.
377 if (fForceDecay == kNoDecay) {
378 if (ParentSelected(kf)) trackIt[i] = 1;
380 if (ParentSelected(kf)) trackIt[i] = 0;
385 } // particle selection loop
387 for (i = 0; i<np-1; i++) {
388 if (!pSelected[i]) continue;
389 TParticle * iparticle = (TParticle *) fParticles->At(i);
390 kf = CheckPDGCode(iparticle->GetPdgCode());
391 p[0]=iparticle->Px();
392 p[1]=iparticle->Py();
393 p[2]=iparticle->Pz();
394 origin[0]=origin0[0]+iparticle->Vx()/10.;
395 origin[1]=origin0[1]+iparticle->Vy()/10.;
396 origin[2]=origin0[2]+iparticle->Vz()/10.;
397 Float_t tof=kconv*iparticle->T();
398 Int_t ipa = iparticle->GetFirstMother()-1;
399 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
400 // printf("\n SetTrack %d %d %d %d", i, iparent, kf, trackIt[i]);
401 SetTrack(fTrackIt*trackIt[i] ,
402 iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
413 if (jev >= fNpart || fNpart == -1) {
414 fKineBias=Float_t(fNpart)/Float_t(fTrials);
415 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
420 SetHighWaterMark(nt);
421 // adjust weight due to kinematic selection
424 fXsection=fPythia->GetPARI(1);
427 Int_t AliGenPythia::GenerateMB()
429 Int_t i, kf, nt, iparent;
432 Float_t polar[3] = {0,0,0};
433 Float_t origin[3] = {0,0,0};
434 Float_t origin0[3] = {0,0,0};
435 // converts from mm/c to s
436 const Float_t kconv=0.001/2.999792458e8;
438 Int_t np = fParticles->GetEntriesFast();
439 Int_t* pParent = new Int_t[np];
440 for (i=0; i< np-1; i++) pParent[i] = -1;
442 for (i = 0; i<np-1; i++) {
444 TParticle * iparticle = (TParticle *) fParticles->At(i);
445 kf = CheckPDGCode(iparticle->GetPdgCode());
446 Int_t ks = iparticle->GetStatusCode();
447 Int_t km = iparticle->GetFirstMother();
448 // printf("\n Particle: %d %d %d", i, kf, ks);
450 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
452 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
454 if (ks == 1) trackIt = 1;
455 Int_t ipa = iparticle->GetFirstMother()-1;
457 iparent = (ipa > -1) ? pParent[ipa] : -1;
460 // store track information
461 p[0]=iparticle->Px();
462 p[1]=iparticle->Py();
463 p[2]=iparticle->Pz();
464 origin[0]=origin0[0]+iparticle->Vx()/10.;
465 origin[1]=origin0[1]+iparticle->Vy()/10.;
466 origin[2]=origin0[2]+iparticle->Vz()/10.;
467 Float_t tof=kconv*iparticle->T();
468 SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
475 if (pParent) delete[] pParent;
477 printf("\n I've put %i particles on the stack \n",nc);
483 void AliGenPythia::FinishRun()
485 // Print x-section summary
489 void AliGenPythia::AdjustWeights()
491 // Adjust the weights after generation of all events
494 Int_t ntrack=gAlice->GetNtrack();
495 for (Int_t i=0; i<ntrack; i++) {
496 part= gAlice->Particle(i);
497 part->SetWeight(part->GetWeight()*fKineBias);
501 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
503 // Treat protons as inside nuclei with mass numbers a1 and a2
509 void AliGenPythia::MakeHeader()
511 // Builds the event header, to be called after each event
512 AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
513 ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
514 gAlice->SetGenEventHeader(header);
518 AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
520 // Assignment operator
527 void AliGenPythia::Streamer(TBuffer &R__b)
529 // Stream an object of class AliGenPythia.
531 if (R__b.IsReading()) {
532 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
533 AliGenerator::Streamer(R__b);
534 R__b >> (Int_t&)fProcess;
535 R__b >> (Int_t&)fStrucFunc;
536 R__b >> (Int_t&)fForceDecay;
540 fParentSelect.Streamer(R__b);
541 fChildSelect.Streamer(R__b);
543 // (AliPythia::Instance())->Streamer(R__b);
546 // if (fDecayer) fDecayer->Streamer(R__b);
548 R__b.WriteVersion(AliGenPythia::IsA());
549 AliGenerator::Streamer(R__b);
550 R__b << (Int_t)fProcess;
551 R__b << (Int_t)fStrucFunc;
552 R__b << (Int_t)fForceDecay;
556 fParentSelect.Streamer(R__b);
557 fChildSelect.Streamer(R__b);
562 // fDecayer->Streamer(R__b);