- Add possibility to read partonic event from file and hadronize.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 1 Jul 2004 07:13:32 +0000 (07:13 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 1 Jul 2004 07:13:32 +0000 (07:13 +0000)
- Add quenching information to event header.

PYTHIA6/AliGenPythia.cxx
PYTHIA6/AliGenPythia.h

index 7317beb1585d7a2976d718fd353bb36e49c2c153..15a8dd2afcf66e4d94402b8c29e7255154fff6d1 100644 (file)
@@ -34,6 +34,7 @@
 #include "AliConst.h"
 #include "AliDecayerPythia.h"
 #include "AliGenPythia.h"
+#include "AliHeader.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"
 #include "AliPythiaRndm.h"
@@ -60,7 +61,7 @@ AliGenPythia::AliGenPythia()
   SetGammaEtaRange();
   SetPtKick();
   SetQuench();
-  
+  SetHadronisation();  
   fSetNuclei = kFALSE;
   if (!AliPythiaRndm::GetPythiaRandom()) 
     AliPythiaRndm::SetPythiaRandom(GetRandom());
@@ -98,6 +99,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetGammaEtaRange();
     SetJetReconstructionMode();
     SetQuench();
+    SetHadronisation();
     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->SetMSTJ(1,2);
+
+    if (fReadFromFile) {
+       fRL  =  AliRunLoader::Open(fFileName, "Partons");
+       fRL->LoadKinematics();
+       fRL->LoadHeader();
+    } else {
+       fRL = 0x0;
+    }
+    
+    
  //
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 
@@ -288,7 +299,7 @@ void AliGenPythia::Init()
     }
 
     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)
     {
 //
-// 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 
 //
-           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);
@@ -358,7 +379,9 @@ void AliGenPythia::Generate()
 //
        Int_t i;
        
+
        Int_t np = fParticles->GetEntriesFast();
+       
        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];
-       for (i=0; i< np; i++) {
+       for (i = 0; i < np; i++) {
            pParent[i]   = -1;
            pSelected[i] =  0;
            trackIt[i]   =  0;
@@ -379,7 +402,7 @@ void AliGenPythia::Generate()
            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());
@@ -575,7 +598,11 @@ Int_t  AliGenPythia::GenerateMB()
 //  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) {
@@ -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);
+           //
+           // 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
@@ -694,7 +735,39 @@ void AliGenPythia::MakeHeader()
                                                        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;
 }
 
-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
 //
-//    AliRunLoader* rl = AliRunLoader::GetRunLoader();
+
     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);
        
-       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();
@@ -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[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])
 {
 //
@@ -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++;
-           
+//         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);
        }
index 3729530e28a87675419997cf6e30e12ae8949ec0..9058258918fc87a7b512b9c15291655dafd04961 100644 (file)
@@ -22,6 +22,8 @@ class AliPythia;
 class TParticle;
 class AliGenPythiaEventHeader;
 class AliStack;
+class AliRunLoader;
+
 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;}
-    
+    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;}
@@ -98,7 +101,7 @@ class AliGenPythia : public AliGenMC
                             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;}
@@ -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
+    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
@@ -180,8 +186,10 @@ class AliGenPythia : public AliGenMC
     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