Updates
[u/mrichter/AliRoot.git] / EVGEN / AliDecayerPolarized.cxx
CommitLineData
c1e10a2a 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#include <TMath.h>
17#include <TF1.h>
18#include <TDatabasePDG.h>
19#include <TParticlePDG.h>
20#include <TParticle.h>
21#include <TVector3.h>
22#include <TRandom.h>
23#include <TObjectTable.h>
24
25#include "AliDecayerPolarized.h"
26#include "AliLog.h"
27
706938e6 28ClassImp(AliDecayerPolarized)
1c56e311 29
30
c1e10a2a 31
32//____________________________________________________________
1c56e311 33AliDecayerPolarized::AliDecayerPolarized():
34 fAlpha(0),
35 fSystRef(kHelicity),
36 fDecProd(kMuon),
f0de6c5c 37 fPol(0),
1c56e311 38 fMother(0),
39 fDaughter1(0),
40 fDaughter2(0)
c1e10a2a 41{
42// Default constructor
f0de6c5c 43
c1e10a2a 44}
45
46//____________________________________________________________
1c56e311 47AliDecayerPolarized::AliDecayerPolarized(Double_t alpha, Polar_t systref, FinState_t decprod):
48 fAlpha(alpha),
49 fSystRef(systref),
50 fDecProd(decprod),
51 fPol(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)),
52 fMother(0),
53 fDaughter1(0),
54 fDaughter2(0)
c1e10a2a 55{
56// Another constructor
c1e10a2a 57 fPol->SetParameter(0,fAlpha);
58 if(fDecProd!=kMuon && fDecProd!=kElectron)
1c56e311 59 AliFatal("Only polarized decay into muons or electrons is implemented!");
60}
61
62AliDecayerPolarized::AliDecayerPolarized(const AliDecayerPolarized &decayer):
63 AliDecayer(decayer),
64 fAlpha(0),
65 fSystRef(kHelicity),
66 fDecProd(kMuon),
67 fPol(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)),
68 fMother(0),
69 fDaughter1(0),
70 fDaughter2(0)
71{
72 // Copy Constructor
73 decayer.Copy(*this);
c1e10a2a 74}
75
76//____________________________________________________________
77AliDecayerPolarized::~AliDecayerPolarized()
78{
79// Destructor
78fd4e85 80 delete fPol;
81 delete fMother;
82 delete fDaughter1;
83 delete fDaughter2;
c1e10a2a 84}
85//____________________________________________________________
86void AliDecayerPolarized::Decay(Int_t ipart, TLorentzVector *p)
87{
88// Polarized 2- body decay
89 TDatabasePDG *pDataBase = TDatabasePDG::Instance();
90 if(ipart!= (pDataBase->GetParticle("J/psi")->PdgCode()) &&
91 ipart!= (pDataBase->GetParticle("psi'")->PdgCode()) &&
92 ipart!= (pDataBase->GetParticle("Upsilon")->PdgCode()) &&
93 ipart!= (pDataBase->GetParticle("Upsilon'")->PdgCode()))
94 AliFatal("Polarized decay only implemented for J/psi, psi', Upsilon and Upsilon' !");
95
96 TParticlePDG *d1 = 0, *d2 = 0;
97 if(fDecProd==kMuon){
98 d1 = pDataBase->GetParticle("mu-");
99 d2 = pDataBase->GetParticle("mu+");
100 }
101 else if(fDecProd==kElectron){
102 d1 = pDataBase->GetParticle("e-");
103 d2 = pDataBase->GetParticle("e+");
104 }
105// energies and momenta in lab frame
106 Double_t e1 = p->M() / 2.;
107 Double_t e2 = e1;
108 Double_t p1 = TMath::Sqrt((e1 + d1->Mass())*(e1 - d1->Mass()));
109
110 Double_t costheta = fPol->GetRandom();
111
112 Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
113 Double_t phi = 2. * TMath::Pi() * gRandom->Rndm();
114 Double_t px1 = p1 * sintheta * TMath::Cos(phi);
115 Double_t py1 = p1 * sintheta * TMath::Sin(phi);
116 Double_t pz1 = p1 * costheta;
117
118 TLorentzVector v1,v2, boosted1, boosted2;
119 v1.SetPxPyPzE(-px1,-py1,-pz1,e1); //in the dimuon rest frame
120 v2.SetPxPyPzE(px1,py1,pz1,e2);
121
78fd4e85 122 TLorentzVector pProj, pTarg; // In the CM frame
c1e10a2a 123 Float_t mp = 0.938;
78fd4e85 124 pProj.SetPxPyPzE(0.,0.,-7000.,TMath::Sqrt(7000.*7000.+mp*mp)); // projectile
125 pTarg.SetPxPyPzE(0.,0.,7000.,TMath::Sqrt(7000.*7000.+mp*mp)); // target
c1e10a2a 126
127 TVector3 betajpsicm;
128 betajpsicm = (-1./p->E()*p->Vect());
129
130 if(fSystRef == kHelicity) {
131 // polarization axis: direction of the JPsi in the CM
132 TVector3 v3jpsi = (p->Vect()).Unit();
133 v1.RotateUz(v3jpsi);
134 v2.RotateUz(v3jpsi);
135 } else if (fSystRef == kColSop){
136 // polarization axis: bisector of proj and target in the dimuon rest frame
78fd4e85 137 pProj.Boost(betajpsicm); //boost proj and targ from CM to DIMU rest frame
138 pTarg.Boost(betajpsicm);
c1e10a2a 139
140 TVector3 zaxisCS;
78fd4e85 141 zaxisCS=(((pProj.Vect()).Unit())-((pTarg.Vect()).Unit())).Unit();
c1e10a2a 142 v1.RotateUz(zaxisCS);
143 v2.RotateUz(zaxisCS);
144 }
145
146// printf("v1 components (mu1 with jpsi at rest)%f %f %f %f\n",v1.Px(),v1.Py(),v1.Pz(),v1.E());
147// printf("v2 components (mu2 with jpsi at rest)%f %f %f %f\n",v2.Px(),v2.Py(),v2.Pz(),v2.E());
148
149 v1.Boost(-betajpsicm); //boost muons from DIMUON rest frame to CM
150 v2.Boost(-betajpsicm);
151
152// printf("v1 components (mu1 in polar. ref. syst.) %f %f %f %f\n",v1.Px(),v1.Py(),v1.Pz(),v1.E());
153// printf("v2 components (mu2 in polar. ref. syst.) %f %f %f %f\n",v2.Px(),v2.Py(),v2.Pz(),v2.E());
154
78fd4e85 155 Int_t statusDecayed=11;
156 Int_t statusUndecayed=1;
c1e10a2a 157
158
78fd4e85 159 fMother = new TParticle(ipart,statusDecayed,0,-1,2,3,p->Px(),p->Py(),p->Pz(),p->E(),0.,0.,0.,0);
160 fDaughter1 = new TParticle(d1->PdgCode(),statusUndecayed,1,-1,0,0,v1.Px(),v1.Py(),v1.Pz(),v1.E(),0.,0.,0.,0);
161 fDaughter2 = new TParticle(d2->PdgCode(),statusUndecayed,1,-1,0,0,v2.Px(),v2.Py(),v2.Pz(),v2.E(),0.,0.,0.,0);
c1e10a2a 162
163}
164
165//____________________________________________________________
166Int_t AliDecayerPolarized::ImportParticles(TClonesArray *part)
167{
168// Return array of particles
169 TClonesArray &cloneparticles = *part;
170 cloneparticles.Clear();
171
172 new(cloneparticles[0])TParticle(*fMother);
173 new(cloneparticles[1])TParticle(*fDaughter1);
174 new(cloneparticles[2])TParticle(*fDaughter2);
175
176 return part->GetEntries();
177}
8a878ecf 178
179void AliDecayerPolarized::Copy(TObject &) const
180{
181 //
182 // Copy *this onto AliDecayerPolarized -- not implemented
183 //
184 Fatal("Copy","Not implemented!\n");
185}
186
3cad937a 187void AliDecayerPolarized::SetForceDecay(Int_t)
188{
189 // This method is dummy
3cad937a 190}
191
192void AliDecayerPolarized::ForceDecay()
193{
194 // This method is dummy
195 AliWarning("Method not implemented for this class !\n");
196}
197
198Float_t AliDecayerPolarized::GetPartialBranchingRatio(Int_t)
199{
200 // This method is dummy
7442901c 201 return 1.;
3cad937a 202}
203
204Float_t AliDecayerPolarized::GetLifetime(Int_t)
205{
206 // This method is dummy
207 AliWarning("Method not implemented for this class !\n");
208 return -1.;
209}
210
211void AliDecayerPolarized::ReadDecayTable()
212{
213 // This method is dummy
214 AliWarning("Method not implemented for this class !\n");
215}
216