PWG2/FEMTOSCOPY/AliFemtoUser/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 "AliFemtoModelHiddenInfo.h"
  15 #include "AliFemtoLorentzVector.h"
  16
  17 //_______________________
  18 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
  19   fSizeInv(0),
  20   fSelectPrimary(false)
  21 {
  22   // Default constructor
  23   fRandom = new TRandom2();
  24 }
  25
  26 //_______________________
  27 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
  28   AliFemtoModelFreezeOutGenerator(),
  29   fSizeInv(0),
  30   fSelectPrimary(false)
  31 {
  32   // Copy constructor
  33   fRandom = new TRandom2();
  34   SetSizeInv(aModel.GetSizeInv());
  35 }
  36 //_______________________
  37 AliFemtoModelGausRinvFreezeOutGenerator::~AliFemtoModelGausRinvFreezeOutGenerator()
  38 {
  39   if (fRandom) delete fRandom;
  40 }
  41 //_______________________
  42 void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
  43 {
  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   if (fSelectPrimary) {
  50     // assume the emission point is in [cm] and try to judge if
  51     // both particles are primary
  52     Double_t dist1 = inf1->GetEmissionPoint()->perp();
  53     Double_t dist2 = inf2->GetEmissionPoint()->perp();
  54
  55     if ((dist1 > 0.005) && (dist2 > 0.005)) {
  56       // At least one particle is non primary
  57       if (!(inf1->GetEmissionPoint())) {
  58         AliFemtoLorentzVector tPos(-1000,1000,-500,0);
  59         inf1->SetEmissionPoint(&tPos);
  60       }
  61       else
  62         inf1->SetEmissionPoint(-1000,1000,-500,0);
  63       if (!(inf2->GetEmissionPoint())) {
  64         AliFemtoLorentzVector tPos(fRandom->Gaus(0,1000.0),fRandom->Gaus(0,1000),fRandom->Gaus(0,1000),0.0);
  65         inf2->SetEmissionPoint(&tPos);
  66       }
  67       else
  68         inf2->SetEmissionPoint(fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), 0.0);
  69
  70       return;
  71     }
  72   }
  73
  74   // Generate two particle emission points with respect
  75   // to their pair momentum
  76   // The source is the 3D Gaussian ellipsoid in the LCMS frame
  77
  78   // Calculate sum momenta
  79   Double_t tPx = inf1->GetTrueMomentum()->x() + inf2->GetTrueMomentum()->x();
  80   Double_t tPy = inf1->GetTrueMomentum()->y() + inf2->GetTrueMomentum()->y();
  81   Double_t tPz = inf1->GetTrueMomentum()->z() + inf2->GetTrueMomentum()->z();
  82   Double_t tM1 = inf1->GetMass();
  83   Double_t tM2 = inf2->GetMass();
  84   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->mag2());
  85   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->mag2());
  86   Double_t tEs = tE1 + tE2;
  87
  88   Double_t tPt = sqrt(tPx*tPx + tPy*tPy);
  89   Double_t tMt = sqrt(tEs*tEs - tPz*tPz);
  90
  91   // Generate positions in PRF from a Gaussian
  92   Double_t tROutS =  fRandom->Gaus(0,fSizeInv); // reuse of long
  93   Double_t tRSideS = fRandom->Gaus(0,fSizeInv);
  94   Double_t tRLongS = fRandom->Gaus(0,fSizeInv);
  95   Double_t tRTimeS = 0;
  96
  97   Double_t tBetat  = tPt/tMt;
  98   Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
  99
 100   Double_t tBetaz  = tPz/tEs;
 101   Double_t tGammaz = 1.0/sqrt(1.0-tBetaz*tBetaz);
 102
 103   Double_t tROut = tGammat * (tROutS + tBetat * tRTimeS);
 104   Double_t tDtL  = tGammat * (tRTimeS + tBetat * tROutS);
 105   Double_t tRSide = tRSideS;
 106
 107   Double_t tRLong = tGammaz * (tRLongS + tBetaz * tDtL);
 108   Double_t tDt    = tGammaz * (tDtL + tBetaz * tRLongS);
 109
 110   tPx /= tPt;
 111   tPy /= tPt;
 112
 113   Double_t tXout  = tROut*tPx-tRSide*tPy;
 114   Double_t tXside = tROut*tPy+tRSide*tPx;
 115   Double_t tXlong = tRLong;
 116   Double_t tXtime = tDt;
 117
 118   if (!(inf1->GetEmissionPoint())) {
 119     AliFemtoLorentzVector tPos(0,0,0,0);
 120     inf1->SetEmissionPoint(&tPos);
 121   }
 122   else
 123     inf1->SetEmissionPoint(0,0,0,0);
 124   if (!(inf2->GetEmissionPoint())) {
 125     AliFemtoLorentzVector tPos(tXout,tXside,tXlong,tXtime);
 126     inf2->SetEmissionPoint(&tPos);
 127   }
 128   else
 129     inf2->SetEmissionPoint(tXout, tXside, tXlong, tXtime);
 130 }
 131
 132 //_______________________
 133 void AliFemtoModelGausRinvFreezeOutGenerator::SetSizeInv(Double_t aSizeInv)
 134 {
 135   fSizeInv = aSizeInv;
 136 }
 137 //_______________________
 138 Double_t AliFemtoModelGausRinvFreezeOutGenerator::GetSizeInv() const
 139 {
 140   return fSizeInv;
 141 }
 142 //_______________________
 143 AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::Clone() const
 144 {
 145   return GetGenerator();
 146 }
 147 //_______________________
 148 inline AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::GetGenerator() const
 149 {
 150   AliFemtoModelFreezeOutGenerator* tModel = new AliFemtoModelGausRinvFreezeOutGenerator(*this);
 151   return tModel;
 152 }
 153 //_______________________
 154 void AliFemtoModelGausRinvFreezeOutGenerator::SetSelectPrimaryFromHidden(bool aUse)
 155 {
 156   fSelectPrimary = aUse;
 157 }
 158 Bool_t AliFemtoModelGausRinvFreezeOutGenerator::GetSelectPrimaryFromHidden()
 159 {
 160   return fSelectPrimary;
 161 }
 162