Changes for kPyJets
[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$
cc363573 18Revision 1.63 2002/10/14 14:55:35 hristov
19Merging the VirtualMC branch to the main development branch (HEAD)
20
b9d0a01d 21Revision 1.52.4.4 2002/10/10 16:40:08 hristov
22Updating VirtualMC to v3-09-02
23
24Revision 1.62 2002/09/24 10:00:01 morsch
25CheckTrigger() corrected.
26
6ae1960b 27Revision 1.61 2002/07/30 09:52:38 morsch
28Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
29
a7f29b20 30Revision 1.60 2002/07/19 14:49:03 morsch
31Typo corrected.
32
438003bb 33Revision 1.59 2002/07/19 14:35:36 morsch
34Count total number of trials. Print mean Q, x1, x2.
35
8a8f459f 36Revision 1.58 2002/07/17 10:04:09 morsch
37SetYHard method added.
38
35967770 39Revision 1.57 2002/05/22 13:22:53 morsch
40Process kPyMbNonDiffr added.
41
5f42cc16 42Revision 1.56 2002/04/26 10:30:01 morsch
43Option kPyBeautyPbMNR added. (N. Carrer).
44
c0acde67 45Revision 1.55 2002/04/17 10:23:56 morsch
46Coding Rule violations corrected.
47
6d4dd317 48Revision 1.54 2002/03/28 11:49:10 morsch
49Pass status code in SetTrack.
50
5ec84200 51Revision 1.53 2002/03/25 14:51:13 morsch
52New stack-fill and count options introduced (N. Carrer).
53
0d7c70a0 54Revision 1.51 2002/03/06 08:46:57 morsch
55- Loop until np-1
56- delete dyn. alloc. arrays (N. Carrer)
57
b9817493 58Revision 1.50 2002/03/03 13:48:50 morsch
59Option kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
60NLO calculations (Nicola Carrer).
61
fbd1348b 62Revision 1.49 2002/02/08 16:50:50 morsch
63Add name and title in constructor.
64
8b31bfac 65Revision 1.48 2001/12/20 11:44:28 morsch
66Add kinematic bias for direct gamma production.
67
b2bb450d 68Revision 1.47 2001/12/19 14:45:00 morsch
69Store number of trials in header.
70
fd8fb861 71Revision 1.46 2001/12/19 10:36:19 morsch
72Add possibility if jet kinematic biasing.
73
988b118b 74Revision 1.45 2001/11/28 08:06:52 morsch
75Use fMaxLifeTime parameter.
76
8cdfd303 77Revision 1.44 2001/11/27 13:13:07 morsch
78Maximum lifetime for long-lived particles to be put on the stack is parameter.
79It can be set via SetMaximumLifetime(..).
80
47fc6bd5 81Revision 1.43 2001/10/21 18:35:56 hristov
82Several pointers were set to zero in the default constructors to avoid memory management problems
83
2685bf00 84Revision 1.42 2001/10/15 08:21:55 morsch
85Vertex truncation settings moved to AliGenMC.
86
4c07486d 87Revision 1.41 2001/10/08 08:45:42 morsch
88Possibility of vertex cut added.
89
a0b2b296 90Revision 1.40 2001/09/25 11:30:23 morsch
91Pass event vertex to header.
92
c5d36e84 93Revision 1.39 2001/07/27 17:09:36 morsch
94Use local SetTrack, KeepTrack and SetHighWaterMark methods
95to delegate either to local stack or to stack owned by AliRun.
96(Piotr Skowronski, A.M.)
97
a99cf51f 98Revision 1.38 2001/07/13 10:58:54 morsch
99- Some coded moved to AliGenMC
100- Improved handling of secondary vertices.
101
fff02fee 102Revision 1.37 2001/06/28 11:17:28 morsch
103SetEventListRange setter added. Events in specified range are listed for
104debugging. (Yuri Kharlov)
105
11dfaf8d 106Revision 1.36 2001/03/30 07:05:49 morsch
107Final print-out in finish run.
108Write parton system for jet-production (preliminary solution).
109
9be6226f 110Revision 1.35 2001/03/09 13:03:40 morsch
111Process_t and Struc_Func_t moved to AliPythia.h
112
f1a48a38 113Revision 1.34 2001/02/14 15:50:40 hristov
114The last particle in event marked using SetHighWaterMark
115
28adac24 116Revision 1.33 2001/01/30 09:23:12 hristov
117Streamers removed (R.Brun)
118
a8a6107b 119Revision 1.32 2001/01/26 19:55:51 hristov
120Major upgrade of AliRoot code
121
2ab0c725 122Revision 1.31 2001/01/17 10:54:31 hristov
123Better protection against FPE
124
0f9e52a3 125Revision 1.30 2000/12/18 08:55:35 morsch
126Make AliPythia dependent generartors work with new scheme of random number generation
127
3356c022 128Revision 1.29 2000/12/04 11:22:03 morsch
129Init of sRandom as in 1.15
130
097ec7c1 131Revision 1.28 2000/12/02 11:41:39 morsch
132Use SetRandom() to initialize random number generator in constructor.
133
1ed25e11 134Revision 1.27 2000/11/30 20:29:02 morsch
135Initialise static variable sRandom in constructor: sRandom = fRandom;
136
ad1df918 137Revision 1.26 2000/11/30 07:12:50 alibrary
138Introducing new Rndm and QA classes
139
65fb704d 140Revision 1.25 2000/10/18 19:11:27 hristov
141Division by zero fixed
142
919c2096 143Revision 1.24 2000/09/18 10:41:35 morsch
144Add possibility to use nuclear structure functions from PDF library V8.
145
811826d8 146Revision 1.23 2000/09/14 14:05:40 morsch
147dito
148
819f84ad 149Revision 1.22 2000/09/14 14:02:22 morsch
150- Correct conversion from mm to cm when passing particle vertex to MC.
151- Correct handling of fForceDecay == all.
152
4451ef92 153Revision 1.21 2000/09/12 14:14:55 morsch
154Call fDecayer->ForceDecay() at the beginning of Generate().
155
00d6ce7d 156Revision 1.20 2000/09/06 14:29:33 morsch
157Use AliPythia for event generation an AliDecayPythia for decays.
158Correct handling of "nodecay" option
159
18edb254 160Revision 1.19 2000/07/11 18:24:56 fca
161Coding convention corrections + few minor bug fixes
162
aee8290b 163Revision 1.18 2000/06/30 12:40:34 morsch
164Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
165Geant units (cm).
166
1512e357 167Revision 1.17 2000/06/09 20:34:07 morsch
168All coding rule violations except RS3 corrected
169
f87cfe57 170Revision 1.16 2000/05/15 15:04:20 morsch
171The full event is written for fNtrack = -1
172Coding rule violations corrected.
173
5ddeb374 174Revision 1.15 2000/04/26 10:14:24 morsch
175Particles array has one entry more than pythia particle list. Upper bound of
176particle loop changed to np-1 (R. Guernane, AM)
177
2a0e6f5b 178Revision 1.14 2000/04/05 08:36:13 morsch
179Check status code of particles in Pythia event
180to avoid double counting as partonic state and final state particle.
181
23211d3e 182Revision 1.13 1999/11/09 07:38:48 fca
183Changes for compatibility with version 2.23 of ROOT
184
084c1b4a 185Revision 1.12 1999/11/03 17:43:20 fca
186New version from G.Martinez & A.Morsch
187
886b6f73 188Revision 1.11 1999/09/29 09:24:14 fca
189Introduction of the Copyright and cvs Log
4c039060 190*/
191
6d4dd317 192//
193// Generator using the TPythia interface (via AliPythia)
194// to generate pp collisions.
195// Using SetNuclei() also nuclear modifications to the structure functions
196// can be taken into account. This makes, of course, only sense for the
197// generation of the products of hard processes (heavy flavor, jets ...)
198//
199// andreas.morsch@cern.ch
200//
201
202
fe4da5cc 203#include "AliGenPythia.h"
fff02fee 204#include "AliGenPythiaEventHeader.h"
18edb254 205#include "AliDecayerPythia.h"
fe4da5cc 206#include "AliRun.h"
207#include "AliPythia.h"
18edb254 208#include "AliPDG.h"
cc363573 209#include "AliConst.h"
1578254f 210#include <TParticle.h>
18edb254 211#include <TSystem.h>
cc363573 212#include <TTree.h>
213#include <TDatabasePDG.h>
f87cfe57 214
fe4da5cc 215 ClassImp(AliGenPythia)
216
217AliGenPythia::AliGenPythia()
fff02fee 218 :AliGenMC()
fe4da5cc 219{
18edb254 220// Default Constructor
2685bf00 221 fParticles = 0;
222 fPythia = 0;
988b118b 223 fDecayer = new AliDecayerPythia();
11dfaf8d 224 SetEventListRange();
988b118b 225 SetJetPhiRange();
226 SetJetEtaRange();
cc363573 227 SetJetEtRange();
a7f29b20 228 SetGammaPhiRange();
229 SetGammaEtaRange();
fe4da5cc 230}
231
232AliGenPythia::AliGenPythia(Int_t npart)
fff02fee 233 :AliGenMC(npart)
fe4da5cc 234{
235// default charm production at 5. 5 TeV
236// semimuonic decay
237// structure function GRVHO
238//
8b31bfac 239 fName = "Pythia";
240 fTitle= "Particle Generator using PYTHIA";
fe4da5cc 241 fXsection = 0.;
811826d8 242 fNucA1=0;
243 fNucA2=0;
fe4da5cc 244 SetProcess();
245 SetStrucFunc();
886b6f73 246 SetForceDecay();
fe4da5cc 247 SetPtHard();
8a8f459f 248 SetYHard();
fe4da5cc 249 SetEnergyCMS();
18edb254 250 fDecayer = new AliDecayerPythia();
1ed25e11 251 // Set random number generator
097ec7c1 252 sRandom=fRandom;
fff02fee 253 fFlavorSelect = 0;
254 // Produced particles
255 fParticles = new TClonesArray("TParticle",1000);
c5d36e84 256 fEventVertex.Set(3);
988b118b 257 SetEventListRange();
258 SetJetPhiRange();
259 SetJetEtaRange();
cc363573 260 SetJetEtRange();
a7f29b20 261 SetGammaPhiRange();
262 SetGammaEtaRange();
719b6b8a 263 // Options determining what to keep in the stack (Heavy flavour generation)
264 fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
265 fFeedDownOpt = kTRUE; // allow feed down from higher family
266 // Fragmentation on/off
267 fFragmentation = kTRUE;
268 // Default counting mode
269 fCountMode = kCountAll;
fe4da5cc 270}
271
f87cfe57 272AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
273{
274// copy constructor
6d4dd317 275 Pythia.Copy(*this);
f87cfe57 276}
277
fe4da5cc 278AliGenPythia::~AliGenPythia()
279{
5ddeb374 280// Destructor
fe4da5cc 281}
282
11dfaf8d 283void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
284{
285 // Set a range of event numbers, for which a table
286 // of generated particle will be printed
287 fDebugEventFirst = eventFirst;
288 fDebugEventLast = eventLast;
289 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
290}
291
fe4da5cc 292void AliGenPythia::Init()
293{
5ddeb374 294// Initialisation
cc363573 295
296 SetMC(AliPythia::Instance());
fe4da5cc 297 fPythia=(AliPythia*) fgMCEvGen;
298//
299 fParentWeight=1./Float_t(fNpart);
300//
301// Forward Paramters to the AliPythia object
00d6ce7d 302 fDecayer->SetForceDecay(fForceDecay);
18edb254 303 fDecayer->Init();
00d6ce7d 304
18edb254 305
fe4da5cc 306 fPythia->SetCKIN(3,fPtHardMin);
35967770 307 fPythia->SetCKIN(4,fPtHardMax);
308 fPythia->SetCKIN(7,fYHardMin);
309 fPythia->SetCKIN(8,fYHardMax);
310
811826d8 311 if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
719b6b8a 312 // Fragmentation?
313 if (fFragmentation) {
314 fPythia->SetMSTP(111,1);
315 } else {
316 fPythia->SetMSTP(111,0);
317 }
318
fe4da5cc 319 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
18edb254 320
cc363573 321// initial state radiation
322 fPythia->SetMSTP(61,fGinit);
323// final state radiation
324 fPythia->SetMSTP(71,fGfinal);
325//
326 fPythia->SetMSTP(91,1);
327 fPythia->SetMSTU(16,1);
328 fPythia->SetMSTJ(1,1);
329//
fe4da5cc 330// Parent and Children Selection
331 switch (fProcess)
332 {
f1a48a38 333 case kPyCharm:
f1a48a38 334 case kPyCharmUnforced:
fbd1348b 335 case kPyCharmPbMNR:
b2bb450d 336 fParentSelect[0] = 411;
337 fParentSelect[1] = 421;
338 fParentSelect[2] = 431;
339 fParentSelect[3] = 4122;
8b31bfac 340 fFlavorSelect = 4;
fe4da5cc 341 break;
719b6b8a 342 case kPyD0PbMNR:
343 fParentSelect[0] = 421;
5f42cc16 344 fFlavorSelect = 4;
719b6b8a 345 break;
f1a48a38 346 case kPyBeauty:
c0acde67 347 case kPyBeautyPbMNR:
fff02fee 348 fParentSelect[0]= 511;
349 fParentSelect[1]= 521;
350 fParentSelect[2]= 531;
351 fParentSelect[3]= 5122;
352 fParentSelect[4]= 5132;
353 fParentSelect[5]= 5232;
354 fParentSelect[6]= 5332;
355 fFlavorSelect = 5;
fe4da5cc 356 break;
f1a48a38 357 case kPyBeautyUnforced:
fff02fee 358 fParentSelect[0] = 511;
359 fParentSelect[1] = 521;
360 fParentSelect[2] = 531;
361 fParentSelect[3] = 5122;
362 fParentSelect[4] = 5132;
363 fParentSelect[5] = 5232;
364 fParentSelect[6] = 5332;
365 fFlavorSelect = 5;
fe4da5cc 366 break;
f1a48a38 367 case kPyJpsiChi:
368 case kPyJpsi:
fff02fee 369 fParentSelect[0] = 443;
fe4da5cc 370 break;
f1a48a38 371 case kPyMb:
5f42cc16 372 case kPyMbNonDiffr:
f1a48a38 373 case kPyJets:
374 case kPyDirectGamma:
5be1fe76 375 break;
fe4da5cc 376 }
8a8f459f 377//
378// This counts the total number of calls to Pyevnt() per run.
379 fTrialsRun = 0;
380 fQ = 0.;
381 fX1 = 0.;
382 fX2 = 0.;
438003bb 383 fNev = 0 ;
8a8f459f 384//
fff02fee 385 AliGenMC::Init();
fe4da5cc 386}
387
388void AliGenPythia::Generate()
389{
5ddeb374 390// Generate one event
00d6ce7d 391 fDecayer->ForceDecay();
a0b2b296 392
9be6226f 393 Float_t polar[3] = {0,0,0};
394 Float_t origin[3] = {0,0,0};
fff02fee 395 Float_t p[3];
5be1fe76 396// converts from mm/c to s
397 const Float_t kconv=0.001/2.999792458e8;
fe4da5cc 398//
fe4da5cc 399 Int_t nt=0;
400 Int_t jev=0;
7921755b 401 Int_t j, kf;
fe4da5cc 402 fTrials=0;
fff02fee 403
fd8fb861 404 // Set collision vertex position
aee8290b 405 if(fVertexSmear==kPerEvent) {
1512e357 406 fPythia->SetMSTP(151,1);
fe4da5cc 407 for (j=0;j<3;j++) {
fff02fee 408 fPythia->SetPARP(151+j, fOsigma[j]*10.);
fe4da5cc 409 }
aee8290b 410 } else if (fVertexSmear==kPerTrack) {
fe4da5cc 411 fPythia->SetMSTP(151,0);
fe4da5cc 412 }
fff02fee 413// event loop
fe4da5cc 414 while(1)
415 {
084c1b4a 416 fPythia->Pyevnt();
11dfaf8d 417 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
418 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
fe4da5cc 419 fTrials++;
47fc6bd5 420
fff02fee 421 fPythia->ImportParticles(fParticles,"All");
422
423//
424//
425//
426 Int_t i;
427
428 Int_t np = fParticles->GetEntriesFast();
a0b2b296 429 if (np == 0 ) continue;
430// Get event vertex and discard the event if the Z coord. is too big
4c07486d 431 TParticle *iparticle = (TParticle *) fParticles->At(0);
432 Float_t distz = iparticle->Vz()/10.;
433 if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
a0b2b296 434//
4c07486d 435 fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
436 fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
437 fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
a0b2b296 438//
fff02fee 439 Int_t* pParent = new Int_t[np];
440 Int_t* pSelected = new Int_t[np];
441 Int_t* trackIt = new Int_t[np];
b9817493 442 for (i=0; i< np; i++) {
fff02fee 443 pParent[i] = -1;
444 pSelected[i] = 0;
719b6b8a 445 trackIt[i] = 0;
fff02fee 446 }
8a8f459f 447
719b6b8a 448 Int_t nc = 0; // Total n. of selected particles
449 Int_t nParents = 0; // Selected parents
450 Int_t nTkbles = 0; // Trackable particles
a7f29b20 451 if (fProcess != kPyMb && fProcess != kPyJets &&
452 fProcess != kPyDirectGamma &&
5f42cc16 453 fProcess != kPyMbNonDiffr) {
fd8fb861 454
b9817493 455 for (i = 0; i<np; i++) {
c5d36e84 456 iparticle = (TParticle *) fParticles->At(i);
23211d3e 457 Int_t ks = iparticle->GetStatusCode();
a8228d85 458 kf = CheckPDGCode(iparticle->GetPdgCode());
fff02fee 459// No initial state partons
2a0e6f5b 460 if (ks==21) continue;
fe4da5cc 461//
fff02fee 462// Heavy Flavor Selection
5be1fe76 463//
fff02fee 464 // quark ?
465 kf = TMath::Abs(kf);
466 Int_t kfl = kf;
467 // meson ?
468 if (kfl > 10) kfl/=100;
469 // baryon
470 if (kfl > 10) kfl/=10;
471 if (kfl > 10) kfl/=10;
472
473 Int_t ipa = iparticle->GetFirstMother()-1;
474 Int_t kfMo = 0;
475
476 if (ipa > -1) {
477 TParticle * mother = (TParticle *) fParticles->At(ipa);
478 kfMo = TMath::Abs(mother->GetPdgCode());
479 }
719b6b8a 480 // What to keep in Stack?
481 Bool_t flavorOK = kFALSE;
482 Bool_t selectOK = kFALSE;
483 if (fFeedDownOpt) {
484 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
485 } else {
486 if (kfl > fFlavorSelect) {
487 nc = -1;
488 break;
489 }
490 if (kfl == fFlavorSelect) flavorOK = kTRUE;
491 }
492 switch (fStackFillOpt) {
493 case kFlavorSelection:
494 selectOK = kTRUE;
495 break;
496 case kParentSelection:
497 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
498 break;
499 }
500 if (flavorOK && selectOK) {
5be1fe76 501//
fff02fee 502// Heavy flavor hadron or quark
fe4da5cc 503//
fff02fee 504// Kinematic seletion on final state heavy flavor mesons
505 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
506 {
719b6b8a 507 continue;
fff02fee 508 }
509 pSelected[i] = 1;
719b6b8a 510 if (ParentSelected(kf)) ++nParents; // Update parent count
fff02fee 511// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
512 } else {
513// Kinematic seletion on decay products
514 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
515 && !KinematicSelection(iparticle, 1))
516 {
719b6b8a 517 continue;
fff02fee 518 }
fe4da5cc 519//
fff02fee 520// Decay products
521// Select if mother was selected and is not tracked
522
523 if (pSelected[ipa] &&
524 !trackIt[ipa] && // mother will be tracked ?
525 kfMo != 5 && // mother is b-quark, don't store fragments
526 kfMo != 4 && // mother is c-quark, don't store fragments
527 kf != 92) // don't store string
528 {
fe4da5cc 529//
fff02fee 530// Semi-stable or de-selected: diselect decay products:
fe4da5cc 531//
5be1fe76 532//
47fc6bd5 533 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
fff02fee 534 {
535 Int_t ipF = iparticle->GetFirstDaughter();
536 Int_t ipL = iparticle->GetLastDaughter();
537 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
538 }
539// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
540 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
541 }
542 }
543 if (pSelected[i] == -1) pSelected[i] = 0;
544 if (!pSelected[i]) continue;
719b6b8a 545 // Count quarks only if you did not include fragmentation
546 if (fFragmentation && kf <= 10) continue;
fff02fee 547 nc++;
548// Decision on tracking
549 trackIt[i] = 0;
550//
551// Track final state particle
552 if (ks == 1) trackIt[i] = 1;
553// Track semi-stable particles
8cdfd303 554 if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
fff02fee 555// Track particles selected by process if undecayed.
556 if (fForceDecay == kNoDecay) {
557 if (ParentSelected(kf)) trackIt[i] = 1;
558 } else {
559 if (ParentSelected(kf)) trackIt[i] = 0;
560 }
719b6b8a 561 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
fff02fee 562//
563//
564
565 } // particle selection loop
719b6b8a 566 if (nc > 0) {
b9817493 567 for (i = 0; i<np; i++) {
fff02fee 568 if (!pSelected[i]) continue;
569 TParticle * iparticle = (TParticle *) fParticles->At(i);
570 kf = CheckPDGCode(iparticle->GetPdgCode());
5ec84200 571 Int_t ks = iparticle->GetStatusCode();
a0b2b296 572 p[0] = iparticle->Px();
573 p[1] = iparticle->Py();
574 p[2] = iparticle->Pz();
575 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
576 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
577 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
578 Float_t tof = kconv*iparticle->T();
579 Int_t ipa = iparticle->GetFirstMother()-1;
fff02fee 580 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a99cf51f 581 SetTrack(fTrackIt*trackIt[i] ,
5ec84200 582 iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
fff02fee 583 pParent[i] = nt;
a99cf51f 584 KeepTrack(nt);
fff02fee 585 } // SetTrack loop
586 }
587 } else {
588 nc = GenerateMB();
5be1fe76 589 } // mb ?
fff02fee 590
b9817493 591 if (pParent) delete[] pParent;
592 if (pSelected) delete[] pSelected;
593 if (trackIt) delete[] trackIt;
594
fe4da5cc 595 if (nc > 0) {
719b6b8a 596 switch (fCountMode) {
597 case kCountAll:
598 // printf(" Count all \n");
599 jev += nc;
600 break;
601 case kCountParents:
602 // printf(" Count parents \n");
603 jev += nParents;
604 break;
605 case kCountTrackables:
606 // printf(" Count trackable \n");
607 jev += nTkbles;
608 break;
609 }
5ddeb374 610 if (jev >= fNpart || fNpart == -1) {
fe4da5cc 611 fKineBias=Float_t(fNpart)/Float_t(fTrials);
23211d3e 612 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
8a8f459f 613
614 fQ += fPythia->GetVINT(51);
615 fX1 += fPythia->GetVINT(41);
616 fX2 += fPythia->GetVINT(42);
617 fTrialsRun += fTrials;
618 fNev++;
fd8fb861 619 MakeHeader();
fe4da5cc 620 break;
621 }
622 }
623 } // event loop
a99cf51f 624 SetHighWaterMark(nt);
fe4da5cc 625// adjust weight due to kinematic selection
626 AdjustWeights();
627// get cross-section
628 fXsection=fPythia->GetPARI(1);
629}
630
fff02fee 631Int_t AliGenPythia::GenerateMB()
fe4da5cc 632{
5ec84200 633//
634// Min Bias selection and other global selections
635//
fff02fee 636 Int_t i, kf, nt, iparent;
637 Int_t nc = 0;
638 Float_t p[3];
639 Float_t polar[3] = {0,0,0};
640 Float_t origin[3] = {0,0,0};
fff02fee 641// converts from mm/c to s
642 const Float_t kconv=0.001/2.999792458e8;
4451ef92 643
fff02fee 644 Int_t np = fParticles->GetEntriesFast();
645 Int_t* pParent = new Int_t[np];
b9817493 646 for (i=0; i< np; i++) pParent[i] = -1;
b2bb450d 647 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
988b118b 648 TParticle* jet1 = (TParticle *) fParticles->At(6);
649 TParticle* jet2 = (TParticle *) fParticles->At(7);
650 if (!CheckTrigger(jet1, jet2)) return 0;
651 }
652
b9817493 653 for (i = 0; i<np; i++) {
fff02fee 654 Int_t trackIt = 0;
655 TParticle * iparticle = (TParticle *) fParticles->At(i);
656 kf = CheckPDGCode(iparticle->GetPdgCode());
657 Int_t ks = iparticle->GetStatusCode();
658 Int_t km = iparticle->GetFirstMother();
fff02fee 659 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
660 (ks != 1) ||
661 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
662 nc++;
663 if (ks == 1) trackIt = 1;
664 Int_t ipa = iparticle->GetFirstMother()-1;
35967770 665
fff02fee 666 iparent = (ipa > -1) ? pParent[ipa] : -1;
35967770 667
fe4da5cc 668//
fff02fee 669// store track information
a0b2b296 670 p[0] = iparticle->Px();
671 p[1] = iparticle->Py();
672 p[2] = iparticle->Pz();
673 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
674 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
675 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
fff02fee 676 Float_t tof=kconv*iparticle->T();
a99cf51f 677 SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
35967770 678 tof, kPPrimary, nt, 1., ks);
a99cf51f 679 KeepTrack(nt);
fff02fee 680 pParent[i] = nt;
681 } // select particle
682 } // particle loop
683
684 if (pParent) delete[] pParent;
fe4da5cc 685
fff02fee 686 printf("\n I've put %i particles on the stack \n",nc);
fff02fee 687 return nc;
688}
fe4da5cc 689
fe4da5cc 690
fff02fee 691void AliGenPythia::FinishRun()
692{
693// Print x-section summary
694 fPythia->Pystat(1);
8a8f459f 695 fQ /= fNev;
696 fX1 /= fNev;
697 fX2 /= fNev;
698 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
699 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
700
701
fe4da5cc 702}
fff02fee 703
fe4da5cc 704void AliGenPythia::AdjustWeights()
705{
5ddeb374 706// Adjust the weights after generation of all events
707//
5ddeb374 708 TParticle *part;
fe4da5cc 709 Int_t ntrack=gAlice->GetNtrack();
710 for (Int_t i=0; i<ntrack; i++) {
2ab0c725 711 part= gAlice->Particle(i);
5ddeb374 712 part->SetWeight(part->GetWeight()*fKineBias);
fe4da5cc 713 }
714}
811826d8 715
716void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
717{
718// Treat protons as inside nuclei with mass numbers a1 and a2
719 fNucA1 = a1;
720 fNucA2 = a2;
721}
fff02fee 722
723
cc363573 724void AliGenPythia::MakeHeader()
fff02fee 725{
726// Builds the event header, to be called after each event
727 AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
cc363573 728//
729// Event type
c5d36e84 730 ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
cc363573 731//
732// Number of trials
fd8fb861 733 ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
cc363573 734//
735// Event Vertex
c5d36e84 736 header->SetPrimaryVertex(fEventVertex);
cc363573 737//
738// Jets that have triggered
739 if (fProcess == kPyJets)
740 {
741 Int_t ntrig, njet;
742 Float_t jets[4][10];
743 GetJets(5.0, 1, njet, ntrig, jets);
744 for (Int_t i = 0; i < ntrig; i++) {
745 ((AliGenPythiaEventHeader*) header)->AddJet(jets[0][i], jets[1][i], jets[2][i],
746 jets[3][i]);
747 }
748 }
c5d36e84 749 gAlice->SetGenEventHeader(header);
fff02fee 750}
811826d8 751
988b118b 752
cc363573 753Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
988b118b 754{
755// Check the kinematic trigger condition
756//
b2bb450d 757 Double_t eta[2];
758 eta[0] = jet1->Eta();
759 eta[1] = jet2->Eta();
760 Double_t phi[2];
761 phi[0] = jet1->Phi();
762 phi[1] = jet2->Phi();
763 Int_t pdg[2];
764 pdg[0] = jet1->GetPdgCode();
765 pdg[1] = jet2->GetPdgCode();
988b118b 766 Bool_t triggered = kFALSE;
b2bb450d 767
768 if (fProcess == kPyJets) {
cc363573 769 Int_t njets = 0;
770 Int_t ntrig = 0;
771 Float_t jets[4][10];
772//
773// Use Pythia clustering on parton level to determine jet axis
774//
775 GetJets(5.0, 1, njets, ntrig, jets);
776 if (ntrig) triggered = kTRUE;
777//
b2bb450d 778 } else {
779 Int_t ij = 0;
780 Int_t ig = 1;
781 if (pdg[0] == kGamma) {
782 ij = 1;
783 ig = 0;
784 }
785 //Check eta range first...
786 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
787 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
788 {
789 //Eta is okay, now check phi range
790 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
791 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
792 {
793 triggered = kTRUE;
794 }
795 }
988b118b 796 }
797 return triggered;
798}
811826d8 799
f87cfe57 800AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
801{
802// Assignment operator
803 return *this;
804}
fe4da5cc 805
cc363573 806void AliGenPythia::LoadEvent()
807{
808//
809// Load event into Pythia Common Block
810//
811
812
813 Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries();
814 (fPythia->GetPyjets())->N = npart;
815
816 for (Int_t part = 0; part < npart; part++) {
817 TParticle *MPart = gAlice->Particle(part);
818 Int_t kf = MPart->GetPdgCode();
819 Int_t ks = MPart->GetStatusCode();
820 Float_t px = MPart->Px();
821 Float_t py = MPart->Py();
822 Float_t pz = MPart->Pz();
823 Float_t e = MPart->Energy();
824 Float_t p = TMath::Sqrt(px * px + py * py + pz * pz);
825 Float_t m = TMath::Sqrt(e * e - p * p);
826
827
828 (fPythia->GetPyjets())->P[0][part] = px;
829 (fPythia->GetPyjets())->P[1][part] = py;
830 (fPythia->GetPyjets())->P[2][part] = pz;
831 (fPythia->GetPyjets())->P[3][part] = e;
832 (fPythia->GetPyjets())->P[4][part] = m;
833
834 (fPythia->GetPyjets())->K[1][part] = kf;
835 (fPythia->GetPyjets())->K[0][part] = ks;
836 }
837}
838
839
840void AliGenPythia::GetJets(Float_t dist, Int_t part, Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
841{
842//
843// Calls the Pythia clustering algorithm to find jets in the current event
844//
845 Int_t n = fPythia->GetN();
846 nJets = 0;
847 nJetsTrig = 0;
848//
849// Configure cluster algorithm
850//
851 fPythia->SetPARU(43, dist);
852 fPythia->SetMSTU(41, part);
853//
854// Call cluster algorithm
855//
856 fPythia->Pyclus(nJets);
857//
858// Loading jets from common block
859//
860 Int_t i;
861 for (i = 0; i < nJets; i++) {
862 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
863 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
864 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
865 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
866 Float_t pt = TMath::Sqrt(px * px + py * py);
867 Float_t phi = TMath::ATan2(py,px);
868 Float_t theta = TMath::ATan2(pt,pz);
869 Float_t et = e * TMath::Sin(theta);
870 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
871
872 if (
873 eta > fEtaMinJet && eta < fEtaMaxJet &&
874 phi > fPhiMinJet && eta < fPhiMaxJet &&
875 et > fEtMinJet && et < fEtMaxJet
876 )
877 {
878// printf("\n*Trigger-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
879 jets[0][nJetsTrig] = px;
880 jets[1][nJetsTrig] = py;
881 jets[2][nJetsTrig] = pz;
882 jets[3][nJetsTrig] = e;
883 nJetsTrig++;
884
885 } else {
886// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
887 }
888 }
889}
fe4da5cc 890
811826d8 891
a8a6107b 892#ifdef never
18edb254 893void AliGenPythia::Streamer(TBuffer &R__b)
894{
895 // Stream an object of class AliGenPythia.
896
897 if (R__b.IsReading()) {
898 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
899 AliGenerator::Streamer(R__b);
900 R__b >> (Int_t&)fProcess;
901 R__b >> (Int_t&)fStrucFunc;
902 R__b >> (Int_t&)fForceDecay;
903 R__b >> fEnergyCMS;
904 R__b >> fKineBias;
905 R__b >> fTrials;
906 fParentSelect.Streamer(R__b);
907 fChildSelect.Streamer(R__b);
908 R__b >> fXsection;
909// (AliPythia::Instance())->Streamer(R__b);
910 R__b >> fPtHardMin;
911 R__b >> fPtHardMax;
912// if (fDecayer) fDecayer->Streamer(R__b);
913 } else {
914 R__b.WriteVersion(AliGenPythia::IsA());
915 AliGenerator::Streamer(R__b);
916 R__b << (Int_t)fProcess;
917 R__b << (Int_t)fStrucFunc;
918 R__b << (Int_t)fForceDecay;
919 R__b << fEnergyCMS;
920 R__b << fKineBias;
921 R__b << fTrials;
922 fParentSelect.Streamer(R__b);
923 fChildSelect.Streamer(R__b);
924 R__b << fXsection;
925// R__b << fPythia;
926 R__b << fPtHardMin;
927 R__b << fPtHardMax;
928 // fDecayer->Streamer(R__b);
929 }
930}
a8a6107b 931#endif
18edb254 932