[u/mrichter/AliRoot.git] / EVGEN / AliOmegaDalitz.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


/* $Id$ */


#include "AliOmegaDalitz.h"
#include <TMath.h>
#include <AliLog.h>
#include <TH1.h>
#include <TPDGCode.h>
#include <TRandom.h>
#include <TParticle.h>
#include <TDatabasePDG.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>

ClassImp(AliOmegaDalitz)


//-----------------------------------------------------------------------------
//
// Generate lepton-pair mass distributions for Dalitz decays according
// to the Kroll-Wada parametrization: N. Kroll, W. Wada: Phys. Rev 98(1955)1355
//
// For the electromagnetic form factor the parameterization from
// Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
//
//-----------------------------------------------------------------------------
//
// Adaption for AliRoot
//
// R. Averbeck
// A. Morsch
//
AliOmegaDalitz::AliOmegaDalitz():
	AliDecayer(),
	fLPMass(0)
{
    // Constructor
}

void AliOmegaDalitz::Init()
{
  // Initialisation
  Int_t    ibin, nbins = 1000;
  Double_t min, max, binwidth;
  Double_t pmass, lmass, omass;
  Double_t epsilon, delta, mLL, q, kwHelp, krollWada, formFactor, weight;

  // Get the particle masses
  // electron
  lmass = (TDatabasePDG::Instance()->GetParticle(kElectron))->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++ )
  {
    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 )
    {
	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 = 
	(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);
  }
}


void AliOmegaDalitz::Decay(Int_t /*idpart*/, TLorentzVector* pparent)
{
//-----------------------------------------------------------------------------
//
//  Generate Omega Dalitz decay
//
//-----------------------------------------------------------------------------

  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();
  // omega
  pmass = (TDatabasePDG::Instance()->GetParticle(223))      ->Mass();
  // 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;
  }
  
  // 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);
  costheta = (2.0 * gRandom->Rndm()) - 1.;
  sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
  phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
  sinphi   = TMath::Sin(phi);
  cosphi   = TMath::Cos(phi);
  // momentum vectors of leptons in virtual photon rest frame
  Double_t pProd1[3] = {p1 * sintheta * cosphi, 
		       p1 * sintheta * sinphi, 
		       p1 * costheta};
  Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi, 
		       -1.0 * p1 * sintheta * sinphi, 
		       -1.0 * p1 * costheta};
  
  // pizero kinematics in omega rest frame
  e3       = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
  p3       = TMath::Sqrt((e3 + omass)  * (e3 - omass));
  costheta = (2.0 * gRandom->Rndm()) - 1.;
  sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
  phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
  sinphi   = TMath::Sin(phi);
  cosphi   = TMath::Cos(phi);	    
  // pizero 4-vector in omega rest frame
  fProducts[2].SetPx(p3 * sintheta * cosphi);
  fProducts[2].SetPy(p3 * sintheta * sinphi);
  fProducts[2].SetPz(p3 * costheta);
  fProducts[2].SetE(e3);

  // lepton 4-vectors in properly rotated virtual photon rest frame
  Double_t pRot1[3] = {0.};
  Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
  Double_t pRot2[3] = {0.};
  Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi); 
  fProducts[0].SetPx(pRot1[0]);
  fProducts[0].SetPy(pRot1[1]);
  fProducts[0].SetPz(pRot1[2]);
  fProducts[0].SetE(e1);
  fProducts[1].SetPx(pRot2[0]);
  fProducts[1].SetPy(pRot2[1]);
  fProducts[1].SetPz(pRot2[2]);
  fProducts[1].SetE(e1);

  // boost the dilepton into the omega's rest frame 
  Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
  TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent, 
		      -1.0 * fProducts[2].Py() / eLPparent,
		      -1.0 * fProducts[2].Pz() / eLPparent);
  fProducts[0].Boost(boostPair);
  fProducts[1].Boost(boostPair);

  // boost all decay products into the lab frame 
  TVector3 boostLab( pparent->Px() / pparent->E(), 
		     pparent->Py() / pparent->E(),
		     pparent->Pz() / pparent->E());
  fProducts[0].Boost(boostLab);
  fProducts[1].Boost(boostLab);
  fProducts[2].Boost(boostLab);
  
  return;

}

Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
{
    // Import TParticles in array particles
    // e+ e- 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};
    for (Int_t i = 2; i > -1; i--) {
	Double_t px = fProducts[i].Px();
	Double_t py = fProducts[i].Py();
	Double_t pz = fProducts[i].Pz();
	Double_t e  = fProducts[i].E();
	//
	new(clonesParticles[2 - i]) TParticle(pdg[i], 1, parent[i], -1, d1[i], d2[i], px, py, pz, e, 0., 0., 0., 0.);
    }
    return (3);
}


void AliOmegaDalitz::
Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
    Double_t cosphi, Double_t sinphi)
{
  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;
}
Commit	Line	Data
8c90f5a2	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16
	17	/* $Id$ */
	18
	19
	20	#include "AliOmegaDalitz.h"
	21	#include <TMath.h>
	22	#include <AliLog.h>
	23	#include <TH1.h>
	24	#include <TPDGCode.h>
	25	#include <TRandom.h>
	26	#include <TParticle.h>
	27	#include <TDatabasePDG.h>
	28	#include <TLorentzVector.h>
	29	#include <TClonesArray.h>
	30
	31	ClassImp(AliOmegaDalitz)
	32
	33
	34	//-----------------------------------------------------------------------------
	35	//
	36	// Generate lepton-pair mass distributions for Dalitz decays according
	37	// to the Kroll-Wada parametrization: N. Kroll, W. Wada: Phys. Rev 98(1955)1355
	38	//
	39	// For the electromagnetic form factor the parameterization from
	40	// Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
	41	//
	42	//-----------------------------------------------------------------------------
	43	//
	44	// Adaption for AliRoot
	45	//
	46	// R. Averbeck
	47	// A. Morsch
	48	//
	49	AliOmegaDalitz::AliOmegaDalitz():
	50	AliDecayer(),
	51	fLPMass(0)
	52	{
	53	// Constructor
	54	}
	55
	56	void AliOmegaDalitz::Init()
	57	{
	58	// Initialisation
	59	Int_t ibin, nbins = 1000;
	60	Double_t min, max, binwidth;
	61	Double_t pmass, lmass, omass;
	62	Double_t epsilon, delta, mLL, q, kwHelp, krollWada, formFactor, weight;
	63
	64	// Get the particle masses
65	// electron
66	lmass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
67	// omega
68	pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
69	// pi0
70	omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
71
72	min = 2.0 * lmass;
73	max = pmass - omass;
74	binwidth = (max - min) / (Double_t)nbins;
75	fLPMass = new TH1F("fLPMass", "Dalitz", nbins, min, max);
76
77	epsilon = (lmass / pmass) * (lmass / pmass);
78	delta = (omass / pmass) * (omass / pmass);
79
80	for ( ibin = 1; ibin <= nbins; ibin++ )
81	{
82	mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
83	q = (mLL / pmass) * (mLL / pmass);
84	if ( q <= 4.0 * epsilon )
85	{
86	AliFatal("Error in calculating Dalitz mass histogram binning!");
87	fLPMass = 0;
88	}
89
90	kwHelp = (1.0 + q / (1.0 - delta)) * (1.0 + q / (1.0 - delta))
91	- 4.0 * q / ((1.0 - delta) * (1.0 - delta));
92	if ( kwHelp <= 0.0 )
93	{
94	AliFatal("Error in calculating Dalitz mass histogram binning!");
95	fLPMass = 0;
96	}
97	krollWada = (2.0 / mLL) * TMath::Exp(1.5 * TMath::Log(kwHelp))
98	* TMath::Sqrt(1.0 - 4.0 * epsilon / q)
99	* (1.0 + 2.0 * epsilon / q);
100	formFactor =
101	(TMath::Power(TMath::Power(0.6519,2),2))
102	/ (TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL, 2), 2)
103	+ TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
104	weight = krollWada * formFactor;
105	printf(" %5d %13.3f \n", ibin, weight);
106
107	fLPMass->AddBinContent(ibin, weight);
108	}
109	}
110
111
112	void AliOmegaDalitz::Decay(Int_t /idpart/, TLorentzVector* pparent)
113	{
114	//-----------------------------------------------------------------------------
115	//
116	// Generate Omega Dalitz decay
117	//
118	//-----------------------------------------------------------------------------
119
120	Double_t pmass, lmass, omass, lpmass;
121	Double_t e1, p1, e3, p3;
122	Double_t betaSquare, lambda;
123	Double_t costheta, sintheta, cosphi, sinphi, phi;
124
125	// Get the particle masses
126	// electron
127	lmass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
128	// omega
129	pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
130	// pi0
131	omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
132
133	// Sample the lepton pair mass from a histogram
134	for( ;; )
135	{
136	lpmass = fLPMass->GetRandom();
137	if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
138	}
139
140	// lepton pair kinematics in virtual photon rest frame
141	e1 = lpmass / 2.;
142	p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
143	betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
144	lambda = betaSquare / (2.0 - betaSquare);
145	costheta = (2.0 * gRandom->Rndm()) - 1.;
146	sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
147	phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
148	sinphi = TMath::Sin(phi);
149	cosphi = TMath::Cos(phi);
150	// momentum vectors of leptons in virtual photon rest frame
151	Double_t pProd1[3] = {p1 * sintheta * cosphi,
152	p1 * sintheta * sinphi,
153	p1 * costheta};
154	Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi,
155	-1.0 * p1 * sintheta * sinphi,
156	-1.0 * p1 * costheta};
157
158	// pizero kinematics in omega rest frame
159	e3 = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
160	p3 = TMath::Sqrt((e3 + omass) * (e3 - omass));
161	costheta = (2.0 * gRandom->Rndm()) - 1.;
162	sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
163	phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
164	sinphi = TMath::Sin(phi);
165	cosphi = TMath::Cos(phi);
166	// pizero 4-vector in omega rest frame
167	fProducts[2].SetPx(p3 * sintheta * cosphi);
168	fProducts[2].SetPy(p3 * sintheta * sinphi);
169	fProducts[2].SetPz(p3 * costheta);
170	fProducts[2].SetE(e3);
171
172	// lepton 4-vectors in properly rotated virtual photon rest frame
173	Double_t pRot1[3] = {0.};
174	Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
175	Double_t pRot2[3] = {0.};
176	Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi);
177	fProducts[0].SetPx(pRot1[0]);
178	fProducts[0].SetPy(pRot1[1]);
179	fProducts[0].SetPz(pRot1[2]);
180	fProducts[0].SetE(e1);
181	fProducts[1].SetPx(pRot2[0]);
182	fProducts[1].SetPy(pRot2[1]);
183	fProducts[1].SetPz(pRot2[2]);
184	fProducts[1].SetE(e1);
185
186	// boost the dilepton into the omega's rest frame
187	Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
188	TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent,
189	-1.0 * fProducts[2].Py() / eLPparent,
190	-1.0 * fProducts[2].Pz() / eLPparent);
191	fProducts[0].Boost(boostPair);
192	fProducts[1].Boost(boostPair);
193
194	// boost all decay products into the lab frame
195	TVector3 boostLab( pparent->Px() / pparent->E(),
196	pparent->Py() / pparent->E(),
197	pparent->Pz() / pparent->E());
198	fProducts[0].Boost(boostLab);
199	fProducts[1].Boost(boostLab);
200	fProducts[2].Boost(boostLab);
201
202	return;
203
204	}
205
206	Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
207	{
208	// Import TParticles in array particles
209	// e+ e- pi0
210
211	TClonesArray &clonesParticles = *particles;
212
213	Int_t pdg [3] = {kPi0, kElectron, -kElectron};
214	Int_t parent[3] = {-1, 0, 0};
215	Int_t d1 [3] = {1, -1, -1};
216	Int_t d2 [3] = {2, -1, -1};
217	for (Int_t i = 2; i > -1; i--) {
218	Double_t px = fProducts[i].Px();
219	Double_t py = fProducts[i].Py();
220	Double_t pz = fProducts[i].Pz();
221	Double_t e = fProducts[i].E();
222	//
223	new(clonesParticles[2 - i]) TParticle(pdg[i], 1, parent[i], -1, d1[i], d2[i], px, py, pz, e, 0., 0., 0., 0.);
224	}
225	return (3);
226	}
227
228
229
230	void AliOmegaDalitz::
231	Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
232	Double_t cosphi, Double_t sinphi)
233	{
234	pout[0] = pin[0]costhetacosphi-pin[1]sinphi+pin[2]sintheta*cosphi;
235	pout[1] = pin[0]costhetasinphi+pin[1]cosphi+pin[2]sintheta*sinphi;
236	pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
237	return;
238	}
239