[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoModelGausRinvFreezeOutGenerator.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoModelGausRinvFreezeOutGenerator - freeze-out                     ///
/// coordinates generator, generating a 3D gaussian ellipsoid in LCMS        ///
/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////
#ifdef __ROOT__
  ClassImp(AliFemtoModelGausRinvFreezeOutGenerator, 1)
#endif

#include "math.h"
#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
#include "AliFemtoModelHiddenInfo.h"
#include "AliFemtoModelGlobalHiddenInfo.h"
#include "AliFemtoLorentzVector.h"

//_______________________
AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
  fSizeInv(0),
  fSelectPrimary(false)
{
  // Default constructor
  fRandom = new TRandom2();
}

//_______________________
AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
  AliFemtoModelFreezeOutGenerator(),
  fSizeInv(0),
  fSelectPrimary(false)
{
  // Copy constructor
  fRandom = new TRandom2();
  SetSizeInv(aModel.GetSizeInv());
}
//_______________________
AliFemtoModelGausRinvFreezeOutGenerator::~AliFemtoModelGausRinvFreezeOutGenerator()
{
  if (fRandom) delete fRandom;
}
//_______________________
AliFemtoModelGausRinvFreezeOutGenerator& AliFemtoModelGausRinvFreezeOutGenerator::operator=(const AliFemtoModelGausRinvFreezeOutGenerator &aModel)
{
  if (this != &aModel) {
    fRandom = new TRandom2();
    SetSizeInv(aModel.GetSizeInv());
  }

  return *this;

}
//_______________________
void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
{
  AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
  AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();

  if ((!inf1) || (!inf2)) { cout << "Hidden info not created! "  << endl; exit(kFALSE); }

  if (fSelectPrimary) {
    const AliFemtoModelGlobalHiddenInfo *infg1 = dynamic_cast<const AliFemtoModelGlobalHiddenInfo *>(aPair->Track1()->HiddenInfo());
    const AliFemtoModelGlobalHiddenInfo *infg2 = dynamic_cast<const AliFemtoModelGlobalHiddenInfo *>(aPair->Track2()->HiddenInfo());
    
    if ((infg1) && (infg2)) {
      // assume the emission point is in [cm] and try to judge if
      // both particles are primary
      Double_t dist1 = infg1->GetGlobalEmissionPoint()->Perp();
      Double_t dist2 = infg2->GetGlobalEmissionPoint()->Perp();

      if ((dist1 > 0.05) && (dist2 > 0.05)) {
	// At least one particle is non primary
	if (!(inf1->GetEmissionPoint())) {
	  AliFemtoLorentzVector tPos(-1000,1000,-500,0);
	  inf1->SetEmissionPoint(&tPos);
	}
	else
	  inf1->SetEmissionPoint(-1000,1000,-500,0);
	if (!(inf2->GetEmissionPoint())) {
	  AliFemtoLorentzVector tPos(fRandom->Gaus(0,1000.0),fRandom->Gaus(0,1000),fRandom->Gaus(0,1000),0.0);
	  inf2->SetEmissionPoint(&tPos);
	}
	else
	  inf2->SetEmissionPoint(fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), 0.0);
	
	return;
      }
    }
  }

  // Generate two particle emission points with respect
  // to their pair momentum 
  // The source is the 3D Gaussian ellipsoid in the LCMS frame

  // Calculate sum momenta
  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_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 tEs = tE1 + tE2;

  Double_t tPt = sqrt(tPx*tPx + tPy*tPy);
  Double_t tMt = sqrt(tEs*tEs - tPz*tPz);

  // Generate positions in PRF from a Gaussian
  Double_t tROutS =  fRandom->Gaus(0,fSizeInv); // reuse of long
  Double_t tRSideS = fRandom->Gaus(0,fSizeInv);
  Double_t tRLongS = fRandom->Gaus(0,fSizeInv);
  Double_t tRTimeS = 0;
      
  Double_t tBetat  = tPt/tMt;
  Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);

  Double_t tBetaz  = tPz/tEs;
  Double_t tGammaz = 1.0/sqrt(1.0-tBetaz*tBetaz);

  Double_t tROut = tGammat * (tROutS + tBetat * tRTimeS);
  Double_t tDtL  = tGammat * (tRTimeS + tBetat * tROutS);
  Double_t tRSide = tRSideS;

  Double_t tRLong = tGammaz * (tRLongS + tBetaz * tDtL);
  Double_t tDt    = tGammaz * (tDtL + tBetaz * tRLongS);
	  
  tPx /= tPt;
  tPy /= tPt;
	  
  Double_t tXout  = tROut*tPx-tRSide*tPy;
  Double_t tXside = tROut*tPy+tRSide*tPx;
  Double_t tXlong = tRLong;
  Double_t tXtime = tDt;
  
  if (!(inf1->GetEmissionPoint())) {
    AliFemtoLorentzVector tPos(0,0,0,0);
    inf1->SetEmissionPoint(&tPos);
  }
  else
    inf1->SetEmissionPoint(0,0,0,0);
  if (!(inf2->GetEmissionPoint())) {
    AliFemtoLorentzVector tPos(tXout,tXside,tXlong,tXtime);
    inf2->SetEmissionPoint(&tPos);
  }
  else
    inf2->SetEmissionPoint(tXout, tXside, tXlong, tXtime);
}

//_______________________
void AliFemtoModelGausRinvFreezeOutGenerator::SetSizeInv(Double_t aSizeInv)
{
  fSizeInv = aSizeInv;
}
//_______________________
Double_t AliFemtoModelGausRinvFreezeOutGenerator::GetSizeInv() const
{
  return fSizeInv;
}
//_______________________
AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::Clone() const
{ 
  return GetGenerator(); 
}
//_______________________
inline AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::GetGenerator() const 
{ 
  AliFemtoModelFreezeOutGenerator* tModel = new AliFemtoModelGausRinvFreezeOutGenerator(*this); 
  return tModel; 
}
//_______________________
void AliFemtoModelGausRinvFreezeOutGenerator::SetSelectPrimaryFromHidden(bool aUse)
{
  fSelectPrimary = aUse;
}
Bool_t AliFemtoModelGausRinvFreezeOutGenerator::GetSelectPrimaryFromHidden()
{
  return fSelectPrimary;
}
Commit	Line	Data
ac830a3d	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoModelGausRinvFreezeOutGenerator - freeze-out ///
	4	/// coordinates generator, generating a 3D gaussian ellipsoid in LCMS ///
	5	/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
	6	/// ///
	7	////////////////////////////////////////////////////////////////////////////////
	8	#ifdef __ROOT__
	9	ClassImp(AliFemtoModelGausRinvFreezeOutGenerator, 1)
	10	#endif
	11
	12	#include "math.h"
	13	#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
65423af9	14	#include "AliFemtoModelHiddenInfo.h"
4c399116	15	#include "AliFemtoModelGlobalHiddenInfo.h"
ac830a3d	16	#include "AliFemtoLorentzVector.h"
	17
	18	//_______________________
	19	AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
1f445b96	20	fSizeInv(0),
1f445b96	21	fSelectPrimary(false)
ac830a3d	22	{
	23	// Default constructor
	24	fRandom = new TRandom2();
	25	}
	26
	27	//_______________________
	28	AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
db2d5b4b	29	AliFemtoModelFreezeOutGenerator(),
1f445b96	30	fSizeInv(0),
1f445b96	31	fSelectPrimary(false)
ac830a3d	32	{
	33	// Copy constructor
	34	fRandom = new TRandom2();
	35	SetSizeInv(aModel.GetSizeInv());
	36	}
	37	//_______________________
	38	AliFemtoModelGausRinvFreezeOutGenerator::~AliFemtoModelGausRinvFreezeOutGenerator()
	39	{
	40	if (fRandom) delete fRandom;
	41	}
	42	//_______________________
5b3c2338	43	AliFemtoModelGausRinvFreezeOutGenerator& AliFemtoModelGausRinvFreezeOutGenerator::operator=(const AliFemtoModelGausRinvFreezeOutGenerator &aModel)
	44	{
	45	if (this != &aModel) {
	46	fRandom = new TRandom2();
	47	SetSizeInv(aModel.GetSizeInv());
	48	}
	49
	50	return *this;
	51
	52	}
	53	//_______________________
ac830a3d	54	void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
ac830a3d	55	{
65423af9	56	AliFemtoModelHiddenInfo inf1 = (AliFemtoModelHiddenInfo ) aPair->Track1()->HiddenInfo();
65423af9	57	AliFemtoModelHiddenInfo inf2 = (AliFemtoModelHiddenInfo ) aPair->Track2()->HiddenInfo();
ac830a3d	58
	59	if ((!inf1) \|\| (!inf2)) { cout << "Hidden info not created! " << endl; exit(kFALSE); }
	60
1f445b96	61	if (fSelectPrimary) {
4c399116	62	const AliFemtoModelGlobalHiddenInfo infg1 = dynamic_cast<const AliFemtoModelGlobalHiddenInfo >(aPair->Track1()->HiddenInfo());
	63	const AliFemtoModelGlobalHiddenInfo infg2 = dynamic_cast<const AliFemtoModelGlobalHiddenInfo >(aPair->Track2()->HiddenInfo());
	64
	65	if ((infg1) && (infg2)) {
	66	// assume the emission point is in [cm] and try to judge if
	67	// both particles are primary
f6944668	68	Double_t dist1 = infg1->GetGlobalEmissionPoint()->Perp();
f6944668	69	Double_t dist2 = infg2->GetGlobalEmissionPoint()->Perp();
1f445b96	70
4c399116	71	if ((dist1 > 0.05) && (dist2 > 0.05)) {
	72	// At least one particle is non primary
	73	if (!(inf1->GetEmissionPoint())) {
	74	AliFemtoLorentzVector tPos(-1000,1000,-500,0);
	75	inf1->SetEmissionPoint(&tPos);
	76	}
	77	else
	78	inf1->SetEmissionPoint(-1000,1000,-500,0);
	79	if (!(inf2->GetEmissionPoint())) {
	80	AliFemtoLorentzVector tPos(fRandom->Gaus(0,1000.0),fRandom->Gaus(0,1000),fRandom->Gaus(0,1000),0.0);
	81	inf2->SetEmissionPoint(&tPos);
	82	}
	83	else
	84	inf2->SetEmissionPoint(fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), 0.0);
	85
	86	return;
1f445b96	87	}
1f445b96	88	}
	89	}
	90
	91	// Generate two particle emission points with respect
	92	// to their pair momentum
	93	// The source is the 3D Gaussian ellipsoid in the LCMS frame
	94
ac830a3d	95	// Calculate sum momenta
	96	Double_t tPx = inf1->GetTrueMomentum()->x() + inf2->GetTrueMomentum()->x();
	97	Double_t tPy = inf1->GetTrueMomentum()->y() + inf2->GetTrueMomentum()->y();
	98	Double_t tPz = inf1->GetTrueMomentum()->z() + inf2->GetTrueMomentum()->z();
	99	Double_t tM1 = inf1->GetMass();
	100	Double_t tM2 = inf2->GetMass();
f6944668	101	Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->Mag2());
f6944668	102	Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->Mag2());
ac830a3d	103	Double_t tEs = tE1 + tE2;
	104
	105	Double_t tPt = sqrt(tPxtPx + tPytPy);
	106	Double_t tMt = sqrt(tEstEs - tPztPz);
	107
	108	// Generate positions in PRF from a Gaussian
	109	Double_t tROutS = fRandom->Gaus(0,fSizeInv); // reuse of long
	110	Double_t tRSideS = fRandom->Gaus(0,fSizeInv);
	111	Double_t tRLongS = fRandom->Gaus(0,fSizeInv);
	112	Double_t tRTimeS = 0;
	113
	114	Double_t tBetat = tPt/tMt;
	115	Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
	116
	117	Double_t tBetaz = tPz/tEs;
	118	Double_t tGammaz = 1.0/sqrt(1.0-tBetaz*tBetaz);
	119
	120	Double_t tROut = tGammat * (tROutS + tBetat * tRTimeS);
	121	Double_t tDtL = tGammat * (tRTimeS + tBetat * tROutS);
	122	Double_t tRSide = tRSideS;
	123
	124	Double_t tRLong = tGammaz * (tRLongS + tBetaz * tDtL);
	125	Double_t tDt = tGammaz * (tDtL + tBetaz * tRLongS);
	126
	127	tPx /= tPt;
	128	tPy /= tPt;
	129
	130	Double_t tXout = tROuttPx-tRSidetPy;
	131	Double_t tXside = tROuttPy+tRSidetPx;
	132	Double_t tXlong = tRLong;
	133	Double_t tXtime = tDt;
	134
	135	if (!(inf1->GetEmissionPoint())) {
	136	AliFemtoLorentzVector tPos(0,0,0,0);
	137	inf1->SetEmissionPoint(&tPos);
	138	}
	139	else
	140	inf1->SetEmissionPoint(0,0,0,0);
	141	if (!(inf2->GetEmissionPoint())) {
	142	AliFemtoLorentzVector tPos(tXout,tXside,tXlong,tXtime);
	143	inf2->SetEmissionPoint(&tPos);
	144	}
	145	else
	146	inf2->SetEmissionPoint(tXout, tXside, tXlong, tXtime);
	147	}
	148
	149	//_______________________
	150	void AliFemtoModelGausRinvFreezeOutGenerator::SetSizeInv(Double_t aSizeInv)
	151	{
	152	fSizeInv = aSizeInv;
	153	}
	154	//_______________________
	155	Double_t AliFemtoModelGausRinvFreezeOutGenerator::GetSizeInv() const
	156	{
	157	return fSizeInv;
	158	}
	159	//_______________________
	160	AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::Clone() const
	161	{
	162	return GetGenerator();
	163	}
	164	//_______________________
	165	inline AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::GetGenerator() const
	166	{
167	AliFemtoModelFreezeOutGenerator* tModel = new AliFemtoModelGausRinvFreezeOutGenerator(*this);
168	return tModel;
169	}
1f445b96	170	//_______________________
	171	void AliFemtoModelGausRinvFreezeOutGenerator::SetSelectPrimaryFromHidden(bool aUse)
	172	{
	173	fSelectPrimary = aUse;
	174	}
	175	Bool_t AliFemtoModelGausRinvFreezeOutGenerator::GetSelectPrimaryFromHidden()
	176	{
	177	return fSelectPrimary;
	178	}
	179