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