PYTHIA6/AliGenPythiaPlus.cxx

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