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