- Add possibility to read partonic event from file and hadronize.

- Add quenching information to event header.

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 7317beb1585d7a2976d718fd353bb36e49c2c153..15a8dd2afcf66e4d94402b8c29e7255154fff6d1 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -34,6 +34,7 @@
 #include "AliConst.h"
 #include "AliDecayerPythia.h"
 #include "AliGenPythia.h"
 #include "AliConst.h"
 #include "AliDecayerPythia.h"
 #include "AliGenPythia.h"

+#include "AliHeader.h"

 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"
 #include "AliPythiaRndm.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"
 #include "AliPythiaRndm.h"


@@ -60,7 +61,7 @@ AliGenPythia::AliGenPythia()
   SetGammaEtaRange();
   SetPtKick();
   SetQuench();
   SetGammaEtaRange();
   SetPtKick();
   SetQuench();

-  
+  SetHadronisation();  

   fSetNuclei = kFALSE;
   if (!AliPythiaRndm::GetPythiaRandom()) 
     AliPythiaRndm::SetPythiaRandom(GetRandom());
   fSetNuclei = kFALSE;
   if (!AliPythiaRndm::GetPythiaRandom()) 
     AliPythiaRndm::SetPythiaRandom(GetRandom());


@@ -98,6 +99,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetGammaEtaRange();
     SetJetReconstructionMode();
     SetQuench();
     SetGammaEtaRange();
     SetJetReconstructionMode();
     SetQuench();

+    SetHadronisation();

     SetPtKick();
     // Options determining what to keep in the stack (Heavy flavour generation)
     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
     SetPtKick();
     // Options determining what to keep in the stack (Heavy flavour generation)
     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor


@@ -188,7 +190,16 @@ void AliGenPythia::Init()
        fPythia->SetMSTP(91,0);
     }
 
        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);
 
  //
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 


@@ -288,7 +299,7 @@ void AliGenPythia::Init()
     }
 
     if (fQuench) {
     }
 
     if (fQuench) {

-       fPythia->InitQuenching(0., 0.1, 6.e5, 0);
+       fPythia->InitQuenching(0., 0.1, 0.6e6, 0);

     }
     
 }
     }
     
 }


@@ -318,33 +329,43 @@ void AliGenPythia::Generate()
     while(1)
     {
 //
     while(1)
     {
 //

-// If quenching option has been selected switch off fragmentation first
+// Produce event

 //
 //

-       if (fQuench) {
-           fPythia->SetMSTJ(1, 0);
-       }

 //
 //

-// Produce event
+// Switch hadronisation off
+//
+       fPythia->SetMSTJ(1, 0);

 //
 //

-       fPythia->Pyevnt();
-       if (fQuench) {
+// 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);
+       }
+       

 //
 //  Run quenching routine 
 //
 //
 //  Run quenching routine 
 //

-           if (fQuench == 1) {
-               fPythia->Quench();
-           } else {
-               fPythia->Pyquen(208., 0, 0.);
-           }
-            
+       if (fQuench == 1) {
+           fPythia->Quench();
+       } else if (fQuench == 2){
+           fPythia->Pyquen(208., 0, 0.);
+       }

 //
 //

-// Switch fragmentation on
-           fPythia->SetMSTJ(1, 1);
+// Switch hadronisation on

 //
 //

-// .. and perform fragmentation
-           fPythia->Pyexec();
-       }
-       
+       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);
        if (gAlice) {
            if (gAlice->GetEvNumber()>=fDebugEventFirst &&
                gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);


@@ -358,7 +379,9 @@ void AliGenPythia::Generate()
 //
        Int_t i;
        
 //
        Int_t i;
        

+

        Int_t np = fParticles->GetEntriesFast();
        Int_t np = fParticles->GetEntriesFast();

+       

        if (np == 0 ) continue;
 //
        
        if (np == 0 ) continue;
 //
        


@@ -366,7 +389,7 @@ void AliGenPythia::Generate()
        Int_t* pParent   = new Int_t[np];
        Int_t* pSelected = new Int_t[np];
        Int_t* trackIt   = new Int_t[np];
        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;
            pParent[i]   = -1;
            pSelected[i] =  0;
            trackIt[i]   =  0;


@@ -379,7 +402,7 @@ void AliGenPythia::Generate()
            fProcess != kPyDirectGamma &&
            fProcess != kPyMbNonDiffr) {
            
            fProcess != kPyDirectGamma &&
            fProcess != kPyMbNonDiffr) {
            

-           for (i = 0; i<np; i++) {
+           for (i = 0; i < np; i++) {

                TParticle* iparticle = (TParticle *) fParticles->At(i);
                Int_t ks = iparticle->GetStatusCode();
                kf = CheckPDGCode(iparticle->GetPdgCode());
                TParticle* iparticle = (TParticle *) fParticles->At(i);
                Int_t ks = iparticle->GetStatusCode();
                kf = CheckPDGCode(iparticle->GetPdgCode());


@@ -575,7 +598,11 @@ Int_t  AliGenPythia::GenerateMB()
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
     
 //  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) {
     Int_t* pParent = new Int_t[np];
     for (i=0; i< np; i++) pParent[i] = -1;
     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {


@@ -618,6 +645,20 @@ Int_t  AliGenPythia::GenerateMB()
                      origin[0], origin[1], origin[2], tof, 
                      polar[0], polar[1], polar[2],
                      kPPrimary, nt, 1., ks);
                      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
            KeepTrack(nt);
            pParent[i] = nt;
        } // select particle


@@ -694,7 +735,39 @@ void AliGenPythia::MakeHeader()
                                                        jets[3][i]);
        }
     }
                                                        jets[3][i]);
        }
     }

-    if (gAlice) gAlice->SetGenEventHeader(fHeader);
+//
+// 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);   

 }
        
 
 }
        
 


@@ -753,13 +826,12 @@ AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
     return *this;
 }
 
     return *this;
 }
 

-void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag)
+void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)

 {
 //
 // Load event into Pythia Common Block
 //
 {
 //
 // Load event into Pythia Common Block
 //

- 
-//    AliRunLoader* rl = AliRunLoader::GetRunLoader();
+

     Int_t npart = stack -> GetNprimary();
     Int_t n0 = 0;
     
     Int_t npart = stack -> GetNprimary();
     Int_t n0 = 0;
     


@@ -774,8 +846,18 @@ void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag)
     for (Int_t part = 0; part < npart; part++) {
        TParticle *MPart = stack->Particle(part);
        
     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 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 px = MPart->Px();
        Float_t py = MPart->Py();


@@ -792,9 +874,13 @@ void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag)
        
        (fPythia->GetPyjets())->K[1][part+n0] = kf;
        (fPythia->GetPyjets())->K[0][part+n0] = 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(Int_t& njets, Float_t jets [4][50])
 {
 //
 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
 {
 //


@@ -875,7 +961,7 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
            jets[2][nJetsTrig] = pz;
            jets[3][nJetsTrig] = e;
            nJetsTrig++;
            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);
        }
        } else {
 //         printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
        }

diff --git a/PYTHIA6/AliGenPythia.h b/PYTHIA6/AliGenPythia.h
index 3729530e28a87675419997cf6e30e12ae8949ec0..9058258918fc87a7b512b9c15291655dafd04961 100644 (file)
--- a/PYTHIA6/AliGenPythia.h
+++ b/PYTHIA6/AliGenPythia.h
@@ -22,6 +22,8 @@ class AliPythia;
 class TParticle;
 class AliGenPythiaEventHeader;
 class AliStack;
 class TParticle;
 class AliGenPythiaEventHeader;
 class AliStack;

+class AliRunLoader;
+

 class AliGenPythia : public AliGenMC
 {
  public:
 class AliGenPythia : public AliGenMC
 {
  public:


@@ -87,7 +89,8 @@ class AliGenPythia : public AliGenMC
     }
     // Set quenching mode 0 = no, 1 = AM, 2 = IL
     virtual void SetQuench(Int_t flag = 0) {fQuench = flag;}
     }
     // 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;}    

            
     // get cross section of process
     virtual Float_t GetXsection() const {return fXsection;}
            
     // get cross section of process
     virtual Float_t GetXsection() const {return fXsection;}


@@ -98,7 +101,7 @@ class AliGenPythia : public AliGenMC
                             Float_t thresh = 0., Float_t etseed = 4.,
                             Float_t minet = 10., Float_t r = 1.);
     
                             Float_t thresh = 0., Float_t etseed = 4.,
                             Float_t minet = 10., Float_t r = 1.);
     

-    void LoadEvent(AliStack* stack, Int_t flag = 0);
+    void LoadEvent(AliStack* stack, Int_t flag = 0, Int_t reHadr = 0);

     // Getters
     virtual Process_t    GetProcess() {return fProcess;}
     virtual StrucFunc_t  GetStrucFunc() {return fStrucFunc;}
     // Getters
     virtual Process_t    GetProcess() {return fProcess;}
     virtual StrucFunc_t  GetStrucFunc() {return fStrucFunc;}


@@ -147,6 +150,9 @@ class AliGenPythia : public AliGenMC
     Float_t     fYHardMax;          //higher y-hard cut
     Int_t       fGinit;             //initial state gluon radiation
     Int_t       fGfinal;            //final state gluon radiation
     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
     Int_t       fQuench;            //Flag for quenching
     Float_t     fPtKick;            //Transverse momentum kick
     Bool_t      fFullEvent;         //!Write Full event if true


@@ -180,8 +186,10 @@ class AliGenPythia : public AliGenMC
     Bool_t fSetNuclei;              // Flag indicating that SetNuclei has been called
     //
 
     Bool_t fSetNuclei;              // Flag indicating that SetNuclei has been called
     //
 

-    CountMode_t fCountMode;         // Options for counting when the event will be finished.
-    AliGenPythiaEventHeader* fHeader; //! Event header
+    CountMode_t fCountMode;            // Options for counting when the event will be finished.
+    AliGenPythiaEventHeader* fHeader;  //! Event header
+    AliRunLoader*            fRL;      //! Run Loader
+    const Text_t* fFileName;           //Name of file to read from

      
     // fCountMode = kCountAll         --> All particles that end up in the
     //                                    stack are counted
      
     // fCountMode = kCountAll         --> All particles that end up in the
     //                                    stack are counted