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 //  Forward Paramters to the AliPythia object
 316     fDecayer->SetForceDecay(fForceDecay);
 317     fDecayer->Init();
 318
 319
 320     fPythia->SetCKIN(3,fPtHardMin);
 321     fPythia->SetCKIN(4,fPtHardMax);
 322     fPythia->SetCKIN(7,fYHardMin);
 323     fPythia->SetCKIN(8,fYHardMax);
 324
 325     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
 326     // Fragmentation?
 327     if (fFragmentation) {
 328       fPythia->SetMSTP(111,1);
 329     } else {
 330       fPythia->SetMSTP(111,0);
 331     }
 332
 333
 334 //  initial state radiation
 335     fPythia->SetMSTP(61,fGinit);
 336 //  final state radiation
 337     fPythia->SetMSTP(71,fGfinal);
 338 //  pt - kick
 339     if (fPtKick > 0.) {
 340         fPythia->SetMSTP(91,1);
 341         fPythia->SetPARP(91,fPtKick);
 342     } else {
 343         fPythia->SetMSTP(91,0);
 344     }
 345
 346
 347     if (fReadFromFile) {
 348         fRL  =  AliRunLoader::Open(fFileName, "Partons");
 349         fRL->LoadKinematics();
 350         fRL->LoadHeader();
 351     } else {
 352         fRL = 0x0;
 353     }
 354 // Switch off Heavy Flavors on request
 355     if (fHFoff) {
 356         fPythia->SetMSTP(58, 3);
 357         fPythia->SetMSTJ(45, 3);
 358         for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
 359     }
 360  //
 361     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 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
 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                 } //  PushTrack loop
 760             }
 761         } else {
 762             nc = GenerateMB();
 763         } // mb ?
 764
 765         GetSubEventTime();
 766
 767         delete[] pParent;
 768         delete[] pSelected;
 769         delete[] trackIt;
 770
 771         if (nc > 0) {
 772           switch (fCountMode) {
 773           case kCountAll:
 774             // printf(" Count all \n");
 775             jev += nc;
 776             break;
 777           case kCountParents:
 778             // printf(" Count parents \n");
 779             jev += nParents;
 780             break;
 781           case kCountTrackables:
 782             // printf(" Count trackable \n");
 783             jev += nTkbles;
 784             break;
 785           }
 786             if (jev >= fNpart || fNpart == -1) {
 787                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 788
 789                 fQ  += fPythia->GetVINT(51);
 790                 fX1 += fPythia->GetVINT(41);
 791                 fX2 += fPythia->GetVINT(42);
 792                 fTrialsRun += fTrials;
 793                 fNev++;
 794                 MakeHeader();
 795                 break;
 796             }
 797         }
 798     } // event loop
 799     SetHighWaterMark(nt);
 800 //  adjust weight due to kinematic selection
 801 //    AdjustWeights();
 802 //  get cross-section
 803     fXsection=fPythia->GetPARI(1);
 804 }
 805
 806 Int_t  AliGenPythia::GenerateMB()
 807 {
 808 //
 809 // Min Bias selection and other global selections
 810 //
 811     Int_t i, kf, nt, iparent;
 812     Int_t nc = 0;
 813     Float_t p[4];
 814     Float_t polar[3]   =   {0,0,0};
 815     Float_t origin[3]  =   {0,0,0};
 816 //  converts from mm/c to s
 817     const Float_t kconv=0.001/2.999792458e8;
 818
 819
 820
 821     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 822
 823
 824
 825     Int_t* pParent = new Int_t[np];
 826     for (i=0; i< np; i++) pParent[i] = -1;
 827     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 828         TParticle* jet1 = (TParticle *) fParticles->At(6);
 829         TParticle* jet2 = (TParticle *) fParticles->At(7);
 830         if (!CheckTrigger(jet1, jet2)) {
 831           delete [] pParent;
 832           return 0;
 833         }
 834     }
 835
 836     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
 837     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
 838
 839       Bool_t ok = kFALSE;
 840
 841       Int_t pdg  = 0;
 842       if (fFragPhotonInCalo) pdg = 22   ; // Photon
 843       else if (fPi0InCalo) pdg = 111 ; // Pi0
 844
 845       for (i=0; i< np; i++) {
 846         TParticle* iparticle = (TParticle *) fParticles->At(i);
 847         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
 848            iparticle->Pt() > fFragPhotonOrPi0MinPt){
 849           Int_t imother = iparticle->GetFirstMother() - 1;
 850           TParticle* pmother = (TParticle *) fParticles->At(imother);
 851           if(pdg == 111 ||
 852              (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
 853             {
 854               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 855               Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
 856               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
 857                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
 858                 ok =kTRUE;
 859             }
 860         }
 861       }
 862       if(!ok)
 863         return 0;
 864     }
 865
 866     if (fTriggerParticle) {
 867         Bool_t triggered = kFALSE;
 868         for (i = 0; i < np; i++) {
 869             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 870             kf = CheckPDGCode(iparticle->GetPdgCode());
 871             if (kf != fTriggerParticle) continue;
 872             if (iparticle->Pt() == 0.) continue;
 873             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 874             triggered = kTRUE;
 875             break;
 876         }
 877         if (!triggered) {
 878           delete [] pParent;
 879           return 0;
 880         }
 881     }
 882
 883
 884     // Check if there is a ccbar or bbbar pair with at least one of the two
 885     // in fYMin < y < fYMax
 886     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 887       TParticle *hvq;
 888       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 889       Float_t yQ;
 890       Int_t   pdgQ;
 891       for(i=0; i<np; i++) {
 892         hvq = (TParticle*)fParticles->At(i);
 893         pdgQ = hvq->GetPdgCode();
 894         if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
 895         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 896         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
 897                             (hvq->Energy()-hvq->Pz()+1.e-13));
 898         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
 899       }
 900       if (!theQ || !theQbar || !inYcut) {
 901         delete[] pParent;
 902         return 0;
 903       }
 904     }
 905
 906     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 907     if ( (fProcess == kPyW ||
 908           fProcess == kPyZ ||
 909           fProcess == kPyMbDefault ||
 910           fProcess == kPyMb ||
 911           fProcess == kPyMbNonDiffr)
 912          && (fCutOnChild == 1) ) {
 913       if ( !CheckKinematicsOnChild() ) {
 914         delete[] pParent;
 915         return 0;
 916       }
 917     }
 918
 919
 920     for (i = 0; i < np; i++) {
 921         Int_t trackIt = 0;
 922         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 923         kf = CheckPDGCode(iparticle->GetPdgCode());
 924         Int_t ks = iparticle->GetStatusCode();
 925         Int_t km = iparticle->GetFirstMother();
 926         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 927             (ks != 1) ||
 928             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 929             nc++;
 930             if (ks == 1) trackIt = 1;
 931             Int_t ipa = iparticle->GetFirstMother()-1;
 932
 933             iparent = (ipa > -1) ? pParent[ipa] : -1;
 934
 935 //
 936 // store track information
 937             p[0] = iparticle->Px();
 938             p[1] = iparticle->Py();
 939             p[2] = iparticle->Pz();
 940             p[3] = iparticle->Energy();
 941
 942
 943             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 944             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 945             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 946
 947             Float_t tof = fEventTime + kconv * iparticle->T();
 948
 949             PushTrack(fTrackIt*trackIt, iparent, kf,
 950                       p[0], p[1], p[2], p[3],
 951                       origin[0], origin[1], origin[2], tof,
 952                       polar[0], polar[1], polar[2],
 953                       kPPrimary, nt, 1., ks);
 954             //
 955             // Special Treatment to store color-flow
 956             //
 957             if (ks == 3 || ks == 13 || ks == 14) {
 958                 TParticle* particle = 0;
 959                 if (fStack) {
 960                     particle = fStack->Particle(nt);
 961                 } else {
 962                     particle = gAlice->Stack()->Particle(nt);
 963                 }
 964                 particle->SetFirstDaughter(fPythia->GetK(2, i));
 965                 particle->SetLastDaughter(fPythia->GetK(3, i));
 966             }
 967
 968             KeepTrack(nt);
 969             pParent[i] = nt;
 970             SetHighWaterMark(nt);
 971
 972         } // select particle
 973     } // particle loop
 974
 975     delete[] pParent;
 976
 977     return 1;
 978 }
 979
 980
 981 void AliGenPythia::FinishRun()
 982 {
 983 // Print x-section summary
 984     fPythia->Pystat(1);
 985
 986     if (fNev > 0.) {
 987         fQ  /= fNev;
 988         fX1 /= fNev;
 989         fX2 /= fNev;
 990     }
 991
 992     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 993     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 994 }
 995
 996 void AliGenPythia::AdjustWeights() const
 997 {
 998 // Adjust the weights after generation of all events
 999 //
1000     if (gAlice) {
1001         TParticle *part;
1002         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1003         for (Int_t i=0; i<ntrack; i++) {
1004             part= gAlice->GetMCApp()->Particle(i);
1005             part->SetWeight(part->GetWeight()*fKineBias);
1006         }
1007     }
1008 }
1009
1010 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1011 {
1012 // Treat protons as inside nuclei with mass numbers a1 and a2
1013
1014     fAProjectile = a1;
1015     fATarget     = a2;
1016     fSetNuclei   = kTRUE;
1017 }
1018
1019
1020 void AliGenPythia::MakeHeader()
1021 {
1022 //
1023 // Make header for the simulated event
1024 //
1025   if (gAlice) {
1026     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1027         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1028   }
1029
1030 // Builds the event header, to be called after each event
1031     if (fHeader) delete fHeader;
1032     fHeader = new AliGenPythiaEventHeader("Pythia");
1033 //
1034 // Event type
1035     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1036 //
1037 // Number of trials
1038     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1039 //
1040 // Event Vertex
1041     fHeader->SetPrimaryVertex(fVertex);
1042 //
1043 // Jets that have triggered
1044
1045     if (fProcess == kPyJets)
1046     {
1047         Int_t ntrig, njet;
1048         Float_t jets[4][10];
1049         GetJets(njet, ntrig, jets);
1050
1051
1052         for (Int_t i = 0; i < ntrig; i++) {
1053             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1054                                                         jets[3][i]);
1055         }
1056     }
1057 //
1058 // Copy relevant information from external header, if present.
1059 //
1060     Float_t uqJet[4];
1061
1062     if (fRL) {
1063         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1064         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1065         {
1066             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1067
1068
1069             exHeader->TriggerJet(i, uqJet);
1070             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1071         }
1072     }
1073 //
1074 // Store quenching parameters
1075 //
1076     if (fQuench){
1077         Double_t z[4];
1078         Double_t xp, yp;
1079         if (fQuench == 1) {
1080             // Pythia::Quench()
1081             fPythia->GetQuenchingParameters(xp, yp, z);
1082         } else {
1083             // Pyquen
1084             Double_t r1 = PARIMP.rb1;
1085             Double_t r2 = PARIMP.rb2;
1086             Double_t b  = PARIMP.b1;
1087             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1088             Double_t phi = PARIMP.psib1;
1089             xp = r * TMath::Cos(phi);
1090             yp = r * TMath::Sin(phi);
1091
1092         }
1093             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1094             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1095         }
1096 //
1097 // Store pt^hard
1098     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1099 //
1100 //  Pass header
1101 //
1102     AddHeader(fHeader);
1103 }
1104
1105 void AliGenPythia::AddHeader(AliGenEventHeader* header)
1106 {
1107     // Add header to container or runloader
1108     if (fContainer) {
1109         fContainer->AddHeader(header);
1110     } else {
1111         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1112     }
1113 }
1114
1115
1116 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1117 {
1118 // Check the kinematic trigger condition
1119 //
1120     Double_t eta[2];
1121     eta[0] = jet1->Eta();
1122     eta[1] = jet2->Eta();
1123     Double_t phi[2];
1124     phi[0] = jet1->Phi();
1125     phi[1] = jet2->Phi();
1126     Int_t    pdg[2];
1127     pdg[0] = jet1->GetPdgCode();
1128     pdg[1] = jet2->GetPdgCode();
1129     Bool_t   triggered = kFALSE;
1130
1131     if (fProcess == kPyJets) {
1132         Int_t njets = 0;
1133         Int_t ntrig = 0;
1134         Float_t jets[4][10];
1135 //
1136 // Use Pythia clustering on parton level to determine jet axis
1137 //
1138         GetJets(njets, ntrig, jets);
1139
1140         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1141 //
1142     } else {
1143         Int_t ij = 0;
1144         Int_t ig = 1;
1145         if (pdg[0] == kGamma) {
1146             ij = 1;
1147             ig = 0;
1148         }
1149         //Check eta range first...
1150         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1151             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1152         {
1153             //Eta is okay, now check phi range
1154             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1155                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1156             {
1157                 triggered = kTRUE;
1158             }
1159         }
1160     }
1161     return triggered;
1162 }
1163
1164
1165
1166 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1167 //
1168 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1169 //
1170     Bool_t checking = kFALSE;
1171     Int_t j, kcode, ks, km;
1172     Int_t nPartAcc = 0; //number of particles in the acceptance range
1173     Int_t numberOfAcceptedParticles = 1;
1174     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1175     Int_t npart = fParticles->GetEntriesFast();
1176
1177     for (j = 0; j<npart; j++) {
1178         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1179         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1180         ks = jparticle->GetStatusCode();
1181         km = jparticle->GetFirstMother();
1182
1183         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1184             nPartAcc++;
1185         }
1186         if( numberOfAcceptedParticles <= nPartAcc){
1187           checking = kTRUE;
1188           break;
1189         }
1190     }
1191
1192     return checking;
1193 }
1194
1195 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1196 {
1197 //
1198 // Load event into Pythia Common Block
1199 //
1200
1201     Int_t npart = stack -> GetNprimary();
1202     Int_t n0 = 0;
1203
1204     if (!flag) {
1205         (fPythia->GetPyjets())->N = npart;
1206     } else {
1207         n0 = (fPythia->GetPyjets())->N;
1208         (fPythia->GetPyjets())->N = n0 + npart;
1209     }
1210
1211
1212     for (Int_t part = 0; part < npart; part++) {
1213         TParticle *mPart = stack->Particle(part);
1214
1215         Int_t kf     =  mPart->GetPdgCode();
1216         Int_t ks     =  mPart->GetStatusCode();
1217         Int_t idf    =  mPart->GetFirstDaughter();
1218         Int_t idl    =  mPart->GetLastDaughter();
1219
1220         if (reHadr) {
1221             if (ks == 11 || ks == 12) {
1222                 ks  -= 10;
1223                 idf  = -1;
1224                 idl  = -1;
1225             }
1226         }
1227
1228         Float_t px = mPart->Px();
1229         Float_t py = mPart->Py();
1230         Float_t pz = mPart->Pz();
1231         Float_t e  = mPart->Energy();
1232         Float_t m  = mPart->GetCalcMass();
1233
1234
1235         (fPythia->GetPyjets())->P[0][part+n0] = px;
1236         (fPythia->GetPyjets())->P[1][part+n0] = py;
1237         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1238         (fPythia->GetPyjets())->P[3][part+n0] = e;
1239         (fPythia->GetPyjets())->P[4][part+n0] = m;
1240
1241         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1242         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1243         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1244         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1245         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1246     }
1247 }
1248
1249
1250 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1251 {
1252 //
1253 //  Calls the Pythia jet finding algorithm to find jets in the current event
1254 //
1255 //
1256 //
1257 //  Save jets
1258     Int_t n     = fPythia->GetN();
1259
1260 //
1261 //  Run Jet Finder
1262     fPythia->Pycell(njets);
1263     Int_t i;
1264     for (i = 0; i < njets; i++) {
1265         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1266         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1267         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1268         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1269
1270         jets[0][i] = px;
1271         jets[1][i] = py;
1272         jets[2][i] = pz;
1273         jets[3][i] = e;
1274     }
1275 }
1276
1277
1278
1279 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1280 {
1281 //
1282 //  Calls the Pythia clustering algorithm to find jets in the current event
1283 //
1284     Int_t n     = fPythia->GetN();
1285     nJets       = 0;
1286     nJetsTrig   = 0;
1287     if (fJetReconstruction == kCluster) {
1288 //
1289 //  Configure cluster algorithm
1290 //
1291         fPythia->SetPARU(43, 2.);
1292         fPythia->SetMSTU(41, 1);
1293 //
1294 //  Call cluster algorithm
1295 //
1296         fPythia->Pyclus(nJets);
1297 //
1298 //  Loading jets from common block
1299 //
1300     } else {
1301
1302 //
1303 //  Run Jet Finder
1304         fPythia->Pycell(nJets);
1305     }
1306
1307     Int_t i;
1308     for (i = 0; i < nJets; i++) {
1309         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1310         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1311         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1312         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1313         Float_t pt    = TMath::Sqrt(px * px + py * py);
1314         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1315         Float_t theta = TMath::ATan2(pt,pz);
1316         Float_t et    = e * TMath::Sin(theta);
1317         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1318         if (
1319             eta > fEtaMinJet && eta < fEtaMaxJet &&
1320             phi > fPhiMinJet && phi < fPhiMaxJet &&
1321             et  > fEtMinJet  && et  < fEtMaxJet
1322             )
1323         {
1324             jets[0][nJetsTrig] = px;
1325             jets[1][nJetsTrig] = py;
1326             jets[2][nJetsTrig] = pz;
1327             jets[3][nJetsTrig] = e;
1328             nJetsTrig++;
1329 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1330         } else {
1331 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1332         }
1333     }
1334 }
1335
1336 void AliGenPythia::GetSubEventTime()
1337 {
1338   // Calculates time of the next subevent
1339   fEventTime = 0.;
1340   if (fEventsTime) {
1341     TArrayF &array = *fEventsTime;
1342     fEventTime = array[fCurSubEvent++];
1343   }
1344   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1345   return;
1346 }
1347
1348 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1349 {
1350   // Is particle in EMCAL acceptance?
1351   // phi in degrees, etamin=-etamax
1352   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1353      eta < fEMCALEta  )
1354     return kTRUE;
1355   else
1356     return kFALSE;
1357 }
1358
1359 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1360 {
1361   // Is particle in PHOS acceptance?
1362   // Acceptance slightly larger considered.
1363   // phi in degrees, etamin=-etamax
1364   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1365      eta < fPHOSEta  )
1366     return kTRUE;
1367   else
1368     return kFALSE;
1369 }
1370
1371
1372
1373 #ifdef never
1374 void AliGenPythia::Streamer(TBuffer &R__b)
1375 {
1376    // Stream an object of class AliGenPythia.
1377
1378    if (R__b.IsReading()) {
1379       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1380       AliGenerator::Streamer(R__b);
1381       R__b >> (Int_t&)fProcess;
1382       R__b >> (Int_t&)fStrucFunc;
1383       R__b >> (Int_t&)fForceDecay;
1384       R__b >> fEnergyCMS;
1385       R__b >> fKineBias;
1386       R__b >> fTrials;
1387       fParentSelect.Streamer(R__b);
1388       fChildSelect.Streamer(R__b);
1389       R__b >> fXsection;
1390 //      (AliPythia::Instance())->Streamer(R__b);
1391       R__b >> fPtHardMin;
1392       R__b >> fPtHardMax;
1393 //      if (fDecayer) fDecayer->Streamer(R__b);
1394    } else {
1395       R__b.WriteVersion(AliGenPythia::IsA());
1396       AliGenerator::Streamer(R__b);
1397       R__b << (Int_t)fProcess;
1398       R__b << (Int_t)fStrucFunc;
1399       R__b << (Int_t)fForceDecay;
1400       R__b << fEnergyCMS;
1401       R__b << fKineBias;
1402       R__b << fTrials;
1403       fParentSelect.Streamer(R__b);
1404       fChildSelect.Streamer(R__b);
1405       R__b << fXsection;
1406 //      R__b << fPythia;
1407       R__b << fPtHardMin;
1408       R__b << fPtHardMax;
1409       //     fDecayer->Streamer(R__b);
1410    }
1411 }
1412 #endif
1413
1414
1415