Option kBJpsi added.

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

diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx
index a87887722ccbfcc5fdc184395e802d60136cf452..4ac8d197931d437378d7a021587115305e1361a5 100644 (file)
--- a/PYTHIA6/AliDecayerPythia.cxx
+++ b/PYTHIA6/AliDecayerPythia.cxx
@@ -13,73 +13,7 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 

-/*
-$Log$
-Revision 1.1  2003/03/15 15:00:47  morsch
-Classed imported from EVGEN.
-
-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.
 
 // Implementation of AliDecayer using Pythia
 // Method forwarding to the AliPythia instance.


@@ -99,10 +33,12 @@ ClassImp(AliDecayerPythia)
 #ifndef WIN32
 # define py1ent py1ent_
 # define opendecaytable opendecaytable_
 #ifndef WIN32
 # define py1ent py1ent_
 # define opendecaytable opendecaytable_

+# define closedecaytable closedecaytable_

 # define type_of_call
 #else
 # define lu1ent PY1ENT
 # define opendecaytable OPENDECAYTABLE
 # define type_of_call
 #else
 # define lu1ent PY1ENT
 # define opendecaytable OPENDECAYTABLE

+# define closedecaytable CLOSEDECAYTABLE

 # define type_of_call _stdcall
 #endif
 
 # define type_of_call _stdcall
 #endif
 


@@ -112,6 +48,9 @@ extern "C" void type_of_call
 extern "C" void type_of_call 
           opendecaytable(Int_t&);
 
 extern "C" void type_of_call 
           opendecaytable(Int_t&);
 

+extern "C" void type_of_call 
+          closedecaytable(Int_t&);
+

 Bool_t AliDecayerPythia::fgInit = kFALSE;
 
 
 Bool_t AliDecayerPythia::fgInit = kFALSE;
 
 


@@ -185,9 +124,39 @@ void AliDecayerPythia::ForceDecay()
 
 //
 // select mode    
 
 //
 // select mode    

-
+    Int_t products[2];
+    Int_t mult[2];
+    Int_t products1[3];
+    Int_t mult1[3];
+    printf("Forcedecay %d %d \n", decay, kAll);
+    

     switch (decay) 
     {
     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     
     case kSemiMuonic:
        ForceParticleDecay(  411,13,1); // D+/-     
        ForceParticleDecay(  421,13,1); // D0     


@@ -209,10 +178,10 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,13,2); // omega
        ForceParticleDecay(  333,13,2); // phi
        ForceParticleDecay(  443,13,2); // J/Psi
        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(  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+/-     
        break;
     case kSemiElectronic:
        ForceParticleDecay(  411,11,1); // D+/-     


@@ -235,17 +204,24 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  221,11,2); // eta
        ForceParticleDecay(  223,11,2); // omega
        ForceParticleDecay(  443,11,2); // J/Psi
        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(  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:
        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     
        break;
     case kBPsiPrimeDiMuon:
        ForceParticleDecay(  511,30443,1); // B0     


@@ -261,12 +237,18 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay( 5122,443,1); // Lambda_b
        ForceParticleDecay(  443,11,2);  // J/Psi    
        break;
        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 
     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(100443,11,2);   // Psi'   

        break;
     case kPiToMu:
        ForceParticleDecay(211,13,1); // pi->mu     
        break;
     case kPiToMu:
        ForceParticleDecay(211,13,1); // pi->mu     


@@ -358,7 +340,7 @@ void AliDecayerPythia::ForceHadronicD()
 void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 {
 //
 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);
 
     Int_t kc=fPythia->Pycomp(particle);
     fPythia->SetMDCY(kc,1,1);


@@ -377,6 +359,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()
 {
 //
 void AliDecayerPythia::DefineParticles()
 {
 //


@@ -487,7 +495,9 @@ void AliDecayerPythia::ReadDecayTable()
 // Read the decay table
     Int_t lun = 15;
     opendecaytable(lun);
 // Read the decay table
     Int_t lun = 15;
     opendecaytable(lun);

-    fPythia->Pyupda(2,lun);
+    fPythia->Pyupda(3,lun);
+    closedecaytable(lun);
+    

 }
 
 #ifdef never
 }
 
 #ifdef never


@@ -511,7 +521,7 @@ void AliDecayerPythia::Streamer(TBuffer &R__b)
 }
 #endif
 
 }
 #endif
 

-void AliDecayerPythia::Copy(AliDecayerPythia &decayer) const
+void AliDecayerPythia::Copy(TObject &) const

 {
   //
   // Copy *this onto AliDecayerPythia -- not implemented
 {
   //
   // Copy *this onto AliDecayerPythia -- not implemented