]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - THerwig/AliGenHerwig.cxx
Removing the lib files
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx
index 9578c55ab6af27c9bb85644134ef74ea135ce62d..f33112a4dc7d17384c9ec1529ddeaffc13acecaf 100644 (file)
@@ -1,3 +1,4 @@
+
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
 #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)
 
 
@@ -54,13 +59,21 @@ ClassImp(AliGenHerwig)
     fHerwig(0x0),
     fProcess(0),
     fPtHardMin(0.),
+    fPtHardMax(9999.),
     fPtRMS(0.),
     fMaxPr(10),
     fMaxErrors(1000),
     fEnSoft(1),
-    fEv1Pr(0),
-    fEv2Pr(0),
-    fFileName(0)
+    fFileName(0),
+    fEtaMinParton(-20.),     
+    fEtaMaxParton(20.),     
+    fPhiMinParton(0.),     
+    fPhiMaxParton(2.* TMath::Pi()),     
+    fEtaMinGamma(-20.),      
+    fEtaMaxGamma(20.),      
+    fPhiMinGamma(0.),      
+    fPhiMaxGamma(2. * TMath::Pi()),
+    fHeader(0)
 {
 // Constructor
   fEnergyCMS = 14000;
@@ -84,14 +97,23 @@ AliGenHerwig::AliGenHerwig(Int_t npart)
     fHerwig(0x0),
     fProcess(0),
     fPtHardMin(0.),
+    fPtHardMax(9999.),
     fPtRMS(0.),
     fMaxPr(10),
     fMaxErrors(1000),
     fEnSoft(1),
-    fEv1Pr(0),
-    fEv2Pr(0),
-    fFileName(0)
+    fFileName(0),
+    fEtaMinParton(-20.),     
+    fEtaMaxParton(20.),     
+    fPhiMinParton(0.),     
+    fPhiMaxParton(2.* TMath::Pi()),     
+    fEtaMinGamma(-20.),      
+    fEtaMaxGamma(20.),      
+    fPhiMinGamma(0.),
+    fPhiMaxGamma(2. * TMath::Pi()),
+    fHeader(0)
 {
+// Constructor
     fEnergyCMS = 14000;
     SetTarget();
     SetProjectile();
@@ -104,13 +126,6 @@ AliGenHerwig::~AliGenHerwig()
 // Destructor
 }
 
-void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
-{
-  fEv1Pr = ++eventFirst;
-  fEv2Pr = ++eventLast;
-  if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
-}
-
 void AliGenHerwig::Init()
 {
 // Initialisation
@@ -119,18 +134,16 @@ 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);
+  printf("SetMAXPR %15d \n", fMaxPr);
   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
@@ -143,9 +156,9 @@ 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     ");
+  //fHerwig->Hwusta("PI0     ");
 
   // compute parameter dependent constants
   fHerwig->PrepareRun();
@@ -168,9 +181,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
@@ -183,9 +193,9 @@ 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     ");
+  //  fHerwig->Hwusta("PI0     ");
 
   // compute parameter dependent constants
   fHerwig->PrepareRunJimmy();
@@ -193,6 +203,7 @@ void AliGenHerwig::InitJimmy()
 
 void AliGenHerwig::InitPDF()
 {
+// Initialize PDF
   switch(fStrucFunc)
     {
 // ONLY USES LHAPDF STRUCTURE FUNCTIONS
@@ -255,10 +266,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
@@ -266,20 +276,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)
     {
@@ -289,8 +294,17 @@ void AliGenHerwig::Generate()
        Int_t np = particles->GetEntriesFast()-1;
        if (np == 0 ) continue;
 
-       Int_t nc=0;
+       //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 ;
+       } 
 
+       Int_t nc = 0;
+       fNprimaries = 0;
+       
        Int_t * newPos = new Int_t[np];
        for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
 
@@ -307,10 +321,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,
@@ -321,10 +337,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) {
@@ -333,11 +349,61 @@ void AliGenHerwig::Generate()
            }
        }
     }
+//
+  MakeHeader();
+//  
   SetHighWaterMark(nt);
 //  adjust weight due to kinematic selection
   AdjustWeights();
 //  get cross-section
   fXsection=fHerwig->GetAVWGT();
+  //printf(">> trials << %d\n",fTrials);
+}
+
+Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const TParticle* parton2)
+{
+// Check the kinematic trigger condition
+//
+//Select events with parton max energy
+    if(fPtHardMax < parton1->Pt()) return kFALSE;
+
+// Select events within angular window
+    Double_t eta[2];
+    eta[0] = parton1->Eta();
+    eta[1] = parton2->Eta();
+    Double_t phi[2];
+    phi[0] = parton1->Phi();
+    phi[1] = parton2->Phi();
+    Int_t    pdg[2]; 
+    pdg[0] = parton1->GetPdgCode();
+    pdg[1] = parton2->GetPdgCode();   
+    printf("min %f, max %f\n",fPtHardMin, fPtHardMax);
+    printf("Parton 1: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton1->GetName(),parton1->Pt(), eta[0],phi[0]*TMath::RadToDeg());
+    printf("Parton 2: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton2->GetName(),parton2->Pt(), eta[1],phi[1]*TMath::RadToDeg());
+    
+    if (fProcess == kHeJets) {
+      //Check if one of the 2 outgoing partons are in the eta-phi window
+      for(Int_t i = 0; i < 2; i++)
+       if ((eta[i] < fEtaMaxParton  && eta[i] > fEtaMinParton) &&
+           (phi[i] < fPhiMaxParton  && phi[i] > fPhiMinParton)) return  kTRUE ;
+    }
+    
+    else {
+      //Check if the gamma and the jet  are in the eta-phi window
+      Int_t igj = 0;
+      Int_t ijj = 0;
+      if(pdg[0] == 22) ijj=1;
+      else igj=1;
+      if ((eta[ijj] < fEtaMaxParton   && eta[ijj] > fEtaMinParton) &&
+         (phi[ijj] < fPhiMaxParton   && phi[ijj] > fPhiMinParton)) {
+       
+       if ((eta[igj] < fEtaMaxGamma   && eta[igj] > fEtaMinGamma) &&
+           (phi[igj] < fPhiMaxGamma   && phi[igj] > fPhiMinGamma)) return  kTRUE;
+       
+      }
+    }
+
+    return kFALSE ;
 }
 
 void AliGenHerwig::AdjustWeights()
@@ -357,7 +423,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
@@ -387,7 +453,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
@@ -400,7 +466,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
 //
@@ -427,3 +493,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;  
+}