PYTHIA6/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 /* $Id$ */
  17
  18 //
  19 // Generator using the TPythia interface (via AliPythia)
  20 // to generate pp collisions.
  21 // Using SetNuclei() also nuclear modifications to the structure functions
  22 // can be taken into account. This makes, of course, only sense for the
  23 // generation of the products of hard processes (heavy flavor, jets ...)
  24 //
  25 // andreas.morsch@cern.ch
  26 //
  27
  28 #include <TDatabasePDG.h>
  29 #include <TParticle.h>
  30 #include <TPDGCode.h>
  31 #include <TSystem.h>
  32 #include <TTree.h>
  33 #include "AliConst.h"
  34 #include "AliDecayerPythia.h"
  35 #include "AliGenPythia.h"
  36 #include "AliHeader.h"
  37 #include "AliGenPythiaEventHeader.h"
  38 #include "AliPythia.h"
  39 #include "AliPythiaRndm.h"
  40 #include "AliRun.h"
  41 #include "AliStack.h"
  42 #include "AliRunLoader.h"
  43 #include "AliMC.h"
  44 #include "pyquenCommon.h"
  45
  46 ClassImp(AliGenPythia)
  47
  48 AliGenPythia::AliGenPythia()
  49                  :AliGenMC()
  50 {
  51 // Default Constructor
  52   fParticles = 0;
  53   fPythia    = 0;
  54   fHeader = 0;
  55   fReadFromFile = 0;
  56   fEventTime = 0.;
  57   fDecayer   = new AliDecayerPythia();
  58   SetEventListRange();
  59   SetJetPhiRange();
  60   SetJetEtaRange();
  61   SetJetEtRange();
  62   SetGammaPhiRange();
  63   SetGammaEtaRange();
  64   SetPtKick();
  65   SetQuench();
  66   SetHadronisation();
  67   fSetNuclei = kFALSE;
  68   if (!AliPythiaRndm::GetPythiaRandom())
  69     AliPythiaRndm::SetPythiaRandom(GetRandom());
  70 }
  71
  72 AliGenPythia::AliGenPythia(Int_t npart)
  73                  :AliGenMC(npart)
  74 {
  75 // default charm production at 5. 5 TeV
  76 // semimuonic decay
  77 // structure function GRVHO
  78 //
  79     fName = "Pythia";
  80     fTitle= "Particle Generator using PYTHIA";
  81     fXsection  = 0.;
  82     fReadFromFile = 0;
  83     fEventTime = 0.;
  84     SetProcess();
  85     SetStrucFunc();
  86     SetForceDecay();
  87     SetPtHard();
  88     SetYHard();
  89     SetEnergyCMS();
  90     fDecayer = new AliDecayerPythia();
  91     // Set random number generator
  92     if (!AliPythiaRndm::GetPythiaRandom())
  93       AliPythiaRndm::SetPythiaRandom(GetRandom());
  94     fFlavorSelect   = 0;
  95     // Produced particles
  96     fParticles = new TClonesArray("TParticle",1000);
  97     fHeader = 0;
  98     SetEventListRange();
  99     SetJetPhiRange();
 100     SetJetEtaRange();
 101     SetJetEtRange();
 102     SetGammaPhiRange();
 103     SetGammaEtaRange();
 104     SetJetReconstructionMode();
 105     SetQuench();
 106     SetHadronisation();
 107     SetPtKick();
 108     // Options determining what to keep in the stack (Heavy flavour generation)
 109     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
 110     fFeedDownOpt = kTRUE;             // allow feed down from higher family
 111     // Fragmentation on/off
 112     fFragmentation = kTRUE;
 113     // Default counting mode
 114     fCountMode = kCountAll;
 115     // Pycel
 116     SetPycellParameters();
 117     fSetNuclei = kFALSE;
 118 }
 119
 120 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 121     :AliGenMC(Pythia)
 122 {
 123 // copy constructor
 124     Pythia.Copy(*this);
 125 }
 126
 127 AliGenPythia::~AliGenPythia()
 128 {
 129 // Destructor
 130 }
 131
 132 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 133                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 134 {
 135 // Set pycell parameters
 136     fPycellEtaMax    =  etamax;
 137     fPycellNEta      =  neta;
 138     fPycellNPhi      =  nphi;
 139     fPycellThreshold =  thresh;
 140     fPycellEtSeed    =  etseed;
 141     fPycellMinEtJet  =  minet;
 142     fPycellMaxRadius =  r;
 143 }
 144
 145
 146
 147 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 148 {
 149   // Set a range of event numbers, for which a table
 150   // of generated particle will be printed
 151   fDebugEventFirst = eventFirst;
 152   fDebugEventLast  = eventLast;
 153   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 154 }
 155
 156 void AliGenPythia::Init()
 157 {
 158 // Initialisation
 159
 160     SetMC(AliPythia::Instance());
 161     fPythia=(AliPythia*) fMCEvGen;
 162
 163 //
 164     fParentWeight=1./Float_t(fNpart);
 165 //
 166 //  Forward Paramters to the AliPythia object
 167     fDecayer->SetForceDecay(fForceDecay);
 168     fDecayer->Init();
 169
 170
 171     fPythia->SetCKIN(3,fPtHardMin);
 172     fPythia->SetCKIN(4,fPtHardMax);
 173     fPythia->SetCKIN(7,fYHardMin);
 174     fPythia->SetCKIN(8,fYHardMax);
 175
 176     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
 177     // Fragmentation?
 178     if (fFragmentation) {
 179       fPythia->SetMSTP(111,1);
 180     } else {
 181       fPythia->SetMSTP(111,0);
 182     }
 183
 184
 185 //  initial state radiation
 186     fPythia->SetMSTP(61,fGinit);
 187 //  final state radiation
 188     fPythia->SetMSTP(71,fGfinal);
 189 //  pt - kick
 190     if (fPtKick > 0.) {
 191         fPythia->SetMSTP(91,1);
 192         fPythia->SetPARP(91,fPtKick);
 193     } else {
 194         fPythia->SetMSTP(91,0);
 195     }
 196
 197
 198     if (fReadFromFile) {
 199         fRL  =  AliRunLoader::Open(fFileName, "Partons");
 200         fRL->LoadKinematics();
 201         fRL->LoadHeader();
 202     } else {
 203         fRL = 0x0;
 204     }
 205
 206
 207  //
 208     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 209
 210 //  Parent and Children Selection
 211     switch (fProcess)
 212     {
 213     case kPyCharm:
 214     case kPyCharmUnforced:
 215     case kPyCharmPbPbMNR:
 216     case kPyCharmppMNR:
 217     case kPyCharmpPbMNR:
 218         fParentSelect[0] =   411;
 219         fParentSelect[1] =   421;
 220         fParentSelect[2] =   431;
 221         fParentSelect[3] =  4122;
 222         fFlavorSelect    =  4;
 223         break;
 224     case kPyD0PbPbMNR:
 225     case kPyD0pPbMNR:
 226     case kPyD0ppMNR:
 227         fParentSelect[0] =   421;
 228         fFlavorSelect    =   4;
 229         break;
 230     case kPyDPlusPbPbMNR:
 231     case kPyDPluspPbMNR:
 232     case kPyDPlusppMNR:
 233         fParentSelect[0] =   411;
 234         fFlavorSelect    =   4;
 235         break;
 236     case kPyBeauty:
 237     case kPyBeautyPbPbMNR:
 238     case kPyBeautypPbMNR:
 239     case kPyBeautyppMNR:
 240         fParentSelect[0]=  511;
 241         fParentSelect[1]=  521;
 242         fParentSelect[2]=  531;
 243         fParentSelect[3]= 5122;
 244         fParentSelect[4]= 5132;
 245         fParentSelect[5]= 5232;
 246         fParentSelect[6]= 5332;
 247         fFlavorSelect   = 5;
 248         break;
 249     case kPyBeautyUnforced:
 250         fParentSelect[0] =  511;
 251         fParentSelect[1] =  521;
 252         fParentSelect[2] =  531;
 253         fParentSelect[3] = 5122;
 254         fParentSelect[4] = 5132;
 255         fParentSelect[5] = 5232;
 256         fParentSelect[6] = 5332;
 257         fFlavorSelect    = 5;
 258         break;
 259     case kPyJpsiChi:
 260     case kPyJpsi:
 261         fParentSelect[0] = 443;
 262         break;
 263     case kPyMb:
 264     case kPyMbNonDiffr:
 265     case kPyJets:
 266     case kPyDirectGamma:
 267         break;
 268     case kPyW:
 269         break;
 270     }
 271 //
 272 //
 273 //  JetFinder for Trigger
 274 //
 275 //  Configure detector (EMCAL like)
 276 //
 277         fPythia->SetPARU(51, fPycellEtaMax);
 278         fPythia->SetMSTU(51, fPycellNEta);
 279         fPythia->SetMSTU(52, fPycellNPhi);
 280 //
 281 //  Configure Jet Finder
 282 //
 283         fPythia->SetPARU(58,  fPycellThreshold);
 284         fPythia->SetPARU(52,  fPycellEtSeed);
 285         fPythia->SetPARU(53,  fPycellMinEtJet);
 286         fPythia->SetPARU(54,  fPycellMaxRadius);
 287         fPythia->SetMSTU(54,  2);
 288 //
 289 //  This counts the total number of calls to Pyevnt() per run.
 290     fTrialsRun = 0;
 291     fQ         = 0.;
 292     fX1        = 0.;
 293     fX2        = 0.;
 294     fNev       = 0 ;
 295 //
 296 //
 297 //
 298     AliGenMC::Init();
 299 //
 300 //
 301 //
 302     if (fSetNuclei) {
 303         fDyBoost = 0;
 304         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 305     }
 306
 307     if (fQuench) {
 308         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 309     }
 310
 311 }
 312
 313 void AliGenPythia::Generate()
 314 {
 315 // Generate one event
 316
 317     fDecayer->ForceDecay();
 318
 319     Float_t polar[3]   =   {0,0,0};
 320     Float_t origin[3]  =   {0,0,0};
 321     Float_t p[4];
 322 //  converts from mm/c to s
 323     const Float_t kconv=0.001/2.999792458e8;
 324 //
 325     Int_t nt=0;
 326     Int_t jev=0;
 327     Int_t j, kf;
 328     fTrials=0;
 329     fEventTime = 0.;
 330
 331
 332
 333     //  Set collision vertex position
 334     if (fVertexSmear == kPerEvent) Vertex();
 335
 336 //  event loop
 337     while(1)
 338     {
 339 //
 340 // Produce event
 341 //
 342 //
 343 // Switch hadronisation off
 344 //
 345         fPythia->SetMSTJ(1, 0);
 346 //
 347 // Either produce new event or read partons from file
 348 //
 349         if (!fReadFromFile) {
 350             fPythia->Pyevnt();
 351             fNpartons = fPythia->GetN();
 352         } else {
 353             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
 354             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
 355             fPythia->SetN(0);
 356             LoadEvent(fRL->Stack(), 0 , 1);
 357             fPythia->Pyedit(21);
 358         }
 359
 360 //
 361 //  Run quenching routine
 362 //
 363         if (fQuench == 1) {
 364             fPythia->Quench();
 365         } else if (fQuench == 2){
 366             fPythia->Pyquen(208., 0, 0.);
 367         }
 368 //
 369 // Switch hadronisation on
 370 //
 371         fPythia->SetMSTJ(1, 1);
 372 //
 373 // .. and perform hadronisation
 374 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 375         fPythia->Pyexec();
 376         fTrials++;
 377         fPythia->ImportParticles(fParticles,"All");
 378         Boost();
 379 //
 380 //
 381 //
 382         Int_t i;
 383
 384
 385         Int_t np = fParticles->GetEntriesFast();
 386
 387         if (np == 0) continue;
 388 //
 389
 390 //
 391         Int_t* pParent   = new Int_t[np];
 392         Int_t* pSelected = new Int_t[np];
 393         Int_t* trackIt   = new Int_t[np];
 394         for (i = 0; i < np; i++) {
 395             pParent[i]   = -1;
 396             pSelected[i] =  0;
 397             trackIt[i]   =  0;
 398         }
 399
 400         Int_t nc = 0;        // Total n. of selected particles
 401         Int_t nParents = 0;  // Selected parents
 402         Int_t nTkbles = 0;   // Trackable particles
 403         if (fProcess != kPyMb && fProcess != kPyJets &&
 404             fProcess != kPyDirectGamma &&
 405             fProcess != kPyMbNonDiffr  &&
 406             fProcess != kPyW) {
 407
 408             for (i = 0; i < np; i++) {
 409                 TParticle* iparticle = (TParticle *) fParticles->At(i);
 410                 Int_t ks = iparticle->GetStatusCode();
 411                 kf = CheckPDGCode(iparticle->GetPdgCode());
 412 // No initial state partons
 413                 if (ks==21) continue;
 414 //
 415 // Heavy Flavor Selection
 416 //
 417                 // quark ?
 418                 kf = TMath::Abs(kf);
 419                 Int_t kfl = kf;
 420                 // Resonance
 421
 422                 if (kfl > 100000) kfl %= 100000;
 423                 if (kfl > 10000)  kfl %= 10000;
 424                 // meson ?
 425                 if  (kfl > 10) kfl/=100;
 426                 // baryon
 427                 if (kfl > 10) kfl/=10;
 428                 Int_t ipa = iparticle->GetFirstMother()-1;
 429                 Int_t kfMo = 0;
 430 //
 431 // Establish mother daughter relation between heavy quarks and mesons
 432 //
 433                 if (kf >= fFlavorSelect && kf <= 6) {
 434                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 435                     if (idau > -1) {
 436                         TParticle* daughter = (TParticle *) fParticles->At(idau);
 437                         Int_t pdgD = daughter->GetPdgCode();
 438                         if (pdgD == 91 || pdgD == 92) {
 439                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 440                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 441                             for (Int_t j = jmin; j <= jmax; j++)
 442                                 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
 443                         } // is string or cluster
 444                     } // has daughter
 445                 } // heavy quark
 446
 447
 448                 if (ipa > -1) {
 449                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 450                     kfMo = TMath::Abs(mother->GetPdgCode());
 451                 }
 452
 453                 // What to keep in Stack?
 454                 Bool_t flavorOK = kFALSE;
 455                 Bool_t selectOK = kFALSE;
 456                 if (fFeedDownOpt) {
 457                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 458                 } else {
 459                     if (kfl > fFlavorSelect) {
 460                         nc = -1;
 461                         break;
 462                     }
 463                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 464                 }
 465                 switch (fStackFillOpt) {
 466                 case kFlavorSelection:
 467                     selectOK = kTRUE;
 468                     break;
 469                 case kParentSelection:
 470                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 471                     break;
 472                 }
 473                 if (flavorOK && selectOK) {
 474 //
 475 // Heavy flavor hadron or quark
 476 //
 477 // Kinematic seletion on final state heavy flavor mesons
 478                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 479                     {
 480                         continue;
 481                     }
 482                     pSelected[i] = 1;
 483                     if (ParentSelected(kf)) ++nParents; // Update parent count
 484 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 485                 } else {
 486 // Kinematic seletion on decay products
 487                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 488                         && !KinematicSelection(iparticle, 1))
 489                     {
 490                         continue;
 491                     }
 492 //
 493 // Decay products
 494 // Select if mother was selected and is not tracked
 495
 496                     if (pSelected[ipa] &&
 497                         !trackIt[ipa]  &&     // mother will be  tracked ?
 498                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 499                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 500                         kf   != 92)           // don't store string
 501                     {
 502 //
 503 // Semi-stable or de-selected: diselect decay products:
 504 //
 505 //
 506                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 507                         {
 508                             Int_t ipF = iparticle->GetFirstDaughter();
 509                             Int_t ipL = iparticle->GetLastDaughter();
 510                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 511                         }
 512 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 513                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 514                     }
 515                 }
 516                 if (pSelected[i] == -1) pSelected[i] = 0;
 517                 if (!pSelected[i]) continue;
 518                 // Count quarks only if you did not include fragmentation
 519                 if (fFragmentation && kf <= 10) continue;
 520
 521                 nc++;
 522 // Decision on tracking
 523                 trackIt[i] = 0;
 524 //
 525 // Track final state particle
 526                 if (ks == 1) trackIt[i] = 1;
 527 // Track semi-stable particles
 528                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 529 // Track particles selected by process if undecayed.
 530                 if (fForceDecay == kNoDecay) {
 531                     if (ParentSelected(kf)) trackIt[i] = 1;
 532                 } else {
 533                     if (ParentSelected(kf)) trackIt[i] = 0;
 534                 }
 535                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 536 //
 537 //
 538
 539             } // particle selection loop
 540             if (nc > 0) {
 541                 for (i = 0; i<np; i++) {
 542                     if (!pSelected[i]) continue;
 543                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 544                     kf = CheckPDGCode(iparticle->GetPdgCode());
 545                     Int_t ks = iparticle->GetStatusCode();
 546                     p[0] = iparticle->Px();
 547                     p[1] = iparticle->Py();
 548                     p[2] = iparticle->Pz();
 549                     p[3] = iparticle->Energy();
 550
 551                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 552                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 553                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 554
 555                     Float_t tof   = kconv*iparticle->T();
 556                     Int_t ipa     = iparticle->GetFirstMother()-1;
 557                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 558
 559                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 560                               p[0], p[1], p[2], p[3],
 561                               origin[0], origin[1], origin[2], tof,
 562                               polar[0], polar[1], polar[2],
 563                               kPPrimary, nt, 1., ks);
 564                     pParent[i] = nt;
 565                     KeepTrack(nt);
 566                 } //  PushTrack loop
 567             }
 568         } else {
 569             nc = GenerateMB();
 570         } // mb ?
 571
 572         GetSubEventTime();
 573
 574         if (pParent)   delete[] pParent;
 575         if (pSelected) delete[] pSelected;
 576         if (trackIt)   delete[] trackIt;
 577
 578         if (nc > 0) {
 579           switch (fCountMode) {
 580           case kCountAll:
 581             // printf(" Count all \n");
 582             jev += nc;
 583             break;
 584           case kCountParents:
 585             // printf(" Count parents \n");
 586             jev += nParents;
 587             break;
 588           case kCountTrackables:
 589             // printf(" Count trackable \n");
 590             jev += nTkbles;
 591             break;
 592           }
 593             if (jev >= fNpart || fNpart == -1) {
 594                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 595                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 596
 597                 fQ  += fPythia->GetVINT(51);
 598                 fX1 += fPythia->GetVINT(41);
 599                 fX2 += fPythia->GetVINT(42);
 600                 fTrialsRun += fTrials;
 601                 fNev++;
 602                 MakeHeader();
 603                 break;
 604             }
 605         }
 606     } // event loop
 607     SetHighWaterMark(nt);
 608 //  adjust weight due to kinematic selection
 609 //    AdjustWeights();
 610 //  get cross-section
 611     fXsection=fPythia->GetPARI(1);
 612 }
 613
 614 Int_t  AliGenPythia::GenerateMB()
 615 {
 616 //
 617 // Min Bias selection and other global selections
 618 //
 619     Int_t i, kf, nt, iparent;
 620     Int_t nc = 0;
 621     Float_t p[4];
 622     Float_t polar[3]   =   {0,0,0};
 623     Float_t origin[3]  =   {0,0,0};
 624 //  converts from mm/c to s
 625     const Float_t kconv=0.001/2.999792458e8;
 626
 627
 628
 629     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 630
 631
 632     Int_t* pParent = new Int_t[np];
 633     for (i=0; i< np; i++) pParent[i] = -1;
 634     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 635         TParticle* jet1 = (TParticle *) fParticles->At(6);
 636         TParticle* jet2 = (TParticle *) fParticles->At(7);
 637         if (!CheckTrigger(jet1, jet2)) return 0;
 638     }
 639
 640     for (i = 0; i < np; i++) {
 641         Int_t trackIt = 0;
 642         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 643         kf = CheckPDGCode(iparticle->GetPdgCode());
 644         Int_t ks = iparticle->GetStatusCode();
 645         Int_t km = iparticle->GetFirstMother();
 646         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 647             (ks != 1) ||
 648             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 649             nc++;
 650             if (ks == 1) trackIt = 1;
 651             Int_t ipa = iparticle->GetFirstMother()-1;
 652
 653             iparent = (ipa > -1) ? pParent[ipa] : -1;
 654
 655 //
 656 // store track information
 657             p[0] = iparticle->Px();
 658             p[1] = iparticle->Py();
 659             p[2] = iparticle->Pz();
 660             p[3] = iparticle->Energy();
 661
 662
 663             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 664             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 665             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 666
 667             Float_t tof = fEventTime + kconv * iparticle->T();
 668
 669             PushTrack(fTrackIt*trackIt, iparent, kf,
 670                       p[0], p[1], p[2], p[3],
 671                       origin[0], origin[1], origin[2], tof,
 672                       polar[0], polar[1], polar[2],
 673                       kPPrimary, nt, 1., ks);
 674             //
 675             // Special Treatment to store color-flow
 676             //
 677             if (ks == 3 || ks == 13 || ks == 14) {
 678                 TParticle* particle = 0;
 679                 if (fStack) {
 680                     particle = fStack->Particle(nt);
 681                 } else {
 682                     particle = gAlice->Stack()->Particle(nt);
 683                 }
 684                 particle->SetFirstDaughter(fPythia->GetK(2, i));
 685                 particle->SetLastDaughter(fPythia->GetK(3, i));
 686             }
 687
 688             KeepTrack(nt);
 689             pParent[i] = nt;
 690             SetHighWaterMark(nt);
 691
 692         } // select particle
 693     } // particle loop
 694
 695     if (pParent) delete[] pParent;
 696
 697     printf("\n I've put %i particles on the stack \n",nc);
 698     return 1;
 699 }
 700
 701
 702 void AliGenPythia::FinishRun()
 703 {
 704 // Print x-section summary
 705     fPythia->Pystat(1);
 706
 707     if (fNev > 0.) {
 708         fQ  /= fNev;
 709         fX1 /= fNev;
 710         fX2 /= fNev;
 711     }
 712
 713     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 714     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 715 }
 716
 717 void AliGenPythia::AdjustWeights()
 718 {
 719 // Adjust the weights after generation of all events
 720 //
 721     if (gAlice) {
 722         TParticle *part;
 723         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 724         for (Int_t i=0; i<ntrack; i++) {
 725             part= gAlice->GetMCApp()->Particle(i);
 726             part->SetWeight(part->GetWeight()*fKineBias);
 727         }
 728     }
 729 }
 730
 731 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 732 {
 733 // Treat protons as inside nuclei with mass numbers a1 and a2
 734
 735     fAProjectile = a1;
 736     fATarget     = a2;
 737     fSetNuclei   = kTRUE;
 738 }
 739
 740
 741 void AliGenPythia::MakeHeader()
 742 {
 743   if (gAlice) {
 744     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 745         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
 746   }
 747
 748 // Builds the event header, to be called after each event
 749     if (fHeader) delete fHeader;
 750     fHeader = new AliGenPythiaEventHeader("Pythia");
 751 //
 752 // Event type
 753     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
 754 //
 755 // Number of trials
 756     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
 757 //
 758 // Event Vertex
 759     fHeader->SetPrimaryVertex(fVertex);
 760 //
 761 // Jets that have triggered
 762
 763     if (fProcess == kPyJets)
 764     {
 765         Int_t ntrig, njet;
 766         Float_t jets[4][10];
 767         GetJets(njet, ntrig, jets);
 768
 769
 770         for (Int_t i = 0; i < ntrig; i++) {
 771             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
 772                                                         jets[3][i]);
 773         }
 774     }
 775 //
 776 // Copy relevant information from external header, if present.
 777 //
 778     Float_t uqJet[4];
 779
 780     if (fRL) {
 781         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
 782         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
 783         {
 784             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
 785
 786
 787             exHeader->TriggerJet(i, uqJet);
 788             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
 789         }
 790     }
 791 //
 792 // Store quenching parameters
 793 //
 794     if (fQuench){
 795         Double_t z[4];
 796         Double_t xp, yp;
 797         if (fQuench == 1) {
 798             // Pythia::Quench()
 799             fPythia->GetQuenchingParameters(xp, yp, z);
 800         } else {
 801             // Pyquen
 802             Double_t r1 = PARIMP.rb1;
 803             Double_t r2 = PARIMP.rb2;
 804             Double_t b  = PARIMP.b1;
 805             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
 806             Double_t phi = PARIMP.psib1;
 807             xp = r * TMath::Cos(phi);
 808             yp = r * TMath::Sin(phi);
 809
 810         }
 811             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
 812             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
 813         }
 814
 815 //
 816 //  Pass header to RunLoader
 817 //
 818     AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(fHeader);
 819 }
 820
 821
 822 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 823 {
 824 // Check the kinematic trigger condition
 825 //
 826     Double_t eta[2];
 827     eta[0] = jet1->Eta();
 828     eta[1] = jet2->Eta();
 829     Double_t phi[2];
 830     phi[0] = jet1->Phi();
 831     phi[1] = jet2->Phi();
 832     Int_t    pdg[2];
 833     pdg[0] = jet1->GetPdgCode();
 834     pdg[1] = jet2->GetPdgCode();
 835     Bool_t   triggered = kFALSE;
 836
 837     if (fProcess == kPyJets) {
 838         Int_t njets = 0;
 839         Int_t ntrig = 0;
 840         Float_t jets[4][10];
 841 //
 842 // Use Pythia clustering on parton level to determine jet axis
 843 //
 844         GetJets(njets, ntrig, jets);
 845
 846         if (ntrig) triggered = kTRUE;
 847 //
 848     } else {
 849         Int_t ij = 0;
 850         Int_t ig = 1;
 851         if (pdg[0] == kGamma) {
 852             ij = 1;
 853             ig = 0;
 854         }
 855         //Check eta range first...
 856         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
 857             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
 858         {
 859             //Eta is okay, now check phi range
 860             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
 861                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
 862             {
 863                 triggered = kTRUE;
 864             }
 865         }
 866     }
 867     return triggered;
 868 }
 869
 870 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 871 {
 872 // Assignment operator
 873     rhs.Copy(*this);
 874     return *this;
 875 }
 876
 877 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
 878 {
 879 //
 880 // Load event into Pythia Common Block
 881 //
 882
 883     Int_t npart = stack -> GetNprimary();
 884     Int_t n0 = 0;
 885
 886     if (!flag) {
 887         (fPythia->GetPyjets())->N = npart;
 888     } else {
 889         n0 = (fPythia->GetPyjets())->N;
 890         (fPythia->GetPyjets())->N = n0 + npart;
 891     }
 892
 893
 894     for (Int_t part = 0; part < npart; part++) {
 895         TParticle *MPart = stack->Particle(part);
 896
 897         Int_t kf     =  MPart->GetPdgCode();
 898         Int_t ks     =  MPart->GetStatusCode();
 899         Int_t idf    =  MPart->GetFirstDaughter();
 900         Int_t idl    =  MPart->GetLastDaughter();
 901
 902         if (reHadr) {
 903             if (ks == 11 || ks == 12) {
 904                 ks  -= 10;
 905                 idf  = -1;
 906                 idl  = -1;
 907             }
 908         }
 909
 910         Float_t px = MPart->Px();
 911         Float_t py = MPart->Py();
 912         Float_t pz = MPart->Pz();
 913         Float_t e  = MPart->Energy();
 914         Float_t m  = MPart->GetCalcMass();
 915
 916
 917         (fPythia->GetPyjets())->P[0][part+n0] = px;
 918         (fPythia->GetPyjets())->P[1][part+n0] = py;
 919         (fPythia->GetPyjets())->P[2][part+n0] = pz;
 920         (fPythia->GetPyjets())->P[3][part+n0] = e;
 921         (fPythia->GetPyjets())->P[4][part+n0] = m;
 922
 923         (fPythia->GetPyjets())->K[1][part+n0] = kf;
 924         (fPythia->GetPyjets())->K[0][part+n0] = ks;
 925         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
 926         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
 927         (fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
 928     }
 929 }
 930
 931
 932 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
 933 {
 934 //
 935 //  Calls the Pythia jet finding algorithm to find jets in the current event
 936 //
 937 //
 938 //
 939 //  Save jets
 940     Int_t n     = fPythia->GetN();
 941
 942 //
 943 //  Run Jet Finder
 944     fPythia->Pycell(njets);
 945     Int_t i;
 946     for (i = 0; i < njets; i++) {
 947         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
 948         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
 949         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 950         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 951
 952         jets[0][i] = px;
 953         jets[1][i] = py;
 954         jets[2][i] = pz;
 955         jets[3][i] = e;
 956     }
 957 }
 958
 959
 960
 961 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 962 {
 963 //
 964 //  Calls the Pythia clustering algorithm to find jets in the current event
 965 //
 966     Int_t n     = fPythia->GetN();
 967     nJets       = 0;
 968     nJetsTrig   = 0;
 969     if (fJetReconstruction == kCluster) {
 970 //
 971 //  Configure cluster algorithm
 972 //
 973         fPythia->SetPARU(43, 2.);
 974         fPythia->SetMSTU(41, 1);
 975 //
 976 //  Call cluster algorithm
 977 //
 978         fPythia->Pyclus(nJets);
 979 //
 980 //  Loading jets from common block
 981 //
 982     } else {
 983
 984 //
 985 //  Run Jet Finder
 986         fPythia->Pycell(nJets);
 987     }
 988
 989     Int_t i;
 990     for (i = 0; i < nJets; i++) {
 991         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
 992         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
 993         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 994         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 995         Float_t pt    = TMath::Sqrt(px * px + py * py);
 996         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
 997         Float_t theta = TMath::ATan2(pt,pz);
 998         Float_t et    = e * TMath::Sin(theta);
 999         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1000         if (
1001             eta > fEtaMinJet && eta < fEtaMaxJet &&
1002             phi > fPhiMinJet && phi < fPhiMaxJet &&
1003             et  > fEtMinJet  && et  < fEtMaxJet
1004             )
1005         {
1006             jets[0][nJetsTrig] = px;
1007             jets[1][nJetsTrig] = py;
1008             jets[2][nJetsTrig] = pz;
1009             jets[3][nJetsTrig] = e;
1010             nJetsTrig++;
1011 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1012         } else {
1013 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1014         }
1015     }
1016 }
1017
1018 void AliGenPythia::GetSubEventTime()
1019 {
1020   // Calculates time of the next subevent
1021   fEventTime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
1022 }
1023
1024 #ifdef never
1025 void AliGenPythia::Streamer(TBuffer &R__b)
1026 {
1027    // Stream an object of class AliGenPythia.
1028
1029    if (R__b.IsReading()) {
1030       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1031       AliGenerator::Streamer(R__b);
1032       R__b >> (Int_t&)fProcess;
1033       R__b >> (Int_t&)fStrucFunc;
1034       R__b >> (Int_t&)fForceDecay;
1035       R__b >> fEnergyCMS;
1036       R__b >> fKineBias;
1037       R__b >> fTrials;
1038       fParentSelect.Streamer(R__b);
1039       fChildSelect.Streamer(R__b);
1040       R__b >> fXsection;
1041 //      (AliPythia::Instance())->Streamer(R__b);
1042       R__b >> fPtHardMin;
1043       R__b >> fPtHardMax;
1044 //      if (fDecayer) fDecayer->Streamer(R__b);
1045    } else {
1046       R__b.WriteVersion(AliGenPythia::IsA());
1047       AliGenerator::Streamer(R__b);
1048       R__b << (Int_t)fProcess;
1049       R__b << (Int_t)fStrucFunc;
1050       R__b << (Int_t)fForceDecay;
1051       R__b << fEnergyCMS;
1052       R__b << fKineBias;
1053       R__b << fTrials;
1054       fParentSelect.Streamer(R__b);
1055       fChildSelect.Streamer(R__b);
1056       R__b << fXsection;
1057 //      R__b << fPythia;
1058       R__b << fPtHardMin;
1059       R__b << fPtHardMax;
1060       //     fDecayer->Streamer(R__b);
1061    }
1062 }
1063 #endif
1064
1065
1066