#include "AliMC.h"
#include "AliLog.h"
#include "PyquenCommon.h"
-#include "AliLog.h"
ClassImp(AliGenPythia)
fFlavorSelect(0),
fXsection(0.),
fPythia(0),
+ fWeightPower(0.),
fPtHardMin(0.),
fPtHardMax(1.e4),
fYHardMin(-1.e10),
fYHardMax(1.e10),
fGinit(1),
fGfinal(1),
+ fCRoff(0),
fHadronisation(1),
fPatchOmegaDalitz(0),
+ fDecayerExodus(0),
fNpartons(0),
fReadFromFile(0),
+ fReadLHEF(0),
fQuench(0),
fQhat(0.),
fLength(0.),
fFeedDownOpt(kTRUE),
fFragmentation(kTRUE),
fSetNuclei(kFALSE),
+ fUseNuclearPDF(kFALSE),
+ fUseLorentzBoost(kTRUE),
fNewMIS(kFALSE),
fHFoff(kFALSE),
fNucPdf(0),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerY(999.),
+ fTriggerEtaMin(0.9),
fTriggerMinPt(-1),
fTriggerMaxPt(1000),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
+ fTriggerMultiplicityEtaMin(0),
+ fTriggerMultiplicityEtaMax(0),
fTriggerMultiplicityPtMin(0),
fCountMode(kCountAll),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
fFlavorSelect(0),
fXsection(0.),
fPythia(0),
+ fWeightPower(0.),
fPtHardMin(0.),
fPtHardMax(1.e4),
fYHardMin(-1.e10),
fYHardMax(1.e10),
fGinit(kTRUE),
fGfinal(kTRUE),
+ fCRoff(kFALSE),
fHadronisation(kTRUE),
- fPatchOmegaDalitz(0),
+ fPatchOmegaDalitz(0),
+ fDecayerExodus(0),
fNpartons(0),
fReadFromFile(kFALSE),
+ fReadLHEF(0),
fQuench(kFALSE),
fQhat(0.),
fLength(0.),
fFeedDownOpt(kTRUE),
fFragmentation(kTRUE),
fSetNuclei(kFALSE),
+ fUseNuclearPDF(kFALSE),
+ fUseLorentzBoost(kTRUE),
fNewMIS(kFALSE),
fHFoff(kFALSE),
fNucPdf(0),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerY(999.),
+ fTriggerEtaMin(0.9),
fTriggerMinPt(-1),
fTriggerMaxPt(1000),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
+ fTriggerMultiplicityEtaMin(0),
+ fTriggerMultiplicityEtaMax(0),
fTriggerMultiplicityPtMin(0),
fCountMode(kCountAll),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
//
fParentWeight=1./Float_t(fNpart);
//
-
-
+ if (fWeightPower != 0)
+ fPythia->SetWeightPower(fWeightPower);
fPythia->SetCKIN(3,fPtHardMin);
fPythia->SetCKIN(4,fPtHardMax);
fPythia->SetCKIN(7,fYHardMin);
fPythia->SetCKIN(8,fYHardMax);
- if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
+ if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
+
+ if(fUseNuclearPDF)
+ fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
// Fragmentation?
if (fFragmentation) {
fPythia->SetMSTP(111,1);
fPythia->SetMSTP(61,fGinit);
// final state radiation
fPythia->SetMSTP(71,fGfinal);
+ //color reconnection strength
+ if(fCRoff==1)fPythia->SetMSTP(95,0);
// pt - kick
if (fPtKick > 0.) {
fPythia->SetMSTP(91,1);
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;
+ case kPyWPWHG:
case kPyW:
case kPyZ:
+ case kPyZgamma:
+ case kPyMBRSingleDiffraction:
+ case kPyMBRDoubleDiffraction:
+ case kPyMBRCentralDiffraction:
break;
}
//
//
//
AliGenMC::Init();
+
+ // Reset Lorentz boost if demanded
+ if(!fUseLorentzBoost) {
+ fDyBoost = 0;
+ Warning("Init","Demand to discard Lorentz boost.\n");
+ }
//
//
//
if (fSetNuclei) {
- fDyBoost = 0;
- Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
+ fDyBoost = 0;
+ Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
}
-
fPythia->SetPARJ(200, 0.0);
fPythia->SetPARJ(199, 0.0);
fPythia->SetPARJ(198, 0.0);
fPythia->SetMSTJ(44, 2); // option to run alpha_s
fPythia->SetPARJ(82, 1.); // Cut off for parton showers
}
+
+ if ( AliLog::GetDebugLevel("","AliGenPythia") >= 1 ) {
+ fPythia->Pystat(4);
+ fPythia->Pystat(5);
+ }
}
void AliGenPythia::Generate()
fPythia->DalitzDecays();
fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
}
- fPythia->Pyexec();
- }
+
+ else if (fDecayerExodus) {
+
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+ fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+ fPythia->Pyexec();
+ fPythia->OmegaDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+ fPythia->PizeroDalitz();
+ fPythia->PhiDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+ fPythia->EtaDalitz();
+ fPythia->EtaprimeDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+ fPythia->RhoDirect();
+ fPythia->OmegaDirect();
+ fPythia->PhiDirect();
+ fPythia->JPsiDirect();
+ }
+
+ fPythia->Pyexec();
+ }
fTrials++;
fPythia->ImportParticles(&fParticles,"All");
fProcess != kPyMbMSEL1 &&
fProcess != kPyW &&
fProcess != kPyZ &&
+ fProcess != kPyZgamma &&
fProcess != kPyCharmppMNRwmi &&
fProcess != kPyBeautyppMNRwmi &&
- fProcess != kPyBeautyJets) {
+ fProcess != kPyBeautyJets &&
+ fProcess != kPyWPWHG &&
+ 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) {
- TParticle* jet1 = (TParticle *) fParticles.At(6);
+ 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);
- if (!CheckTrigger(jet1, jet2)) {
+
+ if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
delete [] pParent;
return 0;
}
// eta cut
if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
continue;
+ //multiplicity check for a given eta range
+ if ((fTriggerMultiplicityEtaMin != fTriggerMultiplicityEtaMax) &&
+ (iparticle->Eta() < fTriggerMultiplicityEtaMin || iparticle->Eta() > fTriggerMultiplicityEtaMax))
+ continue;
// pt cut
if (iparticle->Pt() < fTriggerMultiplicityPtMin)
continue;
kf = CheckPDGCode(iparticle->GetPdgCode());
if (kf != fTriggerParticle) continue;
if (iparticle->Pt() == 0.) continue;
- if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+ if (TMath::Abs(iparticle->Y()) > fTriggerY) continue;
+ if (fTriggerEtaMin == fTriggerEta) {
+ if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+ } else {
+ if (iparticle->Eta() > fTriggerEta || iparticle->Eta() < fTriggerEtaMin) continue;
+ }
if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
triggered = kTRUE;
break;
Bool_t theChild=kFALSE;
Bool_t triggered=kFALSE;
Float_t y;
- Int_t pdg,mpdg,mpdgUpperFamily;
- for(i=0; i<np; i++) {
+ Int_t pdg, mpdg, mpdgUpperFamily;
+ for(i = 0; i < np; i++) {
partCheck = (TParticle*)fParticles.At(i);
pdg = partCheck->GetPdgCode();
if(TMath::Abs(pdg) == fFlavorSelect) { // quark
}
//Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
- if ( (fProcess == kPyW ||
+ if ( (
+ fProcess == kPyWPWHG ||
+ fProcess == kPyW ||
fProcess == kPyZ ||
+ fProcess == kPyZgamma ||
fProcess == kPyMbDefault ||
fProcess == kPyMb ||
fProcess == kPyMbAtlasTuneMC09 ||
return 0;
}
}
-
+
for (i = 0; i < np; i++) {
Int_t trackIt = 0;
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++;
{
// Treat protons as inside nuclei with mass numbers a1 and a2
- fAProjectile = a1;
- fATarget = a2;
- fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
- fSetNuclei = kTRUE;
+ fAProjectile = a1;
+ fATarget = a2;
+ fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
+ fUseNuclearPDF = kTRUE;
+ fSetNuclei = kTRUE;
}
// Number of primaries
fHeader->SetNProduced(fNprimaries);
//
+// Event weight
+ fHeader->SetEventWeight(fPythia->GetVINT(97));
+//
// 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];
((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
}
//
-// Store pt^hard
+// Store pt^hard and cross-section
((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
+ ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
+
+//
+// Store Event Weight
+ ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7)*fPythia->GetPARI(10));
+ // PARI(7) is 1 or -1, for weighted generation with accept/reject, e.g. POWHEG
+ // PARI(10) is a weight associated with reweighted generation, using Pyevwt
//
// 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];
return;
}
-Bool_t AliGenPythia::IsInBarrel(const Float_t eta) const
+Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
{
// Is particle in Central Barrel acceptance?
// etamin=-etamax
return kFALSE;
}
-Bool_t AliGenPythia::IsInEMCAL(const Float_t phi, const Float_t eta) const
+Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
{
// Is particle in EMCAL acceptance?
// phi in degrees, etamin=-etamax
return kFALSE;
}
-Bool_t AliGenPythia::IsInPHOS(const Float_t phi, const Float_t eta, const Int_t iparticle)
+Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
{
// Is particle in PHOS acceptance?
// Acceptance slightly larger considered.
TParticle * part = (TParticle *) fParticles.At(iPart);
Double_t E= part->Energy();
Double_t P= part->P();
- M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
+ Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
+ if(M2<0) return kFALSE;
+ M= TMath::Sqrt(M2);
}
Double_t Mmin, Mmax, wSD, wDD, wND;
return kFALSE;
}
-Bool_t AliGenPythia::CheckDetectorAcceptance(const Float_t phi, const Float_t eta, const Int_t iparticle)
+Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
{
// check the eta/phi correspond to the detectors acceptance
// iparticle is the index of the particle to be checked, for PHOS rotation case
else return kFALSE;
}
-Bool_t AliGenPythia::TriggerOnSelectedParticles(const Int_t np)
+Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
{
// Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
// implemented primaryly for kPyJets, but extended to any kind of process.
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff