First commit.

diff --git a/EVGEN/AliGenPythiaJets.cxx b/EVGEN/AliGenPythiaJets.cxx
new file mode 100644 (file)
index 0000000..70ca2f8
--- /dev/null
+++ b/EVGEN/AliGenPythiaJets.cxx
@@ -0,0 +1,257 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+//
+// 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 "AliGenPythiaJets.h"
+#include "AliRun.h"
+ClassImp(AliGenPythiaJets)
+
+AliGenPythiaJets::AliGenPythiaJets()
+                 :AliGenPythia()
+{
+// Default Constructor
+}
+
+AliGenPythiaJets::AliGenPythiaJets(Int_t npart)
+                 :AliGenPythia(npart)
+{
+    fName = "PythiaJets";
+    fTitle= "Jet Generator using PYTHIA";
+}
+
+AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia)
+{
+// copy constructor
+    Pythia.Copy(*this);
+}
+
+AliGenPythiaJets::~AliGenPythiaJets()
+{
+// Destructor
+}
+
+void AliGenPythiaJets::Init()
+{
+// Initialization
+//
+    printf("AliGenPythiaJets::Init() \n");
+    
+    AliGenPythia::Init();
+    
+    if (fQuench > 0.) {
+        fEtMinJetQ[0]  = fEtMinJet  * fQuench;  
+        fEtMaxJetQ[0]  = fEtMaxJet  * fQuench;  
+        fEtMinJetQ[1]  = fEtMinJet  * (1. - fQuench);  
+        fEtMaxJetQ[1]  = fEtMaxJet  * (1. - fQuench);
+        fPtHardMinQ[0] = fPtHardMin * fQuench;  
+        fPtHardMaxQ[0] = fPtHardMax * fQuench;  
+        fPtHardMinQ[1] = fPtHardMin * (1. - fQuench);  
+        fPtHardMaxQ[1] = fPtHardMax * (1. - fQuench);  
+    }
+}
+
+void AliGenPythiaJets::Generate()
+{
+// Generate one event
+    fDecayer->ForceDecay();
+
+    Float_t polar[3]   =   {0,0,0};
+    Float_t origin[3]  =   {0,0,0};
+    Float_t p[3];
+//  converts from mm/c to s
+    const Float_t kconv=0.001/2.999792458e8;
+//
+    Int_t nt  = 0;
+    Int_t nc  = 0;
+    Int_t jev = 0;
+    Int_t j, kf, iparent;
+    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);
+    }
+//  Event loop    
+    while(1)
+    {
+       if (fQuench > 0.) {
+           fPythia->SetCKIN(3,fPtHardMinQ[jev]);
+           fPythia->SetCKIN(4,fPtHardMaxQ[jev]);
+           fEtMinJet = fEtMinJetQ[jev];
+           fEtMaxJet = fEtMaxJetQ[jev];            
+       }
+       
+       fPythia->Pyevnt();
+       if (gAlice->GetEvNumber()>=fDebugEventFirst &&
+           gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
+       fTrials++;
+        //
+       // Has this jet triggered
+        //
+       if ((fEtMinJet != -1) && ! CheckTrigger()) continue;
+//
+       fPythia->ImportParticles(fParticles,"All");
+       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];
+       for (i=0; i< np; i++) pParent[i] = -1;
+       
+
+       //
+       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 = iparticle->GetFirstMother() - 1;
+               iparent = (ipa > -1) ? pParent[ipa] : -1;
+//
+// Store track information
+//
+               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.;
+               Float_t tof=kconv*iparticle->T();
+               SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
+                        tof, kPPrimary, nt, 1., ks);
+               KeepTrack(nt);
+               pParent[i] = nt;
+           } // select particle
+       } // particle loop 
+       
+       if (pParent) delete[] pParent;
+       printf("\n AliGenPythiaJets: I've put %i particles on the stack \n",nc);
+       if (nc > 0) {
+           jev += 1;
+           if ((fQuench <= 0.) || (fQuench > 0. && jev == 2)) {
+               fKineBias=Float_t(fNpart)/Float_t(fTrials);
+               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+               fNev++;
+               MakeHeader();
+               break;
+           }
+       }
+    }
+    SetHighWaterMark(nt);
+//  Get cross-section
+    fXsection=fPythia->GetPARI(1);
+}
+
+Bool_t AliGenPythiaJets::CheckTrigger()
+{
+// Check the kinematic trigger condition
+//
+    Bool_t   triggered = kFALSE;
+    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) {
+       triggered = kTRUE;
+       Float_t px   = jets[0][0];
+       Float_t py   = jets[1][0];
+       Float_t pz   = jets[2][0];
+       Float_t e    = jets[3][0];
+       Float_t beta = pz/e;
+       Float_t phi  =  TMath::Pi()+TMath::ATan2(-py,-px);
+       TransformEvent(beta, -TMath::Pi()/2. + phi);
+    }
+    return triggered;
+}
+         
+AliGenPythiaJets& AliGenPythiaJets::operator=(const  AliGenPythiaJets& rhs)
+{
+// Assignment operator
+    return *this;
+}
+
+void  AliGenPythiaJets::TransformEvent(Float_t beta, Float_t phi)
+{
+//
+// Perform Lorentz Transformation and Rotation
+//
+    printf("Transforming \n");
+    
+ 
+    Float_t gamma = 1./TMath::Sqrt(1. - beta * beta);
+    Int_t npart = (fPythia->GetPyjets())->N;
+
+    for (Int_t part = 0; part < npart; part++) {
+       Float_t px   =  (fPythia->GetPyjets())->P[0][part];
+       Float_t py   =  (fPythia->GetPyjets())->P[1][part];
+       Float_t pz   =  (fPythia->GetPyjets())->P[2][part];
+       Float_t e    =  (fPythia->GetPyjets())->P[3][part];
+       //
+       // Lorentz Transform
+       //
+       Float_t pzt =  gamma * pz        - gamma * beta * e;
+       Float_t et  = -gamma * beta * pz + gamma        * e;
+       //
+       // Rotation
+       //
+       Float_t pxt =   TMath::Cos(phi) * px +  TMath::Sin(phi) * py;
+       Float_t pyt = - TMath::Sin(phi) * px +  TMath::Cos(phi) * py;
+       //
+       //
+       (fPythia->GetPyjets())->P[0][part] = pxt;
+       (fPythia->GetPyjets())->P[1][part] = pyt;
+       (fPythia->GetPyjets())->P[2][part] = pzt;
+       (fPythia->GetPyjets())->P[3][part] = et;
+    }
+}
+

diff --git a/EVGEN/AliGenPythiaJets.h b/EVGEN/AliGenPythiaJets.h
new file mode 100644 (file)
index 0000000..ae1c6d7
--- /dev/null
+++ b/EVGEN/AliGenPythiaJets.h
@@ -0,0 +1,48 @@
+#ifndef ALIGENPYTHIAJETS_H
+#define ALIGENPYTHIAJETS_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 jets in pp collisions.
+//
+// andreas.morsch@cern.ch
+//
+
+#include "AliGenMC.h"
+#include "AliGenPythia.h"
+
+class TParticle;
+
+class AliGenPythiaJets : public AliGenPythia
+{
+ public:
+    AliGenPythiaJets();
+    AliGenPythiaJets(Int_t npart);
+    AliGenPythiaJets(const AliGenPythiaJets &Pythia);
+    virtual ~AliGenPythiaJets();
+    virtual void    Init();
+    virtual void    Generate();
+    virtual void    TransformEvent(Float_t beta, Float_t gamma);
+    virtual Bool_t  CheckTrigger();
+    virtual void    SetQuenchingFactor(Float_t quench = -1) {fQuench = quench;}
+    
+    // Assignment Operator
+    AliGenPythiaJets & operator=(const AliGenPythiaJets & rhs);
+ protected:
+    Float_t fQuench;                // Quench factor
+    Float_t fEtMinJetQ[2];          // Minimum et of triggered Jet
+    Float_t fEtMaxJetQ[2];          // Maximum et of triggered Jet
+    Float_t fPtHardMinQ[2];         // Lower pT-hard cut 
+    Float_t fPtHardMaxQ[2];         // Higher pT-hard cut
+    ClassDef(AliGenPythiaJets,1)    // AliGenerator Interface to Pythia Jet Production
+};
+#endif
+
+
+
+
+