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     fBremssInCalo(kFALSE),
 115     fPi0InCalo(kFALSE) ,
 116     fBremssPi0InEMCAL(kFALSE),
 117     fBremssPi0InPHOS(kFALSE),
 118     fBremssPi0MinPt(0),
 119     fPHOSMinPhi(220.),
 120     fPHOSMaxPhi(320.),
 121     fPHOSEta(0.12),
 122     fEMCALMinPhi(80.),
 123     fEMCALMaxPhi(190.),
 124     fEMCALEta(0.7)
 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      fBremssInCalo(kFALSE),
 199      fPi0InCalo(kFALSE) ,
 200      fBremssPi0InEMCAL(kFALSE),
 201      fBremssPi0InPHOS(kFALSE),
 202      fBremssPi0MinPt(0),
 203      fPHOSMinPhi(220.),
 204      fPHOSMaxPhi(320.),
 205      fPHOSEta(0.12),
 206      fEMCALMinPhi(80.),
 207      fEMCALMaxPhi(190.),
 208      fEMCALEta(0.7)
 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      fBremssInCalo(kFALSE),
 290      fPi0InCalo(kFALSE) ,
 291      fBremssPi0InEMCAL(kFALSE),
 292      fBremssPi0InPHOS(kFALSE),
 293      fBremssPi0MinPt(0),
 294      fPHOSMinPhi(220.),
 295      fPHOSMaxPhi(320.),
 296      fPHOSEta(0.12),
 297      fEMCALMinPhi(80.),
 298      fEMCALMaxPhi(190.),
 299      fEMCALEta(0.7)
 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 bremsstrahlung or pi0 going to PHOS or EMCAL
 912     if (fProcess == kPyJets && (fBremssInCalo || fPi0InCalo) ) {
 913
 914       Bool_t ok = kFALSE;
 915
 916       Int_t pdg  = 0;
 917       if (fBremssInCalo) 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() > fBremssPi0MinPt){
 924           Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 925           Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
 926           if((fBremssPi0InEMCAL && IsInEMCAL(phi,eta)) ||
 927               fBremssPi0InPHOS    && IsInPHOS(phi,eta))
 928             ok =kTRUE;
 929         }
 930       }
 931       if(!ok)
 932         return 0;
 933     }
 934
 935     if (fTriggerParticle) {
 936         Bool_t triggered = kFALSE;
 937         for (i = 0; i < np; i++) {
 938             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 939             kf = CheckPDGCode(iparticle->GetPdgCode());
 940             if (kf != fTriggerParticle) continue;
 941             if (iparticle->Pt() == 0.) continue;
 942             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 943             triggered = kTRUE;
 944             break;
 945         }
 946         if (!triggered) {
 947           delete [] pParent;
 948           return 0;
 949         }
 950     }
 951
 952
 953     // Check if there is a ccbar or bbbar pair with at least one of the two
 954     // in fYMin < y < fYMax
 955     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 956       TParticle *hvq;
 957       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 958       Float_t yQ;
 959       Int_t   pdgQ;
 960       for(i=0; i<np; i++) {
 961         hvq = (TParticle*)fParticles->At(i);
 962         pdgQ = hvq->GetPdgCode();
 963         if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
 964         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 965         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
 966                             (hvq->Energy()-hvq->Pz()+1.e-13));
 967         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
 968       }
 969       if (!theQ || !theQbar || !inYcut) {
 970         delete[] pParent;
 971         return 0;
 972       }
 973     }
 974
 975     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 976     if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
 977          && (fCutOnChild == 1) ) {
 978       if ( !CheckKinematicsOnChild() ) {
 979         delete[] pParent;
 980         return 0;
 981       }
 982     }
 983
 984
 985     for (i = 0; i < np; i++) {
 986         Int_t trackIt = 0;
 987         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 988         kf = CheckPDGCode(iparticle->GetPdgCode());
 989         Int_t ks = iparticle->GetStatusCode();
 990         Int_t km = iparticle->GetFirstMother();
 991         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 992             (ks != 1) ||
 993             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 994             nc++;
 995             if (ks == 1) trackIt = 1;
 996             Int_t ipa = iparticle->GetFirstMother()-1;
 997
 998             iparent = (ipa > -1) ? pParent[ipa] : -1;
 999
1000 //
1001 // store track information
1002             p[0] = iparticle->Px();
1003             p[1] = iparticle->Py();
1004             p[2] = iparticle->Pz();
1005             p[3] = iparticle->Energy();
1006
1007
1008             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1009             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1010             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1011
1012             Float_t tof = fEventTime + kconv * iparticle->T();
1013
1014             PushTrack(fTrackIt*trackIt, iparent, kf,
1015                       p[0], p[1], p[2], p[3],
1016                       origin[0], origin[1], origin[2], tof,
1017                       polar[0], polar[1], polar[2],
1018                       kPPrimary, nt, 1., ks);
1019             //
1020             // Special Treatment to store color-flow
1021             //
1022             if (ks == 3 || ks == 13 || ks == 14) {
1023                 TParticle* particle = 0;
1024                 if (fStack) {
1025                     particle = fStack->Particle(nt);
1026                 } else {
1027                     particle = gAlice->Stack()->Particle(nt);
1028                 }
1029                 particle->SetFirstDaughter(fPythia->GetK(2, i));
1030                 particle->SetLastDaughter(fPythia->GetK(3, i));
1031             }
1032
1033             KeepTrack(nt);
1034             pParent[i] = nt;
1035             SetHighWaterMark(nt);
1036
1037         } // select particle
1038     } // particle loop
1039
1040     delete[] pParent;
1041
1042     return 1;
1043 }
1044
1045
1046 void AliGenPythia::FinishRun()
1047 {
1048 // Print x-section summary
1049     fPythia->Pystat(1);
1050
1051     if (fNev > 0.) {
1052         fQ  /= fNev;
1053         fX1 /= fNev;
1054         fX2 /= fNev;
1055     }
1056
1057     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1058     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1059 }
1060
1061 void AliGenPythia::AdjustWeights() const
1062 {
1063 // Adjust the weights after generation of all events
1064 //
1065     if (gAlice) {
1066         TParticle *part;
1067         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1068         for (Int_t i=0; i<ntrack; i++) {
1069             part= gAlice->GetMCApp()->Particle(i);
1070             part->SetWeight(part->GetWeight()*fKineBias);
1071         }
1072     }
1073 }
1074
1075 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1076 {
1077 // Treat protons as inside nuclei with mass numbers a1 and a2
1078
1079     fAProjectile = a1;
1080     fATarget     = a2;
1081     fSetNuclei   = kTRUE;
1082 }
1083
1084
1085 void AliGenPythia::MakeHeader()
1086 {
1087 //
1088 // Make header for the simulated event
1089 //
1090   if (gAlice) {
1091     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1092         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1093   }
1094
1095 // Builds the event header, to be called after each event
1096     if (fHeader) delete fHeader;
1097     fHeader = new AliGenPythiaEventHeader("Pythia");
1098 //
1099 // Event type
1100     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1101 //
1102 // Number of trials
1103     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1104 //
1105 // Event Vertex
1106     fHeader->SetPrimaryVertex(fVertex);
1107 //
1108 // Jets that have triggered
1109
1110     if (fProcess == kPyJets)
1111     {
1112         Int_t ntrig, njet;
1113         Float_t jets[4][10];
1114         GetJets(njet, ntrig, jets);
1115
1116
1117         for (Int_t i = 0; i < ntrig; i++) {
1118             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1119                                                         jets[3][i]);
1120         }
1121     }
1122 //
1123 // Copy relevant information from external header, if present.
1124 //
1125     Float_t uqJet[4];
1126
1127     if (fRL) {
1128         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1129         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1130         {
1131             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1132
1133
1134             exHeader->TriggerJet(i, uqJet);
1135             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1136         }
1137     }
1138 //
1139 // Store quenching parameters
1140 //
1141     if (fQuench){
1142         Double_t z[4];
1143         Double_t xp, yp;
1144         if (fQuench == 1) {
1145             // Pythia::Quench()
1146             fPythia->GetQuenchingParameters(xp, yp, z);
1147         } else {
1148             // Pyquen
1149             Double_t r1 = PARIMP.rb1;
1150             Double_t r2 = PARIMP.rb2;
1151             Double_t b  = PARIMP.b1;
1152             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1153             Double_t phi = PARIMP.psib1;
1154             xp = r * TMath::Cos(phi);
1155             yp = r * TMath::Sin(phi);
1156
1157         }
1158             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1159             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1160         }
1161 //
1162 // Store pt^hard
1163     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1164 //
1165 //  Pass header
1166 //
1167     AddHeader(fHeader);
1168 }
1169
1170 void AliGenPythia::AddHeader(AliGenEventHeader* header)
1171 {
1172     // Add header to container or runloader
1173     if (fContainer) {
1174         fContainer->AddHeader(header);
1175     } else {
1176         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1177     }
1178 }
1179
1180
1181 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1182 {
1183 // Check the kinematic trigger condition
1184 //
1185     Double_t eta[2];
1186     eta[0] = jet1->Eta();
1187     eta[1] = jet2->Eta();
1188     Double_t phi[2];
1189     phi[0] = jet1->Phi();
1190     phi[1] = jet2->Phi();
1191     Int_t    pdg[2];
1192     pdg[0] = jet1->GetPdgCode();
1193     pdg[1] = jet2->GetPdgCode();
1194     Bool_t   triggered = kFALSE;
1195
1196     if (fProcess == kPyJets) {
1197         Int_t njets = 0;
1198         Int_t ntrig = 0;
1199         Float_t jets[4][10];
1200 //
1201 // Use Pythia clustering on parton level to determine jet axis
1202 //
1203         GetJets(njets, ntrig, jets);
1204
1205         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1206 //
1207     } else {
1208         Int_t ij = 0;
1209         Int_t ig = 1;
1210         if (pdg[0] == kGamma) {
1211             ij = 1;
1212             ig = 0;
1213         }
1214         //Check eta range first...
1215         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1216             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1217         {
1218             //Eta is okay, now check phi range
1219             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1220                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1221             {
1222                 triggered = kTRUE;
1223             }
1224         }
1225     }
1226     return triggered;
1227 }
1228
1229
1230
1231 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1232 //
1233 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1234 //
1235     Bool_t checking = kFALSE;
1236     Int_t j, kcode, ks, km;
1237     Int_t nPartAcc = 0; //number of particles in the acceptance range
1238     Int_t numberOfAcceptedParticles = 1;
1239     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1240     Int_t npart = fParticles->GetEntriesFast();
1241
1242     for (j = 0; j<npart; j++) {
1243         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1244         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1245         ks = jparticle->GetStatusCode();
1246         km = jparticle->GetFirstMother();
1247
1248         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1249             nPartAcc++;
1250         }
1251         if( numberOfAcceptedParticles <= nPartAcc){
1252           checking = kTRUE;
1253           break;
1254         }
1255     }
1256
1257     return checking;
1258 }
1259
1260
1261 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
1262 {
1263 // Assignment operator
1264     rhs.Copy(*this);
1265     return *this;
1266 }
1267
1268 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1269 {
1270 //
1271 // Load event into Pythia Common Block
1272 //
1273
1274     Int_t npart = stack -> GetNprimary();
1275     Int_t n0 = 0;
1276
1277     if (!flag) {
1278         (fPythia->GetPyjets())->N = npart;
1279     } else {
1280         n0 = (fPythia->GetPyjets())->N;
1281         (fPythia->GetPyjets())->N = n0 + npart;
1282     }
1283
1284
1285     for (Int_t part = 0; part < npart; part++) {
1286         TParticle *mPart = stack->Particle(part);
1287
1288         Int_t kf     =  mPart->GetPdgCode();
1289         Int_t ks     =  mPart->GetStatusCode();
1290         Int_t idf    =  mPart->GetFirstDaughter();
1291         Int_t idl    =  mPart->GetLastDaughter();
1292
1293         if (reHadr) {
1294             if (ks == 11 || ks == 12) {
1295                 ks  -= 10;
1296                 idf  = -1;
1297                 idl  = -1;
1298             }
1299         }
1300
1301         Float_t px = mPart->Px();
1302         Float_t py = mPart->Py();
1303         Float_t pz = mPart->Pz();
1304         Float_t e  = mPart->Energy();
1305         Float_t m  = mPart->GetCalcMass();
1306
1307
1308         (fPythia->GetPyjets())->P[0][part+n0] = px;
1309         (fPythia->GetPyjets())->P[1][part+n0] = py;
1310         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1311         (fPythia->GetPyjets())->P[3][part+n0] = e;
1312         (fPythia->GetPyjets())->P[4][part+n0] = m;
1313
1314         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1315         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1316         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1317         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1318         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1319     }
1320 }
1321
1322
1323 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1324 {
1325 //
1326 //  Calls the Pythia jet finding algorithm to find jets in the current event
1327 //
1328 //
1329 //
1330 //  Save jets
1331     Int_t n     = fPythia->GetN();
1332
1333 //
1334 //  Run Jet Finder
1335     fPythia->Pycell(njets);
1336     Int_t i;
1337     for (i = 0; i < njets; i++) {
1338         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1339         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1340         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1341         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1342
1343         jets[0][i] = px;
1344         jets[1][i] = py;
1345         jets[2][i] = pz;
1346         jets[3][i] = e;
1347     }
1348 }
1349
1350
1351
1352 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1353 {
1354 //
1355 //  Calls the Pythia clustering algorithm to find jets in the current event
1356 //
1357     Int_t n     = fPythia->GetN();
1358     nJets       = 0;
1359     nJetsTrig   = 0;
1360     if (fJetReconstruction == kCluster) {
1361 //
1362 //  Configure cluster algorithm
1363 //
1364         fPythia->SetPARU(43, 2.);
1365         fPythia->SetMSTU(41, 1);
1366 //
1367 //  Call cluster algorithm
1368 //
1369         fPythia->Pyclus(nJets);
1370 //
1371 //  Loading jets from common block
1372 //
1373     } else {
1374
1375 //
1376 //  Run Jet Finder
1377         fPythia->Pycell(nJets);
1378     }
1379
1380     Int_t i;
1381     for (i = 0; i < nJets; i++) {
1382         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1383         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1384         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1385         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1386         Float_t pt    = TMath::Sqrt(px * px + py * py);
1387         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1388         Float_t theta = TMath::ATan2(pt,pz);
1389         Float_t et    = e * TMath::Sin(theta);
1390         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1391         if (
1392             eta > fEtaMinJet && eta < fEtaMaxJet &&
1393             phi > fPhiMinJet && phi < fPhiMaxJet &&
1394             et  > fEtMinJet  && et  < fEtMaxJet
1395             )
1396         {
1397             jets[0][nJetsTrig] = px;
1398             jets[1][nJetsTrig] = py;
1399             jets[2][nJetsTrig] = pz;
1400             jets[3][nJetsTrig] = e;
1401             nJetsTrig++;
1402 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1403         } else {
1404 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1405         }
1406     }
1407 }
1408
1409 void AliGenPythia::GetSubEventTime()
1410 {
1411   // Calculates time of the next subevent
1412   fEventTime = 0.;
1413   if (fEventsTime) {
1414     TArrayF &array = *fEventsTime;
1415     fEventTime = array[fCurSubEvent++];
1416   }
1417   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1418   return;
1419 }
1420
1421 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1422 {
1423   // Is particle in EMCAL acceptance?
1424   // Acceptance slightly larger considered.
1425   // phi in degrees, etamin=-etamax
1426   if(phi > fEMCALMinPhi-0.1  && phi < fEMCALMaxPhi+0.1 &&
1427      eta < fEMCALEta+0.01   )
1428     return kTRUE;
1429   else
1430     return kFALSE;
1431 }
1432
1433 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1434 {
1435   // Is particle in PHOS acceptance?
1436   // Acceptance slightly larger considered.
1437   // phi in degrees, etamin=-etamax
1438   if(phi > fPHOSMinPhi-0.1  && phi < fPHOSMaxPhi+0.1 &&
1439      eta < fPHOSEta+0.01   )
1440     return kTRUE;
1441   else
1442     return kFALSE;
1443 }
1444
1445
1446
1447 #ifdef never
1448 void AliGenPythia::Streamer(TBuffer &R__b)
1449 {
1450    // Stream an object of class AliGenPythia.
1451
1452    if (R__b.IsReading()) {
1453       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1454       AliGenerator::Streamer(R__b);
1455       R__b >> (Int_t&)fProcess;
1456       R__b >> (Int_t&)fStrucFunc;
1457       R__b >> (Int_t&)fForceDecay;
1458       R__b >> fEnergyCMS;
1459       R__b >> fKineBias;
1460       R__b >> fTrials;
1461       fParentSelect.Streamer(R__b);
1462       fChildSelect.Streamer(R__b);
1463       R__b >> fXsection;
1464 //      (AliPythia::Instance())->Streamer(R__b);
1465       R__b >> fPtHardMin;
1466       R__b >> fPtHardMax;
1467 //      if (fDecayer) fDecayer->Streamer(R__b);
1468    } else {
1469       R__b.WriteVersion(AliGenPythia::IsA());
1470       AliGenerator::Streamer(R__b);
1471       R__b << (Int_t)fProcess;
1472       R__b << (Int_t)fStrucFunc;
1473       R__b << (Int_t)fForceDecay;
1474       R__b << fEnergyCMS;
1475       R__b << fKineBias;
1476       R__b << fTrials;
1477       fParentSelect.Streamer(R__b);
1478       fChildSelect.Streamer(R__b);
1479       R__b << fXsection;
1480 //      R__b << fPythia;
1481       R__b << fPtHardMin;
1482       R__b << fPtHardMax;
1483       //     fDecayer->Streamer(R__b);
1484    }
1485 }
1486 #endif
1487
1488
1489