EMCAL/DCAL Trigger Mapping for Run 2
[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
32 //____________________________________________________________
33 AliDecayerPolarized::AliDecayerPolarized():
34     fAlpha(0),
35     fSystRef(kHelicity),
36     fDecProd(kMuon),
37     fPol(0),
38     fMother(0),
39     fDaughter1(0),
40     fDaughter2(0)
41 {
42 // Default constructor
43
44 }
45
46 //____________________________________________________________
47 AliDecayerPolarized::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)
55 {
56 // Another constructor
57     fPol->SetParameter(0,fAlpha);
58     if(fDecProd!=kMuon && fDecProd!=kElectron)
59         AliFatal("Only polarized decay into muons or electrons is implemented!");
60 }
61
62 AliDecayerPolarized::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);
74 }
75
76 //____________________________________________________________
77 AliDecayerPolarized::~AliDecayerPolarized()
78 {
79 // Destructor
80     delete fPol; 
81     delete fMother;
82     delete fDaughter1;
83     delete fDaughter2;
84 }
85 //____________________________________________________________
86 void 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         
122     TLorentzVector pProj, pTarg; // In the CM frame
123     Float_t mp = 0.938;
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
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
137       pProj.Boost(betajpsicm);   //boost proj and targ from CM to DIMU rest frame
138       pTarg.Boost(betajpsicm);
139
140       TVector3 zaxisCS;
141       zaxisCS=(((pProj.Vect()).Unit())-((pTarg.Vect()).Unit())).Unit();
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
155     Int_t statusDecayed=11;
156     Int_t statusUndecayed=1;
157         
158    
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);
162
163 }
164
165 //____________________________________________________________
166 Int_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 }
178
179 void  AliDecayerPolarized::Copy(TObject &) const
180 {
181     //
182     // Copy *this onto AliDecayerPolarized -- not implemented
183     //
184     Fatal("Copy","Not implemented!\n");
185 }
186
187 void AliDecayerPolarized::SetForceDecay(Int_t)
188 {
189     // This method is dummy
190 }
191
192 void AliDecayerPolarized::ForceDecay()
193 {
194     // This method is dummy
195     AliWarning("Method not implemented for this class !\n");
196 }
197
198 Float_t AliDecayerPolarized::GetPartialBranchingRatio(Int_t)
199 {
200     // This method is dummy
201     return  1.;
202 }
203
204 Float_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
211 void AliDecayerPolarized::ReadDecayTable()
212 {
213     // This method is dummy
214     AliWarning("Method not implemented for this class !\n");
215 }
216