Decayer with alice decay options using Pythia8
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Dec 2008 13:29:42 +0000 (13:29 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Dec 2008 13:29:42 +0000 (13:29 +0000)
PYTHIA8/AliDecayerPythia8.cxx [new file with mode: 0644]
PYTHIA8/AliDecayerPythia8.h [new file with mode: 0644]
PYTHIA8/AliPythia8LinkDef.h
PYTHIA8/libAliPythia8.pkg

diff --git a/PYTHIA8/AliDecayerPythia8.cxx b/PYTHIA8/AliDecayerPythia8.cxx
new file mode 100644 (file)
index 0000000..c3eff97
--- /dev/null
@@ -0,0 +1,479 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+// Implementation of AliDecayer using Pythia8
+// Author: andreas.morsch@cern.ch
+#include <TMath.h>
+#include <TPDGCode.h>
+#include <TPythia8.h>
+#include "AliDecayerPythia8.h"
+#include "ParticleData.h"
+
+ClassImp(AliDecayerPythia8)
+
+AliDecayerPythia8::AliDecayerPythia8():
+    TPythia8Decayer(),
+    fDecay(kAll),
+    fHeavyFlavour(kTRUE)
+{
+    // Constructor
+}
+
+void AliDecayerPythia8::ForceDecay()
+{
+  // 
+// Force a particle decay mode
+// Switch heavy flavour production off if requested
+    if (!fHeavyFlavour) SwitchOffHeavyFlavour();
+//
+    Decay_t decay = fDecay;
+    TPythia8::Instance()->ReadString("HadronLevel:Decay = on");
+    
+    if (decay == kNoDecayHeavy) return;
+
+//
+// select mode    
+    Int_t products[2];
+    Int_t mult[2];
+    Int_t products1[3];
+    Int_t mult1[3];
+    
+    switch (decay) 
+    {
+    case kHardMuons:
+       products1[0] =     13;
+       products1[1] =    443;
+       products1[2] = 100443;
+       mult1[0] = 1;
+       mult1[1] = 1;
+       mult1[2] = 1;
+       ForceParticleDecay(  511, products1, mult1, 3); 
+       ForceParticleDecay(  521, products1, mult1, 3); 
+       ForceParticleDecay(  531, products1, mult1, 3); 
+       ForceParticleDecay( 5122, products1, mult1, 3); 
+       ForceParticleDecay( 5132, products1, mult1, 3); 
+       ForceParticleDecay( 5232, products1, mult1, 3); 
+       ForceParticleDecay( 5332, products1, mult1, 3); 
+       ForceParticleDecay( 100443, 443, 1);  // Psi'  -> J/Psi X    
+       ForceParticleDecay(    443,  13, 2);  // J/Psi -> mu+ mu-
+       ForceParticleDecay(  411,13,1); // D+/-     
+       ForceParticleDecay(  421,13,1); // D0     
+       ForceParticleDecay(  431,13,1); // D_s     
+       ForceParticleDecay( 4122,13,1); // Lambda_c    
+       ForceParticleDecay( 4132,13,1); // Xsi_c     
+       ForceParticleDecay( 4232,13,1); // Sigma_c 
+       ForceParticleDecay( 4332,13,1); // Omega_c     
+       break;
+   case kChiToJpsiGammaToMuonMuon:
+       products[0] =    443;
+       products[1] =     22;
+       mult[0] = 1;
+       mult[1] = 1;
+       ForceParticleDecay( 20443, products, mult, 2);   // Chi_1c  -> J/Psi  Gamma
+       ForceParticleDecay(   445, products, mult, 2);   // Chi_2c  -> J/Psi  Gamma
+       ForceParticleDecay(   443, 13, 2);               // J/Psi -> mu+ mu-
+       break;
+    case kChiToJpsiGammaToElectronElectron:
+       products[0] =    443;
+       products[1] =     22;
+       mult[0] = 1;
+       mult[1] = 1;
+       ForceParticleDecay( 20443, products, mult, 2);   // Chi_1c  -> J/Psi  Gamma
+       ForceParticleDecay(   445, products, mult, 2);   // Chi_2c  -> J/Psi  Gamma
+       ForceParticleDecay(   443, 11, 2);               // J/Psi -> e+ e-
+       break;
+
+    case kBSemiMuonic:
+       ForceParticleDecay(  511,13,1); // B0     
+       ForceParticleDecay(  521,13,1); // B+/-     
+       ForceParticleDecay(  531,13,1); // B_s     
+       ForceParticleDecay( 5122,13,1); // Lambda_b    
+       ForceParticleDecay( 5132,13,1); // Xsi_b    
+       ForceParticleDecay( 5232,13,1); // Sigma_b    
+       ForceParticleDecay( 5332,13,1); // Omega_b    
+       break;
+    case kSemiMuonic:
+       ForceParticleDecay(  411,13,1); // D+/-     
+       ForceParticleDecay(  421,13,1); // D0     
+       ForceParticleDecay(  431,13,1); // D_s     
+       ForceParticleDecay( 4122,13,1); // Lambda_c    
+       ForceParticleDecay( 4132,13,1); // Xsi_c     
+       ForceParticleDecay( 4232,13,1); // Sigma_c 
+       ForceParticleDecay( 4332,13,1); // Omega_c     
+       ForceParticleDecay(  511,13,1); // B0     
+       ForceParticleDecay(  521,13,1); // B+/-     
+       ForceParticleDecay(  531,13,1); // B_s     
+       ForceParticleDecay( 5122,13,1); // Lambda_b    
+       ForceParticleDecay( 5132,13,1); // Xsi_b    
+       ForceParticleDecay( 5232,13,1); // Sigma_b    
+       ForceParticleDecay( 5332,13,1); // Omega_b    
+       break;
+    case kDiMuon:
+       ForceParticleDecay(  113,13,2); // rho
+       ForceParticleDecay(  221,13,2); // eta
+       ForceParticleDecay(  223,13,2); // omega
+       ForceParticleDecay(  333,13,2); // phi
+       ForceParticleDecay(  443,13,2); // J/Psi
+       ForceParticleDecay(100443,13,2);// Psi'
+       ForceParticleDecay(  553,13,2); // Upsilon
+       ForceParticleDecay(100553,13,2);// Upsilon'
+       ForceParticleDecay(200553,13,2);// Upsilon''
+       break;
+    case kBSemiElectronic:
+       ForceParticleDecay(  511,11,1); // B0     
+       ForceParticleDecay(  521,11,1); // B+/-     
+       ForceParticleDecay(  531,11,1); // B_s     
+       ForceParticleDecay( 5122,11,1); // Lambda_b     
+       ForceParticleDecay( 5132,11,1); // Xsi_b    
+       ForceParticleDecay( 5232,11,1); // Sigma_b    
+       ForceParticleDecay( 5332,11,1); // Omega_b    
+       break;
+    case kSemiElectronic:
+       ForceParticleDecay(  411,11,1); // D+/-     
+       ForceParticleDecay(  421,11,1); // D0     
+       ForceParticleDecay(  431,11,1); // D_s     
+       ForceParticleDecay( 4122,11,1); // Lambda_c     
+       ForceParticleDecay( 4132,11,1); // Xsi_c     
+       ForceParticleDecay( 4232,11,1); // Sigma_c 
+       ForceParticleDecay( 4332,11,1); // Omega_c     
+       ForceParticleDecay(  511,11,1); // B0     
+       ForceParticleDecay(  521,11,1); // B+/-     
+       ForceParticleDecay(  531,11,1); // B_s     
+       ForceParticleDecay( 5122,11,1); // Lambda_b     
+       ForceParticleDecay( 5132,11,1); // Xsi_b    
+       ForceParticleDecay( 5232,11,1); // Sigma_b    
+       ForceParticleDecay( 5332,11,1); // Omega_b    
+       break;
+    case kDiElectron:
+       ForceParticleDecay(  113,11,2); // rho
+       ForceParticleDecay(  333,11,2); // phi
+       ForceParticleDecay(  221,11,2); // eta
+       ForceParticleDecay(  223,11,2); // omega
+       ForceParticleDecay(  443,11,2); // J/Psi
+       ForceParticleDecay(100443,11,2);// Psi'
+       ForceParticleDecay(  553,11,2); // Upsilon
+       ForceParticleDecay(100553,11,2);// Upsilon'
+       ForceParticleDecay(200553,11,2);// Upsilon''
+       break;
+    case kBJpsiDiMuon:
+
+       products[0] =    443;
+       products[1] = 100443;
+       mult[0] = 1;
+       mult[1] = 1;
+
+       ForceParticleDecay(  511, products, mult, 2); // B0   -> J/Psi (Psi') X   
+       ForceParticleDecay(  521, products, mult, 2); // B+/- -> J/Psi (Psi') X     
+       ForceParticleDecay(  531, products, mult, 2); // B_s  -> J/Psi (Psi') X     
+       ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi') X 
+       ForceParticleDecay( 100443, 443, 1);          // Psi'  -> J/Psi X    
+       ForceParticleDecay(    443,13,2);             // J/Psi -> mu+ mu-   
+       break;
+    case kBPsiPrimeDiMuon:
+       ForceParticleDecay(  511,100443,1); // B0     
+       ForceParticleDecay(  521,100443,1); // B+/-     
+       ForceParticleDecay(  531,100443,1); // B_s     
+       ForceParticleDecay( 5122,100443,1); // Lambda_b 
+       ForceParticleDecay(100443,13,2);    // Psi'   
+       break;
+    case kBJpsiDiElectron:
+       ForceParticleDecay(  511,443,1); // B0     
+       ForceParticleDecay(  521,443,1); // B+/-     
+       ForceParticleDecay(  531,443,1); // B_s     
+       ForceParticleDecay( 5122,443,1); // Lambda_b
+       ForceParticleDecay(  443,11,2);  // J/Psi    
+       break;
+    case kBJpsi:
+       ForceParticleDecay(  511,443,1); // B0     
+       ForceParticleDecay(  521,443,1); // B+/-     
+       ForceParticleDecay(  531,443,1); // B_s     
+       ForceParticleDecay( 5122,443,1); // Lambda_b
+       break;
+    case kBPsiPrimeDiElectron:
+       ForceParticleDecay(  511,100443,1); // B0     
+       ForceParticleDecay(  521,100443,1); // B+/-     
+       ForceParticleDecay(  531,100443,1); // B_s     
+       ForceParticleDecay( 5122,100443,1); // Lambda_b 
+       ForceParticleDecay(100443,11,2);   // Psi'   
+       break;
+    case kPiToMu:
+       ForceParticleDecay(211,13,1); // pi->mu     
+       break;
+    case kKaToMu:
+       ForceParticleDecay(321,13,1); // K->mu     
+       break;
+    case kAllMuonic:
+       ForceParticleDecay(211,13,1); // pi->mu
+       ForceParticleDecay(321,13,1); // K->mu    
+       break;
+    case kWToMuon:
+        ForceParticleDecay(  24, 13,1); // W -> mu
+       break;
+    case kWToCharm:
+        ForceParticleDecay(   24, 4,1); // W -> c
+       break;
+    case kWToCharmToMuon:
+        ForceParticleDecay(   24, 4,1); // W -> c
+       ForceParticleDecay(  411,13,1); // D+/- -> mu
+       ForceParticleDecay(  421,13,1); // D0  -> mu
+       ForceParticleDecay(  431,13,1); // D_s  -> mu
+       ForceParticleDecay( 4122,13,1); // Lambda_c
+       ForceParticleDecay( 4132,13,1); // Xsi_c
+       ForceParticleDecay( 4232,13,1); // Sigma_c
+       ForceParticleDecay( 4332,13,1); // Omega_c
+       break;
+    case kZDiMuon:
+        ForceParticleDecay(  23, 13,2); // Z -> mu+ mu-
+       break;
+    case kZDiElectron:
+        ForceParticleDecay(  23, 11,2); // Z -> e+ e-
+       break;
+    case kHadronicD:
+       ForceHadronicD(1);
+       break;
+    case kHadronicDWithout4Bodies:
+       ForceHadronicD(0);
+       break;
+    case kPhiKK:
+       ForceParticleDecay(333,321,2); // Phi->K+K-
+       break;
+    case kOmega:
+//     ForceOmega();
+    case kAll:
+       break;
+    case kNoDecay:
+       TPythia8::Instance()->ReadString("HadronLevel:Decay = off");
+       break;
+    case kNoDecayHeavy:
+    case kNeutralPion:
+       break;
+    }
+}
+
+Float_t AliDecayerPythia8::GetPartialBranchingRatio(Int_t ipart)
+{
+    // Get the partial branching ration for the forced decay channels
+    
+    Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
+    Pythia8::ParticleDataTable table = thePythia->particleData;
+    Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(ipart);
+    Pythia8::DecayTable  decays      = pd->decay;
+    
+    
+    Int_t nc = decays.size();
+    Float_t br = 0.;
+//
+//  Loop over decay channels
+    for (Int_t ic = 0; ic < nc; ic++) {
+       Pythia8::DecayChannel& decCh = decays[ic];
+       for (Int_t i = 1; i <= decCh.multiplicity(); i++) {
+           br += decCh.bRatio();
+       }
+    }
+    return (br);
+}
+
+
+Int_t AliDecayerPythia8::CountProducts(Pythia8::DecayChannel& channel , Int_t particle)
+{
+// Count decay products of a given type
+    Int_t np = 0;
+    for (Int_t i = 1; i <= channel.multiplicity(); i++) {
+       if (TMath::Abs(channel.product(i)) == particle) np++;
+    }
+    return np;
+}
+
+
+
+void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
+{
+//
+//  Force decay of particle into products with multiplicity mult
+
+    Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
+    Pythia8::ParticleDataTable table = thePythia->particleData;
+    Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(particle);
+    pd->setMayDecay(true);
+    Pythia8::DecayTable  decays      = pd->decay;
+    
+    
+    Int_t nc = decays.size();
+//
+//  Loop over decay channels
+    for (Int_t ic = 0; ic < nc; ic++) {
+       Pythia8::DecayChannel& decCh = decays[ic];
+       if (CountProducts(decCh, product) >= mult) {
+           decCh.onMode(1);
+       } else {
+           decCh.onMode(0);
+       }
+    }
+}
+
+void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
+{
+//
+//  Force decay of particle into products with multiplicity mult
+
+    Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
+    Pythia8::ParticleDataTable table = thePythia->particleData;
+    Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(particle);
+    pd->setMayDecay(true);
+    Pythia8::DecayTable decays       = pd->decay;
+    
+    Int_t nc = decays.size();
+//
+//  Loop over decay channels
+    for (Int_t ic = 0; ic < nc; ic++) {
+       Int_t nprod = 0;
+       Pythia8::DecayChannel& decCh = decays[ic];
+
+       for (Int_t i = 0; i < npart; i++) {
+           nprod += (CountProducts(decCh, products[i]) >= mult[i]);
+       }
+       
+       if (nprod) {
+           decCh.onMode(1);
+       } else {
+           decCh.onMode(0);
+       }
+    }
+}
+
+
+Float_t AliDecayerPythia8::GetLifetime(Int_t kf)
+{
+    // Return lifetime of particle
+    Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
+    Pythia8::ParticleDataTable table = thePythia->particleData;
+    Float_t tau = table.tau0(kf);
+    return ( tau);
+}
+
+void  AliDecayerPythia8::SwitchOffHeavyFlavour()
+{
+    // Switch off heavy flavour production
+    //
+// Maximum number of quark flavours used in pdf 
+    TPythia8::Instance()->ReadString("PDFinProcess:nQuarkIn = 3");
+// Maximum number of flavors that can be used in showers
+    TPythia8::Instance()->ReadString("SpaceShower:nQuarkIn = 3");
+    TPythia8::Instance()->ReadString("TimeShower:nGammaToQuark = 3");
+    TPythia8::Instance()->ReadString("TimeShower:nGluonToQuark = 3");
+}
+
+
+void AliDecayerPythia8::ForceHadronicD(Int_t optUse4Bodies)
+{
+//
+// Force golden D decay modes
+//
+    const Int_t kNHadrons = 5;
+    Int_t hadron[kNHadrons] = {411,  421, 431, 4112, 4122};
+
+    // for D+ -> K0* (-> K- pi+) pi+
+    Int_t iKstar0    =  313;
+    Int_t iKstarbar0 = -313;
+    Int_t products[2] = {kKPlus, kPiMinus}, mult[2] = {1, 1};
+    ForceParticleDecay(iKstar0, products, mult, 2);
+    // for Ds -> Phi pi+
+    Int_t iPhi = 333;
+    ForceParticleDecay(iPhi, kKPlus, 2); // Phi->K+K-
+    // for D0 -> rho0 pi+ k-
+    Int_t iRho0=113;
+    ForceParticleDecay(iRho0, kPiPlus, 2); // Rho0->pi+pi-
+    // for Lambda_c -> Delta++ K-
+    Int_t iDeltaPP = 2224;
+    Int_t productsD[2] = {kProton, kPiPlus}, multD[2] = {1, 1};
+    ForceParticleDecay(iDeltaPP, productsD, multD, 2);
+
+
+    Int_t decayP1[kNHadrons][4] = 
+       { 
+           {kKMinus, kPiPlus,    kPiPlus, 0},
+           {kKMinus, kPiPlus,    0      , 0},
+           {kKPlus , iKstarbar0, 0      , 0},
+           {-1     , -1        , -1     , -1},
+           {kProton, iKstarbar0, 0      , 0}
+       };
+    Int_t decayP2[kNHadrons][4] = 
+       { 
+           {iKstarbar0, kPiPlus, 0      , 0},
+           {kKMinus   , kPiPlus, kPiPlus, kPiMinus},
+           {iPhi      , kPiPlus, 0      , 0},
+           {-1        , -1     , -1     , -1},
+           {iDeltaPP  , kKMinus, 0      , 0}
+       };
+    Int_t decayP3[kNHadrons][4] = 
+       { 
+           {-1        , -1     , -1      , -1},
+           {kKMinus   , kPiPlus, iRho0   , 0 },
+           {-1        , -1     , -1      , -1},
+           {-1        , -1     , -1      , -1},
+           {kProton   , kKMinus, kPiPlus , 0}
+       };
+    if(optUse4Bodies==0){
+       for(Int_t iDau=0;iDau<4;iDau++){
+           decayP2[1][iDau]=-1;
+           decayP3[1][iDau]=-1;
+       }
+    }  
+
+
+    Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
+    Pythia8::ParticleDataTable table = thePythia->particleData;
+    
+    for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
+    {
+       Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(hadron[ihadron]);
+       pd->setMayDecay(true);
+       Pythia8::DecayTable  decays      = pd->decay;
+       
+       
+       for (Int_t ic = 0; ic < decays.size(); ic++) {
+           Pythia8::DecayChannel& decCh = decays[ic];
+           if ((
+                   decCh.product(0) == decayP1[ihadron][0] &&
+                   decCh.product(1) == decayP1[ihadron][1] &&
+                   decCh.product(2) == decayP1[ihadron][2] &&
+                   decCh.product(3) == decayP1[ihadron][3] &&
+                   decCh.product(4) == 0
+                   ) 
+               || (
+                   decCh.product(0) == decayP2[ihadron][0] &&
+                   decCh.product(1) == decayP2[ihadron][1] &&
+                   decCh.product(2) == decayP2[ihadron][2] &&
+                   decCh.product(3) == decayP2[ihadron][3] &&
+                   decCh.product(4) == 0
+                   ) 
+               || (
+                   decCh.product(0) == decayP3[ihadron][0] &&
+                   decCh.product(1) == decayP3[ihadron][1] &&
+                   decCh.product(2) == decayP3[ihadron][2] &&
+                   decCh.product(3) == decayP3[ihadron][3] &&
+                   decCh.product(4) == 0
+                   ))
+           {
+               decCh.onMode(1);
+           } else {
+               decCh.onMode(0);
+           } // selected channel ?
+       } // decay channels
+    } // hadrons
+}
+
diff --git a/PYTHIA8/AliDecayerPythia8.h b/PYTHIA8/AliDecayerPythia8.h
new file mode 100644 (file)
index 0000000..54c9708
--- /dev/null
@@ -0,0 +1,43 @@
+#ifndef ALIDECAYERPYTHIA8_H
+#define ALIDECAYERPYTHIA8_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$*/
+
+// Implementation of TVirtualMCDecayer using Pythia8
+// Author: andreas.morsch@cern.ch
+
+#include <TPythia8Decayer.h>
+#include <ParticleData.h>
+#include "AliDecayer.h"
+
+class AliDecayerPythia8 : public TPythia8Decayer {
+ public:
+  AliDecayerPythia8();
+  virtual ~AliDecayerPythia8(){;}
+  virtual void    SetForceDecay(Decay_t decay) {fDecay=decay;}
+  virtual void    SetForceDecay(Int_t decay)   {SetForceDecay((Decay_t) decay);}
+  virtual void    ForceDecay();
+  virtual Float_t GetPartialBranchingRatio(Int_t ipart);
+  virtual void    HeavyFlavourOff() {fHeavyFlavour = kFALSE;}
+  virtual Float_t GetLifetime(Int_t kf);
+ private:
+  Int_t    CountProducts(Pythia8::DecayChannel& decCh, Int_t particle);
+  void     ForceParticleDecay(Int_t particle, Int_t product, Int_t mult);
+  void     ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart);
+  void     SwitchOffHeavyFlavour();
+  void     ForceHadronicD(Int_t optUser4Bodies = 1);
+ protected:
+  Decay_t     fDecay;           //  Forced decay mode
+  Bool_t      fHeavyFlavour;    //! Flag for heavy flavors
+  ClassDef(AliDecayerPythia8, 1) // Particle Decayer using Pythia8
+};
+#endif
+
+
+
+
+
+
+
index 46111e8..089c96a 100644 (file)
@@ -3,4 +3,5 @@
 #pragma link off all classes;
 #pragma link off all functions;
 #pragma link C++ class AliPythia8+;
+#pragma link C++ class AliDecayerPythia8+;
 #endif
index 1072fbd..ab6c42b 100644 (file)
@@ -1,4 +1,4 @@
-SRCS= AliPythia8.cxx
+SRCS= AliPythia8.cxx AliDecayerPythia8.cxx
 HDRS:= $(SRCS:.cxx=.h) 
 DHDR:= AliPythia8LinkDef.h     
 EINCLUDE:= $(PYTHIA8)/include PYTHIA6