PWG2/FEMTOSCOPY/AliFemtoUser/Model/AliFemtoModelGausRinvFreezeOutGenerator.cxx

   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"
  14 #include "Model/AliFemtoModelHiddenInfo.h"
  15 #include "AliFemtoLorentzVector.h"
  16
  17 //_______________________
  18 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
  19   fSizeInv(0)
  20 {
  21   // Default constructor
  22   fRandom = new TRandom2();
  23 }
  24
  25 //_______________________
  26 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
  27   fSizeInv(0)
  28 {
  29   // Copy constructor
  30   fRandom = new TRandom2();
  31   SetSizeInv(aModel.GetSizeInv());
  32 }
  33 //_______________________
  34 AliFemtoModelGausRinvFreezeOutGenerator::~AliFemtoModelGausRinvFreezeOutGenerator()
  35 {
  36   if (fRandom) delete fRandom;
  37 }
  38 //_______________________
  39 void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
  40 {
  41   // Generate two particle emission points with respect
  42   // to their pair momentum
  43   // The source is the 3D Gaussian ellipsoid in the LCMS frame
  44   AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->track1()->HiddenInfo();
  45   AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->track2()->HiddenInfo();
  46
  47   if ((!inf1) || (!inf2)) { cout << "Hidden info not created! "  << endl; exit(kFALSE); }
  48
  49   // Calculate sum momenta
  50   Double_t tPx = inf1->GetTrueMomentum()->x() + inf2->GetTrueMomentum()->x();
  51   Double_t tPy = inf1->GetTrueMomentum()->y() + inf2->GetTrueMomentum()->y();
  52   Double_t tPz = inf1->GetTrueMomentum()->z() + inf2->GetTrueMomentum()->z();
  53   Double_t tM1 = inf1->GetMass();
  54   Double_t tM2 = inf2->GetMass();
  55   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->mag2());
  56   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->mag2());
  57   Double_t tEs = tE1 + tE2;
  58
  59   Double_t tPt = sqrt(tPx*tPx + tPy*tPy);
  60   Double_t tMt = sqrt(tEs*tEs - tPz*tPz);
  61
  62   // Generate positions in PRF from a Gaussian
  63   Double_t tROutS =  fRandom->Gaus(0,fSizeInv); // reuse of long
  64   Double_t tRSideS = fRandom->Gaus(0,fSizeInv);
  65   Double_t tRLongS = fRandom->Gaus(0,fSizeInv);
  66   Double_t tRTimeS = 0;
  67
  68   Double_t tBetat  = tPt/tMt;
  69   Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
  70
  71   Double_t tBetaz  = tPz/tEs;
  72   Double_t tGammaz = 1.0/sqrt(1.0-tBetaz*tBetaz);
  73
  74   Double_t tROut = tGammat * (tROutS + tBetat * tRTimeS);
  75   Double_t tDtL  = tGammat * (tRTimeS + tBetat * tROutS);
  76   Double_t tRSide = tRSideS;
  77
  78   Double_t tRLong = tGammaz * (tRLongS + tBetaz * tDtL);
  79   Double_t tDt    = tGammaz * (tDtL + tBetaz * tRLongS);
  80
  81   tPx /= tPt;
  82   tPy /= tPt;
  83
  84   Double_t tXout  = tROut*tPx-tRSide*tPy;
  85   Double_t tXside = tROut*tPy+tRSide*tPx;
  86   Double_t tXlong = tRLong;
  87   Double_t tXtime = tDt;
  88
  89   if (!(inf1->GetEmissionPoint())) {
  90     AliFemtoLorentzVector tPos(0,0,0,0);
  91     inf1->SetEmissionPoint(&tPos);
  92   }
  93   else
  94     inf1->SetEmissionPoint(0,0,0,0);
  95   if (!(inf2->GetEmissionPoint())) {
  96     AliFemtoLorentzVector tPos(tXout,tXside,tXlong,tXtime);
  97     inf2->SetEmissionPoint(&tPos);
  98   }
  99   else
 100     inf2->SetEmissionPoint(tXout, tXside, tXlong, tXtime);
 101 }
 102
 103 //_______________________
 104 void AliFemtoModelGausRinvFreezeOutGenerator::SetSizeInv(Double_t aSizeInv)
 105 {
 106   fSizeInv = aSizeInv;
 107 }
 108 //_______________________
 109 Double_t AliFemtoModelGausRinvFreezeOutGenerator::GetSizeInv() const
 110 {
 111   return fSizeInv;
 112 }
 113 //_______________________
 114 AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::Clone() const
 115 {
 116   return GetGenerator();
 117 }
 118 //_______________________
 119 inline AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::GetGenerator() const
 120 {
 121   AliFemtoModelFreezeOutGenerator* tModel = new AliFemtoModelGausRinvFreezeOutGenerator(*this);
 122   return tModel;
 123 }