[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(),
        fEPMass(0),
	fMPMass(0)
{
    // Constructor
}

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

  // Get the particle masses
  // electron
  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();

  for ( idecay = 1; idecay<=2; idecay++ )
  {
    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!");
	if ( idecay == 1 ) 
	  fEPMass = 0;
	else
	  fMPMass = 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!");
	if ( idecay == 1 ) 
	  fEPMass = 0;
	else
	  fMPMass = 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;   
      if ( idecay == 1 ) 
	fEPMass->AddBinContent(ibin, weight);
      else
	fMPMass->AddBinContent(ibin, weight);
    }
  }
}


void AliOmegaDalitz::Decay(Int_t idlepton, 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
  // lepton
  lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
  // omega
  pmass = pparent->M();
  // pi0
  omass = (TDatabasePDG::Instance()->GetParticle(kPi0))     ->Mass();

  // Sample the lepton pair mass from a histogram
  for( ;; ) 
  {
    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);
  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
    // l+ l- pi0
    
    TClonesArray &clonesParticles = *particles;

    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();
	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;
}

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)
{
  // Copy constructor
  dalitz.Copy(*this);
}
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(),
0439b2b7	51	fEPMass(0),
0439b2b7	52	fMPMass(0)
8c90f5a2	53	{
	54	// Constructor
	55	}
	56
	57	void AliOmegaDalitz::Init()
	58	{
	59	// Initialisation
0439b2b7	60	Int_t idecay, ibin, nbins = 1000;
8c90f5a2	61	Double_t min, max, binwidth;
0439b2b7	62	Double_t pmass, lmass, omass, emass, mmass;
8c90f5a2	63	Double_t epsilon, delta, mLL, q, kwHelp, krollWada, formFactor, weight;
	64
	65	// Get the particle masses
	66	// electron
0439b2b7	67	emass = (TDatabasePDG::Instance()->GetParticle(kElectron))->Mass();
	68	// muon
	69	mmass = (TDatabasePDG::Instance()->GetParticle(kMuonPlus))->Mass();
8c90f5a2	70	// omega
	71	pmass = (TDatabasePDG::Instance()->GetParticle(223)) ->Mass();
	72	// pi0
	73	omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
	74
0439b2b7	75	for ( idecay = 1; idecay<=2; idecay++ )
8c90f5a2	76	{
0439b2b7	77	if ( idecay == 1 )
	78	lmass = emass;
	79	else
	80	lmass = mmass;
	81
	82	min = 2.0 * lmass;
	83	max = pmass - omass;
	84	binwidth = (max - min) / (Double_t)nbins;
	85	if ( idecay == 1 )
	86	fEPMass = new TH1F("fEPMass", "Dalitz electron pair", nbins, min, max);
	87	else
	88	fMPMass = new TH1F("fMPMass", "Dalitz muon pair", nbins, min, max);
	89
	90	epsilon = (lmass / pmass) * (lmass / pmass);
	91	delta = (omass / pmass) * (omass / pmass);
	92
	93	for ( ibin = 1; ibin <= nbins; ibin++ )
8c90f5a2	94	{
0439b2b7	95	mLL = min + (Double_t)(ibin - 1) * binwidth + binwidth / 2.0;
	96	q = (mLL / pmass) * (mLL / pmass);
	97	if ( q <= 4.0 * epsilon )
	98	{
8c90f5a2	99	AliFatal("Error in calculating Dalitz mass histogram binning!");
0439b2b7	100	if ( idecay == 1 )
	101	fEPMass = 0;
	102	else
	103	fMPMass = 0;
	104	}
	105
	106	kwHelp = (1.0 + q / (1.0 - delta)) * (1.0 + q / (1.0 - delta))
	107	- 4.0 * q / ((1.0 - delta) * (1.0 - delta));
	108	if ( kwHelp <= 0.0 )
	109	{
8c90f5a2	110	AliFatal("Error in calculating Dalitz mass histogram binning!");
0439b2b7	111	if ( idecay == 1 )
	112	fEPMass = 0;
	113	else
	114	fMPMass = 0;
	115	}
	116	krollWada = (2.0 / mLL) * TMath::Exp(1.5 * TMath::Log(kwHelp))
	117	* TMath::Sqrt(1.0 - 4.0 * epsilon / q)
	118	* (1.0 + 2.0 * epsilon / q);
	119	formFactor =
8c90f5a2	120	(TMath::Power(TMath::Power(0.6519,2),2))
	121	/ (TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL, 2), 2)
	122	+ TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
0439b2b7	123	weight = krollWada * formFactor;
	124	if ( idecay == 1 )
	125	fEPMass->AddBinContent(ibin, weight);
	126	else
	127	fMPMass->AddBinContent(ibin, weight);
	128	}
8c90f5a2	129	}
	130	}
	131
	132
0439b2b7	133	void AliOmegaDalitz::Decay(Int_t idlepton, TLorentzVector* pparent)
8c90f5a2	134	{
	135	//-----------------------------------------------------------------------------
	136	//
	137	// Generate Omega Dalitz decay
	138	//
	139	//-----------------------------------------------------------------------------
	140
	141	Double_t pmass, lmass, omass, lpmass;
	142	Double_t e1, p1, e3, p3;
	143	Double_t betaSquare, lambda;
	144	Double_t costheta, sintheta, cosphi, sinphi, phi;
	145
	146	// Get the particle masses
0439b2b7	147	// lepton
0439b2b7	148	lmass = (TDatabasePDG::Instance()->GetParticle(idlepton))->Mass();
8c90f5a2	149	// omega
d060703c	150	pmass = pparent->M();
8c90f5a2	151	// pi0
	152	omass = (TDatabasePDG::Instance()->GetParticle(kPi0)) ->Mass();
	153
	154	// Sample the lepton pair mass from a histogram
	155	for( ;; )
	156	{
0439b2b7	157	if ( idlepton == kElectron )
	158	lpmass = fEPMass->GetRandom();
	159	else
	160	lpmass = fMPMass->GetRandom();
	161	if ( pmass - omass > lpmass && lpmass / 2. > lmass ) break;
8c90f5a2	162	}
	163
	164	// lepton pair kinematics in virtual photon rest frame
	165	e1 = lpmass / 2.;
	166	p1 = TMath::Sqrt((e1 + lmass) * (e1 - lmass));
	167	betaSquare = 1.0 - 4.0 * (lmass * lmass) / (lpmass * lpmass);
	168	lambda = betaSquare / (2.0 - betaSquare);
	169	costheta = (2.0 * gRandom->Rndm()) - 1.;
	170	sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
	171	phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
	172	sinphi = TMath::Sin(phi);
	173	cosphi = TMath::Cos(phi);
	174	// momentum vectors of leptons in virtual photon rest frame
	175	Double_t pProd1[3] = {p1 * sintheta * cosphi,
	176	p1 * sintheta * sinphi,
	177	p1 * costheta};
	178	Double_t pProd2[3] = {-1.0 * p1 * sintheta * cosphi,
	179	-1.0 * p1 * sintheta * sinphi,
	180	-1.0 * p1 * costheta};
	181
	182	// pizero kinematics in omega rest frame
	183	e3 = (pmass * pmass + omass * omass - lpmass * lpmass)/(2. * pmass);
	184	p3 = TMath::Sqrt((e3 + omass) * (e3 - omass));
	185	costheta = (2.0 * gRandom->Rndm()) - 1.;
	186	sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
	187	phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
	188	sinphi = TMath::Sin(phi);
	189	cosphi = TMath::Cos(phi);
	190	// pizero 4-vector in omega rest frame
	191	fProducts[2].SetPx(p3 * sintheta * cosphi);
	192	fProducts[2].SetPy(p3 * sintheta * sinphi);
	193	fProducts[2].SetPz(p3 * costheta);
	194	fProducts[2].SetE(e3);
	195
	196	// lepton 4-vectors in properly rotated virtual photon rest frame
	197	Double_t pRot1[3] = {0.};
	198	Rot(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
	199	Double_t pRot2[3] = {0.};
	200	Rot(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi);
	201	fProducts[0].SetPx(pRot1[0]);
	202	fProducts[0].SetPy(pRot1[1]);
	203	fProducts[0].SetPz(pRot1[2]);
	204	fProducts[0].SetE(e1);
	205	fProducts[1].SetPx(pRot2[0]);
	206	fProducts[1].SetPy(pRot2[1]);
	207	fProducts[1].SetPz(pRot2[2]);
	208	fProducts[1].SetE(e1);
	209
	210	// boost the dilepton into the omega's rest frame
	211	Double_t eLPparent = TMath::Sqrt(p3 * p3 + lpmass * lpmass);
	212	TVector3 boostPair( -1.0 * fProducts[2].Px() / eLPparent,
	213	-1.0 * fProducts[2].Py() / eLPparent,
	214	-1.0 * fProducts[2].Pz() / eLPparent);
	215	fProducts[0].Boost(boostPair);
	216	fProducts[1].Boost(boostPair);
	217
	218	// boost all decay products into the lab frame
	219	TVector3 boostLab( pparent->Px() / pparent->E(),
	220	pparent->Py() / pparent->E(),
	221	pparent->Pz() / pparent->E());
	222	fProducts[0].Boost(boostLab);
	223	fProducts[1].Boost(boostLab);
	224	fProducts[2].Boost(boostLab);
	225
226	return;
227
228	}
229
230	Int_t AliOmegaDalitz::ImportParticles(TClonesArray *particles)
231	{
232	// Import TParticles in array particles
0439b2b7	233	// l+ l- pi0
8c90f5a2	234
	235	TClonesArray &clonesParticles = *particles;
	236
0439b2b7	237	Int_t pdg [3] = {kElectron, -kElectron, kPi0};
	238	if ( fProducts[1].M() > 0.001 )
	239	{
	240	pdg[0] = kMuonPlus;
	241	pdg[1] = -kMuonPlus;
	242	}
	243	Int_t parent[3] = {0, 0, -1};
	244	Int_t d1 [3] = {-1, -1, 1};
	245	Int_t d2 [3] = {-1, -1, 2};
8c90f5a2	246	for (Int_t i = 2; i > -1; i--) {
	247	Double_t px = fProducts[i].Px();
	248	Double_t py = fProducts[i].Py();
	249	Double_t pz = fProducts[i].Pz();
	250	Double_t e = fProducts[i].E();
	251	//
	252	new(clonesParticles[2 - i]) TParticle(pdg[i], 1, parent[i], -1, d1[i], d2[i], px, py, pz, e, 0., 0., 0., 0.);
	253	}
	254	return (3);
	255	}
	256
	257
8c90f5a2	258	void AliOmegaDalitz::
	259	Rot(Double_t pin[3], Double_t pout[3], Double_t costheta, Double_t sintheta,
	260	Double_t cosphi, Double_t sinphi)
	261	{
	262	pout[0] = pin[0]costhetacosphi-pin[1]sinphi+pin[2]sintheta*cosphi;
	263	pout[1] = pin[0]costhetasinphi+pin[1]cosphi+pin[2]sintheta*sinphi;
	264	pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta;
	265	return;
	266	}
	267
05b093f3	268	void AliOmegaDalitz::Decay(TClonesArray* array)
	269	{
	270	//
	271	// Replace all omega dalitz decays with the correct matrix element decays
	272	//
	273	printf("-->Correcting Dalitz \n");
	274	Int_t nt = array->GetEntriesFast();
	275	TParticle* dp[3];
	276	for (Int_t i = 0; i < nt; i++) {
	277	TParticle* part = (TParticle*) (array->At(i));
	278	if (part->GetPdgCode() != 223) continue;
	279
	280	Int_t fd = part->GetFirstDaughter() - 1;
	281	Int_t ld = part->GetLastDaughter() - 1;
	282
	283	if (fd < 0) continue;
	284	if ((ld - fd) != 2) continue;
	285
	286	for (Int_t j = 0; j < 3; j++) dp[j] = (TParticle*) (array->At(fd+j));
	287	if ((dp[0]->GetPdgCode() != kPi0) \|\|
0439b2b7	288	((TMath::Abs(dp[1]->GetPdgCode()) != kElectron) &&
0439b2b7	289	(TMath::Abs(dp[1]->GetPdgCode()) != kMuonPlus))) continue;
05b093f3	290	TLorentzVector omega(part->Px(), part->Py(), part->Pz(), part->Energy());
0439b2b7	291	if ( TMath::Abs(dp[1]->GetPdgCode()) == kElectron )
	292	Decay(kElectron, &omega);
	293	else
	294	Decay(kMuonPlus, &omega);
05b093f3	295	for (Int_t j = 0; j < 3; j++) dp[j]->SetMomentum(fProducts[2-j]);
	296	printf("original %13.3f %13.3f %13.3f %13.3f \n",
	297	part->Px(), part->Py(), part->Pz(), part->Energy());
	298	printf("new %13.3f %13.3f %13.3f %13.3f \n",
	299	dp[0]->Px()+dp[1]->Px()+dp[2]->Px(),
	300	dp[0]->Py()+dp[1]->Py()+dp[2]->Py(),
	301	dp[0]->Pz()+dp[1]->Pz()+dp[2]->Pz(),
	302	dp[0]->Energy()+dp[1]->Energy()+dp[2]->Energy());
	303
	304
	305	}
	306	}
	307	AliOmegaDalitz& AliOmegaDalitz::operator=(const AliOmegaDalitz& rhs)
	308	{
	309	// Assignment operator
	310	rhs.Copy(*this);
	311	return *this;
	312	}
	313
	314	void AliOmegaDalitz::Copy(TObject&) const
	315	{
	316	//
	317	// Copy
	318	//
	319	Fatal("Copy","Not implemented!\n");
	320	}
	321
	322	AliOmegaDalitz::AliOmegaDalitz(const AliOmegaDalitz &dalitz)
	323	: AliDecayer(),
0439b2b7	324	fEPMass(0),
0439b2b7	325	fMPMass(0)
05b093f3	326	{
	327	// Copy constructor
	328	dalitz.Copy(*this);
	329	}