]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliDecayerPythia.cxx
Removed the table directory component parameter for the TPC
[u/mrichter/AliRoot.git] / PYTHIA6 / AliDecayerPythia.cxx
index 1d59f966a96cb69d3ab761b059e4c904edeb2d3b..6a7d662f5498c4837518e84fbd5a9663eebe8550 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.18  2003/02/12 09:09:19  morsch
-Pylist removed.
-
-Revision 1.17  2003/01/31 16:54:38  morsch
-Use TPDCCode.h instead of AliPDG.
-
-Revision 1.16  2003/01/31 15:56:42  morsch
-Forcing of decay channels independent of order in decay table.
-
-Revision 1.15  2002/10/14 14:55:35  hristov
-Merging the VirtualMC branch to the main development branch (HEAD)
-
-Revision 1.14  2002/10/11 10:05:18  morsch
-pdg code for psi' corrected.
-
-Revision 1.10.6.3  2002/10/10 16:40:08  hristov
-Updating VirtualMC to v3-09-02
-
-Revision 1.13  2002/09/16 10:40:48  morsch
-Use correct pdg codes for Upsilon(2S) = 100553 and Upsilon(3S) = 200553.
-
-Revision 1.12  2002/06/05 14:05:46  morsch
-Decayer option kPhiKK for forced phi->K+K- decay added.
-
-Revision 1.11  2002/04/26 10:32:59  morsch
-Option kNoDecayHeavy added.
-
-Revision 1.10  2002/02/22 17:28:05  morsch
-ReadDecayTable() and WriteDecayTable() methods added.
-
-Revision 1.9  2001/12/20 10:37:13  morsch
-- Add omega forced decay.
-- Semileptonic decays for some more B and D baryons.
-
-Revision 1.8  2001/07/04 10:28:20  morsch
-Introduce GetLifetime(Int_T kf) method until functionality provided by
-TParticlePDG.
-
-Revision 1.7  2001/04/12 07:23:28  morsch
-Reactivate forcing option for dimuon and dielectron decay channels of phi (333).
-
-Revision 1.6  2001/03/27 10:53:26  morsch
-Save pythia default decay table at first initialization. Reload at each
-following Init() call.
-
-Revision 1.5  2001/03/09 13:04:06  morsch
-Decay_t moved to AliDecayer.h
-
-Revision 1.4  2001/01/30 09:23:11  hristov
-Streamers removed (R.Brun)
-
-Revision 1.3  2000/12/21 16:24:06  morsch
-Coding convention clean-up
-
-Revision 1.2  2000/09/12 13:58:45  morsch
-SetForceDcay(..) sets the member data storing the forced decay information.
-ForceDecay() executes the change of the decay table.
-
-Revision 1.1  2000/09/06 14:23:43  morsch
-Realisation of AliDecayer using Pythia6
-
-*/
+/* $Id$ */
 
 // Implementation of AliDecayer using Pythia
 // Method forwarding to the AliPythia instance.
@@ -96,10 +33,12 @@ ClassImp(AliDecayerPythia)
 #ifndef WIN32
 # define py1ent py1ent_
 # define opendecaytable opendecaytable_
+# define closedecaytable closedecaytable_
 # define type_of_call
 #else
 # define lu1ent PY1ENT
 # define opendecaytable OPENDECAYTABLE
+# define closedecaytable CLOSEDECAYTABLE
 # define type_of_call _stdcall
 #endif
 
@@ -109,6 +48,9 @@ extern "C" void type_of_call
 extern "C" void type_of_call 
           opendecaytable(Int_t&);
 
+extern "C" void type_of_call 
+          closedecaytable(Int_t&);
+
 Bool_t AliDecayerPythia::fgInit = kFALSE;
 
 
@@ -117,6 +59,7 @@ AliDecayerPythia::AliDecayerPythia()
 // Default Constructor
     fPythia=AliPythia::Instance();
     for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
+    ReadDecayTable();
 }
 
 void AliDecayerPythia::Init()
@@ -181,9 +124,38 @@ void AliDecayerPythia::ForceDecay()
 
 //
 // 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 kSemiMuonic:
        ForceParticleDecay(  411,13,1); // D+/-     
        ForceParticleDecay(  421,13,1); // D0     
@@ -205,10 +177,10 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,13,2); // omega
        ForceParticleDecay(  333,13,2); // phi
        ForceParticleDecay(  443,13,2); // J/Psi
-       ForceParticleDecay(100443,13,2); // Psi'
+       ForceParticleDecay(100443,13,2);// Psi'
        ForceParticleDecay(  553,13,2); // Upsilon
-       ForceParticleDecay(100553,13,2); // Upsilon'
-       ForceParticleDecay(200553,13,2); // Upsilon''
+       ForceParticleDecay(100553,13,2);// Upsilon'
+       ForceParticleDecay(200553,13,2);// Upsilon''
        break;
     case kSemiElectronic:
        ForceParticleDecay(  411,11,1); // D+/-     
@@ -231,23 +203,30 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  221,11,2); // eta
        ForceParticleDecay(  223,11,2); // omega
        ForceParticleDecay(  443,11,2); // J/Psi
-       ForceParticleDecay(100443,11,2); // Psi'
+       ForceParticleDecay(100443,11,2);// Psi'
        ForceParticleDecay(  553,11,2); // Upsilon
-       ForceParticleDecay(100553,11,2); // Upsilon'
-       ForceParticleDecay(200553,11,2); // Upsilon''
+       ForceParticleDecay(100553,11,2);// Upsilon'
+       ForceParticleDecay(200553,11,2);// Upsilon''
        break;
     case kBJpsiDiMuon:
-       ForceParticleDecay(  511,443,1); // B0     
-       ForceParticleDecay(  521,443,1); // B+/-     
-       ForceParticleDecay(  531,443,1); // B_s     
-       ForceParticleDecay( 5122,443,1); // Lambda_b
-       ForceParticleDecay(  443,13,2);  // J/Psi    
+
+       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,30443,1); // B0     
-       ForceParticleDecay(  521,30443,1); // B+/-     
-       ForceParticleDecay(  531,30443,1); // B_s     
-       ForceParticleDecay( 5122,30443,1); // Lambda_b 
+       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:
@@ -257,19 +236,41 @@ void AliDecayerPythia::ForceDecay()
        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,30443,1); // B0     
-       ForceParticleDecay(  521,30443,1); // B+/-     
-       ForceParticleDecay(  531,30443,1); // B_s     
-       ForceParticleDecay( 5122,30443,1); // Lambda_b 
-       ForceParticleDecay(100443,11,2);    // Psi'   
+       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 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 kHadronicD:
        ForceHadronicD();
        break;
@@ -354,7 +355,7 @@ void AliDecayerPythia::ForceHadronicD()
 void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 {
 //
-//  force decay of particle into products with multiplicity mult
+//  Force decay of particle into products with multiplicity mult
 
     Int_t kc=fPythia->Pycomp(particle);
     fPythia->SetMDCY(kc,1,1);
@@ -373,6 +374,32 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m
     }
 }
 
+void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
+{
+//
+//  Force decay of particle into products with multiplicity mult
+
+    Int_t kc=fPythia->Pycomp(particle);
+    fPythia->SetMDCY(kc,1,1);
+    Int_t ifirst=fPythia->GetMDCY(kc,2);
+    Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
+    fBraPart[kc] = 1;
+//
+//  Loop over decay channels
+    for (Int_t channel = ifirst; channel <= ilast; channel++) {
+       Int_t nprod = 0;
+       for (Int_t i = 0; i < npart; i++) {
+           nprod += (CountProducts(channel, products[i]) >= mult[i]);
+       }
+       if (nprod) {
+           fPythia->SetMDME(channel,1,1);
+       } else {
+           fPythia->SetMDME(channel,1,0);
+           fBraPart[kc]-=fPythia->GetBRAT(channel);
+       }
+    }
+}
+
 void AliDecayerPythia::DefineParticles()
 {
 //
@@ -481,11 +508,11 @@ void AliDecayerPythia::ReadDecayTable()
 {
 //
 // Read the decay table
-    printf("Reading Deacy table\n \n ");
-    
     Int_t lun = 15;
     opendecaytable(lun);
-    fPythia->Pyupda(2,lun);
+    fPythia->Pyupda(3,lun);
+    closedecaytable(lun);
+    
 }
 
 #ifdef never
@@ -509,7 +536,7 @@ void AliDecayerPythia::Streamer(TBuffer &R__b)
 }
 #endif
 
-void AliDecayerPythia::Copy(AliDecayerPythia &decayer) const
+void AliDecayerPythia::Copy(TObject &) const
 {
   //
   // Copy *this onto AliDecayerPythia -- not implemented
@@ -1489,6 +1516,7 @@ void AliDecayerPythia::Copy(AliDecayerPythia &decayer) const
            797    0   42    0.001000    e+              nu_e            pi0                                             
            798    0   42    0.001000    e+              nu_e            eta                                             
            799    0   42    0.001000    e+              nu_e            eta'                                            
+
            800    0   42    0.001000    e+              nu_e            rho0                                            
            801    0   42    0.001000    e+              nu_e            omega                                           
            802    1   42    0.070000    mu+             nu_mu           Kbar0                                           
@@ -3669,3 +3697,4 @@ void AliDecayerPythia::Copy(AliDecayerPythia &decayer) const
           2553    1    0    0.610139    W+              e-                                                              
 
 */
+