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