AliGenerator interface class to HIJING using THijing (test version)

diff --git a/EVGEN/AliGenHijing.cxx b/EVGEN/AliGenHijing.cxx
new file mode 100644 (file)
index 0000000..073da1c
--- /dev/null
+++ b/EVGEN/AliGenHijing.cxx
@@ -0,0 +1,333 @@
+/**************************************************************************
+ * 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$
+*/
+
+#include "AliGenHijing.h"
+#include "AliRun.h"
+
+#include <TArrayI.h>
+#include <TParticle.h>
+#include <THijing.h>
+
+ ClassImp(AliGenHijing)
+
+AliGenHijing::AliGenHijing()
+                 :AliGenerator()
+{
+// Constructor
+}
+
+AliGenHijing::AliGenHijing(Int_t npart)
+                 :AliGenerator(npart)
+{
+// Default PbPb collisions at 5. 5 TeV
+//
+    SetEnergyCMS();
+    SetImpactParameterRange();
+    SetTarget();
+    SetProjectile();
+    fKeep=0;
+    fQuench=1;
+    fShadowing=1;
+    fTrigger=0;
+    fDecaysOff=1;
+    fEvaluate=0;
+    fSelectAll=0;
+}
+
+AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
+{
+// copy constructor
+}
+
+
+AliGenHijing::~AliGenHijing()
+{
+// Destructor
+}
+
+void AliGenHijing::Init()
+{
+// Initialisation
+    SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
+                     fAProjectile, fZProjectile, fATarget, fZTarget, 
+                     fMinImpactParam, fMaxImpactParam));
+
+    fHijing=(THijing*) fgMCEvGen;
+    fHijing->Initialize();
+    fHijing->SetIHPR2(3,  fTrigger);
+    fHijing->SetIHPR2(4,  fQuench);
+    fHijing->SetIHPR2(6,  fShadowing);
+    fHijing->SetIHPR2(12, fDecaysOff);    
+    fHijing->SetIHPR2(21, fKeep);
+//
+    if (fEvaluate) EvaluateCrossSections();
+}
+
+void AliGenHijing::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];
+    Float_t tof;
+
+    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;
+
+    if(!particles) particles=new TClonesArray("TParticle",10000);
+    
+    fTrials=0;
+    for (j=0;j<3;j++) origin0[j]=fOrigin[j];
+    if(fVertexSmear==perEvent) {
+       gMC->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]));
+//         fHijing->SetMSTP(151,0);
+       }
+    } else if (fVertexSmear==perTrack) {
+//     fHijing->SetMSTP(151,0);
+       for (j=0;j<3;j++) {
+//         fHijing->SetPARP(151+j, fOsigma[j]*10.);
+       }
+    }
+    while(1)
+    {
+
+       fHijing->GenerateEvent();
+       fTrials++;
+       fHijing->ImportParticles(particles,"All");
+       Int_t np = particles->GetEntriesFast();
+       printf("\n **************************************************%d\n",np);
+       Int_t nc=0;
+       if (np == 0 ) continue;
+       Int_t i;
+       Int_t * newPos = new Int_t[np];
+       for (i = 0; i<np-1; i++) *(newPos+i)=i;
+       
+       for (i = 0; i<np-1; i++) {
+           TParticle *  iparticle       = (TParticle *) particles->At(i);
+
+           Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
+           Bool_t  hasDaughter          =  (iparticle->GetFirstDaughter() >=0);
+           Bool_t  selected             =  kTRUE;
+           Bool_t  hasSelectedDaughters =  kFALSE;
+
+           if (!fSelectAll) selected = KinematicSelection(iparticle);
+           kf        = iparticle->GetPdgCode();
+//         Int_t id1=iparticle->GetFirstDaughter();
+//         Int_t id2=iparticle->GetLastDaughter();         
+//         printf("\n particle  %d %d %d %d %d %d\n",i, kf, id1, id2, hasDaughter, selected);
+
+           if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
+//
+// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
+//
+
+           if (selected || hasSelectedDaughters) {
+               nc++;
+               ks        = iparticle->GetStatusCode();
+               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;
+               tof=kconv*iparticle->T();
+               imo=-1;
+               if (hasMother) {
+                   imo=iparticle->GetFirstMother();
+                   imo=*(newPos+imo);
+               }
+               
+//             printf("\n selected iparent %d %d %d \n",i, kf, imo);
+               if (hasDaughter) {
+                   gAlice->SetTrack(0,imo,kf,p,origin,polar,
+                                    tof,"Primary",nt);
+               } else {
+                   gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
+                                    tof,"Secondary",nt);
+               }
+               *(newPos+i)=nt;
+           } // selected
+       } // particle loop 
+       delete newPos;
+       printf("\n I've put %i particles on the stack \n",nc);
+       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;
+           }
+       }
+    } // event loop
+}
+
+Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
+{
+// Perform kinematic selection
+    Float_t px=particle->Px();
+    Float_t py=particle->Py();
+    Float_t pz=particle->Pz();
+    Float_t  e=particle->Energy();
+
+//
+//  transverse momentum cut    
+    Float_t pt=TMath::Sqrt(px*px+py*py);
+    if (pt > fPtMax || pt < fPtMin) 
+    {
+//     printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
+       return kFALSE;
+    }
+//
+// momentum cut
+    Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
+    if (p > fPMax || p < fPMin) 
+    {
+//     printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
+       return kFALSE;
+    }
+    
+//
+// theta cut
+    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
+    if (theta > fThetaMax || theta < fThetaMin) 
+    {
+       
+//     printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
+       return kFALSE;
+    }
+
+//
+// rapidity cut
+    Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
+    if (y > fYMax || y < fYMin)
+    {
+//     printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
+       return kFALSE;
+    }
+
+//
+// phi cut
+    Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
+    if (phi > fPhiMax || phi < fPhiMin)
+    {
+//     printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
+       return kFALSE;
+    }
+
+    return kTRUE;
+}
+
+void AliGenHijing::KeepFullEvent()
+{
+    fKeep=1;
+}
+
+void AliGenHijing::EvaluateCrossSections()
+{
+//     Glauber Calculation of geometrical x-section
+//
+    Float_t xTot=0.;          // barn
+    Float_t xTotHard=0.;      // barn 
+    Float_t xPart=0.;         // barn
+    Float_t xPartHard=0.;     // barn 
+    Float_t sigmaHard=0.1;    // mbarn
+    Float_t bMin=0.;
+    Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
+    const Float_t kdib=0.2;
+    Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;
+
+
+    printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
+    printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
+    Int_t i;
+    Float_t oldvalue=0.;
+    
+    for (i=0; i<kMax; i++)
+    {
+       Float_t xb=bMin+i*kdib;
+       Float_t ov;
+       ov=fHijing->Profile(xb);
+       Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
+       Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
+       xTot+=gb;
+       xTotHard+=gbh;
+       if (xb > fMinImpactParam && xb < fMaxImpactParam)
+       {
+           xPart+=gb;
+           xPartHard+=gbh;
+       }
+       
+       if ((xTot-oldvalue)/oldvalue<0.0001) break;
+       oldvalue=xTot;
+       printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
+       printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
+    }
+    printf("\n Total cross section (barn): %f \n",xTot);
+    printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
+    printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
+    printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
+}
+
+Bool_t AliGenHijing::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();       
+//     printf("\n Has daughters %d %d:", imin, imax);
+       for (i=imin; i<= imax; i++){
+           TParticle *  jparticle       = (TParticle *) particles->At(i);      
+//         Int_t ip=jparticle->GetPdgCode();
+//         printf("\n consider daughter %d %d", i,ip);
+           
+           if (KinematicSelection(jparticle)) {selected=kTRUE; break;}
+           if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
+       }
+    } else {
+       return kFALSE;
+    }
+
+    return selected;
+}
+
+
+
+AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
+{
+// Assignment operator
+    return *this;
+}

diff --git a/EVGEN/AliGenHijing.h b/EVGEN/AliGenHijing.h
new file mode 100644 (file)
index 0000000..23e4c89
--- /dev/null
+++ b/EVGEN/AliGenHijing.h
@@ -0,0 +1,104 @@
+#ifndef ALIGENHIJING_H
+#define ALIGENHIJING_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+
+#include "AliGenerator.h"
+#include "GenTypeDefs.h"
+#include <THijing.h>
+#include <TArrayI.h>
+
+class TArrayI;
+class TParticle;
+
+class AliGenHijing : public AliGenerator
+{
+    enum {kNoTrigger, kHardProcesses, kDirectPhotons};
+
+ public:
+    AliGenHijing();
+    AliGenHijing(Int_t npart);
+    AliGenHijing(const AliGenHijing &Hijing);
+    virtual ~AliGenHijing();
+    virtual void    Generate();
+    virtual void    Init();
+    // set centre of mass energy
+    virtual void    SetEnergyCMS(Float_t energy=5500) {fEnergyCMS=energy;}
+    virtual void    SetReferenceFrame(char* frame="CMS     ") {fFrame=frame;}
+    virtual void    SetProjectile(char* proj="A       ", Int_t a=208, Int_t z=97)
+       {fProjectile = proj;
+       fAProjectile = a;
+       fZProjectile = z;}    
+    virtual void    SetTarget(char* tar="A       ", Int_t a=208, Int_t z=97)
+       {fTarget = tar;
+       fATarget = a;
+       fZTarget =z;}    
+    virtual void    SetImpactParameterRange(Float_t bmin = 0, Float_t bmax =15.)
+       { fMinImpactParam=bmin;
+       fMaxImpactParam=bmax;
+       }
+    virtual void    KeepFullEvent();
+    virtual void    SetJetQuenching(Int_t flag=1)     {fQuench    = flag;}
+    virtual void    SetShadowing(Int_t flag=1)        {fShadowing = flag;}
+    virtual void    SetDecaysOff(Int_t flag=1)        {fDecaysOff = flag;}
+    virtual void    SetTrigger(Int_t flag=kNoTrigger) {fTrigger   = flag;}
+    virtual void    SetEvaluate(Int_t flag=0)         {fEvaluate  = flag;}
+    virtual void    SetSelectAll(Int_t flag=0)        {fSelectAll = flag;}    
+    AliGenHijing & operator=(const AliGenHijing & rhs);
+// Physics Routines        
+    virtual void EvaluateCrossSections();
+    
+    
+ protected:
+    char       *fFrame;         // Reference frame 
+    char       *fProjectile;    // Projectile
+    char       *fTarget;        // Target
+    Int_t       fAProjectile;   // Projectile A
+    Int_t       fZProjectile;   // Projectile Z
+    Int_t       fATarget;       // Target A
+    Int_t       fZTarget;       // Target Z
+    Float_t     fMinImpactParam; // minimum impact parameter
+    Float_t     fMaxImpactParam; // maximum impact parameter   
+    Int_t       fKeep;          // Flag to keep full event information
+    Int_t       fQuench;        // Flag to switch on jet quenching
+    Int_t       fShadowing;     // Flag to switch on voclear effects on parton distribution function
+    Int_t       fDecaysOff;     // Flag to turn off decays of pi0, K_s, D, Lambda, sigma
+    Int_t       fTrigger;       // Trigger type
+    Int_t       fEvaluate;      // Evaluate total and partial cross-sections
+    Int_t       fSelectAll;     // Flag to write the full event
+    
+    Process_t   fProcess;       // Process type
+    StrucFunc_t fStrucFunc;     // Structure Function
+    Decay_t     fForceDecay;    // Decay channel  are forced
+    Float_t     fEnergyCMS;     // Centre of mass energy
+    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
+    THijing    *fHijing;        // Hijing
+    Float_t     fPtHardMin;     // lower pT-hard cut 
+    Float_t     fPtHardMax;     // higher pT-hard cut
+
+ 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);
+    // all kinematic selection cuts go here 
+    Bool_t KinematicSelection(TParticle *particle);
+    // adjust the weight from kinematic cuts
+    void   AdjustWeights();
+    // check seleted daughters
+    Bool_t DaughtersSelection(TParticle* iparticle, TClonesArray* particles);
+    ClassDef(AliGenHijing,1) // AliGenerator interface to Hijing
+};
+#endif
+
+
+
+
+