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