#include <TPDGCode.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
+#include <TParticle.h>
ClassImp(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.;
fPythia(0),
fDecay(kAll),
fHeavyFlavour(kTRUE),
- fLongLived(kFALSE)
+ fLongLived(kFALSE),
+ fPatchOmegaDalitz(0),
+ fPi0(1)
{
// Copy Constructor
decayer.Copy(*this);
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;
{
// 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;
}
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+/-
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+/-
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
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
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;
}
}
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)
{
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
};
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},
// 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] =
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);
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
}
}
if (norm > 0.) fBraPart[kc] /= norm;
+
+
+
}
void AliDecayerPythia::DefineParticles()
}
-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++)
{
}
-#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
{
//