First commit.

diff --git a/EVGEN/AliGenHerwig.cxx b/EVGEN/AliGenHerwig.cxx
new file mode 100644 (file)
index 0000000..b255018
--- /dev/null
+++ b/EVGEN/AliGenHerwig.cxx
@@ -0,0 +1,311 @@
+/**************************************************************************
+ * 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 Herwig as an external generator
+// The main Herwig options are accessable for the user through this interface.
+// Uses the THerwig implementation of TGenerator.
+
+#include <stream.h>
+#include "AliGenHerwig.h"
+#include "AliRun.h"
+
+#include <TParticle.h>
+#include "THerwig6.h"
+
+
+ ClassImp(AliGenHerwig)
+
+AliGenHerwig::AliGenHerwig()
+{
+// Constructor
+}
+
+AliGenHerwig::AliGenHerwig(Int_t npart)
+{
+    SetBeamMomenta();
+    SetTarget();
+    SetProjectile();
+    SetStrucFunc(kGRVLO98);
+    fKeep=0;
+    fTrigger=0;
+    fDecaysOff=1;
+    fSelectAll=0;
+    fFlavor=0;
+    fPtHardMin=10.;
+    fPtRMS=0.0;
+    fEnSoft=1.0;
+    fMaxPr=1;
+    fMaxErrors=1000;
+//  Set random number
+    if (!sRandom) sRandom=fRandom;
+}
+
+AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
+{
+// copy constructor
+}
+
+
+AliGenHerwig::~AliGenHerwig()
+{
+// Destructor
+}
+
+void AliGenHerwig::Init()
+{
+// Initialisation
+  fTarget.Resize(8);
+  fProjectile.Resize(8);
+  SetMC(new THerwig6());
+  fHerwig=(THerwig6*) fgMCEvGen;
+  // initialize common blocks
+  fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
+  // reset parameters according to user needs
+  InitPDF();
+  fHerwig->SetPTMIN(fPtHardMin);
+  fHerwig->SetPTRMS(fPtRMS);
+  fHerwig->SetMAXPR(fMaxPr);
+  fHerwig->SetMAXER(fMaxErrors);
+  fHerwig->SetENSOF(fEnSoft);
+  // compute parameter dependent constants
+  fHerwig->PrepareRun();
+}
+
+void AliGenHerwig::InitPDF()
+{
+  switch(fStrucFunc)
+    {
+    case kGRVLO:
+      fModPDF=5;
+      fAutPDF="GRV";
+      break;
+    case kGRVHO:
+      fModPDF=6;
+      fAutPDF="GRV";
+      break;
+    case kGRVLO98:
+      fModPDF=12;
+      fAutPDF="GRV";
+      break;
+    case kMRSDminus:
+      fModPDF=31;
+      fAutPDF="MRS";
+      break;
+    case kMRSD0:
+      fModPDF=30;
+      fAutPDF="MRS";
+      break;
+    case kMRSG:
+      fModPDF=41;
+      fAutPDF="MRS";
+      break;
+    case kMRSTcgLO:
+      fModPDF=72;
+      fAutPDF="MRS";
+      break;
+    case kCTEQ4M:
+      fModPDF=34;
+      fAutPDF="CTEQ";
+      break;
+    case kCTEQ5L:
+      fModPDF=46;
+      fAutPDF="CTEQ";
+      break;
+    }
+  fAutPDF.Resize(20);      
+  fHerwig->SetMODPDF(1,fModPDF);
+  fHerwig->SetMODPDF(2,fModPDF);
+  fHerwig->SetAUTPDF(1,fAutPDF);
+  fHerwig->SetAUTPDF(2,fAutPDF);
+}
+
+void AliGenHerwig::Generate()
+{
+  // Generate one event
+
+  Float_t polar[3] =   {0,0,0};
+  Float_t origin[3]=   {0,0,0};
+  Float_t origin0[3]=  {0,0,0};
+  Float_t p[3], random[6];
+  
+  static TClonesArray *particles;
+  //  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, ks, imo;
+  kf=0;
+  
+  if(!particles) particles=new TClonesArray("TParticle",10000);
+  
+  fTrials=0;
+  for (j=0;j<3;j++) origin0[j]=fOrigin[j];
+  if(fVertexSmear==kPerEvent) {
+    Rndm(random,6);
+    for (j=0;j<3;j++) {
+      origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
+       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
+    }
+  }
+
+  while(1)
+    {
+      fHerwig->GenerateEvent();
+      fTrials++;
+      fHerwig->ImportParticles(particles,"All");
+      Int_t np = particles->GetEntriesFast()-1;
+      if (np == 0 ) continue;
+
+      Int_t nc=0;
+
+      Int_t * newPos = new Int_t[np];
+      for (Int_t i = 0; i<np; i++) *(newPos+i)=i;
+      
+      for (Int_t i = 0; i<np; i++) {
+       TParticle *  iparticle       = (TParticle *) particles->At(i);
+       imo = iparticle->GetFirstMother();
+       kf        = iparticle->GetPdgCode();
+       ks        = iparticle->GetStatusCode();
+       if (ks == 1 && kf != 0 && KinematicSelection(iparticle,0)) {
+         nc++;
+         p[0]=iparticle->Px();
+         p[1]=iparticle->Py();
+         p[2]=iparticle->Pz();
+         origin[0]=origin0[0]+iparticle->Vx()/10;
+         origin[1]=origin0[1]+iparticle->Vy()/10;
+         origin[2]=origin0[2]+iparticle->Vz()/10;
+         Float_t tof=kconv*iparticle->T();
+         SetTrack(fTrackIt,-1,kf,p,origin,polar,tof,kPPrimary,nt);
+         KeepTrack(nt);
+         newPos[i]=nt;
+       } // end of if: selection of particle
+      } // end of for: particle loop
+      if (newPos) delete[] newPos;
+      printf("\n I've put %i particles on the stack \n",nc);
+      //      MakeHeader();
+
+      if (nc > 0) {
+       jev+=nc;
+       if (jev >= fNpart || fNpart == -1) {
+         fKineBias=Float_t(fNpart)/Float_t(fTrials);
+         printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+         break;
+       }
+      }
+    }
+    SetHighWaterMark(nt);
+//  adjust weight due to kinematic selection
+    AdjustWeights();
+//  get cross-section
+    fXsection=fHerwig->GetAVWGT();
+}
+
+void AliGenHerwig::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);
+    }
+}
+ 
+
+void AliGenHerwig::KeepFullEvent()
+{
+    fKeep=1;
+}
+
+Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
+{
+//
+// Looks recursively if one of the daughters has been selected
+//
+//    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
+    Int_t imin=-1;
+    Int_t imax=-1;
+    Int_t i;
+    Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
+    Bool_t selected=kFALSE;
+    if (hasDaughters) {
+       imin=iparticle->GetFirstDaughter();
+       imax=iparticle->GetLastDaughter();       
+       for (i=imin; i<= imax; i++){
+           TParticle *  jparticle       = (TParticle *) particles->At(i);      
+           Int_t ip=jparticle->GetPdgCode();
+           if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
+               selected=kTRUE; break;
+           }
+           if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
+       }
+    } else {
+       return kFALSE;
+    }
+
+    return selected;
+}
+
+
+Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
+{
+// Select flavor of particle
+// 0: all
+// 4: charm and beauty
+// 5: beauty
+    if (fFlavor == 0) return kTRUE;
+    
+    Int_t ifl=TMath::Abs(pid/100);
+    if (ifl > 10) ifl/=10;
+    return (fFlavor == ifl);
+}
+
+Bool_t AliGenHerwig::Stable(TParticle*  particle)
+{
+// Return true for a stable particle
+//
+    Int_t kf = TMath::Abs(particle->GetPdgCode());
+    
+    if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
+        
+    {
+       return kTRUE;
+    } else {
+       return kFALSE;
+    }
+}
+
+void AliGenHerwig::FinishRun()
+{
+  fHerwig->Hwefin();
+}
+
+AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
+{
+// Assignment operator
+    return *this;
+}
+
+
+extern "C" {
+  Double_t hwr_() {return sRandom->Rndm();}
+}

diff --git a/EVGEN/AliGenHerwig.h b/EVGEN/AliGenHerwig.h
new file mode 100644 (file)
index 0000000..8051f56
--- /dev/null
+++ b/EVGEN/AliGenHerwig.h
@@ -0,0 +1,107 @@
+#ifndef ALIGENHERWIG_H
+#define ALIGENHERWIG_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+// Generator using HERWIG as an external generator
+// The main HERWIG options are accessable for the user through this interface.
+
+
+#include "AliGenMC.h"
+#include <TString.h>
+#include <TArrayI.h>
+#include <AliRndm.h>
+#include <AliStructFuncType.h>
+
+class THerwig6;
+class TArrayI;
+class TParticle;
+class TClonesArray;
+
+
+class AliGenHerwig : public AliGenMC
+
+{
+    enum {kNoTrigger, kHardProcesses, kDirectPhotons};
+
+ public:
+    AliGenHerwig();
+    AliGenHerwig(Int_t npart);
+    AliGenHerwig(const AliGenHerwig &Herwig);
+    virtual ~AliGenHerwig();
+    virtual void    Generate();
+    virtual void    Init();
+    // set centre of mass energy
+    virtual void    SetBeamMomenta(Float_t p1=7000., Float_t p2=7000.)
+       {fMomentum1 = p1; fMomentum2 = p2;}
+    virtual void    SetProjectile(TString proj="P")  {fProjectile = proj;}    
+    virtual void    SetTarget(TString tar="P")       {fTarget = tar;}
+    virtual void    SetProcess(Int_t proc)       {fProcess = proc;}    
+    virtual void    KeepFullEvent();
+    virtual void    SetDecaysOff(Int_t flag=1)        {fDecaysOff = flag;}
+    virtual void    SetTrigger(Int_t flag=kNoTrigger) {fTrigger   = flag;}
+    virtual void    SetFlavor(Int_t flag=0)           {fFlavor    = flag;}    
+    virtual void    SetSelectAll(Int_t flag=0)        {fSelectAll = flag;}    
+    AliGenHerwig &  operator=(const AliGenHerwig & rhs);
+    virtual void    SetStrucFunc(StrucFunc_t func = kGRVHO) 
+      {fStrucFunc = func;}
+    virtual void    SetPtHardMin(Double_t pt) {fPtHardMin=pt;}
+    virtual void    SetPtRMS(Double_t pt) {fPtRMS=pt;}
+    virtual void    SetMaxPr(Int_t i) {fMaxPr=i;}
+    virtual void    SetMaxErrors(Int_t i) {fMaxErrors=i;}
+    virtual void    FinishRun();
+    virtual void    SetEnSoft(Double_t e) {fEnSoft=e;}
+ protected:
+    Bool_t SelectFlavor(Int_t pid);
+
+ protected:
+    TString     fProjectile;     // Projectile
+    TString     fTarget;         // Target
+    TString     fAutPDF;         // PDF group
+    Int_t       fModPDF;         // PDF set
+    StrucFunc_t fStrucFunc;      //Structure Function
+    Int_t       fKeep;           // Flag to keep full event information
+    Int_t       fDecaysOff;      // Flag to turn off decays of pi0, K_s, D, Lambda, sigma
+    Int_t       fTrigger;        // Trigger type
+    Int_t       fSelectAll;      // Flag to write the full event
+    Int_t       fFlavor;         // Selected particle flavor 4: charm+beauty 5: beauty
+    Float_t     fEnergyCMS;      // Centre of mass energy
+    Float_t     fMomentum1;      // Momentum of projectile
+    Float_t     fMomentum2;      // Momentum of target
+    Float_t     fKineBias;       // Bias from kinematic selection
+    Int_t       fTrials;         // Number of trials
+    TArrayI     fParentSelect;   // Parent particles to be selected 
+    TArrayI     fChildSelect;    // Decay products to be selected
+    Float_t     fXsection;       // Cross-section
+    THerwig6    *fHerwig;        // Herwig
+    Int_t       fProcess;        // Process number
+    Double_t    fPtHardMin;      // lower pT-hard cut 
+    Double_t    fPtRMS;          // intrinsic pt of incoming hadrons
+    Int_t       fMaxPr;          // maximum number of events to print out
+    Int_t       fMaxErrors;      // maximum number of errors allowed
+    Double_t    fEnSoft;          // change on soft energy distribution
+      
+ private:
+    // check if particle is selected as parent particle
+    Bool_t ParentSelected(Int_t ip);
+    // check if particle is selected as child particle
+    Bool_t ChildSelected(Int_t ip);
+    // adjust the weight from kinematic cuts
+    void   AdjustWeights();
+    // check seleted daughters
+    Bool_t DaughtersSelection(TParticle* iparticle, TClonesArray* particles);
+    // check if stable
+    Bool_t Stable(TParticle*  particle);
+    
+    void InitPDF();
+
+    ClassDef(AliGenHerwig,1) // AliGenerator interface to Herwig
+};
+#endif
+
+
+
+
+