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
  49 AliGenPythia::AliGenPythia():
  50     AliGenMC(),
  51     fProcess(kPyCharm),
  52     fStrucFunc(kCTEQ5L),
  53     fEnergyCMS(5500.),
  54     fKineBias(0.),
  55     fTrials(0),
  56     fTrialsRun(0),
  57     fQ(0.),
  58     fX1(0.),
  59     fX2(0.),
  60     fEventTime(0.),
  61     fInteractionRate(0.),
  62     fTimeWindow(0.),
  63     fCurSubEvent(0),
  64     fEventsTime(0),
  65     fNev(0),
  66     fFlavorSelect(0),
  67     fXsection(0.),
  68     fPythia(0),
  69     fPtHardMin(0.),
  70     fPtHardMax(1.e4),
  71     fYHardMin(-1.e10),
  72     fYHardMax(1.e10),
  73     fGinit(1),
  74     fGfinal(1),
  75     fHadronisation(1),
  76     fNpartons(0),
  77     fReadFromFile(0),
  78     fQuench(0),
  79     fPtKick(1.),
  80     fFullEvent(kTRUE),
  81     fDecayer(new AliDecayerPythia()),
  82     fDebugEventFirst(-1),
  83     fDebugEventLast(-1),
  84     fEtMinJet(0.),
  85     fEtMaxJet(1.e4),
  86     fEtaMinJet(-20.),
  87     fEtaMaxJet(20.),
  88     fPhiMinJet(0.),
  89     fPhiMaxJet(2.* TMath::Pi()),
  90     fJetReconstruction(kCell),
  91     fEtaMinGamma(-20.),
  92     fEtaMaxGamma(20.),
  93     fPhiMinGamma(0.),
  94     fPhiMaxGamma(2. * TMath::Pi()),
  95     fPycellEtaMax(2.),
  96     fPycellNEta(274),
  97     fPycellNPhi(432),
  98     fPycellThreshold(0.),
  99     fPycellEtSeed(4.),
 100     fPycellMinEtJet(10.),
 101     fPycellMaxRadius(1.),
 102     fStackFillOpt(kFlavorSelection),
 103     fFeedDownOpt(kTRUE),
 104     fFragmentation(kTRUE),
 105     fSetNuclei(kFALSE),
 106     fNewMIS(kFALSE),
 107     fHFoff(kFALSE),
 108     fTriggerParticle(0),
 109     fTriggerEta(0.9),
 110     fCountMode(kCountAll),
 111     fHeader(0),
 112     fRL(0),
 113     fFileName(0),
 114     fFragPhotonInCalo(kFALSE),
 115     fPi0InCalo(kFALSE) ,
 116     fCheckEMCAL(kFALSE),
 117     fCheckPHOS(kFALSE),
 118     fFragPhotonOrPi0MinPt(0),
 119     fPHOSMinPhi(219.),
 120     fPHOSMaxPhi(321.),
 121     fPHOSEta(0.13),
 122     fEMCALMinPhi(79.),
 123     fEMCALMaxPhi(191.),
 124     fEMCALEta(0.71)
 125
 126 {
 127 // Default Constructor
 128   SetNuclei(0,0);
 129   if (!AliPythiaRndm::GetPythiaRandom())
 130       AliPythiaRndm::SetPythiaRandom(GetRandom());
 131 }
 132
 133 AliGenPythia::AliGenPythia(Int_t npart)
 134     :AliGenMC(npart),
 135      fProcess(kPyCharm),
 136      fStrucFunc(kCTEQ5L),
 137      fEnergyCMS(5500.),
 138      fKineBias(0.),
 139      fTrials(0),
 140      fTrialsRun(0),
 141      fQ(0.),
 142      fX1(0.),
 143      fX2(0.),
 144      fEventTime(0.),
 145      fInteractionRate(0.),
 146      fTimeWindow(0.),
 147      fCurSubEvent(0),
 148      fEventsTime(0),
 149      fNev(0),
 150      fFlavorSelect(0),
 151      fXsection(0.),
 152      fPythia(0),
 153      fPtHardMin(0.),
 154      fPtHardMax(1.e4),
 155      fYHardMin(-1.e10),
 156      fYHardMax(1.e10),
 157      fGinit(kTRUE),
 158      fGfinal(kTRUE),
 159      fHadronisation(kTRUE),
 160      fNpartons(0),
 161      fReadFromFile(kFALSE),
 162      fQuench(kFALSE),
 163      fPtKick(1.),
 164      fFullEvent(kTRUE),
 165      fDecayer(new AliDecayerPythia()),
 166      fDebugEventFirst(-1),
 167      fDebugEventLast(-1),
 168      fEtMinJet(0.),
 169      fEtMaxJet(1.e4),
 170      fEtaMinJet(-20.),
 171      fEtaMaxJet(20.),
 172      fPhiMinJet(0.),
 173      fPhiMaxJet(2.* TMath::Pi()),
 174      fJetReconstruction(kCell),
 175      fEtaMinGamma(-20.),
 176      fEtaMaxGamma(20.),
 177      fPhiMinGamma(0.),
 178      fPhiMaxGamma(2. * TMath::Pi()),
 179      fPycellEtaMax(2.),
 180      fPycellNEta(274),
 181      fPycellNPhi(432),
 182      fPycellThreshold(0.),
 183      fPycellEtSeed(4.),
 184      fPycellMinEtJet(10.),
 185      fPycellMaxRadius(1.),
 186      fStackFillOpt(kFlavorSelection),
 187      fFeedDownOpt(kTRUE),
 188      fFragmentation(kTRUE),
 189      fSetNuclei(kFALSE),
 190      fNewMIS(kFALSE),
 191      fHFoff(kFALSE),
 192      fTriggerParticle(0),
 193      fTriggerEta(0.9),
 194      fCountMode(kCountAll),
 195      fHeader(0),
 196      fRL(0),
 197      fFileName(0),
 198      fFragPhotonInCalo(kFALSE),
 199      fPi0InCalo(kFALSE) ,
 200      fCheckEMCAL(kFALSE),
 201      fCheckPHOS(kFALSE),
 202      fFragPhotonOrPi0MinPt(0),
 203      fPHOSMinPhi(219.),
 204      fPHOSMaxPhi(321.),
 205      fPHOSEta(0.13),
 206      fEMCALMinPhi(79.),
 207      fEMCALMaxPhi(191.),
 208      fEMCALEta(0.71)
 209 {
 210 // default charm production at 5. 5 TeV
 211 // semimuonic decay
 212 // structure function GRVHO
 213 //
 214     fName = "Pythia";
 215     fTitle= "Particle Generator using PYTHIA";
 216     SetForceDecay();
 217     // Set random number generator
 218     if (!AliPythiaRndm::GetPythiaRandom())
 219       AliPythiaRndm::SetPythiaRandom(GetRandom());
 220     fParticles = new TClonesArray("TParticle",1000);
 221     SetNuclei(0,0);
 222  }
 223
 224 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 225     :AliGenMC(Pythia),
 226      fProcess(kPyCharm),
 227      fStrucFunc(kCTEQ5L),
 228      fEnergyCMS(5500.),
 229      fKineBias(0.),
 230      fTrials(0),
 231      fTrialsRun(0),
 232      fQ(0.),
 233      fX1(0.),
 234      fX2(0.),
 235      fEventTime(0.),
 236      fInteractionRate(0.),
 237      fTimeWindow(0.),
 238      fCurSubEvent(0),
 239      fEventsTime(0),
 240      fNev(0),
 241      fFlavorSelect(0),
 242      fXsection(0.),
 243      fPythia(0),
 244      fPtHardMin(0.),
 245      fPtHardMax(1.e4),
 246      fYHardMin(-1.e10),
 247      fYHardMax(1.e10),
 248      fGinit(kTRUE),
 249      fGfinal(kTRUE),
 250      fHadronisation(kTRUE),
 251      fNpartons(0),
 252      fReadFromFile(kFALSE),
 253      fQuench(kFALSE),
 254      fPtKick(1.),
 255      fFullEvent(kTRUE),
 256      fDecayer(new AliDecayerPythia()),
 257      fDebugEventFirst(-1),
 258      fDebugEventLast(-1),
 259      fEtMinJet(0.),
 260      fEtMaxJet(1.e4),
 261      fEtaMinJet(-20.),
 262      fEtaMaxJet(20.),
 263      fPhiMinJet(0.),
 264      fPhiMaxJet(2.* TMath::Pi()),
 265      fJetReconstruction(kCell),
 266      fEtaMinGamma(-20.),
 267      fEtaMaxGamma(20.),
 268      fPhiMinGamma(0.),
 269      fPhiMaxGamma(2. * TMath::Pi()),
 270      fPycellEtaMax(2.),
 271      fPycellNEta(274),
 272      fPycellNPhi(432),
 273      fPycellThreshold(0.),
 274      fPycellEtSeed(4.),
 275      fPycellMinEtJet(10.),
 276      fPycellMaxRadius(1.),
 277      fStackFillOpt(kFlavorSelection),
 278      fFeedDownOpt(kTRUE),
 279      fFragmentation(kTRUE),
 280      fSetNuclei(kFALSE),
 281      fNewMIS(kFALSE),
 282      fHFoff(kFALSE),
 283      fTriggerParticle(0),
 284      fTriggerEta(0.9),
 285      fCountMode(kCountAll),
 286      fHeader(0),
 287      fRL(0),
 288      fFileName(0),
 289      fFragPhotonInCalo(kFALSE),
 290      fPi0InCalo(kFALSE) ,
 291      fCheckEMCAL(kFALSE),
 292      fCheckPHOS(kFALSE),
 293      fFragPhotonOrPi0MinPt(0),
 294      fPHOSMinPhi(219.),
 295      fPHOSMaxPhi(321.),
 296      fPHOSEta(0.13),
 297      fEMCALMinPhi(79.),
 298      fEMCALMaxPhi(191.),
 299      fEMCALEta(0.71)
 300 {
 301 // copy constructor
 302     Pythia.Copy(*this);
 303 }
 304
 305 AliGenPythia::~AliGenPythia()
 306 {
 307 // Destructor
 308   if(fEventsTime) delete fEventsTime;
 309 }
 310
 311 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
 312 {
 313 // Generate pileup using user specified rate
 314     fInteractionRate = rate;
 315     fTimeWindow = timewindow;
 316     GeneratePileup();
 317 }
 318
 319 void AliGenPythia::GeneratePileup()
 320 {
 321 // Generate sub events time for pileup
 322     fEventsTime = 0;
 323     if(fInteractionRate == 0.) {
 324       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
 325       return;
 326     }
 327
 328     Int_t npart = NumberParticles();
 329     if(npart < 0) {
 330       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
 331       return;
 332     }
 333
 334     if(fEventsTime) delete fEventsTime;
 335     fEventsTime = new TArrayF(npart);
 336     TArrayF &array = *fEventsTime;
 337     for(Int_t ipart = 0; ipart < npart; ipart++)
 338       array[ipart] = 0.;
 339
 340     Float_t eventtime = 0.;
 341     while(1)
 342       {
 343         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 344         if(eventtime > fTimeWindow) break;
 345         array.Set(array.GetSize()+1);
 346         array[array.GetSize()-1] = eventtime;
 347       }
 348
 349     eventtime = 0.;
 350     while(1)
 351       {
 352         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 353         if(TMath::Abs(eventtime) > fTimeWindow) break;
 354         array.Set(array.GetSize()+1);
 355         array[array.GetSize()-1] = eventtime;
 356       }
 357
 358     SetNumberParticles(fEventsTime->GetSize());
 359 }
 360
 361 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 362                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 363 {
 364 // Set pycell parameters
 365     fPycellEtaMax    =  etamax;
 366     fPycellNEta      =  neta;
 367     fPycellNPhi      =  nphi;
 368     fPycellThreshold =  thresh;
 369     fPycellEtSeed    =  etseed;
 370     fPycellMinEtJet  =  minet;
 371     fPycellMaxRadius =  r;
 372 }
 373
 374
 375
 376 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 377 {
 378   // Set a range of event numbers, for which a table
 379   // of generated particle will be printed
 380   fDebugEventFirst = eventFirst;
 381   fDebugEventLast  = eventLast;
 382   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 383 }
 384
 385 void AliGenPythia::Init()
 386 {
 387 // Initialisation
 388
 389     SetMC(AliPythia::Instance());
 390     fPythia=(AliPythia*) fMCEvGen;
 391
 392 //
 393     fParentWeight=1./Float_t(fNpart);
 394 //
 395 //  Forward Paramters to the AliPythia object
 396     fDecayer->SetForceDecay(fForceDecay);
 397     fDecayer->Init();
 398
 399
 400     fPythia->SetCKIN(3,fPtHardMin);
 401     fPythia->SetCKIN(4,fPtHardMax);
 402     fPythia->SetCKIN(7,fYHardMin);
 403     fPythia->SetCKIN(8,fYHardMax);
 404
 405     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
 406     // Fragmentation?
 407     if (fFragmentation) {
 408       fPythia->SetMSTP(111,1);
 409     } else {
 410       fPythia->SetMSTP(111,0);
 411     }
 412
 413
 414 //  initial state radiation
 415     fPythia->SetMSTP(61,fGinit);
 416 //  final state radiation
 417     fPythia->SetMSTP(71,fGfinal);
 418 //  pt - kick
 419     if (fPtKick > 0.) {
 420         fPythia->SetMSTP(91,1);
 421         fPythia->SetPARP(91,fPtKick);
 422     } else {
 423         fPythia->SetMSTP(91,0);
 424     }
 425
 426
 427     if (fReadFromFile) {
 428         fRL  =  AliRunLoader::Open(fFileName, "Partons");
 429         fRL->LoadKinematics();
 430         fRL->LoadHeader();
 431     } else {
 432         fRL = 0x0;
 433     }
 434 // Switch off Heavy Flavors on request
 435     if (fHFoff) {
 436         fPythia->SetMSTP(58, 3);
 437         fPythia->SetMSTJ(45, 3);
 438         for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
 439     }
 440  //
 441     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 442
 443 //  Parent and Children Selection
 444     switch (fProcess)
 445     {
 446     case kPyOldUEQ2ordered:
 447     case kPyOldUEQ2ordered2:
 448     case kPyOldPopcorn:
 449       break;
 450     case kPyCharm:
 451     case kPyCharmUnforced:
 452     case kPyCharmPbPbMNR:
 453     case kPyCharmpPbMNR:
 454     case kPyCharmppMNR:
 455     case kPyCharmppMNRwmi:
 456         fParentSelect[0] =   411;
 457         fParentSelect[1] =   421;
 458         fParentSelect[2] =   431;
 459         fParentSelect[3] =  4122;
 460         fFlavorSelect    =  4;
 461         break;
 462     case kPyD0PbPbMNR:
 463     case kPyD0pPbMNR:
 464     case kPyD0ppMNR:
 465         fParentSelect[0] =   421;
 466         fFlavorSelect    =   4;
 467         break;
 468     case kPyDPlusPbPbMNR:
 469     case kPyDPluspPbMNR:
 470     case kPyDPlusppMNR:
 471         fParentSelect[0] =   411;
 472         fFlavorSelect    =   4;
 473         break;
 474     case kPyDPlusStrangePbPbMNR:
 475     case kPyDPlusStrangepPbMNR:
 476     case kPyDPlusStrangeppMNR:
 477         fParentSelect[0] =   431;
 478         fFlavorSelect    =   4;
 479         break;
 480     case kPyBeauty:
 481     case kPyBeautyPbPbMNR:
 482     case kPyBeautypPbMNR:
 483     case kPyBeautyppMNR:
 484     case kPyBeautyppMNRwmi:
 485         fParentSelect[0]=  511;
 486         fParentSelect[1]=  521;
 487         fParentSelect[2]=  531;
 488         fParentSelect[3]= 5122;
 489         fParentSelect[4]= 5132;
 490         fParentSelect[5]= 5232;
 491         fParentSelect[6]= 5332;
 492         fFlavorSelect   = 5;
 493         break;
 494     case kPyBeautyUnforced:
 495         fParentSelect[0] =  511;
 496         fParentSelect[1] =  521;
 497         fParentSelect[2] =  531;
 498         fParentSelect[3] = 5122;
 499         fParentSelect[4] = 5132;
 500         fParentSelect[5] = 5232;
 501         fParentSelect[6] = 5332;
 502         fFlavorSelect    = 5;
 503         break;
 504     case kPyJpsiChi:
 505     case kPyJpsi:
 506         fParentSelect[0] = 443;
 507         break;
 508     case kPyMb:
 509     case kPyMbNonDiffr:
 510     case kPyMbMSEL1:
 511     case kPyJets:
 512     case kPyDirectGamma:
 513     case kPyLhwgMb:
 514         break;
 515     case kPyW:
 516     case kPyZ:
 517         break;
 518     }
 519 //
 520 //
 521 //  JetFinder for Trigger
 522 //
 523 //  Configure detector (EMCAL like)
 524 //
 525     fPythia->SetPARU(51, fPycellEtaMax);
 526     fPythia->SetMSTU(51, fPycellNEta);
 527     fPythia->SetMSTU(52, fPycellNPhi);
 528 //
 529 //  Configure Jet Finder
 530 //
 531     fPythia->SetPARU(58,  fPycellThreshold);
 532     fPythia->SetPARU(52,  fPycellEtSeed);
 533     fPythia->SetPARU(53,  fPycellMinEtJet);
 534     fPythia->SetPARU(54,  fPycellMaxRadius);
 535     fPythia->SetMSTU(54,  2);
 536 //
 537 //  This counts the total number of calls to Pyevnt() per run.
 538     fTrialsRun = 0;
 539     fQ         = 0.;
 540     fX1        = 0.;
 541     fX2        = 0.;
 542     fNev       = 0 ;
 543 //
 544 //
 545 //
 546     AliGenMC::Init();
 547 //
 548 //
 549 //
 550     if (fSetNuclei) {
 551         fDyBoost = 0;
 552         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 553     }
 554
 555     if (fQuench) {
 556         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 557     }
 558     fPythia->SetPARJ(200, 0.0);
 559
 560     if (fQuench == 3) {
 561         // Nestor's change of the splittings
 562         fPythia->SetPARJ(200, 0.8);
 563         fPythia->SetMSTJ(41, 1);  // QCD radiation only
 564         fPythia->SetMSTJ(42, 2);  // angular ordering
 565         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 566         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 567         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 568         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 569     }
 570 }
 571
 572 void AliGenPythia::Generate()
 573 {
 574 // Generate one event
 575
 576     fDecayer->ForceDecay();
 577
 578     Float_t polar[3]   =   {0,0,0};
 579     Float_t origin[3]  =   {0,0,0};
 580     Float_t p[4];
 581 //  converts from mm/c to s
 582     const Float_t kconv=0.001/2.999792458e8;
 583 //
 584     Int_t nt=0;
 585     Int_t jev=0;
 586     Int_t j, kf;
 587     fTrials=0;
 588     fEventTime = 0.;
 589
 590
 591
 592     //  Set collision vertex position
 593     if (fVertexSmear == kPerEvent) Vertex();
 594
 595 //  event loop
 596     while(1)
 597     {
 598 //
 599 // Produce event
 600 //
 601 //
 602 // Switch hadronisation off
 603 //
 604         fPythia->SetMSTJ(1, 0);
 605 //
 606 // Either produce new event or read partons from file
 607 //
 608         if (!fReadFromFile) {
 609             if (!fNewMIS) {
 610                 fPythia->Pyevnt();
 611             } else {
 612                 fPythia->Pyevnw();
 613             }
 614             fNpartons = fPythia->GetN();
 615         } else {
 616             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
 617             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
 618             fPythia->SetN(0);
 619             LoadEvent(fRL->Stack(), 0 , 1);
 620             fPythia->Pyedit(21);
 621         }
 622
 623 //
 624 //  Run quenching routine
 625 //
 626         if (fQuench == 1) {
 627             fPythia->Quench();
 628         } else if (fQuench == 2){
 629             fPythia->Pyquen(208., 0, 0.);
 630         } else if (fQuench == 3) {
 631             // Quenching is via multiplicative correction of the splittings
 632         }
 633
 634 //
 635 // Switch hadronisation on
 636 //
 637         fPythia->SetMSTJ(1, 1);
 638 //
 639 // .. and perform hadronisation
 640 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 641         fPythia->Pyexec();
 642         fTrials++;
 643         fPythia->ImportParticles(fParticles,"All");
 644         Boost();
 645 //
 646 //
 647 //
 648         Int_t i;
 649
 650
 651         Int_t np = fParticles->GetEntriesFast();
 652
 653         if (np == 0) continue;
 654 //
 655
 656 //
 657         Int_t* pParent   = new Int_t[np];
 658         Int_t* pSelected = new Int_t[np];
 659         Int_t* trackIt   = new Int_t[np];
 660         for (i = 0; i < np; i++) {
 661             pParent[i]   = -1;
 662             pSelected[i] =  0;
 663             trackIt[i]   =  0;
 664         }
 665
 666         Int_t nc = 0;        // Total n. of selected particles
 667         Int_t nParents = 0;  // Selected parents
 668         Int_t nTkbles = 0;   // Trackable particles
 669         if (fProcess != kPyMb && fProcess != kPyJets &&
 670             fProcess != kPyDirectGamma &&
 671             fProcess != kPyMbNonDiffr  &&
 672             fProcess != kPyMbMSEL1     &&
 673             fProcess != kPyW && fProcess != kPyZ &&
 674             fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
 675
 676             for (i = 0; i < np; i++) {
 677                 TParticle* iparticle = (TParticle *) fParticles->At(i);
 678                 Int_t ks = iparticle->GetStatusCode();
 679                 kf = CheckPDGCode(iparticle->GetPdgCode());
 680 // No initial state partons
 681                 if (ks==21) continue;
 682 //
 683 // Heavy Flavor Selection
 684 //
 685                 // quark ?
 686                 kf = TMath::Abs(kf);
 687                 Int_t kfl = kf;
 688                 // Resonance
 689
 690                 if (kfl > 100000) kfl %= 100000;
 691                 if (kfl > 10000)  kfl %= 10000;
 692                 // meson ?
 693                 if  (kfl > 10) kfl/=100;
 694                 // baryon
 695                 if (kfl > 10) kfl/=10;
 696                 Int_t ipa = iparticle->GetFirstMother()-1;
 697                 Int_t kfMo = 0;
 698 //
 699 // Establish mother daughter relation between heavy quarks and mesons
 700 //
 701                 if (kf >= fFlavorSelect && kf <= 6) {
 702                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 703                     if (idau > -1) {
 704                         TParticle* daughter = (TParticle *) fParticles->At(idau);
 705                         Int_t pdgD = daughter->GetPdgCode();
 706                         if (pdgD == 91 || pdgD == 92) {
 707                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 708                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 709                             for (Int_t j = jmin; j <= jmax; j++)
 710                                 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
 711                         } // is string or cluster
 712                     } // has daughter
 713                 } // heavy quark
 714
 715
 716                 if (ipa > -1) {
 717                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 718                     kfMo = TMath::Abs(mother->GetPdgCode());
 719                 }
 720
 721                 // What to keep in Stack?
 722                 Bool_t flavorOK = kFALSE;
 723                 Bool_t selectOK = kFALSE;
 724                 if (fFeedDownOpt) {
 725                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 726                 } else {
 727                     if (kfl > fFlavorSelect) {
 728                         nc = -1;
 729                         break;
 730                     }
 731                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 732                 }
 733                 switch (fStackFillOpt) {
 734                 case kFlavorSelection:
 735                     selectOK = kTRUE;
 736                     break;
 737                 case kParentSelection:
 738                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 739                     break;
 740                 }
 741                 if (flavorOK && selectOK) {
 742 //
 743 // Heavy flavor hadron or quark
 744 //
 745 // Kinematic seletion on final state heavy flavor mesons
 746                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 747                     {
 748                         continue;
 749                     }
 750                     pSelected[i] = 1;
 751                     if (ParentSelected(kf)) ++nParents; // Update parent count
 752 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 753                 } else {
 754 // Kinematic seletion on decay products
 755                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 756                         && !KinematicSelection(iparticle, 1))
 757                     {
 758                         continue;
 759                     }
 760 //
 761 // Decay products
 762 // Select if mother was selected and is not tracked
 763
 764                     if (pSelected[ipa] &&
 765                         !trackIt[ipa]  &&     // mother will be  tracked ?
 766                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 767                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 768                         kf   != 92)           // don't store string
 769                     {
 770 //
 771 // Semi-stable or de-selected: diselect decay products:
 772 //
 773 //
 774                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 775                         {
 776                             Int_t ipF = iparticle->GetFirstDaughter();
 777                             Int_t ipL = iparticle->GetLastDaughter();
 778                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 779                         }
 780 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 781                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 782                     }
 783                 }
 784                 if (pSelected[i] == -1) pSelected[i] = 0;
 785                 if (!pSelected[i]) continue;
 786                 // Count quarks only if you did not include fragmentation
 787                 if (fFragmentation && kf <= 10) continue;
 788
 789                 nc++;
 790 // Decision on tracking
 791                 trackIt[i] = 0;
 792 //
 793 // Track final state particle
 794                 if (ks == 1) trackIt[i] = 1;
 795 // Track semi-stable particles
 796                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 797 // Track particles selected by process if undecayed.
 798                 if (fForceDecay == kNoDecay) {
 799                     if (ParentSelected(kf)) trackIt[i] = 1;
 800                 } else {
 801                     if (ParentSelected(kf)) trackIt[i] = 0;
 802                 }
 803                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 804 //
 805 //
 806
 807             } // particle selection loop
 808             if (nc > 0) {
 809                 for (i = 0; i<np; i++) {
 810                     if (!pSelected[i]) continue;
 811                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 812                     kf = CheckPDGCode(iparticle->GetPdgCode());
 813                     Int_t ks = iparticle->GetStatusCode();
 814                     p[0] = iparticle->Px();
 815                     p[1] = iparticle->Py();
 816                     p[2] = iparticle->Pz();
 817                     p[3] = iparticle->Energy();
 818
 819                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 820                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 821                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 822
 823                     Float_t tof   = kconv*iparticle->T();
 824                     Int_t ipa     = iparticle->GetFirstMother()-1;
 825                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 826
 827                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 828                               p[0], p[1], p[2], p[3],
 829                               origin[0], origin[1], origin[2], tof,
 830                               polar[0], polar[1], polar[2],
 831                               kPPrimary, nt, 1., ks);
 832                     pParent[i] = nt;
 833                     KeepTrack(nt);
 834                 } //  PushTrack loop
 835             }
 836         } else {
 837             nc = GenerateMB();
 838         } // mb ?
 839
 840         GetSubEventTime();
 841
 842         delete[] pParent;
 843         delete[] pSelected;
 844         delete[] trackIt;
 845
 846         if (nc > 0) {
 847           switch (fCountMode) {
 848           case kCountAll:
 849             // printf(" Count all \n");
 850             jev += nc;
 851             break;
 852           case kCountParents:
 853             // printf(" Count parents \n");
 854             jev += nParents;
 855             break;
 856           case kCountTrackables:
 857             // printf(" Count trackable \n");
 858             jev += nTkbles;
 859             break;
 860           }
 861             if (jev >= fNpart || fNpart == -1) {
 862                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 863
 864                 fQ  += fPythia->GetVINT(51);
 865                 fX1 += fPythia->GetVINT(41);
 866                 fX2 += fPythia->GetVINT(42);
 867                 fTrialsRun += fTrials;
 868                 fNev++;
 869                 MakeHeader();
 870                 break;
 871             }
 872         }
 873     } // event loop
 874     SetHighWaterMark(nt);
 875 //  adjust weight due to kinematic selection
 876 //    AdjustWeights();
 877 //  get cross-section
 878     fXsection=fPythia->GetPARI(1);
 879 }
 880
 881 Int_t  AliGenPythia::GenerateMB()
 882 {
 883 //
 884 // Min Bias selection and other global selections
 885 //
 886     Int_t i, kf, nt, iparent;
 887     Int_t nc = 0;
 888     Float_t p[4];
 889     Float_t polar[3]   =   {0,0,0};
 890     Float_t origin[3]  =   {0,0,0};
 891 //  converts from mm/c to s
 892     const Float_t kconv=0.001/2.999792458e8;
 893
 894
 895
 896     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 897
 898
 899
 900     Int_t* pParent = new Int_t[np];
 901     for (i=0; i< np; i++) pParent[i] = -1;
 902     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 903         TParticle* jet1 = (TParticle *) fParticles->At(6);
 904         TParticle* jet2 = (TParticle *) fParticles->At(7);
 905         if (!CheckTrigger(jet1, jet2)) {
 906           delete [] pParent;
 907           return 0;
 908         }
 909     }
 910
 911     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
 912     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
 913
 914       Bool_t ok = kFALSE;
 915
 916       Int_t pdg  = 0;
 917       if (fFragPhotonInCalo) pdg = 22   ; // Photon
 918       else if (fPi0InCalo) pdg = 111 ; // Pi0
 919
 920       for (i=0; i< np; i++) {
 921         TParticle* iparticle = (TParticle *) fParticles->At(i);
 922         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
 923            iparticle->Pt() > fFragPhotonOrPi0MinPt){
 924           if(pdg == 111 ||
 925              (pdg == 22 && iparticle->GetFirstMother() != 11))//No photon from hadron decay
 926             {
 927               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 928               Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
 929               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
 930                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
 931                 ok =kTRUE;
 932             }
 933         }
 934       }
 935       if(!ok)
 936         return 0;
 937     }
 938
 939     if (fTriggerParticle) {
 940         Bool_t triggered = kFALSE;
 941         for (i = 0; i < np; i++) {
 942             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 943             kf = CheckPDGCode(iparticle->GetPdgCode());
 944             if (kf != fTriggerParticle) continue;
 945             if (iparticle->Pt() == 0.) continue;
 946             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 947             triggered = kTRUE;
 948             break;
 949         }
 950         if (!triggered) {
 951           delete [] pParent;
 952           return 0;
 953         }
 954     }
 955
 956
 957     // Check if there is a ccbar or bbbar pair with at least one of the two
 958     // in fYMin < y < fYMax
 959     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 960       TParticle *hvq;
 961       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 962       Float_t yQ;
 963       Int_t   pdgQ;
 964       for(i=0; i<np; i++) {
 965         hvq = (TParticle*)fParticles->At(i);
 966         pdgQ = hvq->GetPdgCode();
 967         if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
 968         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 969         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
 970                             (hvq->Energy()-hvq->Pz()+1.e-13));
 971         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
 972       }
 973       if (!theQ || !theQbar || !inYcut) {
 974         delete[] pParent;
 975         return 0;
 976       }
 977     }
 978
 979     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 980     if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
 981          && (fCutOnChild == 1) ) {
 982       if ( !CheckKinematicsOnChild() ) {
 983         delete[] pParent;
 984         return 0;
 985       }
 986     }
 987
 988
 989     for (i = 0; i < np; i++) {
 990         Int_t trackIt = 0;
 991         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 992         kf = CheckPDGCode(iparticle->GetPdgCode());
 993         Int_t ks = iparticle->GetStatusCode();
 994         Int_t km = iparticle->GetFirstMother();
 995         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 996             (ks != 1) ||
 997             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 998             nc++;
 999             if (ks == 1) trackIt = 1;
1000             Int_t ipa = iparticle->GetFirstMother()-1;
1001
1002             iparent = (ipa > -1) ? pParent[ipa] : -1;
1003
1004 //
1005 // store track information
1006             p[0] = iparticle->Px();
1007             p[1] = iparticle->Py();
1008             p[2] = iparticle->Pz();
1009             p[3] = iparticle->Energy();
1010
1011
1012             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1013             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1014             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1015
1016             Float_t tof = fEventTime + kconv * iparticle->T();
1017
1018             PushTrack(fTrackIt*trackIt, iparent, kf,
1019                       p[0], p[1], p[2], p[3],
1020                       origin[0], origin[1], origin[2], tof,
1021                       polar[0], polar[1], polar[2],
1022                       kPPrimary, nt, 1., ks);
1023             //
1024             // Special Treatment to store color-flow
1025             //
1026             if (ks == 3 || ks == 13 || ks == 14) {
1027                 TParticle* particle = 0;
1028                 if (fStack) {
1029                     particle = fStack->Particle(nt);
1030                 } else {
1031                     particle = gAlice->Stack()->Particle(nt);
1032                 }
1033                 particle->SetFirstDaughter(fPythia->GetK(2, i));
1034                 particle->SetLastDaughter(fPythia->GetK(3, i));
1035             }
1036
1037             KeepTrack(nt);
1038             pParent[i] = nt;
1039             SetHighWaterMark(nt);
1040
1041         } // select particle
1042     } // particle loop
1043
1044     delete[] pParent;
1045
1046     return 1;
1047 }
1048
1049
1050 void AliGenPythia::FinishRun()
1051 {
1052 // Print x-section summary
1053     fPythia->Pystat(1);
1054
1055     if (fNev > 0.) {
1056         fQ  /= fNev;
1057         fX1 /= fNev;
1058         fX2 /= fNev;
1059     }
1060
1061     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1062     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1063 }
1064
1065 void AliGenPythia::AdjustWeights() const
1066 {
1067 // Adjust the weights after generation of all events
1068 //
1069     if (gAlice) {
1070         TParticle *part;
1071         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1072         for (Int_t i=0; i<ntrack; i++) {
1073             part= gAlice->GetMCApp()->Particle(i);
1074             part->SetWeight(part->GetWeight()*fKineBias);
1075         }
1076     }
1077 }
1078
1079 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1080 {
1081 // Treat protons as inside nuclei with mass numbers a1 and a2
1082
1083     fAProjectile = a1;
1084     fATarget     = a2;
1085     fSetNuclei   = kTRUE;
1086 }
1087
1088
1089 void AliGenPythia::MakeHeader()
1090 {
1091 //
1092 // Make header for the simulated event
1093 //
1094   if (gAlice) {
1095     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1096         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1097   }
1098
1099 // Builds the event header, to be called after each event
1100     if (fHeader) delete fHeader;
1101     fHeader = new AliGenPythiaEventHeader("Pythia");
1102 //
1103 // Event type
1104     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1105 //
1106 // Number of trials
1107     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1108 //
1109 // Event Vertex
1110     fHeader->SetPrimaryVertex(fVertex);
1111 //
1112 // Jets that have triggered
1113
1114     if (fProcess == kPyJets)
1115     {
1116         Int_t ntrig, njet;
1117         Float_t jets[4][10];
1118         GetJets(njet, ntrig, jets);
1119
1120
1121         for (Int_t i = 0; i < ntrig; i++) {
1122             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1123                                                         jets[3][i]);
1124         }
1125     }
1126 //
1127 // Copy relevant information from external header, if present.
1128 //
1129     Float_t uqJet[4];
1130
1131     if (fRL) {
1132         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1133         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1134         {
1135             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1136
1137
1138             exHeader->TriggerJet(i, uqJet);
1139             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1140         }
1141     }
1142 //
1143 // Store quenching parameters
1144 //
1145     if (fQuench){
1146         Double_t z[4];
1147         Double_t xp, yp;
1148         if (fQuench == 1) {
1149             // Pythia::Quench()
1150             fPythia->GetQuenchingParameters(xp, yp, z);
1151         } else {
1152             // Pyquen
1153             Double_t r1 = PARIMP.rb1;
1154             Double_t r2 = PARIMP.rb2;
1155             Double_t b  = PARIMP.b1;
1156             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1157             Double_t phi = PARIMP.psib1;
1158             xp = r * TMath::Cos(phi);
1159             yp = r * TMath::Sin(phi);
1160
1161         }
1162             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1163             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1164         }
1165 //
1166 // Store pt^hard
1167     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1168 //
1169 //  Pass header
1170 //
1171     AddHeader(fHeader);
1172 }
1173
1174 void AliGenPythia::AddHeader(AliGenEventHeader* header)
1175 {
1176     // Add header to container or runloader
1177     if (fContainer) {
1178         fContainer->AddHeader(header);
1179     } else {
1180         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1181     }
1182 }
1183
1184
1185 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1186 {
1187 // Check the kinematic trigger condition
1188 //
1189     Double_t eta[2];
1190     eta[0] = jet1->Eta();
1191     eta[1] = jet2->Eta();
1192     Double_t phi[2];
1193     phi[0] = jet1->Phi();
1194     phi[1] = jet2->Phi();
1195     Int_t    pdg[2];
1196     pdg[0] = jet1->GetPdgCode();
1197     pdg[1] = jet2->GetPdgCode();
1198     Bool_t   triggered = kFALSE;
1199
1200     if (fProcess == kPyJets) {
1201         Int_t njets = 0;
1202         Int_t ntrig = 0;
1203         Float_t jets[4][10];
1204 //
1205 // Use Pythia clustering on parton level to determine jet axis
1206 //
1207         GetJets(njets, ntrig, jets);
1208
1209         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1210 //
1211     } else {
1212         Int_t ij = 0;
1213         Int_t ig = 1;
1214         if (pdg[0] == kGamma) {
1215             ij = 1;
1216             ig = 0;
1217         }
1218         //Check eta range first...
1219         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1220             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1221         {
1222             //Eta is okay, now check phi range
1223             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1224                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1225             {
1226                 triggered = kTRUE;
1227             }
1228         }
1229     }
1230     return triggered;
1231 }
1232
1233
1234
1235 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1236 //
1237 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1238 //
1239     Bool_t checking = kFALSE;
1240     Int_t j, kcode, ks, km;
1241     Int_t nPartAcc = 0; //number of particles in the acceptance range
1242     Int_t numberOfAcceptedParticles = 1;
1243     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1244     Int_t npart = fParticles->GetEntriesFast();
1245
1246     for (j = 0; j<npart; j++) {
1247         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1248         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1249         ks = jparticle->GetStatusCode();
1250         km = jparticle->GetFirstMother();
1251
1252         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1253             nPartAcc++;
1254         }
1255         if( numberOfAcceptedParticles <= nPartAcc){
1256           checking = kTRUE;
1257           break;
1258         }
1259     }
1260
1261     return checking;
1262 }
1263
1264
1265 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
1266 {
1267 // Assignment operator
1268     rhs.Copy(*this);
1269     return *this;
1270 }
1271
1272 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1273 {
1274 //
1275 // Load event into Pythia Common Block
1276 //
1277
1278     Int_t npart = stack -> GetNprimary();
1279     Int_t n0 = 0;
1280
1281     if (!flag) {
1282         (fPythia->GetPyjets())->N = npart;
1283     } else {
1284         n0 = (fPythia->GetPyjets())->N;
1285         (fPythia->GetPyjets())->N = n0 + npart;
1286     }
1287
1288
1289     for (Int_t part = 0; part < npart; part++) {
1290         TParticle *mPart = stack->Particle(part);
1291
1292         Int_t kf     =  mPart->GetPdgCode();
1293         Int_t ks     =  mPart->GetStatusCode();
1294         Int_t idf    =  mPart->GetFirstDaughter();
1295         Int_t idl    =  mPart->GetLastDaughter();
1296
1297         if (reHadr) {
1298             if (ks == 11 || ks == 12) {
1299                 ks  -= 10;
1300                 idf  = -1;
1301                 idl  = -1;
1302             }
1303         }
1304
1305         Float_t px = mPart->Px();
1306         Float_t py = mPart->Py();
1307         Float_t pz = mPart->Pz();
1308         Float_t e  = mPart->Energy();
1309         Float_t m  = mPart->GetCalcMass();
1310
1311
1312         (fPythia->GetPyjets())->P[0][part+n0] = px;
1313         (fPythia->GetPyjets())->P[1][part+n0] = py;
1314         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1315         (fPythia->GetPyjets())->P[3][part+n0] = e;
1316         (fPythia->GetPyjets())->P[4][part+n0] = m;
1317
1318         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1319         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1320         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1321         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1322         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1323     }
1324 }
1325
1326
1327 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1328 {
1329 //
1330 //  Calls the Pythia jet finding algorithm to find jets in the current event
1331 //
1332 //
1333 //
1334 //  Save jets
1335     Int_t n     = fPythia->GetN();
1336
1337 //
1338 //  Run Jet Finder
1339     fPythia->Pycell(njets);
1340     Int_t i;
1341     for (i = 0; i < njets; i++) {
1342         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1343         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1344         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1345         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1346
1347         jets[0][i] = px;
1348         jets[1][i] = py;
1349         jets[2][i] = pz;
1350         jets[3][i] = e;
1351     }
1352 }
1353
1354
1355
1356 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1357 {
1358 //
1359 //  Calls the Pythia clustering algorithm to find jets in the current event
1360 //
1361     Int_t n     = fPythia->GetN();
1362     nJets       = 0;
1363     nJetsTrig   = 0;
1364     if (fJetReconstruction == kCluster) {
1365 //
1366 //  Configure cluster algorithm
1367 //
1368         fPythia->SetPARU(43, 2.);
1369         fPythia->SetMSTU(41, 1);
1370 //
1371 //  Call cluster algorithm
1372 //
1373         fPythia->Pyclus(nJets);
1374 //
1375 //  Loading jets from common block
1376 //
1377     } else {
1378
1379 //
1380 //  Run Jet Finder
1381         fPythia->Pycell(nJets);
1382     }
1383
1384     Int_t i;
1385     for (i = 0; i < nJets; i++) {
1386         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1387         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1388         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1389         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1390         Float_t pt    = TMath::Sqrt(px * px + py * py);
1391         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1392         Float_t theta = TMath::ATan2(pt,pz);
1393         Float_t et    = e * TMath::Sin(theta);
1394         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1395         if (
1396             eta > fEtaMinJet && eta < fEtaMaxJet &&
1397             phi > fPhiMinJet && phi < fPhiMaxJet &&
1398             et  > fEtMinJet  && et  < fEtMaxJet
1399             )
1400         {
1401             jets[0][nJetsTrig] = px;
1402             jets[1][nJetsTrig] = py;
1403             jets[2][nJetsTrig] = pz;
1404             jets[3][nJetsTrig] = e;
1405             nJetsTrig++;
1406 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1407         } else {
1408 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1409         }
1410     }
1411 }
1412
1413 void AliGenPythia::GetSubEventTime()
1414 {
1415   // Calculates time of the next subevent
1416   fEventTime = 0.;
1417   if (fEventsTime) {
1418     TArrayF &array = *fEventsTime;
1419     fEventTime = array[fCurSubEvent++];
1420   }
1421   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1422   return;
1423 }
1424
1425 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1426 {
1427   // Is particle in EMCAL acceptance?
1428   // phi in degrees, etamin=-etamax
1429   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1430      eta < fEMCALEta  )
1431     return kTRUE;
1432   else
1433     return kFALSE;
1434 }
1435
1436 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1437 {
1438   // Is particle in PHOS acceptance?
1439   // Acceptance slightly larger considered.
1440   // phi in degrees, etamin=-etamax
1441   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1442      eta < fPHOSEta  )
1443     return kTRUE;
1444   else
1445     return kFALSE;
1446 }
1447
1448
1449
1450 #ifdef never
1451 void AliGenPythia::Streamer(TBuffer &R__b)
1452 {
1453    // Stream an object of class AliGenPythia.
1454
1455    if (R__b.IsReading()) {
1456       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1457       AliGenerator::Streamer(R__b);
1458       R__b >> (Int_t&)fProcess;
1459       R__b >> (Int_t&)fStrucFunc;
1460       R__b >> (Int_t&)fForceDecay;
1461       R__b >> fEnergyCMS;
1462       R__b >> fKineBias;
1463       R__b >> fTrials;
1464       fParentSelect.Streamer(R__b);
1465       fChildSelect.Streamer(R__b);
1466       R__b >> fXsection;
1467 //      (AliPythia::Instance())->Streamer(R__b);
1468       R__b >> fPtHardMin;
1469       R__b >> fPtHardMax;
1470 //      if (fDecayer) fDecayer->Streamer(R__b);
1471    } else {
1472       R__b.WriteVersion(AliGenPythia::IsA());
1473       AliGenerator::Streamer(R__b);
1474       R__b << (Int_t)fProcess;
1475       R__b << (Int_t)fStrucFunc;
1476       R__b << (Int_t)fForceDecay;
1477       R__b << fEnergyCMS;
1478       R__b << fKineBias;
1479       R__b << fTrials;
1480       fParentSelect.Streamer(R__b);
1481       fChildSelect.Streamer(R__b);
1482       R__b << fXsection;
1483 //      R__b << fPythia;
1484       R__b << fPtHardMin;
1485       R__b << fPtHardMax;
1486       //     fDecayer->Streamer(R__b);
1487    }
1488 }
1489 #endif
1490
1491
1492