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