//
AliOmegaDalitz::AliOmegaDalitz():
AliDecayer(),
- fLPMass(0)
+ fEPMass(0),
+ fMPMass(0),
+ fInit(0)
{
// Constructor
}
void AliOmegaDalitz::Init()
{
// Initialisation
- Int_t ibin, nbins = 1000;
+ Int_t idecay, ibin, nbins = 1000;
Double_t min, max, binwidth;
- Double_t pmass, lmass, omass;
+ Double_t pmass, lmass, omass, emass, mmass;
Double_t epsilon, delta, mLL, q, kwHelp, krollWada, formFactor, weight;
// Get the particle masses
// electron
- lmass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
+ emass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
+ // muon
+ mmass = (TDatabasePDG::Instance()->GetParticle(kMuonPlus))->Mass();
// omega
pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
// pi0
omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
- min = 2.0 * lmass;
- max = pmass - omass;
- binwidth = (max - min) / (Double_t)nbins;
- fLPMass = new TH1F("fLPMass", "Dalitz", nbins, min, max);
-
- epsilon = (lmass / pmass) * (lmass / pmass);
- delta = (omass / pmass) * (omass / pmass);
-
- for ( ibin = 1; ibin <= nbins; ibin++ )
+ for ( idecay = 1; idecay<=2; idecay++ )
{
- mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
- q = (mLL / pmass) * (mLL / pmass);
- if ( q <= 4.0 * epsilon )
+ if ( idecay == 1 )
+ lmass = emass;
+ else
+ lmass = mmass;
+
+ min = 2.0 * lmass;
+ max = pmass - omass;
+ binwidth = (max - min) / (Double_t)nbins;
+ if ( idecay == 1 )
+ fEPMass = new TH1F("fEPMass", "Dalitz electron pair", nbins, min, max);
+ else
+ fMPMass = new TH1F("fMPMass", "Dalitz muon pair", nbins, min, max);
+
+ epsilon = (lmass / pmass) * (lmass / pmass);
+ delta = (omass / pmass) * (omass / pmass);
+
+ for ( ibin = 1; ibin <= nbins; ibin++ )
{
+ mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
+ q = (mLL / pmass) * (mLL / pmass);
+ if ( q <= 4.0 * epsilon )
+ {
AliFatal("Error in calculating Dalitz mass histogram binning!");
- fLPMass = 0;
- }
+ }
- kwHelp = (1.0 + q / (1.0 - delta)) * (1.0 + q / (1.0 - delta))
- - 4.0 * q / ((1.0 - delta) * (1.0 - delta));
- if ( kwHelp <= 0.0 )
- {
+ kwHelp = (1.0 + q / (1.0 - delta)) * (1.0 + q / (1.0 - delta))
+ - 4.0 * q / ((1.0 - delta) * (1.0 - delta));
+ if ( kwHelp <= 0.0 )
+ {
AliFatal("Error in calculating Dalitz mass histogram binning!");
- fLPMass = 0;
- }
- krollWada = (2.0 / mLL) * TMath::Exp(1.5 * TMath::Log(kwHelp))
- * TMath::Sqrt(1.0 - 4.0 * epsilon / q)
- * (1.0 + 2.0 * epsilon / q);
- formFactor =
+ }
+ krollWada = (2.0 / mLL) * TMath::Exp(1.5 * TMath::Log(kwHelp))
+ * TMath::Sqrt(1.0 - 4.0 * epsilon / q)
+ * (1.0 + 2.0 * epsilon / q);
+ formFactor =
(TMath::Power(TMath::Power(0.6519,2),2))
/ (TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL, 2), 2)
+ TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
- weight = krollWada * formFactor;
- fLPMass->AddBinContent(ibin, weight);
+ weight = krollWada * formFactor;
+ if ( idecay == 1 )
+ fEPMass->AddBinContent(ibin, weight);
+ else
+ fMPMass->AddBinContent(ibin, weight);
+ }
}
}
-void AliOmegaDalitz::Decay(Int_t /*idpart*/, TLorentzVector* pparent)
+void AliOmegaDalitz::Decay(Int_t idlepton, TLorentzVector* pparent)
{
//-----------------------------------------------------------------------------
//
//
//-----------------------------------------------------------------------------
+ if (!fInit) {
+ Init();
+ fInit=1;
+ }
+
Double_t pmass, lmass, omass, lpmass;
Double_t e1, p1, e3, p3;
- Double_t betaSquare, lambda;
Double_t costheta, sintheta, cosphi, sinphi, phi;
// Get the particle masses
- // electron
- lmass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
+ // lepton
+ lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
// omega
- pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
+ pmass = pparent->M();
// pi0
omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
// Sample the lepton pair mass from a histogram
for( ;; )
{
- lpmass = fLPMass->GetRandom();
- if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
+ if ( idlepton == kElectron )
+ lpmass = fEPMass->GetRandom();
+ else
+ lpmass = fMPMass->GetRandom();
+ if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
}
// lepton pair kinematics in virtual photon rest frame
e1 = lpmass / 2.;
p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
- betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
- lambda = betaSquare / (2.0 - betaSquare);
+ // betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
+ // lambda = betaSquare / (2.0 - betaSquare);
costheta = (2.0 * gRandom->Rndm()) - 1.;
sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
{
// Import TParticles in array particles
- // e+ e- pi0
+ // l+ l- pi0
TClonesArray &clonesParticles = *particles;
- Int_t pdg [3] = {kPi0, kElectron, -kElectron};
- Int_t parent[3] = {-1, 0, 0};
- Int_t d1 [3] = {1, -1, -1};
- Int_t d2 [3] = {2, -1, -1};
+ Int_t pdg [3] = {kElectron, -kElectron, kPi0};
+ if ( fProducts[1].M() > 0.001 )
+ {
+ pdg[0] = kMuonPlus;
+ pdg[1] = -kMuonPlus;
+ }
+ Int_t parent[3] = {0, 0, -1};
+ Int_t d1 [3] = {-1, -1, 1};
+ Int_t d2 [3] = {-1, -1, 2};
for (Int_t i = 2; i > -1; i--) {
Double_t px = fProducts[i].Px();
Double_t py = fProducts[i].Py();
}
-
void AliOmegaDalitz::
Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
- Double_t cosphi, Double_t sinphi)
+ Double_t cosphi, Double_t sinphi) const
{
+// Perform rotation
pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
for (Int_t j = 0; j < 3; j++) dp[j] = (TParticle*) (array->At(fd+j));
if ((dp[0]->GetPdgCode() != kPi0) ||
- (TMath::Abs(dp[1]->GetPdgCode()) != kElectron)) continue;
+ ((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
+ (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
- Decay(223, &omega);
+ if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron )
+ Decay(kElectron, &omega);
+ else
+ Decay(kMuonPlus, &omega);
for (Int_t j = 0; j < 3; j++) dp[j]->SetMomentum(fProducts[2-j]);
printf("original %13.3f %13.3f %13.3f %13.3f \n",
part->Px(), part->Py(), part->Pz(), part->Energy());
AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
: AliDecayer(),
- fLPMass(0)
+ fEPMass(0),
+ fMPMass(0),
+ fInit(0)
{
// Copy constructor
dalitz.Copy(*this);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff