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