]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - THerwig/AliGenHerwig.cxx
Adding extra check for GPU_FORCE_64BIT_PTR env var
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx
index a6779ce75be0150eb1e26bfb15500ef0a0c20194..17c5c7a1abfe93beb58e2b6faaf7593d79a78b08 100644 (file)
 #include <THerwig6.h>
 
 #include "AliGenHerwig.h"
+#include "AliGenHerwigEventHeader.h"
 #include "AliHerwigRndm.h"
 #include "AliMC.h"
 #include "AliRun.h"
 #include "driver.h"
 
+using std::cerr;
+using std::endl;
+
 ClassImp(AliGenHerwig)
 
 
@@ -69,7 +73,8 @@ ClassImp(AliGenHerwig)
     fEtaMinGamma(-20.),      
     fEtaMaxGamma(20.),      
     fPhiMinGamma(0.),      
-    fPhiMaxGamma(2. * TMath::Pi())  
+    fPhiMaxGamma(2. * TMath::Pi()),
+    fHeader(0)
 {
 // Constructor
   fEnergyCMS = 14000;
@@ -107,9 +112,11 @@ AliGenHerwig::AliGenHerwig(Int_t npart)
     fPhiMaxParton(2.* TMath::Pi()),     
     fEtaMinGamma(-20.),      
     fEtaMaxGamma(20.),      
-    fPhiMinGamma(0.),      
-    fPhiMaxGamma(2. * TMath::Pi())  
+    fPhiMinGamma(0.),
+    fPhiMaxGamma(2. * TMath::Pi()),
+    fHeader(0)
 {
+// Constructor
     fEnergyCMS = 14000;
     SetTarget();
     SetProjectile();
@@ -124,9 +131,9 @@ AliGenHerwig::~AliGenHerwig()
 
 void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 {
-  fEv1Pr = ++eventFirst;
-  fEv2Pr = ++eventLast;
-  if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
+  fEv1Pr = eventFirst;
+  fEv2Pr = eventLast;
+  if ( fEv2Pr == -1 ) fEv2Pr = fEv1Pr;
 }
 
 void AliGenHerwig::Init()
@@ -137,18 +144,15 @@ void AliGenHerwig::Init()
   SetMC(new THerwig6());
   fHerwig=(THerwig6*) fMCEvGen;
   // initialize common blocks
-  fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
+  fHerwig->Initialize(fProjectile.Data(), fTarget.Data(), fMomentum1, fMomentum2, fProcess); // 
   // reset parameters according to user needs
   InitPDF();
   fHerwig->SetPTMIN(fPtHardMin);
+  fHerwig->SetPTMAX(fPtHardMax);
   fHerwig->SetPTRMS(fPtRMS);
   fHerwig->SetMAXPR(fMaxPr);
   fHerwig->SetMAXER(fMaxErrors);
   fHerwig->SetENSOF(fEnSoft);
-
-  fHerwig->SetEV1PR(fEv1Pr);
-  fHerwig->SetEV2PR(fEv2Pr);
-
 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
 //       RMASS(1)=0.32
 //       RMASS(2)=0.32
@@ -161,7 +165,7 @@ void AliGenHerwig::Init()
   fHerwig->SetRMASS(4,1.2);
   fHerwig->SetRMASS(5,4.75);
 
-  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
+  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
 
   //fHerwig->Hwusta("PI0     ");
 
@@ -186,9 +190,6 @@ void AliGenHerwig::InitJimmy()
   fHerwig->SetMAXER(fMaxErrors);
   fHerwig->SetENSOF(fEnSoft);
 
-  fHerwig->SetEV1PR(fEv1Pr);
-  fHerwig->SetEV2PR(fEv2Pr);
-
 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
 //       RMASS(1)=0.32
 //       RMASS(2)=0.32
@@ -201,7 +202,7 @@ void AliGenHerwig::InitJimmy()
   fHerwig->SetRMASS(4,1.2);
   fHerwig->SetRMASS(5,4.75);
 
-  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
+  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
 
   //  fHerwig->Hwusta("PI0     ");
 
@@ -211,6 +212,7 @@ void AliGenHerwig::InitJimmy()
 
 void AliGenHerwig::InitPDF()
 {
+// Initialize PDF
   switch(fStrucFunc)
     {
 // ONLY USES LHAPDF STRUCTURE FUNCTIONS
@@ -273,10 +275,9 @@ 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[4], random[6];
+  Float_t polar[3]  =   {0,0,0};
+  Float_t origin[3] =   {0,0,0};
+  Float_t p[4];
 
   static TClonesArray *particles;
   //  converts from mm/c to s
@@ -284,20 +285,15 @@ void AliGenHerwig::Generate()
   //
   Int_t nt=0;
   Int_t jev=0;
-  Int_t j, kf, ks, imo;
+  Int_t 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]));
-    }
-  }
+
+  //  Set collision vertex position 
+  if (fVertexSmear == kPerEvent) Vertex();
 
   while(1)
     {
@@ -310,13 +306,22 @@ void AliGenHerwig::Generate()
        //Check hard partons or direct gamma in kine range
 
        if (fProcess == kHeJets || fProcess == kHeDirectGamma) {
-         TParticle* parton1 = (TParticle *) particles->At(6);
-         TParticle* parton2 = (TParticle *) particles->At(7);
-         if (!CheckParton(parton1, parton2))  continue ;
+           TParticle* parton1 = (TParticle *) particles->At(6);
+           TParticle* parton2 = (TParticle *) particles->At(7);
+           if (!CheckParton(parton1, parton2))  continue ;
        } 
 
-       Int_t nc=0;
+       // 
+       if (gAlice) {
+           if (gAlice->GetEvNumber()>=fEv1Pr &&
+               gAlice->GetEvNumber()<=fEv2Pr) fHerwig->PrintEvt();
+
+       }
 
+       
+       Int_t nc = 0;
+       fNprimaries = 0;
+       
        Int_t * newPos = new Int_t[np];
        for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
 
@@ -333,10 +338,12 @@ void AliGenHerwig::Generate()
                p[1]=iparticle->Py();
                p[2]=iparticle->Pz();
                p[3]=iparticle->Energy();
-               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();
+
+               origin[0] = fVertex[0] + iparticle->Vx()/10; // [cm]
+               origin[1] = fVertex[1] + iparticle->Vy()/10; // [cm]
+               origin[2] = fVertex[2] + iparticle->Vz()/10; // [cm]
+
+               Float_t tof = fTime + kconv*iparticle->T();
                Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
                Int_t   trackIt = (ks == 1) && fTrackIt;
                PushTrack(trackIt, iparent, kf,
@@ -347,10 +354,10 @@ void AliGenHerwig::Generate()
                          kPPrimary, nt, fHerwig->GetEVWGT(), ks);
                KeepTrack(nt);
                newPos[i]=nt;
+               fNprimaries++;
            } // end of if: selection of particle
        } // end of for: particle loop
        if (newPos) delete[] newPos;
-       //      MakeHeader();
        if (nc > 0) {
            jev+=nc;
            if (jev >= fNpart || fNpart == -1) {
@@ -359,6 +366,9 @@ void AliGenHerwig::Generate()
            }
        }
     }
+//
+  MakeHeader();
+//  
   SetHighWaterMark(nt);
 //  adjust weight due to kinematic selection
   AdjustWeights();
@@ -367,7 +377,7 @@ void AliGenHerwig::Generate()
   //printf(">> trials << %d\n",fTrials);
 }
 
-Bool_t AliGenHerwig::CheckParton(TParticle* parton1, TParticle* parton2)
+Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const TParticle* parton2)
 {
 // Check the kinematic trigger condition
 //
@@ -430,7 +440,7 @@ void AliGenHerwig::KeepFullEvent()
     fKeep=1;
 }
 
-Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
+Bool_t AliGenHerwig::DaughtersSelection(const TParticle* iparticle, const TClonesArray* particles)
 {
 //
 // Looks recursively if one of the daughters has been selected
@@ -460,7 +470,7 @@ Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* part
 }
 
 
-Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
+Bool_t AliGenHerwig::SelectFlavor(Int_t pid) const
 {
 // Select flavor of particle
 // 0: all
@@ -473,7 +483,7 @@ Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
     return (fFlavor == ifl);
 }
 
-Bool_t AliGenHerwig::Stable(TParticle*  particle)
+Bool_t AliGenHerwig::Stable(const TParticle*  particle) const
 {
 // Return true for a stable particle
 //
@@ -500,3 +510,32 @@ void AliGenHerwig::FinishRunJimmy()
 
 }
 
+
+void AliGenHerwig::MakeHeader()
+{
+//
+// Make header for the simulated event
+// 
+  if (fHeader) delete fHeader;
+  fHeader = new AliGenHerwigEventHeader("Herwig");
+//
+// Event type  
+    ((AliGenHerwigEventHeader*) fHeader)->SetProcessType(fHerwig->GetIHPRO());
+//
+// Number of trials
+    ((AliGenHerwigEventHeader*) fHeader)->SetTrials(fTrials);
+//
+// Event weight (cross section)
+    ((AliGenHerwigEventHeader*) fHeader)->SetWeight(fHerwig->GetEVWGT());
+//
+// Event Vertex 
+    fHeader->SetPrimaryVertex(fVertex);
+    fHeader->SetInteractionTime(fTime);
+//
+// Number of primaries
+    fHeader->SetNProduced(fNprimaries);
+//  Pass header
+//
+    AddHeader(fHeader);
+    fHeader = 0x0;  
+}