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