Fix for the loophole in the magnets currents check: the L3Off/DipON was passing the...

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

diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx
index dd8f190a54220ee020233e7b3050b9b912fdfb67..affa1648ace9cdd94c71c17156fbfae6ce467e71 100644 (file)
--- a/PYTHIA6/AliDecayerPythia.cxx
+++ b/PYTHIA6/AliDecayerPythia.cxx
@@ -24,6 +24,7 @@
 
 #include "AliDecayerPythia.h"
 #include "AliPythia.h"
+#include "AliLog.h"
 #include <TPDGCode.h>
 #include <TLorentzVector.h>
 #include <TClonesArray.h>
@@ -103,13 +104,16 @@ void AliDecayerPythia::Init()
        }
     }
 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
-    
+
+/*
     if (fDecay != kNeutralPion) {
        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
     } else {
        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
     }
+*/
 
+    fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
 
     fPythia->SetMDCY(fPythia->Pycomp(310) ,1,0);
     fPythia->SetMDCY(fPythia->Pycomp(3122),1,0);
@@ -437,6 +441,8 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
   ForceParticleDecay(iLambda_1520, productsL, multL, 2);
   // for Lambda_c -> Lambda pi+
   Int_t iLambda=3122;
+  //for Lambda_c -> antiK0 p
+  Int_t iK0bar=-311;
 
 
   Int_t decayP1[kNHadrons][4] = 
@@ -482,19 +488,38 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
       {iLambda   , kPiPlus,  0      ,  0}
     };
 
+    // for Lambda_c -> K0bar p
+    Int_t decayP6[kNHadrons][4] =
+    {
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {kProton    , iK0bar,  0      ,  0}
+    };
+
   if(optUse4Bodies==0){
     for(Int_t iDau=0;iDau<4;iDau++){
       decayP2[1][iDau]=-1;
       decayP3[1][iDau]=-1;
+      decayP4[1][iDau]=-1; 
     }
   }  
   
   for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
     {
       Int_t kc = fPythia->Pycomp(hadron[ihadron]);     
-      fPythia->SetMDCY(kc,1,1);        
       Int_t ifirst = fPythia->GetMDCY(kc,2);
       Int_t ilast  = ifirst + fPythia->GetMDCY(kc,3)-1;
+      Double_t norm  = 0.;
+      for (channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+      if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
+        char pName[16];
+        fPythia->Pyname(hadron[ihadron],pName);
+        AliWarning(Form("Total branching ratio of %s (PDG code = %d) not equal to 1 (= %f)",pName,hadron[ihadron],norm));
+      }
+      fBraPart[kc] = norm;
+      fPythia->SetMDCY(kc,1,1);        
 
       for (channel = ifirst; channel <= ilast; channel++) {
        if ((
@@ -527,6 +552,13 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
               fPythia->GetKFDP(channel,3) == decayP5[ihadron][2] &&
               fPythia->GetKFDP(channel,4) == decayP5[ihadron][3] &&
               fPythia->GetKFDP(channel,5) == 0
+        ) || (
+              fPythia->GetKFDP(channel,1) == decayP6[ihadron][0] &&
+              fPythia->GetKFDP(channel,2) == decayP6[ihadron][1] &&
+              fPythia->GetKFDP(channel,3) == decayP6[ihadron][2] &&
+              fPythia->GetKFDP(channel,4) == decayP6[ihadron][3] &&
+              fPythia->GetKFDP(channel,5) == 0
+
              ))
          
          {
@@ -536,6 +568,7 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
            fBraPart[kc] -= fPythia->GetBRAT(channel);
          } // selected channel ?
       } // decay channels
+      if (norm > 0) fBraPart[kc] /= norm;
     } // hadrons
 }
 
@@ -550,7 +583,14 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m
     fPythia->SetMDCY(kc,1,1);
     Int_t ifirst=fPythia->GetMDCY(kc,2);
     Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
-    fBraPart[kc] = 1;
+    Double_t norm = 0.;
+    for (Int_t channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+    if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
+        char pName[16];
+        fPythia->Pyname(particle,pName);
+        AliWarning(Form("Total branching ratio of %s (PDG code = %d) not equal to 1 (= %f)",pName,particle,norm));
+    }
+    fBraPart[kc] = norm;
 //
 //  Loop over decay channels
     for (Int_t channel=ifirst; channel<=ilast;channel++) {
@@ -558,9 +598,10 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m
            fPythia->SetMDME(channel,1,1);
        } else {
            fPythia->SetMDME(channel,1,0);
-           fBraPart[kc]-=fPythia->GetBRAT(channel);
+           fBraPart[kc] -= fPythia->GetBRAT(channel);
        }
     }
+    if (norm > 0.) fBraPart[kc] /= norm;
 }
 
 void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
@@ -572,7 +613,14 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t
     fPythia->SetMDCY(kc,1,1);
     Int_t ifirst=fPythia->GetMDCY(kc,2);
     Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
-    fBraPart[kc] = 1;
+    Double_t norm = 0.;
+    for (Int_t channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+    if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
+        char pName[16];
+        fPythia->Pyname(particle,pName);
+        AliWarning(Form("Total branching ratio of %s (PDG code = %d) not equal to 1 (= %f)",pName,particle,norm));
+    }
+    fBraPart[kc] = norm;
 //
 //  Loop over decay channels
     for (Int_t channel = ifirst; channel <= ilast; channel++) {
@@ -584,9 +632,10 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t
            fPythia->SetMDME(channel,1,1);
        } else {
            fPythia->SetMDME(channel,1,0);
-           fBraPart[kc]-=fPythia->GetBRAT(channel);
+           fBraPart[kc] -= fPythia->GetBRAT(channel);
        }
     }
+    if (norm > 0.) fBraPart[kc] /= norm;
 }
 
 void AliDecayerPythia::DefineParticles()