[u/mrichter/AliRoot.git] / EVGEN / AliDecayerPolarized.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include <TMath.h>
#include <TF1.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TParticle.h>
#include <TVector3.h>
#include <TRandom.h>
#include <TObjectTable.h>

#include "AliDecayerPolarized.h"
#include "AliLog.h"

ClassImp(AliDecayerPolarized)

//____________________________________________________________
AliDecayerPolarized::AliDecayerPolarized()
{
// Default constructor
    fAlpha=0;
    fSystRef=kHelicity;
    fDecProd=kMuon;
    fMother=0;
    fDaughter1=0;
    fDaughter2=0;
    fPol = new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.);
    fPol->SetParameter(0,fAlpha);
}

//____________________________________________________________
AliDecayerPolarized::AliDecayerPolarized(Double_t alpha, Polar_t systref, FinState_t decprod)
{
// Another constructor
    fAlpha=alpha;
    fSystRef=systref;
    fDecProd=decprod;
    fMother=0;
    fDaughter1=0;
    fDaughter2=0;
    fPol = new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.);
    fPol->SetParameter(0,fAlpha);
    if(fDecProd!=kMuon && fDecProd!=kElectron)
    AliFatal("Only polarized decay into muons or electrons is implemented!");
}

//____________________________________________________________
AliDecayerPolarized::~AliDecayerPolarized()
{
// Destructor
    if(fPol) delete fPol; fPol=0;
    if(fMother) delete fMother; fMother=0;
    if(fDaughter1) delete fDaughter1; fDaughter1=0;
    if(fDaughter2) delete fDaughter2; fDaughter2=0;
}
//____________________________________________________________
void AliDecayerPolarized::Decay(Int_t ipart, TLorentzVector *p)
{
// Polarized 2- body decay 
    TDatabasePDG *pDataBase = TDatabasePDG::Instance();
    if(ipart!= (pDataBase->GetParticle("J/psi")->PdgCode()) &&
       ipart!= (pDataBase->GetParticle("psi'")->PdgCode())  &&
       ipart!= (pDataBase->GetParticle("Upsilon")->PdgCode())  &&
       ipart!= (pDataBase->GetParticle("Upsilon'")->PdgCode()))  
       AliFatal("Polarized decay only implemented for J/psi, psi', Upsilon and Upsilon' !");

    TParticlePDG *d1 = 0, *d2 = 0;
    if(fDecProd==kMuon){
      d1 = pDataBase->GetParticle("mu-");
      d2 = pDataBase->GetParticle("mu+");
    }
    else if(fDecProd==kElectron){
      d1 = pDataBase->GetParticle("e-");
      d2 = pDataBase->GetParticle("e+");
    }
// energies and momenta in lab frame 
    Double_t e1 = p->M() / 2.;
    Double_t e2 = e1; 
    Double_t p1 = TMath::Sqrt((e1 + d1->Mass())*(e1 - d1->Mass())); 

    Double_t costheta = fPol->GetRandom();
    
    Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
    Double_t phi      = 2. * TMath::Pi() * gRandom->Rndm(); 
    Double_t px1 = p1 * sintheta * TMath::Cos(phi); 
    Double_t py1 = p1 * sintheta * TMath::Sin(phi); 
    Double_t pz1 = p1 * costheta; 

    TLorentzVector v1,v2, boosted1, boosted2;	
    v1.SetPxPyPzE(-px1,-py1,-pz1,e1);   //in the dimuon rest frame
    v2.SetPxPyPzE(px1,py1,pz1,e2); 
	
    TLorentzVector PProj, PTarg; // In the CM frame
    Float_t mp = 0.938;
    PProj.SetPxPyPzE(0.,0.,-7000.,TMath::Sqrt(7000.*7000.+mp*mp)); // projectile
    PTarg.SetPxPyPzE(0.,0.,7000.,TMath::Sqrt(7000.*7000.+mp*mp)); // target
 
    TVector3 betajpsicm;
    betajpsicm = (-1./p->E()*p->Vect());
      
    if(fSystRef == kHelicity) {
      //  polarization axis: direction of the JPsi in the CM 
      TVector3 v3jpsi = (p->Vect()).Unit();  
      v1.RotateUz(v3jpsi);  
      v2.RotateUz(v3jpsi); 	
    } else if (fSystRef == kColSop){
      //  polarization axis: bisector of proj and target in the dimuon rest frame
      PProj.Boost(betajpsicm);   //boost proj and targ from CM to DIMU rest frame
      PTarg.Boost(betajpsicm);

      TVector3 zaxisCS;
      zaxisCS=(((PProj.Vect()).Unit())-((PTarg.Vect()).Unit())).Unit();
      v1.RotateUz(zaxisCS);
      v2.RotateUz(zaxisCS);
    }
 
//    printf("v1 components (mu1 with jpsi at rest)%f %f %f %f\n",v1.Px(),v1.Py(),v1.Pz(),v1.E());
//    printf("v2 components (mu2 with jpsi at rest)%f %f %f %f\n",v2.Px(),v2.Py(),v2.Pz(),v2.E());

    v1.Boost(-betajpsicm);  //boost muons from DIMUON rest frame to CM
    v2.Boost(-betajpsicm); 

//    printf("v1 components (mu1 in polar. ref. syst.) %f %f %f %f\n",v1.Px(),v1.Py(),v1.Pz(),v1.E());
//    printf("v2 components (mu2 in polar. ref. syst.) %f %f %f %f\n",v2.Px(),v2.Py(),v2.Pz(),v2.E());

    Int_t status_decayed=11;
    Int_t status_undecayed=1;
        
   
    fMother = new TParticle(ipart,status_decayed,0,-1,2,3,p->Px(),p->Py(),p->Pz(),p->E(),0.,0.,0.,0);   
    fDaughter1 = new TParticle(d1->PdgCode(),status_undecayed,1,-1,0,0,v1.Px(),v1.Py(),v1.Pz(),v1.E(),0.,0.,0.,0);   
    fDaughter2 = new TParticle(d2->PdgCode(),status_undecayed,1,-1,0,0,v2.Px(),v2.Py(),v2.Pz(),v2.E(),0.,0.,0.,0);

}

//____________________________________________________________
Int_t AliDecayerPolarized::ImportParticles(TClonesArray *part)
{
// Return array of particles  
  TClonesArray &cloneparticles = *part;
  cloneparticles.Clear();
  
  new(cloneparticles[0])TParticle(*fMother);
  new(cloneparticles[1])TParticle(*fDaughter1);
  new(cloneparticles[2])TParticle(*fDaughter2);

  return part->GetEntries();
}

void  AliDecayerPolarized::Copy(TObject &) const
{
    //
    // Copy *this onto AliDecayerPolarized -- not implemented
    //
    Fatal("Copy","Not implemented!\n");
}

void AliDecayerPolarized::SetForceDecay(Int_t)
{
    // This method is dummy
}

void AliDecayerPolarized::ForceDecay()
{
    // This method is dummy
    AliWarning("Method not implemented for this class !\n");
}

Float_t AliDecayerPolarized::GetPartialBranchingRatio(Int_t)
{
    // This method is dummy
    return  1.;
}

Float_t AliDecayerPolarized::GetLifetime(Int_t)
{
    // This method is dummy
    AliWarning("Method not implemented for this class !\n");
    return -1.;
}

void AliDecayerPolarized::ReadDecayTable()
{
    // This method is dummy
    AliWarning("Method not implemented for this class !\n");
}
Commit	Line	Data
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
	28	ClassImp(AliDecayerPolarized)
	29
	30	//____________________________________________________________
3cad937a	31	AliDecayerPolarized::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]xx",-1.,1.);
	41	fPol->SetParameter(0,fAlpha);
	42	}
	43
	44	//____________________________________________________________
3cad937a	45	AliDecayerPolarized::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]xx",-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.+mpmp)); // projectile
	109	PTarg.SetPxPyPzE(0.,0.,7000.,TMath::Sqrt(7000.7000.+mpmp)); // 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	}
8a878ecf	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
3cad937a	171	void AliDecayerPolarized::SetForceDecay(Int_t)
	172	{
	173	// This method is dummy
3cad937a	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
7442901c	185	return 1.;
3cad937a	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