PYTHIA6/AliGenPythia.cxx

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