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