X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PYTHIA6%2FAliGenPythia.cxx;h=c73145bbc05de3308010f2b7e10fb1f6d245d1eb;hb=31d94b145092dc24f6c2c4fed82a5a324be08e03;hp=c7cd57f958ac92f965f1c5abbef5729c002a1c33;hpb=adf4d898e02daf4e16db15bc3efb5186d8be83c5;p=u%2Fmrichter%2FAliRoot.git

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index c7cd57f958a..c73145bbc05 100644
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -13,205 +13,7 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.1  2003/03/15 15:00:48  morsch
-Classed imported from EVGEN.
-
-Revision 1.69  2003/01/14 10:50:19  alibrary
-Cleanup of STEER coding conventions
-
-Revision 1.68  2002/12/11 09:16:16  morsch
-Use GetJets to fill header.
-
-Revision 1.67  2002/12/09 15:24:09  morsch
-Same trigger routine can use Pycell or Pyclus.
-
-Revision 1.66  2002/12/09 08:22:56  morsch
-UA1 jet finder (Pycell) for software triggering added.
-
-Revision 1.65  2002/11/19 08:57:10  morsch
-Configuration of pt-kick added.
-
-Revision 1.64  2002/11/15 00:43:06  morsch
-Changes for kPyJets
-- initial and final state g-radiation + pt-kick default
-- trigger based on parton clusters (using pyclus)
-- trigger jets are stored in header.
-
-Revision 1.63  2002/10/14 14:55:35  hristov
-Merging the VirtualMC branch to the main development branch (HEAD)
-
-Revision 1.52.4.4  2002/10/10 16:40:08  hristov
-Updating VirtualMC to v3-09-02
-
-Revision 1.62  2002/09/24 10:00:01  morsch
-CheckTrigger() corrected.
-
-Revision 1.61  2002/07/30 09:52:38  morsch
-Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
-
-Revision 1.60  2002/07/19 14:49:03  morsch
-Typo corrected.
-
-Revision 1.59  2002/07/19 14:35:36  morsch
-Count total number of trials. Print mean Q, x1, x2.
-
-Revision 1.58  2002/07/17 10:04:09  morsch
-SetYHard method added.
-
-Revision 1.57  2002/05/22 13:22:53  morsch
-Process kPyMbNonDiffr added.
-
-Revision 1.56  2002/04/26 10:30:01  morsch
-Option kPyBeautyPbMNR added. (N. Carrer).
-
-Revision 1.55  2002/04/17 10:23:56  morsch
-Coding Rule violations corrected.
-
-Revision 1.54  2002/03/28 11:49:10  morsch
-Pass status code in SetTrack.
-
-Revision 1.53  2002/03/25 14:51:13  morsch
-New stack-fill and count options introduced (N. Carrer).
-
-Revision 1.51  2002/03/06 08:46:57  morsch
-- Loop until np-1
-- delete dyn. alloc. arrays (N. Carrer)
-
-Revision 1.50  2002/03/03 13:48:50  morsch
-Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
-NLO calculations (Nicola Carrer).
-
-Revision 1.49  2002/02/08 16:50:50  morsch
-Add name and title in constructor.
-
-Revision 1.48  2001/12/20 11:44:28  morsch
-Add kinematic bias for direct gamma production.
-
-Revision 1.47  2001/12/19 14:45:00  morsch
-Store number of trials in header.
-
-Revision 1.46  2001/12/19 10:36:19  morsch
-Add possibility if jet kinematic biasing.
-
-Revision 1.45  2001/11/28 08:06:52  morsch
-Use fMaxLifeTime parameter.
-
-Revision 1.44  2001/11/27 13:13:07  morsch
-Maximum lifetime for long-lived particles to be put on the stack is parameter.
-It can be set via SetMaximumLifetime(..).
-
-Revision 1.43  2001/10/21 18:35:56  hristov
-Several pointers were set to zero in the default constructors to avoid memory management problems
-
-Revision 1.42  2001/10/15 08:21:55  morsch
-Vertex truncation settings moved to AliGenMC.
-
-Revision 1.41  2001/10/08 08:45:42  morsch
-Possibility of vertex cut added.
-
-Revision 1.40  2001/09/25 11:30:23  morsch
-Pass event vertex to header.
-
-Revision 1.39  2001/07/27 17:09:36  morsch
-Use local SetTrack, KeepTrack and SetHighWaterMark methods
-to delegate either to local stack or to stack owned by AliRun.
-(Piotr Skowronski, A.M.)
-
-Revision 1.38  2001/07/13 10:58:54  morsch
-- Some coded moved to AliGenMC
-- Improved handling of secondary vertices.
-
-Revision 1.37  2001/06/28 11:17:28  morsch
-SetEventListRange setter added. Events in specified range are listed for
-debugging. (Yuri Kharlov)
-
-Revision 1.36  2001/03/30 07:05:49  morsch
-Final print-out in finish run.
-Write parton system for jet-production (preliminary solution).
-
-Revision 1.35  2001/03/09 13:03:40  morsch
-Process_t and Struc_Func_t moved to AliPythia.h
-
-Revision 1.34  2001/02/14 15:50:40  hristov
-The last particle in event marked using SetHighWaterMark
-
-Revision 1.33  2001/01/30 09:23:12  hristov
-Streamers removed (R.Brun)
-
-Revision 1.32  2001/01/26 19:55:51  hristov
-Major upgrade of AliRoot code
-
-Revision 1.31  2001/01/17 10:54:31  hristov
-Better protection against FPE
-
-Revision 1.30  2000/12/18 08:55:35  morsch
-Make AliPythia dependent generartors work with new scheme of random number generation
-
-Revision 1.29  2000/12/04 11:22:03  morsch
-Init of sRandom as in 1.15
-
-Revision 1.28  2000/12/02 11:41:39  morsch
-Use SetRandom() to initialize random number generator in constructor.
-
-Revision 1.27  2000/11/30 20:29:02  morsch
-Initialise static variable sRandom in constructor: sRandom = fRandom;
-
-Revision 1.26  2000/11/30 07:12:50  alibrary
-Introducing new Rndm and QA classes
-
-Revision 1.25  2000/10/18 19:11:27  hristov
-Division by zero fixed
-
-Revision 1.24  2000/09/18 10:41:35  morsch
-Add possibility to use nuclear structure functions from PDF library V8.
-
-Revision 1.23  2000/09/14 14:05:40  morsch
-dito
-
-Revision 1.22  2000/09/14 14:02:22  morsch
-- Correct conversion from mm to cm when passing particle vertex to MC.
-- Correct handling of fForceDecay == all.
-
-Revision 1.21  2000/09/12 14:14:55  morsch
-Call fDecayer->ForceDecay() at the beginning of Generate().
-
-Revision 1.20  2000/09/06 14:29:33  morsch
-Use AliPythia for event generation an AliDecayPythia for decays.
-Correct handling of "nodecay" option
-
-Revision 1.19  2000/07/11 18:24:56  fca
-Coding convention corrections + few minor bug fixes
-
-Revision 1.18  2000/06/30 12:40:34  morsch
-Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
-Geant units (cm).
-
-Revision 1.17  2000/06/09 20:34:07  morsch
-All coding rule violations except RS3 corrected
-
-Revision 1.16  2000/05/15 15:04:20  morsch
-The full event is written for fNtrack = -1
-Coding rule violations corrected.
-
-Revision 1.15  2000/04/26 10:14:24  morsch
-Particles array has one entry more than pythia particle list. Upper bound of
-particle loop changed to np-1 (R. Guernane, AM)
-
-Revision 1.14  2000/04/05 08:36:13  morsch
-Check status code of particles in Pythia event
-to avoid double counting as partonic state and final state particle.
-
-Revision 1.13  1999/11/09 07:38:48  fca
-Changes for compatibility with version 2.23 of ROOT
-
-Revision 1.12  1999/11/03 17:43:20  fca
-New version from G.Martinez & A.Morsch
-
-Revision 1.11  1999/09/29 09:24:14  fca
-Introduction of the Copyright and cvs Log
-*/
+/* $Id$ */
 
 //
 // Generator using the TPythia interface (via AliPythia)
@@ -232,11 +34,16 @@ Introduction of the Copyright and cvs Log
 #include "AliConst.h"
 #include "AliDecayerPythia.h"
 #include "AliGenPythia.h"
+#include "AliHeader.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"
+#include "AliPythiaRndm.h"
 #include "AliRun.h"
+#include "AliStack.h"
+#include "AliRunLoader.h"
+#include "AliMC.h"
 
- ClassImp(AliGenPythia)
+ClassImp(AliGenPythia)
 
 AliGenPythia::AliGenPythia()
                  :AliGenMC()
@@ -244,6 +51,8 @@ AliGenPythia::AliGenPythia()
 // Default Constructor
   fParticles = 0;
   fPythia    = 0;
+  fHeader = 0;
+  fReadFromFile = 0;
   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();
@@ -252,6 +61,11 @@ AliGenPythia::AliGenPythia()
   SetGammaPhiRange();
   SetGammaEtaRange();
   SetPtKick();
+  SetQuench();
+  SetHadronisation();  
+  fSetNuclei = kFALSE;
+  if (!AliPythiaRndm::GetPythiaRandom()) 
+    AliPythiaRndm::SetPythiaRandom(GetRandom());
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
@@ -264,8 +78,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fName = "Pythia";
     fTitle= "Particle Generator using PYTHIA";
     fXsection  = 0.;
-    fNucA1=0;
-    fNucA2=0;
+    fReadFromFile = 0;
     SetProcess();
     SetStrucFunc();
     SetForceDecay();
@@ -274,11 +87,12 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetEnergyCMS();
     fDecayer = new AliDecayerPythia();
     // Set random number generator 
-    sRandom=fRandom;
+    if (!AliPythiaRndm::GetPythiaRandom()) 
+      AliPythiaRndm::SetPythiaRandom(GetRandom());
     fFlavorSelect   = 0;
     // Produced particles  
     fParticles = new TClonesArray("TParticle",1000);
-    fEventVertex.Set(3);
+    fHeader = 0;
     SetEventListRange();
     SetJetPhiRange();
     SetJetEtaRange();
@@ -286,6 +100,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetGammaPhiRange();
     SetGammaEtaRange();
     SetJetReconstructionMode();
+    SetQuench();
+    SetHadronisation();
     SetPtKick();
     // Options determining what to keep in the stack (Heavy flavour generation)
     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
@@ -294,9 +110,13 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fFragmentation = kTRUE;
     // Default counting mode
     fCountMode = kCountAll;
+    // Pycel
+    SetPycellParameters();
+    fSetNuclei = kFALSE;
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
+    :AliGenMC(Pythia)
 {
 // copy constructor
     Pythia.Copy(*this);
@@ -307,6 +127,21 @@ AliGenPythia::~AliGenPythia()
 // Destructor
 }
 
+void AliGenPythia::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 AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 {
   // Set a range of event numbers, for which a table
@@ -321,7 +156,8 @@ void AliGenPythia::Init()
 // Initialisation
     
     SetMC(AliPythia::Instance());
-    fPythia=(AliPythia*) fgMCEvGen;
+    fPythia=(AliPythia*) fMCEvGen;
+    
 //
     fParentWeight=1./Float_t(fNpart);
 //
@@ -335,7 +171,7 @@ void AliGenPythia::Init()
     fPythia->SetCKIN(7,fYHardMin);
     fPythia->SetCKIN(8,fYHardMax);
     
-    if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
+    if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);  
     // Fragmentation?
     if (fFragmentation) {
       fPythia->SetMSTP(111,1);
@@ -356,7 +192,16 @@ void AliGenPythia::Init()
 	fPythia->SetMSTP(91,0);
     }
 
- //   fPythia->SetMSTJ(1,2);
+
+    if (fReadFromFile) {
+	fRL  =  AliRunLoader::Open(fFileName, "Partons");
+	fRL->LoadKinematics();
+	fRL->LoadHeader();
+    } else {
+	fRL = 0x0;
+    }
+    
+    
  //
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 
@@ -366,6 +211,8 @@ void AliGenPythia::Init()
     case kPyCharm:
     case kPyCharmUnforced:
     case kPyCharmPbPbMNR:
+    case kPyCharmppMNR:
+    case kPyCharmpPbMNR:
 	fParentSelect[0] =   411;
 	fParentSelect[1] =   421;
 	fParentSelect[2] =   431;
@@ -378,6 +225,12 @@ void AliGenPythia::Init()
 	fParentSelect[0] =   421;
 	fFlavorSelect    =   4;	
 	break;
+    case kPyDPlusPbPbMNR:
+    case kPyDPluspPbMNR:
+    case kPyDPlusppMNR:
+	fParentSelect[0] =   411;
+	fFlavorSelect    =   4;	
+	break;
     case kPyBeauty:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
@@ -412,6 +265,23 @@ void AliGenPythia::Init()
 	break;
     }
 //
+//
+//  JetFinder for Trigger
+//
+//  Configure detector (EMCAL like)
+//
+	fPythia->SetPARU(51, fPycellEtaMax);
+	fPythia->SetMSTU(51, fPycellNEta);
+	fPythia->SetMSTU(52, fPycellNPhi);
+//
+//  Configure Jet Finder
+//  
+	fPythia->SetPARU(58,  fPycellThreshold);
+	fPythia->SetPARU(52,  fPycellEtSeed);
+	fPythia->SetPARU(53,  fPycellMinEtJet);
+	fPythia->SetPARU(54,  fPycellMaxRadius);
+	fPythia->SetMSTU(54,  2);
+//
 //  This counts the total number of calls to Pyevnt() per run.
     fTrialsRun = 0;
     fQ         = 0.;
@@ -419,17 +289,32 @@ void AliGenPythia::Init()
     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);
+    }
+    
 }
 
 void AliGenPythia::Generate()
 {
 // Generate one event
+    
     fDecayer->ForceDecay();
 
     Float_t polar[3]   =   {0,0,0};
     Float_t origin[3]  =   {0,0,0};
-    Float_t p[3];
+    Float_t p[4];
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
 //
@@ -437,46 +322,76 @@ void AliGenPythia::Generate()
     Int_t jev=0;
     Int_t j, kf;
     fTrials=0;
+    
 
     //  Set collision vertex position 
-    if(fVertexSmear==kPerEvent) {
-	fPythia->SetMSTP(151,1);
-	for (j=0;j<3;j++) {
-	    fPythia->SetPARP(151+j, fOsigma[j]*10.);
-	}
-    } else if (fVertexSmear==kPerTrack) {
-	fPythia->SetMSTP(151,0);
-    }
+    if (fVertexSmear == kPerEvent) Vertex();
+    
 //  event loop    
     while(1)
     {
-	fPythia->Pyevnt();
-	if (gAlice->GetEvNumber()>=fDebugEventFirst &&
-	    gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
-	fTrials++;
+//
+// Produce event
+//
+//
+// Switch hadronisation off
+//
+	fPythia->SetMSTJ(1, 0);
+//
+// Either produce new event or read partons from file
+//	
+	if (!fReadFromFile) {
+	    fPythia->Pyevnt();
+	    fNpartons = fPythia->GetN();
+	} else {
+	    printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
+	    fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
+	    fPythia->SetN(0);
+	    LoadEvent(fRL->Stack(), 0 , 1);
+	    fPythia->Pyedit(21);
+	}
 	
-	fPythia->ImportParticles(fParticles,"All");
+//
+//  Run quenching routine 
+//
+	if (fQuench == 1) {
+	    fPythia->Quench();
+	} else if (fQuench == 2){
+	    fPythia->Pyquen(208., 0, 0.);
+	}
+//
+// Switch hadronisation on
+//
+	fPythia->SetMSTJ(1, 1);
+//
+// .. and perform hadronisation
+//	printf("Calling hadronisation %d\n", fPythia->GetN());
+	fPythia->Pyexec();	
 
+	if (gAlice) {
+	    if (gAlice->GetEvNumber()>=fDebugEventFirst &&
+		gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
+	}
+	
+	fTrials++;
+	fPythia->ImportParticles(fParticles,"All");
+	Boost();
 //
 //
 //
 	Int_t i;
 	
+
 	Int_t np = fParticles->GetEntriesFast();
+	
 	if (np == 0 ) continue;
-// Get event vertex and discard the event if the Z coord. is too big	
-	TParticle *iparticle = (TParticle *) fParticles->At(0);
-	Float_t distz = iparticle->Vz()/10.;
-	if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 //
-	fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
-	fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
-	fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
+	
 //
 	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++) {
+	for (i = 0; i < np; i++) {
 	    pParent[i]   = -1;
 	    pSelected[i] =  0;
 	    trackIt[i]   =  0;
@@ -489,8 +404,8 @@ void AliGenPythia::Generate()
 	    fProcess != kPyDirectGamma &&
 	    fProcess != kPyMbNonDiffr) {
 	    
-	    for (i = 0; i<np; i++) {
-		iparticle = (TParticle *) fParticles->At(i);
+	    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
@@ -518,21 +433,21 @@ void AliGenPythia::Generate()
 		Bool_t flavorOK = kFALSE;
 		Bool_t selectOK = kFALSE;
 		if (fFeedDownOpt) {
-		  if (kfl >= fFlavorSelect) flavorOK = kTRUE;
+		    if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 		} else {
-		  if (kfl > fFlavorSelect) {
-		    nc = -1;
-		    break;
-		  }
-		  if (kfl == fFlavorSelect) flavorOK = kTRUE;
+		    if (kfl > fFlavorSelect) {
+			nc = -1;
+			break;
+		    }
+		    if (kfl == fFlavorSelect) flavorOK = kTRUE;
 		}
 		switch (fStackFillOpt) {
 		case kFlavorSelection:
-		  selectOK = kTRUE;
-		  break;
+		    selectOK = kTRUE;
+		    break;
 		case kParentSelection:
-		  if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
-		  break;
+		    if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
+		    break;
 		}
 		if (flavorOK && selectOK) { 
 //
@@ -588,7 +503,7 @@ void AliGenPythia::Generate()
 // Track final state particle
 		if (ks == 1) trackIt[i] = 1;
 // Track semi-stable particles
-		if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
+		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;
@@ -609,17 +524,24 @@ void AliGenPythia::Generate()
 		    p[0] = iparticle->Px();
 		    p[1] = iparticle->Py();
 		    p[2] = iparticle->Pz();
-		    origin[0] = fOrigin[0]+iparticle->Vx()/10.;
-		    origin[1] = fOrigin[1]+iparticle->Vy()/10.;
-		    origin[2] = fOrigin[2]+iparticle->Vz()/10.;
+		    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     = iparticle->GetFirstMother()-1;
 		    Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
-		    SetTrack(fTrackIt*trackIt[i] ,
-				     iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
+ 
+		    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); 
-		} //  SetTrack loop
+		} //  PushTrack loop
 	    }
   	} else {
 	    nc = GenerateMB();
@@ -660,7 +582,7 @@ void AliGenPythia::Generate()
     } // event loop
     SetHighWaterMark(nt);
 //  adjust weight due to kinematic selection
-    AdjustWeights();
+//    AdjustWeights();
 //  get cross-section
     fXsection=fPythia->GetPARI(1);
 }
@@ -672,13 +594,17 @@ Int_t  AliGenPythia::GenerateMB()
 //
     Int_t i, kf, nt, iparent;
     Int_t nc = 0;
-    Float_t p[3];
+    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 = fParticles->GetEntriesFast();
+
+    
+    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) {
@@ -707,12 +633,34 @@ Int_t  AliGenPythia::GenerateMB()
 	    p[0] = iparticle->Px();
 	    p[1] = iparticle->Py();
 	    p[2] = iparticle->Pz();
-	    origin[0] = fOrigin[0]+iparticle->Vx()/10.;
-	    origin[1] = fOrigin[1]+iparticle->Vy()/10.;
-	    origin[2] = fOrigin[2]+iparticle->Vz()/10.;
+	    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();
-	    SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
-		     tof, kPPrimary, nt, 1., ks);
+
+	    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);
+	    //
+	    // 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;
 	} // select particle
@@ -729,48 +677,55 @@ void AliGenPythia::FinishRun()
 {
 // Print x-section summary
     fPythia->Pystat(1);
-    fQ  /= fNev;
-    fX1 /= fNev;
-    fX2 /= fNev;    
+
+    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 AliGenPythia::AdjustWeights()
 {
 // Adjust the weights after generation of all events
 //
-    TParticle *part;
-    Int_t ntrack=gAlice->GetNtrack();
-    for (Int_t i=0; i<ntrack; i++) {
-	part= gAlice->Particle(i);
-	part->SetWeight(part->GetWeight()*fKineBias);
+    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 AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 {
 // Treat protons as inside nuclei with mass numbers a1 and a2  
-    fNucA1 = a1;
-    fNucA2 = a2;
+
+    fAProjectile = a1;
+    fATarget     = a2;
+    fSetNuclei   = kTRUE;
 }
 
 
 void AliGenPythia::MakeHeader()
 {
 // Builds the event header, to be called after each event
-    AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
+    if (fHeader) delete fHeader;
+    fHeader = new AliGenPythiaEventHeader("Pythia");
 //
 // Event type  
-    ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+    ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
 //
 // Number of trials
-    ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
+    ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
 //
 // Event Vertex 
-    header->SetPrimaryVertex(fEventVertex);
+    fHeader->SetPrimaryVertex(fVertex);
 //
 // Jets that have triggered
     if (fProcess == kPyJets)
@@ -780,11 +735,43 @@ void AliGenPythia::MakeHeader()
 	GetJets(njet, ntrig, jets);
 	
 	for (Int_t i = 0; i < ntrig; i++) {
-	    ((AliGenPythiaEventHeader*) header)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
+	    ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
 							jets[3][i]);
 	}
     }
-    gAlice->SetGenEventHeader(header);
+//
+// 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;
+	
+	fPythia->GetQuenchingParameters(xp, yp, z);
+	
+	((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
+	((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
+    }
+    
+//
+//  Pass header to RunLoader
+//
+    AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(fHeader);   
 }
 	
 
@@ -839,62 +826,71 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 {
 // Assignment operator
+    rhs.Copy(*this);
     return *this;
 }
 
-void  AliGenPythia::LoadEvent()
+void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
 {
 //
 // Load event into Pythia Common Block
 //
- 
-
-    Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries(); 
-   (fPythia->GetPyjets())->N = npart;
 
+    Int_t npart = stack -> GetNprimary();
+    Int_t n0 = 0;
+    
+    if (!flag) {
+	(fPythia->GetPyjets())->N = npart;
+    } else {
+	n0 = (fPythia->GetPyjets())->N;
+	(fPythia->GetPyjets())->N = n0 + npart;
+    }
+    
+    
     for (Int_t part = 0; part < npart; part++) {
-	TParticle *MPart = gAlice->Particle(part);
-	Int_t kf     = MPart->GetPdgCode();
-	Int_t ks     = MPart->GetStatusCode();
+	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 p  = TMath::Sqrt(px * px + py * py + pz * pz);
-	Float_t m  = TMath::Sqrt(e * e - p * p);
+	Float_t m  = MPart->GetCalcMass();
 	
 	
-	(fPythia->GetPyjets())->P[0][part] = px;
-	(fPythia->GetPyjets())->P[1][part] = py;
-	(fPythia->GetPyjets())->P[2][part] = pz;
-	(fPythia->GetPyjets())->P[3][part] = e;
-	(fPythia->GetPyjets())->P[4][part] = m;
+	(fPythia->GetPyjets())->P[0][part+n0] = px;
+	(fPythia->GetPyjets())->P[1][part+n0] = py;
+	(fPythia->GetPyjets())->P[2][part+n0] = pz;
+	(fPythia->GetPyjets())->P[3][part+n0] = e;
+	(fPythia->GetPyjets())->P[4][part+n0] = m;
 	
-	(fPythia->GetPyjets())->K[1][part] = kf;
-	(fPythia->GetPyjets())->K[0][part] = ks;
+	(fPythia->GetPyjets())->K[1][part+n0] = kf;
+	(fPythia->GetPyjets())->K[0][part+n0] = ks;
+	(fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
+	(fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
+	(fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
     }
 }
 
-void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax, 
-			      Int_t& njets, Float_t jets [4][50])
+
+void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
 {
 //
 //  Calls the Pythia jet finding algorithm to find jets in the current event
 //
 //
-//  Configure detector (EMCAL like)
-//
-    fPythia->SetPARU(51,2.);
-    fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
-    fPythia->SetMSTU(52,3 * 144);
-//
-//  Configure Jet Finder
-//  
-    fPythia->SetPARU(58, eCellMin);
-    fPythia->SetPARU(52, eCellSeed);
-    fPythia->SetPARU(53, eMin);
-    fPythia->SetPARU(54, rMax);
-    fPythia->SetMSTU(54, 2);
 //
 //  Save jets
     Int_t n     = fPythia->GetN();
@@ -940,20 +936,7 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 //  Loading jets from common block
 //
     } else {
-//
-//  Configure detector (EMCAL like)
-//
-	fPythia->SetPARU(51,2.);
-	fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
-	fPythia->SetMSTU(52,3 * 144);
-//
-//  Configure Jet Finder
-//  
-	fPythia->SetPARU(58,  0.0);
-	fPythia->SetPARU(52,  4.0);
-	fPythia->SetPARU(53, 10.0);
-	fPythia->SetPARU(54,  1.0);
-	fPythia->SetMSTU(54,  2);
+
 //
 //  Run Jet Finder
 	fPythia->Pycell(nJets);
@@ -966,14 +949,14 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 	Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 	Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 	Float_t pt    = TMath::Sqrt(px * px + py * py);
-	Float_t phi   = TMath::ATan2(py,px);
+	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 && eta < fPhiMaxJet &&
+	    phi > fPhiMinJet && phi < fPhiMaxJet &&
 	    et  > fEtMinJet  && et  < fEtMaxJet     
 	    ) 
 	{
@@ -982,7 +965,7 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 	    jets[2][nJetsTrig] = pz;
 	    jets[3][nJetsTrig] = e;
 	    nJetsTrig++;
-	    
+//	    printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
 	} else {
 //	    printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
 	}
@@ -1031,3 +1014,4 @@ void AliGenPythia::Streamer(TBuffer &R__b)
 }
 #endif
 
+