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