[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoModelWeightGeneratorBasic.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoModelWeightGeneratorBasic -  basic femtoscopic weight generator  ///
/// only return a simple                                                          ///
/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////
#ifdef __ROOT__
  ClassImp(AliFemtoModelWeightGeneratorBasic, 1)
#endif

#include "AliFemtoModelWeightGeneratorBasic.h"
#include "AliFemtoModelHiddenInfo.h"

//________________________
AliFemtoModelWeightGeneratorBasic::AliFemtoModelWeightGeneratorBasic():
  AliFemtoModelWeightGenerator()
{
  /* no-op */
}
//________________________
AliFemtoModelWeightGeneratorBasic::AliFemtoModelWeightGeneratorBasic(const AliFemtoModelWeightGeneratorBasic &aModel) :
  AliFemtoModelWeightGenerator(aModel)
{
  /* no-op */
}

//________________________
AliFemtoModelWeightGeneratorBasic::~AliFemtoModelWeightGeneratorBasic()
{
  /* no-op */
}

AliFemtoModelWeightGeneratorBasic& AliFemtoModelWeightGeneratorBasic::operator=(const AliFemtoModelWeightGeneratorBasic &aModel)
{
  if (this != &aModel) {
    AliFemtoModelWeightGenerator::operator=(aModel);
  }
  
  return *this;
}

//________________________
Double_t AliFemtoModelWeightGeneratorBasic::GenerateWeight(AliFemtoPair *aPair)
{
  // Generate a simple femtoscopic weight coming only from 
  // quantum statistics - symmetrization or anti-symmetrization
  // of the pair wave function

  // Get hidden information pointers
  AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
  AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();

  // Calculate pair variables
  Double_t tPx = inf1->GetTrueMomentum()->x()+inf2->GetTrueMomentum()->x();
  Double_t tPy = inf1->GetTrueMomentum()->y()+inf2->GetTrueMomentum()->y();
  Double_t tPz = inf1->GetTrueMomentum()->z()+inf2->GetTrueMomentum()->z();
  //  double tE  = inf1->GetTrueMomentum()->e +inf2->GetTrueMomentum()->.e;
  Double_t tM1 = inf1->GetMass();
  Double_t tM2 = inf2->GetMass();
  Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->Mag2());
  Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->Mag2());
  Double_t tE  = tE1 + tE2;
  Double_t tPt = tPx*tPx + tPy*tPy;
  Double_t tMt = tE*tE - tPz*tPz;//mCVK;
  Double_t tM  = sqrt(tMt - tPt);
  tMt = sqrt(tMt);
  tPt = sqrt(tPt);
  Double_t tBetat = tPt/tMt;

  // Boost to LCMS
  Double_t tBeta = tPz/tE;
  Double_t tGamma = tE/tMt;	    
  fKStarLong = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
  Double_t tE1L = tGamma * (tE1  - tBeta * inf1->GetTrueMomentum()->z());
  
  // Transform positions to LCMS
//   Double_t tP1zl = tGamma * (inf1->GetEmissionPoint()->z() - tBeta * inf1->GetEmissionPoint()->t());
//   Double_t tP1tl = tGamma * (inf1->GetEmissionPoint()->t() - tBeta * inf1->GetEmissionPoint()->z());
  
//   Double_t tP2zl = tGamma * (inf2->GetEmissionPoint()->z() - tBeta * inf2->GetEmissionPoint()->t());
//   Double_t tP2tl = tGamma * (inf2->GetEmissionPoint()->t() - tBeta * inf2->GetEmissionPoint()->z());
  
//   Double_t tP1pzl = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
//   Double_t tP1el  = tGamma * (tE1  - tBeta * inf1->GetTrueMomentum()->z());
  
//   Double_t tP2pzl = tGamma * (inf2->GetTrueMomentum()->z() - tBeta * tE2);
//   Double_t tP2el  = tGamma * (tE2  - tBeta * inf2->GetTrueMomentum()->z());
  
  // Rotate in transverse plane
  fKStarOut  = ( inf1->GetTrueMomentum()->x()*tPx + inf1->GetTrueMomentum()->y()*tPy)/tPt;
  fKStarSide = (-inf1->GetTrueMomentum()->x()*tPy + inf1->GetTrueMomentum()->y()*tPx)/tPt;
      
//   Double_t tP1pxl = fKStarOut;
//   Double_t tP1pyl = fKStarSide;
  
//   Double_t tP2pxl = (inf2->GetTrueMomentum()->x()*tPx + inf2->GetTrueMomentum()->y()*tPy)/tPt;
//   Double_t tP2pyl = (inf2->GetTrueMomentum()->y()*tPx - inf2->GetTrueMomentum()->x()*tPy)/tPt;;
  
//   Double_t tKO = tP1pxl - tP2pxl;
//   Double_t tKS = tP1pyl - tP2pyl;
//   Double_t tKL = tP1pzl - tP2pzl;
//   Double_t tDE = tP1el  - tP2el;
  
  // save the rotated coordinates in LCMS variables
//   Double_t tP1xl = ( inf1->GetEmissionPoint()->x()*tPx + inf1->GetEmissionPoint()->y()*tPy)/tPt;
//   Double_t tP1yl = (-inf1->GetEmissionPoint()->x()*tPy + inf1->GetEmissionPoint()->y()*tPx)/tPt;

//   Double_t tP2xl = ( inf2->GetEmissionPoint()->x()*tPx + inf2->GetEmissionPoint()->y()*tPy)/tPt;
//   Double_t tP2yl = (-inf2->GetEmissionPoint()->x()*tPy + inf2->GetEmissionPoint()->y()*tPx)/tPt;
  
  // Boost to pair cms
  fKStarOut = tMt/tM * (fKStarOut - tPt/tMt * tE1L);
  
  tBetat = tPt/tMt;
//   Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
  
//   Double_t tP1xp = tGammat*(tP1xl - tBetat*tP1tl);
//   Double_t tP1tp = tGammat*(tP1tl - tBetat*tP1xl);
  
//   Double_t tP2xp = tGammat*(tP2xl - tBetat*tP2tl);
//   Double_t tP2tp = tGammat*(tP2tl - tBetat*tP2xl);
  
//   Double_t tRO = (tP1xl - tP2xl)/0.197327;
//   Double_t tRS = (tP1yl - tP2yl)/0.197327;
//   Double_t tRL = (tP1zl - tP2zl)/0.197327;
//   Double_t tDT = (tP1tl - tP2tl)/0.197327;
  
  Double_t tDX = inf1->GetEmissionPoint()->x()-inf2->GetEmissionPoint()->x();
  Double_t tDY = inf1->GetEmissionPoint()->y()-inf2->GetEmissionPoint()->y();
  Double_t tRLong = inf1->GetEmissionPoint()->z()-inf2->GetEmissionPoint()->z();
  Double_t tDTime = inf1->GetEmissionPoint()->t()-inf2->GetEmissionPoint()->t();

  Double_t tROut = (tDX*tPx + tDY*tPy)/tPt;
  Double_t tRSide = (-tDX*tPy + tDY*tPx)/tPt;

  fRStarSide = tRSide;
  Double_t tRSS = fRStarSide/0.197327;

  fRStarLong = tGamma*(tRLong - tBeta* tDTime);
  Double_t tDTimePairLCMS = tGamma*(tDTime - tBeta* tRLong);

  Double_t tRLS = fRStarLong/0.197327;
  tBeta = tPt/tMt;
  tGamma = tMt/tM;

  fRStarOut = tGamma*(tROut - tBeta* tDTimePairLCMS);
  Double_t tROS = fRStarOut/0.197327;
//   Double_t tDTimePairCMS = tGamma*(tDTimePairLCMS - tBeta* tROut);
  fRStar = ::sqrt(fRStarOut*fRStarOut + fRStarSide*fRStarSide +
			   fRStarLong*fRStarLong);
  fKStar = ::sqrt(fKStarOut*fKStarOut + fKStarSide*fKStarSide + fKStarLong*fKStarLong);
//   Double_t tRSt = fRStar/0.197327;

  if (fPairType != fgkPairTypeNone) {
    if ((fPairType == PionPlusPionPlus()) || (fPairType == KaonPlusKaonPlus()))
      return 1.0 + cos (2*(fKStarOut * tROS + fKStarSide * tRSS + fKStarLong * tRLS));
    else if (fPairType == ProtonProton())
      return 1.0 - 0.5 * cos (2*(fKStarOut * tROS + fKStarSide * tRSS + fKStarLong * tRLS));
    else 
      return 1.0;
  }
  else {
    Int_t tPairType = GetPairTypeFromPair(aPair);
    if ((tPairType == PionPlusPionPlus()) || (tPairType == KaonPlusKaonPlus()))
      return 1.0 + cos (2*(fKStarOut * tROS + fKStarSide * tRSS + fKStarLong * tRLS));
    else if (tPairType == ProtonProton())
      return 1.0 - 0.5 * cos (2*(fKStarOut * tROS + fKStarSide * tRSS + fKStarLong * tRLS));
    else 
      return 1.0;
    
  }
}

//________________________
void     AliFemtoModelWeightGeneratorBasic::SetPairType(Int_t aPairType)
{
  AliFemtoModelWeightGenerator::SetPairType(aPairType);
}
//________________________
void     AliFemtoModelWeightGeneratorBasic::SetPairTypeFromPair(AliFemtoPair *aPair)
{
  AliFemtoModelWeightGenerator::SetPairTypeFromPair(aPair);
}
//________________________
Int_t    AliFemtoModelWeightGeneratorBasic::GetPairType() const
{
  return AliFemtoModelWeightGenerator::GetPairType();
}
//________________________
AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorBasic::Clone() const
{
  return GetGenerator();
}
//________________________
AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorBasic::GetGenerator() const
{
  AliFemtoModelWeightGeneratorBasic *tGen = new AliFemtoModelWeightGeneratorBasic(*this);
  return tGen;
}
Commit	Line	Data
75c432a7	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoModelWeightGeneratorBasic - basic femtoscopic weight generator ///
	4	/// only return a simple ///
	5	/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
	6	/// ///
	7	////////////////////////////////////////////////////////////////////////////////
	8	#ifdef __ROOT__
	9	ClassImp(AliFemtoModelWeightGeneratorBasic, 1)
	10	#endif
	11
	12	#include "AliFemtoModelWeightGeneratorBasic.h"
	13	#include "AliFemtoModelHiddenInfo.h"
	14
	15	//________________________
	16	AliFemtoModelWeightGeneratorBasic::AliFemtoModelWeightGeneratorBasic():
	17	AliFemtoModelWeightGenerator()
	18	{
	19	/* no-op */
	20	}
	21	//________________________
	22	AliFemtoModelWeightGeneratorBasic::AliFemtoModelWeightGeneratorBasic(const AliFemtoModelWeightGeneratorBasic &aModel) :
	23	AliFemtoModelWeightGenerator(aModel)
	24	{
	25	/* no-op */
	26	}
	27
	28	//________________________
	29	AliFemtoModelWeightGeneratorBasic::~AliFemtoModelWeightGeneratorBasic()
	30	{
	31	/* no-op */
	32	}
	33
c66f18ec	34	AliFemtoModelWeightGeneratorBasic& AliFemtoModelWeightGeneratorBasic::operator=(const AliFemtoModelWeightGeneratorBasic &aModel)
	35	{
	36	if (this != &aModel) {
	37	AliFemtoModelWeightGenerator::operator=(aModel);
	38	}
	39
	40	return *this;
	41	}
	42
75c432a7	43	//________________________
	44	Double_t AliFemtoModelWeightGeneratorBasic::GenerateWeight(AliFemtoPair *aPair)
	45	{
	46	// Generate a simple femtoscopic weight coming only from
	47	// quantum statistics - symmetrization or anti-symmetrization
	48	// of the pair wave function
	49
	50	// Get hidden information pointers
d0e92d9a	51	AliFemtoModelHiddenInfo inf1 = (AliFemtoModelHiddenInfo ) aPair->Track1()->HiddenInfo();
d0e92d9a	52	AliFemtoModelHiddenInfo inf2 = (AliFemtoModelHiddenInfo ) aPair->Track2()->HiddenInfo();
75c432a7	53
	54	// Calculate pair variables
	55	Double_t tPx = inf1->GetTrueMomentum()->x()+inf2->GetTrueMomentum()->x();
	56	Double_t tPy = inf1->GetTrueMomentum()->y()+inf2->GetTrueMomentum()->y();
	57	Double_t tPz = inf1->GetTrueMomentum()->z()+inf2->GetTrueMomentum()->z();
	58	// double tE = inf1->GetTrueMomentum()->e +inf2->GetTrueMomentum()->.e;
	59	Double_t tM1 = inf1->GetMass();
	60	Double_t tM2 = inf2->GetMass();
69c1c8ff	61	Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->Mag2());
69c1c8ff	62	Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->Mag2());
75c432a7	63	Double_t tE = tE1 + tE2;
	64	Double_t tPt = tPxtPx + tPytPy;
	65	Double_t tMt = tEtE - tPztPz;//mCVK;
	66	Double_t tM = sqrt(tMt - tPt);
	67	tMt = sqrt(tMt);
	68	tPt = sqrt(tPt);
	69	Double_t tBetat = tPt/tMt;
	70
	71	// Boost to LCMS
	72	Double_t tBeta = tPz/tE;
	73	Double_t tGamma = tE/tMt;
	74	fKStarLong = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
	75	Double_t tE1L = tGamma * (tE1 - tBeta * inf1->GetTrueMomentum()->z());
	76
	77	// Transform positions to LCMS
	78	// Double_t tP1zl = tGamma * (inf1->GetEmissionPoint()->z() - tBeta * inf1->GetEmissionPoint()->t());
	79	// Double_t tP1tl = tGamma * (inf1->GetEmissionPoint()->t() - tBeta * inf1->GetEmissionPoint()->z());
	80
	81	// Double_t tP2zl = tGamma * (inf2->GetEmissionPoint()->z() - tBeta * inf2->GetEmissionPoint()->t());
	82	// Double_t tP2tl = tGamma * (inf2->GetEmissionPoint()->t() - tBeta * inf2->GetEmissionPoint()->z());
	83
	84	// Double_t tP1pzl = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
	85	// Double_t tP1el = tGamma * (tE1 - tBeta * inf1->GetTrueMomentum()->z());
	86
	87	// Double_t tP2pzl = tGamma * (inf2->GetTrueMomentum()->z() - tBeta * tE2);
	88	// Double_t tP2el = tGamma * (tE2 - tBeta * inf2->GetTrueMomentum()->z());
	89
	90	// Rotate in transverse plane
	91	fKStarOut = ( inf1->GetTrueMomentum()->x()tPx + inf1->GetTrueMomentum()->y()tPy)/tPt;
	92	fKStarSide = (-inf1->GetTrueMomentum()->x()tPy + inf1->GetTrueMomentum()->y()tPx)/tPt;
	93
	94	// Double_t tP1pxl = fKStarOut;
	95	// Double_t tP1pyl = fKStarSide;
	96
	97	// Double_t tP2pxl = (inf2->GetTrueMomentum()->x()tPx + inf2->GetTrueMomentum()->y()tPy)/tPt;
	98	// Double_t tP2pyl = (inf2->GetTrueMomentum()->y()tPx - inf2->GetTrueMomentum()->x()tPy)/tPt;;
	99
	100	// Double_t tKO = tP1pxl - tP2pxl;
	101	// Double_t tKS = tP1pyl - tP2pyl;
	102	// Double_t tKL = tP1pzl - tP2pzl;
	103	// Double_t tDE = tP1el - tP2el;
	104
	105	// save the rotated coordinates in LCMS variables
	106	// Double_t tP1xl = ( inf1->GetEmissionPoint()->x()tPx + inf1->GetEmissionPoint()->y()tPy)/tPt;
	107	// Double_t tP1yl = (-inf1->GetEmissionPoint()->x()tPy + inf1->GetEmissionPoint()->y()tPx)/tPt;
	108
	109	// Double_t tP2xl = ( inf2->GetEmissionPoint()->x()tPx + inf2->GetEmissionPoint()->y()tPy)/tPt;
	110	// Double_t tP2yl = (-inf2->GetEmissionPoint()->x()tPy + inf2->GetEmissionPoint()->y()tPx)/tPt;
	111
	112	// Boost to pair cms
	113	fKStarOut = tMt/tM * (fKStarOut - tPt/tMt * tE1L);
	114
	115	tBetat = tPt/tMt;
	116	// Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
	117
	118	// Double_t tP1xp = tGammat(tP1xl - tBetattP1tl);
	119	// Double_t tP1tp = tGammat(tP1tl - tBetattP1xl);
	120
	121	// Double_t tP2xp = tGammat(tP2xl - tBetattP2tl);
	122	// Double_t tP2tp = tGammat(tP2tl - tBetattP2xl);
	123
	124	// Double_t tRO = (tP1xl - tP2xl)/0.197327;
	125	// Double_t tRS = (tP1yl - tP2yl)/0.197327;
	126	// Double_t tRL = (tP1zl - tP2zl)/0.197327;
127	// Double_t tDT = (tP1tl - tP2tl)/0.197327;
128
129	Double_t tDX = inf1->GetEmissionPoint()->x()-inf2->GetEmissionPoint()->x();
130	Double_t tDY = inf1->GetEmissionPoint()->y()-inf2->GetEmissionPoint()->y();
131	Double_t tRLong = inf1->GetEmissionPoint()->z()-inf2->GetEmissionPoint()->z();
132	Double_t tDTime = inf1->GetEmissionPoint()->t()-inf2->GetEmissionPoint()->t();
133
134	Double_t tROut = (tDXtPx + tDYtPy)/tPt;
135	Double_t tRSide = (-tDXtPy + tDYtPx)/tPt;
136
137	fRStarSide = tRSide;
138	Double_t tRSS = fRStarSide/0.197327;
139
140	fRStarLong = tGamma(tRLong - tBeta tDTime);
141	Double_t tDTimePairLCMS = tGamma(tDTime - tBeta tRLong);
142
143	Double_t tRLS = fRStarLong/0.197327;
144	tBeta = tPt/tMt;
145	tGamma = tMt/tM;
146
147	fRStarOut = tGamma(tROut - tBeta tDTimePairLCMS);
148	Double_t tROS = fRStarOut/0.197327;
149	// Double_t tDTimePairCMS = tGamma(tDTimePairLCMS - tBeta tROut);
150	fRStar = ::sqrt(fRStarOutfRStarOut + fRStarSidefRStarSide +
151	fRStarLong*fRStarLong);
152	fKStar = ::sqrt(fKStarOutfKStarOut + fKStarSidefKStarSide + fKStarLong*fKStarLong);
153	// Double_t tRSt = fRStar/0.197327;
154
0b3bd1ac	155	if (fPairType != fgkPairTypeNone) {
	156	if ((fPairType == PionPlusPionPlus()) \|\| (fPairType == KaonPlusKaonPlus()))
	157	return 1.0 + cos (2(fKStarOut tROS + fKStarSide * tRSS + fKStarLong * tRLS));
	158	else if (fPairType == ProtonProton())
	159	return 1.0 - 0.5 * cos (2(fKStarOut tROS + fKStarSide * tRSS + fKStarLong * tRLS));
	160	else
	161	return 1.0;
	162	}
	163	else {
	164	Int_t tPairType = GetPairTypeFromPair(aPair);
	165	if ((tPairType == PionPlusPionPlus()) \|\| (tPairType == KaonPlusKaonPlus()))
	166	return 1.0 + cos (2(fKStarOut tROS + fKStarSide * tRSS + fKStarLong * tRLS));
	167	else if (tPairType == ProtonProton())
	168	return 1.0 - 0.5 * cos (2(fKStarOut tROS + fKStarSide * tRSS + fKStarLong * tRLS));
	169	else
	170	return 1.0;
	171
	172	}
75c432a7	173	}
	174
	175	//________________________
	176	void AliFemtoModelWeightGeneratorBasic::SetPairType(Int_t aPairType)
	177	{
	178	AliFemtoModelWeightGenerator::SetPairType(aPairType);
	179	}
	180	//________________________
	181	void AliFemtoModelWeightGeneratorBasic::SetPairTypeFromPair(AliFemtoPair *aPair)
	182	{
	183	AliFemtoModelWeightGenerator::SetPairTypeFromPair(aPair);
	184	}
	185	//________________________
fcda1d4e	186	Int_t AliFemtoModelWeightGeneratorBasic::GetPairType() const
75c432a7	187	{
	188	return AliFemtoModelWeightGenerator::GetPairType();
	189	}
	190	//________________________
	191	AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorBasic::Clone() const
	192	{
	193	return GetGenerator();
	194	}
	195	//________________________
	196	AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorBasic::GetGenerator() const
	197	{
	198	AliFemtoModelWeightGeneratorBasic tGen = new AliFemtoModelWeightGeneratorBasic(this);
	199	return tGen;
	200	}