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