]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliDecayerPythia.cxx
fix spline selection for MC for train tests - Extension
[u/mrichter/AliRoot.git] / PYTHIA6 / AliDecayerPythia.cxx
index 2b429d185fcd3863ce57d437f8204101e4518d10..0b9be557115e33c6e218160b46d2925e8a879f65 100644 (file)
@@ -28,6 +28,7 @@
 #include <TPDGCode.h>
 #include <TLorentzVector.h>
 #include <TClonesArray.h>
+#include <TParticle.h>
 
 ClassImp(AliDecayerPythia)
 
@@ -59,7 +60,9 @@ AliDecayerPythia::AliDecayerPythia():
     fPythia(AliPythia::Instance()),
     fDecay(kAll),
     fHeavyFlavour(kTRUE),
-    fLongLived(kFALSE)
+    fLongLived(kFALSE),
+    fPatchOmegaDalitz(0),
+    fPi0(1)
 {
 // Default Constructor
     for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
@@ -71,7 +74,9 @@ AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer):
     fPythia(0),
     fDecay(kAll),
     fHeavyFlavour(kTRUE),
-    fLongLived(kFALSE)
+    fLongLived(kFALSE),
+    fPatchOmegaDalitz(0),
+    fPi0(1)
 {
     // Copy Constructor
     decayer.Copy(*this);
@@ -115,8 +120,8 @@ void AliDecayerPythia::Init()
        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
     }
 */
+    if (fPi0) fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
 
-    fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
     Int_t isw = 0;
     if (fLongLived) isw = 1;
     
@@ -137,19 +142,74 @@ void AliDecayerPythia::Decay(Int_t idpart, TLorentzVector* p)
 {
 //  Decay a particle
 // 
-    Float_t energy = p->Energy();
-    Float_t theta  = p->Theta();
-    Float_t phi    = p->Phi();
-    
+  Float_t energy = p->Energy();
+  Float_t theta  = p->Theta();
+  Float_t phi    = p->Phi();
+
+  if(!fDecayerExodus) {
     Lu1Ent(0, idpart, energy, theta, phi);
-    fPythia->GetPrimaries();
+  } else {
+
+  // EXODUS decayer
+    if(idpart == 111){
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->PizeroDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+    }
+    else if(idpart == 221){
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->EtaDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+    }
+    else if(idpart == 113){
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->RhoDirect();
+    }
+    else if(idpart == 223){
+      fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->OmegaDirect();
+      fPythia->OmegaDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+    }
+    else if(idpart == 331){
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->EtaprimeDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+    }
+    else if(idpart == 333){
+      fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->PhiDirect();
+      fPythia->PhiDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+    }
+    else if(idpart == 443){
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->JPsiDirect();
+    }  
+  }
+
+  fPythia->GetPrimaries();
 }
 
 Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
 {
 // Import the decay products
 //
-    return fPythia->ImportParticles(particles, "All");
+  const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
+  const Float_t kconvL=1./10; // mm to cm conversion
+  int np = fPythia->ImportParticles(particles, "All");
+  // pythia assigns decay time in mm/c, convert to seconds
+  for (int ip=1;ip<np;ip++) {
+    TParticle* prod = (TParticle*)particles->At(ip);
+    if (!prod) continue;
+    prod->SetProductionVertex(prod->Vx()*kconvL,prod->Vy()*kconvL,prod->Vz()*kconvL,kconvT*prod->T());
+  }
+  return np;
 }
 
 
@@ -250,6 +310,9 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(100553,13,2);// Upsilon'
        ForceParticleDecay(200553,13,2);// Upsilon''
        break;
+    case kJpsiDiMuon:
+       ForceParticleDecay(  443,13,2); // J/Psi
+       break;
     case kBSemiElectronic:
        ForceParticleDecay(  511,11,1); // B0     
        ForceParticleDecay(  521,11,1); // B+/-     
@@ -328,6 +391,13 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  531,443,1); // B_s     
        ForceParticleDecay( 5122,443,1); // Lambda_b
        break;
+    case kBJpsiUndecayed:
+        ForceParticleDecay(  511,443,1); // B0     
+        ForceParticleDecay(  521,443,1); // B+/-     
+        ForceParticleDecay(  531,443,1); // B_s     
+        ForceParticleDecay( 5122,443,1); // Lambda_b
+        fPythia->SetMDCY(fPythia->Pycomp(443),1,0); // switch-off J/psi 
+        break;
     case kBPsiPrimeDiElectron:
        ForceParticleDecay(  511,100443,1); // B0     
        ForceParticleDecay(  521,100443,1); // B+/-     
@@ -400,6 +470,7 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,11,1); // omega     
        ForceParticleDecay(  331,11,1); // etaprime     
        ForceParticleDecay(  333,11,1); // phi     
+       ForceParticleDecay(  443,11,1); // jpsi     
        break;
     case kDiElectronEM:
        ForceParticleDecay(  111,11,2); // pi^0     
@@ -408,6 +479,7 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,11,2); // omega     
        ForceParticleDecay(  331,11,2); // etaprime     
        ForceParticleDecay(  333,11,2); // phi     
+       ForceParticleDecay(  443,11,2); // jpsi     
        break;
     case kGammaEM:
        ForceParticleDecay(  111,22,1); // pi^0     
@@ -416,7 +488,11 @@ void AliDecayerPythia::ForceDecay()
        ForceParticleDecay(  223,22,1); // omega     
        ForceParticleDecay(  331,22,1); // etaprime     
        ForceParticleDecay(  333,22,1); // phi     
+       ForceParticleDecay(  443,22,1); // jpsi     
        break;
+    case kBeautyUpgrade:
+        ForceBeautyUpgrade();
+      break;
     }
 }
 
@@ -432,6 +508,34 @@ void  AliDecayerPythia::SwitchOffHeavyFlavour()
     for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
 }
 
+void  AliDecayerPythia::ForceBeautyUpgrade()
+{
+ //
+ // Force dedicated decay channels of signals ineresting
+ // for the ITS upgrade (Lb, Lc, Xi_c, B)
+ //
+
+   // Lb: 50% of them in Lc pi+ and the rest with a Lc in final state
+   if(gRandom->Rndm()<0.50) {
+     const Int_t prod3[3]={4122,211,0};
+     Int_t mult3[3]={1,1,1};
+     ForceParticleDecay(5122,prod3,mult3,3,1);
+   } else {
+     ForceParticleDecay( 5122, 4122, 1);
+   }
+   // Xi_c
+   ForceParticleDecay( 4232, 3312, 1);
+   // B+ -> D0pi+
+   const Int_t prod[2]={421,211};
+   Int_t mult[2]={1,1};
+   ForceParticleDecay(521,prod,mult,2,1);
+   // B0 -> D*-pi+
+   const Int_t prod2[2]={413,211};
+   ForceParticleDecay(511,prod2,mult,2,1);
+   // force charm hadronic decays (D mesons and Lc)
+   ForceHadronicD(0);
+}
+
 void  AliDecayerPythia::Lu1Ent(Int_t flag, Int_t idpart, 
                      Double_t mom, Double_t theta, Double_t phi)
 {
@@ -479,9 +583,9 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
   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 iLambda1520 = 3124;
   Int_t productsL[2] = {kProton, kKMinus}, multL[2] = {1, 1};
-  ForceParticleDecay(iLambda_1520, productsL, multL, 2);
+  ForceParticleDecay(iLambda1520, productsL, multL, 2);
   // for Lambda_c -> Lambda pi+
   Int_t iLambda=3122;
   //for Lambda_c -> antiK0 p
@@ -506,7 +610,7 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
     };
   Int_t decayP3[kNHadrons][4] = 
     { 
-      {-1        , -1     , -1      , -1},
+      {kPiPlus   , iPhi   , 0       , 0},
       {kKMinus   , kPiPlus, iRho0   , 0 },
       {-1        , -1     , -1      , -1},
       {-1        , -1     , -1      , -1},
@@ -515,11 +619,11 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
   // 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  , kKPlus  , 0       ,  0},
       {iKstarbar0  , kPiPlus , kPiMinus,  0},
       {-1          , -1      , -1      , -1},
       {-1          , -1      , -1      , -1},
-      {iLambda_1520, kPiPlus ,  0      ,  0}
+      {iLambda1520 , kPiPlus ,  0      ,  0}
     };
   // for Lambda_c -> Lambda pi+ 
     Int_t decayP5[kNHadrons][4] =
@@ -627,7 +731,9 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m
     Int_t ifirst=fPythia->GetMDCY(kc,2);
     Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
     Double_t norm = 0.;
-    for (Int_t channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+    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);
@@ -647,7 +753,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, Bool_t flag)
+void AliDecayerPythia::ForceParticleDecay(Int_t particle, const Int_t* products, Int_t* mult, Int_t npart, Bool_t flag)
 {
 //
 //  Force decay of particle into products with multiplicity mult
@@ -679,6 +785,9 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t
        }
     }
     if (norm > 0.) fBraPart[kc] /= norm;
+
+
+
 }
 
 void AliDecayerPythia::DefineParticles()
@@ -784,7 +893,9 @@ void  AliDecayerPythia::ForceLambda()
 }
 
 
-void AliDecayerPythia::SwitchOffBDecay(){
+void AliDecayerPythia::SwitchOffBDecay()
+{
+// Switch off B-decays
   Int_t heavyB[]={511,521,531,5122,5132,5232,5332};
   for(int i=0;i<4;i++)
     {
@@ -824,27 +935,6 @@ void AliDecayerPythia::ReadDecayTable()
     
 }
 
-#ifdef never
-void AliDecayerPythia::Streamer(TBuffer &R__b)
-{
-   // Stream an object of class AliDecayerPythia.
-
-   if (R__b.IsReading()) {
-      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
-      AliDecayer::Streamer(R__b);
-      (AliPythia::Instance())->Streamer(R__b);
-      R__b >> (Int_t&)fDecay;
-      R__b.ReadStaticArray(fBraPart);
-   } else {
-      R__b.WriteVersion(AliDecayerPythia::IsA());
-      AliDecayer::Streamer(R__b);
-      R__b << fPythia;
-      R__b << (Int_t)fDecay;
-      R__b.WriteArray(fBraPart, 501);
-   }
-}
-#endif
-
 void AliDecayerPythia::Copy(TObject &) const
 {
   //