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