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