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