X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PYTHIA6%2FAliDecayerPythia.cxx;h=1465fdf9865ed16cd5b5f016fe828eec06a0a412;hb=1e2a45d476f5ebdc72da16cf15591845b378f51a;hp=a433adc5ecd658f00abbe0ffe6bbabb96a179062;hpb=b64fad39394f6631d38e78dfb2d1a4b66bac4c4a;p=u%2Fmrichter%2FAliRoot.git diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx index a433adc5ecd..1465fdf9865 100644 --- a/PYTHIA6/AliDecayerPythia.cxx +++ b/PYTHIA6/AliDecayerPythia.cxx @@ -24,9 +24,11 @@ #include "AliDecayerPythia.h" #include "AliPythia.h" +#include "AliLog.h" #include #include #include +#include ClassImp(AliDecayerPythia) @@ -57,7 +59,9 @@ Bool_t AliDecayerPythia::fgInit = kFALSE; AliDecayerPythia::AliDecayerPythia(): fPythia(AliPythia::Instance()), fDecay(kAll), - fHeavyFlavour(kTRUE) + fHeavyFlavour(kTRUE), + fLongLived(kFALSE), + fPatchOmegaDalitz(0) { // Default Constructor for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.; @@ -68,10 +72,13 @@ AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer): AliDecayer(decayer), fPythia(0), fDecay(kAll), - fHeavyFlavour(kTRUE) + fHeavyFlavour(kTRUE), + fLongLived(kFALSE), + fPatchOmegaDalitz(0) { // Copy Constructor decayer.Copy(*this); + for (Int_t i = 0; i < 501; i++) fBraPart[i] = 0.; } void AliDecayerPythia::Init() @@ -113,16 +120,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(); } @@ -134,8 +143,15 @@ void AliDecayerPythia::Decay(Int_t idpart, TLorentzVector* p) Float_t energy = p->Energy(); Float_t theta = p->Theta(); Float_t phi = p->Phi(); - - Lu1Ent(0, idpart, energy, theta, phi); + if (!fPatchOmegaDalitz) { + Lu1Ent(0, idpart, energy, theta, phi); + } else { + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0); + Lu1Ent(0, idpart, energy, theta, phi); + fPythia->DalitzDecays(); + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1); + fPythia->Pyexec(); + } fPythia->GetPrimaries(); } @@ -143,7 +159,16 @@ 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;ipAt(ip); + if (!prod) continue; + prod->SetProductionVertex(prod->Vx()*kconvL,prod->Vy()*kconvL,prod->Vz()*kconvL,kconvT*prod->T()); + } + return np; } @@ -155,7 +180,6 @@ void AliDecayerPythia::ForceDecay() // Decay_t decay=fDecay; fPythia->SetMSTJ(21,2); - if (decay == kNoDecayHeavy) return; // // select mode @@ -245,6 +269,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+/- @@ -281,6 +308,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; @@ -315,6 +350,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+/- @@ -365,6 +407,10 @@ void AliDecayerPythia::ForceDecay() break; case kOmega: ForceOmega(); + break; + case kLambda: + ForceLambda(); + break; case kAll: break; case kNoDecay: @@ -373,6 +419,36 @@ 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; + case kBeautyUpgrade: + ForceBeautyUpgrade(); + break; } } @@ -388,6 +464,21 @@ 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) + // + + ForceParticleDecay( 5122, 4122, 1); + ForceParticleDecay( 4232, 3312, 1); + const Int_t prod[2]={421,211}; + Int_t mult[2]={1,1}; + ForceParticleDecay(521,prod,mult,2,1); + ForceHadronicD(1); +} + void AliDecayerPythia::Lu1Ent(Int_t flag, Int_t idpart, Double_t mom, Double_t theta, Double_t phi) { @@ -435,11 +526,13 @@ 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 + Int_t iK0bar=-311; Int_t decayP1[kNHadrons][4] = @@ -460,7 +553,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}, @@ -469,11 +562,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] = @@ -485,19 +578,38 @@ 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; } } 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 (( @@ -530,6 +642,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 + )) { @@ -539,6 +658,7 @@ void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies) fBraPart[kc] -= fPythia->GetBRAT(channel); } // selected channel ? } // decay channels + if (norm > 0) fBraPart[kc] /= norm; } // hadrons } @@ -553,7 +673,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++) { @@ -561,12 +688,13 @@ 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) +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 @@ -575,7 +703,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++) { @@ -583,13 +718,14 @@ 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); + fBraPart[kc] -= fPythia->GetBRAT(channel); } } + if (norm > 0.) fBraPart[kc] /= norm; } void AliDecayerPythia::DefineParticles() @@ -673,7 +809,37 @@ 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() +{ +// Switch off B-decays + 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) { @@ -707,27 +873,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 { //