fPatchOmegaDalitz(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),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
fPatchOmegaDalitz(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),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
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(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:
//
//
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
}
+
+ StdoutToAliDebug(1, fPythia->Pystat(4));
+ StdoutToAliDebug(1, fPythia->Pystat(5));
}
void AliGenPythia::Generate()
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) {
+ 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);
}
//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 ||
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));
+
//
// 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];
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;