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