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