Switch off decays of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-

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

diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx
index 89dfddd53209bec94a95ae2e25bec99606045e04..f0a98cddddd6f2b76771644bad985fba0b80d3c2 100644 (file)
--- a/PYTHIA6/AliDecayerPythia.cxx
+++ b/PYTHIA6/AliDecayerPythia.cxx
@@ -89,7 +89,17 @@ void AliDecayerPythia::Init()
          }
        }
     }
-
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+    fPythia->SetMDCY(fPythia->Pycomp(111) ,1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(310) ,1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3122),1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3112),1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3212),1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3222),1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3312),1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3322),1,0);
+    fPythia->SetMDCY(fPythia->Pycomp(3334),1,0);
+// .. Force decay channels
     ForceDecay();
 }
 
@@ -256,7 +266,7 @@ void AliDecayerPythia::ForceDecay()
        break;
     case kKaToMu:
        ForceParticleDecay(321,13,1); // K->mu     
-       break;     
+       break;
     case kWToMuon:
         ForceParticleDecay(  24, 13,1); // W -> mu
        break;
@@ -264,15 +274,18 @@ void AliDecayerPythia::ForceDecay()
         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;
+        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 kHadronicD:
        ForceHadronicD();
        break;
@@ -288,9 +301,6 @@ void AliDecayerPythia::ForceDecay()
        break;
     case kNoDecayHeavy:
        break;
-    case kNewTest:
-       ForceParticleDecay(  511,13,1); // B0 -> mu
-       break;
     }
 }
 
@@ -324,12 +334,28 @@ void AliDecayerPythia::ForceHadronicD()
     const Int_t kNHadrons = 4;
     Int_t channel;
     Int_t hadron[kNHadrons] = {411,  421, 431, 4112};
-    Int_t decayP[kNHadrons][3] = 
+    // 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-
+    Int_t decayP1[kNHadrons][3] = 
+
+       { 
+           {kKMinus, kPiPlus,    kPiPlus},
+           {kKMinus, kPiPlus,    0      },
+           {kKPlus , iKstarbar0, 0     },
+           {-1     , -1        , -1        }
+       };
+    Int_t decayP2[kNHadrons][3] = 
        { 
-           {kKMinus, kPiPlus, kPiPlus},
-           {kKMinus, kPiPlus, 0      },
-           {-1     , -1     , -1     },
-           {-1     , -1     , -1     }
+           {iKstarbar0, kPiPlus, 0   },
+           {-1        , -1     , -1  },
+           {iPhi      , kPiPlus, 0  },
+           {-1        , -1     , -1  }
        };
     
 
@@ -341,12 +367,18 @@ void AliDecayerPythia::ForceHadronicD()
        Int_t ilast  = ifirst + fPythia->GetMDCY(kc,3)-1;
 
        for (channel = ifirst; channel <= ilast; channel++) {
-           if (
-               fPythia->GetKFDP(channel,1) == decayP[ihadron][0] &&
-               fPythia->GetKFDP(channel,2) == decayP[ihadron][1] &&
-               fPythia->GetKFDP(channel,3) == decayP[ihadron][2] &&
+           if ((
+               fPythia->GetKFDP(channel,1) == decayP1[ihadron][0] &&
+               fPythia->GetKFDP(channel,2) == decayP1[ihadron][1] &&
+               fPythia->GetKFDP(channel,3) == decayP1[ihadron][2] &&
                fPythia->GetKFDP(channel,4) == 0
-               )
+               ) || (
+               fPythia->GetKFDP(channel,1) == decayP2[ihadron][0] &&
+               fPythia->GetKFDP(channel,2) == decayP2[ihadron][1] &&
+               fPythia->GetKFDP(channel,3) == decayP2[ihadron][2] &&
+               fPythia->GetKFDP(channel,4) == 0
+               ))
+
            {
                fPythia->SetMDME(channel,1,1);
            } else {
@@ -357,6 +389,7 @@ void AliDecayerPythia::ForceHadronicD()
     } // hadrons
 }
 
+
 void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 {
 //