Base class for common usage if Pythia6 and Pythia8.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 Mar 2008 20:04:17 +0000 (20:04 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 Mar 2008 20:04:17 +0000 (20:04 +0000)
Generator for Pythia6 and Pythia8 event generation.

PYTHIA6/AliGenPythiaPlus.cxx [new file with mode: 0644]
PYTHIA6/AliGenPythiaPlus.h [new file with mode: 0644]
PYTHIA6/AliPythia6.cxx [new file with mode: 0644]
PYTHIA6/AliPythia6.h [new file with mode: 0644]
PYTHIA6/AliPythia6LinkDef.h
PYTHIA6/AliPythiaBase.cxx [new file with mode: 0644]
PYTHIA6/AliPythiaBase.h [new file with mode: 0644]
PYTHIA6/libAliPythia6.pkg

diff --git a/PYTHIA6/AliGenPythiaPlus.cxx b/PYTHIA6/AliGenPythiaPlus.cxx
new file mode 100644 (file)
index 0000000..35ba99a
--- /dev/null
@@ -0,0 +1,1422 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//
+// Generator using the TPythia interface (via AliPythia)
+// to generate pp collisions.
+// Using SetNuclei() also nuclear modifications to the structure functions
+// can be taken into account. This makes, of course, only sense for the
+// generation of the products of hard processes (heavy flavor, jets ...)
+//
+// andreas.morsch@cern.ch
+//
+
+#include <TClonesArray.h>
+#include <TDatabasePDG.h>
+#include <TParticle.h>
+#include <TPDGCode.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include "AliConst.h"
+#include "AliDecayerPythia.h"
+#include "AliGenPythiaPlus.h"
+#include "AliHeader.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliPythiaBase.h"
+#include "AliPythiaRndm.h"
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliRunLoader.h"
+#include "AliMC.h"
+#include "PyquenCommon.h"
+
+ClassImp(AliGenPythiaPlus)
+
+
+AliGenPythiaPlus::AliGenPythiaPlus():
+    AliGenMC(),
+    fPythia(0),
+    fProcess(kPyCharm),          
+    fStrucFunc(kCTEQ5L), 
+    fEnergyCMS(5500.),
+    fKineBias(0.),
+    fTrials(0),
+    fTrialsRun(0),
+    fQ(0.),
+    fX1(0.),
+    fX2(0.),
+    fEventTime(0.),
+    fInteractionRate(0.),
+    fTimeWindow(0.),
+    fCurSubEvent(0),
+    fEventsTime(0),
+    fNev(0),
+    fFlavorSelect(0),
+    fXsection(0.),
+    fPtHardMin(0.),
+    fPtHardMax(1.e4),
+    fYHardMin(-1.e10),
+    fYHardMax(1.e10),
+    fGinit(1),
+    fGfinal(1),
+    fHadronisation(1),
+    fNpartons(0),
+    fReadFromFile(0),
+    fQuench(0),
+    fPtKick(1.),
+    fFullEvent(kTRUE),
+    fDecayer(new AliDecayerPythia()),
+    fDebugEventFirst(-1),
+    fDebugEventLast(-1),
+    fEtMinJet(0.),      
+    fEtMaxJet(1.e4),      
+    fEtaMinJet(-20.),     
+    fEtaMaxJet(20.),     
+    fPhiMinJet(0.),     
+    fPhiMaxJet(2.* TMath::Pi()),     
+    fJetReconstruction(kCell),
+    fEtaMinGamma(-20.),      
+    fEtaMaxGamma(20.),      
+    fPhiMinGamma(0.),      
+    fPhiMaxGamma(2. * TMath::Pi()),      
+    fPycellEtaMax(2.),     
+    fPycellNEta(274),       
+    fPycellNPhi(432),       
+    fPycellThreshold(0.),  
+    fPycellEtSeed(4.),     
+    fPycellMinEtJet(10.),  
+    fPycellMaxRadius(1.), 
+    fStackFillOpt(kFlavorSelection),   
+    fFeedDownOpt(kTRUE),    
+    fFragmentation(kTRUE),
+    fSetNuclei(kFALSE),
+    fNewMIS(kFALSE),   
+    fHFoff(kFALSE),    
+    fTriggerParticle(0),
+    fTriggerEta(0.9),     
+    fCountMode(kCountAll),      
+    fHeader(0),  
+    fRL(0),      
+    fFileName(0),
+    fFragPhotonInCalo(kFALSE),
+    fPi0InCalo(kFALSE) ,
+    fPhotonInCalo(kFALSE),
+    fCheckEMCAL(kFALSE),
+    fCheckPHOS(kFALSE),
+    fCheckPHOSeta(kFALSE),
+    fFragPhotonOrPi0MinPt(0), 
+    fPhotonMinPt(0), 
+    fPHOSMinPhi(219.),
+    fPHOSMaxPhi(321.),
+    fPHOSEta(0.13),
+    fEMCALMinPhi(79.),
+    fEMCALMaxPhi(191.),
+    fEMCALEta(0.71)
+
+{
+// Default Constructor
+  SetNuclei(0,0);
+  if (!AliPythiaRndm::GetPythiaRandom()) 
+      AliPythiaRndm::SetPythiaRandom(GetRandom());
+}
+
+AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia)
+    :AliGenMC(-1),
+     fPythia(pythia),
+     fProcess(kPyCharm),          
+     fStrucFunc(kCTEQ5L), 
+     fEnergyCMS(5500.),
+     fKineBias(0.),
+     fTrials(0),
+     fTrialsRun(0),
+     fQ(0.),
+     fX1(0.),
+     fX2(0.),
+     fEventTime(0.),
+     fInteractionRate(0.),
+     fTimeWindow(0.),
+     fCurSubEvent(0),
+     fEventsTime(0),
+     fNev(0),
+     fFlavorSelect(0),
+     fXsection(0.),
+     fPtHardMin(0.),
+     fPtHardMax(1.e4),
+     fYHardMin(-1.e10),
+     fYHardMax(1.e10),
+     fGinit(kTRUE),
+     fGfinal(kTRUE),
+     fHadronisation(kTRUE),
+     fNpartons(0),
+     fReadFromFile(kFALSE),
+     fQuench(kFALSE),
+     fPtKick(1.),
+     fFullEvent(kTRUE),
+     fDecayer(new AliDecayerPythia()),
+     fDebugEventFirst(-1),
+     fDebugEventLast(-1),
+     fEtMinJet(0.),      
+     fEtMaxJet(1.e4),      
+     fEtaMinJet(-20.),     
+     fEtaMaxJet(20.),     
+     fPhiMinJet(0.),     
+     fPhiMaxJet(2.* TMath::Pi()),     
+     fJetReconstruction(kCell),
+     fEtaMinGamma(-20.),      
+     fEtaMaxGamma(20.),      
+     fPhiMinGamma(0.),      
+     fPhiMaxGamma(2. * TMath::Pi()),      
+     fPycellEtaMax(2.),     
+     fPycellNEta(274),       
+     fPycellNPhi(432),       
+     fPycellThreshold(0.),  
+     fPycellEtSeed(4.),     
+     fPycellMinEtJet(10.),  
+     fPycellMaxRadius(1.), 
+     fStackFillOpt(kFlavorSelection),   
+     fFeedDownOpt(kTRUE),    
+     fFragmentation(kTRUE),
+     fSetNuclei(kFALSE),
+     fNewMIS(kFALSE),   
+     fHFoff(kFALSE),    
+     fTriggerParticle(0),
+     fTriggerEta(0.9),     
+     fCountMode(kCountAll),      
+     fHeader(0),  
+     fRL(0),      
+     fFileName(0),
+     fFragPhotonInCalo(kFALSE),
+     fPi0InCalo(kFALSE) ,
+     fPhotonInCalo(kFALSE),
+     fCheckEMCAL(kFALSE),
+     fCheckPHOS(kFALSE),
+     fCheckPHOSeta(kFALSE),
+     fFragPhotonOrPi0MinPt(0),
+     fPhotonMinPt(0),
+     fPHOSMinPhi(219.),
+     fPHOSMaxPhi(321.),
+     fPHOSEta(0.13),
+     fEMCALMinPhi(79.),
+     fEMCALMaxPhi(191.),
+     fEMCALEta(0.71)
+{
+// default charm production at 5. 5 TeV
+// semimuonic decay
+// structure function GRVHO
+//
+    fName = "Pythia";
+    fTitle= "Particle Generator using PYTHIA";
+    SetForceDecay();
+    // Set random number generator 
+    if (!AliPythiaRndm::GetPythiaRandom()) 
+      AliPythiaRndm::SetPythiaRandom(GetRandom());
+    fParticles = new TClonesArray("TParticle",1000);
+    SetNuclei(0,0);
+ }
+
+AliGenPythiaPlus::~AliGenPythiaPlus()
+{
+// Destructor
+  if(fEventsTime) delete fEventsTime;
+}
+
+void AliGenPythiaPlus::SetInteractionRate(Float_t rate,Float_t timewindow)
+{
+// Generate pileup using user specified rate
+    fInteractionRate = rate;
+    fTimeWindow = timewindow;
+    GeneratePileup();
+}
+
+void AliGenPythiaPlus::GeneratePileup()
+{
+// Generate sub events time for pileup
+    fEventsTime = 0;
+    if(fInteractionRate == 0.) {
+      Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
+      return;
+    }
+
+    Int_t npart = NumberParticles();
+    if(npart < 0) {
+      Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
+      return;
+    }
+
+    if(fEventsTime) delete fEventsTime;
+    fEventsTime = new TArrayF(npart);
+    TArrayF &array = *fEventsTime;
+    for(Int_t ipart = 0; ipart < npart; ipart++)
+      array[ipart] = 0.;
+
+    Float_t eventtime = 0.;
+    while(1)
+      {
+       eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+       if(eventtime > fTimeWindow) break;
+       array.Set(array.GetSize()+1);
+       array[array.GetSize()-1] = eventtime;
+      }
+
+    eventtime = 0.;
+    while(1)
+      {
+       eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+       if(TMath::Abs(eventtime) > fTimeWindow) break;
+       array.Set(array.GetSize()+1);
+       array[array.GetSize()-1] = eventtime;
+      }
+
+    SetNumberParticles(fEventsTime->GetSize());
+}
+
+void AliGenPythiaPlus::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+                                      Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
+{
+// Set pycell parameters
+    fPycellEtaMax    =  etamax;
+    fPycellNEta      =  neta;
+    fPycellNPhi      =  nphi;
+    fPycellThreshold =  thresh;
+    fPycellEtSeed    =  etseed;
+    fPycellMinEtJet  =  minet;
+    fPycellMaxRadius =  r;
+}
+
+
+
+void AliGenPythiaPlus::SetEventListRange(Int_t eventFirst, Int_t eventLast)
+{
+  // Set a range of event numbers, for which a table
+  // of generated particle will be printed
+  fDebugEventFirst = eventFirst;
+  fDebugEventLast  = eventLast;
+  if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
+}
+
+void AliGenPythiaPlus::Init()
+{
+// Initialisation
+    
+//    SetMC(AliPythia::Instance());
+//    fPythia=(AliPythia*) fMCEvGen;
+    
+//
+    fParentWeight=1./Float_t(fNpart);
+//
+
+    
+    fPythia->SetPtHardRange(fPtHardMin, fPtHardMax);
+    fPythia->SetYHardRange(fYHardMin, fYHardMax);
+    
+    if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);  
+    // Fragmentation?
+    if (fFragmentation) {
+      fPythia->SetFragmentation(1);
+    } else {
+      fPythia->SetFragmentation(0);
+    }
+
+
+//  initial state radiation   
+    fPythia->SetInitialAndFinalStateRadiation(fGinit, fGfinal);
+
+//  pt - kick
+    fPythia->SetIntrinsicKt(fPtKick);
+
+    if (fReadFromFile) {
+       fRL  =  AliRunLoader::Open(fFileName, "Partons");
+       fRL->LoadKinematics();
+       fRL->LoadHeader();
+    } else {
+       fRL = 0x0;
+    }
+ //
+    fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc);
+    //  Forward Paramters to the AliPythia object
+    fDecayer->SetForceDecay(fForceDecay);    
+// Switch off Heavy Flavors on request  
+    if (fHFoff) {
+       fPythia->SwitchHFOff();
+       // Switch off g->QQbar splitting in decay table
+       ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
+    }
+
+    fDecayer->Init();
+
+
+//  Parent and Children Selection
+    switch (fProcess) 
+    {
+    case kPyOldUEQ2ordered:
+    case kPyOldUEQ2ordered2:
+    case kPyOldPopcorn:
+      break;
+    case kPyCharm:
+    case kPyCharmUnforced:
+    case kPyCharmPbPbMNR:
+    case kPyCharmpPbMNR:
+    case kPyCharmppMNR:
+    case kPyCharmppMNRwmi:
+       fParentSelect[0] =   411;
+       fParentSelect[1] =   421;
+       fParentSelect[2] =   431;
+       fParentSelect[3] =  4122;
+       fFlavorSelect    =  4;  
+       break;
+    case kPyD0PbPbMNR:
+    case kPyD0pPbMNR:
+    case kPyD0ppMNR:
+       fParentSelect[0] =   421;
+       fFlavorSelect    =   4; 
+       break;
+    case kPyDPlusPbPbMNR:
+    case kPyDPluspPbMNR:
+    case kPyDPlusppMNR:
+       fParentSelect[0] =   411;
+       fFlavorSelect    =   4; 
+       break;
+    case kPyDPlusStrangePbPbMNR:
+    case kPyDPlusStrangepPbMNR:
+    case kPyDPlusStrangeppMNR:
+       fParentSelect[0] =   431;
+       fFlavorSelect    =   4; 
+       break;
+    case kPyBeauty:
+    case kPyBeautyPbPbMNR:
+    case kPyBeautypPbMNR:
+    case kPyBeautyppMNR:
+    case kPyBeautyppMNRwmi:
+       fParentSelect[0]=  511;
+       fParentSelect[1]=  521;
+       fParentSelect[2]=  531;
+       fParentSelect[3]= 5122;
+       fParentSelect[4]= 5132;
+       fParentSelect[5]= 5232;
+       fParentSelect[6]= 5332;
+       fFlavorSelect   = 5;    
+       break;
+    case kPyBeautyUnforced:
+       fParentSelect[0] =  511;
+       fParentSelect[1] =  521;
+       fParentSelect[2] =  531;
+       fParentSelect[3] = 5122;
+       fParentSelect[4] = 5132;
+       fParentSelect[5] = 5232;
+       fParentSelect[6] = 5332;
+       fFlavorSelect    = 5;   
+       break;
+    case kPyJpsiChi:
+    case kPyJpsi:
+       fParentSelect[0] = 443;
+       break;
+    case kPyMbDefault:
+    case kPyMb:
+    case kPyMbNonDiffr:
+    case kPyMbMSEL1:
+    case kPyJets:
+    case kPyDirectGamma:
+    case kPyLhwgMb:    
+       break;
+    case kPyW:
+    case kPyZ:
+        break;
+    }
+//
+//
+//  JetFinder for Trigger
+//
+//  Configure detector (EMCAL like)
+//
+    fPythia->SetPycellParameters(fPycellEtaMax,fPycellNEta, fPycellNPhi,
+                                fPycellThreshold, fPycellEtSeed, 
+                                fPycellMinEtJet, fPycellMaxRadius);
+//
+//  This counts the total number of calls to Pyevnt() per run.
+    fTrialsRun = 0;
+    fQ         = 0.;
+    fX1        = 0.;
+    fX2        = 0.;    
+    fNev       = 0 ;
+//    
+//
+//
+    AliGenMC::Init();
+//
+//
+//  
+    if (fSetNuclei) {
+       fDyBoost = 0;
+       Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
+    }
+    
+    if (fQuench) {
+       fPythia->InitQuenching(0., 0.1, 0.6e6, 0, 0.97, 30);
+    }
+
+//    fPythia->SetPARJ(200, 0.0);
+
+//    if (fQuench == 3) {
+//     // Nestor's change of the splittings
+//     fPythia->SetPARJ(200, 0.8);
+//     fPythia->SetMSTJ(41, 1);  // QCD radiation only
+//     fPythia->SetMSTJ(42, 2);  // angular ordering
+//     fPythia->SetMSTJ(44, 2);  // option to run alpha_s
+//     fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
+//     fPythia->SetMSTJ(50, 0);  // No coherence in first branching
+//     fPythia->SetPARJ(82, 1.); // Cut off for parton showers
+//    }
+}
+
+void AliGenPythiaPlus::Generate()
+{
+// Generate one event
+    
+    fDecayer->ForceDecay();
+
+    Float_t polar[3]   =   {0,0,0};
+    Float_t origin[3]  =   {0,0,0};
+    Float_t p[4];
+//  converts from mm/c to s
+    const Float_t kconv=0.001/2.999792458e8;
+//
+    Int_t nt=0;
+    Int_t jev=0;
+    Int_t j, kf;
+    fTrials=0;
+    fEventTime = 0.;
+    
+    
+
+    //  Set collision vertex position 
+    if (fVertexSmear == kPerEvent) Vertex();
+    
+//  event loop    
+    while(1)
+    {
+//
+// Produce event
+//
+//
+// Switch hadronisation off
+//
+//     fPythia->SwitchHadronisationOff();
+//
+// Either produce new event or read partons from file
+//     
+       if (!fReadFromFile) {
+           if (!fNewMIS) {
+               fPythia->GenerateEvent();
+           } else {
+               fPythia->GenerateMIEvent();
+           }
+           fNpartons = fPythia->GetNumberOfParticles();
+       } else {
+           printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
+           fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
+           fPythia->SetNumberOfParticles(0);
+           fPythia->LoadEvent(fRL->Stack(), 0 , 1);
+           fPythia->EditEventList(21);
+       }
+       
+//
+//  Run quenching routine 
+//
+       if (fQuench == 1) {
+           fPythia->Quench();
+       } else if (fQuench == 2){
+           fPythia->Pyquen(208., 0, 0.);
+       } else if (fQuench == 3) {
+           // Quenching is via multiplicative correction of the splittings
+       }
+       
+//
+// Switch hadronisation on
+//
+//     fPythia->SwitchHadronisationOn();
+//
+// .. and perform hadronisation
+//     printf("Calling hadronisation %d\n", fPythia->GetN());
+//     fPythia->HadronizeEvent();      
+       fTrials++;
+       fPythia->GetParticles(fParticles);
+       Boost();
+//
+//
+//
+       Int_t i;
+       
+       fNprimaries = 0;
+       Int_t np = fParticles->GetEntriesFast();
+       
+       if (np == 0) continue;
+//
+       
+//
+       Int_t* pParent   = new Int_t[np];
+       Int_t* pSelected = new Int_t[np];
+       Int_t* trackIt   = new Int_t[np];
+       for (i = 0; i < np; i++) {
+           pParent[i]   = -1;
+           pSelected[i] =  0;
+           trackIt[i]   =  0;
+       }
+
+       Int_t nc = 0;        // Total n. of selected particles
+       Int_t nParents = 0;  // Selected parents
+       Int_t nTkbles = 0;   // Trackable particles
+       if (fProcess != kPyMbDefault && 
+           fProcess != kPyMb && 
+           fProcess != kPyJets && 
+           fProcess != kPyDirectGamma &&
+           fProcess != kPyMbNonDiffr  &&
+           fProcess != kPyMbMSEL1     &&
+           fProcess != kPyW && 
+           fProcess != kPyZ &&
+           fProcess != kPyCharmppMNRwmi && 
+           fProcess != kPyBeautyppMNRwmi) {
+           
+           for (i = 0; i < np; i++) {
+               TParticle* iparticle = (TParticle *) fParticles->At(i);
+               Int_t ks = iparticle->GetStatusCode();
+               kf = CheckPDGCode(iparticle->GetPdgCode());
+// No initial state partons
+               if (ks==21) continue;
+//
+// Heavy Flavor Selection
+//
+               // quark ?
+               kf = TMath::Abs(kf);
+               Int_t kfl = kf;
+               // Resonance
+
+               if (kfl > 100000) kfl %= 100000;
+               if (kfl > 10000)  kfl %= 10000;
+               // meson ?
+               if  (kfl > 10) kfl/=100;
+               // baryon
+               if (kfl > 10) kfl/=10;
+               Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+               Int_t kfMo = 0;
+//
+// Establish mother daughter relation between heavy quarks and mesons
+//
+               if (kf >= fFlavorSelect && kf <= 6) {
+                   Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
+                   if (idau > -1) {
+                       TParticle* daughter = (TParticle *) fParticles->At(idau);
+                       Int_t pdgD = daughter->GetPdgCode();
+                       if (pdgD == 91 || pdgD == 92) {
+                           Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
+                           Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1)  : (daughter->GetLastDaughter());
+
+                           for (Int_t j = jmin; j <= jmax; j++)
+                               ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
+                       } // is string or cluster
+                   } // has daughter
+               } // heavy quark
+               
+
+               if (ipa > -1) {
+                   TParticle *  mother = (TParticle *) fParticles->At(ipa);
+                   kfMo = TMath::Abs(mother->GetPdgCode());
+               }
+               
+               // What to keep in Stack?
+               Bool_t flavorOK = kFALSE;
+               Bool_t selectOK = kFALSE;
+               if (fFeedDownOpt) {
+                   if (kfl >= fFlavorSelect) flavorOK = kTRUE;
+               } else {
+                   if (kfl > fFlavorSelect) {
+                       nc = -1;
+                       break;
+                   }
+                   if (kfl == fFlavorSelect) flavorOK = kTRUE;
+               }
+               switch (fStackFillOpt) {
+               case kFlavorSelection:
+                   selectOK = kTRUE;
+                   break;
+               case kParentSelection:
+                   if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
+                   break;
+               }
+               if (flavorOK && selectOK) { 
+//
+// Heavy flavor hadron or quark
+//
+// Kinematic seletion on final state heavy flavor mesons
+                   if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
+                   {
+                       continue;
+                   }
+                   pSelected[i] = 1;
+                   if (ParentSelected(kf)) ++nParents; // Update parent count
+//                 printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
+               } else {
+// Kinematic seletion on decay products
+                   if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
+                       && !KinematicSelection(iparticle, 1)) 
+                   {
+                       continue;
+                   }
+//
+// Decay products 
+// Select if mother was selected and is not tracked
+
+                   if (pSelected[ipa] && 
+                       !trackIt[ipa]  &&     // mother will be  tracked ?
+                       kfMo !=  5 &&         // mother is b-quark, don't store fragments          
+                       kfMo !=  4 &&         // mother is c-quark, don't store fragments 
+                       kf   != 92)           // don't store string
+                   {
+//
+// Semi-stable or de-selected: diselect decay products:
+// 
+//
+                       if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
+                       {
+                           Int_t ipF = iparticle->GetFirstDaughter();
+                           Int_t ipL = iparticle->GetLastDaughter();   
+                           if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
+                       }
+//                     printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
+                       pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
+                   }
+               }
+               if (pSelected[i] == -1) pSelected[i] = 0;
+               if (!pSelected[i]) continue;
+               // Count quarks only if you did not include fragmentation
+               if (fFragmentation && kf <= 10) continue;
+
+               nc++;
+// Decision on tracking
+               trackIt[i] = 0;
+//
+// Track final state particle
+               if (ks == 1) trackIt[i] = 1;
+// Track semi-stable particles
+               if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
+// Track particles selected by process if undecayed. 
+               if (fForceDecay == kNoDecay) {
+                   if (ParentSelected(kf)) trackIt[i] = 1;
+               } else {
+                   if (ParentSelected(kf)) trackIt[i] = 0;
+               }
+               if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
+//
+//
+
+           } // particle selection loop
+           if (nc > 0) {
+               for (i = 0; i < np; i++) {
+                   if (!pSelected[i]) continue;
+                   TParticle *  iparticle = (TParticle *) fParticles->At(i);
+                   kf = CheckPDGCode(iparticle->GetPdgCode());
+                   Int_t ks = iparticle->GetStatusCode();  
+                   p[0] = iparticle->Px();
+                   p[1] = iparticle->Py();
+                   p[2] = iparticle->Pz();
+                   p[3] = iparticle->Energy();
+                   
+                   origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
+                   origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
+                   origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
+                   
+                   Float_t tof   = kconv*iparticle->T();
+                   Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+                   Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
+                   PushTrack(fTrackIt*trackIt[i], iparent, kf, 
+                             p[0], p[1], p[2], p[3], 
+                             origin[0], origin[1], origin[2], tof, 
+                             polar[0], polar[1], polar[2],
+                             kPPrimary, nt, 1., ks);
+                   pParent[i] = nt;
+                   KeepTrack(nt);
+                   fNprimaries++;
+               } //  PushTrack loop
+           }
+       } else {
+           nc = GenerateMB();
+       } // mb ?
+       
+       GetSubEventTime();
+
+       delete[] pParent;
+       delete[] pSelected;
+       delete[] trackIt;
+
+       if (nc > 0) {
+         switch (fCountMode) {
+         case kCountAll:
+           // printf(" Count all \n");
+           jev += nc;
+           break;
+         case kCountParents:
+           // printf(" Count parents \n");
+           jev += nParents;
+           break;
+         case kCountTrackables:
+           // printf(" Count trackable \n");
+           jev += nTkbles;
+           break;
+         }
+           if (jev >= fNpart || fNpart == -1) {
+               fKineBias=Float_t(fNpart)/Float_t(fTrials);
+               fPythia->GetXandQ(fQ, fX1, fX2);
+               fTrialsRun += fTrials;
+               fNev++;
+               MakeHeader();
+               break;
+           }
+       }
+    } // event loop
+    SetHighWaterMark(nt);
+//  Adjust weight due to kinematic selection
+//  AdjustWeights();
+//  get cross-section
+    fXsection = fPythia->GetXSection();
+}
+
+Int_t  AliGenPythiaPlus::GenerateMB()
+{
+//
+// Min Bias selection and other global selections
+//
+    Int_t i, kf, nt, iparent;
+    Int_t nc = 0;
+    Float_t p[4];
+    Float_t polar[3]   =   {0,0,0};
+    Float_t origin[3]  =   {0,0,0};
+//  converts from mm/c to s
+    const Float_t kconv = 0.001 / 2.999792458e8;
+    
+    Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+    
+    Int_t* pParent = new Int_t[np];
+    for (i=0; i< np; i++) pParent[i] = -1;
+    if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+       TParticle* jet1 = (TParticle *) fParticles->At(6);
+       TParticle* jet2 = (TParticle *) fParticles->At(7);
+       if (!CheckTrigger(jet1, jet2)) {
+         delete [] pParent;
+         return 0;
+       }
+    }
+
+    // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
+    if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+
+      Bool_t ok = kFALSE;
+
+      Int_t pdg  = 0; 
+      if (fFragPhotonInCalo) pdg = 22   ; // Photon
+      else if (fPi0InCalo) pdg = 111 ; // Pi0
+
+      for (i=0; i< np; i++) {
+       TParticle* iparticle = (TParticle *) fParticles->At(i);
+       if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg && 
+          iparticle->Pt() > fFragPhotonOrPi0MinPt){
+           Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+         TParticle* pmother = (TParticle *) fParticles->At(imother);
+         if(pdg == 111 || 
+            (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
+           {
+             Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+             Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax         
+             if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
+                (fCheckPHOS    && IsInPHOS(phi,eta)) )
+               ok =kTRUE;
+           }
+       }
+      }
+      if(!ok)
+       return 0;
+    }
+    
+    
+     // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
+    if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+
+      Bool_t okd = kFALSE;
+
+      Int_t pdg  = 22; 
+      Int_t iphcand = -1;
+      for (i=0; i< np; i++) {
+        TParticle* iparticle = (TParticle *) fParticles->At(i);
+        Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+        Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax 
+        
+        if(iparticle->GetStatusCode() == 1 
+           && iparticle->GetPdgCode() == pdg   
+           && iparticle->Pt() > fPhotonMinPt    
+           && eta < fPHOSEta){                 
+           
+           // first check if the photon is in PHOS phi
+           if(IsInPHOS(phi,eta)){ 
+               okd = kTRUE;
+               break;
+           } 
+           if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
+            
+        }
+      }
+      
+      if(!okd && iphcand != -1) // execute rotation in phi 
+          RotatePhi(iphcand,okd);
+      
+      if(!okd)
+       return 0;
+    }
+    
+    if (fTriggerParticle) {
+       Bool_t triggered = kFALSE;
+       for (i = 0; i < np; i++) {
+           TParticle *  iparticle = (TParticle *) fParticles->At(i);
+           kf = CheckPDGCode(iparticle->GetPdgCode());
+           if (kf != fTriggerParticle) continue;
+           if (iparticle->Pt() == 0.) continue;
+           if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+           triggered = kTRUE;
+           break;
+       }
+       if (!triggered) {
+         delete [] pParent;
+         return 0;
+       }
+    }
+       
+
+    // Check if there is a ccbar or bbbar pair with at least one of the two
+    // in fYMin < y < fYMax
+    if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
+      TParticle *hvq;
+      Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
+      Float_t yQ;  
+      Int_t   pdgQ;
+      for(i=0; i<np; i++) {
+       hvq = (TParticle*)fParticles->At(i);
+       pdgQ = hvq->GetPdgCode();  
+       if(TMath::Abs(pdgQ) != fFlavorSelect) continue; 
+       if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
+       yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
+                           (hvq->Energy()-hvq->Pz()+1.e-13));
+       if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
+      }
+      if (!theQ || !theQbar || !inYcut) {
+       delete[] pParent;
+       return 0;
+      }
+    }
+
+    //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
+    if ( (fProcess == kPyW ||
+         fProcess == kPyZ ||
+         fProcess == kPyMbDefault ||
+         fProcess == kPyMb ||
+         fProcess == kPyMbNonDiffr)  
+        && (fCutOnChild == 1) ) {
+      if ( !CheckKinematicsOnChild() ) {
+       delete[] pParent;
+       return 0;
+      }
+    }
+  
+
+    for (i = 0; i < np; i++) {
+       Int_t trackIt = 0;
+       TParticle *  iparticle = (TParticle *) fParticles->At(i);
+       kf = CheckPDGCode(iparticle->GetPdgCode());
+       Int_t ks = iparticle->GetStatusCode();
+       Int_t km = iparticle->GetFirstMother();
+       if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
+           (ks != 1) ||
+           (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+           nc++;
+           if (ks == 1) trackIt = 1;
+
+           Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
+           iparent = (ipa > -1) ? pParent[ipa] : -1;
+           if (ipa >= np) fPythia->EventListing();
+           
+//
+// store track information
+           p[0] = iparticle->Px();
+           p[1] = iparticle->Py();
+           p[2] = iparticle->Pz();
+           p[3] = iparticle->Energy();
+
+           
+           origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
+           origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
+           origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
+           
+           Float_t tof = fEventTime + kconv * iparticle->T();
+
+           PushTrack(fTrackIt*trackIt, iparent, kf, 
+                     p[0], p[1], p[2], p[3], 
+                     origin[0], origin[1], origin[2], tof, 
+                     polar[0], polar[1], polar[2],
+                     kPPrimary, nt, 1., ks);
+           fNprimaries++;
+
+           
+           //
+           // Special Treatment to store color-flow
+           //
+           if (ks == 3 || ks == 13 || ks == 14) {
+               TParticle* particle = 0;
+               if (fStack) {
+                   particle = fStack->Particle(nt);
+               } else {
+                   particle = gAlice->Stack()->Particle(nt);
+               }
+//             particle->SetFirstDaughter(fPythia->GetK(2, i));
+//             particle->SetLastDaughter(fPythia->GetK(3, i));         
+           }
+           
+           KeepTrack(nt);
+           pParent[i] = nt;
+           SetHighWaterMark(nt);
+           
+       } // select particle
+    } // particle loop 
+
+    delete[] pParent;
+    
+    return 1;
+}
+
+
+void AliGenPythiaPlus::FinishRun()
+{
+// Print x-section summary
+    fPythia->PrintStatistics();
+
+    if (fNev > 0.) {
+       fQ  /= fNev;
+       fX1 /= fNev;
+       fX2 /= fNev;    
+    }
+    
+    printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
+    printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
+}
+
+void AliGenPythiaPlus::AdjustWeights() const
+{
+// Adjust the weights after generation of all events
+//
+    if (gAlice) {
+       TParticle *part;
+       Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
+       for (Int_t i=0; i<ntrack; i++) {
+           part= gAlice->GetMCApp()->Particle(i);
+           part->SetWeight(part->GetWeight()*fKineBias);
+       }
+    }
+}
+    
+void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
+{
+// Treat protons as inside nuclei with mass numbers a1 and a2  
+
+    fAProjectile = a1;
+    fATarget     = a2;
+    fSetNuclei   = kTRUE;
+}
+
+
+void AliGenPythiaPlus::MakeHeader()
+{
+//
+// Make header for the simulated event
+// 
+  if (gAlice) {
+    if (gAlice->GetEvNumber()>=fDebugEventFirst &&
+       gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
+  }
+
+// Builds the event header, to be called after each event
+    if (fHeader) delete fHeader;
+    fHeader = new AliGenPythiaEventHeader("Pythia");
+//
+// Event type  
+    ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
+//
+// Number of trials
+    ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
+//
+// Event Vertex 
+    fHeader->SetPrimaryVertex(fVertex);
+    
+//
+// Number of primaries
+    fHeader->SetNProduced(fNprimaries);
+//
+// Jets that have triggered
+
+    if (fProcess == kPyJets)
+    {
+       Int_t ntrig, njet;
+       Float_t jets[4][10];
+       GetJets(njet, ntrig, jets);
+
+       
+       for (Int_t i = 0; i < ntrig; i++) {
+           ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
+                                                       jets[3][i]);
+       }
+    }
+//
+// Copy relevant information from external header, if present.
+//
+    Float_t uqJet[4];
+    
+    if (fRL) {
+       AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
+       for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
+       {
+           printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
+           
+           
+           exHeader->TriggerJet(i, uqJet);
+           ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
+       }
+    }
+//
+// Store quenching parameters
+//
+    if (fQuench){
+       Double_t z[4];
+       Double_t xp, yp;
+       if (fQuench == 1) {
+           // Pythia::Quench()
+           fPythia->GetQuenchingParameters(xp, yp, z);
+       } else {
+           // Pyquen
+           Double_t r1 = PARIMP.rb1;
+           Double_t r2 = PARIMP.rb2;
+           Double_t b  = PARIMP.b1;
+           Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
+           Double_t phi = PARIMP.psib1;
+           xp = r * TMath::Cos(phi);
+           yp = r * TMath::Sin(phi);
+           
+       }
+           ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
+           ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
+       }
+//
+// Store pt^hard 
+    ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
+//
+//  Pass header
+//
+    AddHeader(fHeader);
+    fHeader = 0x0;
+}
+
+Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
+{
+// Check the kinematic trigger condition
+//
+    Double_t eta[2];
+    eta[0] = jet1->Eta();
+    eta[1] = jet2->Eta();
+    Double_t phi[2];
+    phi[0] = jet1->Phi();
+    phi[1] = jet2->Phi();
+    Int_t    pdg[2]; 
+    pdg[0] = jet1->GetPdgCode();
+    pdg[1] = jet2->GetPdgCode();    
+    Bool_t   triggered = kFALSE;
+
+    if (fProcess == kPyJets) {
+       Int_t njets = 0;
+       Int_t ntrig = 0;
+       Float_t jets[4][10];
+//
+// Use Pythia clustering on parton level to determine jet axis
+//
+       GetJets(njets, ntrig, jets);
+       
+       if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
+//
+    } else {
+       Int_t ij = 0;
+       Int_t ig = 1;
+       if (pdg[0] == kGamma) {
+           ij = 1;
+           ig = 0;
+       }
+       //Check eta range first...
+       if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
+           (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
+       {
+           //Eta is okay, now check phi range
+           if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
+               (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
+           {
+               triggered = kTRUE;
+           }
+       }
+    }
+    return triggered;
+}
+
+
+
+Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
+//
+//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
+//
+    Bool_t checking = kFALSE;
+    Int_t j, kcode, ks, km;
+    Int_t nPartAcc = 0; //number of particles in the acceptance range
+    Int_t numberOfAcceptedParticles = 1;
+    if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
+    Int_t npart = fParticles->GetEntriesFast();
+    
+    for (j = 0; j<npart; j++) {
+       TParticle *  jparticle = (TParticle *) fParticles->At(j);
+       kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
+       ks = jparticle->GetStatusCode();
+       km = jparticle->GetFirstMother(); 
+       
+       if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
+           nPartAcc++;
+       }
+       if( numberOfAcceptedParticles <= nPartAcc){
+         checking = kTRUE;
+         break;
+       }
+    }
+
+    return checking;
+}
+
+void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
+{
+//
+//  Calls the Pythia jet finding algorithm to find jets in the current event
+//
+//
+//
+//  Save jets
+//
+//  Run Jet Finder
+    fPythia->Pycell(njets);
+    Int_t i;
+    for (i = 0; i < njets; i++) {
+       Float_t px, py, pz, e;
+       fPythia->GetJet(i, px, py, pz, e);
+       jets[0][i] = px;
+       jets[1][i] = py;
+       jets[2][i] = pz;
+       jets[3][i] = e;
+    }
+}
+
+
+void  AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
+{
+//
+//  Calls the Pythia clustering algorithm to find jets in the current event
+//
+    nJets       = 0;
+    nJetsTrig   = 0;
+    if (fJetReconstruction == kCluster) {
+//
+//  Configure cluster algorithm
+//    
+//     fPythia->SetPARU(43, 2.);
+//     fPythia->SetMSTU(41, 1);
+//
+//  Call cluster algorithm
+//    
+       fPythia->Pyclus(nJets);
+//
+//  Loading jets from common block
+//
+    } else {
+
+//
+//  Run Jet Finder
+       fPythia->Pycell(nJets);
+    }
+
+    Int_t i;
+    for (i = 0; i < nJets; i++) {
+       Float_t px, py, pz, e;
+       fPythia->GetJet(i, px, py, pz, e);
+       Float_t pt    = TMath::Sqrt(px * px + py * py);
+       Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);  
+       Float_t theta = TMath::ATan2(pt,pz);
+       Float_t et    = e * TMath::Sin(theta);
+       Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
+       if (
+           eta > fEtaMinJet && eta < fEtaMaxJet && 
+           phi > fPhiMinJet && phi < fPhiMaxJet &&
+           et  > fEtMinJet  && et  < fEtMaxJet     
+           ) 
+       {
+           jets[0][nJetsTrig] = px;
+           jets[1][nJetsTrig] = py;
+           jets[2][nJetsTrig] = pz;
+           jets[3][nJetsTrig] = e;
+           nJetsTrig++;
+       } else {
+       }
+    }
+}
+
+void AliGenPythiaPlus::GetSubEventTime()
+{
+  // Calculates time of the next subevent
+  fEventTime = 0.;
+  if (fEventsTime) {
+    TArrayF &array = *fEventsTime;
+    fEventTime = array[fCurSubEvent++];
+  }
+  //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
+  return;
+}
+
+
+
+
+Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
+{
+  // Is particle in EMCAL acceptance? 
+  // phi in degrees, etamin=-etamax
+  if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi && 
+     eta < fEMCALEta  ) 
+    return kTRUE;
+  else 
+    return kFALSE;
+}
+
+Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
+{
+  // Is particle in PHOS acceptance? 
+  // Acceptance slightly larger considered.
+  // phi in degrees, etamin=-etamax
+  if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi && 
+     eta < fPHOSEta  ) 
+    return kTRUE;
+  else 
+    return kFALSE;
+}
+
+void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
+{
+  //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi 
+  Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
+  Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
+  Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
+  
+  //calculate deltaphi
+  TParticle* ph = (TParticle *) fParticles->At(iphcand);
+  Double_t phphi = ph->Phi();
+  Double_t deltaphi = phiPHOS - phphi;
+
+  
+  
+  //loop for all particles and produce the phi rotation
+  Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+  Double_t oldphi, newphi;
+  Double_t newVx, newVy, R, Vz, time; 
+  Double_t newPx, newPy, pt, Pz, e;
+  for(Int_t i=0; i< np; i++) {
+      TParticle* iparticle = (TParticle *) fParticles->At(i);
+      oldphi = iparticle->Phi();
+      newphi = oldphi + deltaphi;
+      if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle 
+      if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
+      
+      R = iparticle->R();
+      newVx = R*TMath::Cos(newphi);
+      newVy = R*TMath::Sin(newphi);
+      Vz = iparticle->Vz(); // don't transform
+      time = iparticle->T(); // don't transform
+      
+      pt = iparticle->Pt();
+      newPx = pt*TMath::Cos(newphi);
+      newPy = pt*TMath::Sin(newphi);
+      Pz = iparticle->Pz(); // don't transform
+      e = iparticle->Energy(); // don't transform
+      
+      // apply rotation 
+      iparticle->SetProductionVertex(newVx, newVy, Vz, time);
+      iparticle->SetMomentum(newPx, newPy, Pz, e);
+      
+  } //end particle loop 
+  
+   // now let's check that we put correctly the candidate photon in PHOS
+   Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
+   Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax 
+   if(IsInPHOS(phi,eta)) 
+      okdd = kTRUE;
+}
+
+
+#ifdef never
+void AliGenPythiaPlus::Streamer(TBuffer &R__b)
+{
+   // Stream an object of class AliGenPythia.
+
+   if (R__b.IsReading()) {
+      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+      AliGenerator::Streamer(R__b);
+      R__b >> (Int_t&)fProcess;
+      R__b >> (Int_t&)fStrucFunc;
+      R__b >> (Int_t&)fForceDecay;
+      R__b >> fEnergyCMS;
+      R__b >> fKineBias;
+      R__b >> fTrials;
+      fParentSelect.Streamer(R__b);
+      fChildSelect.Streamer(R__b);
+      R__b >> fXsection;
+//      (AliPythia::Instance())->Streamer(R__b);
+      R__b >> fPtHardMin;
+      R__b >> fPtHardMax;
+//      if (fDecayer) fDecayer->Streamer(R__b);
+   } else {
+      R__b.WriteVersion(AliGenPythiaPlus::IsA());
+      AliGenerator::Streamer(R__b);
+      R__b << (Int_t)fProcess;
+      R__b << (Int_t)fStrucFunc;
+      R__b << (Int_t)fForceDecay;
+      R__b << fEnergyCMS;
+      R__b << fKineBias;
+      R__b << fTrials;
+      fParentSelect.Streamer(R__b);
+      fChildSelect.Streamer(R__b);
+      R__b << fXsection;
+//      R__b << fPythia;
+      R__b << fPtHardMin;
+      R__b << fPtHardMax;
+      //     fDecayer->Streamer(R__b);
+   }
+}
+#endif
+
+
+
diff --git a/PYTHIA6/AliGenPythiaPlus.h b/PYTHIA6/AliGenPythiaPlus.h
new file mode 100644 (file)
index 0000000..366e96b
--- /dev/null
@@ -0,0 +1,286 @@
+#ifndef ALIGENPYTHIAPLUS_H
+#define ALIGENPYTHIAPLUS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+/* $Id$ */
+
+//
+// Generator using the TPythia interface (via AliPythia)
+// to generate pp collisions.
+// Using SetNuclei() also nuclear modifications to the structure functions
+// can be taken into account. This makes, of course, only sense for the
+// generation of the products of hard processes (heavy flavor, jets ...)
+//
+// andreas.morsch@cern.ch
+//
+
+#include "AliGenMC.h"
+#include "AliPythia.h"
+
+class AliPythiaBase;
+class TParticle;
+class AliGenPythiaEventHeader;
+class AliGenEventHeader;
+class AliStack;
+class AliRunLoader;
+
+class AliGenPythiaPlus : public AliGenMC
+{
+ public:
+
+    typedef enum {kFlavorSelection, kParentSelection} StackFillOpt_t;
+    typedef enum {kCountAll, kCountParents, kCountTrackables} CountMode_t;
+    typedef enum {kCluster, kCell} JetRecMode_t;
+         
+    AliGenPythiaPlus();
+    AliGenPythiaPlus(AliPythiaBase* pythia);
+    virtual ~AliGenPythiaPlus();
+    virtual void    Generate();
+    virtual void    Init();
+    // Range of events to be printed
+    virtual void    SetEventListRange(Int_t eventFirst=-1, Int_t eventLast=-1);
+    // Select process type
+    virtual void    SetProcess(Process_t proc = kPyCharm) {fProcess = proc;}
+    // Select structure function
+    virtual void    SetStrucFunc(StrucFunc_t func =  kCTEQ5L) {fStrucFunc = func;}
+    // Select pt of hard scattering 
+    virtual void    SetPtHard(Float_t ptmin = 0, Float_t ptmax = 1.e10)
+       {fPtHardMin = ptmin; fPtHardMax = ptmax; }
+    // y of hard scattering
+    virtual void    SetYHard(Float_t ymin = -1.e10, Float_t ymax = 1.e10)
+       {fYHardMin = ymin; fYHardMax = ymax; }
+    // Set initial and final state gluon radiation
+    virtual void    SetGluonRadiation(Int_t iIn, Int_t iFin)
+       {fGinit = iIn; fGfinal = iFin;}
+    // Intrinsic kT
+    virtual void    SetPtKick(Float_t kt = 1.)
+       {fPtKick = kt;}
+    // Use the Pythia 6.3 new multiple interations scenario
+    virtual void    UseNewMultipleInteractionsScenario() {fNewMIS = kTRUE;}
+    // Switch off heavy flavors
+    virtual void    SwitchHFOff() {fHFoff = kTRUE;}
+    // Set centre of mass energy
+    virtual void    SetEnergyCMS(Float_t energy = 5500) {fEnergyCMS = energy;}
+    // Treat protons as inside nuclei with mass numbers a1 and a2
+    virtual void    SetNuclei(Int_t a1, Int_t a2);
+    //
+    // Trigger options
+    //
+    // Energy range for jet trigger
+    virtual void    SetJetEtRange(Float_t etmin = 0., Float_t etmax = 1.e4)
+       {fEtMinJet = etmin; fEtMaxJet = etmax;}
+    // Eta range for jet trigger
+    virtual void    SetJetEtaRange(Float_t etamin = -20., Float_t etamax = 20.)
+       {fEtaMinJet = etamin; fEtaMaxJet = etamax;}
+    // Phi range for jet trigger
+    virtual void    SetJetPhiRange(Float_t phimin = 0., Float_t phimax = 360.)
+       {fPhiMinJet = TMath::Pi()*phimin/180.; fPhiMaxJet = TMath::Pi()*phimax/180.;}
+    // Jet reconstruction mode; default is cone algorithm
+    virtual void    SetJetReconstructionMode(Int_t mode = kCell) {fJetReconstruction = mode;}
+    // Eta range for gamma trigger 
+    virtual void    SetGammaEtaRange(Float_t etamin = -20., Float_t etamax = 20.)
+       {fEtaMinGamma = etamin; fEtaMaxGamma = etamax;}
+    // Phi range for gamma trigger
+    virtual void    SetGammaPhiRange(Float_t phimin = 0., Float_t phimax = 360.)
+       {fPhiMinGamma = TMath::Pi()*phimin/180.; fPhiMaxGamma = TMath::Pi()*phimax/180.;}
+   // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
+    virtual void  SetFragPhotonInCalo(Bool_t b)  {fFragPhotonInCalo = b;}
+    virtual void  SetPi0InCalo       (Bool_t b)  {fPi0InCalo    = b;}
+    virtual void  SetPhotonInCalo(Bool_t b)      {fPhotonInCalo = b;}
+    virtual void  SetCheckPHOS (Bool_t b)        {fCheckPHOS    = b;}
+    virtual void  SetCheckEMCAL(Bool_t b)        {fCheckEMCAL   = b;}
+    virtual void  SetFragPhotonInEMCAL(Bool_t b) {fCheckEMCAL   = b; fFragPhotonInCalo = b;}
+    virtual void  SetFragPhotonInPHOS(Bool_t b)  {fCheckPHOS    = b; fFragPhotonInCalo = b;}
+    virtual void  SetPi0InEMCAL(Bool_t b)        {fCheckEMCAL   = b; fPi0InCalo        = b;}
+    virtual void  SetPi0InPHOS(Bool_t b)         {fCheckPHOS    = b; fPi0InCalo        = b;}
+    virtual void  SetPhotonInEMCAL(Bool_t b)     {fCheckEMCAL   = b; fPhotonInCalo     = b;}
+    virtual void  SetPhotonInPHOS(Bool_t b)      {fCheckPHOS    = b; fPhotonInCalo     = b;}
+    virtual void  SetPhotonInPHOSeta(Bool_t b)   {fCheckPHOSeta = b; fPhotonInCalo     = b;}
+    virtual void  SetFragPhotonOrPi0MinPt(Float_t pt)      {fFragPhotonOrPi0MinPt = pt;}
+    virtual void  SetPhotonMinPt(Float_t pt)     {fPhotonMinPt = pt;}
+    // Trigger and rotate event 
+    void RotatePhi(Int_t iphcand, Bool_t& okdd);
+    // Trigger on a single particle
+    virtual void    SetTriggerParticle(Int_t particle = 0, Float_t etamax = 0.9) 
+       {fTriggerParticle = particle; fTriggerEta = etamax;}
+    //
+    // Heavy flavor options
+    //
+    // Set option for feed down from higher family
+    virtual void SetFeedDownHigherFamily(Bool_t opt) {
+       fFeedDownOpt = opt;
+    }
+    // Set option for selecting particles kept in stack according to flavor
+    // or to parent selection
+    virtual void SetStackFillOpt(StackFillOpt_t opt) {
+       fStackFillOpt = opt;
+    }
+    // Set fragmentation option
+    virtual void SetFragmentation(Bool_t opt) {
+       fFragmentation = opt;
+    }
+    // Set counting mode
+    virtual void SetCountMode(CountMode_t mode) {
+       fCountMode = mode;
+    }
+    //
+    // Quenching
+    //
+    // Set quenching mode 0 = no, 1 = AM, 2 = IL
+    virtual void SetQuench(Int_t flag = 0) {fQuench = flag;}
+    virtual void SetHadronisation(Int_t flag = 1) {fHadronisation = flag;}
+    virtual void SetReadFromFile(const Text_t *filname) {fFileName = filname;  fReadFromFile = 1;}    
+
+    //
+    // Pile-up
+    //
+    // Get interaction rate for pileup studies
+    virtual void    SetInteractionRate(Float_t rate,Float_t timewindow = 90.e-6);
+    virtual Float_t GetInteractionRate() const {return fInteractionRate;}
+    // Get cross section of process
+    virtual Float_t GetXsection() const {return fXsection;}
+    // Get triggered jets
+    void GetJets(Int_t& njets, Int_t& ntrig, Float_t jets[4][10]);
+    void RecJetsUA1(Int_t& njets, Float_t jets[4][50]);
+    void SetPycellParameters(Float_t etamax = 2., Int_t neta = 274, Int_t nphi = 432,
+                            Float_t thresh = 0., Float_t etseed = 4.,
+                            Float_t minet = 10., Float_t r = 1.);
+    
+    // Getters
+    virtual Process_t    GetProcess() const {return fProcess;}
+    virtual StrucFunc_t  GetStrucFunc() const {return fStrucFunc;}
+    virtual void         GetPtHard(Float_t& ptmin, Float_t& ptmax) const
+       {ptmin = fPtHardMin; ptmax = fPtHardMax;}
+    virtual Float_t      GetEnergyCMS() const {return fEnergyCMS;}
+    virtual void         GetNuclei(Int_t&  a1, Int_t& a2) const
+       {a1 = fAProjectile; a2 = fATarget;}
+    virtual void         GetJetEtRange(Float_t& etamin, Float_t& etamax) const
+       {etamin = fEtaMinJet; etamax = fEtaMaxJet;}
+    virtual void         GetJetPhiRange(Float_t& phimin, Float_t& phimax) const
+       {phimin = fPhiMinJet*180./TMath::Pi(); phimax = fPhiMaxJet*180/TMath::Pi();}
+    virtual void         GetGammaEtaRange(Float_t& etamin, Float_t& etamax) const
+       {etamin = fEtaMinGamma; etamax = fEtaMaxGamma;}
+    virtual void         GetGammaPhiRange(Float_t& phimin, Float_t& phimax) const
+       {phimin = fPhiMinGamma*180./TMath::Pi(); phimax = fPhiMaxGamma*180./TMath::Pi();}
+    //
+    Bool_t IsInEMCAL(Float_t phi, Float_t eta);
+    Bool_t IsInPHOS(Float_t phi, Float_t eta);
+    //
+    virtual void FinishRun();
+    Bool_t CheckTrigger(TParticle* jet1, TParticle* jet2);
+    //Used in some processes to selected child properties
+    Bool_t CheckKinematicsOnChild();
+    void     GetSubEventTime();
+
+ protected:
+    // adjust the weight from kinematic cuts
+    void     AdjustWeights() const;
+    Int_t    GenerateMB();
+    void     MakeHeader();    
+    void     GeneratePileup();
+    AliPythiaBase *fPythia;         //!Pythia 
+    Process_t   fProcess;           //Process type
+    StrucFunc_t fStrucFunc;         //Structure Function
+    Float_t     fEnergyCMS;         //Centre of mass energy
+    Float_t     fKineBias;          //!Bias from kinematic selection
+    Int_t       fTrials;            //!Number of trials for current event
+    Int_t       fTrialsRun;         //!Number of trials for run
+    Float_t     fQ;                 //Mean Q
+    Float_t     fX1;                //Mean x1
+    Float_t     fX2;                //Mean x2
+    Float_t     fEventTime;         //Time of the subevent
+    Float_t     fInteractionRate;   //Interaction rate (set by user)
+    Float_t     fTimeWindow;        //Time window for pileup events (set by user)
+    Int_t       fCurSubEvent;       //Index of the current sub-event
+    TArrayF     *fEventsTime;       //Subevents time for pileup
+    Int_t       fNev;               //Number of events 
+    Int_t       fFlavorSelect;      //Heavy Flavor Selection
+    Float_t     fXsection;          //Cross-section
+    Float_t     fPtHardMin;         //lower pT-hard cut 
+    Float_t     fPtHardMax;         //higher pT-hard cut
+    Float_t     fYHardMin;          //lower  y-hard cut 
+    Float_t     fYHardMax;          //higher y-hard cut
+    Int_t       fGinit;             //initial state gluon radiation
+    Int_t       fGfinal;            //final state gluon radiation
+    Int_t       fHadronisation;     //hadronisation
+    Int_t       fNpartons;          //Number of partons before hadronisation
+    Int_t       fReadFromFile;      //read partons from file
+    Int_t       fQuench;            //Flag for quenching
+    Float_t     fPtKick;            //Transverse momentum kick
+    Bool_t      fFullEvent;         //!Write Full event if true
+    AliDecayer  *fDecayer;          //!Pointer to the decayer instance
+    Int_t       fDebugEventFirst;   //!First event to debug
+    Int_t       fDebugEventLast;    //!Last  event to debug
+    Float_t     fEtMinJet;          //Minimum et of triggered Jet
+    Float_t     fEtMaxJet;          //Maximum et of triggered Jet
+    Float_t     fEtaMinJet;         //Minimum eta of triggered Jet
+    Float_t     fEtaMaxJet;         //Maximum eta of triggered Jet
+    Float_t     fPhiMinJet;         //Minimum phi of triggered Jet
+    Float_t     fPhiMaxJet;         //Maximum phi of triggered Jet
+    Int_t       fJetReconstruction; //Jet Reconstruction mode 
+    Float_t     fEtaMinGamma;       // Minimum eta of triggered gamma
+    Float_t     fEtaMaxGamma;       // Maximum eta of triggered gamma
+    Float_t     fPhiMinGamma;       // Minimum phi of triggered gamma
+    Float_t     fPhiMaxGamma;       // Maximum phi of triggered gamma
+    Float_t     fPycellEtaMax;      // Max. eta for Pycell 
+    Int_t       fPycellNEta;        // Number of eta bins for Pycell 
+    Int_t       fPycellNPhi;        // Number of phi bins for Pycell
+    Float_t    fPycellThreshold;   // Pycell threshold
+    Float_t    fPycellEtSeed;      // Pycell seed
+    Float_t    fPycellMinEtJet;    // Pycell min. jet et
+    Float_t    fPycellMaxRadius;   // Pycell cone radius
+    StackFillOpt_t fStackFillOpt;   // Stack filling with all particles with
+                                    // that flavour or only with selected
+                                    // parents and their decays
+    Bool_t fFeedDownOpt;            // Option to set feed down from higher
+                                    // quark families (e.g. b->c)
+    Bool_t  fFragmentation;         // Option to activate fragmentation by Pythia
+    Bool_t  fSetNuclei;             // Flag indicating that SetNuclei has been called
+    Bool_t  fNewMIS;                // Flag for the new multipple interactions scenario
+    Bool_t  fHFoff;                 // Flag for switching heafy flavor production off
+    Int_t   fTriggerParticle;       // Trigger on this particle ...
+    Float_t fTriggerEta;            // .. within |eta| < fTriggerEta
+    CountMode_t fCountMode;         // Options for counting when the event will be finished.     
+    // fCountMode = kCountAll         --> All particles that end up in the
+    //                                    stack are counted
+    // fCountMode = kCountParents     --> Only selected parents are counted
+    // fCountMode = kCountTrackabless --> Only particles flagged for tracking
+    //                                     are counted
+    //
+    //
+
+    AliGenPythiaEventHeader* fHeader;  //! Event header
+    AliRunLoader*            fRL;      //! Run Loader
+    const Text_t* fFileName;           //! Name of file to read from
+
+
+    Bool_t fFragPhotonInCalo; // Option to ask for Fragmentation Photon in calorimeters acceptance
+    Bool_t fPi0InCalo;        // Option to ask for Pi0 in calorimeters acceptance
+    Bool_t fPhotonInCalo;     // Option to ask for Decay Photon in calorimeter acceptance
+    Bool_t fCheckEMCAL;       // Option to ask for FragPhoton or Pi0 in calorimeters EMCAL acceptance
+    Bool_t fCheckPHOS;        // Option to ask for FragPhoton or Pi0 in calorimeters PHOS acceptance
+    Bool_t fCheckPHOSeta;     // Option to ask for PHOS eta acceptance
+    Float_t fFragPhotonOrPi0MinPt; // Minimum momentum of Fragmentation Photon or Pi0
+    Float_t fPhotonMinPt;          // Minimum momentum of Photon 
+    //Calorimeters eta-phi acceptance 
+    Float_t fPHOSMinPhi;           // Minimum phi PHOS
+    Float_t fPHOSMaxPhi;           // Maximum phi PHOS
+    Float_t fPHOSEta;              // Minimum eta PHOS
+    Float_t fEMCALMinPhi;          // Minimum phi EMCAL
+    Float_t fEMCALMaxPhi;          // Maximum phi EMCAL
+    Float_t fEMCALEta;             // Maximum eta EMCAL
+
+ private:
+    AliGenPythiaPlus(const AliGenPythiaPlus &Pythia);
+    AliGenPythiaPlus & operator=(const AliGenPythiaPlus & rhs);
+
+    ClassDef(AliGenPythiaPlus, 1) // AliGenerator interface to Pythia
+};
+#endif
+
+
+
+
+
diff --git a/PYTHIA6/AliPythia6.cxx b/PYTHIA6/AliPythia6.cxx
new file mode 100644 (file)
index 0000000..da77784
--- /dev/null
@@ -0,0 +1,1499 @@
+
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliPythia.cxx,v 1.40 2007/10/09 08:43:24 morsch Exp $ */
+
+#include "AliPythia6.h"
+#include "AliStack.h"
+#include "AliPythiaRndm.h"
+#include "AliFastGlauber.h"
+#include "AliQuenchingWeights.h"
+
+#include "TVector3.h"
+#include "TParticle.h"
+#include "PyquenCommon.h"
+
+ClassImp(AliPythia6)
+
+#ifndef WIN32
+# define pyclus pyclus_
+# define pycell pycell_
+# define pyshow pyshow_
+# define pyrobo pyrobo_
+# define pyquen pyquen_
+# define pyevnw pyevnw_
+# define type_of_call
+#else
+# define pyclus PYCLUS
+# define pycell PYCELL
+# define pyrobo PYROBO
+# define pyquen PYQUEN
+# define pyevnw PYEVNW
+# define type_of_call _stdcall
+#endif
+
+extern "C" void type_of_call pyclus(Int_t & );
+extern "C" void type_of_call pycell(Int_t & );
+extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
+extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyevnw();
+
+
+//_____________________________________________________________________________
+
+AliPythia6* AliPythia6::fgAliPythia=NULL;
+
+AliPythia6::AliPythia6():
+    TPythia6(),
+    AliPythiaBase(),
+    fProcess(kPyMb),
+    fEcms(0.),
+    fStrucFunc(kCTEQ5L),
+    fXJet(0.),
+    fYJet(0.),
+    fNGmax(30),
+    fZmax(0.97),
+    fGlauber(0),
+    fQuenchingWeights(0)
+{
+// Default Constructor
+//
+//  Set random number
+    if (!AliPythiaRndm::GetPythiaRandom()) 
+      AliPythiaRndm::SetPythiaRandom(GetRandom());
+    fGlauber          = 0;
+    fQuenchingWeights = 0;
+}
+
+AliPythia6::AliPythia6(const AliPythia6& pythia):
+    TPythia6(),
+    AliPythiaBase(),
+    fProcess(kPyMb),
+    fEcms(0.),
+    fStrucFunc(kCTEQ5L),
+    fXJet(0.),
+    fYJet(0.),
+    fNGmax(30),
+    fZmax(0.97),
+    fGlauber(0),
+    fQuenchingWeights(0)
+{
+    // Copy Constructor
+    pythia.Copy(*this);
+}
+
+void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
+{
+// Initialise the process to generate 
+    if (!AliPythiaRndm::GetPythiaRandom()) 
+      AliPythiaRndm::SetPythiaRandom(GetRandom());
+    
+    fProcess = process;
+    fEcms = energy;
+    fStrucFunc = strucfunc;
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+    SetMDCY(Pycomp(111) ,1,0);
+    SetMDCY(Pycomp(310) ,1,0);
+    SetMDCY(Pycomp(3122),1,0);
+    SetMDCY(Pycomp(3112),1,0);
+    SetMDCY(Pycomp(3212),1,0);
+    SetMDCY(Pycomp(3222),1,0);
+    SetMDCY(Pycomp(3312),1,0);
+    SetMDCY(Pycomp(3322),1,0);
+    SetMDCY(Pycomp(3334),1,0);
+    // Select structure function 
+    SetMSTP(52,2);
+    SetMSTP(51,AliStructFuncType::PDFsetIndex(strucfunc));
+    // Particles produced in string fragmentation point directly to either of the two endpoints
+    // of the string (depending in the side they were generated from).
+    SetMSTU(16,2);
+
+//
+// Pythia initialisation for selected processes//
+//
+// Make MSEL clean
+//
+    for (Int_t i=1; i<= 200; i++) {
+       SetMSUB(i,0);
+    }
+//  select charm production
+    switch (process) 
+    {
+    case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
+//        Multiple interactions on.
+       SetMSTP(81,1);
+// Double Gaussian matter distribution.
+       SetMSTP(82,4);
+       SetPARP(83,0.5);
+       SetPARP(84,0.4);
+//  pT0.
+       SetPARP(82,2.0);
+//  Reference energy for pT0 and energy rescaling pace.
+       SetPARP(89,1800);
+       SetPARP(90,0.25);
+//  String drawing almost completely minimizes string length.
+       SetPARP(85,0.9);
+       SetPARP(86,0.95);
+// ISR and FSR activity.
+       SetPARP(67,4);
+       SetPARP(71,4);
+// Lambda_FSR scale.
+       SetPARJ(81,0.29);
+       break;
+    case kPyOldUEQ2ordered2:   
+// Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+       SetMSTP(81,1);
+// Double Gaussian matter distribution.
+       SetMSTP(82,4);
+       SetPARP(83,0.5);
+       SetPARP(84,0.4);
+// pT0.
+       SetPARP(82,2.0);
+// Reference energy for pT0 and energy rescaling pace.
+       SetPARP(89,1800);
+       SetPARP(90,0.16);  // here is the difference with  kPyOldUEQ2ordered
+// String drawing almost completely minimizes string length.
+       SetPARP(85,0.9);
+       SetPARP(86,0.95);
+// ISR and FSR activity.
+       SetPARP(67,4);
+       SetPARP(71,4);
+// Lambda_FSR scale.
+       SetPARJ(81,0.29);       
+       break;
+    case kPyOldPopcorn:  
+// Old production mechanism: Old Popcorn
+       SetMSEL(1);
+       SetMSTJ(12,3); 
+// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
+       SetMSTP(88,2); 
+// (D=1)see can be used to form  baryons (BARYON JUNCTION)
+       SetMSTJ(1,1);  
+       AtlasTuning();
+       break;
+    case kPyCharm:
+       SetMSEL(4);
+//  heavy quark masses
+
+       SetPMAS(4,1,1.2);
+//
+//    primordial pT
+       SetMSTP(91,1);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
+//
+       break;
+    case kPyBeauty:
+       SetMSEL(5);
+       SetPMAS(5,1,4.75);
+       break;
+    case kPyJpsi:
+       SetMSEL(0);
+// gg->J/Psi g
+       SetMSUB(86,1);
+       break;
+    case kPyJpsiChi:
+       SetMSEL(0);
+// gg->J/Psi g
+       SetMSUB(86,1);
+// gg-> chi_0c g
+       SetMSUB(87,1);
+// gg-> chi_1c g
+       SetMSUB(88,1);
+// gg-> chi_2c g
+       SetMSUB(89,1);  
+       break;
+    case kPyCharmUnforced:
+       SetMSEL(0);
+// gq->qg   
+       SetMSUB(28,1);
+// gg->qq
+       SetMSUB(53,1);
+// gg->gg
+       SetMSUB(68,1);
+       break;
+    case kPyBeautyUnforced:
+       SetMSEL(0);
+// gq->qg   
+       SetMSUB(28,1);
+// gg->qq
+       SetMSUB(53,1);
+// gg->gg
+       SetMSUB(68,1);
+       break;
+    case kPyMb:
+// Minimum Bias pp-Collisions
+//
+//   
+//      select Pythia min. bias model
+       SetMSEL(0);
+       SetMSUB(92,1);             // single diffraction AB-->XB
+       SetMSUB(93,1);             // single diffraction AB-->AX
+       SetMSUB(94,1);             // double diffraction
+       SetMSUB(95,1);             // low pt production
+
+       AtlasTuning();
+       break;
+    case kPyMbDefault:
+// Minimum Bias pp-Collisions
+//
+//   
+//      select Pythia min. bias model
+       SetMSEL(0);
+       SetMSUB(92,1);             // single diffraction AB-->XB
+       SetMSUB(93,1);             // single diffraction AB-->AX
+       SetMSUB(94,1);             // double diffraction
+       SetMSUB(95,1);             // low pt production
+
+       break;
+    case kPyLhwgMb:
+// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
+//  -> Pythia 6.3 or above is needed
+//   
+       SetMSEL(0);
+       SetMSUB(92,1);             // single diffraction AB-->XB
+       SetMSUB(93,1);             // single diffraction AB-->AX
+       SetMSUB(94,1);             // double diffraction
+       SetMSUB(95,1);             // low pt production
+       SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));
+       SetMSTP(52,2);
+       SetMSTP(68,1);
+       SetMSTP(70,2);
+       SetMSTP(81,1);             // Multiple Interactions ON
+       SetMSTP(82,4);             // Double Gaussian Model
+       SetMSTP(88,1);
+
+       SetPARP(82,2.3);           // [GeV]    PT_min at Ref. energy
+       SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
+       SetPARP(84,0.5);           // Core radius
+       SetPARP(85,0.9);           // Regulates gluon prod. mechanism
+       SetPARP(90,0.2);           // 2*epsilon (exponent in power law)
+
+       break;
+    case kPyMbNonDiffr:
+// Minimum Bias pp-Collisions
+//
+//   
+//      select Pythia min. bias model
+       SetMSEL(0);
+       SetMSUB(95,1);             // low pt production
+
+       AtlasTuning();
+       break;
+    case kPyMbMSEL1:
+       ConfigHeavyFlavor();
+// Intrinsic <kT^2>
+        SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
+        SetPARP(91,1.);     // <kT^2> = PARP(91,1.)^2
+        SetPARP(93,5.);     // Upper cut-off
+// Set Q-quark mass
+        SetPMAS(4,1,1.2);   // Charm quark mass
+        SetPMAS(5,1,4.78);  // Beauty quark mass
+       SetPARP(71,4.);     // Defaut value
+// Atlas Tuning
+       AtlasTuning();
+       break;
+    case kPyJets:
+//
+//  QCD Jets
+//
+       SetMSEL(1);
+ // Pythia Tune A (CDF)
+ //
+       SetPARP(67,2.5);           // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
+       SetMSTP(82,4);             // Double Gaussian Model
+       SetPARP(82,2.0);           // [GeV]    PT_min at Ref. energy
+       SetPARP(84,0.4);           // Core radius
+       SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
+       SetPARP(86,0.95);          // Regulates gluon prod. mechanism
+       SetPARP(89,1800.);         // [GeV]   Ref. energy
+       SetPARP(90,0.25);          // 2*epsilon (exponent in power law)
+       break;
+    case kPyDirectGamma:
+       SetMSEL(10);
+       break;
+    case kPyCharmPbPbMNR:
+    case kPyD0PbPbMNR:
+    case kPyDPlusPbPbMNR:
+    case kPyDPlusStrangePbPbMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // c-cbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with Pb-Pb collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+       ConfigHeavyFlavor();
+      // Intrinsic <kT>
+      SetMSTP(91,1);
+      SetPARP(91,1.304);
+      SetPARP(93,6.52);
+      // Set c-quark mass
+      SetPMAS(4,1,1.2);
+      break;
+    case kPyCharmpPbMNR:
+    case kPyD0pPbMNR:
+    case kPyDPluspPbMNR:
+    case kPyDPlusStrangepPbMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // c-cbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with p-Pb collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+       ConfigHeavyFlavor();
+      // Intrinsic <kT>
+       SetMSTP(91,1);
+       SetPARP(91,1.16);
+       SetPARP(93,5.8);
+       
+      // Set c-quark mass
+       SetPMAS(4,1,1.2);
+      break;
+    case kPyCharmppMNR:
+    case kPyD0ppMNR:
+    case kPyDPlusppMNR:
+    case kPyDPlusStrangeppMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // c-cbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with pp collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+       ConfigHeavyFlavor();
+      // Intrinsic <kT^2>
+       SetMSTP(91,1);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
+       
+      // Set c-quark mass
+       SetPMAS(4,1,1.2);
+      break;
+    case kPyCharmppMNRwmi:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // c-cbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with pp collisions
+      // and with kCTEQ5L PDFs.
+      // Added multiple interactions according to ATLAS tune settings.
+      // To get a "reasonable" agreement with MNR results, events have to be 
+      // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+      // set to 2.76 GeV.
+      // To get a "perfect" agreement with MNR results, events have to be 
+      // generated in four ptHard bins with the following relative 
+      // normalizations:
+      // 2.76-3 GeV: 25%
+      //    3-4 GeV: 40%
+      //    4-8 GeV: 29%
+      //     >8 GeV:  6%
+       ConfigHeavyFlavor();
+      // Intrinsic <kT^2>
+       SetMSTP(91,1);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
+
+      // Set c-quark mass
+       SetPMAS(4,1,1.2);
+       AtlasTuning();
+       break;
+    case kPyBeautyPbPbMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // b-bbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with Pb-Pb collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+       ConfigHeavyFlavor();
+      // QCD scales
+       SetPARP(67,1.0);
+       SetPARP(71,1.0);
+      // Intrinsic <kT>
+       SetMSTP(91,1);
+       SetPARP(91,2.035);
+       SetPARP(93,10.17);
+      // Set b-quark mass
+       SetPMAS(5,1,4.75);
+      break;
+    case kPyBeautypPbMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // b-bbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with p-Pb collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+       ConfigHeavyFlavor();
+      // QCD scales
+       SetPARP(67,1.0);
+       SetPARP(71,1.0);
+      // Intrinsic <kT>
+       SetMSTP(91,1);
+       SetPARP(91,1.60);
+       SetPARP(93,8.00);
+      // Set b-quark mass
+       SetPMAS(5,1,4.75);
+      break;
+    case kPyBeautyppMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // b-bbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with pp collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+       ConfigHeavyFlavor();
+      // QCD scales
+       SetPARP(67,1.0);
+       SetPARP(71,1.0);
+       
+       // Intrinsic <kT>
+       SetMSTP(91,1);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
+       
+       // Set b-quark mass
+       SetPMAS(5,1,4.75);
+      break;
+     case kPyBeautyppMNRwmi:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // b-bbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with pp collisions
+      // and with kCTEQ5L PDFs.
+      // Added multiple interactions according to ATLAS tune settings.
+      // To get a "reasonable" agreement with MNR results, events have to be 
+      // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+      // set to 2.76 GeV.
+      // To get a "perfect" agreement with MNR results, events have to be 
+      // generated in four ptHard bins with the following relative 
+      // normalizations:
+      // 2.76-4 GeV:  5% 
+      //    4-6 GeV: 31%
+      //    6-8 GeV: 28%
+      //     >8 GeV: 36%
+        ConfigHeavyFlavor();
+      // QCD scales
+        SetPARP(67,1.0);
+        SetPARP(71,1.0);
+        
+        // Intrinsic <kT>
+        SetMSTP(91,1);
+        SetPARP(91,1.);
+        SetPARP(93,5.);
+
+      // Set b-quark mass
+        SetPMAS(5,1,4.75);
+
+        AtlasTuning();
+        break; 
+    case kPyW:
+
+      //Inclusive production of W+/-
+      SetMSEL(0);
+      //f fbar -> W+ 
+      SetMSUB(2,1);
+      //       //f fbar -> g W+
+      //       SetMSUB(16,1);
+      //       //f fbar -> gamma W+
+      //       SetMSUB(20,1);
+      //       //f g -> f W+  
+      //       SetMSUB(31,1);
+      //       //f gamma -> f W+
+      //       SetMSUB(36,1);
+      
+      // Initial/final parton shower on (Pythia default)
+      // With parton showers on we are generating "W inclusive process"
+      SetMSTP(61,1); //Initial QCD & QED showers on
+      SetMSTP(71,1); //Final QCD & QED showers on
+      
+      break;  
+
+    case kPyZ:
+
+      //Inclusive production of Z
+      SetMSEL(0);
+      //f fbar -> Z/gamma
+      SetMSUB(1,1);
+      
+      //       // f fbar -> g Z/gamma
+      //       SetMSUB(15,1);
+      //       // f fbar -> gamma Z/gamma
+      //       SetMSUB(19,1);
+      //       // f g -> f Z/gamma
+      //       SetMSUB(30,1);
+      //       // f gamma -> f Z/gamma
+      //       SetMSUB(35,1);
+      
+      //only Z included, not gamma
+      SetMSTP(43,2);
+      
+      // Initial/final parton shower on (Pythia default)
+      // With parton showers on we are generating "Z inclusive process"
+      SetMSTP(61,1); //Initial QCD & QED showers on
+      SetMSTP(71,1); //Final QCD & QED showers on
+      
+      break;  
+
+    }
+//
+//  Initialize PYTHIA
+    SetMSTP(41,1);   // all resonance decays switched on
+    Initialize("CMS","p","p",fEcms);
+    
+}
+
+Int_t AliPythia6::CheckedLuComp(Int_t kf)
+{
+// Check Lund particle code (for debugging)
+    Int_t kc=Pycomp(kf);
+    return kc;
+}
+
+void AliPythia6::SetNuclei(Int_t a1, Int_t a2)
+{
+// Treat protons as inside nuclei with mass numbers a1 and a2  
+//    The MSTP array in the PYPARS common block is used to enable and 
+//    select the nuclear structure functions. 
+//    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
+//            =1: internal PYTHIA acording to MSTP(51) 
+//            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
+//    If the following mass number both not equal zero, nuclear corrections of the stf are used.
+//    MSTP(192) : Mass number of nucleus side 1
+//    MSTP(193) : Mass number of nucleus side 2
+    SetMSTP(52,2);
+    SetMSTP(192, a1);
+    SetMSTP(193, a2);  
+}
+       
+
+AliPythia6* AliPythia6::Instance()
+{ 
+// Set random number generator 
+    if (fgAliPythia) {
+       return fgAliPythia;
+    } else {
+       fgAliPythia = new AliPythia6();
+       return fgAliPythia;
+    }
+}
+
+void AliPythia6::PrintParticles()
+{ 
+// Print list of particl properties
+    Int_t np = 0;
+    char*   name = new char[16];    
+    for (Int_t kf=0; kf<1000000; kf++) {
+       for (Int_t c = 1;  c > -2; c-=2) {
+           Int_t kc = Pycomp(c*kf);
+           if (kc) {
+               Float_t mass  = GetPMAS(kc,1);
+               Float_t width = GetPMAS(kc,2);  
+               Float_t tau   = GetPMAS(kc,4);
+
+               Pyname(kf,name);
+       
+               np++;
+               
+               printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
+                      c*kf, name, mass, width, tau);
+           }
+       }
+    }
+    printf("\n Number of particles %d \n \n", np);
+}
+
+void  AliPythia6::ResetDecayTable()
+{
+//  Set default values for pythia decay switches
+    Int_t i;
+    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
+    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
+}
+
+void  AliPythia6::SetDecayTable()
+{
+//  Set default values for pythia decay switches
+//
+    Int_t i;
+    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
+    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
+}
+
+void  AliPythia6::Pyclus(Int_t& njet)
+{
+//  Call Pythia clustering algorithm
+//
+    pyclus(njet);
+}
+
+void  AliPythia6::Pycell(Int_t& njet)
+{
+//  Call Pythia jet reconstruction algorithm
+//
+    pycell(njet);
+}
+
+void AliPythia6::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
+{
+    // Get jet number i
+    Int_t n = GetN();
+    px    = GetPyjets()->P[0][n+i];
+    py    = GetPyjets()->P[1][n+i];
+    pz    = GetPyjets()->P[2][n+i];
+    e     = GetPyjets()->P[3][n+i];
+}
+
+void  AliPythia6::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
+{
+//  Call Pythia showering
+//
+    pyshow(ip1, ip2, qmax);
+}
+
+void AliPythia6::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
+{
+    pyrobo(imi, ima, the, phi, bex, bey, bez);
+}
+
+
+
+void AliPythia6::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
+{
+// Initializes 
+// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
+// (2) The nuclear geometry using the Glauber Model
+//     
+    
+    fGlauber = new AliFastGlauber();
+    fGlauber->Init(2);
+    fGlauber->SetCentralityClass(cMin, cMax); 
+
+    fQuenchingWeights = new AliQuenchingWeights();
+    fQuenchingWeights->InitMult();
+    fQuenchingWeights->SetK(k);
+    fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
+    fNGmax = ngmax;
+    fZmax  = zmax;
+    
+}
+
+
+void  AliPythia6::Quench()
+{
+//
+//
+//  Simple Jet Quenching routine:
+//  =============================
+//  The jet formed by all final state partons radiated by the parton created 
+//  in the hard collisions is quenched by a factor (1-z) using light cone variables in 
+//  the initial parton reference frame:
+//  (E + p_z)new = (1-z) (E + p_z)old
+//
+//
+//
+//
+//   The lost momentum is first balanced by one gluon with virtuality > 0.   
+//   Subsequently the gluon splits to yield two gluons with E = p.
+//
+//
+// 
+    static Float_t eMean = 0.;
+    static Int_t   icall = 0;
+    
+    Double_t p0[4][5];
+    Double_t p1[4][5];
+    Double_t p2[4][5];
+    Int_t   klast[4] = {-1, -1, -1, -1};
+
+    Int_t numpart   = fPyjets->N;
+    Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
+    Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
+    Bool_t  quenched[4];
+    Double_t wjtKick[4];
+    Int_t nGluon[4];
+    Int_t qPdg[4];
+    Int_t   imo, kst, pdg;
+    
+//
+//  Sore information about Primary partons
+//
+//  j =
+//  0, 1 partons from hard scattering
+//  2, 3 partons from initial state radiation
+// 
+    for (Int_t i = 2; i <= 7; i++) {
+       Int_t j = 0;
+       // Skip gluons that participate in hard scattering
+       if (i == 4 || i == 5) continue;
+       // Gluons from hard Scattering
+       if (i == 6 || i == 7) {
+           j = i - 6;
+           pxq[j]    = fPyjets->P[0][i];
+           pyq[j]    = fPyjets->P[1][i];
+           pzq[j]    = fPyjets->P[2][i];
+           eq[j]     = fPyjets->P[3][i];
+           mq[j]     = fPyjets->P[4][i];
+       } else {
+           // Gluons from initial state radiation
+           //
+           // Obtain 4-momentum vector from difference between original parton and parton after gluon 
+           // radiation. Energy is calculated independently because initial state radition does not 
+           // conserve strictly momentum and energy for each partonic system independently.
+           //
+           // Not very clean. Should be improved !
+           //
+           //
+           j = i;
+           pxq[j]    = fPyjets->P[0][i] - fPyjets->P[0][i+2];
+           pyq[j]    = fPyjets->P[1][i] - fPyjets->P[1][i+2];
+           pzq[j]    = fPyjets->P[2][i] - fPyjets->P[2][i+2];
+           mq[j]     = fPyjets->P[4][i];
+           eq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
+       }
+//
+//  Calculate some kinematic variables
+//
+       yq[j]     = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
+       pq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
+       phiq[j]   = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
+       ptq[j]    = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
+       thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
+       qPdg[j]   =  fPyjets->K[1][i];
+    }
+  
+    Double_t int0[4];
+    Double_t int1[4];
+    
+    fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
+
+    for (Int_t j = 0; j < 4; j++) {
+       //
+       // Quench only central jets and with E > 10.
+       //
+
+
+       Int_t itype = (qPdg[j] == 21) ? 2 : 1;
+       Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
+
+       if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
+           fZQuench[j] = 0.;
+       } else {
+           if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
+               icall ++;
+               eMean += eloss;
+           }
+           //
+           // Extra pt
+           Double_t l =   fQuenchingWeights->CalcLk(int0[j], int1[j]);     
+           wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->CalcQk(int0[j], int1[j]));
+           //
+           // Fractional energy loss
+           fZQuench[j] = eloss / eq[j];
+           //
+           // Avoid complete loss
+           //
+           if (fZQuench[j] == 1.) fZQuench[j] = fZmax;
+           //
+           // Some debug printing
+
+           
+//         printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n", 
+//                j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
+           
+//         fZQuench[j] = 0.8;
+//         while (fZQuench[j] >= 0.95)  fZQuench[j] = gRandom->Exp(0.2);
+       }
+       
+       quenched[j] = (fZQuench[j] > 0.01);
+    } // primary partons
+    
+    
+
+    Double_t pNew[1000][4];
+    Int_t    kNew[1000];
+    Int_t icount = 0;
+    Double_t zquench[4];
+    
+//
+//  System Loop    
+    for (Int_t isys = 0; isys < 4; isys++) {
+//      Skip to next system if not quenched.
+       if (!quenched[isys]) continue;
+       
+       nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
+       if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
+       zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
+       wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
+
+
+       
+       Int_t igMin = -1;
+       Int_t igMax = -1;
+       Double_t pg[4] = {0., 0., 0., 0.};
+       
+//
+// Loop on radiation events
+
+       for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
+           while (1) {
+               icount = 0;
+               for (Int_t k = 0; k < 4; k++)
+               {
+                   p0[isys][k] = 0.;
+                   p1[isys][k] = 0.;
+                   p2[isys][k] = 0.;
+               }
+//      Loop over partons
+               for (Int_t i = 0; i < numpart; i++)
+               {
+                   imo =  fPyjets->K[2][i];
+                   kst =  fPyjets->K[0][i];
+                   pdg =  fPyjets->K[1][i];
+                   
+               
+               
+//      Quarks and gluons only
+                   if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
+//      Particles from hard scattering only
+                   
+                   if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
+                   Int_t imom = imo % 1000;
+                   if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
+                   if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;               
+                   
+                   
+//      Skip comment lines
+                   if (kst != 1 && kst != 2) continue;
+//
+//      Parton kinematic
+                   px    = fPyjets->P[0][i];
+                   py    = fPyjets->P[1][i];
+                   pz    = fPyjets->P[2][i];
+                   e     = fPyjets->P[3][i];
+                   m     = fPyjets->P[4][i];
+                   pt    = TMath::Sqrt(px * px + py * py);
+                   p     = TMath::Sqrt(px * px + py * py + pz * pz); 
+                   phi   = TMath::Pi() + TMath::ATan2(-py, -px);
+                   theta = TMath::ATan2(pt, pz);
+               
+//
+//      Save 4-momentum sum for balancing
+                   Int_t index = isys;
+                   
+                   p0[index][0] += px;
+                   p0[index][1] += py;
+                   p0[index][2] += pz;
+                   p0[index][3] += e;
+               
+                   klast[index] = i;
+                   
+//
+//      Fractional energy loss
+                   Double_t z = zquench[index];
+                   
+                   
+//      Don't fully quench radiated gluons
+//
+                   if (imo > 1000) {
+//      This small factor makes sure that the gluons are not too close in phase space to avoid recombination
+//
+
+                       z = 0.02;
+                   }
+//                 printf("z: %d %f\n", imo, z);
+                   
+
+//
+                   
+                   //
+                   //
+                   //      Transform into frame in which initial parton is along z-axis
+                   //
+                   TVector3 v(px, py, pz);
+                   v.RotateZ(-phiq[index]);  v.RotateY(-thetaq[index]);
+                   Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl  = v.Z();
+
+                   Double_t jt  = TMath::Sqrt(pxs * pxs + pys * pys);
+                   Double_t mt2 = jt * jt + m * m;
+                   Double_t zmax = 1.;     
+                   //
+                   // Kinematic limit on z
+                   //
+                   if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
+                   //
+                   // Change light-cone kinematics rel. to initial parton
+                   //  
+                   Double_t eppzOld = e + pl;
+                   Double_t empzOld = e - pl;
+                   
+                   Double_t eppzNew = (1. - z) * eppzOld;
+                   Double_t empzNew = empzOld - mt2 * z / eppzOld;
+                   Double_t eNew    = 0.5 * (eppzNew + empzNew);
+                   Double_t plNew   = 0.5 * (eppzNew - empzNew);
+                   
+                   Double_t jtNew;
+                   //
+                   // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
+                   Double_t mt2New = eppzNew * empzNew;
+                   if (mt2New < 1.e-8) mt2New = 0.;
+                   if (z < zmax) {
+                       if (m * m > mt2New) {
+                           //
+                           // This should not happen 
+                           //
+                           Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
+                           jtNew = 0;
+                       } else {
+                           jtNew    = TMath::Sqrt(mt2New - m * m);
+                       }
+                   } else {
+                       // If pT is to small (probably a leading massive particle) we scale only the energy
+                       // This can cause negative masses of the radiated gluon
+                       // Let's hope for the best ...
+                       jtNew = jt;
+                       eNew  = TMath::Sqrt(plNew * plNew + mt2);
+                       
+                   }
+                   //
+                   //     Calculate new px, py
+                   //
+                   Double_t pxNew   = jtNew / jt * pxs;
+                   Double_t pyNew   = jtNew / jt * pys;        
+                   
+//                 Double_t dpx = pxs - pxNew;
+//                 Double_t dpy = pys - pyNew;
+//                 Double_t dpz = pl  - plNew;
+//                 Double_t de  = e   - eNew;
+//                 Double_t dmass2 = de * de  - dpx * dpx - dpy * dpy - dpz * dpz;
+//                 printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
+//                 printf("New mass (2) %e %e \n", pxNew, pyNew);
+                   //
+                   //      Rotate back
+                   //  
+                   TVector3 w(pxNew, pyNew, plNew);
+                   w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
+                   pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
+               
+                   p1[index][0] += pxNew;
+                   p1[index][1] += pyNew;
+                   p1[index][2] += plNew;
+                   p1[index][3] += eNew;       
+                   //
+                   // Updated 4-momentum vectors
+                   //
+                   pNew[icount][0]  = pxNew;
+                   pNew[icount][1]  = pyNew;
+                   pNew[icount][2]  = plNew;
+                   pNew[icount][3]  = eNew;
+                   kNew[icount]     = i;
+                   icount++;
+               } // parton loop
+               //
+               // Check if there was phase-space for quenching
+               //
+
+               if (icount == 0) quenched[isys] = kFALSE;
+               if (!quenched[isys]) break;
+               
+               for (Int_t j = 0; j < 4; j++) 
+               {
+                   p2[isys][j] = p0[isys][j] - p1[isys][j];
+               }
+               p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
+               if (p2[isys][4] > 0.) {
+                   p2[isys][4] = TMath::Sqrt(p2[isys][4]);
+                   break;
+               } else {
+                   printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
+                   printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
+                   if (p2[isys][4] < -0.01) {
+                       printf("Negative mass squared !\n");
+                       // Here we have to put the gluon back to mass shell
+                       // This will lead to a small energy imbalance
+                       p2[isys][4]  = 0.;
+                       p2[isys][3]  = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
+                       break;
+                   } else {
+                       p2[isys][4] = 0.;
+                       break;
+                   }
+               }
+               /*
+               zHeavy *= 0.98;
+               printf("zHeavy lowered to %f\n", zHeavy);
+               if (zHeavy < 0.01) {
+                   printf("No success ! \n");
+                   icount = 0;
+                   quenched[isys] = kFALSE;
+                   break;
+               }
+               */
+           } // iteration on z (while)
+           
+//         Update  event record
+           for (Int_t k = 0; k < icount; k++) {
+//             printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
+               fPyjets->P[0][kNew[k]] = pNew[k][0];
+               fPyjets->P[1][kNew[k]] = pNew[k][1];
+               fPyjets->P[2][kNew[k]] = pNew[k][2];
+               fPyjets->P[3][kNew[k]] = pNew[k][3];
+           }
+           //
+           // Add the gluons
+           //
+           Int_t ish = 0;    
+           Int_t iGlu;
+           if (!quenched[isys]) continue;
+//
+//      Last parton from shower i
+           Int_t in = klast[isys];
+//
+//      Continue if no parton in shower i selected
+           if (in == -1) continue;
+//  
+//      If this is the second initial parton and it is behind the first move pointer by previous ish
+           if (isys == 1 && klast[1] > klast[0]) in += ish;
+//
+//      Starting index
+           
+//         jmin = in - 1;
+// How many additional gluons will be generated
+           ish  = 1;
+           if (p2[isys][4] > 0.05) ish = 2;
+//
+//      Position of gluons
+           iGlu = numpart;
+           if (iglu == 0) igMin = iGlu;
+           igMax = iGlu;
+           numpart += ish;
+           (fPyjets->N) += ish;
+           
+           if (ish == 1) {
+               fPyjets->P[0][iGlu] = p2[isys][0];
+               fPyjets->P[1][iGlu] = p2[isys][1];
+               fPyjets->P[2][iGlu] = p2[isys][2];
+               fPyjets->P[3][iGlu] = p2[isys][3];
+               fPyjets->P[4][iGlu] = p2[isys][4];
+               
+               fPyjets->K[0][iGlu] = 1;
+               if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
+               fPyjets->K[1][iGlu] = 21;       
+               fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
+               fPyjets->K[3][iGlu] = -1;       
+               fPyjets->K[4][iGlu] = -1;
+               
+               pg[0] += p2[isys][0];
+               pg[1] += p2[isys][1];
+               pg[2] += p2[isys][2];
+               pg[3] += p2[isys][3];
+           } else {
+               //
+               // Split gluon in rest frame.
+               //
+               Double_t bx   =  p2[isys][0] / p2[isys][3];
+               Double_t by   =  p2[isys][1] / p2[isys][3];
+               Double_t bz   =  p2[isys][2] / p2[isys][3];
+               Double_t pst  =  p2[isys][4] / 2.;
+               //
+               // Isotropic decay ????
+               Double_t cost = 2. * gRandom->Rndm() - 1.;
+               Double_t sint = TMath::Sqrt(1. - cost * cost);
+               Double_t phi =  2. * TMath::Pi() * gRandom->Rndm();
+               
+               Double_t pz1 =   pst * cost;
+               Double_t pz2 =  -pst * cost;
+               Double_t pt1 =   pst * sint;
+               Double_t pt2 =  -pst * sint;
+               Double_t px1 =   pt1 * TMath::Cos(phi);
+               Double_t py1 =   pt1 * TMath::Sin(phi);     
+               Double_t px2 =   pt2 * TMath::Cos(phi);
+               Double_t py2 =   pt2 * TMath::Sin(phi);     
+               
+               fPyjets->P[0][iGlu] = px1;
+               fPyjets->P[1][iGlu] = py1;
+               fPyjets->P[2][iGlu] = pz1;
+               fPyjets->P[3][iGlu] = pst;
+               fPyjets->P[4][iGlu] = 0.;
+               
+               fPyjets->K[0][iGlu] = 1 ;
+               fPyjets->K[1][iGlu] = 21;       
+               fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
+               fPyjets->K[3][iGlu] = -1;       
+               fPyjets->K[4][iGlu] = -1;
+               
+               fPyjets->P[0][iGlu+1] = px2;
+               fPyjets->P[1][iGlu+1] = py2;
+               fPyjets->P[2][iGlu+1] = pz2;
+               fPyjets->P[3][iGlu+1] = pst;
+               fPyjets->P[4][iGlu+1] = 0.;
+               
+               fPyjets->K[0][iGlu+1] = 1;
+               if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
+               fPyjets->K[1][iGlu+1] = 21;     
+               fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
+               fPyjets->K[3][iGlu+1] = -1;     
+               fPyjets->K[4][iGlu+1] = -1;
+               SetMSTU(1,0);
+               SetMSTU(2,0);
+               //
+               // Boost back
+               //
+               Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
+           }
+/*    
+           for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
+               Double_t px, py, pz;
+               px = fPyjets->P[0][ig]; 
+               py = fPyjets->P[1][ig]; 
+               pz = fPyjets->P[2][ig]; 
+               TVector3 v(px, py, pz);
+               v.RotateZ(-phiq[isys]);
+               v.RotateY(-thetaq[isys]);
+               Double_t pxs     = v.X(); Double_t pys = v.Y(); Double_t pzs  = v.Z();     
+               Double_t r       = AliPythiaRndm::GetPythiaRandom()->Rndm();
+               Double_t jtKick  = 0.3 * TMath::Sqrt(-TMath::Log(r));
+               if (ish == 2)   jtKick  = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
+               Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
+               pxs += jtKick * TMath::Cos(phiKick);
+               pys += jtKick * TMath::Sin(phiKick);
+               TVector3 w(pxs, pys, pzs);
+               w.RotateY(thetaq[isys]);
+               w.RotateZ(phiq[isys]);
+               fPyjets->P[0][ig] = w.X(); 
+               fPyjets->P[1][ig] = w.Y(); 
+               fPyjets->P[2][ig] = w.Z(); 
+               fPyjets->P[2][ig] = w.Mag();
+           }
+*/
+       } // kGluon         
+       
+       
+    // Check energy conservation
+       Double_t pxs = 0.;
+       Double_t pys = 0.;
+       Double_t pzs = 0.;      
+       Double_t es  = 14000.;
+       
+       for (Int_t i = 0; i < numpart; i++)
+       {
+           kst =  fPyjets->K[0][i];
+           if (kst != 1 && kst != 2) continue;
+           pxs += fPyjets->P[0][i];
+           pys += fPyjets->P[1][i];
+           pzs += fPyjets->P[2][i];        
+           es  -= fPyjets->P[3][i];        
+       }
+       if (TMath::Abs(pxs) > 1.e-2 ||
+           TMath::Abs(pys) > 1.e-2 ||
+           TMath::Abs(pzs) > 1.e-1) {
+           printf("%e %e %e %e\n", pxs, pys, pzs, es);
+//             Fatal("Quench()", "4-Momentum non-conservation");
+       }
+       
+    } // end quenching loop (systems)
+// Clean-up
+    for (Int_t i = 0; i < numpart; i++)
+    {
+       imo =  fPyjets->K[2][i];
+       if (imo > 1000) {
+           fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
+       }
+    }
+//     this->Pylist(1);
+} // end quench
+
+
+void AliPythia6::Pyquen(Double_t a, Int_t ibf, Double_t b)
+{
+    // Igor Lokthine's quenching routine
+    // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
+
+    pyquen(a, ibf, b);
+}
+
+void AliPythia6::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
+{
+    // Set the parameters for the PYQUEN package.
+    // See comments in PyquenCommon.h
+    
+    
+    PYQPAR.t0    = t0;
+    PYQPAR.tau0  = tau0;
+    PYQPAR.nf    = nf;
+    PYQPAR.iengl = iengl;
+    PYQPAR.iangl = iangl;
+}
+
+void  AliPythia6::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
+{
+//
+// Load event into Pythia Common Block
+//
+
+    Int_t npart = stack -> GetNprimary();
+    Int_t n0 = 0;
+    
+    if (!flag) {
+       GetPyjets()->N = npart;
+    } else {
+       n0 = GetPyjets()->N;
+       GetPyjets()->N = n0 + npart;
+    }
+    
+    
+    for (Int_t part = 0; part < npart; part++) {
+       TParticle *mPart = stack->Particle(part);
+       
+       Int_t kf     =  mPart->GetPdgCode();
+       Int_t ks     =  mPart->GetStatusCode();
+       Int_t idf    =  mPart->GetFirstDaughter();
+       Int_t idl    =  mPart->GetLastDaughter();
+       
+       if (reHadr) {
+           if (ks == 11 || ks == 12) {
+               ks  -= 10;
+               idf  = -1;
+               idl  = -1;
+           }
+       }
+       
+       Float_t px = mPart->Px();
+       Float_t py = mPart->Py();
+       Float_t pz = mPart->Pz();
+       Float_t e  = mPart->Energy();
+       Float_t m  = mPart->GetCalcMass();
+       
+       
+       (GetPyjets())->P[0][part+n0] = px;
+       (GetPyjets())->P[1][part+n0] = py;
+       (GetPyjets())->P[2][part+n0] = pz;
+       (GetPyjets())->P[3][part+n0] = e;
+       (GetPyjets())->P[4][part+n0] = m;
+       
+       (GetPyjets())->K[1][part+n0] = kf;
+       (GetPyjets())->K[0][part+n0] = ks;
+       (GetPyjets())->K[3][part+n0] = idf + 1;
+       (GetPyjets())->K[4][part+n0] = idl + 1;
+       (GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
+    }
+}
+
+
+void AliPythia6::Pyevnw()
+{
+    // New multiple interaction scenario
+    pyevnw();
+}
+
+void AliPythia6::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
+{
+    // Return event specific quenching parameters
+    xp = fXJet;
+    yp = fYJet;
+    for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
+
+}
+
+void AliPythia6::ConfigHeavyFlavor()
+{
+    //
+    // Default configuration for Heavy Flavor production
+    //
+    // All QCD processes
+    //
+    SetMSEL(1);
+    
+    // No multiple interactions
+    SetMSTP(81,0);
+    SetPARP(81, 0.);
+    SetPARP(82, 0.);    
+    // Initial/final parton shower on (Pythia default)
+    SetMSTP(61,1);
+    SetMSTP(71,1);
+    
+    // 2nd order alpha_s
+    SetMSTP(2,2);
+    
+    // QCD scales
+    SetMSTP(32,2);
+    SetPARP(34,1.0);
+}
+
+void AliPythia6::AtlasTuning()
+{
+    //
+    // Configuration for the ATLAS tuning
+       SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ5L));
+       SetMSTP(81,1);             // Multiple Interactions ON
+       SetMSTP(82,4);             // Double Gaussian Model
+       SetPARP(81,1.9);           // Min. pt for multiple interactions (default in 6.2-14) 
+       SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
+       SetPARP(89,1000.);         // [GeV]   Ref. energy
+       SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
+       SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
+       SetPARP(84,0.5);           // Core radius
+       SetPARP(85,0.33);          // Regulates gluon prod. mechanism
+       SetPARP(86,0.66);          // Regulates gluon prod. mechanism
+       SetPARP(67,1);             // Regulates Initial State Radiation
+}
+
+void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
+{
+    // Set the pt hard range
+    SetCKIN(3, ptmin);
+    SetCKIN(4, ptmax);
+}
+
+void AliPythia6::SetYHardRange(Float_t ymin, Float_t ymax)
+{
+    // Set the y hard range
+    SetCKIN(7, ymin);
+    SetCKIN(8, ymax);
+}
+
+
+void AliPythia6::SetFragmentation(Int_t flag)
+{
+    // Switch fragmentation on/off
+    SetMSTP(111, flag);
+}
+
+void AliPythia6::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
+{
+//  initial state radiation    
+    SetMSTP(61, flag1);
+//  final state radiation
+    SetMSTP(71, flag2);
+}
+
+void AliPythia6::SetIntrinsicKt(Float_t kt)
+{
+    // Set the inreinsic kt
+    if (kt > 0.) {
+       SetMSTP(91,1);
+       SetPARP(91,kt); 
+       SetPARP(93, 4. * kt);
+    } else {
+       SetMSTP(91,0);
+    }
+}
+
+void AliPythia6::SwitchHFOff()
+{
+    // Switch off heavy flavor
+    // Maximum number of quark flavours used in pdf 
+    SetMSTP(58, 3);
+    // Maximum number of flavors that can be used in showers
+    SetMSTJ(45, 3);    
+}
+
+void AliPythia6::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+                                      Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
+{
+// Set pycell parameters
+    SetPARU(51, etamax);
+    SetMSTU(51, neta);
+    SetMSTU(52, nphi);
+    SetPARU(58, thresh);
+    SetPARU(52, etseed);
+    SetPARU(53, minet);
+    SetPARU(54, r);
+    SetMSTU(54,  2);
+}
+
+void AliPythia6::ModifiedSplitting()
+{
+    // Modified splitting probability as a model for quenching
+    SetPARJ(200, 0.8);
+    SetMSTJ(41, 1);  // QCD radiation only
+    SetMSTJ(42, 2);  // angular ordering
+    SetMSTJ(44, 2);  // option to run alpha_s
+    SetMSTJ(47, 0);  // No correction back to hard scattering element
+    SetMSTJ(50, 0);  // No coherence in first branching
+    SetPARJ(82, 1.); // Cut off for parton showers
+}
+
+void AliPythia6::SwitchHadronisationOff()
+{
+    // Switch off hadronisarion
+    SetMSTJ(1, 0);
+}
+
+void AliPythia6::SwitchHadronisationOn()
+{
+    // Switch on hadronisarion
+    SetMSTJ(1, 1);
+}
+
+
+void AliPythia6::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
+{
+    // Get x1, x2 and Q for this event
+    
+    q  = GetVINT(51);
+    x1 = GetVINT(41);
+    x2 = GetVINT(42);
+}
+
+Float_t AliPythia6::GetXSection()
+{
+    // Get the total cross-section
+    return (GetPARI(1));
+}
+
+Float_t AliPythia6::GetPtHard()
+{
+    // Get the pT hard for this event
+    return GetVINT(47);
+}
+
+Int_t AliPythia6::ProcessCode()
+{
+    // Get the subprocess code
+    return GetMSTI(1);
+}
+
+void AliPythia6::PrintStatistics()
+{
+    // End of run statistics
+    Pystat(1);
+}
+
+void AliPythia6::EventListing()
+{
+    // End of run statistics
+    Pylist(2);
+}
+
+AliPythia6& AliPythia6::operator=(const  AliPythia6& rhs)
+{
+// Assignment operator
+    rhs.Copy(*this);
+    return *this;
+}
+
+ void AliPythia6::Copy(TObject&) const
+{
+    //
+    // Copy 
+    //
+    Fatal("Copy","Not implemented!\n");
+}
diff --git a/PYTHIA6/AliPythia6.h b/PYTHIA6/AliPythia6.h
new file mode 100644 (file)
index 0000000..8ea59d4
--- /dev/null
@@ -0,0 +1,108 @@
+#ifndef ALIPYTHIA6_H
+#define ALIPYTHIA6_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliPythia.h,v 1.22 2007/10/09 08:43:24 morsch Exp $ */
+
+#include <TPythia6.h>
+#include "AliPythiaBase.h"
+
+
+class AliFastGlauber;
+class AliQuenchingWeights;
+class AliStack;
+
+class AliPythia6 : public TPythia6, public AliPythiaBase
+{
+
+ public:
+    virtual ~AliPythia6(){;}
+    virtual Int_t Version() {return (6);}
+    // convert to compressed code and print result (for debugging only)
+    virtual Int_t CheckedLuComp(Int_t kf);
+    // Pythia initialisation for selected processes
+    virtual void ProcInit
+       (Process_t process, Float_t energy, StrucFunc_t strucfunc);
+    virtual void  GenerateEvent()   {Pyevnt();}
+    virtual void  GenerateMIEvent() {Pyevnw();}
+    virtual void  HadronizeEvent()  {Pyexec();}
+    virtual Int_t GetNumberOfParticles() {return GetN();}
+    virtual void  SetNumberOfParticles(Int_t i) {SetN(i);}
+    virtual void  EditEventList(Int_t i) {Pyedit(i);}
+    virtual void  PrintStatistics();
+    virtual void  EventListing();    
+    virtual Int_t GetParticles(TClonesArray *particles) {return ImportParticles(particles, "All");}
+    // Treat protons as inside nuclei
+    virtual void  SetNuclei(Int_t a1, Int_t a2);
+    // Print particle properties
+    virtual void PrintParticles();
+    // Reset the decay table
+    virtual void ResetDecayTable();
+    //
+    // Common Physics Configuration
+    virtual void SetPtHardRange(Float_t ptmin, Float_t ptmax);
+    virtual void SetYHardRange(Float_t ymin, Float_t ymax);
+    virtual void SetFragmentation(Int_t flag);
+    virtual void SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2);
+    virtual void SetIntrinsicKt(Float_t kt);
+    virtual void SwitchHFOff();
+    virtual void SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+                                    Float_t thresh, Float_t etseed, Float_t minet, Float_t r);
+    virtual void ModifiedSplitting();
+    virtual void SwitchHadronisationOff();
+    virtual void SwitchHadronisationOn();
+    //
+    // Common Getters
+    virtual void    GetXandQ(Float_t& x1, Float_t& x2, Float_t& q);
+    virtual Float_t GetXSection();
+    virtual Int_t   ProcessCode();
+    virtual Float_t GetPtHard();
+    //
+    //
+    virtual void SetDecayTable();
+    virtual void Pyevnw();
+    virtual void Pycell(Int_t& nclus);
+    virtual void Pyclus(Int_t& nclus);
+    virtual void GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e);
+    virtual void LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr);
+    virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax);
+    virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez);
+    virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k, Int_t iECMethod, Float_t zmax = 0.97, Int_t ngmax = 30);
+    virtual void SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl);
+    virtual void Pyquen(Double_t a, Int_t ibf, Double_t b);
+    virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]);
+    // return instance of the singleton
+    static  AliPythia6* Instance();
+    virtual void Quench();
+    // Assignment Operator
+    AliPythia6 & operator=(const AliPythia6 & rhs);
+    void Copy(TObject&) const;
+ protected:
+    Process_t             fProcess;           // Process type
+    Float_t               fEcms;              // Centre of mass energy
+    StrucFunc_t           fStrucFunc;         // Structure function
+    Int_t                 fDefMDCY[501];      //  ! Default decay switches per particle
+    Int_t                 fDefMDME[2001];     //  ! Default decay switches per mode
+    Double_t              fZQuench[4];        //  ! Quenching fractions for this even
+    Double_t              fXJet;              //  ! Jet production point X
+    Double_t              fYJet;              //  ! Jet production point Y
+    Int_t                 fNGmax;             //    Maximum number of radiated gluons in quenching
+    Float_t               fZmax;              //    Maximum energy loss in quenching
+    AliFastGlauber*       fGlauber;           //  ! The Glauber model
+    AliQuenchingWeights*  fQuenchingWeights;  //  ! The Quenching Weights model
+    static AliPythia6*    fgAliPythia;        //  Pointer to single instance
+ private: 
+    AliPythia6();
+    AliPythia6(const AliPythia6& pythia);
+    void ConfigHeavyFlavor();
+    void AtlasTuning();
+    ClassDef(AliPythia6,1) //ALICE UI to PYTHIA
+};
+
+#endif
+
+
+
+
+
index 638298c..c71c3b2 100644 (file)
@@ -5,8 +5,10 @@
 #pragma link off all functions;
  
 #pragma link C++ enum  Process_t;
+#pragma link C++ class AliPythiaBase+;
 #pragma link C++ class AliPythia+;
+#pragma link C++ class AliPythia6+;
 #pragma link C++ class AliGenPythia+;
+#pragma link C++ class AliGenPythiaPlus+;
 #pragma link C++ class AliDecayerPythia+;
-#pragma link C++ class AliGenPythiaJets+;
 #endif
diff --git a/PYTHIA6/AliPythiaBase.cxx b/PYTHIA6/AliPythiaBase.cxx
new file mode 100644 (file)
index 0000000..898d195
--- /dev/null
@@ -0,0 +1,29 @@
+
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+#include "AliPythiaBase.h"
+
+ClassImp(AliPythiaBase)
+
+AliPythiaBase::AliPythiaBase()
+{
+    // Default constructor
+}
+
+
+
diff --git a/PYTHIA6/AliPythiaBase.h b/PYTHIA6/AliPythiaBase.h
new file mode 100644 (file)
index 0000000..6b79896
--- /dev/null
@@ -0,0 +1,86 @@
+#ifndef ALIPYTHIABASE_H
+#define ALIPYTHIABASE_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+#include "AliRndm.h"
+#include <TObject.h>
+#include "AliStructFuncType.h"
+#include "PythiaProcesses.h"
+
+class AliFastGlauber;
+class AliQuenchingWeights;
+class AliStack;
+class TClonesArray;
+
+class AliPythiaBase : public AliRndm 
+{
+
+ public:
+    AliPythiaBase();
+    virtual ~AliPythiaBase(){;}
+    virtual Int_t Version()                                                                                         = 0;
+    // Convert to compressed code and print result (for debugging only)
+    virtual Int_t CheckedLuComp(Int_t kf)                                                                           = 0;
+    // Pythia initialisation for selected processes
+    virtual void  ProcInit (Process_t process, Float_t energy, StrucFunc_t strucfunc)                               = 0;
+    virtual void  GenerateEvent()                                                                                   = 0;
+    virtual void  GenerateMIEvent()                                                                                 = 0;
+    virtual Int_t GetNumberOfParticles()                                                                            = 0;
+    virtual void  SetNumberOfParticles(Int_t i)                                                                     = 0;
+    virtual void  EditEventList(Int_t i)                                                                            = 0;
+    virtual void  HadronizeEvent()                                                                                  = 0;
+    virtual Int_t GetParticles(TClonesArray *particles)                                                             = 0;
+    virtual void  PrintStatistics()                                                                                 = 0;
+    virtual void  EventListing()                                                                                    = 0;    
+    // Treat protons as inside nuclei
+    virtual void  SetNuclei(Int_t a1, Int_t a2)                                                                     = 0;
+    // Print particle properties
+    virtual void PrintParticles()                                                                                   = 0;
+    // Reset the decay table
+    virtual void ResetDecayTable()                                                                                  = 0;
+    //
+    // Common Physics Configuration
+    virtual void SetPtHardRange(Float_t ptmin, Float_t ptmax)                                                       = 0;
+    virtual void SetYHardRange(Float_t ymin, Float_t ymax)                                                          = 0;
+    virtual void SetFragmentation(Int_t flag)                                                                       = 0;
+    virtual void SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)                                         = 0;
+    virtual void SetIntrinsicKt(Float_t kt)                                                                         = 0;
+    virtual void SwitchHFOff()                                                                                      = 0;
+    virtual void SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+                                    Float_t thresh, Float_t etseed, Float_t minet, Float_t r)                      = 0;
+    virtual void ModifiedSplitting()                                                                                = 0;
+    virtual void SwitchHadronisationOff()                                                                           = 0;
+    virtual void SwitchHadronisationOn()                                                                            = 0;
+    //
+    // Common Getters
+    virtual void    GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)                                                  = 0;
+    virtual Float_t GetXSection()                                                                                   = 0;
+    virtual Int_t   ProcessCode()                                                                                   = 0;
+    virtual Float_t GetPtHard()                                                                                     = 0;
+    //
+    //
+    virtual void SetDecayTable()                                                                                    = 0;
+    virtual void Pycell(Int_t& nclus)                                                                               = 0;
+    virtual void Pyclus(Int_t& nclus)                                                                               = 0;
+    virtual void GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)                                 = 0;
+    virtual void LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)                                               = 0;
+    virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)                                                        = 0;
+    virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez) = 0;
+    virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)   = 0;
+    virtual void Pyquen(Double_t a, Int_t ibf, Double_t b)                                                          = 0;
+    virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])                                  = 0;
+    // return instance of the singleton
+    virtual void Quench()                                                                                           = 0;
+    virtual void ConfigHeavyFlavor()                                                                                = 0;
+    virtual void AtlasTuning()                                                                                      = 0;
+    ClassDef(AliPythiaBase, 1) //ALICE UI to PYTHIA
+};
+
+#endif
+
+
+
+
+
index 168d322..a2d4497 100644 (file)
@@ -1,18 +1,9 @@
-#-*- Mode: Makefile -*-
-
-SRCS= AliPythia.cxx AliGenPythia.cxx AliDecayerPythia.cxx \
-AliGenPythiaJets.cxx
+SRCS= AliPythiaBase.cxx AliPythia.cxx AliPythia6.cxx AliGenPythia.cxx AliGenPythiaPlus.cxx AliDecayerPythia.cxx
 
 HDRS= $(SRCS:.cxx=.h) 
 
 DHDR:=AliPythia6LinkDef.h
 
-EXPORT:=AliPythia.h AliDecayerPythia.h
+EXPORT:=AliPythiaBase.h AliPythia.h AliDecayerPythia.h
 
 EINCLUDE:=FASTSIM
-
-ifeq (win32gcc,$(ALICE_TARGET))
-PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) \
-                         -lEVGEN -lSTEER -lpythia6 -lSTEERBase -lFASTSIM \
-                         -L$(shell root-config --libdir) -lVMC -lEG -lEGPythia6
-endif