Modifications needed to use EXODUS Decayer

[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx

diff --git a/PYTHIA6/AliPythia.cxx b/PYTHIA6/AliPythia.cxx
index 360cd8e59948cae1653b481bd63b2ccdbb5f80f7..d8a267940827c7d656afbdffa6931eb5e725907c 100644 (file)
--- a/PYTHIA6/AliPythia.cxx
+++ b/PYTHIA6/AliPythia.cxx
@@ -20,6 +20,7 @@
 #include "AliFastGlauber.h"
 #include "AliQuenchingWeights.h"
 #include "AliOmegaDalitz.h"
+#include "AliDecayerExodus.h"
 #include "AliLog.h"
 #include "TVector3.h"
 #include "TLorentzVector.h"
@@ -79,7 +80,8 @@ AliPythia::AliPythia():
     fGlauber(0),
     fQuenchingWeights(0),
     fItune(-1),
-    fOmegaDalitz()
+    fOmegaDalitz(),
+    fExodus()
 {
 // Default Constructor
 //
@@ -109,7 +111,8 @@ AliPythia::AliPythia(const AliPythia& pythia):
     fGlauber(0),
     fQuenchingWeights(0),
     fItune(-1),
-    fOmegaDalitz()
+    fOmegaDalitz(),
+    fExodus()
 {
     // Copy Constructor
     Int_t i;
@@ -674,6 +677,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       Initialize("CMS",fProjectile,fTarget,fEcms);
     }
     fOmegaDalitz.Init();
+    fExodus.Init();
 }
 
 Int_t AliPythia::CheckedLuComp(Int_t kf)
@@ -1548,3 +1552,235 @@ void AliPythia::DalitzDecays()
     }
   }
 }
+
+
+//
+// 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];
+        }
+     }
+  }
+}
+
+