]> 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
Files consistent with the new physics selection and event cuts.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliDecayerPythia.cxx
diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx
index f0eb6495cd4acc769c002b562900ba8d4951373b..2b429d185fcd3863ce57d437f8204101e4518d10 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>
@@ -57,7 +58,8 @@ Bool_t AliDecayerPythia::fgInit = kFALSE;
 AliDecayerPythia::AliDecayerPythia():
     fPythia(AliPythia::Instance()),
     fDecay(kAll),
-    fHeavyFlavour(kTRUE)
+    fHeavyFlavour(kTRUE),
+    fLongLived(kFALSE)
 {
 // Default Constructor
     for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
@@ -68,10 +70,12 @@ AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer):
     AliDecayer(decayer),
     fPythia(0),
     fDecay(kAll),
-    fHeavyFlavour(kTRUE)
+    fHeavyFlavour(kTRUE),
+    fLongLived(kFALSE)
 {
     // Copy Constructor
     decayer.Copy(*this);
+    for (Int_t i = 0; i < 501; i++) fBraPart[i] = 0.;
 }
 
 void AliDecayerPythia::Init()
@@ -113,16 +117,18 @@ void AliDecayerPythia::Init()
 */
 
     fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
-
-    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);
-
+    Int_t isw = 0;
+    if (fLongLived) isw = 1;
+    
+    fPythia->SetMDCY(fPythia->Pycomp(310) ,1, isw);
+    fPythia->SetMDCY(fPythia->Pycomp(3122),1, isw);
+    fPythia->SetMDCY(fPythia->Pycomp(3112),1, isw);
+//    fPythia->SetMDCY(fPythia->Pycomp(3212),1, isw); // Sigma0 decays elem.
+    fPythia->SetMDCY(fPythia->Pycomp(3222),1, isw);
+    fPythia->SetMDCY(fPythia->Pycomp(3312),1, isw);
+    fPythia->SetMDCY(fPythia->Pycomp(3322),1, isw);
+    fPythia->SetMDCY(fPythia->Pycomp(3334),1, isw);
+    
 // .. Force decay channels
     ForceDecay();
 }
@@ -155,7 +161,6 @@ void AliDecayerPythia::ForceDecay()
 //
     Decay_t decay=fDecay;
     fPythia->SetMSTJ(21,2);
-    if (decay == kNoDecayHeavy) return;
 
 //
 // select mode    
@@ -281,6 +286,14 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(100553,11,2);// Upsilon'
        ForceParticleDecay(200553,11,2);// Upsilon''
        break;
+    case kPsiPrimeJpsiDiElectron:
+       products[0] =    443;
+        products[1] =    211;
+        mult[0] = 1;
+        mult[1] = 2;
+        ForceParticleDecay( 100443, products, mult, 2, 1);
+        ForceParticleDecay(    443,11,2);             
+        break; 
     case kBJpsiDiMuon:
 
        products[0] =    443;
@@ -365,6 +378,10 @@ void AliDecayerPythia::ForceDecay()
        break;
     case kOmega:
        ForceOmega();
+       break;
+    case kLambda:
+       ForceLambda();
+       break;
     case kAll:
        break;
     case kNoDecay:
@@ -373,6 +390,33 @@ void AliDecayerPythia::ForceDecay()
     case kNoDecayHeavy:
     case kNeutralPion:
        break;
+    case kNoDecayBeauty:                                                    
+        SwitchOffBDecay(); 
+       break;
+    case kElectronEM:
+        ForceParticleDecay(  111,11,1); // pi^0     
+       ForceParticleDecay(  221,11,1); // eta     
+       ForceParticleDecay(  113,11,1); // rho     
+       ForceParticleDecay(  223,11,1); // omega     
+       ForceParticleDecay(  331,11,1); // etaprime     
+       ForceParticleDecay(  333,11,1); // phi     
+       break;
+    case kDiElectronEM:
+       ForceParticleDecay(  111,11,2); // pi^0     
+       ForceParticleDecay(  221,11,2); // eta     
+       ForceParticleDecay(  113,11,2); // rho     
+       ForceParticleDecay(  223,11,2); // omega     
+       ForceParticleDecay(  331,11,2); // etaprime     
+       ForceParticleDecay(  333,11,2); // phi     
+       break;
+    case kGammaEM:
+       ForceParticleDecay(  111,22,1); // pi^0     
+       ForceParticleDecay(  221,22,1); // eta     
+       ForceParticleDecay(  113,22,1); // rho     
+       ForceParticleDecay(  223,22,1); // omega     
+       ForceParticleDecay(  331,22,1); // etaprime     
+       ForceParticleDecay(  333,22,1); // phi     
+       break;
     }
 }
 
@@ -440,6 +484,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] = 
@@ -485,10 +531,21 @@ 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; 
     }
   }  
   
@@ -497,8 +554,13 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
       Int_t kc = fPythia->Pycomp(hadron[ihadron]);     
       Int_t ifirst = fPythia->GetMDCY(kc,2);
       Int_t ilast  = ifirst + fPythia->GetMDCY(kc,3)-1;
-      Float_t norm  = 0.;
+      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);        
 
@@ -533,6 +595,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
+
              ))
          
          {
@@ -557,8 +626,13 @@ 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;
-    Float_t norm = 0.;
+    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
@@ -573,7 +647,7 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m
     if (norm > 0.) fBraPart[kc] /= norm;
 }
 
-void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
+void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart, Bool_t flag)
 {
 //
 //  Force decay of particle into products with multiplicity mult
@@ -582,8 +656,13 @@ 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;
-    Float_t norm = 0.;
+    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
@@ -592,9 +671,9 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t
        for (Int_t i = 0; i < npart; i++) {
            nprod += (CountProducts(channel, products[i]) >= mult[i]);
        }
-       if (nprod) {
+       if ((nprod && !flag) || ((nprod == npart) && flag)) {
            fPythia->SetMDME(channel,1,1);
-       } else {
+       } else {                // 
            fPythia->SetMDME(channel,1,0);
            fBraPart[kc] -= fPythia->GetBRAT(channel);
        }
@@ -683,7 +762,35 @@ void  AliDecayerPythia::ForceOmega()
     } // decay channels
 }
 
+void  AliDecayerPythia::ForceLambda()
+{
+    // Force Lambda -> p pi-
+    Int_t kc=fPythia->Pycomp(3122);
+    fPythia->SetMDCY(kc,1,1);  
+    Int_t ifirst = fPythia->GetMDCY(kc,2);
+    Int_t ilast  = ifirst + fPythia->GetMDCY(kc,3)-1;
+    for (Int_t channel = ifirst; channel <= ilast; channel++) {
+       if (
+           fPythia->GetKFDP(channel,1) == kProton   &&
+           fPythia->GetKFDP(channel,2) == kPiMinus  &&
+           fPythia->GetKFDP(channel,3) == 0
+           )
+       {
+           fPythia->SetMDME(channel,1,1);
+       } else {
+           fPythia->SetMDME(channel,1,0);
+       } // selected channel ?
+    } // decay channels
+}
+
 
+void AliDecayerPythia::SwitchOffBDecay(){
+  Int_t heavyB[]={511,521,531,5122,5132,5232,5332};
+  for(int i=0;i<4;i++)
+    {
+      fPythia->SetMDCY(fPythia->Pycomp(heavyB[i]),1,0);
+    }
+}
 
 Float_t  AliDecayerPythia::GetPartialBranchingRatio(Int_t kf)
 {