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(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)),
  38     fMother(0),
  39     fDaughter1(0),
  40     fDaughter2(0)
  41 {
  42 // Default constructor
  43     fPol->SetParameter(0,fAlpha);
  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     if(fPol) delete fPol; fPol=0;
  81     if(fMother) delete fMother; fMother=0;
  82     if(fDaughter1) delete fDaughter1; fDaughter1=0;
  83     if(fDaughter2) delete fDaughter2; fDaughter2=0;
  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 status_decayed=11;
 156     Int_t status_undecayed=1;
 157
 158
 159     fMother = new TParticle(ipart,status_decayed,0,-1,2,3,p->Px(),p->Py(),p->Pz(),p->E(),0.,0.,0.,0);
 160     fDaughter1 = new TParticle(d1->PdgCode(),status_undecayed,1,-1,0,0,v1.Px(),v1.Py(),v1.Pz(),v1.E(),0.,0.,0.,0);
 161     fDaughter2 = new TParticle(d2->PdgCode(),status_undecayed,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