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