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