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