]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliPythia.cxx
ALIROOT-5575 Missing protection in AliAnalysisTaskFiltereTree
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
index 9b5e6819bac3627f8709f9d3e0b8d9c1fcb94915..d8a267940827c7d656afbdffa6931eb5e725907c 100644 (file)
@@ -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;
@@ -375,15 +378,17 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        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;
@@ -613,6 +618,16 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       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:
@@ -635,6 +650,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     
       case kPyCharmPWHG:
       case kPyBeautyPWHG:
+      case kPyWPWHG:
       //    number of warnings printed on the shell
       SetMSTU(26,20);
             
@@ -655,12 +671,13 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     }
 //  
     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)
@@ -1535,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];
+        }
+     }
+  }
+}
+
+