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