#include "AliFastGlauber.h"
#include "AliQuenchingWeights.h"
#include "AliOmegaDalitz.h"
+#include "AliDecayerExodus.h"
#include "AliLog.h"
#include "TVector3.h"
#include "TLorentzVector.h"
fGlauber(0),
fQuenchingWeights(0),
fItune(-1),
- fOmegaDalitz()
+ fOmegaDalitz(),
+ fExodus()
{
// Default Constructor
//
fGlauber(0),
fQuenchingWeights(0),
fItune(-1),
- fOmegaDalitz()
+ fOmegaDalitz(),
+ fExodus()
{
// Copy Constructor
Int_t i;
SetMSEL(1);
// Pythia Tune A (CDF)
//
- SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
- SetMSTP(82,4); // Double Gaussian Model
- SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
- SetPARP(84,0.4); // Core radius
- SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
- SetPARP(86,0.95); // Regulates gluon prod. mechanism
- SetPARP(89,1800.); // [GeV] Ref. energy
- SetPARP(90,0.25); // 2*epsilon (exponent in power law)
- break;
+ if (fItune < 0) {
+ SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
+ SetMSTP(82,4); // Double Gaussian Model
+ SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
+ SetPARP(84,0.4); // Core radius
+ SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
+ SetPARP(86,0.95); // Regulates gluon prod. mechanism
+ SetPARP(89,1800.); // [GeV] Ref. energy
+ SetPARP(90,0.25); // 2*epsilon (exponent in power law)
+ }
+ break;
case kPyDirectGamma:
SetMSEL(10);
break;
SetMSTP(71,1); //Final QCD & QED showers on
break;
+ case kPyZgamma:
+ //Inclusive production of Z
+ SetMSEL(0);
+ //f fbar -> Z/gamma
+ SetMSUB(1,1);
+ // Initial/final parton shower on (Pythia default)
+ // With parton showers on we are generating "Z inclusive process"
+ SetMSTP(61,1); //Initial QCD & QED showers on
+ SetMSTP(71,1); //Final QCD & QED showers on
+ break;
case kPyMBRSingleDiffraction:
case kPyMBRDoubleDiffraction:
case kPyMBRCentralDiffraction:
case kPyCharmPWHG:
case kPyBeautyPWHG:
+ case kPyWPWHG:
// number of warnings printed on the shell
SetMSTU(26,20);
}
//
SetMSTP(41,1); // all resonance decays switched on
- if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG) {
+ if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG || process == kPyWPWHG) {
Initialize("USER","","",0.);
} else {
Initialize("CMS",fProjectile,fTarget,fEcms);
}
fOmegaDalitz.Init();
+ fExodus.Init();
}
Int_t AliPythia::CheckedLuComp(Int_t kf)
}
}
}
+
+
+//
+// Replace all dalitz(pi0,eta,omega,eta',phi) and resonance(rho,omega,phi,jpsi) decays with the correct matrix element decays
+// for di-electron cocktail calculations
+//
+
+
+void AliPythia::PizeroDalitz()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 111) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+ if ((fPyjets->K[1][fd] != 22) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11) )
+ continue;
+ TLorentzVector pizero(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &pizero);
+ for (Int_t j = 0; j < 3; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_pion())[2-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::EtaDalitz()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 221) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+ if ((fPyjets->K[1][fd] != 22) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+ continue;
+ TLorentzVector eta(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &eta);
+ for (Int_t j = 0; j < 3; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_eta())[2-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::RhoDirect()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 113) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 1) continue;
+ if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+ continue;
+ TLorentzVector rho(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &rho);
+ for (Int_t j = 0; j < 2; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_rho())[1-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::OmegaDalitz()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 223) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+ if ((fPyjets->K[1][fd] != 111) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+ continue;
+ TLorentzVector omegadalitz(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &omegadalitz);
+ for (Int_t j = 0; j < 3; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_omega_dalitz())[2-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::OmegaDirect()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 223) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 1) continue;
+ if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+ continue;
+ TLorentzVector omegadirect(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &omegadirect);
+ for (Int_t j = 0; j < 2; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_omega())[1-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::EtaprimeDalitz()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 331) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+ if ((fPyjets->K[1][fd] != 22) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+ continue;
+ TLorentzVector etaprime(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &etaprime);
+ for (Int_t j = 0; j < 3; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_etaprime())[2-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::PhiDalitz()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 333) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+ if ((fPyjets->K[1][fd] != 221) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+ continue;
+ TLorentzVector phidalitz(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &phidalitz);
+ for (Int_t j = 0; j < 3; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_phi_dalitz())[2-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+
+void AliPythia::PhiDirect()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 333) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 1) continue;
+ if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+ continue;
+ TLorentzVector phi(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &phi);
+ for (Int_t j = 0; j < 2; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_phi())[1-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+void AliPythia::JPsiDirect()
+{
+
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 443) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 1) continue;
+ if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+ continue;
+ TLorentzVector jpsi(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+ fExodus.Decay(pdg, &jpsi);
+ for (Int_t j = 0; j < 2; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fExodus.Products_jpsi())[1-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}
+
+