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