--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+#include <TMath.h>
+#include <TF1.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+#include <TParticle.h>
+#include <TVector3.h>
+#include <TRandom.h>
+#include <TObjectTable.h>
+
+#include "AliDecayerPolarized.h"
+#include "AliLog.h"
+
+ClassImp(AliDecayerPolarized)
+
+//____________________________________________________________
+AliDecayerPolarized::AliDecayerPolarized():AliDecayer()
+{
+// Default constructor
+ fAlpha=0;
+ fSystRef=kHelicity;
+ fDecProd=kMuon;
+ fMother=0;
+ fDaughter1=0;
+ fDaughter2=0;
+ fPol = new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.);
+ fPol->SetParameter(0,fAlpha);
+}
+
+//____________________________________________________________
+AliDecayerPolarized::AliDecayerPolarized(Double_t alpha, Polar_t systref, FinState_t decprod):AliDecayer()
+{
+// Another constructor
+ fAlpha=alpha;
+ fSystRef=systref;
+ fDecProd=decprod;
+ fMother=0;
+ fDaughter1=0;
+ fDaughter2=0;
+ fPol = new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.);
+ fPol->SetParameter(0,fAlpha);
+ if(fDecProd!=kMuon && fDecProd!=kElectron)
+ AliFatal("Only polarized decay into muons or electrons is implemented!");
+}
+
+//____________________________________________________________
+AliDecayerPolarized::~AliDecayerPolarized()
+{
+// Destructor
+ if(fPol) delete fPol; fPol=0;
+ if(fMother) delete fMother; fMother=0;
+ if(fDaughter1) delete fDaughter1; fDaughter1=0;
+ if(fDaughter2) delete fDaughter2; fDaughter2=0;
+}
+//____________________________________________________________
+void AliDecayerPolarized::Decay(Int_t ipart, TLorentzVector *p)
+{
+// Polarized 2- body decay
+ TDatabasePDG *pDataBase = TDatabasePDG::Instance();
+ if(ipart!= (pDataBase->GetParticle("J/psi")->PdgCode()) &&
+ ipart!= (pDataBase->GetParticle("psi'")->PdgCode()) &&
+ ipart!= (pDataBase->GetParticle("Upsilon")->PdgCode()) &&
+ ipart!= (pDataBase->GetParticle("Upsilon'")->PdgCode()))
+ AliFatal("Polarized decay only implemented for J/psi, psi', Upsilon and Upsilon' !");
+
+ TParticlePDG *d1 = 0, *d2 = 0;
+ if(fDecProd==kMuon){
+ d1 = pDataBase->GetParticle("mu-");
+ d2 = pDataBase->GetParticle("mu+");
+ }
+ else if(fDecProd==kElectron){
+ d1 = pDataBase->GetParticle("e-");
+ d2 = pDataBase->GetParticle("e+");
+ }
+// energies and momenta in lab frame
+ Double_t e1 = p->M() / 2.;
+ Double_t e2 = e1;
+ Double_t p1 = TMath::Sqrt((e1 + d1->Mass())*(e1 - d1->Mass()));
+
+ Double_t costheta = fPol->GetRandom();
+
+ Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
+ Double_t phi = 2. * TMath::Pi() * gRandom->Rndm();
+ Double_t px1 = p1 * sintheta * TMath::Cos(phi);
+ Double_t py1 = p1 * sintheta * TMath::Sin(phi);
+ Double_t pz1 = p1 * costheta;
+
+ TLorentzVector v1,v2, boosted1, boosted2;
+ v1.SetPxPyPzE(-px1,-py1,-pz1,e1); //in the dimuon rest frame
+ v2.SetPxPyPzE(px1,py1,pz1,e2);
+
+ TLorentzVector PProj, PTarg; // In the CM frame
+ Float_t mp = 0.938;
+ PProj.SetPxPyPzE(0.,0.,-7000.,TMath::Sqrt(7000.*7000.+mp*mp)); // projectile
+ PTarg.SetPxPyPzE(0.,0.,7000.,TMath::Sqrt(7000.*7000.+mp*mp)); // target
+
+ TVector3 betajpsicm;
+ betajpsicm = (-1./p->E()*p->Vect());
+
+ if(fSystRef == kHelicity) {
+ // polarization axis: direction of the JPsi in the CM
+ TVector3 v3jpsi = (p->Vect()).Unit();
+ v1.RotateUz(v3jpsi);
+ v2.RotateUz(v3jpsi);
+ } else if (fSystRef == kColSop){
+ // polarization axis: bisector of proj and target in the dimuon rest frame
+ PProj.Boost(betajpsicm); //boost proj and targ from CM to DIMU rest frame
+ PTarg.Boost(betajpsicm);
+
+ TVector3 zaxisCS;
+ zaxisCS=(((PProj.Vect()).Unit())-((PTarg.Vect()).Unit())).Unit();
+ v1.RotateUz(zaxisCS);
+ v2.RotateUz(zaxisCS);
+ }
+
+// printf("v1 components (mu1 with jpsi at rest)%f %f %f %f\n",v1.Px(),v1.Py(),v1.Pz(),v1.E());
+// printf("v2 components (mu2 with jpsi at rest)%f %f %f %f\n",v2.Px(),v2.Py(),v2.Pz(),v2.E());
+
+ v1.Boost(-betajpsicm); //boost muons from DIMUON rest frame to CM
+ v2.Boost(-betajpsicm);
+
+// printf("v1 components (mu1 in polar. ref. syst.) %f %f %f %f\n",v1.Px(),v1.Py(),v1.Pz(),v1.E());
+// printf("v2 components (mu2 in polar. ref. syst.) %f %f %f %f\n",v2.Px(),v2.Py(),v2.Pz(),v2.E());
+
+ Int_t status_decayed=11;
+ Int_t status_undecayed=1;
+
+
+ fMother = new TParticle(ipart,status_decayed,0,-1,2,3,p->Px(),p->Py(),p->Pz(),p->E(),0.,0.,0.,0);
+ fDaughter1 = new TParticle(d1->PdgCode(),status_undecayed,1,-1,0,0,v1.Px(),v1.Py(),v1.Pz(),v1.E(),0.,0.,0.,0);
+ fDaughter2 = new TParticle(d2->PdgCode(),status_undecayed,1,-1,0,0,v2.Px(),v2.Py(),v2.Pz(),v2.E(),0.,0.,0.,0);
+
+}
+
+//____________________________________________________________
+Int_t AliDecayerPolarized::ImportParticles(TClonesArray *part)
+{
+// Return array of particles
+ TClonesArray &cloneparticles = *part;
+ cloneparticles.Clear();
+
+ new(cloneparticles[0])TParticle(*fMother);
+ new(cloneparticles[1])TParticle(*fDaughter1);
+ new(cloneparticles[2])TParticle(*fDaughter2);
+
+ return part->GetEntries();
+}
--- /dev/null
+#ifndef ALIDECAYERPOLARIZED_H
+#define ALIDECAYERPOLARIZED_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Class to generate decay products for polarized heavy quarkonia
+
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include <TF1.h>
+
+#include "AliDecayer.h"
+
+
+class AliDecayerPolarized : public AliDecayer
+{
+ public:
+ typedef enum { kNoPol = 0, kColSop = 1, kHelicity = 2} Polar_t;
+ typedef enum { kElectron = 1, kMuon = 2} FinState_t;
+ AliDecayerPolarized();
+ AliDecayerPolarized(Double_t alpha, Polar_t systref, FinState_t decprod);
+ virtual ~AliDecayerPolarized();
+ void SetPolDec(Double_t alpha=0) {fAlpha=alpha;}
+ void SetPolRefSys(Polar_t systref=kColSop) {fSystRef=systref;}
+ void SetDecProd(FinState_t decprod=kMuon) {fDecProd=decprod;}
+ virtual void Init(){;}
+ virtual void Decay(Int_t ipart, TLorentzVector *p);
+ virtual Int_t ImportParticles(TClonesArray *part);
+ AliDecayerPolarized(const AliDecayerPolarized &decayer):AliDecayer(decayer)
+ {decayer.Copy(*this);}
+ virtual AliDecayerPolarized &operator=(const AliDecayerPolarized &decayer)
+ {decayer.Copy(*this);return(*this);}
+
+ protected:
+ Double_t fAlpha; // Polarization parameter
+ Polar_t fSystRef; // Reference system for polarization
+ FinState_t fDecProd; // Choice of decay products
+ TF1 *fPol; // ! Angular distribution for decay products
+ TParticle *fMother; // ! Particle that has to be decayed
+ TParticle *fDaughter1; // ! Decay product no. 1
+ TParticle *fDaughter2; // ! Decay product no. 2
+
+ ClassDef(AliDecayerPolarized,1) // Polarized 2-body quarkonium decay
+};
+#endif
+
+