Fixes to remove compilation warnings
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx
index ff8b535..e2ce8fc 100644 (file)
 // The main Herwig options are accessable for the user through this interface.
 // Uses the THerwig implementation of TGenerator.
 
-#include "AliGenHerwig.h"
-#include "AliHerwigRndm.h"
-#include "AliRun.h"
 
+#include <Riostream.h>
+#include <TClonesArray.h>
 #include <TParticle.h>
-#include "THerwig6.h"
 
-#include "Riostream.h"
-#include "AliMC.h"
+#include <THerwig6.h>
 
+#include "AliGenHerwig.h"
+#include "AliHerwigRndm.h"
+#include "AliMC.h"
+#include "AliRun.h"
 #include "driver.h"
 
 ClassImp(AliGenHerwig)
@@ -37,15 +38,14 @@ ClassImp(AliGenHerwig)
 
   AliGenHerwig::AliGenHerwig() :
     AliGenMC(),
-    fAutPDF("GRV"),
-    fModPDF(5),
-    fStrucFunc(kGRVHO),
+    fAutPDF("LHAPDF"),
+    fModPDF(19070),
+    fStrucFunc(kCTEQ5L),
     fKeep(0),
     fDecaysOff(1),
     fTrigger(0),
     fSelectAll(0),
     fFlavor(0),
-    fEnergyCMS(14000),
     fMomentum1(7000),
     fMomentum2(7000),
     fKineBias(1),
@@ -53,48 +53,70 @@ ClassImp(AliGenHerwig)
     fXsection(0),
     fHerwig(0x0),
     fProcess(0),
-    fPtHardMin(0),
-    fPtRMS(0),
+    fPtHardMin(0.),
+    fPtHardMax(9999.),
+    fPtRMS(0.),
     fMaxPr(10),
     fMaxErrors(1000),
     fEnSoft(1),
     fEv1Pr(0),
-    fEv2Pr(0)
+    fEv2Pr(0),
+    fFileName(0),
+    fEtaMinParton(-20.),     
+    fEtaMaxParton(20.),     
+    fPhiMinParton(0.),     
+    fPhiMaxParton(2.* TMath::Pi()),     
+    fEtaMinGamma(-20.),      
+    fEtaMaxGamma(20.),      
+    fPhiMinGamma(0.),      
+    fPhiMaxGamma(2. * TMath::Pi())  
 {
 // Constructor
+  fEnergyCMS = 14000;
 }
 
 AliGenHerwig::AliGenHerwig(Int_t npart)
-    :AliGenMC(npart)
+    :AliGenMC(npart),
+    fAutPDF("LHAPDF"),
+    fModPDF(19070),
+    fStrucFunc(kCTEQ5L),
+    fKeep(0),
+    fDecaysOff(1),
+    fTrigger(0),
+    fSelectAll(0),
+    fFlavor(0),
+    fMomentum1(7000),
+    fMomentum2(7000),
+    fKineBias(1),
+    fTrials(0),
+    fXsection(0),
+    fHerwig(0x0),
+    fProcess(0),
+    fPtHardMin(0.),
+    fPtHardMax(9999.),
+    fPtRMS(0.),
+    fMaxPr(10),
+    fMaxErrors(1000),
+    fEnSoft(1),
+    fEv1Pr(0),
+    fEv2Pr(0),
+    fFileName(0),
+    fEtaMinParton(-20.),     
+    fEtaMaxParton(20.),     
+    fPhiMinParton(0.),     
+    fPhiMaxParton(2.* TMath::Pi()),     
+    fEtaMinGamma(-20.),      
+    fEtaMaxGamma(20.),      
+    fPhiMinGamma(0.),      
+    fPhiMaxGamma(2. * TMath::Pi())  
 {
-    SetBeamMomenta();
+    fEnergyCMS = 14000;
     SetTarget();
     SetProjectile();
-    SetStrucFunc(kGRVLO98);
-    fKeep=0;
-    fTrigger=0;
-    fDecaysOff=1;
-    fSelectAll=0;
-    fFlavor=0;
-    fPtHardMin=0.;
-    fPtRMS=0.0;
-    fEnSoft=1.0;
-    fMaxPr=10;
-    fMaxErrors=1000;
-    fEv1Pr=0;
-    fEv2Pr=0;
-    // Set random number generator   
+    // Set random number generator
     AliHerwigRndm::SetHerwigRandom(GetRandom());
 }
 
-AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
-    :AliGenMC(Herwig)
-{
-// Copy constructor
-    Herwig.Copy(*this);
-}
-
-
 AliGenHerwig::~AliGenHerwig()
 {
 // Destructor
@@ -119,11 +141,12 @@ void AliGenHerwig::Init()
   // 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);
 
@@ -132,70 +155,115 @@ void AliGenHerwig::Init()
 //       RMASS(2)=0.32
 //       RMASS(3)=0.5
 //       RMASS(4)=1.55
-//       RMASS(5)=4.95
+//       RMASS(5)=4.75
 //       RMASS(6)=174.3
 //       RMASS(13)=0.75
 
   fHerwig->SetRMASS(4,1.2);
   fHerwig->SetRMASS(5,4.75);
-  
+
   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
 
-  fHerwig->Hwusta("PI0     ");
+  //fHerwig->Hwusta("PI0     ");
 
   // compute parameter dependent constants
   fHerwig->PrepareRun();
 }
 
+void AliGenHerwig::InitJimmy()
+{
+// Initialisation
+  fTarget.Resize(8);
+  fProjectile.Resize(8);
+  SetMC(new THerwig6());
+  fHerwig=(THerwig6*) fMCEvGen;
+  // initialize common blocks
+  fHerwig->InitializeJimmy(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);
+
+  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
+//       RMASS(3)=0.5
+//       RMASS(4)=1.55
+//       RMASS(5)=4.75
+//       RMASS(6)=174.3
+//       RMASS(13)=0.75
+
+  fHerwig->SetRMASS(4,1.2);
+  fHerwig->SetRMASS(5,4.75);
+
+  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
+
+  //  fHerwig->Hwusta("PI0     ");
+
+  // compute parameter dependent constants
+  fHerwig->PrepareRunJimmy();
+}
+
 void AliGenHerwig::InitPDF()
 {
   switch(fStrucFunc)
     {
-    case kGRVLO:
-      fModPDF=5;
-      fAutPDF="GRV";
-      break;
-    case kGRVHO:
-      fModPDF=6;
-      fAutPDF="GRV";
-      break;
+// ONLY USES LHAPDF STRUCTURE FUNCTIONS
     case kGRVLO98:
-      fModPDF=12;
-      fAutPDF="GRV";
+      fModPDF=80060;
+      fAutPDF="HWLHAPDF";
       break;
-    case kMRSDminus:
-      fModPDF=31;
-      fAutPDF="MRS";
+    case kCTEQ6:
+      fModPDF=10040;
+      fAutPDF="HWLHAPDF";
       break;
-    case kMRSD0:
-      fModPDF=30;
-      fAutPDF="MRS";
+    case kCTEQ61:
+      fModPDF=10100;
+      fAutPDF="HWLHAPDF";
       break;
-    case kMRSG:
-      fModPDF=41;
-      fAutPDF="MRS";
+    case kCTEQ6m:
+      fModPDF=10050;
+      fAutPDF="HWLHAPDF";
       break;
-    case kMRSTcgLO:
-      fModPDF=72;
-      fAutPDF="MRS";
+    case kCTEQ6l:
+      fModPDF=10041;
+      fAutPDF="HWLHAPDF";
       break;
-    case kCTEQ4M:
-      fModPDF=34;
-      fAutPDF="CTEQ";
+    case kCTEQ6ll:
+      fModPDF=10042;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ5M:
+      fModPDF=19050;
+      fAutPDF="HWLHAPDF";
       break;
     case kCTEQ5L:
-      fModPDF=46;
-      fAutPDF="CTEQ";
+      fModPDF=19070;
+      fAutPDF="HWLHAPDF";
       break;
-    case kCTEQ5M:
-      fModPDF=48;
-      fAutPDF="CTEQ";
+    case kCTEQ4M:
+      fModPDF=19150;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ4L:
+      fModPDF=19170;
+      fAutPDF="HWLHAPDF";
       break;
+//    case kMRST2004nlo:
+//      fModPDF=20400;
+//      fAutPDF="HWLHAPDF";
+//      break;
     default:
       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
       break;
     }
-  fAutPDF.Resize(20);      
+  fAutPDF.Resize(20);
   fHerwig->SetMODPDF(1,fModPDF);
   fHerwig->SetMODPDF(2,fModPDF);
   fHerwig->SetAUTPDF(1,fAutPDF);
@@ -210,7 +278,7 @@ void AliGenHerwig::Generate()
   Float_t origin[3]=   {0,0,0};
   Float_t origin0[3]=  {0,0,0};
   Float_t p[4], random[6];
-  
+
   static TClonesArray *particles;
   //  converts from mm/c to s
   const Float_t kconv=0.001/2.999792458e8;
@@ -219,9 +287,9 @@ void AliGenHerwig::Generate()
   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) {
@@ -231,7 +299,7 @@ void AliGenHerwig::Generate()
        TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
     }
   }
-  
+
   while(1)
     {
        fHerwig->GenerateEvent();
@@ -239,19 +307,27 @@ void AliGenHerwig::Generate()
        fHerwig->ImportParticles(particles,"All");
        Int_t np = particles->GetEntriesFast()-1;
        if (np == 0 ) continue;
-       
+
+       //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;
-       
+
        Int_t * newPos = new Int_t[np];
        for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
-       
+
        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 != 3 && 
-               KinematicSelection(iparticle,0)) 
+           if (ks != 3 &&
+               KinematicSelection(iparticle,0))
            {
                nc++;
                p[0]=iparticle->Px();
@@ -266,7 +342,7 @@ void AliGenHerwig::Generate()
                Int_t   trackIt = (ks == 1) && fTrackIt;
                PushTrack(trackIt, iparent, kf,
                          p[0], p[1], p[2], p[3],
-                         origin[0], origin[1], origin[2], 
+                         origin[0], origin[1], origin[2],
                          tof,
                          polar[0], polar[1], polar[2],
                          kPPrimary, nt, fHerwig->GetEVWGT(), ks);
@@ -289,6 +365,53 @@ void AliGenHerwig::Generate()
   AdjustWeights();
 //  get cross-section
   fXsection=fHerwig->GetAVWGT();
+  //printf(">> trials << %d\n",fTrials);
+}
+
+Bool_t AliGenHerwig::CheckParton(TParticle* parton1, 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()
@@ -301,7 +424,7 @@ void AliGenHerwig::AdjustWeights()
         part->SetWeight(part->GetWeight()*fKineBias);
     }
 }
+
 
 void AliGenHerwig::KeepFullEvent()
 {
@@ -321,9 +444,9 @@ Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* part
     Bool_t selected=kFALSE;
     if (hasDaughters) {
        imin=iparticle->GetFirstDaughter();
-       imax=iparticle->GetLastDaughter();       
+       imax=iparticle->GetLastDaughter();
        for (i=imin; i<= imax; i++){
-           TParticle *  jparticle       = (TParticle *) particles->At(i);      
+           TParticle *  jparticle       = (TParticle *) particles->At(i);
            Int_t ip=jparticle->GetPdgCode();
            if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
                selected=kTRUE; break;
@@ -345,7 +468,7 @@ Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
 // 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);
@@ -356,9 +479,9 @@ 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 {
@@ -371,12 +494,10 @@ void AliGenHerwig::FinishRun()
   fHerwig->Hwefin();
 }
 
-
-AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
+void AliGenHerwig::FinishRunJimmy()
 {
-// Assignment operator
-    rhs.Copy(*this);
-    return (*this);
-}
+  fHerwig->Hwefin();
+  fHerwig->Jmefin();
 
+}