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