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 "AliFemtoModelGlobalHiddenInfo.h"
  16 #include "AliFemtoLorentzVector.h"
  17
  18 //_______________________
  19 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
  20   fSizeInv(0),
  21   fSelectPrimary(false)
  22 {
  23   // Default constructor
  24   fRandom = new TRandom2();
  25 }
  26
  27 //_______________________
  28 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
  29   AliFemtoModelFreezeOutGenerator(),
  30   fSizeInv(0),
  31   fSelectPrimary(false)
  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 //_______________________
  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 //_______________________
  54 void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
  55 {
  56   AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
  57   AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
  58
  59   if ((!inf1) || (!inf2)) { cout << "Hidden info not created! "  << endl; exit(kFALSE); }
  60
  61   if (fSelectPrimary) {
  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
  68       Double_t dist1 = infg1->GetGlobalEmissionPoint()->Perp();
  69       Double_t dist2 = infg2->GetGlobalEmissionPoint()->Perp();
  70
  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;
  87       }
  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
  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();
 101   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->Mag2());
 102   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->Mag2());
 103   Double_t tEs = tE1 + tE2;
 104
 105   Double_t tPt = sqrt(tPx*tPx + tPy*tPy);
 106   Double_t tMt = sqrt(tEs*tEs - tPz*tPz);
 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  = tROut*tPx-tRSide*tPy;
 131   Double_t tXside = tROut*tPy+tRSide*tPx;
 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 }
 170 //_______________________
 171 void AliFemtoModelGausRinvFreezeOutGenerator::SetSelectPrimaryFromHidden(bool aUse)
 172 {
 173   fSelectPrimary = aUse;
 174 }
 175 Bool_t AliFemtoModelGausRinvFreezeOutGenerator::GetSelectPrimaryFromHidden()
 176 {
 177   return fSelectPrimary;
 178 }
 179