Warning messages removed.
[u/mrichter/AliRoot.git] / EVGEN / AliDecayerPolarized.cxx
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
28 ClassImp(AliDecayerPolarized)
29
30 //____________________________________________________________
31 AliDecayerPolarized::AliDecayerPolarized()
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 //____________________________________________________________
45 AliDecayerPolarized::AliDecayerPolarized(Double_t alpha, Polar_t systref, FinState_t decprod)
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 //____________________________________________________________
61 AliDecayerPolarized::~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 //____________________________________________________________
70 void 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 //____________________________________________________________
150 Int_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 }
162
163 void  AliDecayerPolarized::Copy(TObject &) const
164 {
165     //
166     // Copy *this onto AliDecayerPolarized -- not implemented
167     //
168     Fatal("Copy","Not implemented!\n");
169 }
170
171 void AliDecayerPolarized::SetForceDecay(Int_t)
172 {
173     // This method is dummy
174 }
175
176 void AliDecayerPolarized::ForceDecay()
177 {
178     // This method is dummy
179     AliWarning("Method not implemented for this class !\n");
180 }
181
182 Float_t AliDecayerPolarized::GetPartialBranchingRatio(Int_t)
183 {
184     // This method is dummy
185     return  1.;
186 }
187
188 Float_t AliDecayerPolarized::GetLifetime(Int_t)
189 {
190     // This method is dummy
191     AliWarning("Method not implemented for this class !\n");
192     return -1.;
193 }
194
195 void AliDecayerPolarized::ReadDecayTable()
196 {
197     // This method is dummy
198     AliWarning("Method not implemented for this class !\n");
199 }
200