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