]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliDecayerPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Some cleanup in the makefiles
[u/mrichter/AliRoot.git] / PYTHIA6 / AliDecayerPythia.cxx
diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx
index d4f3432f8483a92d40facd1d289461287b397b9b..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>
@@ -76,6 +77,7 @@ AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer):
 
 void AliDecayerPythia::Init()
 {
+    
 // Initialisation
 //
     if (!fgInit) {
@@ -102,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);
@@ -177,8 +182,7 @@ void AliDecayerPythia::ForceDecay()
        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(    443,  13, 2);  // J/Psi -> mu+ mu-
        ForceParticleDecay(  411,13,1); // D+/-     
        ForceParticleDecay(  421,13,1); // D0     
        ForceParticleDecay(  431,13,1); // D_s     
@@ -187,6 +191,25 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay( 4232,13,1); // Sigma_c 
        ForceParticleDecay( 4332,13,1); // Omega_c     
        break;
+   case kChiToJpsiGammaToMuonMuon:
+       products[0] =    443;
+       products[1] =     22;
+       mult[0] = 1;
+       mult[1] = 1;
+       ForceParticleDecay( 20443, products, mult, 2);   // Chi_1c  -> J/Psi  Gamma
+       ForceParticleDecay(   445, products, mult, 2);   // Chi_2c  -> J/Psi  Gamma
+       ForceParticleDecay(   443, 13, 2);               // J/Psi -> mu+ mu-
+       break;
+    case kChiToJpsiGammaToElectronElectron:
+       products[0] =    443;
+       products[1] =     22;
+       mult[0] = 1;
+       mult[1] = 1;
+       ForceParticleDecay( 20443, products, mult, 2);   // Chi_1c  -> J/Psi  Gamma
+       ForceParticleDecay(   445, products, mult, 2);   // Chi_2c  -> J/Psi  Gamma
+       ForceParticleDecay(   443, 11, 2);               // J/Psi -> e+ e-
+       break;
+
     case kBSemiMuonic:
        ForceParticleDecay(  511,13,1); // B0     
        ForceParticleDecay(  521,13,1); // B+/-     
@@ -412,6 +435,14 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
   Int_t iDeltaPP = 2224;
   Int_t productsD[2] = {kProton, kPiPlus}, multD[2] = {1, 1};
   ForceParticleDecay(iDeltaPP, productsD, multD, 2);
+  // for Lambda_c -> Lambda(1520) pi+ -> p K- pi+ 
+  Int_t iLambda_1520 = 3124;
+  Int_t productsL[2] = {kProton, kKMinus}, multL[2] = {1, 1};
+  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] = 
@@ -438,19 +469,57 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
       {-1        , -1     , -1      , -1},
       {kProton   , kKMinus, kPiPlus , 0}
     };
+  // for Lambda_c -> Lambda_1520 pi+ -> p K- pi+, D0-> K*0 pi+ pi- -> K3pi
+    Int_t decayP4[kNHadrons][4] =
+    {
+      {-1          , -1      , -1      , -1},
+      {iKstarbar0  , kPiPlus , kPiMinus,  0},
+      {-1          , -1      , -1      , -1},
+      {-1          , -1      , -1      , -1},
+      {iLambda_1520, kPiPlus ,  0      ,  0}
+    };
+  // for Lambda_c -> Lambda pi+ 
+    Int_t decayP5[kNHadrons][4] =
+    {
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {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 ((
@@ -471,6 +540,25 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
              fPythia->GetKFDP(channel,3) == decayP3[ihadron][2] &&
              fPythia->GetKFDP(channel,4) == decayP3[ihadron][3] &&
              fPythia->GetKFDP(channel,5) == 0
+        ) || (
+              fPythia->GetKFDP(channel,1) == decayP4[ihadron][0] &&
+              fPythia->GetKFDP(channel,2) == decayP4[ihadron][1] &&
+              fPythia->GetKFDP(channel,3) == decayP4[ihadron][2] &&
+              fPythia->GetKFDP(channel,4) == decayP4[ihadron][3] &&
+              fPythia->GetKFDP(channel,5) == 0
+        ) || (
+              fPythia->GetKFDP(channel,1) == decayP5[ihadron][0] &&
+              fPythia->GetKFDP(channel,2) == decayP5[ihadron][1] &&
+              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
+
              ))
          
          {
@@ -480,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
 }
 
@@ -494,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++) {
@@ -502,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)
@@ -516,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++) {
@@ -528,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()