]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliDecayerPythia.cxx
All the necessary for fast Chi_c generation. (P. Ladron de Guevara)
[u/mrichter/AliRoot.git] / PYTHIA6 / AliDecayerPythia.cxx
index 6712c618ede14480462401554d09b2aa06243bd1..7579cfecfecc7da0816197d96084c79398cf3e55 100644 (file)
@@ -56,7 +56,8 @@ Bool_t AliDecayerPythia::fgInit = kFALSE;
 
 AliDecayerPythia::AliDecayerPythia():
     fPythia(AliPythia::Instance()),
-    fDecay(kAll)
+    fDecay(kAll),
+    fHeavyFlavour(kTRUE)
 {
 // Default Constructor
     for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
@@ -66,7 +67,8 @@ AliDecayerPythia::AliDecayerPythia():
 AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer):
     AliDecayer(decayer),
     fPythia(0),
-    fDecay(kAll)
+    fDecay(kAll),
+    fHeavyFlavour(kTRUE)
 {
     // Copy Constructor
     decayer.Copy(*this);
@@ -100,14 +102,14 @@ void AliDecayerPythia::Init()
        }
     }
 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
-
+    
     if (fDecay != kNeutralPion) {
-       fPythia->SetMDCY(fPythia->Pycomp(111) ,1,0);
+       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);
     fPythia->SetMDCY(fPythia->Pycomp(3112),1,0);
@@ -116,6 +118,7 @@ void AliDecayerPythia::Init()
     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();
 }
@@ -130,8 +133,6 @@ void AliDecayerPythia::Decay(Int_t idpart, TLorentzVector* p)
     
     Lu1Ent(0, idpart, energy, theta, phi);
     fPythia->GetPrimaries();
-//    fPythia->Pylist(1);
-    
 }
 
 Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
@@ -145,6 +146,9 @@ Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
 void AliDecayerPythia::ForceDecay()
 {
 // Force a particle decay mode
+// Switch heavy flavour production off if requested
+    if (!fHeavyFlavour) SwitchOffHeavyFlavour();
+//
     Decay_t decay=fDecay;
     fPythia->SetMSTJ(21,2);
     if (decay == kNoDecayHeavy) return;
@@ -173,8 +177,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     
@@ -183,6 +186,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+/-     
@@ -302,6 +324,10 @@ void AliDecayerPythia::ForceDecay()
     case kKaToMu:
        ForceParticleDecay(321,13,1); // K->mu     
        break;
+    case kAllMuonic:
+       ForceParticleDecay(211,13,1); // pi->mu
+       ForceParticleDecay(321,13,1); // K->mu    
+       break;
     case kWToMuon:
         ForceParticleDecay(  24, 13,1); // W -> mu
        break;
@@ -321,8 +347,14 @@ void AliDecayerPythia::ForceDecay()
     case kZDiMuon:
         ForceParticleDecay(  23, 13,2); // Z -> mu+ mu-
        break;
+    case kZDiElectron:
+        ForceParticleDecay(  23, 11,2); // Z -> e+ e-
+       break;
     case kHadronicD:
-       ForceHadronicD();
+       ForceHadronicD(1);
+       break;
+    case kHadronicDWithout4Bodies:
+       ForceHadronicD(0);
        break;
     case kPhiKK:
        ForceParticleDecay(333,321,2); // Phi->K+K-
@@ -337,10 +369,21 @@ void AliDecayerPythia::ForceDecay()
     case kNoDecayHeavy:
     case kNeutralPion:
        break;
-
     }
 }
 
+void  AliDecayerPythia::SwitchOffHeavyFlavour()
+{
+    // Switch off heavy flavour production
+    //
+    // Maximum number of quark flavours used in pdf 
+    fPythia->SetMSTP(58, 3);
+    // Maximum number of flavors that can be used in showers
+    fPythia->SetMSTJ(45, 3);   
+    // Switch off g->QQbar splitting in decay table
+    for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
+}
+
 void  AliDecayerPythia::Lu1Ent(Int_t flag, Int_t idpart, 
                      Double_t mom, Double_t theta, Double_t phi)
 {
@@ -363,70 +406,103 @@ Int_t AliDecayerPythia::CountProducts(Int_t channel, Int_t particle)
 }
 
 
-void AliDecayerPythia::ForceHadronicD()
+void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
 {
 //
 // Force golden D decay modes
 //
-    const Int_t kNHadrons = 4;
-    Int_t channel;
-    Int_t hadron[kNHadrons] = {411,  421, 431, 4112};
-    // 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] = 
-       { 
-           {iKstarbar0, kPiPlus, 0   },
-           {-1        , -1     , -1  },
-           {iPhi      , kPiPlus, 0  },
-           {-1        , -1     , -1  }
-       };
-    
-
-    for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
+  const Int_t kNHadrons = 5;
+  Int_t channel;
+  Int_t hadron[kNHadrons] = {411,  421, 431, 4112, 4122};
+
+  // 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-
+  // for D0 -> rho0 pi+ k-
+  Int_t iRho0=113;
+  ForceParticleDecay(iRho0,kPiPlus,2); // Rho0->pi+pi-
+  // for Lambda_c -> Delta++ K-
+  Int_t iDeltaPP = 2224;
+  Int_t productsD[2] = {kProton, kPiPlus}, multD[2] = {1, 1};
+  ForceParticleDecay(iDeltaPP, productsD, multD, 2);
+
+
+  Int_t decayP1[kNHadrons][4] = 
+    { 
+      {kKMinus, kPiPlus,    kPiPlus, 0},
+      {kKMinus, kPiPlus,    0      , 0},
+      {kKPlus , iKstarbar0, 0      , 0},
+      {-1     , -1        , -1     , -1},
+      {kProton, iKstarbar0, 0      , 0}
+    };
+  Int_t decayP2[kNHadrons][4] = 
+    { 
+      {iKstarbar0, kPiPlus, 0      , 0},
+      {kKMinus   , kPiPlus, kPiPlus, kPiMinus},
+      {iPhi      , kPiPlus, 0      , 0},
+      {-1        , -1     , -1     , -1},
+      {iDeltaPP  , kKMinus, 0      , 0}
+    };
+  Int_t decayP3[kNHadrons][4] = 
+    { 
+      {-1        , -1     , -1      , -1},
+      {kKMinus   , kPiPlus, iRho0   , 0 },
+      {-1        , -1     , -1      , -1},
+      {-1        , -1     , -1      , -1},
+      {kProton   , kKMinus, kPiPlus , 0}
+    };
+  if(optUse4Bodies==0){
+    for(Int_t iDau=0;iDau<4;iDau++){
+      decayP2[1][iDau]=-1;
+      decayP3[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;
-
-       for (channel = ifirst; channel <= ilast; channel++) {
-           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 {
-               fPythia->SetMDME(channel,1,0);
-               fBraPart[kc] -= fPythia->GetBRAT(channel);
-           } // selected channel ?
-       } // decay channels
+      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;
+
+      for (channel = ifirst; channel <= ilast; channel++) {
+       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) == decayP1[ihadron][3] &&
+            fPythia->GetKFDP(channel,5) == 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) == decayP2[ihadron][3] &&
+            fPythia->GetKFDP(channel,5) == 0
+       ) || (
+             fPythia->GetKFDP(channel,1) == decayP3[ihadron][0] &&
+             fPythia->GetKFDP(channel,2) == decayP3[ihadron][1] &&
+             fPythia->GetKFDP(channel,3) == decayP3[ihadron][2] &&
+             fPythia->GetKFDP(channel,4) == decayP3[ihadron][3] &&
+             fPythia->GetKFDP(channel,5) == 0
+             ))
+         
+         {
+           fPythia->SetMDME(channel,1,1);
+         } else {
+           fPythia->SetMDME(channel,1,0);
+           fBraPart[kc] -= fPythia->GetBRAT(channel);
+         } // selected channel ?
+      } // decay channels
     } // hadrons
 }
 
 
+
 void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 {
 //