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