fPatchOmegaDalitz(0),
fNpartons(0),
fReadFromFile(0),
+ fReadLHEF(0),
+ fPartonShower(1),
fQuench(0),
fQhat(0.),
fLength(0.),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
fPatchOmegaDalitz(0),
fNpartons(0),
fReadFromFile(kFALSE),
+ fReadLHEF(0),
+ fPartonShower(1),
fQuench(kFALSE),
fQhat(0.),
fLength(0.),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
fPythia->SetMSTP(91,0);
}
-
+ if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
+
if (fReadFromFile) {
fRL = AliRunLoader::Open(fkFileName, "Partons");
fRL->LoadKinematics();
case kPyCharmpPbMNR:
case kPyCharmppMNR:
case kPyCharmppMNRwmi:
+ case kPyCharmPWHG:
fParentSelect[0] = 411;
fParentSelect[1] = 421;
fParentSelect[2] = 431;
case kPyBeautypPbMNR:
case kPyBeautyppMNR:
case kPyBeautyppMNRwmi:
+ case kPyBeautyPWHG:
fParentSelect[0]= 511;
fParentSelect[1]= 521;
fParentSelect[2]= 531;
case kPyMbNonDiffr:
case kPyMbMSEL1:
case kPyJets:
+ case kPyJetsPWHG:
case kPyDirectGamma:
case kPyLhwgMb:
break;
fProcess != kPyZ &&
fProcess != kPyCharmppMNRwmi &&
fProcess != kPyBeautyppMNRwmi &&
- fProcess != kPyBeautyJets) {
+ fProcess != kPyBeautyJets &&
+ fProcess != kPyJetsPWHG &&
+ fProcess != kPyCharmPWHG &&
+ fProcess != kPyBeautyPWHG) {
for (i = 0; i < np; i++) {
TParticle* iparticle = (TParticle *) fParticles.At(i);
Int_t* pParent = new Int_t[np];
for (i=0; i< np; i++) pParent[i] = -1;
- if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
+ if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
&& fEtMinJet > 0.) {
TParticle* jet1 = (TParticle *) fParticles.At(6);
TParticle* jet2 = (TParticle *) fParticles.At(7);
Int_t km = iparticle->GetFirstMother();
if (
(((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
- ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
+ ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
)
{
nc++;
// Jets that have triggered
//Need to store jets for b-jet studies too!
- if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
{
Int_t ntrig, njet;
Float_t jets[4][10];
//
// Store pt^hard
((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
+
+//
+// Store Event Weight
+ ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7));
+
//
// Pass header
//
pdg[1] = jet2->GetPdgCode();
Bool_t triggered = kFALSE;
- if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
+ if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
Int_t njets = 0;
Int_t ntrig = 0;
Float_t jets[4][10];
virtual void SetHadronisation(Int_t flag = 1) {fHadronisation = flag;}
virtual void SetPatchOmegaDalitz(Int_t flag = 1) {fPatchOmegaDalitz = flag;}
virtual void SetReadFromFile(const Text_t *filname) {fkFileName = filname; fReadFromFile = 1;}
+ virtual void SetReadLHEF(const Text_t *filename) {fkNameLHEF = filename; fReadLHEF = 1;}
//
// Pile-up
Bool_t fPatchOmegaDalitz; //flag for omega dalitz decay patch
Int_t fNpartons; //Number of partons before hadronisation
Int_t fReadFromFile; //read partons from file
+ Int_t fReadLHEF; //read lhef file
Int_t fQuench; //Flag for quenching
Float_t fQhat; //Transport coefficient (GeV^2/fm)
Float_t fLength; //Medium length (fm)
AliGenPythiaEventHeader* fHeader; //! Event header
AliRunLoader* fRL; //! Run Loader
const Text_t* fkFileName; //! Name of file to read from
-
+ const Text_t* fkNameLHEF; //! Name of lhef file to read from
+ Int_t fPartonShower;
Bool_t fFragPhotonInCalo; // Option to ask for Fragmentation Photon in calorimeters acceptance
Bool_t fHadronInCalo; // Option to ask for hadron (not pi0) in calorimeters acceptance
case kPyCharmpPbMNR:
case kPyCharmppMNR:
case kPyCharmppMNRwmi:
+ case kPyCharmPWHG:
fParentSelect[0] = 411;
fParentSelect[1] = 421;
fParentSelect[2] = 431;
case kPyBeautypPbMNR:
case kPyBeautyppMNR:
case kPyBeautyppMNRwmi:
+ case kPyBeautyPWHG:
fParentSelect[0]= 511;
fParentSelect[1]= 521;
fParentSelect[2]= 531;
case kPyMbNonDiffr:
case kPyMbMSEL1:
case kPyJets:
+ case kPyJetsPWHG:
case kPyDirectGamma:
case kPyLhwgMb:
break;
fProcess != kPyW &&
fProcess != kPyZ &&
fProcess != kPyCharmppMNRwmi &&
- fProcess != kPyBeautyppMNRwmi) {
+ fProcess != kPyBeautyppMNRwmi &&
+ fProcess != kPyJetsPWHG &&
+ fProcess != kPyCharmPWHG &&
+ fProcess != kPyBeautyPWHG) {
for (i = 0; i < np; i++) {
TParticle* iparticle = (TParticle *) fParticles.At(i);
Int_t* pParent = new Int_t[np];
for (i=0; i< np; i++) pParent[i] = -1;
- if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG ) {
TParticle* jet1 = (TParticle *) fParticles.At(6);
TParticle* jet2 = (TParticle *) fParticles.At(7);
if (!CheckTrigger(jet1, jet2)) {
}
// Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
- if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+ if ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && (fFragPhotonInCalo || fPi0InCalo) ) {
Bool_t ok = kFALSE;
// Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
- if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+ if ((fProcess == kPyJets || fProcess == kPyJetsPWHG || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
Bool_t okd = kFALSE;
Int_t km = iparticle->GetFirstMother();
if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
(ks != 1) ||
- (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+ ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && ks == 21 && km == 0 && i>1)) {
nc++;
if (ks == 1) trackIt = 1;
//
// Jets that have triggered
- if (fProcess == kPyJets)
+ if (fProcess == kPyJets || fProcess == kPyJetsPWHG)
{
Int_t ntrig, njet;
Float_t jets[4][10];
pdg[1] = jet2->GetPdgCode();
Bool_t triggered = kFALSE;
- if (fProcess == kPyJets) {
+ if (fProcess == kPyJets || fProcess == kPyJetsPWHG) {
Int_t njets = 0;
Int_t ntrig = 0;
Float_t jets[4][10];
case kPyMBRDoubleDiffraction:
case kPyMBRCentralDiffraction:
break;
+ case kPyJetsPWHG:
+ // N.B.
+ // ====
+ // For the case of jet production the following parameter setting
+ // limits the transverse momentum of secondary scatterings, due
+ // to multiple parton interactions, to be less than that of the
+ // primary interaction (see POWHEG Dijet paper arXiv:1012.3380
+ // [hep-ph] sec. 4.1 and also the PYTHIA Manual).
+ SetMSTP(86,1);
+
+ // maximum number of errors before pythia aborts (def=10)
+ SetMSTU(22,10);
+ // number of warnings printed on the shell
+ SetMSTU(26,20);
+ break;
+
+ case kPyCharmPWHG:
+ case kPyBeautyPWHG:
+ // number of warnings printed on the shell
+ SetMSTU(26,20);
+
+ break;
}
//
// Initialize PYTHIA
//
SetMSTP(41,1); // all resonance decays switched on
- Initialize("CMS",fProjectile,fTarget,fEcms);
+ if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG) {
+ Initialize("USER","","",0.);
+ } else {
+ Initialize("CMS",fProjectile,fTarget,fEcms);
+ }
fOmegaDalitz.Init();
}
case kPyMBRDoubleDiffraction:
case kPyMBRCentralDiffraction:
break;
-
+ case kPyJetsPWHG:
+ // N.B.
+ // ====
+ // For the case of jet production the following parameter setting
+ // limits the transverse momentum of secondary scatterings, due
+ // to multiple parton interactions, to be less than that of the
+ // primary interaction (see POWHEG Dijet paper arXiv:1012.3380
+ // [hep-ph] sec. 4.1 and also the PYTHIA Manual).
+ SetMSTP(86,1);
+ // maximum number of errors before pythia aborts (def=10)
+ SetMSTU(22,10);
+ // number of warnings printed on the shell
+ SetMSTU(26,20);
+ break;
}
//
// Initialize PYTHIA
SetMSTP(41,1); // all resonance decays switched on
- Initialize("CMS",fProjectile,fTarget,fEcms);
-
+ if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG) {
+ Initialize("USER","","",0.);
+ } else {
+ Initialize("CMS",fProjectile,fTarget,fEcms);
+ }
}
Int_t AliPythia6::CheckedLuComp(Int_t kf)
set ( CSRCS pythia6.4.25/main.c pythia6.4.25/pythia6_common_address.c)
-set ( FSRCS pythia6.4.25/pythia6_common_block_address.F pythia6.4.25/tpythia6_called_from_cc.F pythia6.4.25/pythia-6.4.25.f pythia6.4.25/pydummy_6.4.25.f)
+set ( FSRCS pythia6.4.25/pythia6_common_block_address.F pythia6.4.25/tpythia6_called_from_cc.F pythia6.4.25/pythia-6.4.25.f pythia6.4.25/pydummy_6.4.25.f POWHEG-BOX/setup-PYTHIA-lhef.f POWHEG-BOX/lhefread.f)
if( ALICE_TARGET STREQUAL "win32gcc")
--- /dev/null
+c -*-Fortran-*-
+
+ integer maxpup
+ parameter(maxpup=100)
+ integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
+ double precision ebmup,xsecup,xerrup,xmaxup
+ common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
+ & idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
+ & xmaxup(maxpup),lprup(maxpup)
+ integer maxnup
+ parameter (maxnup=500)
+ integer nup,idprup,idup,istup,mothup,icolup
+ double precision xwgtup,scalup,aqedup,aqcdup,pup,vtimup,spinup
+ common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,
+ & idup(maxnup),istup(maxnup),mothup(2,maxnup),
+ & icolup(2,maxnup),pup(5,maxnup),vtimup(maxnup),
+ & spinup(maxnup)
+ save /hepeup/
--- /dev/null
+c -*-Fortran-*-
+
+ integer NMXHEP,NEVHEP,NHEP,ISTHEP,IDHEP,
+ & JMOHEP,JDAHEP
+ double precision phep,vhep
+ PARAMETER (NMXHEP=4000)
+ COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
+ & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
+ save /HEPEVT/
--- /dev/null
+c...lhefheader(nlf)
+c...reads initialization information from a les houches events file on unit nlf.
+ subroutine lhefreadhdr(nlf)
+ implicit none
+ integer nlf
+ character * 100 string
+ integer ipr
+ include 'LesHouches.h'
+ 1 read(nlf,fmt='(a)',err=998,end=998) string
+ if(string(1:5).eq.'<init') then
+ read(nlf,*) idbmup(1),idbmup(2),ebmup(1),ebmup(2),
+ & pdfgup(1),pdfgup(2),pdfsup(1),pdfsup(2),idwtup,nprup
+ do ipr=1,nprup
+ read(nlf,*) xsecup(ipr),xerrup(ipr),xmaxup(ipr),
+ & lprup(ipr)
+ enddo
+ goto 999
+ else
+ goto 1
+ endif
+ 998 write(*,*) 'lhefreadhdr: could not find <init> data'
+ call exit(1)
+ 999 end
+
+
+c...reads event information from a les houches events file on unit nlf.
+ subroutine lhefreadev(nlf)
+ implicit none
+ integer nlf
+ character * 100 string
+ include 'LesHouches.h'
+ integer i,j
+ 1 continue
+ string=' '
+ read(nlf,fmt='(a)',err=777,end=666) string
+ if(string.eq.'</LesHouchesEvents>') then
+ goto 998
+ endif
+ if(string(1:6).eq.'<event') then
+c on error try next event. The error may be cause by merging
+c truncated event files. On EOF return with no event found
+ read(nlf,*,end=998,err=1)nup,idprup,xwgtup,scalup,aqedup,aqcdup
+ do i=1,nup
+ read(nlf,*,end=998,err=1) idup(i),istup(i),mothup(1,i),
+ & mothup(2,i),icolup(1,i),icolup(2,i),(pup(j,i),j=1,5),
+ & vtimup(i),spinup(i)
+ enddo
+ goto 999
+ else
+ goto 1
+ endif
+c no event found:
+ 777 continue
+ print *,"Error in reading"
+ print *,string
+ stop
+ 666 continue
+ print *,"reached EOF"
+ print *,string
+ stop
+ 998 continue
+ print *,"read </LesHouchesEvents>"
+ nup=0
+ 999 end
--- /dev/null
+ subroutine UPINIT
+ implicit none
+ include 'hepevt.h'
+ include 'LesHouches.h'
+ double precision parp,pari
+ integer mstp,msti
+ common/pypars/mstp(200),parp(200),msti(200),pari(200)
+ integer MSTU,MSTJ
+ double precision PARU,PARJ
+ COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
+ integer MDCY,MDME,KFDP
+ double precision brat
+ COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
+ integer pycomp
+ external pycomp
+ integer maxev
+ common/mcmaxev/maxev
+ nevhep=0
+c read the header first, so lprup is set
+ call lhefreadhdr(97)
+ end
+
+ subroutine UPEVNT
+ implicit none
+ call lhefreadev(97)
+ end
+
+ subroutine UPVETO
+c pythia routine to abort event
+ end
kPyCharmppMNR, kPyCharmppMNRwmi, kPyD0ppMNR, kPyDPlusppMNR, kPyDPlusStrangeppMNR,
kPyBeautyppMNR, kPyBeautyppMNRwmi, kPyBeautyJets, kPyW, kPyZ, kPyLambdacppMNR, kPyMbMSEL1,
kPyOldUEQ2ordered, kPyOldUEQ2ordered2, kPyOldPopcorn,
- kPyLhwgMb, kPyMbDefault, kPyMbAtlasTuneMC09, kPyMBRSingleDiffraction, kPyMBRDoubleDiffraction, kPyMBRCentralDiffraction
+ kPyLhwgMb, kPyMbDefault, kPyMbAtlasTuneMC09, kPyMBRSingleDiffraction, kPyMBRDoubleDiffraction, kPyMBRCentralDiffraction, kPyJetsPWHG, kPyCharmPWHG, kPyBeautyPWHG
}
Process_t;
C...New example: handles a standard Les Houches Events File.
- SUBROUTINE UPINIT
+ SUBROUTINE UPINITD
C...Double precision and integer declarations.
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
C...New example: handles a standard Les Houches Events File.
- SUBROUTINE UPEVNT
+ SUBROUTINE UPEVNTD
C...Double precision and integer declarations.
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
C...IVETO = 0 : retain current event and generate in full;
C... = 1 : abort generation of current event and move to next.
- SUBROUTINE UPVETO(IVETO)
+ SUBROUTINE UPVETOD(IVETO)
C...HEPEVT commonblock.
PARAMETER (NMXHEP=4000)