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