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