Updates in SDD preprocessor: new infos from DA, HLT mode no longer needed (F. Prino)

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

diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx
index a433adc5ecd658f00abbe0ffe6bbabb96a179062..125956ee2cd5e50157ca8e9e948010791e0ea330 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>
@@ -489,15 +490,24 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
     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 ((
@@ -539,6 +549,7 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
            fBraPart[kc] -= fPythia->GetBRAT(channel);
          } // selected channel ?
       } // decay channels
+      if (norm > 0) fBraPart[kc] /= norm;
     } // hadrons
 }
 
@@ -553,7 +564,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++) {
@@ -561,9 +579,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)
@@ -575,7 +594,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++) {
@@ -587,9 +613,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()