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