//
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);
+ for ( idecay = 1; idecay<=2; idecay++ )
+ {
+ if ( idecay == 1 )
+ lmass = emass;
+ else
+ lmass = mmass;
- epsilon = (lmass / pmass) * (lmass / pmass);
- delta = (omass / pmass) * (omass / pmass);
+ 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);
- 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 )
+ 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;
- printf(" %5d %13.3f \n", ibin, weight);
-
- 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;
return;
}
+void AliOmegaDalitz::Decay(TClonesArray* array)
+{
+ //
+ // Replace all omega dalitz decays with the correct matrix element decays
+ //
+ printf("-->Correcting Dalitz \n");
+ Int_t nt = array->GetEntriesFast();
+ TParticle* dp[3];
+ for (Int_t i = 0; i < nt; i++) {
+ TParticle* part = (TParticle*) (array->At(i));
+ if (part->GetPdgCode() != 223) continue;
+
+ Int_t fd = part->GetFirstDaughter() - 1;
+ Int_t ld = part->GetLastDaughter() - 1;
+
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+
+ 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) &&
+ (TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
+ TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
+ 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());
+ printf("new %13.3f %13.3f %13.3f %13.3f \n",
+ dp[0]->Px()+dp[1]->Px()+dp[2]->Px(),
+ dp[0]->Py()+dp[1]->Py()+dp[2]->Py(),
+ dp[0]->Pz()+dp[1]->Pz()+dp[2]->Pz(),
+ dp[0]->Energy()+dp[1]->Energy()+dp[2]->Energy());
+
+
+ }
+}
+AliOmegaDalitz& AliOmegaDalitz::operator=(const AliOmegaDalitz& rhs)
+{
+// Assignment operator
+ rhs.Copy(*this);
+ return *this;
+}
+
+void AliOmegaDalitz::Copy(TObject&) const
+{
+ //
+ // Copy
+ //
+ Fatal("Copy","Not implemented!\n");
+}
+
+AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
+ : AliDecayer(),
+ fEPMass(0),
+ fMPMass(0),
+ fInit(0)
+{
+ // Copy constructor
+ dalitz.Copy(*this);
+}